2 #include <JANA/JApplication.h>
27 #include <TDirectory.h>
31 #include <TFitResult.h>
32 #include <TFitResultPtr.h>
36 #include <JANA/JApplication.h>
37 #include <JANA/JFactory.h>
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");
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)");
97 if(CALC_NEW_CONSTANTS_BEAM) {
106 TDirectory *
main = gDirectory;
107 gDirectory->mkdir(
"FCAL_LED_shifts")->cd();
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);
115 m_totalEnergy =
new TH1I(
"fcal_total_energy",
"Total Energy", 500, 0, 10000);
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) {
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)) );
164 vector< const DFCALGeometry* > geomVec;
165 eventLoop->Get( geomVec );
167 if( geomVec.size() != 1 ){
169 cerr <<
"No geometry accessbile." << endl;
170 return RESOURCE_UNAVAILABLE;
173 m_fcalGeom = geomVec[0];
177 vector< const DTranslationTable* > ttabVec;
178 eventLoop->Get(ttabVec);
184 m_runnumber = runnumber;
186 eventLoop->GetCalib(
"/FCAL/ADC_Offsets", old_ADCoffsets);
195 uint64_t eventnumber)
200 vector< const DFCALHit* > hits;
201 eventLoop->Get( hits );
205 map< int, pair<uint32_t,uint32_t> > crate_slot_map;
209 double total_energy = 0.;
210 for(
auto hit : hits ) {
211 total_energy += hit->E;
214 for(
auto hit : hits ) {
217 int channel_index = m_fcalGeom->channel(hit->row, hit->column);
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);
229 crate_slot_map[channel_index] = crate_slot;
238 japp->RootFillLock(
this);
240 m_totalEnergy->Fill(total_energy);
243 for(
auto hit : hits ) {
246 int channel_index = m_fcalGeom->channel(hit->row, hit->column);
248 pair<uint32_t,uint32_t> crate_slot = crate_slot_map[channel_index];
251 if(crate_slot.first == 0) {
256 m_crateTimes[crate_slot.first]->Fill(hit->t);
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));
265 m_slotTimes[crate_slot]->Fill(hit->t);
266 m_channelTimes[channel_index]->Fill(hit->t);
273 japp->RootFillUnLock(
this);
284 if( fabs((time - reference_time)-4) < 1. )
286 else if( fabs((time - reference_time)+4) < 1. )
303 if(CALC_NEW_CONSTANTS_BEAM) {
306 auto ref_file =
new TFile(REFERENCE_FILE_NAME.c_str());
308 ofstream outf(Form(
"fcal_adc_offsets_r%i.txt", m_runnumber));
311 auto fgaus =
new TF1(
"fgaus",
"gaus", -200, 200);
315 map<uint32_t, double> crate_shifts;
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;
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);
329 if(m_crateTimes.find(crate) == m_crateTimes.end())
continue;
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);
336 crate_shifts[crate] = adc_shift;
338 cout <<
"crate " << crate <<
" ADC shift = " << (int)adc_shift
339 <<
" reference time = " << reference_time <<
" time = " << time << endl;
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;
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;
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);
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);
364 slot_shifts[slotTimeEntry.first] = adc_shift;
366 cout <<
"crate " << crate <<
"slot " << slot <<
" ADC shift = " << (int)adc_shift
367 <<
" reference time = " << reference_time <<
" time = " << time << endl;
374 int row = m_fcalGeom->row(i);
375 int col = m_fcalGeom->column(i);
382 auto daq_index = m_ttab->GetDAQIndex(channel_info);
383 pair<uint32_t,uint32_t> crate_slot(daq_index.rocid,daq_index.slot);
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];
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;
404 double max = ref_hist->GetBinCenter(ref_hist->GetMaximumBin());
405 TFitResultPtr ref_fr = ref_hist->Fit(fgaus,
"SQ",
"", max - 2.5, max + 2.5);
407 double reference_time = ref_fr->Parameter(1);
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);
412 double time = fr->Parameter(1);
420 m_fadcShifts->Fill(col,row,adc_shift);
422 outf << (old_ADCoffsets[i] + adc_shift) << endl;
432 if(CALC_NEW_CONSTANTS_LED) {
435 auto ref_file =
new TFile(REFERENCE_FILE_NAME.c_str());
437 ofstream outf(Form(
"fcal_adc_offsets_r%i.txt", m_runnumber));
440 map<uint32_t, double> crate_shifts;
444 auto ref_hist = (TH1I*)ref_file->Get(Form(
"FCAL_LED_shifts/crate_times_%i",crate));
446 double reference_time = ref_hist->GetMean();
447 double time = m_crateTimes[crate]->GetMean();
450 crate_shifts[crate] = adc_shift;
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;
462 auto ref_hist = (TH1I*)ref_file->Get(Form(
"FCAL_LED_shifts/slot_times_c%i_s%i",crate,slot));
464 double reference_time = ref_hist->GetMean();
465 double time = slotTimeEntry.second->GetMean();
468 slot_shifts[slotTimeEntry.first] = adc_shift;
478 int row = m_fcalGeom->row(i);
479 int col = m_fcalGeom->column(i);
486 auto daq_index = m_ttab->GetDAQIndex(channel_info);
487 pair<uint32_t,uint32_t> crate_slot(daq_index.rocid,daq_index.slot);
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];
501 auto ref_hist = (TH1I*)ref_file->Get(Form(
"FCAL_LED_shifts/channel_times_%i",i));
503 double reference_time = ref_hist->GetMean();
504 double time = m_channelTimes[i]->GetMean();
509 m_fadcShifts->Fill(col,row,adc_shift);
511 outf << (old_ADCoffsets[i] + adc_shift) << endl;
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.
~JEventProcessor_FCAL_LED_shifts()
jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
JEventProcessor_FCAL_LED_shifts()
jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
int main(int argc, char *argv[])