Belle II Software  release-05-01-25
eclMuMuECollectorModule.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/eclMuMuECollector/eclMuMuECollectorModule.h>
12 
13 //ROOT
14 #include <TH2F.h>
15 
16 //Framework
17 #include <framework/gearbox/Const.h>
18 #include <framework/dataobjects/EventMetaData.h>
19 
20 //Tracking
21 #include <tracking/dataobjects/ExtHit.h>
22 
23 //ECL
24 #include <ecl/dataobjects/ECLDigit.h>
25 #include <ecl/dataobjects/ECLCalDigit.h>
26 #include <ecl/geometry/ECLNeighbours.h>
27 #include <ecl/dbobjects/ECLCrystalCalib.h>
28 
29 //MDST
30 #include <mdst/dataobjects/Track.h>
31 #include <mdst/dataobjects/TRGSummary.h>
32 
33 using namespace std;
34 using namespace Belle2;
35 using namespace ECL;
36 
37 
38 //-----------------------------------------------------------------
39 // Register the Module
40 //-----------------------------------------------------------------
41 REG_MODULE(eclMuMuECollector)
42 
43 //-----------------------------------------------------------------
44 // Implementation
45 //-----------------------------------------------------------------
46 
47 //-----------------------------------------------------------------------------------------------------
48 
49 eclMuMuECollectorModule::eclMuMuECollectorModule() : CalibrationCollectorModule(), m_ECLExpMuMuE("ECLExpMuMuE"),
50  m_ElectronicsCalib("ECLCrystalElectronics"), m_MuMuECalib("ECLCrystalEnergyMuMu")
51 {
52  // Set module properties
53  setDescription("Calibration Collector Module for ECL single crystal energy calibration using muons");
54  setPropertyFlags(c_ParallelProcessingCertified);
55  addParam("minPairMass", m_minPairMass, "minimum invariant mass of the muon pair (GeV/c^2)", 9.0);
56  addParam("minTrackLength", m_minTrackLength, "minimum extrapolated track length in the crystal (cm)", 30.);
57  addParam("MaxNeighbourE", m_MaxNeighbourE, "maximum energy allowed in a neighbouring crystal (GeV)", 0.010);
58  addParam("thetaLabMinDeg", m_thetaLabMinDeg, "miniumum muon theta in lab (degrees)", 17.);
59  addParam("thetaLabMaxDeg", m_thetaLabMaxDeg, "maximum muon theta in lab (degrees)", 150.);
60  addParam("measureTrueEnergy", m_measureTrueEnergy, "use MC events to obtain expected energies", false);
61  addParam("requireL1", m_requireL1, "only use events that have a level 1 trigger", true);
62 }
63 
66 void eclMuMuECollectorModule::prepare()
67 {
68 
72  B2INFO("eclMuMuECollector: Experiment = " << m_evtMetaData->getExperiment() << " run = " << m_evtMetaData->getRun());
73 
76  auto EnVsCrysID = new TH2F("EnVsCrysID", "Normalized energy for each crystal;crystal ID;E/Expected", 8736, 0, 8736, 120, 0.1, 2.5);
77  registerObject<TH2F>("EnVsCrysID", EnVsCrysID);
78 
79  auto ExpEvsCrys = new TH1F("ExpEvsCrys", "Sum expected energy vs crystal ID;crystal ID;Energy (GeV)", 8736, 0, 8736);
80  registerObject<TH1F>("ExpEvsCrys", ExpEvsCrys);
81 
82  auto ElecCalibvsCrys = new TH1F("ElecCalibvsCrys", "Sum electronics calib const vs crystal ID;crystal ID;calibration constant",
83  8736, 0, 8736);
84  registerObject<TH1F>("ElecCalibvsCrys", ElecCalibvsCrys);
85 
86  auto InitialCalibvsCrys = new TH1F("InitialCalibvsCrys",
87  "Sum initial muon pair calib const vs crystal ID;crystal ID;calibration constant", 8736, 0, 8736);
88  registerObject<TH1F>("InitialCalibvsCrys", InitialCalibvsCrys);
89 
90  auto CalibEntriesvsCrys = new TH1F("CalibEntriesvsCrys", "Entries in calib vs crys histograms;crystal ID;Entries per crystal", 8736,
91  0, 8736);
92  registerObject<TH1F>("CalibEntriesvsCrys", CalibEntriesvsCrys);
93 
95  auto RawDigitAmpvsCrys = new TH2F("RawDigitAmpvsCrys", "Digit Amplitude vs crystal ID;crystal ID;Amplitude", 8736, 0, 8736, 250, 0,
96  25000);
97  registerObject<TH2F>("RawDigitAmpvsCrys", RawDigitAmpvsCrys);
98 
99  auto RawDigitTimevsCrys = new TH2F("RawDigitTimevsCrys", "Digit Time vs crystal ID;crystal ID;Time", 8736, 0, 8736, 200, -2000,
100  2000);
101  registerObject<TH2F>("RawDigitTimevsCrys", RawDigitTimevsCrys);
102 
103 
104  //------------------------------------------------------------------------
106  myNeighbours4 = new ECLNeighbours("NC", 1);
107  myNeighbours8 = new ECLNeighbours("N", 1);
108 
109 
110  //------------------------------------------------------------------------
112  B2INFO("Input parameters to eclMuMuECollector:");
113  B2INFO("minPairMass: " << m_minPairMass);
114  B2INFO("minTrackLength: " << m_minTrackLength);
115  B2INFO("MaxNeighbourE: " << m_MaxNeighbourE);
116  B2INFO("thetaLabMinDeg: " << m_thetaLabMinDeg);
117  B2INFO("thetaLabMaxDeg: " << m_thetaLabMaxDeg);
118  if (m_thetaLabMinDeg < 1.) {
119  cotThetaLabMax = 9999.;
120  } else if (m_thetaLabMinDeg > 179.) {
121  cotThetaLabMax = -9999.;
122  } else {
123  double thetaLabMin = m_thetaLabMinDeg * TMath::Pi() / 180.;
124  cotThetaLabMax = 1. / tan(thetaLabMin);
125  }
126  if (m_thetaLabMaxDeg < 1.) {
127  cotThetaLabMin = 9999.;
128  } else if (m_thetaLabMaxDeg > 179.) {
129  cotThetaLabMin = -9999.;
130  } else {
131  double thetaLabMax = m_thetaLabMaxDeg * TMath::Pi() / 180.;
132  cotThetaLabMin = 1. / tan(thetaLabMax);
133  }
134  B2INFO("cotThetaLabMin: " << cotThetaLabMin);
135  B2INFO("cotThetaLabMax: " << cotThetaLabMax);
136  B2INFO("measureTrueEnergy: " << m_measureTrueEnergy);
137  B2INFO("requireL1: " << m_requireL1);
138 
140  EperCrys.resize(8736);
141 
144  if (m_ECLExpMuMuE.hasChanged()) {ExpMuMuE = m_ECLExpMuMuE->getCalibVector();}
145  if (m_ElectronicsCalib.hasChanged()) {ElectronicsCalib = m_ElectronicsCalib->getCalibVector();}
146  if (m_MuMuECalib.hasChanged()) {MuMuECalib = m_MuMuECalib->getCalibVector();}
147 
149  for (int ic = 1; ic < 9000; ic += 1000) {
150  B2INFO("DB constants for cellID=" << ic << ": ExpMuMuE = " << ExpMuMuE[ic - 1] << " ElectronicsCalib = " << ElectronicsCalib[ic - 1]
151  << " MuMuECalib = " << MuMuECalib[ic - 1]);
152  }
153 
155  for (int crysID = 0; crysID < 8736; crysID++) {
156  if (ElectronicsCalib[crysID] <= 0) {B2FATAL("eclMuMuECollector: ElectronicsCalib = " << ElectronicsCalib[crysID] << " for crysID = " << crysID);}
157  if (ExpMuMuE[crysID] == 0) {B2FATAL("eclMuMuECollector: ExpMuMuE = 0 for crysID = " << crysID);}
158  if (MuMuECalib[crysID] == 0) {B2FATAL("eclMuMuECollector: MuMuECalib = 0 for crysID = " << crysID);}
159  }
160 
163  m_trackArray.isRequired();
164  m_eclDigitArray.isRequired();
165 
166 }
167 
168 
171 void eclMuMuECollectorModule::collect()
172 {
173 
175  if (iEvent == 0) {
176  for (int crysID = 0; crysID < 8736; crysID++) {
177  getObjectPtr<TH1F>("ExpEvsCrys")->Fill(crysID + 0.001, ExpMuMuE[crysID]);
178  getObjectPtr<TH1F>("ElecCalibvsCrys")->Fill(crysID + 0.001, ElectronicsCalib[crysID]);
179  getObjectPtr<TH1F>("InitialCalibvsCrys")->Fill(crysID + 0.001, MuMuECalib[crysID]);
180  getObjectPtr<TH1F>("CalibEntriesvsCrys")->Fill(crysID + 0.001);
181  }
182  }
183 
184  if (iEvent % 10000 == 0) {B2INFO("eclMuMuECollector: iEvent = " << iEvent);}
185  iEvent++;
186 
189  bool newConst = false;
190  if (m_ECLExpMuMuE.hasChanged()) {
191  newConst = true;
192  B2INFO("ECLExpMuMuE has changed, exp = " << m_evtMetaData->getExperiment() << " run = " << m_evtMetaData->getRun());
193  ExpMuMuE = m_ECLExpMuMuE->getCalibVector();
194  }
195  if (m_ElectronicsCalib.hasChanged()) {
196  newConst = true;
197  B2INFO("ECLCrystalElectronics has changed, exp = " << m_evtMetaData->getExperiment() << " run = " << m_evtMetaData->getRun());
198  ElectronicsCalib = m_ElectronicsCalib->getCalibVector();
199  }
200  if (m_MuMuECalib.hasChanged()) {
201  newConst = true;
202  B2INFO("ECLCrystalEnergyMuMu has changed, exp = " << m_evtMetaData->getExperiment() << " run = " << m_evtMetaData->getRun());
203  MuMuECalib = m_MuMuECalib->getCalibVector();
204  }
205 
206  if (newConst) {
207  for (int ic = 1; ic < 9000; ic += 1000) {
208  B2INFO("DB constants for cellID=" << ic << ": ExpMuMuE = " << ExpMuMuE[ic - 1] << " ElectronicsCalib = " <<
209  ElectronicsCalib[ic - 1]
210  << " MuMuECalib = " << MuMuECalib[ic - 1]);
211  }
212 
214  for (int crysID = 0; crysID < 8736; crysID++) {
215  if (ElectronicsCalib[crysID] <= 0) {B2FATAL("eclMuMuECollector: ElectronicsCalib = " << ElectronicsCalib[crysID] << " for crysID = " << crysID);}
216  if (ExpMuMuE[crysID] == 0) {B2FATAL("eclMuMuECollector: ExpMuMuE = 0 for crysID = " << crysID);}
217  if (MuMuECalib[crysID] == 0) {B2FATAL("eclMuMuECollector: MuMuECalib = 0 for crysID = " << crysID);}
218  }
219  }
220 
221 
222 
225  if (m_requireL1) {
226  unsigned int L1TriggerResults = m_TRGResults->getTRGSummary(0);
227  if (L1TriggerResults == 0) {return;}
228  }
229 
230  //------------------------------------------------------------------------
232  int nTrack = m_trackArray.getEntries();
233  if (nTrack < 2) {return;}
234 
236  double maxpt[2] = {0., 0.};
237  int iTrack[2] = { -1, -1};
238  for (int it = 0; it < nTrack; it++) {
239  const TrackFitResult* temptrackFit = m_trackArray[it]->getTrackFitResult(Const::ChargedStable(211));
240  if (not temptrackFit) {continue;}
241  int imu = 0;
242  if (temptrackFit->getChargeSign() == 1) {imu = 1; }
243 
244  double temppt = temptrackFit->getTransverseMomentum();
245  double cotThetaLab = temptrackFit->getCotTheta();
246  if (temppt > maxpt[imu] && cotThetaLab > cotThetaLabMin && cotThetaLab < cotThetaLabMax) {
247  maxpt[imu] = temppt;
248  iTrack[imu] = it;
249  }
250  }
251 
253  if (iTrack[0] == -1 || iTrack[1] == -1) { return; }
254 
256  TLorentzVector mu0 = m_trackArray[iTrack[0]]->getTrackFitResult(Const::ChargedStable(211))->get4Momentum();
257  TLorentzVector mu1 = m_trackArray[iTrack[1]]->getTrackFitResult(Const::ChargedStable(211))->get4Momentum();
258  if ((mu0 + mu1).M() < m_minPairMass) { return; }
259 
260  //------------------------------------------------------------------------
262  int extCrysID[2] = { -1, -1};
263  Const::EDetector eclID = Const::EDetector::ECL;
264  for (int imu = 0; imu < 2; imu++) {
265  TVector3 temppos[2] = {};
266  int IDEnter = -99;
267  for (auto& extHit : m_trackArray[iTrack[imu]]->getRelationsTo<ExtHit>()) {
268  int pdgCode = extHit.getPdgCode();
269  Const::EDetector detectorID = extHit.getDetectorID(); // subsystem ID
270  int temp0 = extHit.getCopyID(); // ID within that subsystem; for ecl it is crystal ID
271 
273  if (detectorID == eclID && TMath::Abs(pdgCode) == 13 && extHit.getStatus() == EXT_ENTER) {
274  IDEnter = temp0;
275  temppos[0] = extHit.getPosition();
276  }
277 
279  if (detectorID == eclID && TMath::Abs(pdgCode) == 13 && extHit.getStatus() == EXT_EXIT && temp0 == IDEnter) {
280  temppos[1] = extHit.getPosition();
281 
283  double trackLength = (temppos[1] - temppos[0]).Mag();
284  if (trackLength > m_minTrackLength) {extCrysID[imu] = temp0;}
285  break;
286  }
287  }
288  }
289 
291  if (extCrysID[0] == -1 && extCrysID[1] == -1) { return; }
292 
293  //------------------------------------------------------------------------
295  memset(&EperCrys[0], 0, EperCrys.size()*sizeof EperCrys[0]);
296  for (auto& eclDigit : m_eclDigitArray) {
297  int crysID = eclDigit.getCellId() - 1;
298  getObjectPtr<TH2F>("RawDigitAmpvsCrys")->Fill(crysID + 0.001, eclDigit.getAmp());
299 
301  EperCrys[crysID] = eclDigit.getAmp() * abs(MuMuECalib[crysID]) * ElectronicsCalib[crysID];
302  if (EperCrys[crysID] > 0.01) {
303  getObjectPtr<TH2F>("RawDigitTimevsCrys")->Fill(crysID + 0.001, eclDigit.getTimeFit());
304  }
305  }
306 
308  if (m_measureTrueEnergy) {
309 
310  for (auto& eclCalDigit : m_eclCalDigitArray) {
311  int tempCrysID = eclCalDigit.getCellId() - 1;
312  EperCrys[tempCrysID] = eclCalDigit.getEnergy();
313  }
314  }
315 
316  //------------------------------------------------------------------------
318  for (int imu = 0; imu < 2; imu++) {
319  int crysID = extCrysID[imu];
320  int cellID = crysID + 1;
321  if (crysID > -1) {
322 
323  bool noNeighbourSignal = true;
324  if (cellID >= firstcellIDN4 && crysID <= lastcellIDN4) {
325  for (const auto& tempCellID : myNeighbours4->getNeighbours(cellID)) {
326  int tempCrysID = tempCellID - 1;
327  if (tempCellID != cellID && EperCrys[tempCrysID] > m_MaxNeighbourE) {
328  noNeighbourSignal = false;
329  break;
330  }
331  }
332  } else {
333  for (const auto& tempCellID : myNeighbours8->getNeighbours(cellID)) {
334  int tempCrysID = tempCellID - 1;
335  if (tempCellID != cellID && EperCrys[tempCrysID] > m_MaxNeighbourE) {
336  noNeighbourSignal = false;
337  break;
338  }
339  }
340  }
341 
343  if (noNeighbourSignal) {
344 
346  getObjectPtr<TH2F>("EnVsCrysID")->Fill(crysID + 0.001, EperCrys[crysID] / abs(ExpMuMuE[crysID]));
347  getObjectPtr<TH1F>("ExpEvsCrys")->Fill(crysID + 0.001, ExpMuMuE[crysID]);
348  getObjectPtr<TH1F>("ElecCalibvsCrys")->Fill(crysID + 0.001, ElectronicsCalib[crysID]);
349  getObjectPtr<TH1F>("InitialCalibvsCrys")->Fill(crysID + 0.001, MuMuECalib[crysID]);
350  getObjectPtr<TH1F>("CalibEntriesvsCrys")->Fill(crysID + 0.001);
351  }
352  }
353  }
354 }
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::eclMuMuECollectorModule
Calibration collector module that uses muon pairs to do ECL single crystal energy calibration.
Definition: eclMuMuECollectorModule.h:46
Belle2::Const::EDetector
EDetector
Enum for identifying the detector components (detector and subdetector).
Definition: Const.h:44
Belle2::TrackFitResult
Values of the result of a track fit with a given particle hypothesis.
Definition: TrackFitResult.h:59
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::TrackFitResult::getTransverseMomentum
double getTransverseMomentum() const
Getter for the absolute value of the transverse momentum at the perigee.
Definition: TrackFitResult.h:140
Belle2::Const::ChargedStable
Provides a type-safe way to pass members of the chargedStableSet set.
Definition: Const.h:465