Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JEventProcessor_FCAL_LED_shifts.cc
Go to the documentation of this file.
2 #include <JANA/JApplication.h>
3 #include "FCAL/DFCALShower.h"
4 #include "FCAL/DFCALGeometry.h"
5 #include "FCAL/DFCALHit.h"
6 #include "FCAL/DFCALDigiHit.h"
7 #include "FCAL/DFCALCluster.h"
9 #include "PID/DVertex.h"
10 #include "DVector3.h"
12 #include <TTree.h>
13 #include "DVector3.h"
14 #include "PID/DParticleID.h"
15 #include "TRIGGER/DTrigger.h"
16 #include "GlueX.h"
17 #include <vector>
18 #include <map>
19 #include <deque>
20 #include <string>
21 #include <iostream>
22 #include <algorithm>
23 #include <iostream>
24 #include <fstream>
25 #include <stdio.h>
26 #include <stdlib.h>
27 #include <TDirectory.h>
28 #include <TH1I.h>
29 #include <TH2F.h>
30 #include <TFile.h>
31 #include <TFitResult.h>
32 #include <TFitResultPtr.h>
33 
34 using namespace jana;
35 
36 #include <JANA/JApplication.h>
37 #include <JANA/JFactory.h>
38 
39 const int numChannels = 2800;
40 
41 const int numCrates = 12;
42 const int firstCrate = 11;
43 const int numSlots = 16;
44 const int firstSlot = 3;
45 
46 // Define Histograms
47 
48 
49 extern "C"{
50  void InitPlugin(JApplication *app){
51  InitJANAPlugin(app);
52  app->AddProcessor(new JEventProcessor_FCAL_LED_shifts());
53  }
54 } // "C"
55 
56 
57 //------------------
58 // JEventProcessor_FCAL_LED_shifts (Constructor)
59 //------------------
61 {
62 
63 }
64 
65 //------------------
66 // ~JEventProcessor_FCAL_LED_shifts (Destructor)
67 //------------------
69 {
70 
71 }
72 
73 //------------------
74 // init
75 //------------------
77 {
78  // Set parameters
79  CALC_NEW_CONSTANTS_LED = true;
80  gPARMS->SetDefaultParameter("FCAL_SHIFT:CALC_NEW_CONSTANTS_LED", CALC_NEW_CONSTANTS_LED, "True if we should calculate new offsets based on a reference file for LED data");
81  CALC_NEW_CONSTANTS_BEAM = false;
82  gPARMS->SetDefaultParameter("FCAL_SHIFT:CALC_NEW_CONSTANTS_BEAM", CALC_NEW_CONSTANTS_BEAM, "True if we should calculate new offsets based on a reference file for beam data");
83  REFERENCE_FILE_NAME = "hd_root-r42485-led.root";
84  gPARMS->SetDefaultParameter("FCAL_SHIFT:REFERENCE_FILE_NAME", REFERENCE_FILE_NAME, "Reference file for new offsets");
85 
86  FCAL_TOTAL_ENERGY_HI = 9000.;
87  gPARMS->SetDefaultParameter("FCAL_SHIFT:FCAL_TOTAL_ENERGY_HI", FCAL_TOTAL_ENERGY_HI, "Total event energy cut (hi level)");
88  FCAL_TOTAL_ENERGY_LO = 8500.;
89  gPARMS->SetDefaultParameter("FCAL_SHIFT:FCAL_TOTAL_ENERGY_LO", FCAL_TOTAL_ENERGY_LO, "Total event energy cut (lo level)");
90 
91 
92  // the histogram limits should be different for LED and beam data events
93  NBINS_TIME=100;
94  TIME_MIN=60.;
95  TIME_MAX=110.;
96 
97  if(CALC_NEW_CONSTANTS_BEAM) {
98  TIME_MIN=-20.;
99  TIME_MAX=80.;
100  }
101 
102  // This is called once at program startup. If you are creating
103  // and filling historgrams in this plugin, you should lock the
104  // ROOT mutex like this:
105  //
106  TDirectory *main = gDirectory;
107  gDirectory->mkdir("FCAL_LED_shifts")->cd();
108 
109  // for monitoring of results
110  m_fadcShifts = new TH2I("fcal_fadc_shifts", "ADC Shifts", 60, 0, 60, 60, 0, 60);
111  m_fadcShifts->GetXaxis()->SetTitle("column");
112  m_fadcShifts->GetYaxis()->SetTitle("row");
113  m_fadcShifts->GetZaxis()->SetRangeUser(-8,8);
114 
115  m_totalEnergy = new TH1I("fcal_total_energy", "Total Energy", 500, 0, 10000);
116 
117  for (int i = 0; i < numCrates; ++i) {
118  uint32_t crate = firstCrate + i;
119  m_crateTimes[crate] = (new TH1I(Form("crate_times_%i",crate),Form("Hit Times for Crate %i",crate),NBINS_TIME,TIME_MIN,TIME_MAX));
120  for (int i = 0; i < numSlots; ++i) {
121  uint32_t slot = firstSlot + i;
122  //pair<uint32_t,uint32_t> crate_slot(crate,slot);
123  //m_slotTimes[crate_slot] = (new TH1I(Form("slot_times_c%i_s%i",crate,slot),
124  // Form("Hit Times for Crate %i, Slot %i",crate,slot),200,60.,110.));
125  }
126  }
127 
128  for (int i = 0; i < numChannels; ++i) {
129  m_channelTimes.push_back( (new TH1I(Form("channel_times_%i",i),Form("Hit Times for Channel %i",i),NBINS_TIME,TIME_MIN,TIME_MAX)) );
130  }
131 
132 
133  main->cd();
134 
135 
136  return NOERROR;
137 }
138 
139 //------------------
140 // brun
141 //------------------
142 jerror_t JEventProcessor_FCAL_LED_shifts::brun(JEventLoop *eventLoop,
143  int32_t runnumber)
144 {
145  /*
146  // get the FCAL z position from the global geometry interface
147  DApplication *dapp =
148  dynamic_cast<DApplication*>(eventLoop->GetJApplication());
149  const DGeometry *geom = dapp->GetDGeometry(runnumber);
150  if( geom ) {
151 
152  geom->GetFCALZ( m_FCALfront );
153  }
154  else{
155 
156  cerr << "No geometry accessbile." << endl;
157  return RESOURCE_UNAVAILABLE;
158  }
159  */
160 
161  // WARNING: THIS IS SUPER DANGEROUS
162  // FIGURE OUT A WAY TO SAVE THIS INFO FOR erun()
163  // we need an FCAL Geometry object
164  vector< const DFCALGeometry* > geomVec;
165  eventLoop->Get( geomVec );
166 
167  if( geomVec.size() != 1 ){
168 
169  cerr << "No geometry accessbile." << endl;
170  return RESOURCE_UNAVAILABLE;
171  }
172 
173  m_fcalGeom = geomVec[0];
174 
175 
176  // Load the translation table
177  vector< const DTranslationTable* > ttabVec;
178  eventLoop->Get(ttabVec);
179 
180  m_ttab = ttabVec[0];
181 
182  // save this info - not terribly thread safe, but this plugin
183  // should only be used on one run at once
184  m_runnumber = runnumber;
185  // load old offsets
186  eventLoop->GetCalib("/FCAL/ADC_Offsets", old_ADCoffsets);
187 
188  return NOERROR;
189 }
190 
191 //------------------
192 // evnt
193 //------------------
194 jerror_t JEventProcessor_FCAL_LED_shifts::evnt(JEventLoop *eventLoop,
195  uint64_t eventnumber)
196 {
197 
198 
199 
200  vector< const DFCALHit* > hits;
201  eventLoop->Get( hits );
202 
203 
204  // save this mapping outside of the locks below, since the index lookup is pretty slow
205  map< int, pair<uint32_t,uint32_t> > crate_slot_map;
206 
207  // look at the total energy first to pick up a particle LED color
208  // since each color has different timing. pick the most energetic to analyze
209  double total_energy = 0.;
210  for( auto hit : hits ) {
211  total_energy += hit->E;
212  }
213 
214  for( auto hit : hits ) {
215  //if(total_energy > FCAL_TOTAL_ENERGY_LO && total_energy < FCAL_TOTAL_ENERGY_HI) {
216  // convert (row,col) to channel index
217  int channel_index = m_fcalGeom->channel(hit->row, hit->column);
218 
219  // look up crate/slot info
220  DTranslationTable::DChannelInfo channel_info;
221  channel_info.det_sys = DTranslationTable::FCAL;
222  channel_info.fcal.row = hit->row;
223  channel_info.fcal.col = hit->column;
224  auto daq_index = m_ttab->GetDAQIndex(channel_info);
225  pair<uint32_t,uint32_t> crate_slot(daq_index.rocid,daq_index.slot);
226 
227  //cerr << "crate = " << crate_slot.first << " slot = " << crate_slot.second << endl;
228 
229  crate_slot_map[channel_index] = crate_slot;
230  //}
231  }
232 
233  //cerr << total_energy << endl;
234 
235 
236  // FILL HISTOGRAMS
237  // Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock
238  japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
239 
240  m_totalEnergy->Fill(total_energy);
241 
242  //if(total_energy > FCAL_TOTAL_ENERGY_LO && total_energy < FCAL_TOTAL_ENERGY_HI) {
243  for( auto hit : hits ) {
244 
245  // convert (row,col) to channel index
246  int channel_index = m_fcalGeom->channel(hit->row, hit->column);
247 
248  pair<uint32_t,uint32_t> crate_slot = crate_slot_map[channel_index];
249  //cerr << "crate = " << crate_slot.first << " slot = " << crate_slot.second << endl;
250 
251  if(crate_slot.first == 0) { // ignore these hits
252  continue;
253  }
254 
255  // fill histograms
256  m_crateTimes[crate_slot.first]->Fill(hit->t);
257  // make these on demand, since the number of slots is different in each crate
258  if(m_slotTimes.find(crate_slot) == m_slotTimes.end()) {
259  TDirectory *main = gDirectory;
260  gDirectory->cd("FCAL_LED_shifts");
261  m_slotTimes[crate_slot] = (new TH1I(Form("slot_times_c%i_s%i",crate_slot.first,crate_slot.second),
262  Form("Hit Times for Crate %i, Slot %i",crate_slot.first,crate_slot.second),NBINS_TIME/2,TIME_MIN,TIME_MAX));
263  main->cd();
264  }
265  m_slotTimes[crate_slot]->Fill(hit->t);
266  m_channelTimes[channel_index]->Fill(hit->t);
267 
268  //cout << " crate = " << daq_index.rocid << " slot = " << daq_index.slot
269  // << " time = " << hit->t << endl;
270  }
271  //}
272 
273  japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK
274 
275 
276  return NOERROR;
277 }
278 
279 static double CalcADCShift(double reference_time, double time) {
280  // Assume, for now, the simplest case: fADC250 times can shift by +/- 4 ns
281  // Hopefully it stays that way
282  // Use some rough tolerance, return the appropriate correction factor
283 
284  if( fabs((time - reference_time)-4) < 1. )
285  return -4.;
286  else if( fabs((time - reference_time)+4) < 1. )
287  return 4.;
288 
289  return 0.;
290 }
291 
292 
293 //------------------
294 // erun
295 //------------------
297 {
298  // This is called whenever the run number changes, before it is
299  // changed to give you a chance to clean up before processing
300  // events from the next run number.
301 
302 
303  if(CALC_NEW_CONSTANTS_BEAM) {
304  // calculate time shifts
305  //cerr << "opening " << REFERENCE_FILE_NAME << endl;
306  auto ref_file = new TFile(REFERENCE_FILE_NAME.c_str());
307 
308  ofstream outf(Form("fcal_adc_offsets_r%i.txt", m_runnumber));
309 
310  // Have to define an object to pass into TH1::Fit() below since ROOT is being weird
311  auto fgaus = new TF1("fgaus", "gaus", -200, 200);
312 
313 
314  // calculate crate shifts
315  map<uint32_t, double> crate_shifts;
316  for (int i = 0; i < numCrates; ++i) {
317  uint32_t crate = firstCrate + i;
318 
319  auto ref_hist = (TH1I*)ref_file->Get(Form("FCAL_LED_shifts/crate_times_%i",crate));
320  if(ref_hist == NULL) {
321  cout << "skipping crate " << crate << " ..." <<endl;
322  continue;
323  }
324  // reference time
325  double max = ref_hist->GetBinCenter(ref_hist->GetMaximumBin());
326  TFitResultPtr ref_fr = ref_hist->Fit(fgaus, "SQ", "", max - 2.5, max + 2.5);
327  double reference_time = ref_fr->Parameter(1);
328 
329  if(m_crateTimes.find(crate) == m_crateTimes.end()) continue;
330 
331  max = m_crateTimes[crate]->GetBinCenter(m_crateTimes[crate]->GetMaximumBin());
332  TFitResultPtr fr = m_crateTimes[crate]->Fit(fgaus, "SQ", "", max - 2.5, max + 2.5);
333  double time = fr->Parameter(1);
334 
335  double adc_shift = CalcADCShift(reference_time, time);
336  crate_shifts[crate] = adc_shift;
337 
338  cout << "crate " << crate << " ADC shift = " << (int)adc_shift
339  << " reference time = " << reference_time << " time = " << time << endl;
340  }
341 
342  // calculate slot shifts
343  map< pair<uint32_t,uint32_t>, double> slot_shifts;
344  for (auto &slotTimeEntry : m_slotTimes) {
345  uint32_t crate = slotTimeEntry.first.first;
346  uint32_t slot = slotTimeEntry.first.second;
347 
348  auto ref_hist = (TH1I*)ref_file->Get(Form("FCAL_LED_shifts/slot_times_c%i_s%i",crate,slot));
349  if(ref_hist == NULL) {
350  cout << "skipping crate " << crate << " slot " << slot << " ..." <<endl;
351  continue;
352  }
353 
354  // reference time
355  double max = ref_hist->GetBinCenter(ref_hist->GetMaximumBin());
356  TFitResultPtr ref_fr = ref_hist->Fit(fgaus, "SQ", "", max - 2.5, max + 2.5);
357  double reference_time = ref_fr->Parameter(1);
358 
359  max = slotTimeEntry.second->GetBinCenter(slotTimeEntry.second->GetMaximumBin());
360  TFitResultPtr fr = slotTimeEntry.second->Fit(fgaus, "SQ", "", max - 2.5, max + 2.5);
361  double time = fr->Parameter(1);
362 
363  double adc_shift = CalcADCShift(reference_time, time);
364  slot_shifts[slotTimeEntry.first] = adc_shift;
365 
366  cout << "crate " << crate << "slot " << slot << " ADC shift = " << (int)adc_shift
367  << " reference time = " << reference_time << " time = " << time << endl;
368 
369  }
370 
371  // calculate channel shifts
372  for (int i = 0; i < numChannels; ++i) {
373  // figure out detector location and indexing
374  int row = m_fcalGeom->row(i);
375  int col = m_fcalGeom->column(i);
376 
377  // look up crate/slot info
378  DTranslationTable::DChannelInfo channel_info;
379  channel_info.det_sys = DTranslationTable::FCAL;
380  channel_info.fcal.row = row;
381  channel_info.fcal.col = col;
382  auto daq_index = m_ttab->GetDAQIndex(channel_info);
383  pair<uint32_t,uint32_t> crate_slot(daq_index.rocid,daq_index.slot);
384 
385  //m_channelTimes.push_back( (new TH1I(Form("channel_times_%i",i),Form("Hit Times for Channel %i",i),50,60.,110.)) );
386 
387  // determine the shifts for individual channels in the following way:
388  // 1) if the whole crate has shifted, use this value
389  // 2) if the whole slot has shifted, use this value
390  // 3) use the individual channel
391  double adc_shift = 0.;
392  if(static_cast<int>(crate_shifts[daq_index.rocid]) != 0) {
393  adc_shift = crate_shifts[daq_index.rocid];
394  } else if(static_cast<int>(slot_shifts[crate_slot]) != 0) {
395  adc_shift = slot_shifts[crate_slot];
396  } else {
397  auto ref_hist = (TH1I*)ref_file->Get(Form("FCAL_LED_shifts/channel_times_%i",i));
398  if(ref_hist == NULL) {
399  cout << "skipping channel " << i << " ..." <<endl;
400  continue;
401  }
402 
403  // reference time
404  double max = ref_hist->GetBinCenter(ref_hist->GetMaximumBin());
405  TFitResultPtr ref_fr = ref_hist->Fit(fgaus, "SQ", "", max - 2.5, max + 2.5);
406  if(ref_fr>=0) { // make sure the fits converge...
407  double reference_time = ref_fr->Parameter(1);
408 
409  max = m_channelTimes[i]->GetBinCenter(m_channelTimes[i]->GetMaximumBin());
410  TFitResultPtr fr = m_channelTimes[i]->Fit(fgaus, "SQ", "", max - 2.5, max + 2.5);
411  if(fr>=0) {
412  double time = fr->Parameter(1);
413 
414  adc_shift = CalcADCShift(reference_time, time);
415  }
416  }
417  }
418 
419  // histogram the results for monitoring purposes
420  m_fadcShifts->Fill(col,row,adc_shift);
421 
422  outf << (old_ADCoffsets[i] + adc_shift) << endl;
423  }
424 
425 
426  ref_file->Close();
427  outf.close();
428  delete fgaus;
429 
430  }
431 
432  if(CALC_NEW_CONSTANTS_LED) {
433  // calculate time shifts
434  //cerr << "opening " << REFERENCE_FILE_NAME << endl;
435  auto ref_file = new TFile(REFERENCE_FILE_NAME.c_str());
436 
437  ofstream outf(Form("fcal_adc_offsets_r%i.txt", m_runnumber));
438 
439  // calculate crate shifts
440  map<uint32_t, double> crate_shifts;
441  for (int i = 0; i < numCrates; ++i) {
442  uint32_t crate = firstCrate + i;
443 
444  auto ref_hist = (TH1I*)ref_file->Get(Form("FCAL_LED_shifts/crate_times_%i",crate));
445  // reference time
446  double reference_time = ref_hist->GetMean();
447  double time = m_crateTimes[crate]->GetMean();
448 
449  double adc_shift = CalcADCShift(reference_time, time);
450  crate_shifts[crate] = adc_shift;
451 
452  //cout << "crate " << crate << " ADC shift = " << (int)adc_shift
453  // << " reference time = " << reference_time << " time = " << time << endl;
454  }
455 
456  // calculate slot shifts
457  map< pair<uint32_t,uint32_t>, double> slot_shifts;
458  for (auto &slotTimeEntry : m_slotTimes) {
459  uint32_t crate = slotTimeEntry.first.first;
460  uint32_t slot = slotTimeEntry.first.second;
461 
462  auto ref_hist = (TH1I*)ref_file->Get(Form("FCAL_LED_shifts/slot_times_c%i_s%i",crate,slot));
463  // reference time
464  double reference_time = ref_hist->GetMean();
465  double time = slotTimeEntry.second->GetMean();
466 
467  double adc_shift = CalcADCShift(reference_time, time);
468  slot_shifts[slotTimeEntry.first] = adc_shift;
469 
470  //cout << "crate " << crate << "slot " << slot << " ADC shift = " << (int)adc_shift
471  // << " reference time = " << reference_time << " time = " << time << endl;
472 
473  }
474 
475  // calculate channel shifts
476  for (int i = 0; i < numChannels; ++i) {
477  // figure out detector location and indexing
478  int row = m_fcalGeom->row(i);
479  int col = m_fcalGeom->column(i);
480 
481  // look up crate/slot info
482  DTranslationTable::DChannelInfo channel_info;
483  channel_info.det_sys = DTranslationTable::FCAL;
484  channel_info.fcal.row = row;
485  channel_info.fcal.col = col;
486  auto daq_index = m_ttab->GetDAQIndex(channel_info);
487  pair<uint32_t,uint32_t> crate_slot(daq_index.rocid,daq_index.slot);
488 
489  //m_channelTimes.push_back( (new TH1I(Form("channel_times_%i",i),Form("Hit Times for Channel %i",i),50,60.,110.)) );
490 
491  // determine the shifts for individual channels in the following way:
492  // 1) if the whole crate has shifted, use this value
493  // 2) if the whole slot has shifted, use this value
494  // 3) use the individual channel
495  double adc_shift = 0.;
496  if(static_cast<int>(crate_shifts[daq_index.rocid]) != 0) {
497  adc_shift = crate_shifts[daq_index.rocid];
498  } else if(static_cast<int>(slot_shifts[crate_slot]) != 0) {
499  adc_shift = slot_shifts[crate_slot];
500  } else {
501  auto ref_hist = (TH1I*)ref_file->Get(Form("FCAL_LED_shifts/channel_times_%i",i));
502  // reference time
503  double reference_time = ref_hist->GetMean();
504  double time = m_channelTimes[i]->GetMean();
505  adc_shift = CalcADCShift(reference_time, time);
506  }
507 
508  // histogram the results for monitoring purposes
509  m_fadcShifts->Fill(col,row,adc_shift);
510 
511  outf << (old_ADCoffsets[i] + adc_shift) << endl;
512  }
513 
514 
515  ref_file->Close();
516  outf.close();
517 
518  }
519 
520 
521  return NOERROR;
522 }
523 
524 //------------------
525 // fini
526 //------------------
528 {
529 
530  // Called before program exit after event processing is finished.
531  return NOERROR;
532 }
533 
534 
static double CalcADCShift(double reference_time, double time)
jerror_t init(void)
Called once at program start.
jerror_t fini(void)
Called after last event of last event source has been processed.
JApplication * japp
const int firstSlot
const int numChannels
InitPlugin_t InitPlugin
const int numCrates
jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
const int numSlots
const int firstCrate
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
int main(int argc, char *argv[])
Definition: gendoc.cc:6