10 #include <ecl/calibration/eclMergingCrystalTimingAlgorithm.h>
11 #include <ecl/dbobjects/ECLCrystalCalib.h>
12 #include <ecl/dbobjects/ECLReferenceCrystalPerCrateCalib.h>
13 #include <ecl/utility/ECLChannelMapper.h>
14 #include <ecl/digitization/EclConfiguration.h>
18 #include "TDirectory.h"
23 using namespace Calibration;
26 eclMergingCrystalTimingAlgorithm::eclMergingCrystalTimingAlgorithm():
28 readPrevCrysPayload(false),
29 m_ECLCrystalTimeOffsetBhabha(
"ECLCrystalTimeOffsetBhabha"),
30 m_ECLCrystalTimeOffsetBhabhaGamma(
"ECLCrystalTimeOffsetBhabhaGamma"),
31 m_ECLCrystalTimeOffsetCosmic(
"ECLCrystalTimeOffsetCosmic"),
32 m_ECLReferenceCrystalPerCrateCalib(
"ECLReferenceCrystalPerCrateCalib"),
33 m_ECLCrateTimeOffset(
"ECLCrateTimeOffset")
36 "Perform time calibration of ecl crystals by combining previous values from the DB for different calibrations."
50 ExpRun chosenRun = runs.front();
51 B2INFO(
"merging using the ExpRun (" << chosenRun.second <<
"," << chosenRun.first <<
")");
64 bool allObjectsFound =
true;
69 allObjectsFound =
false;
70 B2ERROR(
"No valid DBObject found for 'ECLCrystalTimeOffsetBhabha'");
73 allObjectsFound =
false;
74 B2ERROR(
"No valid DBObject found for 'ECLCrystalTimeOffsetBhabhaGamma'");
77 allObjectsFound =
false;
78 B2ERROR(
"No valid DBObject found for 'ECLCrystalTimeOffsetCosmic'");
83 allObjectsFound =
false;
84 B2ERROR(
"No valid DBObject found for 'ECLReferenceCrystalPerCrateCalib'");
89 allObjectsFound =
false;
90 B2ERROR(
"No valid DBObject found for 'ECLCrateTimeOffset'");
94 if (allObjectsFound) {
95 B2INFO(
"Valid objects found for 'ECLCrystalTimeOffsetBhabha', 'ECLCrystalTimeOffsetBhabhaGamma' and 'ECLCrystalTimeOffsetCosmic')");
96 B2INFO(
"Valid object found for 'ECLReferenceCrystalPerCrateCalib'");
97 B2INFO(
"Valid object found for 'ECLCrateTimeOffset'");
99 B2INFO(
"Exiting with failure");
111 bool minRunNumBool =
false;
112 bool maxRunNumBool =
false;
118 int expNumber = expRun.first;
119 int runNumber = expRun.second;
120 if (!minRunNumBool) {
121 minExpNum = expNumber;
122 minRunNum = runNumber;
123 minRunNumBool =
true;
125 if (!maxRunNumBool) {
126 maxExpNum = expNumber;
127 maxRunNum = runNumber;
128 maxRunNumBool =
true;
130 if (((minRunNum > runNumber) && (minExpNum >= expNumber)) ||
131 (minExpNum > expNumber)) {
132 minExpNum = expNumber;
133 minRunNum = runNumber;
135 if (((maxRunNum < runNumber) && (maxExpNum <= expNumber)) ||
136 (maxExpNum < expNumber))
139 maxExpNum = expNumber;
140 maxRunNum = runNumber;
144 string debugFilenameBase =
"ECLCrystalTimeOffsetMerged";
145 B2INFO(
"debugFilenameBase = " << debugFilenameBase);
146 string runNumsString = string(
"_") + to_string(minExpNum) +
"_" + to_string(minRunNum) + string(
"-") +
147 to_string(maxExpNum) +
"_" + to_string(maxRunNum);
148 string debugFilename = debugFilenameBase + runNumsString + string(
".root");
164 B2INFO(
"Loaded 'ECLCrystalTimeOffset*' calibrations");
169 B2INFO(
"Loaded 'ECLReferenceCrystalPerCrateCalib' calibration");
173 B2INFO(
"Loaded 'ECLCrateTimeOffset' calibration");
178 vector<float> prevValuesCrys(8736);
179 vector<float> prevValuesCrysUnc(8736);
182 vector<float> prevBhabhaValuesCrys(8736);
183 vector<float> prevBhabhaValuesCrysUnc(8736);
187 prevValuesCrys = customPrevCrystalTimeObject->getCalibVector();
188 prevValuesCrysUnc = customPrevCrystalTimeObject->getCalibUncVector();
190 prevBhabhaValuesCrys = customPrevBhabhaCrystalTimeObject->getCalibVector();
191 prevBhabhaValuesCrysUnc = customPrevBhabhaCrystalTimeObject->getCalibUncVector();
194 B2INFO(
"Previous values read from database. Write out for their values for comparison");
195 for (
int ic = 0; ic < 8736; ic += 500) {
196 B2INFO(
"ts custom previous payload: cellID " << ic + 1 <<
" " << prevValuesCrys[ic]);
198 B2INFO(
"Previous bhabha values read from database. Write out for their values for comparison");
199 for (
int ic = 0; ic < 8736; ic += 500) {
200 B2INFO(
"ts custom previous bhabha payload: cellID " << ic + 1 <<
" " << prevBhabhaValuesCrys[ic]);
238 vector<bool> bhabhaGammaCalibGoodQuality(
m_numCrystals,
false);
247 if (bhabhaCalibUnc[ic] != 0
248 && fabs(bhabhaCalibUnc[ic])*TICKS_TO_NS < 2) {
255 bhabhaCalibGoodQuality[ic] = true ;
266 if (bhabhaGammaCalibUnc[ic] != 0
267 && fabs(bhabhaGammaCalibUnc[ic])*TICKS_TO_NS <
269 bhabhaGammaCalibGoodQuality[ic] = true ;
279 vector<float> cosmicCalibShifted = cosmicCalib;
280 vector<float> cosmicCalibShiftedUnc = cosmicCalibUnc;
282 cosmicCalibShifted[ic] -= crateCalib[ic];
283 int crate_id_from_crystal = crystalMapper->
getCrateID(ic + 1);
284 crateCalibShort[crate_id_from_crystal - 1] = crateCalib[ic];
287 for (
int icrate = 0; icrate <
m_numCrates; icrate++) {
288 B2INFO(
"crate time calibration for crate " <<
289 icrate + 1 <<
" = " << crateCalibShort[icrate]);
296 vector<float> meanGoodCrystalCalibPerCrate_bhabha(
m_numCrates, 0);
297 vector<float> meanGoodCrystalCalibPerCrate_cosmic(
m_numCrates, 0);
304 vector<short> numGoodCrystalCalibPerCrate_bhabha(
m_numCrates, 1);
309 for (
int ic_base1 = 1; ic_base1 <=
m_numCrystals; ic_base1++) {
310 int crate_id_from_crystal = crystalMapper->
getCrateID(ic_base1);
312 if (!((bhabhaCalib[ic_base1 - 1] == 0) || (! bhabhaCalibGoodQuality[ic_base1 - 1]))) {
313 meanGoodCrystalCalibPerCrate_bhabha[crate_id_from_crystal - 1] += bhabhaCalib[ic_base1 - 1];
314 meanGoodCrystalCalibPerCrate_cosmic[crate_id_from_crystal - 1] += cosmicCalibShifted[ic_base1 - 1];
315 numGoodCrystalCalibPerCrate_bhabha[crate_id_from_crystal - 1]++ ;
320 for (
int icrate = 0; icrate <
m_numCrates; icrate++) {
321 meanGoodCrystalCalibPerCrate_bhabha[icrate] /= numGoodCrystalCalibPerCrate_bhabha[icrate] ;
322 meanGoodCrystalCalibPerCrate_cosmic[icrate] /= numGoodCrystalCalibPerCrate_bhabha[icrate] ;
323 B2INFO(
"Mean bhabha crystal calibration for crate " << icrate + 1 <<
" = " << meanGoodCrystalCalibPerCrate_bhabha[icrate]);
324 B2INFO(
"Mean cosmic crystal calibration for crate " << icrate + 1 <<
" = " << meanGoodCrystalCalibPerCrate_cosmic[icrate]);
325 B2INFO(
"Number of good crystal calibrations for crate" << icrate + 1 <<
" = " << numGoodCrystalCalibPerCrate_bhabha[icrate]);
333 for (
int ic_base1 = 1; ic_base1 <=
m_numCrystals; ic_base1++) {
334 int crate_id_from_crystal = crystalMapper->
getCrateID(ic_base1);
335 cosmicCalibShifted[ic_base1 - 1] += meanGoodCrystalCalibPerCrate_bhabha[crate_id_from_crystal - 1] -
336 meanGoodCrystalCalibPerCrate_cosmic[crate_id_from_crystal - 1] ;
353 vector<float> newBhabhaMinusCustomPrevCalib__cid(
m_numCrystals);
354 vector<float> newBhabhaMinusCustomPrevBhabhaCalib__cid(
m_numCrystals);
355 vector<float> newBhabhaMinusCustomPrevCalib;
356 vector<float> newBhabhaMinusCustomPrevBhabhaCalib;
361 B2DEBUG(29,
"Merging crystal " << ic + 1);
362 int crate_id_from_crystal = crystalMapper->
getCrateID(ic + 1);
365 bool isRefCrys = false ;
366 for (
int icrate = 0; icrate <
m_numCrates; icrate++) {
367 if (ic + 1 == refCrystalPerCrate[icrate]) {
369 B2INFO(
"Crystal " << ic + 1 <<
" is a reference crystal");
372 B2DEBUG(29,
"Checked reference crystal for crate " << icrate + 1);
374 B2DEBUG(29,
"Checked ref crystal");
376 string whichCalibUsed =
"bhabha";
379 newCalib[ic] = bhabhaCalib[ic];
380 newCalibUnc[ic] = fabs(bhabhaCalibUnc[ic]);
382 double newMinusMerged ;
383 double newMinusBhabha ;
391 ((bhabhaCalib[ic] == 0) || (! bhabhaCalibGoodQuality[ic]))) {
395 if (bhabhaGammaCalib[ic] != 0 || bhabhaGammaCalibGoodQuality[ic]) {
396 newCalib[ic] = bhabhaGammaCalib[ic] ;
397 newCalibUnc[ic] = fabs(bhabhaGammaCalibUnc[ic]);
398 whichCalibUsed =
"bhabhaGamma";
400 newCalib[ic] = cosmicCalibShifted[ic] ;
401 newCalibUnc[ic] = cosmicCalibShiftedUnc[ic];
402 whichCalibUsed =
"cosmic";
405 newMinusMerged = (newCalib[ic] - prevValuesCrys[ic]) * TICKS_TO_NS;
406 newMinusBhabha = (newCalib[ic] - prevBhabhaValuesCrys[ic]) * TICKS_TO_NS;
408 newBhabhaMinusCustomPrevCalib__cid[ic] = newMinusMerged;
409 newBhabhaMinusCustomPrevBhabhaCalib__cid[ic] = newMinusBhabha;
410 newBhabhaMinusCustomPrevCalib.push_back(newMinusMerged);
411 newBhabhaMinusCustomPrevBhabhaCalib.push_back(newMinusBhabha);
414 newMinusMerged = (newCalib[ic] - prevValuesCrys[ic]) * TICKS_TO_NS;
415 newMinusBhabha = (newCalib[ic] - prevBhabhaValuesCrys[ic]) * TICKS_TO_NS;
418 B2INFO(
"Crystal " << ic + 1 <<
", crate " << crate_id_from_crystal <<
419 "::: bhabha = " << bhabhaCalib[ic] <<
"+-" << bhabhaCalibUnc[ic] <<
420 " ticks, rad bhabha = " << bhabhaGammaCalib[ic] <<
"+-" << bhabhaGammaCalibUnc[ic] <<
421 " ticks, cosmic UNshifted = " << cosmicCalib[ic] <<
"+-" << cosmicCalibUnc[ic] <<
422 " ticks, cosmic shifted = " << cosmicCalibShifted[ic] <<
"+-" << cosmicCalibShiftedUnc[ic] <<
423 " ticks. ||| New calib = " << newCalib[ic] <<
" ticks, selected " << whichCalibUsed);
424 B2INFO(
" Crystal " << ic + 1 <<
", pre calib ts = " << prevValuesCrys[ic] <<
425 ", pre calib bhabha ts = " << prevBhabhaValuesCrys[ic] <<
426 " ticks ||| new - pre calib ts = " << newMinusMerged <<
427 ", new - pre calib bhabha ts = " << newMinusBhabha <<
" ns");
436 TDirectory* executeDir = gDirectory;
438 TString fname = debugFilename;
439 TFile hfile(fname,
"recreate");
441 TString htitle =
"ECLCrystalTimeOffsetBhabha;cellID;Bhabha ts values [ticks]";
442 TH1F* bhabhaPayload =
new TH1F(
"bhabhaPayload", htitle,
m_numCrystals, 1, 8737);
444 htitle =
"ECLCrystalTimeOffsetBhabhaGamma;cellID;Radiative bhabha ts values [ticks]";
445 TH1F* bhabhaGammaPayload =
new TH1F(
"bhabhaGammaPayload", htitle,
m_numCrystals, 1, 8737);
447 htitle =
"ECLCrystalTimeOffsetCosmic : unshifted;cellID;Unshifted cosmic ts values [ticks]";
448 TH1F* cosmicPayload =
new TH1F(
"cosmicUnshiftedPayload", htitle,
m_numCrystals, 1, 8737);
450 htitle =
"ECLCrystalTimeOffsetCosmic : shifted;cellID;Shifted cosmic ts values [ticks]";
451 TH1F* cosmicShiftedPayload =
new TH1F(
"cosmicShiftedPayload", htitle,
m_numCrystals, 1, 8737);
453 htitle =
"New ECLCrystalTimeOffset;cellID;New (merged) ts values [ticks]";
454 TH1F* newPayload =
new TH1F(
"newPayload", htitle,
m_numCrystals, 1, 8737);
457 htitle =
"ECLCrystalTimeOffsetBhabha : 'pre-calib' values;cellID;Pre-calibration bhabha ts values [ticks]";
458 TH1F* tsCustomPrevBhabha_payload =
new TH1F(
"tsCustomPrevBhabha_payload", htitle,
m_numCrystals, 1, 8737);
460 htitle =
"ECLCrystalTimeOffset : 'pre-calib' values;cellID;Pre-calibration merged ts values [ticks]";
461 TH1F* tsCustomPrev_payload =
new TH1F(
"tsCustomPrev_payload", htitle,
m_numCrystals, 1, 8737);
465 htitle =
"-Only for crystals where the new bhabha ts value is used-;cellID;ts(new | bhabha) - ts(old = 'pre-calib' | merged) [ns]";
466 TH1F* newBhabhaMinusCustomPrev__cid =
new TH1F(
"newBhabhaMinusCustomPrev__cid", htitle,
m_numCrystals, 1, 8737);
468 htitle =
"-Only for crystals where the new bhabha ts value is used-;cellID;ts(new | bhabha) - ts(old = 'pre-calib' | bhabha) [ns]";
469 TH1F* newBhabhaMinusCustomPrevBhabha__cid =
new TH1F(
"newBhabhaMinusCustomPrevBhabha__cid", htitle,
m_numCrystals, 1, 8737);
471 htitle =
"-Only for crystals where the new bhabha ts value is used-;ts(new | bhabha) - ts(old = 'pre-calib' | merged) [ns];Number of crystals";
472 auto tsNewBhabha_MINUS_tsCustomPrev =
new TH1F(
"TsNewBhabha_MINUS_TsCustomPrev", htitle, 285, -69.5801, 69.5801);
474 htitle =
"-Only for crystals where the new bhabha ts value is used-;ts(new | bhabha) - ts(old = 'pre-calib' | bhabha) [ns];Number of crystals";
475 auto tsNewBhabha_MINUS_tsCustomPrevBhabha =
new TH1F(
"TsNewBhabha_MINUS_TsCustomPrevBhabha", htitle, 285, -69.5801, 69.5801);
479 bhabhaPayload->SetBinContent(cellID, bhabhaCalib[cellID - 1]);
480 bhabhaPayload->SetBinError(cellID, fabs(bhabhaCalibUnc[cellID - 1]));
482 bhabhaGammaPayload->SetBinContent(cellID, bhabhaGammaCalib[cellID - 1]);
483 bhabhaGammaPayload->SetBinError(cellID, fabs(bhabhaGammaCalibUnc[cellID - 1]));
485 cosmicPayload->SetBinContent(cellID, cosmicCalib[cellID - 1]);
486 cosmicPayload->SetBinError(cellID, cosmicCalibUnc[cellID - 1]);
488 cosmicShiftedPayload->SetBinContent(cellID, cosmicCalibShifted[cellID - 1]);
489 cosmicShiftedPayload->SetBinError(cellID, cosmicCalibShiftedUnc[cellID - 1]);
491 newPayload->SetBinContent(cellID, newCalib[cellID - 1]);
492 newPayload->SetBinError(cellID, newCalibUnc[cellID - 1]);
496 newBhabhaMinusCustomPrev__cid->SetBinContent(cellID, newBhabhaMinusCustomPrevCalib__cid[cellID - 1]);
497 newBhabhaMinusCustomPrev__cid->SetBinError(cellID, 0);
499 newBhabhaMinusCustomPrevBhabha__cid->SetBinContent(cellID, newBhabhaMinusCustomPrevBhabhaCalib__cid[cellID - 1]);
500 newBhabhaMinusCustomPrevBhabha__cid->SetBinError(cellID, 0);
503 tsCustomPrev_payload->SetBinContent(cellID, prevValuesCrys[cellID - 1]);
504 tsCustomPrev_payload->SetBinError(cellID, prevValuesCrysUnc[cellID - 1]);
506 tsCustomPrevBhabha_payload->SetBinContent(cellID, prevBhabhaValuesCrys[cellID - 1]);
507 tsCustomPrevBhabha_payload->SetBinError(cellID, fabs(prevBhabhaValuesCrysUnc[cellID - 1]));
512 for (
long unsigned int ib = 0; ib < newBhabhaMinusCustomPrevCalib.size(); ib++) {
513 tsNewBhabha_MINUS_tsCustomPrev->Fill(newBhabhaMinusCustomPrevCalib[ib]);
516 for (
long unsigned int ib = 0; ib < newBhabhaMinusCustomPrevBhabhaCalib.size(); ib++) {
517 tsNewBhabha_MINUS_tsCustomPrevBhabha->Fill(newBhabhaMinusCustomPrevBhabhaCalib[ib]);
525 B2INFO(
"Debugging histograms written to " << fname);