Belle II Software  release-05-02-19
FillTrackFitNtupleModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Stefano Spataro *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <tracking/modules/trackingPerformanceEvaluation/FillTrackFitNtupleModule.h>
12 
13 #include <mdst/dataobjects/HitPatternCDC.h>
14 #include <mdst/dataobjects/HitPatternVXD.h>
15 
16 #include <map>
17 
18 using namespace Belle2;
19 
20 //-----------------------------------------------------------------
21 // Register the Module
22 //-----------------------------------------------------------------
23 REG_MODULE(FillTrackFitNtuple)
24 
25 FillTrackFitNtupleModule::FillTrackFitNtupleModule() :
26  Module()
27 {
28 
29  setDescription("This module fills a ntuple with tracking variables under different hypotheses");
30 
31  addParam("outputFileName", m_rootFileName, "Name of output root file.",
32  std::string("TrackingPerformanceEvaluation_output.root"));
33  addParam("TracksName", m_TracksName, "Name of Track collection.", std::string(""));
34  addParam("RecoTracksName", m_RecoTracksName, "Name of RecoTrack collection.", std::string("RecoTracks"));
35  addParam("ParticleHypothesis", m_ParticleHypothesis, "Particle Hypothesis used in the track fit.", int(211));
36 
37 }
38 
39 void FillTrackFitNtupleModule::initialize()
40 {
41  // Tracks, RecoTracks needed for this module
42  m_RecoTracks.isRequired();
43  m_Tracks.isRequired();
44 
45  //set the ROOT File
46  m_rootFilePtr = new TFile(m_rootFileName.c_str(), "RECREATE");
47 
48  //now create ntuples
49  TString var_list("evt:run:exp:prod:nhits:ncdc:npxd:nsvd:seed_x:");
50  var_list += ("seed_y:seed_z:seed_px:seed_py:seed_pz:seed_p:seed_pt:seed_theta:seed_phi:seed_charge:");
51  var_list += ("nhits_pi:ncdc_pi:npxd_pi:nsvd_pi:nhits_k:ncdc_k:npxd_k:nsvd_k:nhits_p:ncdc_p:npxd_p:nsvd_p:nhits_d:ncdc_d:npxd_d:nsvd_d:");
52  var_list += ("flag_pi:flag_k:flag_p:flag_d:");
53  var_list += ("trk_x_pi:trk_y_pi:trk_z_pi:trk_px_pi:trk_py_pi:trk_pz_pi:trk_p_pi:trk_pt_pi:trk_theta_pi:trk_phi_pi:");
54  var_list += ("trk_charge_pi:trk_chi2_pi:trk_ndf_pi:trk_pvalue_pi:nfailed_pi:");
55  var_list += ("trk_x_k:trk_y_k:trk_z_k:trk_px_k:trk_py_k:trk_pz_k:trk_p_k:trk_pt_k:trk_theta_k:trk_phi_k:");
56  var_list += ("trk_charge_k:trk_chi2_k:trk_ndf_k:trk_pvalue_k:nfailed_k:");
57  var_list += ("trk_x_p:trk_y_p:trk_z_p:trk_px_p:trk_py_p:trk_pz_p:trk_p_p:trk_pt_p:trk_theta_p:trk_phi_p:");
58  var_list += ("trk_charge_p:trk_chi2_p:trk_ndf_p:trk_pvalue_p:nfailed_p:");
59  var_list += ("trk_x_d:trk_y_d:trk_z_d:trk_px_d:trk_py_d:trk_pz_d:trk_p_d:trk_pt_d:trk_theta_d:trk_phi_d:");
60  var_list += ("trk_charge_d:trk_chi2_d:trk_ndf_d:trk_pvalue_d:nfailed_d:");
61  var_list += ("cdcf_pi:cdcl_pi:svdf_pi:svdl_pi:");
62  var_list += ("cdcf_k:cdcl_k:svdf_k:svdl_k:");
63  var_list += ("cdcf_p:cdcl_p:svdf_p:svdl_p:");
64  var_list += ("cdcf_d:cdcl_d:svdf_d:svdl_d");
65  m_n_MultiParticle = new TNtuple("nMultiParticle", "ntuple for multi hyp particle", var_list);
66 }
67 
68 void FillTrackFitNtupleModule::event()
69 {
70 
71  StoreObjPtr<EventMetaData> eventMetaData("EventMetaData", DataStore::c_Event);
72  Float_t event_num = eventMetaData->getEvent();
73  Float_t event_run = eventMetaData->getRun();
74  Float_t event_exp = eventMetaData->getExperiment();
75  Float_t event_prod = eventMetaData->getProduction();
76 
77  B2DEBUG(29, "+++++ Loop on Tracks");
78  StoreArray<Track> tracks(m_Tracks);
79 
80  for (Track& track : tracks) {
81 
82  RecoTrack* recoTrack = track.getRelationsTo<RecoTrack>()[0];
83  if (recoTrack == nullptr) {
84  // if no recoTrack is associated to Track, we skip the track
85  B2WARNING(" the RecoTrack associated to Track is nullptr!");
86  continue;
87  }
88 
89  Float_t nhits = recoTrack->getNumberOfTrackingHits();
90  Float_t ncdc = recoTrack->getNumberOfCDCHits();
91  Float_t npxd = recoTrack->getNumberOfPXDHits();
92  Float_t nsvd = recoTrack->getNumberOfSVDHits();
93 
94  const Const::ChargedStable pdg_list[4] = {Const::pion, Const::kaon, Const::proton, Const::deuteron}; // loop only on these hypotheses
95 
96  std::map <Const::ChargedStable, Float_t> flag; // is the particle hypothesis existing in the track?
97  std::map <Const::ChargedStable, Float_t> trk_x, trk_y, trk_z, trk_px, trk_py, trk_pz, trk_p, trk_pt, trk_theta, trk_phi;
98  std::map <Const::ChargedStable, Float_t> trk_charge, trk_chi2, trk_ndf, trk_pvalue, trk_nfailed;
99  std::map <Const::ChargedStable, Float_t> nhits_pid, ncdc_pid, nsvd_pid, npxd_pid;
100  std::map <Const::ChargedStable, Float_t> first_cdc, last_cdc, first_svd, last_svd;
101 
102  for (const Const::ChargedStable& pdgIter : pdg_list) {
103  trk_x[pdgIter] = 0.; trk_y[pdgIter] = 0.; trk_z[pdgIter] = 0.;
104  trk_px[pdgIter] = 0.; trk_py[pdgIter] = 0.; trk_pz[pdgIter] = 0.;
105  trk_p[pdgIter] = 0.; trk_pt[pdgIter] = 0.; trk_theta[pdgIter] = 0.; trk_phi[pdgIter] = 0.;
106  trk_charge[pdgIter] = 0.;
107  trk_chi2[pdgIter] = 0.; trk_ndf[pdgIter] = 0.; trk_pvalue[pdgIter] = 0.; trk_nfailed[pdgIter] = 0.;
108  nhits_pid[pdgIter] = 0.; ncdc_pid[pdgIter] = 0.; nsvd_pid[pdgIter] = 0.; npxd_pid[pdgIter] = 0.;
109  first_cdc[pdgIter] = -100.; last_cdc [pdgIter] = -100.;
110  first_svd[pdgIter] = -100.; last_svd [pdgIter] = -100.;
111 
112  const TrackFitResult* fitResult = track.getTrackFitResult(pdgIter);
113  if ((fitResult != nullptr) && (fitResult->getParticleType() == pdgIter)) {
114  flag[pdgIter] = kTRUE;
115  } else {
116  flag[pdgIter] = kFALSE;
117  continue;
118  }
119 
120  trk_x[pdgIter] = fitResult->getPosition().X();
121  trk_y[pdgIter] = fitResult->getPosition().Y();
122  trk_z[pdgIter] = fitResult->getPosition().Z();
123  trk_px[pdgIter] = fitResult->getMomentum().X();
124  trk_py[pdgIter] = fitResult->getMomentum().Y();
125  trk_pz[pdgIter] = fitResult->getMomentum().Z();
126  trk_p[pdgIter] = fitResult->getMomentum().Mag();
127  trk_pt[pdgIter] = fitResult->getMomentum().Pt();
128  trk_theta[pdgIter] = fitResult->getMomentum().Theta() * TMath::RadToDeg();
129  trk_phi[pdgIter] = fitResult->getMomentum().Phi() * TMath::RadToDeg();
130  trk_charge[pdgIter] = fitResult->getChargeSign();
131  double chi2 = recoTrack->getTrackFitStatus(recoTrack->getTrackRepresentationForPDG(pdgIter.getPDGCode()))->getChi2();
132  if (isnan(chi2)) chi2 = -10;
133  if (isinf(chi2)) chi2 = -20;
134  trk_chi2[pdgIter] = chi2;
135  trk_ndf[pdgIter] = recoTrack->getTrackFitStatus(recoTrack->getTrackRepresentationForPDG(pdgIter.getPDGCode()))->getNdf();
136  trk_pvalue[pdgIter] = fitResult->getPValue();
137  trk_nfailed[pdgIter] = recoTrack->getTrackFitStatus(recoTrack->getTrackRepresentationForPDG(
138  pdgIter.getPDGCode()))->getNFailedPoints();
139  ncdc_pid[pdgIter] = fitResult->getHitPatternCDC().getNHits();
140  npxd_pid[pdgIter] = fitResult->getHitPatternVXD().getNPXDHits();
141  nsvd_pid[pdgIter] = fitResult->getHitPatternVXD().getNSVDHits();
142  nhits_pid[pdgIter] = ncdc_pid[pdgIter] + npxd_pid[pdgIter] + nsvd_pid[pdgIter];
143  first_cdc[pdgIter] = fitResult->getHitPatternCDC().getFirstLayer();
144  last_cdc[pdgIter] = fitResult->getHitPatternCDC().getLastLayer();
145  first_svd[pdgIter] = fitResult->getHitPatternVXD().getFirstSVDLayer();
146  last_svd[pdgIter] = fitResult->getHitPatternVXD().getLastSVDLayer();
147  }
148 
149  Float_t buffer[] = {event_num, event_run, event_exp, event_prod,
150  nhits, ncdc, npxd, nsvd,
151  (Float_t)recoTrack->getPositionSeed().X(), (Float_t)recoTrack->getPositionSeed().Y(), (Float_t)recoTrack->getPositionSeed().Z(),
152  (Float_t)recoTrack->getMomentumSeed().X(), (Float_t)recoTrack->getMomentumSeed().Y(), (Float_t)recoTrack->getMomentumSeed().Z(), (Float_t)recoTrack->getMomentumSeed().Mag(), (Float_t)recoTrack->getMomentumSeed().Perp(),
153  (Float_t)(recoTrack->getMomentumSeed().Theta()* TMath::RadToDeg()), (Float_t)(recoTrack->getMomentumSeed().Phi()* TMath::RadToDeg()), (Float_t)recoTrack->getChargeSeed(),
154  nhits_pid[Const::pion], ncdc_pid[Const::pion], npxd_pid[Const::pion], nsvd_pid[Const::pion],
155  nhits_pid[Const::kaon], ncdc_pid[Const::kaon], npxd_pid[Const::kaon], nsvd_pid[Const::kaon],
156  nhits_pid[Const::proton], ncdc_pid[Const::proton], npxd_pid[Const::proton], nsvd_pid[Const::proton],
157  nhits_pid[Const::deuteron], ncdc_pid[Const::deuteron], npxd_pid[Const::deuteron], nsvd_pid[Const::deuteron],
158  flag[Const::pion], flag[Const::kaon], flag[Const::proton], flag[Const::deuteron],
159  trk_x[Const::pion], trk_y[Const::pion], trk_z[Const::pion],
160  trk_px[Const::pion], trk_py[Const::pion], trk_pz[Const::pion], trk_p[Const::pion], trk_pt[Const::pion], trk_theta[Const::pion],
161  trk_phi[Const::pion], trk_charge[Const::pion], trk_chi2[Const::pion], trk_ndf[Const::pion], trk_pvalue[Const::pion], trk_nfailed[Const::pion],
162  trk_x[Const::kaon], trk_y[Const::kaon], trk_z[Const::kaon],
163  trk_px[Const::kaon], trk_py[Const::kaon], trk_pz[Const::kaon], trk_p[Const::kaon], trk_pt[Const::kaon], trk_theta[Const::kaon],
164  trk_phi[Const::kaon], trk_charge[Const::kaon], trk_chi2[Const::kaon], trk_ndf[Const::kaon], trk_pvalue[Const::kaon], trk_nfailed[Const::kaon],
165  trk_x[Const::proton], trk_y[Const::proton], trk_z[Const::proton],
166  trk_px[Const::proton], trk_py[Const::proton], trk_pz[Const::proton], trk_p[Const::proton], trk_pt[Const::proton], trk_theta[Const::proton],
167  trk_phi[Const::proton], trk_charge[Const::proton], trk_chi2[Const::proton], trk_ndf[Const::proton], trk_pvalue[Const::proton], trk_nfailed[Const::proton],
168  trk_x[Const::deuteron], trk_y[Const::deuteron], trk_z[Const::deuteron],
169  trk_px[Const::deuteron], trk_py[Const::deuteron], trk_pz[Const::deuteron], trk_p[Const::deuteron], trk_pt[Const::deuteron], trk_theta[Const::deuteron],
170  trk_phi[Const::deuteron], trk_charge[Const::deuteron], trk_chi2[Const::deuteron], trk_ndf[Const::deuteron], trk_pvalue[Const::deuteron], trk_nfailed[Const::deuteron],
171  first_cdc[Const::pion], last_cdc[Const::pion], first_svd[Const::pion], last_svd[Const::pion],
172  first_cdc[Const::kaon], last_cdc[Const::kaon], first_svd[Const::kaon], last_svd[Const::kaon],
173  first_cdc[Const::proton], last_cdc[Const::proton], first_svd[Const::proton], last_svd[Const::proton],
174  first_cdc[Const::deuteron], last_cdc[Const::deuteron], first_svd[Const::deuteron], last_svd[Const::deuteron]
175  };
176  m_n_MultiParticle->Fill(buffer);
177  }
178 
179 
180 
181 }
182 
183 
184 void FillTrackFitNtupleModule::terminate()
185 {
186 
187  if (m_rootFilePtr != nullptr) {
188  m_rootFilePtr->cd();
189  m_n_MultiParticle->Write();
190 
191  m_rootFilePtr->Close();
192  }
193 
194 }
Belle2::RecoTrack::getNumberOfPXDHits
unsigned int getNumberOfPXDHits() const
Return the number of pxd hits.
Definition: RecoTrack.h:423
Belle2::RecoTrack::getTrackFitStatus
const genfit::FitStatus * getTrackFitStatus(const genfit::AbsTrackRep *representation=nullptr) const
Return the track fit status for the given representation or for the cardinal one. You are not allowed...
Definition: RecoTrack.h:537
genfit::FitStatus::getChi2
double getChi2() const
Get chi^2 of the fit.
Definition: FitStatus.h:120
Belle2::HitPatternVXD::getFirstSVDLayer
short getFirstSVDLayer() const
Get the first activated SVD layer index.
Definition: HitPatternVXD.cc:103
Belle2::HitPatternVXD::getNPXDHits
unsigned short getNPXDHits() const
Get total number of hits in the PXD.
Definition: HitPatternVXD.cc:147
Belle2::HitPatternCDC::getFirstLayer
short getFirstLayer() const
Returns the index of the first layer with a hit.
Definition: HitPatternCDC.cc:98
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::RecoTrack::getNumberOfSVDHits
unsigned int getNumberOfSVDHits() const
Return the number of svd hits.
Definition: RecoTrack.h:420
Belle2::RecoTrack::getTrackRepresentationForPDG
genfit::AbsTrackRep * getTrackRepresentationForPDG(int pdgCode)
Return an already created track representation of the given reco track for the PDG.
Definition: RecoTrack.cc:453
Belle2::TrackFitResult::getPValue
double getPValue() const
Getter for Chi2 Probability of the track fit.
Definition: TrackFitResult.h:163
Belle2::FillTrackFitNtupleModule::m_n_MultiParticle
TNtuple * m_n_MultiParticle
Multi particle ntuple.
Definition: FillTrackFitNtupleModule.h:66
Belle2::FillTrackFitNtupleModule
This module takes the Tracks and the RecoTrack input and produce a root file containing an nutple sho...
Definition: FillTrackFitNtupleModule.h:43
Belle2::FillTrackFitNtupleModule::initialize
void initialize() override
Require the store arrays and create the output root file.
Definition: FillTrackFitNtupleModule.cc:39
Belle2::TrackFitResult
Values of the result of a track fit with a given particle hypothesis.
Definition: TrackFitResult.h:59
Belle2::FillTrackFitNtupleModule::m_Tracks
StoreArray< Track > m_Tracks
Track StoreArray.
Definition: FillTrackFitNtupleModule.h:69
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::RecoTrack::getNumberOfTrackingHits
unsigned int getNumberOfTrackingHits() const
Return the number of cdc + svd + pxd hits.
Definition: RecoTrack.h:439
Belle2::TrackFitResult::getHitPatternCDC
HitPatternCDC getHitPatternCDC() const
Getter for the hit pattern in the CDC;.
Definition: TrackFitResult.cc:120
Belle2::Const::kaon
static const ChargedStable kaon
charged kaon particle
Definition: Const.h:536
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::RecoTrack
This is the Reconstruction Event-Data Model Track.
Definition: RecoTrack.h:78
Belle2::Const::pion
static const ChargedStable pion
charged pion particle
Definition: Const.h:535
Belle2::RecoTrack::getPositionSeed
TVector3 getPositionSeed() const
Return the position seed stored in the reco track. ATTENTION: This is not the fitted position.
Definition: RecoTrack.h:477
Belle2::PerformanceEvaluationBaseClass::m_rootFileName
std::string m_rootFileName
root file name
Definition: PerformanceEvaluationBaseClass.h:139
Belle2::TrackFitResult::getParticleType
Const::ParticleType getParticleType() const
Getter for ParticleType of the mass hypothesis of the track fit.
Definition: TrackFitResult.h:154
Belle2::FillTrackFitNtupleModule::m_RecoTracks
StoreArray< RecoTrack > m_RecoTracks
RecoTrack StoreArray.
Definition: FillTrackFitNtupleModule.h:68
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::Const::deuteron
static const ChargedStable deuteron
deuteron particle
Definition: Const.h:538
Belle2::PerformanceEvaluationBaseClass::m_rootFilePtr
TFile * m_rootFilePtr
pointer at root file used for storing histograms
Definition: PerformanceEvaluationBaseClass.h:142
genfit::FitStatus::getNdf
double getNdf() const
Get the degrees of freedom of the fit.
Definition: FitStatus.h:122
Belle2::TrackFitResult::getHitPatternVXD
HitPatternVXD getHitPatternVXD() const
Getter for the hit pattern in the VXD;.
Definition: TrackFitResult.cc:125
Belle2::RecoTrack::getNumberOfCDCHits
unsigned int getNumberOfCDCHits() const
Return the number of cdc hits.
Definition: RecoTrack.h:417
Belle2::RecoTrack::getMomentumSeed
TVector3 getMomentumSeed() const
Return the momentum seed stored in the reco track. ATTENTION: This is not the fitted momentum.
Definition: RecoTrack.h:484
Belle2::Const::proton
static const ChargedStable proton
proton particle
Definition: Const.h:537
Belle2::Const::ChargedStable
Provides a type-safe way to pass members of the chargedStableSet set.
Definition: Const.h:465
Belle2::HitPatternCDC::getLastLayer
short getLastLayer() const
Returns the index of the last layer with a hit.
Definition: HitPatternCDC.cc:106
Belle2::HitPatternCDC::getNHits
unsigned short getNHits() const
Get the total Number of CDC hits in the fit.
Definition: HitPatternCDC.cc:49
Belle2::HitPatternVXD::getNSVDHits
unsigned short getNSVDHits() const
Get total number of hits in the SVD.
Definition: HitPatternVXD.cc:83
Belle2::Track
Class that bundles various TrackFitResults.
Definition: Track.h:35
Belle2::RecoTrack::getChargeSeed
short int getChargeSeed() const
Return the charge seed stored in the reco track. ATTENTION: This is not the fitted charge.
Definition: RecoTrack.h:497
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::FillTrackFitNtupleModule::terminate
void terminate() override
Save output root file with ntuple.
Definition: FillTrackFitNtupleModule.cc:184
Belle2::FillTrackFitNtupleModule::event
void event() override
Loop over Track objects and fill ntuples with tracking parameters.
Definition: FillTrackFitNtupleModule.cc:68
Belle2::DataStore::c_Event
@ c_Event
Different object in each event, all objects/arrays are invalidated after event() function has been ca...
Definition: DataStore.h:61
Belle2::HitPatternVXD::getLastSVDLayer
short getLastSVDLayer() const
Get the last activated SVD layer index.
Definition: HitPatternVXD.cc:112
Belle2::TrackFitResult::getChargeSign
short getChargeSign() const
Return track charge (1 or -1).
Definition: TrackFitResult.h:160