Belle II Software  release-05-02-19
ARICHNtupleModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Luka Santelj, Dino Tahirovic *
7  * *
8  * The purpose of this module is to test the reconstruction of *
9  * the particles with ARICH. *
10  * *
11  * This software is provided "as is" without any warranty. *
12  **************************************************************************/
13 
14 // Own include
15 #include <arich/modules/arichNtuple/ARICHNtupleModule.h>
16 
17 // Hit classes
18 #include <framework/dataobjects/EventMetaData.h>
19 #include <arich/dataobjects/ARICHLikelihood.h>
20 #include <arich/dataobjects/ARICHAeroHit.h>
21 #include <arich/dataobjects/ARICHTrack.h>
22 #include <arich/dataobjects/ARICHPhoton.h>
23 #include <arich/dataobjects/ARICHHit.h>
24 #include <arich/dataobjects/ARICHInfo.h>
25 
26 // framework - DataStore
27 #include <framework/datastore/StoreObjPtr.h>
28 
29 #include <mdst/dataobjects/HitPatternCDC.h>
30 
31 #ifdef ALIGNMENT_USING_BHABHA
32 #include <mdst/dataobjects/ECLCluster.h>
33 #endif
34 
35 // framework aux
36 #include <framework/gearbox/Const.h>
37 #include <framework/logging/Logger.h>
38 
39 // ROOT
40 #include <TVector3.h>
41 #include <vector>
42 
43 using namespace std;
44 
45 namespace Belle2 {
51  //-----------------------------------------------------------------
52  // Register module
53  //-----------------------------------------------------------------
54 
55  REG_MODULE(ARICHNtuple)
56 
57  //-----------------------------------------------------------------
58  // Implementation
59  //-----------------------------------------------------------------
60 
61  ARICHNtupleModule::ARICHNtupleModule() : Module(),
62  m_file(0),
63  m_tree(0)
64  {
65  // set module description (e.g. insert text)
66  setDescription("The module saves variables needed for performance analysis, such as postion and momentum of the hit, likelihoods for hypotheses and number of photons.");
67 
68  // Add parameters
69  addParam("outputFile", m_outputFile, "ROOT output file name", string("ARICHNtuple.root"));
70  }
71 
72  ARICHNtupleModule::~ARICHNtupleModule()
73  {
74  }
75 
76  void ARICHNtupleModule::initialize()
77  {
78 
79  if (m_file) delete m_file;
80  m_file = new TFile(m_outputFile.c_str(), "RECREATE");
81  if (m_file->IsZombie()) {
82  B2FATAL("Couldn't open file '" << m_outputFile << "' for writing!");
83  return;
84  }
85 
86  m_tree = new TTree("arich", "ARICH validation ntuple");
87 
88  m_tree->Branch("evt", &m_arich.evt, "evt/I");
89  m_tree->Branch("run", &m_arich.run, "run/I");
90  m_tree->Branch("exp", &m_arich.exp, "exp/I");
91 
92  m_tree->Branch("charge", &m_arich.charge, "charge/S");
93  m_tree->Branch("pValue", &m_arich.pValue, "pValue/F");
94  m_tree->Branch("d0", &m_arich.z0, "pValue/F");
95  m_tree->Branch("z0", &m_arich.d0, "pValue/F");
96 
97 #ifdef ALIGNMENT_USING_BHABHA
98  m_tree->Branch("eop", &m_arich.eop, "eop/F");
99  m_tree->Branch("e9e21", &m_arich.e9e21, "e9e21/F");
100  m_tree->Branch("etot", &m_arich.etot, "etot/F");
101 #endif
102 
103  m_tree->Branch("PDG", &m_arich.PDG, "PDG/I");
104  m_tree->Branch("motherPDG", &m_arich.motherPDG, "motherPDG/I");
105  m_tree->Branch("primary", &m_arich.primary, "primary/S");
106  m_tree->Branch("seen", &m_arich.seen, "seen/S");
107  m_tree->Branch("scatter", &m_arich.scatter, "scatter/I");
108 
109  m_tree->Branch("rhoProd", &m_arich.rhoProd, "rhoProd/F");
110  m_tree->Branch("zProd", &m_arich.zProd, "zProd/F");
111  m_tree->Branch("phiProd", &m_arich.phiProd, "phiProd/F");
112  m_tree->Branch("rhoDec", &m_arich.rhoDec, "rhoDec/F");
113  m_tree->Branch("zDec", &m_arich.zDec, "zDec/F");
114  m_tree->Branch("phiDec", &m_arich.phiDec, "phiDec/F");
115  m_tree->Branch("status", &m_arich.status, "status/I");
116 
117  m_tree->Branch("detPhot", &m_arich.detPhot, "detPhot/I");
118  m_tree->Branch("numBkg", &m_arich.numBkg, "e/F:mu:pi:K:p:d");
119  m_tree->Branch("expPhot", &m_arich.expPhot, "e/F:mu:pi:K:p:d");
120  m_tree->Branch("logL", &m_arich.logL, "e/F:mu:pi:K:p:d");
121 
122  m_tree->Branch("recHit", &m_arich.recHit, "PDG/I:x/F:y:z:p:theta:phi");
123  m_tree->Branch("mcHit", &m_arich.mcHit, "PDG/I:x/F:y:z:p:theta:phi");
124  m_tree->Branch("winHit", &m_arich.winHit, "x/F:y");
125  m_tree->Branch("nrec", &m_arich.nRec, "nRec/I");
126  m_tree->Branch("nCDC", &m_arich.nCDC, "nCDC/I");
127  m_tree->Branch("inAcceptance", &m_arich.inAcc, "inAcc/O");
128  m_tree->Branch("photons", "std::vector<Belle2::ARICHPhoton>", &m_arich.photons);
129 
130  // required input
131  m_arichTracks.isRequired();
132  m_arichLikelihoods.isRequired();
133 
134  // optional input
135  m_tracks.isOptional();
136  m_arichMCPs.isOptional();
137  m_arichAeroHits.isOptional();
138  m_arichInfo.isOptional();
139 
140  StoreArray<ARICHHit> arichHits;
141  arichHits.isRequired();
142 
143  }
144 
145  void ARICHNtupleModule::beginRun()
146  {
147  }
148 
149  void ARICHNtupleModule::event()
150  {
151 
152  StoreObjPtr<EventMetaData> evtMetaData;
153  StoreArray<ARICHHit> arichHits;
154 
155  if (!m_arichTracks.isValid()) return;
156 
157 #ifdef ALIGNMENT_USING_BHABHA
158  double E_tot = 0.;
159  StoreArray<ECLCluster> eclClusters;
160  for (const auto& trk : m_tracks) {
161  const ECLCluster* eclTrack = trk.getRelated<ECLCluster>();
162  if (eclTrack) {
163  if (eclTrack->getEnergy() > 0.1) E_tot += eclTrack->getEnergy();
164  }
165  }
166  for (const auto& cluster : eclClusters) {
167  if (!cluster.isNeutral()) continue;
168  if (cluster.getHypothesisId() != 5 && cluster.getHypothesisId() != 1)
169  continue;
170  if (cluster.getEnergy() > 0.1) E_tot += cluster.getEnergy();
171  }
172 #endif
173 
174  for (const auto& arichTrack : m_arichTracks) {
175 
176  const ExtHit* extHit = arichTrack.getRelated<ExtHit>();
177 
178  const Track* mdstTrack = NULL;
179  if (extHit) mdstTrack = extHit->getRelated<Track>();
180 
181  const ARICHAeroHit* aeroHit = arichTrack.getRelated<ARICHAeroHit>();
182 
183  const ARICHLikelihood* lkh = NULL;
184  if (mdstTrack) lkh = mdstTrack->getRelated<ARICHLikelihood>();
185  else lkh = arichTrack.getRelated<ARICHLikelihood>();
186 
187  if (!lkh) continue;
188 
189  m_arich.clear();
190 
191  m_arich.evt = evtMetaData->getEvent();
192  m_arich.run = evtMetaData->getRun();
193  m_arich.exp = evtMetaData->getExperiment();
194 
195  // set hapd window hit if available
196  if (arichTrack.hitsWindow()) {
197  TVector2 winHit = arichTrack.windowHitPosition();
198  m_arich.winHit[0] = winHit.X();
199  m_arich.winHit[1] = winHit.Y();
200  }
201 
202  if (lkh->getFlag() == 1) m_arich.inAcc = 1;
203  m_arich.logL.e = lkh->getLogL(Const::electron);
204  m_arich.logL.mu = lkh->getLogL(Const::muon);
205  m_arich.logL.pi = lkh->getLogL(Const::pion);
206  m_arich.logL.K = lkh->getLogL(Const::kaon);
207  m_arich.logL.p = lkh->getLogL(Const::proton);
208  m_arich.logL.d = lkh->getLogL(Const::deuteron);
209 
210  m_arich.expPhot.e = lkh->getExpPhot(Const::electron);
211  m_arich.expPhot.mu = lkh->getExpPhot(Const::muon);
212  m_arich.expPhot.pi = lkh->getExpPhot(Const::pion);
213  m_arich.expPhot.K = lkh->getExpPhot(Const::kaon);
214  m_arich.expPhot.p = lkh->getExpPhot(Const::proton);
215  m_arich.expPhot.d = lkh->getExpPhot(Const::deuteron);
216 
217  m_arich.detPhot = lkh->getDetPhot();
218 
219 #ifdef ALIGNMENT_USING_BHABHA
220  m_arich.etot = E_tot;
221 #endif
222  m_arich.status = 1;
223 
224  m_arich.trgtype = 0;
225  if (m_arichInfo.getEntries() > 0) {
226  auto arichInfo = *m_arichInfo[0];
227  m_arich.trgtype = arichInfo.gettrgtype();
228  }
229 
230  const MCParticle* particle = 0;
231 
232  if (mdstTrack) {
233  const TrackFitResult* fitResult = mdstTrack->getTrackFitResultWithClosestMass(Const::pion);
234  if (fitResult) {
235  m_arich.pValue = fitResult->getPValue();
236  TVector3 trkPos = fitResult->getPosition();
237  m_arich.charge = fitResult->getChargeSign();
238  m_arich.z0 = trkPos.Z();
239  m_arich.d0 = (trkPos.XYvector()).Mod();
240  m_arich.nCDC = fitResult->getHitPatternCDC().getNHits();
241 #ifdef ALIGNMENT_USING_BHABHA
242  TVector3 trkMom = fitResult->getMomentum();
243  const ECLCluster* eclTrack = mdstTrack->getRelated<ECLCluster>();
244  if (eclTrack) {
245  m_arich.eop = eclTrack->getEnergy() / trkMom.Mag();
246  m_arich.e9e21 = eclTrack->getE9oE21();
247  }
248 #endif
249  }
250  m_arich.status += 10;
251  int fromARICHMCP = 0;
252  particle = mdstTrack->getRelated<MCParticle>();
253  if (!particle) { particle = mdstTrack->getRelated<MCParticle>("arichMCParticles"); fromARICHMCP = 1;}
254  if (particle) {
255  m_arich.PDG = particle->getPDG();
256  m_arich.primary = particle->getStatus(MCParticle::c_PrimaryParticle);
257  m_arich.seen = particle->hasSeenInDetector(Const::ARICH);
258  TVector3 prodVertex = particle->getProductionVertex();
259  m_arich.rhoProd = prodVertex.Perp();
260  m_arich.zProd = prodVertex.Z();
261  m_arich.phiProd = prodVertex.Phi();
262  TVector3 decVertex = particle->getDecayVertex();
263  m_arich.rhoDec = decVertex.Perp();
264  m_arich.zDec = decVertex.Z();
265  m_arich.phiDec = decVertex.Phi();
266 
267  if (!fromARICHMCP) {
268  MCParticle* mother = particle->getMother();
269  if (mother) m_arich.motherPDG = mother->getPDG();
270  std::vector<Belle2::MCParticle*> daughs = particle->getDaughters();
271  for (const auto daugh : daughs) {
272  if (daugh == NULL) continue;
273  if (daugh->getPDG() == particle->getPDG()) m_arich.scatter = 1;
274  }
275  }
276  }
277  }
278 
279 
280  // get reconstructed photons associated with track
281  const std::vector<ARICHPhoton>& photons = arichTrack.getPhotons();
282  m_arich.nRec = photons.size();
283  int nphot = 0;
284  for (auto it = photons.begin(); it != photons.end(); ++it) {
285  ARICHPhoton iph = *it;
286  if (iph.getHitID() < arichHits.getEntries()) {
287  auto ihit = *arichHits[iph.getHitID()];
288  int module = ihit.getModule();
289  int channel = ihit.getChannel();
290  int chid = 1 << 16 | module << 8 | channel;
291  iph.setHitID(chid);
292  }
293  m_arich.photons.push_back(iph);
294  nphot++;
295  }
296 
297  TVector3 recPos = arichTrack.getPosition();
298  m_arich.recHit.x = recPos.X();
299  m_arich.recHit.y = recPos.Y();
300  m_arich.recHit.z = recPos.Z();
301 
302  TVector3 recMom = arichTrack.getDirection() * arichTrack.getMomentum();
303  m_arich.recHit.p = recMom.Mag();
304  m_arich.recHit.theta = recMom.Theta();
305  m_arich.recHit.phi = recMom.Phi();
306 
307  if (aeroHit) {
308  TVector3 truePos = aeroHit->getPosition();
309  m_arich.mcHit.x = truePos.X();
310  m_arich.mcHit.y = truePos.Y();
311  m_arich.mcHit.z = truePos.Z();
312 
313  TVector3 trueMom = aeroHit->getMomentum();
314  m_arich.mcHit.p = trueMom.Mag();
315  m_arich.mcHit.theta = trueMom.Theta();
316  m_arich.mcHit.phi = trueMom.Phi();
317 
318  m_arich.mcHit.PDG = aeroHit->getPDG();
319  m_arich.status += 1000;
320 
321  if (!particle) {
322  particle = aeroHit->getRelated<MCParticle>();
323  if (particle) {
324  m_arich.PDG = particle->getPDG();
325  MCParticle* mother = particle->getMother();
326  if (mother) m_arich.motherPDG = mother->getPDG();
327  m_arich.primary = particle->getStatus(MCParticle::c_PrimaryParticle);
328  m_arich.seen = particle->hasSeenInDetector(Const::ARICH);
329  TVector3 prodVertex = particle->getProductionVertex();
330  m_arich.rhoProd = prodVertex.Perp();
331  m_arich.zProd = prodVertex.Z();
332  m_arich.phiProd = prodVertex.Phi();
333  TVector3 decVertex = particle->getDecayVertex();
334  m_arich.rhoDec = decVertex.Perp();
335  m_arich.zDec = decVertex.Z();
336  m_arich.phiDec = decVertex.Phi();
337 
338  std::vector<Belle2::MCParticle*> daughs = particle->getDaughters();
339  for (const auto daugh : daughs) {
340  if (daugh->getPDG() == particle->getPDG()) m_arich.scatter = 1;
341  }
342  }
343  }
344  }
345  m_tree->Fill();
346  }
347  }
348 
349 
350 
351  void ARICHNtupleModule::endRun()
352  {
353  }
354 
355  void ARICHNtupleModule::terminate()
356  {
357  m_file->cd();
358  m_tree->Write();
359  m_file->Close();
360  }
361 
363 } // end Belle2 namespace
364 
Belle2::ARICHLikelihood
This is a class to store ARICH likelihoods in the datastore.
Definition: ARICHLikelihood.h:38
Belle2::ARICHNtupleModule
ARICH reconstruction efficiency test module.
Definition: ARICHNtupleModule.h:73
Belle2::ECLCluster::getE9oE21
double getE9oE21() const
Return E9/E21 (shower shape variable).
Definition: ECLCluster.h:295
Belle2::ARICHLikelihood::getLogL
float getLogL(const Const::ChargedStable &part) const
Return log likelihood for a given particle.
Definition: ARICHLikelihood.h:83
Belle2::TrackFitResult::getMomentum
TVector3 getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
Definition: TrackFitResult.h:116
Belle2::ECLCluster
ECL cluster data.
Definition: ECLCluster.h:39
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::ARICHLikelihood::getFlag
int getFlag() const
Get reconstruction flag.
Definition: ARICHLikelihood.h:76
Belle2::TrackFitResult::getPValue
double getPValue() const
Getter for Chi2 Probability of the track fit.
Definition: TrackFitResult.h:163
Belle2::ARICHAeroHit
Datastore class that holds information on track parameters at the entrance in aerogel.
Definition: ARICHAeroHit.h:37
Belle2::ARICHLikelihood::getExpPhot
float getExpPhot(const Const::ChargedStable &part) const
Return number of expected photons for a given particle.
Definition: ARICHLikelihood.h:102
Belle2::RelationsInterface::getRelated
T * getRelated(const std::string &name="", const std::string &namedRelation="") const
Get the object to or from which this object has a relation.
Definition: RelationsObject.h:280
Belle2::TrackFitResult
Values of the result of a track fit with a given particle hypothesis.
Definition: TrackFitResult.h:59
Belle2::TrackFitResult::getPosition
TVector3 getPosition() const
Getter for vector of position at closest approach of track in r/phi projection.
Definition: TrackFitResult.h:109
Belle2::TrackFitResult::getHitPatternCDC
HitPatternCDC getHitPatternCDC() const
Getter for the hit pattern in the CDC;.
Definition: TrackFitResult.cc:120
Belle2::ARICHPhoton::setHitID
void setHitID(int id)
Set ID of corresponding ARICHHit.
Definition: ARICHPhoton.h:50
Belle2::ARICHLikelihood::getDetPhot
float getDetPhot() const
Return number of detected photons for a given particle.
Definition: ARICHLikelihood.h:92
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::ECLCluster::getEnergy
double getEnergy(const EHypothesisBit &hypothesis) const
Return Energy (GeV).
Definition: ECLCluster.cc:21
Belle2::ExtHit
Store one Ext hit as a ROOT object.
Definition: ExtHit.h:40
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::ARICHPhoton::getHitID
int getHitID()
Get ID of corresponding ARICHHit.
Definition: ARICHPhoton.h:232
Belle2::HitPatternCDC::getNHits
unsigned short getNHits() const
Get the total Number of CDC hits in the fit.
Definition: HitPatternCDC.cc: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
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::TrackFitResult::getChargeSign
short getChargeSign() const
Return track charge (1 or -1).
Definition: TrackFitResult.h:160
Belle2::ARICHPhoton
Struct for ARICH reconstructed photon (hit related to track) information.
Definition: ARICHPhoton.h:27