Belle II Software  release-05-02-19
CDCCosmicAnalysisModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: CDC Group *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include "cdc/modules/cdcCosmicAnalysis/CDCCosmicAnalysisModule.h"
12 #include <framework/geometry/BFieldManager.h>
13 #include <framework/gearbox/Const.h>
14 #include <framework/datastore/StoreObjPtr.h>
15 #include <framework/datastore/StoreArray.h>
16 #include <framework/datastore/RelationArray.h>
17 #include <framework/dataobjects/EventMetaData.h>
18 
19 #include <mdst/dataobjects/TrackFitResult.h>
20 #include <mdst/dataobjects/Track.h>
21 #include <tracking/dataobjects/RecoTrack.h>
22 
23 #include <Math/ProbFuncMathCore.h>
24 #include "algorithm"
25 
26 using namespace std;
27 using namespace Belle2;
28 using namespace CDC;
29 
30 
31 //-----------------------------------------------------------------
32 // Register the Module
33 //-----------------------------------------------------------------
34 REG_MODULE(CDCCosmicAnalysis)
35 
36 
37 //-----------------------------------------------------------------
38 // Implementation
39 //-----------------------------------------------------------------
40 
41 CDCCosmicAnalysisModule::CDCCosmicAnalysisModule() : Module()
42 {
43  setDescription("Module for harvesting parameters of the two half (up/down) of a cosmic track for performance study");
44  setPropertyFlags(c_ParallelProcessingCertified); // specify this flag if you need parallel processing
45  addParam("RecoTracksColName", m_recoTrackArrayName, "Name of collectrion hold RecoTracks", std::string(""));
46  addParam("Output", m_outputFileName, "output file name", string("twotracks.root"));
47  addParam("EventT0Extraction", m_eventT0Extraction, "use event t0 extract t0 or not", true);
48  addParam("treeName", m_treeName, "Output tree name", string("tree"));
49  addParam("phi0InRad", m_phi0InRad, "Phi0 in unit of radian, true: rad, false: deg", true);
50  addParam("StoreTrackParErrors", m_storeTrackParErrors,
51  "Store Track Parameter errors (true) or not (false)", false);
52 }
53 
54 CDCCosmicAnalysisModule::~CDCCosmicAnalysisModule()
55 {
56 }
57 void CDCCosmicAnalysisModule::initialize()
58 {
59  // Register histograms (calls back defineHisto)
60  StoreArray<Belle2::Track> storeTrack(m_trackArrayName);
61  StoreArray<RecoTrack> recoTracks(m_recoTrackArrayName);
62  StoreArray<Belle2::TrackFitResult> storeTrackFitResults(m_trackFitResultArrayName);
63  RelationArray relRecoTrackTrack(recoTracks, storeTrack, m_relRecoTrackTrackName);
64 
65  m_recoTrackArrayName = recoTracks.getName();
66  m_trackFitResultArrayName = storeTrackFitResults.getName();
67  m_relRecoTrackTrackName = relRecoTrackTrack.getName();
68 
69  tfile = new TFile(m_outputFileName.c_str(), "RECREATE");
70  // tree = new TTree("treeTrk", "treeTrk");
71  tree = new TTree(m_treeName.c_str(), m_treeName.c_str());
72  tree->Branch("run", &run, "run/I");
73  tree->Branch("evtT0", &evtT0, "evtT0/D");
74  tree->Branch("charge", &charge, "charge/S");
75 
76  tree->Branch("Pval1", &Pval1, "Pval1/D");
77  tree->Branch("ndf1", &ndf1, "ndf1/D");
78  tree->Branch("Pt1", &Pt1, "Pt1/D");
79  tree->Branch("D01", &D01, "D01/D");
80  tree->Branch("Phi01", &Phi01, "Phi01/D");
81  // tree->Branch("Om1", &Om1, "Om1/D");
82  tree->Branch("Z01", &Z01, "Z01/D");
83  tree->Branch("tanLambda1", &tanLambda1, "tanLambda1/D");
84  tree->Branch("posSeed1", "TVector3", &posSeed1);
85  tree->Branch("Omega1", &Omega1, "Omega1/D");
86  tree->Branch("Mom1", "TVector3", &Mom1);
87 
88  tree->Branch("Pval2", &Pval2, "Pval2/D");
89  tree->Branch("ndf2", &ndf2, "ndf2/D");
90  tree->Branch("Pt2", &Pt2, "Pt2/D");
91  tree->Branch("D02", &D02, "D02/D");
92  tree->Branch("Phi02", &Phi02, "Phi02/D");
93  // tree->Branch("Om2", &Om2, "Om2/D");
94  tree->Branch("Z02", &Z02, "Z02/D");
95  tree->Branch("tanLambda2", &tanLambda2, "tanLambda2/D");
96  tree->Branch("posSeed2", "TVector3", &posSeed2);
97  tree->Branch("Omega2", &Omega2, "Omega2/D");
98  tree->Branch("Mom2", "TVector3", &Mom2);
99 
100  if (m_storeTrackParErrors) {
101  tree->Branch("eD01", &eD01, "eD01/D");
102  tree->Branch("ePhi01", &ePhi01, "ePhi01/D");
103  tree->Branch("eOm1", &eOm1, "eOm1/D");
104  tree->Branch("eZ01", &eZ01, "eZ01/D");
105  tree->Branch("etanL1", &etanL1, "etanL1/D");
106 
107  tree->Branch("eD02", &eD02, "eD02/D");
108  tree->Branch("ePhi02", &ePhi02, "ePhi02/D");
109  tree->Branch("eOm2", &eOm2, "eOm2/D");
110  tree->Branch("eZ02", &eZ02, "eZ02/D");
111  tree->Branch("etanL2", &etanL2, "etanL2/D");
112  }
113 
114 }
115 
116 void CDCCosmicAnalysisModule::beginRun()
117 {
118  B2Vector3D pos(0, 0, 0);
119  B2Vector3D bfield = BFieldManager::getFieldInTesla(pos);
120  if (bfield.Z() > 0.5) {
121  m_bField = true;
122  B2INFO("CDCCosmicAnalysis: Magnetic field is ON");
123  } else {
124  m_bField = false;
125  B2INFO("CDCCosmicAnalysis: Magnetic field is OFF");
126  }
127  B2INFO("CDCCosmicAnalysis: BField at (0,0,0) = " << bfield.Mag());
128 }
129 
130 void CDCCosmicAnalysisModule::event()
131 {
132  const StoreArray<Belle2::Track> storeTrack(m_trackArrayName);
133  const StoreArray<Belle2::TrackFitResult> storeTrackFitResults(m_trackFitResultArrayName);
134  const StoreArray<Belle2::RecoTrack> recoTracks(m_recoTrackArrayName);
135  const RelationArray relTrackTrack(recoTracks, storeTrack, m_relRecoTrackTrackName);
136  const StoreObjPtr<EventMetaData> eventMetaData;
137 
138  if (eventMetaData)
139  run = eventMetaData->getRun();
140 
141  evtT0 = 0;
142  if (m_eventT0Extraction) {
143  // event with is fail to extract event-t0 will be excluded
144  if (m_eventTimeStoreObject.isValid() && m_eventTimeStoreObject->hasEventT0()) {
145  evtT0 = m_eventTimeStoreObject->getEventT0();
146  } else {
147  return;
148  }
149  }
150 
151  // Loop over Tracks
152  int nTr = storeTrack.getEntries();
153 
154  short charge2 = 0;
155  short charge1 = 0;
156  bool up(false), down(false);
157  for (int i = 0; i < nTr; ++i) {
158  const Belle2::Track* b2track = storeTrack[i];
159  const Belle2::TrackFitResult* fitresult;
160  fitresult = b2track->getTrackFitResult(Const::muon);
161  if (!fitresult) {
162  B2WARNING("There is no track fit result for muon hypothesis; try with the closest mass hypothesis...");
163  fitresult = b2track->getTrackFitResultWithClosestMass(Const::muon);
164  if (!fitresult) {
165  B2WARNING("There is also no track fit reslt for the other mass hypothesis");
166  continue;
167  }
168  }
169  const Belle2::RecoTrack* recoTrack = b2track->getRelatedTo<Belle2::RecoTrack>(m_recoTrackArrayName);
170  if (!recoTrack) {
171  B2WARNING("Can not access RecoTrack of this Belle2::Track");
172  continue;
173  }
174 
175  TVector3 posSeed = recoTrack->getPositionSeed();
176  const genfit::FitStatus* fs = recoTrack->getTrackFitStatus();
177 
178  double ndf = fs->getNdf();
179  if (!m_bField) // in case of no magnetic field, #track par=4 instead of 5.
180  ndf += 1;
181 
182  double Chi2 = fs->getChi2();
183  double TrPval = std::max(0., ROOT::Math::chisquared_cdf_c(Chi2, ndf));
184  double Phi0 = fitresult->getPhi0();
185  if (m_phi0InRad == false) { // unit of degrees.
186  Phi0 *= 180 / M_PI;
187  }
188 
189  /*** Two track case.***/
190  if ((posSeed.Y() > 0 && !up) ||
191  (up && posSeed.Y()*posSeed1.Y() > 0 && ndf > ndf1)) {
192  charge1 = fitresult->getChargeSign();
193  posSeed1 = posSeed;
194  D01 = fitresult->getD0();
195  eD01 = sqrt(fitresult->getCovariance5()[0][0]);
196  Phi01 = Phi0;
197  ePhi01 = sqrt(fitresult->getCovariance5()[1][1]);
198  if (m_phi0InRad == false) { // unit of degrees.
199  ePhi01 *= 180 / M_PI;
200  }
201 
202  Omega1 = fitresult->getOmega();
203  Mom1 = fitresult->getMomentum();
204  eOm1 = sqrt(fitresult->getCovariance5()[2][2]);
205  Z01 = fitresult->getZ0();
206  eZ01 = sqrt(fitresult->getCovariance5()[3][3]);
207  tanLambda1 = fitresult->getTanLambda();
208  etanL1 = sqrt(fitresult->getCovariance5()[4][4]);
209  Pt1 = fitresult->getTransverseMomentum();
210  ndf1 = ndf;
211  Pval1 = TrPval;
212  up = true;
213  }
214 
215  if ((posSeed.Y() < 0 && !down) ||
216  (down && posSeed.Y()*posSeed2.Y() > 0 && ndf > ndf2)) {
217  charge2 = fitresult->getChargeSign();
218  posSeed2 = posSeed;
219  D02 = fitresult->getD0();
220  eD02 = sqrt(fitresult->getCovariance5()[0][0]);
221  Phi02 = Phi0;
222  ePhi02 = sqrt(fitresult->getCovariance5()[1][1]);
223  if (m_phi0InRad == false) { // unit of degrees.
224  ePhi02 *= 180 / M_PI;
225  }
226  Omega2 = fitresult->getOmega();
227  Mom2 = fitresult->getMomentum();
228  eOm2 = sqrt(fitresult->getCovariance5()[2][2]);
229  Z02 = fitresult->getZ0();
230  eZ02 = sqrt(fitresult->getCovariance5()[3][3]);
231  tanLambda2 = fitresult->getTanLambda();
232  etanL2 = sqrt(fitresult->getCovariance5()[4][4]);
233  Pt2 = fitresult->getTransverseMomentum();
234  ndf2 = ndf;
235  Pval2 = TrPval;
236  down = true;
237  }
238  }
239  if (m_bField && charge1 * charge2 == 0) return;
240  if (charge1 * charge2 >= 0 && up && down) {
241  charge = charge1;
242  tree->Fill();
243  }
244 
245 }
246 
247 void CDCCosmicAnalysisModule::endRun()
248 {
249 }
250 
251 void CDCCosmicAnalysisModule::terminate()
252 {
253  tfile->cd();
254  tree->Write();
255  tfile->Close();
256 }
257 
258 
Belle2::CDC::CDCCosmicAnalysisModule
Analysis module for CDC CR data.
Definition: CDCCosmicAnalysisModule.h:46
Belle2::RelationArray
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:72
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::Track::getTrackFitResult
const TrackFitResult * getTrackFitResult(const Const::ChargedStable &chargedStable) const
Access to TrackFitResults.
Definition: Track.cc:19
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::getOmega
double getOmega() const
Getter for omega.
Definition: TrackFitResult.h:194
Belle2::StoreAccessorBase::getName
const std::string & getName() const
Return name under which the object is saved in the DataStore.
Definition: StoreAccessorBase.h:130
Belle2::RelationsInterface::getRelatedTo
TO * getRelatedTo(const std::string &name="", const std::string &namedRelation="") const
Get the object to which this object has a relation.
Definition: RelationsObject.h:250
Belle2::TrackFitResult::getTanLambda
double getTanLambda() const
Getter for tanLambda.
Definition: TrackFitResult.h:206
Belle2::B2Vector3::Z
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:434
Belle2::TrackFitResult
Values of the result of a track fit with a given particle hypothesis.
Definition: TrackFitResult.h:59
Belle2::TrackFitResult::getZ0
double getZ0() const
Getter for z0.
Definition: TrackFitResult.h:200
Belle2::B2Vector3< double >
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::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::Track::getTrackFitResultWithClosestMass
const TrackFitResult * getTrackFitResultWithClosestMass(const Const::ChargedStable &requestedType) const
Return the track fit for a fit hypothesis with the closest mass.
Definition: Track.cc:70
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::TrackFitResult::getTransverseMomentum
double getTransverseMomentum() const
Getter for the absolute value of the transverse momentum at the perigee.
Definition: TrackFitResult.h:140
Belle2::TrackFitResult::getPhi0
double getPhi0() const
Getter for phi0.
Definition: TrackFitResult.h:184
genfit::FitStatus::getNdf
double getNdf() const
Get the degrees of freedom of the fit.
Definition: FitStatus.h:122
Belle2::TrackFitResult::getCovariance5
TMatrixDSym getCovariance5() const
Getter for covariance matrix of perigee parameters in matrix form.
Definition: TrackFitResult.cc:106
Belle2::Track
Class that bundles various TrackFitResults.
Definition: Track.h:35
Belle2::StoreArray< Belle2::Track >
genfit::FitStatus
Class where important numbers and properties of a fit can be stored.
Definition: FitStatus.h:80
Belle2::B2Vector3::Mag
DataType Mag() const
The magnitude (rho in spherical coordinate system).
Definition: B2Vector3.h:158
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226
Belle2::TrackFitResult::getChargeSign
short getChargeSign() const
Return track charge (1 or -1).
Definition: TrackFitResult.h:160
Belle2::TrackFitResult::getD0
double getD0() const
Getter for d0.
Definition: TrackFitResult.h:178