9 #include <ecl/modules/eclDigitCalibration/ECLDigitCalibratorModule.h>
12 #include <unordered_map>
19 #include <ecl/dataobjects/ECLDigit.h>
20 #include <ecl/dataobjects/ECLCalDigit.h>
21 #include <ecl/digitization/EclConfiguration.h>
22 #include <ecl/dataobjects/ECLPureCsIInfo.h>
23 #include <ecl/dataobjects/ECLDsp.h>
24 #include <ecl/utility/utilityFunctions.h>
25 #include <ecl/geometry/ECLGeometryPar.h>
26 #include <ecl/dbobjects/ECLCrystalCalib.h>
29 #include <framework/gearbox/Unit.h>
30 #include <framework/logging/Logger.h>
31 #include <framework/utilities/FileSystem.h>
32 #include <framework/geometry/B2Vector3.h>
33 #include <framework/core/Environment.h>
36 #include <mdst/dataobjects/EventLevelClusteringInfo.h>
54 m_calibrationCrystalElectronics("ECLCrystalElectronics"),
55 m_calibrationCrystalEnergy("ECLCrystalEnergy"),
56 m_calibrationCrystalElectronicsTime("ECLCrystalElectronicsTime"),
57 m_calibrationCrystalTimeOffset("ECLCrystalTimeOffset"),
58 m_calibrationCrateTimeOffset("ECLCrateTimeOffset"),
59 m_calibrationCrystalFlightTime("ECLCrystalFlightTime"),
60 m_eclDigits(eclDigitArrayName()),
61 m_eclCalDigits(eclCalDigitArrayName()),
62 m_eclDsps(eclDspArrayName()),
63 m_eclPureCsIInfo(eclPureCsIInfoArrayName())
66 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.");
67 addParam(
"backgroundEnergyCut", m_backgroundEnergyCut,
"Energy cut used to count background digits", 7.0 *
Belle2::Unit::MeV);
68 addParam(
"backgroundTimingCut", m_backgroundTimingCut,
"Timing cut used to count background digits", 110.0 *
Belle2::Unit::ns);
69 addParam(
"fileBackgroundName", m_fileBackgroundName,
"Background filename.",
70 FileSystem::findFile(
"/data/ecl/background_norm.root"));
71 addParam(
"simulatePure", m_simulatePure,
"Flag to simulate pure CsI option",
false);
74 addParam(
"energyDependenceTimeOffsetFitParam_p1", m_energyDependenceTimeOffsetFitParam_p1,
75 "Fit parameter (p1) for applying correction to the time offset as a function of the energy (amplitude)", -999.0);
76 addParam(
"energyDependenceTimeOffsetFitParam_p2", m_energyDependenceTimeOffsetFitParam_p2,
77 "Fit parameter (p2) for applying correction to the time offset as a function of the energy (amplitude)", -999.0);
78 addParam(
"energyDependenceTimeOffsetFitParam_p3", m_energyDependenceTimeOffsetFitParam_p3,
79 "Fit parameter (p3) for applying correction to the time offset as a function of the energy (amplitude)", -999.0);
80 addParam(
"energyDependenceTimeOffsetFitParam_p4", m_energyDependenceTimeOffsetFitParam_p4,
81 "Fit parameter (p4) for applying correction to the time offset as a function of the energy (amplitude)", -999.0);
82 addParam(
"energyDependenceTimeOffsetFitParam_p5", m_energyDependenceTimeOffsetFitParam_p5,
83 "Fit parameter (p5) for applying correction to the time offset as a function of the energy (amplitude)", -999.0);
84 addParam(
"energyDependenceTimeOffsetFitParam_p6", m_energyDependenceTimeOffsetFitParam_p6,
85 "Fit parameter (p6) for applying correction to the time offset as a function of the energy (amplitude)", -999.0);
89 setPropertyFlags(c_ParallelProcessingCertified);
93 m_timeInverseSlope = 0.0;
98 ECLDigitCalibratorModule::~ECLDigitCalibratorModule()
103 void ECLDigitCalibratorModule::initializeCalibration()
107 m_timeInverseSlope = 1.0 / (4.0 * EclConfiguration::m_rf) * 1e3;
108 m_pureCsIEnergyCalib = 0.00005;
109 m_pureCsITimeCalib = 0.1;
110 m_pureCsITimeOffset = 0.31;
112 callbackCalibration(m_calibrationCrystalElectronics, v_calibrationCrystalElectronics, v_calibrationCrystalElectronicsUnc);
113 callbackCalibration(m_calibrationCrystalEnergy, v_calibrationCrystalEnergy, v_calibrationCrystalEnergyUnc);
114 callbackCalibration(m_calibrationCrystalElectronicsTime, v_calibrationCrystalElectronicsTime,
115 v_calibrationCrystalElectronicsTimeUnc);
116 callbackCalibration(m_calibrationCrystalTimeOffset, v_calibrationCrystalTimeOffset, v_calibrationCrystalTimeOffsetUnc);
117 callbackCalibration(m_calibrationCrateTimeOffset, v_calibrationCrateTimeOffset, v_calibrationCrateTimeOffsetUnc);
118 callbackCalibration(m_calibrationCrystalFlightTime, v_calibrationCrystalFlightTime, v_calibrationCrystalFlightTimeUnc);
123 std::vector<float>& constantsUnc)
125 constants = cal->getCalibVector();
126 constantsUnc = cal->getCalibUncVector();
131 void ECLDigitCalibratorModule::initialize()
134 m_eventLevelClusteringInfo.registerInDataStore(eventLevelClusteringInfoName());
137 m_eclDigits.registerInDataStore(eclDigitArrayName());
138 m_eclDsps.registerInDataStore(eclDspArrayName());
139 m_eclCalDigits.registerInDataStore(eclCalDigitArrayName());
140 m_eclCalDigits.registerRelationTo(m_eclDigits);
143 m_eclPureCsIInfo.registerInDataStore(eclPureCsIInfoArrayName());
146 initializeCalibration();
150 m_fileBackground =
new TFile(m_fileBackgroundName.c_str(),
"READ");
151 if (!m_fileBackground) B2FATAL(
"Could not find file: " << m_fileBackgroundName);
152 m_th1fBackground =
dynamic_cast<TH1F*
>(m_fileBackground->Get(
"background"));
153 if (!m_th1fBackground) B2FATAL(
"Could not find m_th1fBackground");
156 m_averageBG = m_th1fBackground->Integral() / m_th1fBackground->GetEntries();
159 if (fabs(c_pol2Var3) > 1e-12) {
160 m_pol2Max = -c_pol2Var2 / (2 * c_pol2Var3);
165 if ((m_energyDependenceTimeOffsetFitParam_p1 != -999) &&
166 (m_energyDependenceTimeOffsetFitParam_p2 != -999) &&
167 (m_energyDependenceTimeOffsetFitParam_p3 != -999) &&
168 (m_energyDependenceTimeOffsetFitParam_p4 != -999) &&
169 (m_energyDependenceTimeOffsetFitParam_p5 != -999) &&
170 (m_energyDependenceTimeOffsetFitParam_p6 != -999)) {
171 B2DEBUG(80,
"m_energyDependenceTimeOffsetFitParam_p1 = " << m_energyDependenceTimeOffsetFitParam_p1);
172 B2DEBUG(80,
"m_energyDependenceTimeOffsetFitParam_p2 = " << m_energyDependenceTimeOffsetFitParam_p2);
173 B2DEBUG(80,
"m_energyDependenceTimeOffsetFitParam_p3 = " << m_energyDependenceTimeOffsetFitParam_p3);
174 B2DEBUG(80,
"m_energyDependenceTimeOffsetFitParam_p4 = " << m_energyDependenceTimeOffsetFitParam_p4);
175 B2DEBUG(80,
"m_energyDependenceTimeOffsetFitParam_p5 = " << m_energyDependenceTimeOffsetFitParam_p5);
176 B2DEBUG(80,
"m_energyDependenceTimeOffsetFitParam_p6 = " << m_energyDependenceTimeOffsetFitParam_p6);
178 ECLTimeUtil->setTimeWalkFuncParams(m_energyDependenceTimeOffsetFitParam_p1,
179 m_energyDependenceTimeOffsetFitParam_p2,
180 m_energyDependenceTimeOffsetFitParam_p3,
181 m_energyDependenceTimeOffsetFitParam_p4,
182 m_energyDependenceTimeOffsetFitParam_p5,
183 m_energyDependenceTimeOffsetFitParam_p6) ;
188 void ECLDigitCalibratorModule::beginRun()
191 if (m_calibrationCrystalElectronics.hasChanged()) {
192 if (m_calibrationCrystalElectronics) {
193 callbackCalibration(m_calibrationCrystalElectronics, v_calibrationCrystalElectronics, v_calibrationCrystalElectronicsUnc);
194 }
else B2ERROR(
"ECLDigitCalibratorModule::beginRun - Couldn't find m_calibrationCrystalElectronics for current run!");
197 if (m_calibrationCrystalEnergy.hasChanged()) {
198 if (m_calibrationCrystalEnergy) {
199 callbackCalibration(m_calibrationCrystalEnergy, v_calibrationCrystalEnergy, v_calibrationCrystalEnergyUnc);
200 }
else B2ERROR(
"ECLDigitCalibratorModule::beginRun - Couldn't find m_calibrationCrystalEnergy for current run!");
203 if (m_calibrationCrystalElectronicsTime.hasChanged()) {
204 if (m_calibrationCrystalElectronicsTime) {
205 callbackCalibration(m_calibrationCrystalElectronicsTime, v_calibrationCrystalElectronicsTime,
206 v_calibrationCrystalElectronicsTimeUnc);
207 }
else B2ERROR(
"ECLDigitCalibratorModule::beginRun - Couldn't find m_calibrationCrystalElectronicsTime for current run!");
210 if (m_calibrationCrystalTimeOffset.hasChanged()) {
211 if (m_calibrationCrystalTimeOffset) {
212 callbackCalibration(m_calibrationCrystalTimeOffset, v_calibrationCrystalTimeOffset, v_calibrationCrystalTimeOffsetUnc);
213 }
else B2ERROR(
"ECLDigitCalibratorModule::beginRun - Couldn't find m_calibrationCrystalTimeOffset for current run!");
216 if (m_calibrationCrateTimeOffset.hasChanged()) {
217 if (m_calibrationCrateTimeOffset) {
218 callbackCalibration(m_calibrationCrateTimeOffset, v_calibrationCrateTimeOffset, v_calibrationCrateTimeOffsetUnc);
219 }
else B2ERROR(
"ECLDigitCalibratorModule::beginRun - Couldn't find m_calibrationCrateTimeOffset for current run!");
222 if (m_calibrationCrystalFlightTime.hasChanged()) {
223 if (m_calibrationCrystalFlightTime) {
224 callbackCalibration(m_calibrationCrystalFlightTime, v_calibrationCrystalFlightTime, v_calibrationCrystalFlightTimeUnc);
225 }
else B2ERROR(
"ECLDigitCalibratorModule::beginRun - Couldn't find m_calibrationCrystalFlightTime for current run!");
231 void ECLDigitCalibratorModule::event()
235 for (
auto& aECLDigit : m_eclDigits) {
237 bool is_pure_csi = 0;
240 const auto aECLCalDigit = m_eclCalDigits.appendNew();
243 const int cellid = aECLDigit.getCellId();
246 if (cellid < 1 or cellid > c_nCrystals) {
247 B2FATAL(
"ECLDigitCalibrationModule::event():" << cellid <<
" out of range!");
251 const int amplitude = aECLDigit.getAmp();
252 double calibratedEnergy = 0;
254 if (m_simulatePure) {
255 if (aECLDigit.getRelated<
ECLPureCsIInfo>(eclPureCsIInfoArrayName()) !=
nullptr) {
256 if (aECLDigit.getRelated<
ECLPureCsIInfo>(eclPureCsIInfoArrayName())->getPureCsI())
262 calibratedEnergy = amplitude * m_pureCsIEnergyCalib;
264 calibratedEnergy = amplitude * v_calibrationCrystalElectronics[cellid - 1] * v_calibrationCrystalEnergy[cellid - 1];
266 if (calibratedEnergy < 0.0) calibratedEnergy = 0.0;
269 const int time = aECLDigit.getTimeFit();
270 const int quality = aECLDigit.getQuality();
271 double calibratedTime = c_timeForFitFailed;
273 aECLCalDigit->addStatus(ECLCalDigit::c_IsFailedFit);
276 calibratedTime = m_pureCsITimeCalib * m_timeInverseSlope * (time - v_calibrationCrystalElectronicsTime[cellid - 1] -
277 v_calibrationCrystalTimeOffset[cellid - 1] -
278 v_calibrationCrateTimeOffset[cellid - 1])
279 - v_calibrationCrystalFlightTime[cellid - 1] + m_pureCsITimeOffset ;
281 calibratedTime = m_timeInverseSlope * (time - v_calibrationCrystalElectronicsTime[cellid - 1] -
282 v_calibrationCrystalTimeOffset[cellid - 1] -
283 v_calibrationCrateTimeOffset[cellid - 1])
284 - v_calibrationCrystalFlightTime[cellid - 1] ;
289 bool m_IsMCFlag = Environment::Instance().isMC();
290 B2DEBUG(35,
"cellid = " << cellid <<
", m_IsMCFlag = " << m_IsMCFlag) ;
293 double energyTimeShift = ECLTimeUtil->energyDependentTimeOffsetElectronic(amplitude * v_calibrationCrystalElectronics[cellid - 1]) *
295 B2DEBUG(35,
"cellid = " << cellid <<
", amplitude = " << amplitude <<
", corrected amplitude = " << amplitude *
296 v_calibrationCrystalElectronics[cellid - 1] <<
", time before t(E) shift = " << calibratedTime <<
", t(E) shift = " <<
297 energyTimeShift <<
" ns") ;
298 calibratedTime -= energyTimeShift ;
302 B2DEBUG(35,
"cellid = " << cellid <<
", amplitude = " << amplitude <<
", calibrated energy = " << calibratedEnergy);
303 B2DEBUG(35,
"cellid = " << cellid <<
", time = " << time <<
", calibratedTime = " << calibratedTime);
306 ECLDsp* aECLDsp = ECLDsp::getByCellID(cellid);
307 aECLCalDigit->setTwoComponentChi2(-1);
308 aECLCalDigit->setTwoComponentSavedChi2(ECLDsp::photonHadron, -1);
309 aECLCalDigit->setTwoComponentSavedChi2(ECLDsp::photonHadronBackgroundPhoton, -1);
310 aECLCalDigit->setTwoComponentSavedChi2(ECLDsp::photonDiodeCrossing, -1);
311 aECLCalDigit->setTwoComponentTotalEnergy(-1);
312 aECLCalDigit->setTwoComponentHadronEnergy(-1);
313 aECLCalDigit->setTwoComponentDiodeEnergy(-1);
319 const double calibratedTwoComponentTotalEnergy = aECLDsp->
getTwoComponentTotalAmp() * v_calibrationCrystalElectronics[cellid - 1] *
320 v_calibrationCrystalEnergy[cellid - 1];
321 const double calibratedTwoComponentHadronEnergy = aECLDsp->
getTwoComponentHadronAmp() * v_calibrationCrystalElectronics[cellid -
323 v_calibrationCrystalEnergy[cellid - 1];
324 const double calibratedTwoComponentDiodeEnergy = aECLDsp->
getTwoComponentDiodeAmp() * v_calibrationCrystalElectronics[cellid - 1] *
325 v_calibrationCrystalEnergy[cellid - 1];
329 aECLCalDigit->setTwoComponentTotalEnergy(calibratedTwoComponentTotalEnergy);
330 aECLCalDigit->setTwoComponentHadronEnergy(calibratedTwoComponentHadronEnergy);
331 aECLCalDigit->setTwoComponentDiodeEnergy(calibratedTwoComponentDiodeEnergy);
332 aECLCalDigit->setTwoComponentChi2(twoComponentChi2);
333 aECLCalDigit->setTwoComponentSavedChi2(ECLDsp::photonHadron, aECLDsp->
getTwoComponentSavedChi2(ECLDsp::photonHadron));
334 aECLCalDigit->setTwoComponentSavedChi2(ECLDsp::photonHadronBackgroundPhoton,
336 aECLCalDigit->setTwoComponentSavedChi2(ECLDsp::photonDiodeCrossing, aECLDsp->
getTwoComponentSavedChi2(ECLDsp::photonDiodeCrossing));
337 aECLCalDigit->setTwoComponentFitType(twoComponentFitType);
343 aECLCalDigit->setCellId(cellid);
345 aECLCalDigit->setEnergy(calibratedEnergy);
346 aECLCalDigit->addStatus(ECLCalDigit::c_IsEnergyCalibrated);
348 aECLCalDigit->setTime(calibratedTime);
349 aECLCalDigit->addStatus(ECLCalDigit::c_IsTimeCalibrated);
352 aECLCalDigit->addRelationTo(&aECLDigit);
356 const int bgCount = determineBackgroundECL();
359 for (
auto& aECLCalDigit : m_eclCalDigits) {
362 const double t99 = getT99(aECLCalDigit.getCellId(),
363 aECLCalDigit.getEnergy(),
364 aECLCalDigit.hasStatus(ECLCalDigit::c_IsFailedFit),
366 aECLCalDigit.setTimeResolution(t99);
368 if (t99 == c_timeResolutionForFitFailed) {
369 aECLCalDigit.addStatus(ECLCalDigit::c_IsFailedTimeResolution);
372 aECLCalDigit.addStatus(ECLCalDigit::c_IsTimeResolutionCalibrated);
377 void ECLDigitCalibratorModule::endRun()
383 void ECLDigitCalibratorModule::terminate()
389 double ECLDigitCalibratorModule::getT99(
const int cellid,
const double energy,
const bool fitfailed,
const int bgcount)
const
393 if (fitfailed)
return c_timeResolutionForFitFailed;
396 const double bglevel = TMath::Min((
double) bgcount / (
double) c_nominalBG * m_th1fBackground->GetBinContent(cellid) / m_averageBG,
400 const double p1 = c_pol2Var1 + c_pol2Var2 * bglevel + c_pol2Var3 * bglevel * bglevel;
407 double t99 = p1 * einv;
410 if (t99 < c_minT99) t99 = c_minT99;
412 B2DEBUG(35,
"ECLDigitCalibratorModule::getCalibratedTimeResolution: dose = " << m_th1fBackground->GetBinContent(
413 cellid) <<
", bglevel = " << bglevel <<
", cellid = " << cellid <<
", t99 = " << t99 <<
", energy = " << energy /
420 int ECLDigitCalibratorModule::determineBackgroundECL()
423 using regionCounter = std::unordered_map<ECL::DetectorRegion, uint>;
425 regionCounter outOfTimeCount {{ECL::DetectorRegion::FWD, 0},
426 {ECL::DetectorRegion::BRL, 0},
427 {ECL::DetectorRegion::BWD, 0}};
432 for (
auto& aECLCalDigit : m_eclCalDigits) {
433 if (abs(aECLCalDigit.getTime()) >= m_backgroundTimingCut) {
434 if (aECLCalDigit.getEnergy() >= m_backgroundEnergyCut) {
436 const B2Vector3D position = geom->GetCrystalPos(aECLCalDigit.getCellId() - 1);
437 const double theta = position.
Theta();
440 const auto detectorRegion = ECL::getDetectorRegion(theta);
443 ++outOfTimeCount.at(detectorRegion);
449 if (!m_eventLevelClusteringInfo) {
450 m_eventLevelClusteringInfo.create();
453 m_eventLevelClusteringInfo->setNECLCalDigitsOutOfTimeFWD(outOfTimeCount.at(ECL::DetectorRegion::FWD));
454 m_eventLevelClusteringInfo->setNECLCalDigitsOutOfTimeBarrel(outOfTimeCount.at(ECL::DetectorRegion::BRL));
455 m_eventLevelClusteringInfo->setNECLCalDigitsOutOfTimeBWD(outOfTimeCount.at(ECL::DetectorRegion::BWD));
457 B2DEBUG(35,
"ECLDigitCalibratorModule::determineBackgroundECL found " << outOfTimeCount.at(ECL::DetectorRegion::FWD) <<
", " <<
458 outOfTimeCount.at(ECL::DetectorRegion::BRL) <<
", " <<
459 outOfTimeCount.at(ECL::DetectorRegion::BWD) <<
" out of time digits in FWD, BRL, BWD");
461 return m_eventLevelClusteringInfo->getNECLCalDigitsOutOfTime();
DataType Theta() const
The polar angle.
Class for accessing objects in the database.
Class to find calibrate digits and convert waveform fit information to physics quantities.
Class to store ECL ShaperDSP waveform ADC data.
double getTwoComponentHadronAmp() const
get two comp hadron amp
double getTwoComponentSavedChi2(TwoComponentFitType FitTypeIn) const
get two comp chi2 for a fit type see enum TwoComponentFitType in ECLDsp.h for description of fit type...
TwoComponentFitType
Offline two component fit type.
double getTwoComponentTotalAmp() const
get two comp total amp
double getTwoComponentChi2() const
get two comp chi2
TwoComponentFitType getTwoComponentFitType() const
get two comp fit type
double getTwoComponentDiodeAmp() const
get two comp diode amp
Class to store ECL crystal type relation to ECLDigit for the simulation pure CsI upgrade option fille...
The Class for ECL Geometry Parameters.
static const double MeV
[megaelectronvolt]
static const double ns
Standard of [time].
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Abstract base class for different kinds of events.