Belle II Software  release-05-01-25
eclee5x5CollectorModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - 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/eclee5x5Collector/eclee5x5CollectorModule.h>
12 
13 //Root
14 #include <TH2F.h>
15 
16 //Analysis
17 #include <analysis/ClusterUtility/ClusterUtils.h>
18 
19 //Framework
20 #include <framework/dataobjects/EventMetaData.h>
21 #include <framework/datastore/RelationVector.h>
22 
23 //MDST
24 #include <mdst/dataobjects/TRGSummary.h>
25 #include <mdst/dataobjects/ECLCluster.h>
26 
27 //ECL
28 #include <ecl/dataobjects/ECLDigit.h>
29 #include <ecl/dataobjects/ECLCalDigit.h>
30 #include <ecl/dbobjects/ECLCrystalCalib.h>
31 
32 
33 using namespace std;
34 using namespace Belle2;
35 
36 
37 //-----------------------------------------------------------------
38 // Register the Module
39 //-----------------------------------------------------------------
40 REG_MODULE(eclee5x5Collector)
41 
42 //-----------------------------------------------------------------
43 // Implementation
44 //-----------------------------------------------------------------
45 
46 //-----------------------------------------------------------------------------------------------------
47 
48 eclee5x5CollectorModule::eclee5x5CollectorModule() : CalibrationCollectorModule(),
49  m_ECLExpee5x5E("ECLExpee5x5E"), m_ElectronicsCalib("ECLCrystalElectronics"), m_ee5x5Calib("ECLCrystalEnergyee5x5"),
50  m_selectdPhiData("ECLeedPhiData"), m_selectdPhiMC("ECLeedPhiMC")
51 {
52  // Set module properties
53  setDescription("Calibration Collector Module for ECL single crystal energy calibration using Bhabha events");
54  setPropertyFlags(c_ParallelProcessingCertified);
55  addParam("thetaLabMinDeg", m_thetaLabMinDeg, "miniumum ecl cluster theta in lab (degrees)", 17.);
56  addParam("thetaLabMaxDeg", m_thetaLabMaxDeg, "maximum ecl cluster theta in lab (degrees)", 150.);
57  addParam("minE0", m_minE0, "minimum energy of cluster 0: E*0/sqrts", 0.45);
58  addParam("minE1", m_minE1, "minimum energy of cluster 1: E*1/sqrts", 0.40);
59  addParam("maxdThetaSum", m_maxdThetaSum, "maximum diff between 180 deg and sum of cluster theta* (deg)", 2.);
60  addParam("dPhiScale", m_dPhiScale, "scale dPhi* cut by this factor", 1.);
61  addParam("maxTime", m_maxTime, "maximum cluster time diff abs(t1-t0)/dt99", 10.);
62  addParam("useCalDigits", m_useCalDigits, "use MC events to obtain expected energies", false);
63  addParam("requireL1", m_requireL1, "only use events that have a level 1 trigger", false);
64 }
65 
66 
67 
70 void eclee5x5CollectorModule::prepare()
71 {
72 
74  m_sqrts = m_boostrotate.getCMSEnergy();
75  B2INFO("eclee5x5Collector: Experiment = " << m_evtMetaData->getExperiment() << " run = " << m_evtMetaData->getRun() << " sqrts = "
76  << m_sqrts);
77 
78 
81  auto EnVsCrysID = new TH2F("EnVsCrysID", "Normalized 5x5 energy for each crystal;crystal ID;E25/Expected", 8736, 0, 8736, 200, 0.7,
82  1.2);
83  registerObject<TH2F>("EnVsCrysID", EnVsCrysID);
84 
85  auto RvsCrysID = new TH1F("RvsCrysID", "E_exp x E_crysID / sigma^2;crysID;sum of E_exp x E_crysID/sigma^2", 8736, 0, 8736);
86  registerObject("RvsCrysID", RvsCrysID);
87 
88  auto NRvsCrysID = new TH1F("NRvsCrysID", "Entries in RvsCrysID vs crysID;crysID;Entries in RvsCrysID", 8736, 0, 8736);
89  registerObject("NRvsCrysID", NRvsCrysID);
90 
91  auto Qmatrix = new TH2F("Qmatrix", "E_i x E_j/sigma^2;crysID i;crysID j", 8736, 0, 8736, 8736, 0, 8736);
92  registerObject("Qmatrix", Qmatrix);
93 
94 
95  auto ElecCalibvsCrys = new TH1F("ElecCalibvsCrys", "Sum electronics calib const vs crystal ID;crystal ID;calibration constant",
96  8736, 0, 8736);
97  registerObject<TH1F>("ElecCalibvsCrys", ElecCalibvsCrys);
98  auto ExpEvsCrys = new TH1F("ExpEvsCrys", "Sum expected energy calib const vs crystalID;crystal ID;calibration constant", 8736, 0,
99  8736);
100  registerObject<TH1F>("ExpEvsCrys", ExpEvsCrys);
101  auto InitialCalibvsCrys = new TH1F("InitialCalibvsCrys", "Sum initial calib const vs crystal ID;crystal ID;calibration constant",
102  8736, 0, 8736);
103  registerObject<TH1F>("InitialCalibvsCrys", InitialCalibvsCrys);
104 
105  auto CalibEntriesvsCrys = new TH1F("CalibEntriesvsCrys", "Entries in calib vs crys histograms;crystal ID;Entries per crystal", 8736,
106  0, 8736);
107  registerObject<TH1F>("CalibEntriesvsCrys", CalibEntriesvsCrys);
108 
109  auto EntriesvsCrys = new TH1F("EntriesvsCrys", "Selected Bhabha clusters vs crystal ID;crystal ID;Entries", 8736, 0, 8736);
110  registerObject<TH1F>("EntriesvsCrys", EntriesvsCrys);
111 
112  auto dPhivsThetaID = new TH2F("dPhivsThetaID",
113  "Phi* vs thetaID forward, pass thetaSum,E0,E1;thetaID of forward cluster;dPhi COM (deg)", 69, 0, 69, 150, 165, 180);
114  registerObject<TH2F>("dPhivsThetaID", dPhivsThetaID);
115 
116  //------------------------------------------------------------------------
118  B2INFO("Input parameters to eclee5x5Collector:");
119  B2INFO("thetaLabMinDeg: " << m_thetaLabMinDeg);
120  B2INFO("thetaLabMaxDeg: " << m_thetaLabMaxDeg);
121  m_thetaLabMin = m_thetaLabMinDeg / TMath::RadToDeg();
122  m_thetaLabMax = m_thetaLabMaxDeg / TMath::RadToDeg();
123  B2INFO("minE0: " << m_minE0);
124  B2INFO("minE1: " << m_minE1);
125  B2INFO("maxdThetaSum: " << m_maxdThetaSum);
126  B2INFO("dPhiScale: " << m_dPhiScale);
127  B2INFO("maxTime: " << m_maxTime);
128  B2INFO("useCalDigits: " << m_useCalDigits);
129  B2INFO("requireL1: " << m_requireL1);
130 
132  EperCrys.resize(8736);
133  m_thetaID.resize(8736);
134 
136  m_eclNeighbours5x5 = new ECL::ECLNeighbours("N", 2);
137 
140  if (m_ECLExpee5x5E.hasChanged()) {
141  Expee5x5E = m_ECLExpee5x5E->getCalibVector();
142  Expee5x5Sigma = m_ECLExpee5x5E->getCalibUncVector();
143  }
144  if (m_ElectronicsCalib.hasChanged()) {ElectronicsCalib = m_ElectronicsCalib->getCalibVector();}
145  if (m_ee5x5Calib.hasChanged()) {ee5x5Calib = m_ee5x5Calib->getCalibVector();}
146  if (m_selectdPhiMC.hasChanged() and m_useCalDigits) {
147  meandPhi = m_selectdPhiMC->getCalibVector();
148  widthdPhi = m_selectdPhiMC->getCalibUncVector();
149  } else if (m_selectdPhiData.hasChanged()) {
150  meandPhi = m_selectdPhiData->getCalibVector();
151  widthdPhi = m_selectdPhiData->getCalibUncVector();
152  }
153 
155  for (int ic = 1; ic < 9000; ic += 1000) {
156  B2INFO("DB constants for cellID=" << ic << ": ee5x5Calib = " << ee5x5Calib[ic - 1] << " Expee5x5E = " << Expee5x5E[ic - 1] <<
157  " ElectronicsCalib = " <<
158  ElectronicsCalib[ic - 1]);
159  }
160 
162  for (int crysID = 0; crysID < 8736; crysID++) {
163  if (ElectronicsCalib[crysID] <= 0) {B2FATAL("eclee5x5Collector: ElectronicsCalib = " << ElectronicsCalib[crysID] << " for crysID = " << crysID);}
164  if (Expee5x5E[crysID] == 0) {B2FATAL("eclee5x5Collector: Expee5x5E = 0 for crysID = " << crysID);}
165  if (ee5x5Calib[crysID] == 0) {B2FATAL("eclee5x5Collector: ee5x5Calib = 0 for crysID = " << crysID);}
166  }
167 
170  m_eclClusterArray.isRequired();
171  m_eclCalDigitArray.isRequired();
172  m_eclDigitArray.isRequired();
173  m_evtMetaData.isRequired();
174 
176  //..Derive ThetaID of each crystal, and cut on dPhi* as a function of thetaID
177  int crysID = 0;
178  for (int it = 0; it < 69; it++) {
179 
180  //..dPhi* cuts are actually a function of thetaID, not crysID
181  m_dPhiMin.push_back(meandPhi[crysID] - m_dPhiScale * widthdPhi[crysID]);
182  m_dPhiMax.push_back(meandPhi[crysID] + m_dPhiScale * widthdPhi[crysID]);
183  for (int ic = 0; ic < m_eclNeighbours5x5->getCrystalsPerRing(it); ic++) {
184  m_thetaID.at(crysID) = it;
185  crysID++;
186  }
187  }
188 
189 }
190 
191 
194 void eclee5x5CollectorModule::collect()
195 {
196 
198  if (storeCalib) {
199  for (int crysID = 0; crysID < 8736; crysID++) {
200  getObjectPtr<TH1F>("ExpEvsCrys")->Fill(crysID + 0.001, Expee5x5E[crysID]);
201  getObjectPtr<TH1F>("ElecCalibvsCrys")->Fill(crysID + 0.001, ElectronicsCalib[crysID]);
202  getObjectPtr<TH1F>("InitialCalibvsCrys")->Fill(crysID + 0.001, ee5x5Calib[crysID]);
203  getObjectPtr<TH1F>("CalibEntriesvsCrys")->Fill(crysID + 0.001);
204  }
205  storeCalib = false;
206  }
207 
210  bool newConst = false;
211  if (m_ECLExpee5x5E.hasChanged()) {
212  newConst = true;
213  B2INFO("ECLExpee5x5E has changed, exp = " << m_evtMetaData->getExperiment() << " run = " << m_evtMetaData->getRun());
214  Expee5x5E = m_ECLExpee5x5E->getCalibVector();
215  Expee5x5Sigma = m_ECLExpee5x5E->getCalibUncVector();
216  }
217  if (m_ElectronicsCalib.hasChanged()) {
218  newConst = true;
219  B2INFO("ECLCrystalElectronics has changed, exp = " << m_evtMetaData->getExperiment() << " run = " << m_evtMetaData->getRun());
220  ElectronicsCalib = m_ElectronicsCalib->getCalibVector();
221  }
222  if (m_ee5x5Calib.hasChanged()) {
223  newConst = true;
224  B2INFO("ECLCrystalEnergyee5x5 has changed, exp = " << m_evtMetaData->getExperiment() << " run = " << m_evtMetaData->getRun());
225  ee5x5Calib = m_ee5x5Calib->getCalibVector();
226  }
227 
228  if (newConst) {
229  for (int ic = 1; ic < 9000; ic += 1000) {
230  B2INFO("DB constants for cellID=" << ic << ": ee5x5Calib = " << ee5x5Calib[ic - 1] << " Expee5x5E = " << Expee5x5E[ic - 1] <<
231  " ElectronicsCalib = " <<
232  ElectronicsCalib[ic - 1]);
233  }
234 
236  for (int crysID = 0; crysID < 8736; crysID++) {
237  if (ElectronicsCalib[crysID] <= 0) {B2FATAL("eclee5x5Collector: ElectronicsCalib = " << ElectronicsCalib[crysID] << " for crysID = " << crysID);}
238  if (Expee5x5E[crysID] == 0) {B2FATAL("eclee5x5Collector: Expee5x5E = 0 for crysID = " << crysID);}
239  if (ee5x5Calib[crysID] == 0) {B2FATAL("eclee5x5Collector: ee5x5Calib = 0 for crysID = " << crysID);}
240  }
241  }
242 
245  if (m_requireL1) {
246  unsigned int L1TriggerResults = m_TRGResults->getTRGSummary(0);
247  if (L1TriggerResults == 0) {return;}
248  }
249 
250  //------------------------------------------------------------------------
252  int icMax[2] = { -1, -1};
253  double maxClustE[2] = { -1., -1.};
254  int nclust = m_eclClusterArray.getEntries();
255  for (int ic = 0; ic < nclust; ic++) {
256  if (m_eclClusterArray[ic]->hasHypothesis(Belle2::ECLCluster::EHypothesisBit::c_nPhotons)) {
257  double eClust = m_eclClusterArray[ic]->getEnergy(Belle2::ECLCluster::EHypothesisBit::c_nPhotons);
258  if (eClust > maxClustE[0]) {
259  maxClustE[1] = maxClustE[0];
260  icMax[1] = icMax[0];
261  maxClustE[0] = eClust;
262  icMax[0] = ic;
263  } else if (eClust > maxClustE[1]) {
264  maxClustE[1] = eClust;
265  icMax[1] = ic;
266  }
267  }
268  }
269 
270  //------------------------------------------------------------------------
273  if (icMax[1] == -1) {return;}
274  double theta0 = m_eclClusterArray[icMax[0]]->getTheta();
275  double theta1 = m_eclClusterArray[icMax[1]]->getTheta();
276  if (theta0 < m_thetaLabMin || theta0 > m_thetaLabMax || theta1 < m_thetaLabMin || theta1 > m_thetaLabMax) {return;}
277 
279  double t0 = m_eclClusterArray[icMax[0]]->getTime();
280  double dt990 = m_eclClusterArray[icMax[0]]->getDeltaTime99();
281  double t1 = m_eclClusterArray[icMax[1]]->getTime();
282  double dt991 = m_eclClusterArray[icMax[1]]->getDeltaTime99();
283  double dt99min = dt990;
284  if (dt991 < dt990) {dt99min = dt991;}
285  if (dt99min <= 0) {dt99min = 0.0001;}
286  if (abs(t1 - t0) > dt99min * m_maxTime) {return;}
287 
288  //------------------------------------------------------------------------
290  ClusterUtils cUtil;
291  const TVector3 clustervertex = cUtil.GetIPPosition();
292 
293  double phi0 = m_eclClusterArray[icMax[0]]->getPhi();
294  TVector3 p30(0., 0., maxClustE[0]);
295  p30.SetTheta(theta0);
296  p30.SetPhi(phi0);
297  const TLorentzVector p40 = cUtil.Get4MomentumFromCluster(m_eclClusterArray[icMax[0]], clustervertex,
298  ECLCluster::EHypothesisBit::c_nPhotons);
299 
300  double phi1 = m_eclClusterArray[icMax[1]]->getPhi();
301  TVector3 p31(0., 0., maxClustE[1]);
302  p31.SetTheta(theta1);
303  p31.SetPhi(phi1);
304  const TLorentzVector p41 = cUtil.Get4MomentumFromCluster(m_eclClusterArray[icMax[1]], clustervertex,
305  ECLCluster::EHypothesisBit::c_nPhotons);
306 
308  TLorentzVector p40COM = m_boostrotate.rotateLabToCms() * p40;
309  TLorentzVector p41COM = m_boostrotate.rotateLabToCms() * p41;
310  double theta01COM = (p41COM.Theta() + p40COM.Theta()) * TMath::RadToDeg();
311  if (abs(theta01COM - 180.) > m_maxdThetaSum) {return;}
312 
314  if (p40COM.E() < m_minE0 * m_sqrts and p41COM.E() < m_minE0 * m_sqrts) {return;}
315  if (p40COM.E() < m_minE1 * m_sqrts or p41COM.E() < m_minE1 * m_sqrts) {return;}
316 
317 
318  //------------------------------------------------------------------------
320  int crysIDMax[2] = { -1, -1};
321  double crysEMax[2] = { -1., -1.};
322  for (int imax = 0; imax < 2; imax++) {
323  auto eclClusterRelations = m_eclClusterArray[icMax[imax]]->getRelationsTo<ECLCalDigit>("ECLCalDigits");
324  for (unsigned int ir = 0; ir < eclClusterRelations.size(); ir++) {
325  const auto calDigit = eclClusterRelations.object(ir);
326  int tempCrysID = calDigit->getCellId() - 1;
327  float tempE = calDigit->getEnergy();
328  if (tempE > crysEMax[imax]) {
329  crysEMax[imax] = tempE;
330  crysIDMax[imax] = tempCrysID;
331  }
332  }
333  }
334 
336  double dphiCOM = abs(p41COM.Phi() - p40COM.Phi()) * TMath::RadToDeg();
337  if (dphiCOM > 180.) {dphiCOM = 360. - dphiCOM;}
338 
339  int thetaIDmin = m_thetaID[crysIDMax[0]];
340  if (m_thetaID[crysIDMax[1]] < m_thetaID[crysIDMax[0]]) {thetaIDmin = m_thetaID[crysIDMax[1]];}
341  getObjectPtr<TH2F>("dPhivsThetaID")->Fill(thetaIDmin + 0.001, dphiCOM);
342  if (dphiCOM<m_dPhiMin.at(thetaIDmin) or dphiCOM>m_dPhiMax.at(thetaIDmin)) {return;}
343 
344 
345  //------------------------------------------------------------------------
347  memset(&EperCrys[0], 0, EperCrys.size()*sizeof EperCrys[0]);
348 
349  if (m_useCalDigits) {
350  for (auto& eclCalDigit : m_eclCalDigitArray) {
351  int tempCrysID = eclCalDigit.getCellId() - 1;
352  EperCrys[tempCrysID] = eclCalDigit.getEnergy();
353  }
354  } else {
355  for (auto& eclDigit : m_eclDigitArray) {
356  int tempCrysID = eclDigit.getCellId() - 1;
358  EperCrys[tempCrysID] = eclDigit.getAmp() * abs(ee5x5Calib[tempCrysID]) * ElectronicsCalib[tempCrysID];
359  }
360  }
361 
362  //------------------------------------------------------------------------
363  //** Quantities needed for the 5x5 calibration */
364  for (int ic = 0; ic < 2; ic++) {
365  int crysMax = crysIDMax[ic];
366  float expE = abs(Expee5x5E[crysMax]);
367  float sigmaExp = Expee5x5Sigma[crysMax];
368  std::vector<short int> neighbours = m_eclNeighbours5x5->getNeighbours(crysMax + 1);
369 
370  //** Energy in 5x5, and expected energy corrected for crystals that will not be calibrated */
371  double reducedExpE = expE;
372  double E25 = 0.;
373  for (auto& cellID : neighbours) {
374  E25 += EperCrys[cellID - 1];
375  if (ee5x5Calib[cellID - 1] < 0.) {
376  reducedExpE -= EperCrys[cellID - 1];
377  }
378  }
379 
380  //** now the vector and matrix used in the calibration */
381  double rexpE = reducedExpE / sigmaExp;
382  for (auto& celli : neighbours) {
383  if (ee5x5Calib[celli - 1] > 0.) {
384  float rEi = EperCrys[celli - 1] / sigmaExp;
385  getObjectPtr<TH1F>("RvsCrysID")->Fill(celli - 0.999, rexpE * rEi);
386  getObjectPtr<TH1F>("NRvsCrysID")->Fill(celli - 0.999);
387  for (auto& cellj : neighbours) {
388  if (ee5x5Calib[cellj - 1] > 0.) {
389  float rEj = EperCrys[cellj - 1] / sigmaExp;
390  getObjectPtr<TH2F>("Qmatrix")->Fill(celli - 0.999, cellj - 0.999, rEi * rEj);
391  }
392  }
393  }
394  }
395 
396  //** Record normalized energy and corresponding calib values. Expee5x5E is negative if the algorithm was unable to calculate a value. In this case, the nominal input value has been stored with a minus sign */
397  getObjectPtr<TH2F>("EnVsCrysID")->Fill(crysMax + 0.001, E25 / expE);
398  getObjectPtr<TH1F>("EntriesvsCrys")->Fill(crysMax + 0.001);
399  getObjectPtr<TH1F>("ExpEvsCrys")->Fill(crysMax + 0.001, Expee5x5E[crysMax]);
400  getObjectPtr<TH1F>("ElecCalibvsCrys")->Fill(crysMax + 0.001, ElectronicsCalib[crysMax]);
401  getObjectPtr<TH1F>("InitialCalibvsCrys")->Fill(crysMax + 0.001, ee5x5Calib[crysMax]);
402  getObjectPtr<TH1F>("CalibEntriesvsCrys")->Fill(crysMax + 0.001);
403  }
404 }
Belle2::CalibrationCollectorModule
Calibration collector module base class.
Definition: CalibrationCollectorModule.h:44
Belle2::ECLCalDigit
Class to store calibrated ECLDigits: ECLCalDigits.
Definition: ECLCalDigit.h:38
Belle2::ECLCalDigit::getCellId
int getCellId() const
Get Cell ID.
Definition: ECLCalDigit.h:129
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::ClusterUtils::Get4MomentumFromCluster
const TLorentzVector Get4MomentumFromCluster(const ECLCluster *cluster, ECLCluster::EHypothesisBit hypo)
Returns four momentum vector.
Definition: ClusterUtils.cc:26
Belle2::ECLCluster::EHypothesisBit::c_nPhotons
@ c_nPhotons
CR is split into n photons (N1)
Belle2::ClusterUtils::GetIPPosition
const TVector3 GetIPPosition()
Returns default IP position from beam parameters.
Definition: ClusterUtils.cc:182
Belle2::eclee5x5CollectorModule
Calibration collector module that uses e+e- --> e+e- to do ECL single crystal energy calibration.
Definition: eclee5x5CollectorModule.h:46
Belle2::ECL::ECLNeighbours
Class to get the neighbours for a given cell id.
Definition: ECLNeighbours.h:38
Belle2::ClusterUtils
Class to provide momentum-related information from ECLClusters.
Definition: ClusterUtils.h:44
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19