Belle II Software  release-05-01-25
eclGammaGammaECollectorModule.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/eclGammaGammaECollector/eclGammaGammaECollectorModule.h>
12 
13 //Root
14 #include <TH2F.h>
15 
16 //Analysis
17 #include <analysis/utility/PCmsLabTransform.h>
18 #include <analysis/ClusterUtility/ClusterUtils.h>
19 
20 //Framework
21 #include <framework/gearbox/Const.h>
22 #include <framework/dataobjects/EventMetaData.h>
23 #include <framework/datastore/RelationVector.h>
24 
25 //MDST
26 #include <mdst/dataobjects/TRGSummary.h>
27 #include <mdst/dataobjects/HitPatternCDC.h>
28 #include <mdst/dataobjects/TrackFitResult.h>
29 #include <mdst/dataobjects/Track.h>
30 #include <mdst/dataobjects/ECLCluster.h>
31 
32 //ECL
33 #include <ecl/dataobjects/ECLDigit.h>
34 #include <ecl/dataobjects/ECLCalDigit.h>
35 #include <ecl/dbobjects/ECLCrystalCalib.h>
36 
37 
38 
39 using namespace std;
40 using namespace Belle2;
41 
42 
43 //-----------------------------------------------------------------
44 // Register the Module
45 //-----------------------------------------------------------------
46 REG_MODULE(eclGammaGammaECollector)
47 
48 //-----------------------------------------------------------------
49 // Implementation
50 //-----------------------------------------------------------------
51 
52 //-----------------------------------------------------------------------------------------------------
53 
54 eclGammaGammaECollectorModule::eclGammaGammaECollectorModule() : CalibrationCollectorModule(),
55  m_ECLExpGammaGammaE("ECLExpGammaGammaE"),
56  m_ElectronicsCalib("ECLCrystalElectronics"), m_GammaGammaECalib("ECLCrystalEnergyGammaGamma")
57 {
58  // Set module properties
59  setDescription("Calibration Collector Module for ECL single crystal energy calibration using gamma gamma events");
60  setPropertyFlags(c_ParallelProcessingCertified);
61  addParam("thetaLabMinDeg", m_thetaLabMinDeg, "miniumum photon theta in lab (degrees)", 0.);
62  addParam("thetaLabMaxDeg", m_thetaLabMaxDeg, "maximum photon theta in lab (degrees)", 180.);
63  addParam("minPairMass", m_minPairMass, "minimum invariant mass of the pair of photons (GeV/c^2)", 9.);
64  addParam("mindPhi", m_mindPhi, "minimum delta phi between clusters (deg)", 179.);
65  addParam("maxTime", m_maxTime, "maximum (time-<t>)/dt99 of photons", 1.);
66  addParam("measureTrueEnergy", m_measureTrueEnergy, "use MC events to obtain expected energies", false);
67  addParam("requireL1", m_requireL1, "only use events that have a level 1 trigger", true);
68 }
69 
70 
71 
74 void eclGammaGammaECollectorModule::prepare()
75 {
76 
78  B2INFO("eclGammaGammaECollector: Experiment = " << m_evtMetaData->getExperiment() << " run = " << m_evtMetaData->getRun());
79 
82  auto EnVsCrysID = new TH2F("EnVsCrysID", "Normalized energy for each crystal;crystal ID;E/Expected", 8736, 0, 8736, 140, 0, 1.4);
83  registerObject<TH2F>("EnVsCrysID", EnVsCrysID);
84 
85  auto ExpEvsCrys = new TH1F("ExpEvsCrys", "Sum expected energy vs crystal ID;crystal ID;Energy (GeV)", 8736, 0, 8736);
86  registerObject<TH1F>("ExpEvsCrys", ExpEvsCrys);
87 
88  auto ElecCalibvsCrys = new TH1F("ElecCalibvsCrys", "Sum electronics calib const vs crystal ID;crystal ID;calibration constant",
89  8736, 0, 8736);
90  registerObject<TH1F>("ElecCalibvsCrys", ElecCalibvsCrys);
91 
92  auto InitialCalibvsCrys = new TH1F("InitialCalibvsCrys",
93  "Sum initial gamma gamma calib const vs crystal ID;crystal ID;calibration constant", 8736, 0, 8736);
94  registerObject<TH1F>("InitialCalibvsCrys", InitialCalibvsCrys);
95 
96  auto CalibEntriesvsCrys = new TH1F("CalibEntriesvsCrys", "Entries in calib vs crys histograms;crystal ID;Entries per crystal", 8736,
97  0, 8736);
98  registerObject<TH1F>("CalibEntriesvsCrys", CalibEntriesvsCrys);
99 
101  auto RawDigitAmpvsCrys = new TH2F("RawDigitAmpvsCrys", "Digit Amplitude vs crystal ID;crystal ID;Amplitude", 8736, 0, 8736, 200, 0,
102  200000);
103  registerObject<TH2F>("RawDigitAmpvsCrys", RawDigitAmpvsCrys);
104 
105  auto RawDigitTimevsCrys = new TH2F("RawDigitTimevsCrys", "Digit Time vs crystal ID;crystal ID;Time", 8736, 0, 8736, 200, -2000,
106  2000);
107  registerObject<TH2F>("RawDigitTimevsCrys", RawDigitTimevsCrys);
108 
109  auto TimeMinusAverage = new TH1F("TimeMinusAverage", "Photon time - average / dt99;(T-T_ave)/dt99 ", 100, -10, 10);
110  registerObject<TH1F>("TimeMinusAverage", TimeMinusAverage);
111 
112 
113  //------------------------------------------------------------------------
115  B2INFO("Input parameters to eclGammaGammaECollector:");
116  B2INFO("thetaLabMinDeg: " << m_thetaLabMinDeg);
117  B2INFO("thetaLabMaxDeg: " << m_thetaLabMaxDeg);
118  thetaLabMin = m_thetaLabMinDeg / TMath::RadToDeg();
119  thetaLabMax = m_thetaLabMaxDeg / TMath::RadToDeg();
120  B2INFO("minPairMass: " << m_minPairMass);
121  B2INFO("mindPhi: " << m_mindPhi);
122  B2INFO("maxTime: " << m_maxTime);
123  B2INFO("measureTrueEnergy: " << m_measureTrueEnergy);
124  B2INFO("requireL1: " << m_requireL1);
125 
127  EperCrys.resize(8736);
128 
131  if (m_ECLExpGammaGammaE.hasChanged()) {ExpGammaGammaE = m_ECLExpGammaGammaE->getCalibVector();}
132  if (m_ElectronicsCalib.hasChanged()) {ElectronicsCalib = m_ElectronicsCalib->getCalibVector();}
133  if (m_GammaGammaECalib.hasChanged()) {GammaGammaECalib = m_GammaGammaECalib->getCalibVector();}
134 
136  for (int ic = 1; ic < 9000; ic += 1000) {
137  B2INFO("DB constants for cellID=" << ic << ": ExpGammaGammaE = " << ExpGammaGammaE[ic - 1] << " ElectronicsCalib = " <<
138  ElectronicsCalib[ic - 1]
139  << " GammaGammaECalib = " << GammaGammaECalib[ic - 1]);
140  }
141 
143  for (int crysID = 0; crysID < 8736; crysID++) {
144  if (ElectronicsCalib[crysID] <= 0) {B2FATAL("eclGammaGammaECollector: ElectronicsCalib = " << ElectronicsCalib[crysID] << " for crysID = " << crysID);}
145  if (ExpGammaGammaE[crysID] == 0) {B2FATAL("eclGammaGammaECollector: ExpGammaGammaE = 0 for crysID = " << crysID);}
146  if (GammaGammaECalib[crysID] == 0) {B2FATAL("eclGammaGammaECollector: GammaGammaECalib = 0 for crysID = " << crysID);}
147  }
148 
151  m_trackArray.isRequired();
152  m_eclClusterArray.isRequired();
153  m_eclCalDigitArray.isRequired();
154  m_eclDigitArray.isRequired();
155 
156 }
157 
158 
161 void eclGammaGammaECollectorModule::collect()
162 {
163 
165  if (storeCalib) {
166  for (int crysID = 0; crysID < 8736; crysID++) {
167  getObjectPtr<TH1F>("ExpEvsCrys")->Fill(crysID + 0.001, ExpGammaGammaE[crysID]);
168  getObjectPtr<TH1F>("ElecCalibvsCrys")->Fill(crysID + 0.001, ElectronicsCalib[crysID]);
169  getObjectPtr<TH1F>("InitialCalibvsCrys")->Fill(crysID + 0.001, GammaGammaECalib[crysID]);
170  getObjectPtr<TH1F>("CalibEntriesvsCrys")->Fill(crysID + 0.001);
171  }
172  storeCalib = false;
173  }
174 
177  bool newConst = false;
178  if (m_ECLExpGammaGammaE.hasChanged()) {
179  newConst = true;
180  B2INFO("ECLExpGammaGammaE has changed, exp = " << m_evtMetaData->getExperiment() << " run = " << m_evtMetaData->getRun());
181  ExpGammaGammaE = m_ECLExpGammaGammaE->getCalibVector();
182  }
183  if (m_ElectronicsCalib.hasChanged()) {
184  newConst = true;
185  B2INFO("ECLCrystalElectronics has changed, exp = " << m_evtMetaData->getExperiment() << " run = " << m_evtMetaData->getRun());
186  ElectronicsCalib = m_ElectronicsCalib->getCalibVector();
187  }
188  if (m_GammaGammaECalib.hasChanged()) {
189  newConst = true;
190  B2INFO("ECLCrystalEnergyGammaGamma has changed, exp = " << m_evtMetaData->getExperiment() << " run = " << m_evtMetaData->getRun());
191  GammaGammaECalib = m_GammaGammaECalib->getCalibVector();
192  }
193 
194  if (newConst) {
195  for (int ic = 1; ic < 9000; ic += 1000) {
196  B2INFO("DB constants for cellID=" << ic << ": ExpGammaGammaE = " << ExpGammaGammaE[ic - 1] << " ElectronicsCalib = " <<
197  ElectronicsCalib[ic - 1]
198  << " GammaGammaECalib = " << GammaGammaECalib[ic - 1]);
199  }
200 
202  for (int crysID = 0; crysID < 8736; crysID++) {
203  if (ElectronicsCalib[crysID] <= 0) {B2FATAL("eclGammaGammaECollector: ElectronicsCalib = " << ElectronicsCalib[crysID] << " for crysID = " << crysID);}
204  if (ExpGammaGammaE[crysID] == 0) {B2FATAL("eclGammaGammaECollector: ExpGammaGammaE = 0 for crysID = " << crysID);}
205  if (GammaGammaECalib[crysID] == 0) {B2FATAL("eclGammaGammaECollector: GammaGammaECalib = 0 for crysID = " << crysID);}
206  }
207  }
208 
211  if (m_requireL1) {
212  unsigned int L1TriggerResults = m_TRGResults->getTRGSummary(0);
213  if (L1TriggerResults == 0) {return;}
214  }
215 
216  //------------------------------------------------------------------------
218  int nGoodTrk = 0;
219  for (auto& track : m_trackArray) {
220  const TrackFitResult* temptrackFit = track.getTrackFitResult(Const::ChargedStable(211));
221  if (temptrackFit) {
222  double pt = temptrackFit->getTransverseMomentum();
223  double z0 = temptrackFit->getZ0();
224  double d0 = temptrackFit->getD0();
225  double pValue = temptrackFit->getPValue();
226  int nCDChits = temptrackFit->getHitPatternCDC().getNHits();
227  if (pt > minTrkpt && abs(z0) < maxZ0 && abs(d0) < maxD0 && pValue > minpValue && nCDChits >= minCDChits) {nGoodTrk++;}
228  }
229  }
230  if (nGoodTrk > 0) {return;}
231 
232  //------------------------------------------------------------------------
234  const ECLCluster::EHypothesisBit usePhotons = ECLCluster::EHypothesisBit::c_nPhotons;
235  int icMax[2] = { -1, -1};
236  double maxClustE[2] = { -1., -1.};
237  int nclust = m_eclClusterArray.getEntries();
238  for (int ic = 0; ic < nclust; ic++) {
239  if (m_eclClusterArray[ic]->hasHypothesis(usePhotons)) {
240  double eClust = m_eclClusterArray[ic]->getEnergy(usePhotons);
241  if (eClust > maxClustE[0]) {
242  maxClustE[1] = maxClustE[0];
243  icMax[1] = icMax[0];
244  maxClustE[0] = eClust;
245  icMax[0] = ic;
246  } else if (eClust > maxClustE[1]) {
247  maxClustE[1] = eClust;
248  icMax[1] = ic;
249  }
250  }
251  }
252 
253  //------------------------------------------------------------------------
256  if (icMax[0] == -1 || icMax[1] == -1) {return;}
257  double theta0 = m_eclClusterArray[icMax[0]]->getTheta();
258  double theta1 = m_eclClusterArray[icMax[1]]->getTheta();
259  if (theta0 < thetaLabMin || theta0 > thetaLabMax || theta1 < thetaLabMin || theta1 > thetaLabMax) {return;}
260 
262  double t0 = m_eclClusterArray[icMax[0]]->getTime();
263  double t990 = m_eclClusterArray[icMax[0]]->getDeltaTime99();
264  double t1 = m_eclClusterArray[icMax[1]]->getTime();
265  double t991 = m_eclClusterArray[icMax[1]]->getDeltaTime99();
266  double taverage = (t0 / (t990 * t990) + t1 / (t991 * t991)) / (1. / (t990 * t990) + 1. / (t991 * t991));
267  if (abs(t0 - taverage) > t990 * m_maxTime || abs(t1 - taverage) > t991 * m_maxTime) {return;}
268 
270  ClusterUtils cUtil;
271  const TVector3 clustervertex = cUtil.GetIPPosition();
272 
273  double phi0 = m_eclClusterArray[icMax[0]]->getPhi();
274  TVector3 p30(0., 0., maxClustE[0]);
275  p30.SetTheta(theta0);
276  p30.SetPhi(phi0);
277  const TLorentzVector p40 = cUtil.Get4MomentumFromCluster(m_eclClusterArray[icMax[0]], clustervertex, usePhotons);
278 
279  double phi1 = m_eclClusterArray[icMax[1]]->getPhi();
280  TVector3 p31(0., 0., maxClustE[1]);
281  p31.SetTheta(theta1);
282  p31.SetPhi(phi1);
283  const TLorentzVector p41 = cUtil.Get4MomentumFromCluster(m_eclClusterArray[icMax[1]], clustervertex, usePhotons);
284 
285  double pairmass = (p40 + p41).M();
286  if (pairmass < m_minPairMass) {return;}
287 
289  PCmsLabTransform boostrotate;
290  TLorentzVector p40COM = boostrotate.rotateLabToCms() * p40;
291  TLorentzVector p41COM = boostrotate.rotateLabToCms() * p41;
292  double dphi = abs(p41COM.Phi() - p40COM.Phi()) * TMath::RadToDeg();
293  if (dphi > 180.) {dphi = 360. - dphi;}
294  if (dphi < m_mindPhi) {return;}
295 
296  //------------------------------------------------------------------------
298  memset(&EperCrys[0], 0, EperCrys.size()*sizeof EperCrys[0]);
299  for (auto& eclDigit : m_eclDigitArray) {
300  int crysID = eclDigit.getCellId() - 1;
301  getObjectPtr<TH2F>("RawDigitAmpvsCrys")->Fill(crysID + 0.001, eclDigit.getAmp());
302 
304  EperCrys[crysID] = eclDigit.getAmp() * abs(GammaGammaECalib[crysID]) * ElectronicsCalib[crysID];
305  if (EperCrys[crysID] > 0.01) {
306  getObjectPtr<TH2F>("RawDigitTimevsCrys")->Fill(crysID + 0.001, eclDigit.getTimeFit());
307  }
308  }
309 
311  if (m_measureTrueEnergy) {
312  for (auto& eclCalDigit : m_eclCalDigitArray) {
313  int tempCrysID = eclCalDigit.getCellId() - 1;
314  double tempE = eclCalDigit.getEnergy();
315  EperCrys[tempCrysID] = tempE;
316  }
317  }
318 
319  //------------------------------------------------------------------------
320  //** Find the most energetic crystal in each photon cluster */
321  int crysIDMax[2] = { -1, -1};
322  double crysEMax[2] = { -1., -1.};
323  for (int imax = 0; imax < 2; imax++) {
324  auto eclClusterRelations = m_eclClusterArray[icMax[imax]]->getRelationsTo<ECLCalDigit>("ECLCalDigits");
325  for (unsigned int ir = 0; ir < eclClusterRelations.size(); ir++) {
326  const auto calDigit = eclClusterRelations.object(ir);
327  int tempCrysID = calDigit->getCellId() - 1;
328  float tempE = EperCrys[tempCrysID];
329  if (tempE > crysEMax[imax]) {
330  crysEMax[imax] = tempE;
331  crysIDMax[imax] = tempCrysID;
332  }
333  }
334  }
335 
336  //------------------------------------------------------------------------
337  //** Store the normalized energies of the two crystals */
338  for (int ic = 0; ic < 2; ic++) {
339  if (crysIDMax[ic] >= 0) {
341  getObjectPtr<TH2F>("EnVsCrysID")->Fill(crysIDMax[ic] + 0.001, EperCrys[crysIDMax[ic]] / abs(ExpGammaGammaE[crysIDMax[ic]]));
342  getObjectPtr<TH1F>("ExpEvsCrys")->Fill(crysIDMax[ic] + 0.001, ExpGammaGammaE[crysIDMax[ic]]);
343  getObjectPtr<TH1F>("ElecCalibvsCrys")->Fill(crysIDMax[ic] + 0.001, ElectronicsCalib[crysIDMax[ic]]);
344  getObjectPtr<TH1F>("InitialCalibvsCrys")->Fill(crysIDMax[ic] + 0.001, GammaGammaECalib[crysIDMax[ic]]);
345  getObjectPtr<TH1F>("CalibEntriesvsCrys")->Fill(crysIDMax[ic] + 0.001);
346  }
347  }
348  if (crysIDMax[0] >= 0 && crysIDMax[1] >= 0) {
349  getObjectPtr<TH1F>("TimeMinusAverage")->Fill((t0 - taverage) / t990);
350  getObjectPtr<TH1F>("TimeMinusAverage")->Fill((t1 - taverage) / t991);
351  }
352 }
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::ECLCluster::EHypothesisBit
EHypothesisBit
The hypothesis bits for this ECLCluster (Connected region (CR) is split using this hypothesis.
Definition: ECLCluster.h:43
Belle2::TrackFitResult::getPValue
double getPValue() const
Getter for Chi2 Probability of the track fit.
Definition: TrackFitResult.h:163
Belle2::ClusterUtils::Get4MomentumFromCluster
const TLorentzVector Get4MomentumFromCluster(const ECLCluster *cluster, ECLCluster::EHypothesisBit hypo)
Returns four momentum vector.
Definition: ClusterUtils.cc:26
Belle2::TrackFitResult
Values of the result of a track fit with a given particle hypothesis.
Definition: TrackFitResult.h:59
Belle2::ClusterUtils::GetIPPosition
const TVector3 GetIPPosition()
Returns default IP position from beam parameters.
Definition: ClusterUtils.cc:182
Belle2::TrackFitResult::getHitPatternCDC
HitPatternCDC getHitPatternCDC() const
Getter for the hit pattern in the CDC;.
Definition: TrackFitResult.cc:120
Belle2::TrackFitResult::getZ0
double getZ0() const
Getter for z0.
Definition: TrackFitResult.h:200
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
Belle2::eclGammaGammaECollectorModule
Calibration collector module that uses e+e- --> gamma gamma to do ECL single crystal energy calibrati...
Definition: eclGammaGammaECollectorModule.h:42
Belle2::TrackFitResult::getTransverseMomentum
double getTransverseMomentum() const
Getter for the absolute value of the transverse momentum at the perigee.
Definition: TrackFitResult.h:140
Belle2::PCmsLabTransform
Class to hold Lorentz transformations from/to CMS and boost vector.
Definition: PCmsLabTransform.h:37
Belle2::PCmsLabTransform::rotateLabToCms
const TLorentzRotation rotateLabToCms() const
Returns Lorentz transformation from Lab to CMS.
Definition: PCmsLabTransform.h:74
Belle2::Const::ChargedStable
Provides a type-safe way to pass members of the chargedStableSet set.
Definition: Const.h:465
Belle2::HitPatternCDC::getNHits
unsigned short getNHits() const
Get the total Number of CDC hits in the fit.
Definition: HitPatternCDC.cc:49
Belle2::TrackFitResult::getD0
double getD0() const
Getter for d0.
Definition: TrackFitResult.h:178