Belle II Software  release-05-02-19
TOPNtupleModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Marko Staric *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 // Own include
12 #include <top/modules/TOPNtuple/TOPNtupleModule.h>
13 
14 #include <top/geometry/TOPGeometryPar.h>
15 
16 // framework - DataStore
17 #include <framework/datastore/StoreArray.h>
18 #include <framework/datastore/StoreObjPtr.h>
19 
20 // framework aux
21 #include <framework/gearbox/Const.h>
22 #include <framework/logging/Logger.h>
23 
24 // DataStore classes
25 #include <framework/dataobjects/EventMetaData.h>
26 #include <mdst/dataobjects/Track.h>
27 #include <mdst/dataobjects/MCParticle.h>
28 #include <top/dataobjects/TOPLikelihood.h>
29 #include <top/dataobjects/TOPBarHit.h>
30 #include <tracking/dataobjects/ExtHit.h>
31 
32 
33 using namespace std;
34 
35 namespace Belle2 {
41  using namespace TOP;
42 
43  //-----------------------------------------------------------------
44  // Register module
45  //-----------------------------------------------------------------
46 
47  REG_MODULE(TOPNtuple)
48 
49  //-----------------------------------------------------------------
50  // Implementation
51  //-----------------------------------------------------------------
52 
53  TOPNtupleModule::TOPNtupleModule() : Module()
54  {
55  // set module description
56  setDescription("Writes ntuple of TOPLikelihoods with tracking info into a root file");
57 
58  // Add parameters
59  addParam("outputFileName", m_outputFileName, "Output file name",
60  string("TOPNtuple.root"));
61 
62  }
63 
64  TOPNtupleModule::~TOPNtupleModule()
65  {
66  }
67 
68  void TOPNtupleModule::initialize()
69  {
70 
71  m_file = TFile::Open(m_outputFileName.c_str(), "RECREATE");
72  if (m_file->IsZombie()) {
73  B2FATAL("Couldn't open file '" << m_outputFileName << "' for writing!");
74  return;
75  }
76 
77  m_tree = new TTree("top", "TOP validation ntuple");
78 
79  m_tree->Branch("evt", &m_top.evt, "evt/I");
80  m_tree->Branch("run", &m_top.run, "run/I");
81 
82  m_tree->Branch("p", &m_top.p, "p/F");
83  m_tree->Branch("cth", &m_top.cth, "cth/F");
84  m_tree->Branch("phi", &m_top.phi, "phi/F");
85  m_tree->Branch("pValue", &m_top.pValue, "pValue/F");
86 
87  m_tree->Branch("PDG", &m_top.PDG, "PDG/I");
88  m_tree->Branch("motherPDG", &m_top.motherPDG, "motherPDG/I");
89  m_tree->Branch("primary", &m_top.primary, "primary/S");
90  m_tree->Branch("seen", &m_top.seen, "seen/S");
91  m_tree->Branch("rhoProd", &m_top.rhoProd, "rhoProd/F");
92  m_tree->Branch("zProd", &m_top.zProd, "zProd/F");
93  m_tree->Branch("phiProd", &m_top.phiProd, "phiProd/F");
94  m_tree->Branch("rhoDec", &m_top.rhoDec, "rhoDec/F");
95  m_tree->Branch("zDec", &m_top.zDec, "zDec/F");
96  m_tree->Branch("phiDec", &m_top.phiDec, "phiDec/F");
97 
98  m_tree->Branch("numPhot", &m_top.numPhot, "numPhot/I");
99  m_tree->Branch("numBkg", &m_top.numBkg, "numBkg/F");
100  m_tree->Branch("phot", &m_top.phot, "e/F:mu:pi:K:p:d");
101  m_tree->Branch("logL", &m_top.logL, "e/F:mu:pi:K:p:d");
102 
103  m_tree->Branch("extHit", &m_top.extHit, "moduleID/I:PDG:x/F:y:z:p:theta:phi:time");
104  m_tree->Branch("barHit", &m_top.barHit, "moduleID/I:PDG:x/F:y:z:p:theta:phi:time");
105 
106  StoreArray<Track> tracks;
107  tracks.isRequired();
108  StoreArray<ExtHit> extHits;
109  extHits.isRequired();
110  StoreArray<TOPLikelihood> likelihoods;
111  likelihoods.isRequired();
112  StoreArray<MCParticle> mcParticles;
113  mcParticles.isOptional();
114  StoreArray<TOPBarHit> barHits;
115  barHits.isOptional();
116 
117  }
118 
119  void TOPNtupleModule::beginRun()
120  {
121  }
122 
123  void TOPNtupleModule::event()
124  {
125 
126  StoreObjPtr<EventMetaData> evtMetaData;
127  StoreArray<Track> tracks;
128 
129  const auto* geo = TOPGeometryPar::Instance()->getGeometry();
130 
131  for (const auto& track : tracks) {
132  const auto* trackFit = track.getTrackFitResultWithClosestMass(Const::pion);
133  if (!trackFit) continue;
134  const TOPLikelihood* top = track.getRelated<TOPLikelihood>();
135  if (!top) continue;
136  if (top->getFlag() != 1) continue;
137 
138  const ExtHit* extHit = top->getRelated<ExtHit>();
139  const TOPBarHit* barHit = top->getRelated<TOPBarHit>();
140  const MCParticle* mcParticle = track.getRelated<MCParticle>();
141  const MCParticle* mother = 0;
142  if (mcParticle) mother = mcParticle->getMother();
143 
144  m_top.clear();
145 
146  m_top.evt = evtMetaData->getEvent();
147  m_top.run = evtMetaData->getRun();
148 
149  TVector3 mom = trackFit->getMomentum();
150  m_top.p = mom.Mag();
151  m_top.cth = mom.CosTheta();
152  m_top.phi = mom.Phi();
153  m_top.pValue = trackFit->getPValue();
154 
155  if (mcParticle) {
156  m_top.PDG = mcParticle->getPDG();
157  if (mother) m_top.motherPDG = mother->getPDG();
158  m_top.primary = mcParticle->getStatus(MCParticle::c_PrimaryParticle);
159  m_top.seen = mcParticle->hasSeenInDetector(Const::TOP);
160  TVector3 prodVertex = mcParticle->getProductionVertex();
161  m_top.rhoProd = prodVertex.Perp();
162  m_top.zProd = prodVertex.Z();
163  m_top.phiProd = prodVertex.Phi();
164  TVector3 decVertex = mcParticle->getDecayVertex();
165  m_top.rhoDec = decVertex.Perp();
166  m_top.zDec = decVertex.Z();
167  m_top.phiDec = decVertex.Phi();
168  }
169 
170  m_top.numPhot = top->getNphot();
171  m_top.numBkg = top->getEstBkg();
172  m_top.phot.e = top->getEstPhot(Const::electron);
173  m_top.phot.mu = top->getEstPhot(Const::muon);
174  m_top.phot.pi = top->getEstPhot(Const::pion);
175  m_top.phot.K = top->getEstPhot(Const::kaon);
176  m_top.phot.p = top->getEstPhot(Const::proton);
177  m_top.phot.d = top->getEstPhot(Const::deuteron);
178 
179  m_top.logL.e = top->getLogL(Const::electron);
180  m_top.logL.mu = top->getLogL(Const::muon);
181  m_top.logL.pi = top->getLogL(Const::pion);
182  m_top.logL.K = top->getLogL(Const::kaon);
183  m_top.logL.p = top->getLogL(Const::proton);
184  m_top.logL.d = top->getLogL(Const::deuteron);
185 
186  if (extHit) {
187  int moduleID = extHit->getCopyID();
188  TVector3 position = extHit->getPosition();
189  TVector3 momentum = extHit->getMomentum();
190  if (geo->isModuleIDValid(moduleID)) {
191  const auto& module = geo->getModule(moduleID);
192  position = module.pointToLocal(position);
193  momentum = module.momentumToLocal(momentum);
194  }
195  m_top.extHit.moduleID = moduleID;
196  m_top.extHit.PDG = extHit->getPdgCode();
197  m_top.extHit.x = position.X();
198  m_top.extHit.y = position.Y();
199  m_top.extHit.z = position.Z();
200  m_top.extHit.p = momentum.Mag();
201  m_top.extHit.theta = momentum.Theta();
202  m_top.extHit.phi = momentum.Phi();
203  m_top.extHit.time = extHit->getTOF();
204  }
205 
206  if (barHit) {
207  int moduleID = barHit->getModuleID();
208  TVector3 position = barHit->getPosition();
209  TVector3 momentum = barHit->getMomentum();
210  if (geo->isModuleIDValid(moduleID)) {
211  const auto& module = geo->getModule(moduleID);
212  position = module.pointToLocal(position);
213  momentum = module.momentumToLocal(momentum);
214  }
215  m_top.barHit.moduleID = moduleID;
216  m_top.barHit.PDG = barHit->getPDG();
217  m_top.barHit.x = position.X();
218  m_top.barHit.y = position.Y();
219  m_top.barHit.z = position.Z();
220  m_top.barHit.p = momentum.Mag();
221  m_top.barHit.theta = momentum.Theta();
222  m_top.barHit.phi = momentum.Phi();
223  m_top.barHit.time = barHit->getTime();
224  }
225 
226  m_tree->Fill();
227  }
228 
229  }
230 
231 
232  void TOPNtupleModule::endRun()
233  {
234  }
235 
236  void TOPNtupleModule::terminate()
237  {
238  m_file->cd();
239  m_tree->Write();
240  m_file->Close();
241  }
242 
243 
245 } // end Belle2 namespace
246 
Belle2::ExtHit::getPdgCode
int getPdgCode() const
Get PDG code of this extrapolation's hypothesis.
Definition: ExtHit.h:126
Belle2::ExtHit::getPosition
TVector3 getPosition() const
Get position of this extrapolation hit.
Definition: ExtHit.h:153
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::ExtHit::getTOF
double getTOF() const
Get time of flight from the point of closest approach near the origin to this hit.
Definition: ExtHit.h:145
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::ExtHit
Store one Ext hit as a ROOT object.
Definition: ExtHit.h:40
Belle2::TOPBarHit
Class to store track parameters of incoming MC particles relation to MCParticle filled in top/simulat...
Definition: TOPBarHit.h:36
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::ExtHit::getCopyID
int getCopyID() const
Get detector-element ID of sensitive element within detector.
Definition: ExtHit.h:134
Belle2::TOPLikelihood
Class to store TOP log likelihoods (output of TOPReconstructor).
Definition: TOPLikelihood.h:37
Belle2::ExtHit::getMomentum
TVector3 getMomentum() const
Get momentum at this extrapolation hit.
Definition: ExtHit.h:157
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:43
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::TOPNtupleModule
Module to write out a ntuple summarizing TOP reconstruction output.
Definition: TOPNtupleModule.h:38