11 #include <ecl/modules/eclCosmicECollector/eclCosmicECollectorModule.h>
17 #include <framework/dataobjects/EventMetaData.h>
20 #include <ecl/dbobjects/ECLCrystalCalib.h>
21 #include <ecl/dataobjects/ECLDigit.h>
22 #include <ecl/geometry/ECLGeometryPar.h>
23 #include <ecl/geometry/ECLNeighbours.h>
26 #include <trg/ecl/TrgEclMapping.h>
45 m_ECLExpCosmicEDifferent("ECLExpCosmicEDifferent"), m_ElectronicsCalib("ECLCrystalElectronics"),
46 m_CosmicECalib("ECLCrystalEnergyCosmic")
49 setDescription(
"Calibration Collector Module for ECL single crystal energy calibration using cosmic rays");
50 setPropertyFlags(c_ParallelProcessingCertified);
51 addParam(
"minCrysE", m_minCrysE,
"Minimum energy in GeV for a crystal to be considered hit", 0.010);
52 addParam(
"mockupL1", m_mockupL1,
"Calculate energy per trigger cell in lieu of trigger simulation",
false);
53 addParam(
"trigThreshold", m_trigThreshold,
"Minimum energy in GeV per trigger cell to mock up L1 trigger", 0.1);
59 void eclCosmicECollectorModule::prepare()
61 B2INFO(
"eclCosmicECollector: Experiment = " << m_evtMetaData->getExperiment() <<
" run = " << m_evtMetaData->getRun());
64 m_eclDigitArray.isRequired();
68 auto EnvsCrysSameRing =
new TH2F(
"EnvsCrysSameRing",
"Normalized energy vs crystal ID, same theta ring;crystal ID;E/Expected", 8736,
69 0, 8736, 99, .025, 2.5);
70 registerObject<TH2F>(
"EnvsCrysSameRing", EnvsCrysSameRing);
72 auto EnvsCrysDifferentRing =
new TH2F(
"EnvsCrysDifferentRing",
73 "Normalized energy vs crystal ID, different theta rings;crystal ID;E/Expected", 8736, 0, 8736, 99, .025, 2.5);
74 registerObject<TH2F>(
"EnvsCrysDifferentRing", EnvsCrysDifferentRing);
76 auto ExpEvsCrysSameRing =
new TH1F(
"ExpEvsCrysSameRing",
77 "Sum expected energy vs crystal ID, same theta ring;crystal ID;Energy (GeV)", 8736, 0, 8736);
78 registerObject<TH1F>(
"ExpEvsCrysSameRing", ExpEvsCrysSameRing);
80 auto ExpEvsCrysDifferentRing =
new TH1F(
"ExpEvsCrysDifferentRing",
81 "Sum expected energy vs crystal ID, different theta rings;crystal ID;Energy (GeV)", 8736, 0, 8736);
82 registerObject<TH1F>(
"ExpEvsCrysDifferentRing", ExpEvsCrysDifferentRing);
84 auto ElecCalibvsCrysSameRing =
new TH1F(
"ElecCalibvsCrysSameRing",
85 "Sum electronics calib const vs crystal ID, same theta ring;crystal ID", 8736, 0, 8736);
86 registerObject<TH1F>(
"ElecCalibvsCrysSameRing", ElecCalibvsCrysSameRing);
88 auto ElecCalibvsCrysDifferentRing =
new TH1F(
"ElecCalibvsCrysDifferentRing",
89 "Sum electronics calib const vs crystal ID, different theta rings;crystal ID", 8736, 0, 8736);
90 registerObject<TH1F>(
"ElecCalibvsCrysDifferentRing", ElecCalibvsCrysDifferentRing);
92 auto InitialCalibvsCrysSameRing =
new TH1F(
"InitialCalibvsCrysSameRing",
93 "Sum initial cosmic calib const vs crystal ID, same theta ring;crystal ID", 8736, 0, 8736);
94 registerObject<TH1F>(
"InitialCalibvsCrysSameRing", InitialCalibvsCrysSameRing);
96 auto InitialCalibvsCrysDifferentRing =
new TH1F(
"InitialCalibvsCrysDifferentRing",
97 "Sum initial cosmic calib const vs crystal ID, different theta rings;crystal ID", 8736, 0, 8736);
98 registerObject<TH1F>(
"InitialCalibvsCrysDifferentRing", InitialCalibvsCrysDifferentRing);
100 auto CalibEntriesvsCrysSameRing =
new TH1F(
"CalibEntriesvsCrysSameRing",
101 "Entries in calib vs crys histograms, same theta ring;crystal ID;Entries per crystal", 8736, 0, 8736);
102 registerObject<TH1F>(
"CalibEntriesvsCrysSameRing", CalibEntriesvsCrysSameRing);
104 auto CalibEntriesvsCrysDifferentRing =
new TH1F(
"CalibEntriesvsCrysDifferentRing",
105 "Entries in calib vs crys histograms, different theta rings;crystal ID;Entries per crystal", 8736, 0, 8736);
106 registerObject<TH1F>(
"CalibEntriesvsCrysDifferentRing", CalibEntriesvsCrysDifferentRing);
108 auto RawDigitAmpvsCrys =
new TH2F(
"RawDigitAmpvsCrys",
"Digit Amplitude vs crystal ID;crystal ID;Amplitude", 8736, 0, 8736, 100, 0,
110 registerObject<TH2F>(
"RawDigitAmpvsCrys", RawDigitAmpvsCrys);
114 B2INFO(
"Input parameters to eclCosmicECollector:");
115 B2INFO(
"trigThreshold: " << m_trigThreshold);
116 B2INFO(
"minCrysE: " << m_minCrysE);
119 FirstSet.resize(8737);
120 EperCrys.resize(8736);
121 HitCrys.resize(8736);
122 EnergyPerTC.resize(576);
123 TCperCrys.resize(8736);
128 if (m_ECLExpCosmicESame.hasChanged()) {ExpCosmicESame = m_ECLExpCosmicESame->getCalibVector();}
129 if (m_ECLExpCosmicEDifferent.hasChanged()) {ExpCosmicEDifferent = m_ECLExpCosmicEDifferent->getCalibVector();}
130 if (m_ElectronicsCalib.hasChanged()) {ElectronicsCalib = m_ElectronicsCalib->getCalibVector();}
131 if (m_CosmicECalib.hasChanged()) {CosmicECalib = m_CosmicECalib->getCalibVector();}
134 for (
int ic = 1; ic < 9000; ic += 1000) {
135 B2INFO(
"DB constants for cellID=" << ic <<
": ExpCosmicESame = " << ExpCosmicESame[ic - 1] <<
" ExpCosmicEDifferent = " <<
136 ExpCosmicEDifferent[ic - 1] <<
" ElectronicsCalib = " << ElectronicsCalib[ic - 1] <<
" CosmicECalib = " << CosmicECalib[ic - 1]);
140 for (
int ic = 0; ic < 8736; ic++) {
141 if (ElectronicsCalib[ic] <= 0) {B2FATAL(
"eclCosmicECollector: ElectronicsCalib = " << ElectronicsCalib[ic] <<
" for crysID = " << ic);}
142 if (ExpCosmicESame[ic] == 0) {B2FATAL(
"eclCosmicECollector: ExpCosmicESame = 0 for crysID = " << ic);}
143 if (ExpCosmicEDifferent[ic] == 0) {B2FATAL(
"eclCosmicECollector: ExpCosmicEDifferent = 0 for crysID = " << ic);}
144 if (CosmicECalib[ic] == 0) {B2FATAL(
"eclCosmicECollector: CosmicECalib = 0 for crysID = " << ic);}
151 const short m_crystalsPerRing[69] = {
152 48, 48, 64, 64, 64, 96, 96, 96, 96, 96, 96, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 144, 96, 96, 96, 96, 96, 64, 64, 64
156 short firstCrystal[69] = {};
157 for (
int it = 1; it < 69; it++) {firstCrystal[it] = firstCrystal[it - 1] + m_crystalsPerRing[it - 1];}
160 std::vector<short int> ThetaIDCrys(8736);
161 std::vector<short int> PhiIDCrys(8736);
163 for (
int it = 0; it < 69; it++) {
164 for (
int ic = 0; ic < m_crystalsPerRing[it]; ic++) {
165 ThetaIDCrys.at(cID) = it;
166 PhiIDCrys.at(cID) = ic;
174 for (
int crysID = 0; crysID < 8736; crysID++) {
184 double n2Overn0 = (double)m_crystalsPerRing[ThetaID + 2] / (
double)m_crystalsPerRing[ThetaID];
185 for (
int phiID = 0; phiID < m_crystalsPerRing[ThetaID]; phiID++) {
186 int crysID = firstCrystal[ThetaID] + phiID;
189 int phiIDA = phiID + 1;
190 if (phiIDA == m_crystalsPerRing[ThetaID]) {phiIDA = 0;}
191 int phiIDB = phiID - 1;
192 if (phiIDB == -1) {phiIDB = m_crystalsPerRing[ThetaID] - 1;}
193 FirstSet[crysID] = CenterCrys.size();
194 CenterCrys.push_back(crysID);
195 NeighbourA.push_back(firstCrystal[ThetaID] + phiIDA);
196 NeighbourB.push_back(firstCrystal[ThetaID] + phiIDB);
199 double dphiB = n2Overn0 * phiID + 0.0001;
201 CenterCrys.push_back(crysID);
202 NeighbourA.push_back(firstCrystal[ThetaID + 1] + phiID);
203 NeighbourB.push_back(firstCrystal[ThetaID + 2] + phiIDB);
206 CenterCrys.push_back(crysID);
207 NeighbourA.push_back(firstCrystal[ThetaID + 1] + phiID);
208 NeighbourB.push_back(firstCrystal[ThetaID + 2] + phiIDB);
217 for (
int crysID = firstCrystal[1]; crysID < firstCrystal[68]; crysID++) {
218 std::vector<short int> neighbours = myNeighbours4->
getNeighbours(crysID + 1);
223 std::vector<short int> nextThetaNeighbours;
224 std::vector<short int> previousThetaNeighbours;
225 for (
auto& ID1 : neighbours) {
227 if (temp0 != crysID && ThetaIDCrys[temp0] == ThetaIDCrys[crysID] && nA == -1) {
229 }
else if (temp0 != crysID && ThetaIDCrys[temp0] == ThetaIDCrys[crysID] && nA != -1) {
234 if (ThetaIDCrys[temp0] == ThetaIDCrys[crysID] + 1) {nextThetaNeighbours.push_back(temp0);}
235 if (ThetaIDCrys[temp0] == ThetaIDCrys[crysID] - 1) {previousThetaNeighbours.push_back(temp0);}
239 if (nA >= 0 && nB >= 0) {
240 FirstSet[crysID] = CenterCrys.size();
241 CenterCrys.push_back(crysID);
242 NeighbourA.push_back(nA);
243 NeighbourB.push_back(nB);
245 B2FATAL(
"No neighbour pair with the same thetaID for crysID = " << crysID);
249 for (
auto& IDn : nextThetaNeighbours) {
250 for (
auto& IDp : previousThetaNeighbours) {
251 CenterCrys.push_back(crysID);
252 NeighbourA.push_back(IDn);
253 NeighbourB.push_back(IDp);
262 for (
int phiID = 0; phiID < m_crystalsPerRing[ThetaID]; phiID++) {
263 int crysID = firstCrystal[ThetaID] + phiID;
266 int phiIDA = phiID + 1;
267 if (phiIDA == m_crystalsPerRing[ThetaID]) {phiIDA = 0;}
268 int phiIDB = phiID - 1;
269 if (phiIDB == -1) {phiIDB = m_crystalsPerRing[ThetaID] - 1;}
270 FirstSet[crysID] = CenterCrys.size();
271 CenterCrys.push_back(crysID);
272 NeighbourA.push_back(firstCrystal[ThetaID] + phiIDA);
273 NeighbourB.push_back(firstCrystal[ThetaID] + phiIDB);
276 CenterCrys.push_back(crysID);
277 NeighbourA.push_back(firstCrystal[ThetaID - 1] + phiID);
278 NeighbourB.push_back(firstCrystal[ThetaID - 2] + phiID);
282 FirstSet[8736] = CenterCrys.size();
288 void eclCosmicECollectorModule::collect()
294 for (
int crysID = 0; crysID < 8736; crysID++) {
295 getObjectPtr<TH1F>(
"ExpEvsCrysSameRing")->Fill(crysID + 0.001, ExpCosmicESame[crysID]);
296 getObjectPtr<TH1F>(
"ElecCalibvsCrysSameRing")->Fill(crysID + 0.001, ElectronicsCalib[crysID]);
297 getObjectPtr<TH1F>(
"InitialCalibvsCrysSameRing")->Fill(crysID + 0.001, CosmicECalib[crysID]);
298 getObjectPtr<TH1F>(
"CalibEntriesvsCrysSameRing")->Fill(crysID + 0.001);
300 getObjectPtr<TH1F>(
"ExpEvsCrysDifferentRing")->Fill(crysID + 0.001, ExpCosmicEDifferent[crysID]);
301 getObjectPtr<TH1F>(
"ElecCalibvsCrysDifferentRing")->Fill(crysID + 0.001, ElectronicsCalib[crysID]);
302 getObjectPtr<TH1F>(
"InitialCalibvsCrysDifferentRing")->Fill(crysID + 0.001, CosmicECalib[crysID]);
303 getObjectPtr<TH1F>(
"CalibEntriesvsCrysDifferentRing")->Fill(crysID + 0.001);
310 if (iEvent % 10000 == 0) {B2INFO(
"eclCosmicECollector: iEvent = " << iEvent);}
313 if (m_ECLExpCosmicESame.hasChanged()) {B2FATAL(
"eclCosmicECollector: ECLExpCosmicESame has changed");}
314 if (m_ECLExpCosmicEDifferent.hasChanged()) {B2FATAL(
"eclCosmicECollector: ECLExpCosmicEDifferent has changed");}
315 if (m_ElectronicsCalib.hasChanged()) {B2FATAL(
"eclCosmicECollector: ElectronicsCalib has changed");}
316 if (m_CosmicECalib.hasChanged()) {
317 B2DEBUG(9,
"eclCosmicECollector: new values for ECLCosmicECalib");
318 CosmicECalib = m_CosmicECalib->getCalibVector();
323 memset(&EperCrys[0], 0, EperCrys.size()*
sizeof EperCrys[0]);
324 memset(&EnergyPerTC[0], 0, EnergyPerTC.size()*
sizeof EnergyPerTC[0]);
327 for (
auto& eclDigit : m_eclDigitArray) {
328 int crysID = eclDigit.getCellId() - 1;
331 EperCrys[crysID] = eclDigit.getAmp() * abs(CosmicECalib[crysID]) * ElectronicsCalib[crysID];
332 EnergyPerTC[TCperCrys[crysID]] += EperCrys[crysID];
333 getObjectPtr<TH2F>(
"RawDigitAmpvsCrys")->Fill(crysID + 0.001, eclDigit.getAmp());
335 for (
int crysID = 0; crysID < 8736; crysID++) {HitCrys[crysID] = EperCrys[crysID] > m_minCrysE;}
341 for (
int tc = 0; tc < 576; tc++) {
342 if (EnergyPerTC[tc] > maxTrig) {maxTrig = EnergyPerTC[tc];}
344 bool ECLtrigger = maxTrig > m_trigThreshold;
345 if (!ECLtrigger) {
return;}
350 for (
int i = 0; i < FirstSet[8736]; i++) {
351 if (HitCrys[NeighbourA[i]] && HitCrys[NeighbourB[i]]) {
352 int crysID = CenterCrys[i];
353 if (i == FirstSet[crysID]) {
356 getObjectPtr<TH2F>(
"EnvsCrysSameRing")->Fill(crysID + 0.001, EperCrys[crysID] / abs(ExpCosmicESame[crysID]));
357 getObjectPtr<TH1F>(
"ExpEvsCrysSameRing")->Fill(crysID + 0.001, ExpCosmicESame[crysID]);
358 getObjectPtr<TH1F>(
"ElecCalibvsCrysSameRing")->Fill(crysID + 0.001, ElectronicsCalib[crysID]);
359 getObjectPtr<TH1F>(
"InitialCalibvsCrysSameRing")->Fill(crysID + 0.001, CosmicECalib[crysID]);
360 getObjectPtr<TH1F>(
"CalibEntriesvsCrysSameRing")->Fill(crysID + 0.001);
362 getObjectPtr<TH2F>(
"EnvsCrysDifferentRing")->Fill(crysID + 0.001, EperCrys[crysID] / abs(ExpCosmicEDifferent[crysID]));
363 getObjectPtr<TH1F>(
"ExpEvsCrysDifferentRing")->Fill(crysID + 0.001, ExpCosmicEDifferent[crysID]);
364 getObjectPtr<TH1F>(
"ElecCalibvsCrysDifferentRing")->Fill(crysID + 0.001, ElectronicsCalib[crysID]);
365 getObjectPtr<TH1F>(
"InitialCalibvsCrysDifferentRing")->Fill(crysID + 0.001, CosmicECalib[crysID]);
366 getObjectPtr<TH1F>(
"CalibEntriesvsCrysDifferentRing")->Fill(crysID + 0.001);