Belle II Software  release-06-02-00
ECLChargedPIDDataAnalysisModule.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 #include <iostream>
10 #include <ecl/modules/eclChargedPIDDataAnalysisExpert/ECLChargedPIDDataAnalysisModule.h>
11 #include <framework/datastore/RelationVector.h>
12 #include <framework/logging/Logger.h>
13 #include <framework/gearbox/Const.h>
14 
15 #include <mdst/dataobjects/MCParticle.h>
16 #include <mdst/dataobjects/Track.h>
17 #include <ecl/dataobjects/ECLShower.h>
18 
19 using namespace std;
20 using namespace Belle2;
21 
22 //-----------------------------------------------------------------
23 // Register the Module
24 //-----------------------------------------------------------------
25 
26 REG_MODULE(ECLChargedPIDDataAnalysis)
27 
28 //-----------------------------------------------------------------
29 // Implementation
30 //-----------------------------------------------------------------
31 
32 ECLChargedPIDDataAnalysisModule::ECLChargedPIDDataAnalysisModule() :
33  Module(),
34  m_rootFilePtr(0),
35  m_writeToRoot(1),
36  m_eclShowers(eclShowerArrayName()),
37 
38  // N1 Hypothesis
39  n1_tree(0),
40  n1_iExperiment(0),
41  n1_iRun(0),
42  n1_iEvent(0),
43 
44  // Shower
45  n1_eclShowerMultip(0),
46  n1_eclShowerEnergy(0),
47  n1_eclShowerTheta(0),
48  n1_eclShowerPhi(0),
49  n1_eclShowerR(0),
50  n1_eclShowerHypothesisId(0),
51  n1_eclShowerAbsZernike40(0),
52  n1_eclShowerAbsZernike51(0),
53 
54  // MC
55  n1_mcMultip(0),
56  n1_mcPdg(0),
57  n1_mcMothPdg(0),
58  n1_mcEnergy(0),
59  n1_mcP(0),
60  n1_mcTheta(0),
61  n1_mcPhi(0),
62 
63  // Tracks
64  n1_trkMultip(0),
65  n1_trkPdg(0),
66  n1_trkCharge(0),
67  n1_trkP(0),
68  n1_trkTheta(0),
69  n1_trkPhi(0),
70 
71  n1_eclEoP(0),
72 
73  // N2 Hypothesis
74  n2_tree(0),
75  n2_iExperiment(0),
76  n2_iRun(0),
77  n2_iEvent(0),
78 
79  // Shower
80  n2_eclShowerMultip(0),
81  n2_eclShowerEnergy(0),
82  n2_eclShowerTheta(0),
83  n2_eclShowerPhi(0),
84  n2_eclShowerR(0),
85  n2_eclShowerHypothesisId(0),
86  n2_eclShowerAbsZernike40(0),
87  n2_eclShowerAbsZernike51(0),
88 
89  // MC
90  n2_mcMultip(0),
91  n2_mcPdg(0),
92  n2_mcMothPdg(0),
93  n2_mcEnergy(0),
94  n2_mcP(0),
95  n2_mcTheta(0),
96  n2_mcPhi(0),
97 
98  // Tracks
99  n2_trkMultip(0),
100  n2_trkPdg(0),
101  n2_trkCharge(0),
102  n2_trkP(0),
103  n2_trkTheta(0),
104  n2_trkPhi(0),
105 
106  n2_eclEoP(0)
107 {
108  // Set module properties
109  setDescription("This module produces an ntuple with ECL-related quantities starting from mdst");
110  addParam("writeToRoot", m_writeToRoot,
111  "set true if you want to save the informations in a root file named by parameter 'rootFileName'",
112  bool(true));
113  addParam("rootFileName", m_rootFileName,
114  "fileName used for root file where info are saved. Will be ignored if parameter 'writeToRoot' is false (standard)",
115  string("eclChargedPID"));
116 }
117 
118 ECLChargedPIDDataAnalysisModule::~ECLChargedPIDDataAnalysisModule()
119 {
120 }
121 
122 void ECLChargedPIDDataAnalysisModule::initialize()
123 {
124  B2INFO("[ECLChargedPIDDataAnalysis Module]: Starting initialization of ECLChargedPIDDataAnalysis Module.");
125 
126  if (m_writeToRoot) {
127  m_rootFilePtr = new TFile(m_rootFileName.c_str(), "RECREATE");
128  } else {
129  m_rootFilePtr = nullptr;
130  }
131  // initialize N1 tree
132  n1_tree = new TTree("n1_tree", "ECL Charged PID tree: N1 Hypothesis");
133 
134  n1_tree->Branch("expNo", &n1_iExperiment, "expNo/I");
135  n1_tree->Branch("runNo", &n1_iRun, "runNo/I");
136  n1_tree->Branch("evtNo", &n1_iEvent, "evtNo/I");
137 
138  // shower
139  n1_tree->Branch("eclShowerMultip", &n1_eclShowerMultip, "eclShowerMultip/I");
140  n1_tree->Branch("eclShowerEnergy", "std::vector<double>", &n1_eclShowerEnergy);
141  n1_tree->Branch("eclShowerTheta", "std::vector<double>", &n1_eclShowerTheta);
142  n1_tree->Branch("eclShowerPhi", "std::vector<double>", &n1_eclShowerPhi);
143  n1_tree->Branch("eclShowerR", "std::vector<double>", &n1_eclShowerR);
144  n1_tree->Branch("eclShowerHypothesisId", "std::vector<int>", &n1_eclShowerHypothesisId);
145  n1_tree->Branch("eclShowerAbsZernike40", "std::vector<double>", &n1_eclShowerAbsZernike40);
146  n1_tree->Branch("eclShowerAbsZernike51", "std::vector<double>", &n1_eclShowerAbsZernike51);
147 
148  // MC particle
149  n1_tree->Branch("mcMultip", &n1_mcMultip, "mcMultip/I");
150  n1_tree->Branch("mcPdg", "std::vector<int>", &n1_mcPdg);
151  n1_tree->Branch("mcMothPdg", "std::vector<int>", &n1_mcMothPdg);
152  n1_tree->Branch("mcEnergy", "std::vector<double>", &n1_mcEnergy);
153  n1_tree->Branch("mcP", "std::vector<double>", &n1_mcP);
154  n1_tree->Branch("mcTheta", "std::vector<double>", &n1_mcTheta);
155  n1_tree->Branch("mcPhi", "std::vector<double>", &n1_mcPhi);
156 
157  // tracks
158  n1_tree->Branch("trkMultip", &n1_trkMultip, "trkMulti/I");
159  n1_tree->Branch("trkPdg", "std::vector<int>", &n1_trkPdg);
160  n1_tree->Branch("trkCharge", "std::vector<int>", &n1_trkCharge);
161  n1_tree->Branch("trkP", "std::vector<double>", &n1_trkP);
162  n1_tree->Branch("trkTheta", "std::vector<double>", &n1_trkTheta);
163  n1_tree->Branch("trkPhi", "std::vector<double>", &n1_trkPhi);
164 
165  n1_tree->Branch("eclEoP", "std::vector<double>", &n1_eclEoP);
166 
167  // initialize N2 tree
168  n2_tree = new TTree("n2_tree", "ECL Charged PID tree: N2 Hypothesis");
169 
170  n2_tree->Branch("expNo", &n2_iExperiment, "expNo/I");
171  n2_tree->Branch("runNo", &n2_iRun, "runNo/I");
172  n2_tree->Branch("evtNo", &n2_iEvent, "evtNo/I");
173 
174  // shower
175  n2_tree->Branch("eclShowerMultip", &n2_eclShowerMultip, "eclShowerMultip/I");
176  n2_tree->Branch("eclShowerEnergy", "std::vector<double>", &n2_eclShowerEnergy);
177  n2_tree->Branch("eclShowerTheta", "std::vector<double>", &n2_eclShowerTheta);
178  n2_tree->Branch("eclShowerPhi", "std::vector<double>", &n2_eclShowerPhi);
179  n2_tree->Branch("eclShowerR", "std::vector<double>", &n2_eclShowerR);
180  n2_tree->Branch("eclShowerHypothesisId", "std::vector<int>", &n2_eclShowerHypothesisId);
181  n2_tree->Branch("eclShowerAbsZernike40", "std::vector<double>", &n2_eclShowerAbsZernike40);
182  n2_tree->Branch("eclShowerAbsZernike51", "std::vector<double>", &n2_eclShowerAbsZernike51);
183 
184  // MC particle
185  n2_tree->Branch("mcMultip", &n2_mcMultip, "mcMultip/I");
186  n2_tree->Branch("mcPdg", "std::vector<int>", &n2_mcPdg);
187  n2_tree->Branch("mcMothPdg", "std::vector<int>", &n2_mcMothPdg);
188  n2_tree->Branch("mcEnergy", "std::vector<double>", &n2_mcEnergy);
189  n2_tree->Branch("mcP", "std::vector<double>", &n2_mcP);
190  n2_tree->Branch("mcTheta", "std::vector<double>", &n2_mcTheta);
191  n2_tree->Branch("mcPhi", "std::vector<double>", &n2_mcPhi);
192 
193  // tracks
194  n2_tree->Branch("trkMultip", &n2_trkMultip, "trkMulti/I");
195  n2_tree->Branch("trkPdg", "std::vector<int>", &n2_trkPdg);
196  n2_tree->Branch("trkCharge", "std::vector<int>", &n2_trkCharge);
197  n2_tree->Branch("trkP", "std::vector<double>", &n2_trkP);
198  n2_tree->Branch("trkTheta", "std::vector<double>", &n2_trkTheta);
199  n2_tree->Branch("trkPhi", "std::vector<double>", &n2_trkPhi);
200 
201  n2_tree->Branch("eclEoP", "std::vector<double>", &n2_eclEoP);
202 
203  B2INFO("[ECLChargedPIDDataAnalysis Module]: Initialization of ECLChargedPIDDataAnalysis Module completed.");
204 }
205 
206 void ECLChargedPIDDataAnalysisModule::beginRun()
207 {
208 }
209 
210 void ECLChargedPIDDataAnalysisModule::event()
211 {
212 
213  B2DEBUG(1, " ++++++++++++++ ECLChargedPIDDataAnalysisModule");
214 
215  // Showers
216  n1_eclShowerMultip = 0;
217  n1_eclShowerEnergy->clear();
218  n1_eclShowerTheta->clear();
219  n1_eclShowerPhi->clear();
220  n1_eclShowerR->clear();
221  n1_eclShowerHypothesisId->clear();
222  n1_eclShowerAbsZernike40->clear();
223  n1_eclShowerAbsZernike51->clear();
224 
225  // MC
226  n1_mcMultip = 0;
227  n1_mcPdg->clear();
228  n1_mcMothPdg->clear();
229  n1_mcEnergy->clear();
230  n1_mcP->clear();
231  n1_mcTheta->clear();
232  n1_mcPhi->clear();
233 
234  // Tracks
235  n1_trkMultip = 0;
236  n1_trkPdg->clear();
237  n1_trkCharge->clear();
238  n1_trkP->clear();
239  n1_trkTheta->clear();
240  n1_trkPhi->clear();
241 
242  n1_eclEoP->clear();
243 
244  // Showers
245  n2_eclShowerMultip = 0;
246  n2_eclShowerEnergy->clear();
247  n2_eclShowerTheta->clear();
248  n2_eclShowerPhi->clear();
249  n2_eclShowerR->clear();
250  n2_eclShowerHypothesisId->clear();
251  n2_eclShowerAbsZernike40->clear();
252  n2_eclShowerAbsZernike51->clear();
253 
254  // MC
255  n2_mcMultip = 0;
256  n2_mcPdg->clear();
257  n2_mcMothPdg->clear();
258  n2_mcEnergy->clear();
259  n2_mcP->clear();
260  n2_mcTheta->clear();
261  n2_mcPhi->clear();
262 
263  // Tracks
264  n2_trkMultip = 0;
265  n2_trkPdg->clear();
266  n2_trkCharge->clear();
267  n2_trkP->clear();
268  n2_trkTheta->clear();
269  n2_trkPhi->clear();
270 
271  n2_eclEoP->clear();
272 
273  if (m_EventMetaData) {
274  n1_iExperiment = m_EventMetaData->getExperiment();
275  n1_iRun = m_EventMetaData->getRun();
276  n1_iEvent = m_EventMetaData->getEvent();
277  n2_iExperiment = m_EventMetaData->getExperiment();
278  n2_iRun = m_EventMetaData->getRun();
279  n2_iEvent = m_EventMetaData->getEvent();
280  } else {
281  n1_iExperiment = -1;
282  n1_iRun = -1;
283  n1_iEvent = -1;
284  n2_iExperiment = -1;
285  n2_iRun = -1;
286  n2_iEvent = -1;
287  }
288 
289  // get the matched MC particle
290  for (const MCParticle& imcpart : m_mcParticles) {
291  if (!imcpart.hasStatus(MCParticle::c_PrimaryParticle)) continue; // only check primaries
292  if (imcpart.hasStatus(MCParticle::c_Initial)) continue; // ignore initial particles
293  if (imcpart.hasStatus(MCParticle::c_IsVirtual)) continue; // ignore virtual particles
294 
295  n1_mcMultip++;
296  n2_mcMultip++;
297 
298  // get mc particle kinematics
299  n1_mcPdg->push_back(imcpart.getPDG());
300  if (imcpart.getMother() != nullptr) n1_mcMothPdg->push_back(imcpart.getMother()->getPDG());
301  else n1_mcMothPdg->push_back(-999);
302  n1_mcEnergy->push_back(imcpart.getEnergy());
303  n1_mcP->push_back(imcpart.getMomentum().Mag());
304  n1_mcTheta->push_back(imcpart.getMomentum().Theta());
305  n1_mcPhi->push_back(imcpart.getMomentum().Phi());
306 
307  n2_mcPdg->push_back(imcpart.getPDG());
308  if (imcpart.getMother() != nullptr) n2_mcMothPdg->push_back(imcpart.getMother()->getPDG());
309  else n2_mcMothPdg->push_back(-999);
310  n2_mcEnergy->push_back(imcpart.getEnergy());
311  n2_mcP->push_back(imcpart.getMomentum().Mag());
312  n2_mcTheta->push_back(imcpart.getMomentum().Theta());
313  n2_mcPhi->push_back(imcpart.getMomentum().Phi());
314 
315  // loop over all matched tracks to find index of max momentum
316  int index = 0;
317  int index_max_mom = -1;
318  double max_mom = -1;
319  for (const auto& itrk : imcpart.getRelationsFrom<Track>()) {
320  // get the track fit results
321  const TrackFitResult* atrkF = itrk.getTrackFitResult(Const::pion);
322  if (atrkF == nullptr) continue; //go to next track if no fit result
323  if (atrkF->getMomentum().Mag() > max_mom) {
324  max_mom = atrkF->getMomentum().Mag();
325  index_max_mom = index;
326  }
327  index++;
328  }
329  if (index_max_mom == -1) continue; // go to next mc part if no track found
330 
331  // get the track w/ max momentum
332  const auto itrack = imcpart.getRelationsFrom<Track>()[index_max_mom];
333  // get the track fit results
334  const TrackFitResult* atrkF = itrack->getTrackFitResult(Const::pion);
335 
336  n1_trkMultip++;
337  n2_trkMultip++;
338 
339  // get trk kinematics
340  n1_trkPdg->push_back(atrkF->getParticleType().getPDGCode());
341  n1_trkCharge->push_back(atrkF->getChargeSign());
342  n1_trkP->push_back(atrkF->getMomentum().Mag());
343  n1_trkTheta->push_back(atrkF->getMomentum().Theta());
344  n1_trkPhi->push_back(atrkF->getMomentum().Phi());
345 
346  n2_trkPdg->push_back(atrkF->getParticleType().getPDGCode());
347  n2_trkCharge->push_back(atrkF->getChargeSign());
348  n2_trkP->push_back(atrkF->getMomentum().Mag());
349  n2_trkTheta->push_back(atrkF->getMomentum().Theta());
350  n2_trkPhi->push_back(atrkF->getMomentum().Phi());
351 
352  // loop over all matched ECLShowers (N1,N2) to find index of max energy
353  int jndex1 = -1;
354  int jndex1_max_e = -1;
355  double max_e1 = -1;
356  for (const auto& i1sh : itrack->getRelationsTo<ECLShower>()) {
357  ++jndex1;
358  // use HypoID 5 (N1 Photon hypothesis)
359  if (i1sh.getHypothesisId() != 5) continue;
360  // look only at showers passing the timing selection
361  if (abs(i1sh.getTime()) > i1sh.getDeltaTime99()) continue;
362  if (i1sh.getEnergy() > max_e1) {
363  max_e1 = i1sh.getEnergy();
364  jndex1_max_e = jndex1;
365  }
366  }
367  int jndex2 = -1;
368  int jndex2_max_e = -1;
369  double max_e2 = -1;
370  for (const auto& i2sh : itrack->getRelationsTo<ECLShower>()) {
371  ++jndex2;
372  // use Hypo ID 6 (N2 neutral hadron hypothesis)
373  if (i2sh.getHypothesisId() != 6) continue;
374  // look only at showers passing the timing selection
375  if (abs(i2sh.getTime()) > i2sh.getDeltaTime99()) continue;
376  if (i2sh.getEnergy() > max_e2) {
377  max_e2 = i2sh.getEnergy();
378  jndex2_max_e = jndex2;
379  }
380  }
381 
382  // get the N1, N2 shower w/ max energy
383  if (jndex1_max_e != -1) {
384  const auto i1shower = itrack->getRelationsTo<ECLShower>()[jndex1_max_e];
385  // get shower kinematics
386  n1_eclShowerEnergy->push_back(i1shower->getEnergy());
387  n1_eclShowerTheta->push_back(i1shower->getTheta());
388  n1_eclShowerPhi->push_back(i1shower->getPhi());
389  n1_eclShowerR->push_back(i1shower->getR());
390  n1_eclShowerHypothesisId->push_back(i1shower->getHypothesisId());
391  // get shower Zernike moments
392  n1_eclShowerAbsZernike40->push_back(i1shower->getAbsZernike40());
393  n1_eclShowerAbsZernike51->push_back(i1shower->getAbsZernike51());
394  // get E/p
395  n1_eclEoP->push_back((i1shower->getEnergy()) / (atrkF->getMomentum().Mag()));
396  n1_eclShowerMultip++;
397  }
398  if (jndex2_max_e != -1) {
399  const auto i2shower = itrack->getRelationsTo<ECLShower>()[jndex2_max_e];
400  // get shower kinematics
401  n2_eclShowerEnergy->push_back(i2shower->getEnergy());
402  n2_eclShowerTheta->push_back(i2shower->getTheta());
403  n2_eclShowerPhi->push_back(i2shower->getPhi());
404  n2_eclShowerR->push_back(i2shower->getR());
405  n2_eclShowerHypothesisId->push_back(i2shower->getHypothesisId());
406  // get shower Zernike moments
407  n2_eclShowerAbsZernike40->push_back(i2shower->getAbsZernike40());
408  n2_eclShowerAbsZernike51->push_back(i2shower->getAbsZernike51());
409  // get E/p
410  n2_eclEoP->push_back((i2shower->getEnergy()) / (atrkF->getMomentum().Mag()));
411  n2_eclShowerMultip++;
412  }
413  }
414 
415  n1_tree->Fill();
416  n2_tree->Fill();
417 
418 }
419 
420 void ECLChargedPIDDataAnalysisModule::endRun()
421 {
422 }
423 
424 void ECLChargedPIDDataAnalysisModule::terminate()
425 {
426  if (m_rootFilePtr != nullptr) {
427  m_rootFilePtr->cd();
428  n1_tree->Write();
429  n2_tree->Write();
430  }
431 
432 }
433 
434 
Belle2::ECLChargedPIDDataAnalysisModule
The ECL Charged PID Data Analysis Module.
Definition: ECLChargedPIDDataAnalysisModule.h:37
Belle2::ECLShower
Class to store ECL Showers.
Definition: ECLShower.h:25
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:32
Belle2::Module
Base class for Modules.
Definition: Module.h:72
Belle2::RelationsInterface::getRelationsTo
RelationVector< TO > getRelationsTo(const std::string &name="", const std::string &namedRelation="") const
Get the relations that point from this object to another store array.
Definition: RelationsObject.h:197
Belle2::TrackFitResult
Values of the result of a track fit with a given particle hypothesis.
Definition: TrackFitResult.h:49
Belle2::TrackFitResult::getEnergy
double getEnergy() const
Getter for the Energy at the closest approach of the track in the r/phi projection.
Definition: TrackFitResult.h:123
Belle2::TrackFitResult::getMomentum
TVector3 getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
Definition: TrackFitResult.h:106
Belle2::Track
Class that bundles various TrackFitResults.
Definition: Track.h:25
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17