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