Belle II Software  release-05-01-25
eclCosmicECollectorModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2018 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Christopher Hearty *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 //This module
11 #include <ecl/modules/eclCosmicECollector/eclCosmicECollectorModule.h>
12 
13 //Root
14 #include <TH2F.h>
15 
16 //Framework
17 #include <framework/dataobjects/EventMetaData.h>
18 
19 // ECL
20 #include <ecl/dbobjects/ECLCrystalCalib.h>
21 #include <ecl/dataobjects/ECLDigit.h>
22 #include <ecl/geometry/ECLGeometryPar.h>
23 #include <ecl/geometry/ECLNeighbours.h>
24 
25 //Trigger
26 #include <trg/ecl/TrgEclMapping.h>
27 
28 using namespace std;
29 using namespace Belle2;
30 using namespace ECL;
31 
32 
34 // Register the Module
36 REG_MODULE(eclCosmicECollector)
37 
38 
39 // Implementation
44 eclCosmicECollectorModule::eclCosmicECollectorModule() : CalibrationCollectorModule(), m_ECLExpCosmicESame("ECLExpCosmicESame"),
45  m_ECLExpCosmicEDifferent("ECLExpCosmicEDifferent"), m_ElectronicsCalib("ECLCrystalElectronics"),
46  m_CosmicECalib("ECLCrystalEnergyCosmic")
47 {
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);
54 }
55 
56 
59 void eclCosmicECollectorModule::prepare()
60 {
61  B2INFO("eclCosmicECollector: Experiment = " << m_evtMetaData->getExperiment() << " run = " << m_evtMetaData->getRun());
62 
64  m_eclDigitArray.isRequired();
65 
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);
71 
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);
75 
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);
79 
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);
83 
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);
87 
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);
91 
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);
95 
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);
99 
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);
103 
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);
107 
108  auto RawDigitAmpvsCrys = new TH2F("RawDigitAmpvsCrys", "Digit Amplitude vs crystal ID;crystal ID;Amplitude", 8736, 0, 8736, 100, 0,
109  2000);
110  registerObject<TH2F>("RawDigitAmpvsCrys", RawDigitAmpvsCrys);
111 
114  B2INFO("Input parameters to eclCosmicECollector:");
115  B2INFO("trigThreshold: " << m_trigThreshold);
116  B2INFO("minCrysE: " << m_minCrysE);
117 
119  FirstSet.resize(8737);
120  EperCrys.resize(8736);
121  HitCrys.resize(8736);
122  EnergyPerTC.resize(576);
123  TCperCrys.resize(8736);
124 
125 
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();}
132 
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]);
137  }
138 
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);}
145  }
146 
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
153  };
154 
156  short firstCrystal[69] = {};
157  for (int it = 1; it < 69; it++) {firstCrystal[it] = firstCrystal[it - 1] + m_crystalsPerRing[it - 1];}
158 
160  std::vector<short int> ThetaIDCrys(8736);
161  std::vector<short int> PhiIDCrys(8736);
162  int cID = 0;
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;
167  cID++;
168  }
169  }
170 
173  TrgEclMapping* trgecl_obj = new TrgEclMapping();
174  for (int crysID = 0; crysID < 8736; crysID++) {
175  TCperCrys[crysID] = trgecl_obj->getTCIdFromXtalId(crysID + 1) - 1;
176  }
183  int ThetaID = 0;
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;
187 
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);
197 
199  double dphiB = n2Overn0 * phiID + 0.0001;
200  phiIDB = dphiB;
201  CenterCrys.push_back(crysID);
202  NeighbourA.push_back(firstCrystal[ThetaID + 1] + phiID);
203  NeighbourB.push_back(firstCrystal[ThetaID + 2] + phiIDB);
204 
205  phiIDB++;
206  CenterCrys.push_back(crysID);
207  NeighbourA.push_back(firstCrystal[ThetaID + 1] + phiID);
208  NeighbourB.push_back(firstCrystal[ThetaID + 2] + phiIDB);
209  }
210 
215  ECLNeighbours* myNeighbours4 = new ECLNeighbours("F", 0.95);
216 
217  for (int crysID = firstCrystal[1]; crysID < firstCrystal[68]; crysID++) {
218  std::vector<short int> neighbours = myNeighbours4->getNeighbours(crysID + 1);
219 
221  int nA = -1;
222  int nB = -1;
223  std::vector<short int> nextThetaNeighbours;
224  std::vector<short int> previousThetaNeighbours;
225  for (auto& ID1 : neighbours) {
226  int temp0 = ID1 - 1;
227  if (temp0 != crysID && ThetaIDCrys[temp0] == ThetaIDCrys[crysID] && nA == -1) {
228  nA = temp0;
229  } else if (temp0 != crysID && ThetaIDCrys[temp0] == ThetaIDCrys[crysID] && nA != -1) {
230  nB = temp0;
231  }
232 
234  if (ThetaIDCrys[temp0] == ThetaIDCrys[crysID] + 1) {nextThetaNeighbours.push_back(temp0);}
235  if (ThetaIDCrys[temp0] == ThetaIDCrys[crysID] - 1) {previousThetaNeighbours.push_back(temp0);}
236  }
237 
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);
244  } else {
245  B2FATAL("No neighbour pair with the same thetaID for crysID = " << crysID);
246  }
247 
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);
254  }
255  }
256  }
257 
261  ThetaID = 68;
262  for (int phiID = 0; phiID < m_crystalsPerRing[ThetaID]; phiID++) {
263  int crysID = firstCrystal[ThetaID] + phiID;
264 
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);
274 
276  CenterCrys.push_back(crysID);
277  NeighbourA.push_back(firstCrystal[ThetaID - 1] + phiID);
278  NeighbourB.push_back(firstCrystal[ThetaID - 2] + phiID);
279  }
280 
282  FirstSet[8736] = CenterCrys.size();
283 }
284 
285 
288 void eclCosmicECollectorModule::collect()
289 {
290 
293  if (iEvent == 0) {
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);
299 
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);
304 
305  }
306  }
307 
309  /* Metadata; check if required DB objects have changed */
310  if (iEvent % 10000 == 0) {B2INFO("eclCosmicECollector: iEvent = " << iEvent);}
311  iEvent++;
312 
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();
319  }
320 
323  memset(&EperCrys[0], 0, EperCrys.size()*sizeof EperCrys[0]);
324  memset(&EnergyPerTC[0], 0, EnergyPerTC.size()*sizeof EnergyPerTC[0]);
325 
326 
327  for (auto& eclDigit : m_eclDigitArray) {
328  int crysID = eclDigit.getCellId() - 1;
329 
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());
334  }
335  for (int crysID = 0; crysID < 8736; crysID++) {HitCrys[crysID] = EperCrys[crysID] > m_minCrysE;}
336 
339  if (m_mockupL1) {
340  float maxTrig = 0.;
341  for (int tc = 0; tc < 576; tc++) {
342  if (EnergyPerTC[tc] > maxTrig) {maxTrig = EnergyPerTC[tc];}
343  }
344  bool ECLtrigger = maxTrig > m_trigThreshold;
345  if (!ECLtrigger) {return;}
346  }
347 
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]) {
354 
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);
361  } else {
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);
367  }
368  }
369  }
370 }
Belle2::CalibrationCollectorModule
Calibration collector module base class.
Definition: CalibrationCollectorModule.h:44
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::TrgEclMapping
A class of TC Mapping.
Definition: TrgEclMapping.h:31
Belle2::ECL::ECLNeighbours
Class to get the neighbours for a given cell id.
Definition: ECLNeighbours.h:38
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::eclCosmicECollectorModule
class eclCosmicECollectorModule.
Definition: eclCosmicECollectorModule.h:41
Belle2::TrgEclMapping::getTCIdFromXtalId
int getTCIdFromXtalId(int)
get [TC ID] from [Xtal ID]
Definition: TrgEclMapping.cc:36
Belle2::ECL::ECLNeighbours::getNeighbours
const std::vector< short int > & getNeighbours(short int cid) const
Return the neighbours for a given cell ID.
Definition: ECLNeighbours.cc:318