20 #include <ecl/modules/eclDigitCalibration/ECLDigitCalibratorModule.h>
23 #include <unordered_map>
30 #include <ecl/dataobjects/ECLDigit.h>
31 #include <ecl/dataobjects/ECLCalDigit.h>
32 #include <ecl/digitization/EclConfiguration.h>
33 #include <ecl/dataobjects/ECLPureCsIInfo.h>
34 #include <ecl/dataobjects/ECLDsp.h>
35 #include <ecl/utility/utilityFunctions.h>
36 #include <ecl/geometry/ECLGeometryPar.h>
37 #include <ecl/dbobjects/ECLCrystalCalib.h>
40 #include <framework/gearbox/Unit.h>
41 #include <framework/logging/Logger.h>
42 #include <framework/utilities/FileSystem.h>
43 #include <framework/geometry/B2Vector3.h>
44 #include <framework/core/Environment.h>
47 #include <mdst/dataobjects/EventLevelClusteringInfo.h>
65 m_calibrationCrystalElectronics("ECLCrystalElectronics"),
66 m_calibrationCrystalEnergy("ECLCrystalEnergy"),
67 m_calibrationCrystalElectronicsTime("ECLCrystalElectronicsTime"),
68 m_calibrationCrystalTimeOffset("ECLCrystalTimeOffset"),
69 m_calibrationCrateTimeOffset("ECLCrateTimeOffset"),
70 m_calibrationCrystalFlightTime("ECLCrystalFlightTime"),
71 m_eclDigits(eclDigitArrayName()),
72 m_eclCalDigits(eclCalDigitArrayName()),
73 m_eclDsps(eclDspArrayName()),
74 m_eclPureCsIInfo(eclPureCsIInfoArrayName())
77 setDescription(
"Applies digit energy, time and time-resolution calibration to each ECL digit. Counts number of out-of-time background digits to determine the event-by-event background level.");
78 addParam(
"backgroundEnergyCut", m_backgroundEnergyCut,
"Energy cut used to count background digits", 7.0 *
Belle2::Unit::MeV);
79 addParam(
"backgroundTimingCut", m_backgroundTimingCut,
"Timing cut used to count background digits", 110.0 *
Belle2::Unit::ns);
80 addParam(
"fileBackgroundName", m_fileBackgroundName,
"Background filename.",
81 FileSystem::findFile(
"/data/ecl/background_norm.root"));
82 addParam(
"simulatePure", m_simulatePure,
"Flag to simulate pure CsI option",
false);
85 addParam(
"energyDependenceTimeOffsetFitParam_p1", m_energyDependenceTimeOffsetFitParam_p1,
86 "Fit parameter (p1) for applying correction to the time offset as a function of the energy (amplitude)", -999.0);
87 addParam(
"energyDependenceTimeOffsetFitParam_p2", m_energyDependenceTimeOffsetFitParam_p2,
88 "Fit parameter (p2) for applying correction to the time offset as a function of the energy (amplitude)", -999.0);
89 addParam(
"energyDependenceTimeOffsetFitParam_p3", m_energyDependenceTimeOffsetFitParam_p3,
90 "Fit parameter (p3) for applying correction to the time offset as a function of the energy (amplitude)", -999.0);
91 addParam(
"energyDependenceTimeOffsetFitParam_p4", m_energyDependenceTimeOffsetFitParam_p4,
92 "Fit parameter (p4) for applying correction to the time offset as a function of the energy (amplitude)", -999.0);
93 addParam(
"energyDependenceTimeOffsetFitParam_p5", m_energyDependenceTimeOffsetFitParam_p5,
94 "Fit parameter (p5) for applying correction to the time offset as a function of the energy (amplitude)", -999.0);
95 addParam(
"energyDependenceTimeOffsetFitParam_p6", m_energyDependenceTimeOffsetFitParam_p6,
96 "Fit parameter (p6) for applying correction to the time offset as a function of the energy (amplitude)", -999.0);
100 setPropertyFlags(c_ParallelProcessingCertified);
104 m_timeInverseSlope = 0.0;
109 ECLDigitCalibratorModule::~ECLDigitCalibratorModule()
114 void ECLDigitCalibratorModule::initializeCalibration()
117 m_timeInverseSlope = 1.0 / (4.0 * EclConfiguration::m_rf) *
119 m_pureCsIEnergyCalib = 0.00005;
120 m_pureCsITimeCalib = 0.1;
121 m_pureCsITimeOffset = 0.31;
123 callbackCalibration(m_calibrationCrystalElectronics, v_calibrationCrystalElectronics, v_calibrationCrystalElectronicsUnc);
124 callbackCalibration(m_calibrationCrystalEnergy, v_calibrationCrystalEnergy, v_calibrationCrystalEnergyUnc);
125 callbackCalibration(m_calibrationCrystalElectronicsTime, v_calibrationCrystalElectronicsTime,
126 v_calibrationCrystalElectronicsTimeUnc);
127 callbackCalibration(m_calibrationCrystalTimeOffset, v_calibrationCrystalTimeOffset, v_calibrationCrystalTimeOffsetUnc);
128 callbackCalibration(m_calibrationCrateTimeOffset, v_calibrationCrateTimeOffset, v_calibrationCrateTimeOffsetUnc);
129 callbackCalibration(m_calibrationCrystalFlightTime, v_calibrationCrystalFlightTime, v_calibrationCrystalFlightTimeUnc);
134 std::vector<float>& constantsUnc)
136 constants = cal->getCalibVector();
137 constantsUnc = cal->getCalibUncVector();
142 void ECLDigitCalibratorModule::initialize()
145 m_eventLevelClusteringInfo.registerInDataStore(eventLevelClusteringInfoName());
148 m_eclDigits.registerInDataStore(eclDigitArrayName());
149 m_eclDsps.registerInDataStore(eclDspArrayName());
150 m_eclCalDigits.registerInDataStore(eclCalDigitArrayName());
151 m_eclCalDigits.registerRelationTo(m_eclDigits);
154 m_eclPureCsIInfo.registerInDataStore(eclPureCsIInfoArrayName());
157 initializeCalibration();
161 m_fileBackground =
new TFile(m_fileBackgroundName.c_str(),
"READ");
162 if (!m_fileBackground) B2FATAL(
"Could not find file: " << m_fileBackgroundName);
163 m_th1fBackground =
dynamic_cast<TH1F*
>(m_fileBackground->Get(
"background"));
164 if (!m_th1fBackground) B2FATAL(
"Could not find m_th1fBackground");
167 m_averageBG = m_th1fBackground->Integral() / m_th1fBackground->GetEntries();
170 if (fabs(c_pol2Var3) > 1e-12) {
171 m_pol2Max = -c_pol2Var2 / (2 * c_pol2Var3);
176 if ((m_energyDependenceTimeOffsetFitParam_p1 != -999) &&
177 (m_energyDependenceTimeOffsetFitParam_p2 != -999) &&
178 (m_energyDependenceTimeOffsetFitParam_p3 != -999) &&
179 (m_energyDependenceTimeOffsetFitParam_p4 != -999) &&
180 (m_energyDependenceTimeOffsetFitParam_p5 != -999) &&
181 (m_energyDependenceTimeOffsetFitParam_p6 != -999)) {
182 B2DEBUG(80,
"m_energyDependenceTimeOffsetFitParam_p1 = " << m_energyDependenceTimeOffsetFitParam_p1);
183 B2DEBUG(80,
"m_energyDependenceTimeOffsetFitParam_p2 = " << m_energyDependenceTimeOffsetFitParam_p2);
184 B2DEBUG(80,
"m_energyDependenceTimeOffsetFitParam_p3 = " << m_energyDependenceTimeOffsetFitParam_p3);
185 B2DEBUG(80,
"m_energyDependenceTimeOffsetFitParam_p4 = " << m_energyDependenceTimeOffsetFitParam_p4);
186 B2DEBUG(80,
"m_energyDependenceTimeOffsetFitParam_p5 = " << m_energyDependenceTimeOffsetFitParam_p5);
187 B2DEBUG(80,
"m_energyDependenceTimeOffsetFitParam_p6 = " << m_energyDependenceTimeOffsetFitParam_p6);
189 ECLTimeUtil->setTimeWalkFuncParams(m_energyDependenceTimeOffsetFitParam_p1,
190 m_energyDependenceTimeOffsetFitParam_p2,
191 m_energyDependenceTimeOffsetFitParam_p3,
192 m_energyDependenceTimeOffsetFitParam_p4,
193 m_energyDependenceTimeOffsetFitParam_p5,
194 m_energyDependenceTimeOffsetFitParam_p6) ;
199 void ECLDigitCalibratorModule::beginRun()
202 if (m_calibrationCrystalElectronics.hasChanged()) {
203 if (m_calibrationCrystalElectronics) {
204 callbackCalibration(m_calibrationCrystalElectronics, v_calibrationCrystalElectronics, v_calibrationCrystalElectronicsUnc);
205 }
else B2ERROR(
"ECLDigitCalibratorModule::beginRun - Couldn't find m_calibrationCrystalElectronics for current run!");
208 if (m_calibrationCrystalEnergy.hasChanged()) {
209 if (m_calibrationCrystalEnergy) {
210 callbackCalibration(m_calibrationCrystalEnergy, v_calibrationCrystalEnergy, v_calibrationCrystalEnergyUnc);
211 }
else B2ERROR(
"ECLDigitCalibratorModule::beginRun - Couldn't find m_calibrationCrystalEnergy for current run!");
214 if (m_calibrationCrystalElectronicsTime.hasChanged()) {
215 if (m_calibrationCrystalElectronicsTime) {
216 callbackCalibration(m_calibrationCrystalElectronicsTime, v_calibrationCrystalElectronicsTime,
217 v_calibrationCrystalElectronicsTimeUnc);
218 }
else B2ERROR(
"ECLDigitCalibratorModule::beginRun - Couldn't find m_calibrationCrystalElectronicsTime for current run!");
221 if (m_calibrationCrystalTimeOffset.hasChanged()) {
222 if (m_calibrationCrystalTimeOffset) {
223 callbackCalibration(m_calibrationCrystalTimeOffset, v_calibrationCrystalTimeOffset, v_calibrationCrystalTimeOffsetUnc);
224 }
else B2ERROR(
"ECLDigitCalibratorModule::beginRun - Couldn't find m_calibrationCrystalTimeOffset for current run!");
227 if (m_calibrationCrateTimeOffset.hasChanged()) {
228 if (m_calibrationCrateTimeOffset) {
229 callbackCalibration(m_calibrationCrateTimeOffset, v_calibrationCrateTimeOffset, v_calibrationCrateTimeOffsetUnc);
230 }
else B2ERROR(
"ECLDigitCalibratorModule::beginRun - Couldn't find m_calibrationCrateTimeOffset for current run!");
233 if (m_calibrationCrystalFlightTime.hasChanged()) {
234 if (m_calibrationCrystalFlightTime) {
235 callbackCalibration(m_calibrationCrystalFlightTime, v_calibrationCrystalFlightTime, v_calibrationCrystalFlightTimeUnc);
236 }
else B2ERROR(
"ECLDigitCalibratorModule::beginRun - Couldn't find m_calibrationCrystalFlightTime for current run!");
242 void ECLDigitCalibratorModule::event()
246 for (
auto& aECLDigit : m_eclDigits) {
248 bool is_pure_csi = 0;
251 const auto aECLCalDigit = m_eclCalDigits.appendNew();
254 const int cellid = aECLDigit.getCellId();
257 if (cellid < 1 or cellid > c_nCrystals) {
258 B2FATAL(
"ECLDigitCalibrationModule::event():" << cellid <<
" out of range!");
262 const int amplitude = aECLDigit.getAmp();
263 double calibratedEnergy = 0;
265 if (m_simulatePure) {
266 if (aECLDigit.getRelated<
ECLPureCsIInfo>(eclPureCsIInfoArrayName()) != NULL) {
267 if (aECLDigit.getRelated<
ECLPureCsIInfo>(eclPureCsIInfoArrayName())->getPureCsI())
273 calibratedEnergy = amplitude * m_pureCsIEnergyCalib;
275 calibratedEnergy = amplitude * v_calibrationCrystalElectronics[cellid - 1] * v_calibrationCrystalEnergy[cellid - 1];
277 if (calibratedEnergy < 0.0) calibratedEnergy = 0.0;
280 const int time = aECLDigit.getTimeFit();
281 const int quality = aECLDigit.getQuality();
282 double calibratedTime = c_timeForFitFailed;
284 aECLCalDigit->addStatus(ECLCalDigit::c_IsFailedFit);
287 calibratedTime = m_pureCsITimeCalib * m_timeInverseSlope * (time - v_calibrationCrystalElectronicsTime[cellid - 1] -
288 v_calibrationCrystalTimeOffset[cellid - 1] -
289 v_calibrationCrateTimeOffset[cellid - 1])
290 - v_calibrationCrystalFlightTime[cellid - 1] + m_pureCsITimeOffset ;
292 calibratedTime = m_timeInverseSlope * (time - v_calibrationCrystalElectronicsTime[cellid - 1] -
293 v_calibrationCrystalTimeOffset[cellid - 1] -
294 v_calibrationCrateTimeOffset[cellid - 1])
295 - v_calibrationCrystalFlightTime[cellid - 1] ;
300 bool m_IsMCFlag = Environment::Instance().isMC();
301 B2DEBUG(35,
"cellid = " << cellid <<
", m_IsMCFlag = " << m_IsMCFlag) ;
304 double energyTimeShift = ECLTimeUtil->energyDependentTimeOffsetElectronic(amplitude * v_calibrationCrystalElectronics[cellid - 1]) *
306 B2DEBUG(35,
"cellid = " << cellid <<
", amplitude = " << amplitude <<
", corrected amplitude = " << amplitude *
307 v_calibrationCrystalElectronics[cellid - 1] <<
", time before t(E) shift = " << calibratedTime <<
", t(E) shift = " <<
308 energyTimeShift <<
" ns") ;
309 calibratedTime -= energyTimeShift ;
313 B2DEBUG(35,
"cellid = " << cellid <<
", amplitude = " << amplitude <<
", calibrated energy = " << calibratedEnergy);
314 B2DEBUG(35,
"cellid = " << cellid <<
", time = " << time <<
", calibratedTime = " << calibratedTime);
319 aECLCalDigit->setTwoComponentSavedChi2(ECLDsp::photonHadron, -1);
320 aECLCalDigit->setTwoComponentSavedChi2(ECLDsp::photonHadronBackgroundPhoton, -1);
321 aECLCalDigit->setTwoComponentSavedChi2(ECLDsp::photonDiodeCrossing, -1);
322 aECLCalDigit->setTwoComponentTotalEnergy(-1);
323 aECLCalDigit->setTwoComponentHadronEnergy(-1);
324 aECLCalDigit->setTwoComponentDiodeEnergy(-1);
330 const double calibratedTwoComponentTotalEnergy = aECLDsp->
getTwoComponentTotalAmp() * v_calibrationCrystalElectronics[cellid - 1] *
331 v_calibrationCrystalEnergy[cellid - 1];
332 const double calibratedTwoComponentHadronEnergy = aECLDsp->
getTwoComponentHadronAmp() * v_calibrationCrystalElectronics[cellid -
334 v_calibrationCrystalEnergy[cellid - 1];
335 const double calibratedTwoComponentDiodeEnergy = aECLDsp->
getTwoComponentDiodeAmp() * v_calibrationCrystalElectronics[cellid - 1] *
336 v_calibrationCrystalEnergy[cellid - 1];
340 aECLCalDigit->setTwoComponentTotalEnergy(calibratedTwoComponentTotalEnergy);
341 aECLCalDigit->setTwoComponentHadronEnergy(calibratedTwoComponentHadronEnergy);
342 aECLCalDigit->setTwoComponentDiodeEnergy(calibratedTwoComponentDiodeEnergy);
343 aECLCalDigit->setTwoComponentChi2(twoComponentChi2);
344 aECLCalDigit->setTwoComponentSavedChi2(ECLDsp::photonHadron, aECLDsp->
getTwoComponentSavedChi2(ECLDsp::photonHadron));
345 aECLCalDigit->setTwoComponentSavedChi2(ECLDsp::photonHadronBackgroundPhoton,
347 aECLCalDigit->setTwoComponentSavedChi2(ECLDsp::photonDiodeCrossing, aECLDsp->
getTwoComponentSavedChi2(ECLDsp::photonDiodeCrossing));
348 aECLCalDigit->setTwoComponentFitType(twoComponentFitType);
354 aECLCalDigit->setCellId(cellid);
356 aECLCalDigit->setEnergy(calibratedEnergy);
357 aECLCalDigit->addStatus(ECLCalDigit::c_IsEnergyCalibrated);
359 aECLCalDigit->setTime(calibratedTime);
360 aECLCalDigit->addStatus(ECLCalDigit::c_IsTimeCalibrated);
363 aECLCalDigit->addRelationTo(&aECLDigit);
367 const int bgCount = determineBackgroundECL();
370 for (
auto& aECLCalDigit : m_eclCalDigits) {
373 const double t99 = getT99(aECLCalDigit.getCellId(),
374 aECLCalDigit.getEnergy(),
375 aECLCalDigit.hasStatus(ECLCalDigit::c_IsFailedFit),
377 aECLCalDigit.setTimeResolution(t99);
379 if (t99 == c_timeResolutionForFitFailed) {
380 aECLCalDigit.addStatus(ECLCalDigit::c_IsFailedTimeResolution);
383 aECLCalDigit.addStatus(ECLCalDigit::c_IsTimeResolutionCalibrated);
388 void ECLDigitCalibratorModule::endRun()
394 void ECLDigitCalibratorModule::terminate()
400 double ECLDigitCalibratorModule::getT99(
const int cellid,
const double energy,
const bool fitfailed,
const int bgcount)
const
404 if (fitfailed)
return c_timeResolutionForFitFailed;
407 const double bglevel = TMath::Min((
double) bgcount / (
double) c_nominalBG * m_th1fBackground->GetBinContent(cellid) / m_averageBG,
411 const double p1 = c_pol2Var1 + c_pol2Var2 * bglevel + c_pol2Var3 * bglevel * bglevel;
418 double t99 = p1 * einv;
421 if (t99 < c_minT99) t99 = c_minT99;
423 B2DEBUG(35,
"ECLDigitCalibratorModule::getCalibratedTimeResolution: dose = " << m_th1fBackground->GetBinContent(
424 cellid) <<
", bglevel = " << bglevel <<
", cellid = " << cellid <<
", t99 = " << t99 <<
", energy = " << energy /
431 int ECLDigitCalibratorModule::determineBackgroundECL()
434 using regionCounter = std::unordered_map<ECL::DetectorRegion, uint>;
436 regionCounter outOfTimeCount {{ECL::DetectorRegion::FWD, 0},
437 {ECL::DetectorRegion::BRL, 0},
438 {ECL::DetectorRegion::BWD, 0}};
443 for (
auto& aECLCalDigit : m_eclCalDigits) {
444 if (abs(aECLCalDigit.getTime()) >= m_backgroundTimingCut) {
445 if (aECLCalDigit.getEnergy() >= m_backgroundEnergyCut) {
447 const B2Vector3D position = geom->GetCrystalPos(aECLCalDigit.getCellId() - 1);
448 const double theta = position.
Theta();
451 const auto detectorRegion = ECL::getDetectorRegion(theta);
454 ++outOfTimeCount.at(detectorRegion);
460 if (!m_eventLevelClusteringInfo) {
461 m_eventLevelClusteringInfo.create();
464 m_eventLevelClusteringInfo->setNECLCalDigitsOutOfTimeFWD(outOfTimeCount.at(ECL::DetectorRegion::FWD));
465 m_eventLevelClusteringInfo->setNECLCalDigitsOutOfTimeBarrel(outOfTimeCount.at(ECL::DetectorRegion::BRL));
466 m_eventLevelClusteringInfo->setNECLCalDigitsOutOfTimeBWD(outOfTimeCount.at(ECL::DetectorRegion::BWD));
468 B2DEBUG(35,
"ECLDigitCalibratorModule::determineBackgroundECL found " << outOfTimeCount.at(ECL::DetectorRegion::FWD) <<
", " <<
469 outOfTimeCount.at(ECL::DetectorRegion::BRL) <<
", " <<
470 outOfTimeCount.at(ECL::DetectorRegion::BWD) <<
" out of time digits in FWD, BRL, BWD");
472 return m_eventLevelClusteringInfo->getNECLCalDigitsOutOfTime();