Belle II Software development
CDCCosmicAnalysisModule.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include "cdc/modules/cdcCosmicAnalysis/CDCCosmicAnalysisModule.h"
10#include <framework/geometry/BFieldManager.h>
11#include <framework/geometry/B2Vector3.h>
12#include <framework/gearbox/Const.h>
13#include <framework/datastore/RelationArray.h>
14
15#include <Math/ProbFuncMathCore.h>
16
17using namespace std;
18using namespace Belle2;
19using namespace CDC;
20
21
22//-----------------------------------------------------------------
23// Register the Module
24//-----------------------------------------------------------------
25REG_MODULE(CDCCosmicAnalysis);
26
27
28//-----------------------------------------------------------------
29// Implementation
30//-----------------------------------------------------------------
31
33{
34 setDescription("Module for harvesting parameters of the two half (up/down) of a cosmic track for performance study");
35 setPropertyFlags(c_ParallelProcessingCertified); // specify this flag if you need parallel processing
36 addParam("RecoTracksColName", m_recoTrackArrayName, "Name of collectrion hold RecoTracks", std::string(""));
37 addParam("Output", m_outputFileName, "output file name", string("twotracks.root"));
38 addParam("EventT0Extraction", m_eventT0Extraction, "use event t0 extract t0 or not", true);
39 addParam("treeName", m_treeName, "Output tree name", string("tree"));
40 addParam("phi0InRad", m_phi0InRad, "Phi0 in unit of radian, true: rad, false: deg", true);
41 addParam("StoreTrackParErrors", m_storeTrackParErrors,
42 "Store Track Parameter errors (true) or not (false)", false);
43}
44
49{
50 // Register histograms (calls back defineHisto)
51 m_Tracks.isRequired(m_trackArrayName);
55
56 m_relRecoTrackTrackName = relRecoTrackTrack.getName();
57
58 tfile = new TFile(m_outputFileName.c_str(), "RECREATE");
59 // tree = new TTree("treeTrk", "treeTrk");
60 tree = new TTree(m_treeName.c_str(), m_treeName.c_str());
61 tree->Branch("run", &run, "run/I");
62 tree->Branch("evtT0", &evtT0, "evtT0/D");
63 tree->Branch("charge", &charge, "charge/S");
64
65 tree->Branch("Pval1", &Pval1, "Pval1/D");
66 tree->Branch("ndf1", &ndf1, "ndf1/D");
67 tree->Branch("Pt1", &Pt1, "Pt1/D");
68 tree->Branch("D01", &D01, "D01/D");
69 tree->Branch("Phi01", &Phi01, "Phi01/D");
70 // tree->Branch("Om1", &Om1, "Om1/D");
71 tree->Branch("Z01", &Z01, "Z01/D");
72 tree->Branch("tanLambda1", &tanLambda1, "tanLambda1/D");
73 tree->Branch("posSeed1", "TVector3", &posSeed1);
74 tree->Branch("Omega1", &Omega1, "Omega1/D");
75 tree->Branch("Mom1", "TVector3", &Mom1);
76
77 tree->Branch("Pval2", &Pval2, "Pval2/D");
78 tree->Branch("ndf2", &ndf2, "ndf2/D");
79 tree->Branch("Pt2", &Pt2, "Pt2/D");
80 tree->Branch("D02", &D02, "D02/D");
81 tree->Branch("Phi02", &Phi02, "Phi02/D");
82 // tree->Branch("Om2", &Om2, "Om2/D");
83 tree->Branch("Z02", &Z02, "Z02/D");
84 tree->Branch("tanLambda2", &tanLambda2, "tanLambda2/D");
85 tree->Branch("posSeed2", "TVector3", &posSeed2);
86 tree->Branch("Omega2", &Omega2, "Omega2/D");
87 tree->Branch("Mom2", "TVector3", &Mom2);
88
90 tree->Branch("eD01", &eD01, "eD01/D");
91 tree->Branch("ePhi01", &ePhi01, "ePhi01/D");
92 tree->Branch("eOm1", &eOm1, "eOm1/D");
93 tree->Branch("eZ01", &eZ01, "eZ01/D");
94 tree->Branch("etanL1", &etanL1, "etanL1/D");
95
96 tree->Branch("eD02", &eD02, "eD02/D");
97 tree->Branch("ePhi02", &ePhi02, "ePhi02/D");
98 tree->Branch("eOm2", &eOm2, "eOm2/D");
99 tree->Branch("eZ02", &eZ02, "eZ02/D");
100 tree->Branch("etanL2", &etanL2, "etanL2/D");
101 }
102
103}
104
106{
107 ROOT::Math::XYZVector pos(0, 0, 0);
108 ROOT::Math::XYZVector bfield = BFieldManager::getFieldInTesla(pos);
109 if (bfield.Z() > 0.5) {
110 m_bField = true;
111 B2INFO("CDCCosmicAnalysis: Magnetic field is ON");
112 } else {
113 m_bField = false;
114 B2INFO("CDCCosmicAnalysis: Magnetic field is OFF");
115 }
116 B2INFO("CDCCosmicAnalysis: BField at (0,0,0) = " << bfield.R());
117}
118
120{
121 if (m_EventMetaData.isValid())
122 run = m_EventMetaData->getRun();
123
124 evtT0 = 0;
126 // event with is fail to extract event-t0 will be excluded
127 if (m_eventTimeStoreObject.isValid() && m_eventTimeStoreObject->hasEventT0()) {
128 evtT0 = m_eventTimeStoreObject->getEventT0();
129 } else {
130 return;
131 }
132 }
133
134 // Loop over Tracks
135 int nTr = m_Tracks.getEntries();
136
137 short charge2 = 0;
138 short charge1 = 0;
139 bool up(false), down(false);
140 for (int i = 0; i < nTr; ++i) {
141 const Belle2::Track* b2track = m_Tracks[i];
142 const Belle2::TrackFitResult* fitresult;
143 fitresult = b2track->getTrackFitResult(Const::muon);
144 if (!fitresult) {
145 B2WARNING("There is no track fit result for muon hypothesis; try with the closest mass hypothesis...");
146 fitresult = b2track->getTrackFitResultWithClosestMass(Const::muon);
147 if (!fitresult) {
148 B2WARNING("There is also no track fit reslt for the other mass hypothesis");
149 continue;
150 }
151 }
153 if (!recoTrack) {
154 B2WARNING("Can not access RecoTrack of this Belle2::Track");
155 continue;
156 }
157
158 B2Vector3D posSeed = recoTrack->getPositionSeed();
159 const genfit::FitStatus* fs = recoTrack->getTrackFitStatus();
160
161 double ndf = fs->getNdf();
162 if (!m_bField) // in case of no magnetic field, #track par=4 instead of 5.
163 ndf += 1;
164
165 double Chi2 = fs->getChi2();
166 double TrPval = std::max(0., ROOT::Math::chisquared_cdf_c(Chi2, ndf));
167 double Phi0 = fitresult->getPhi0();
168 if (m_phi0InRad == false) { // unit of degrees.
169 Phi0 *= 180 / M_PI;
170 }
171
172 /*** Two track case.***/
173 if ((posSeed.Y() > 0 && !up) ||
174 (up && posSeed.Y()*posSeed1.Y() > 0 && ndf > ndf1)) {
175 charge1 = fitresult->getChargeSign();
176 posSeed1 = posSeed;
177 D01 = fitresult->getD0();
178 eD01 = sqrt(fitresult->getCovariance5()[0][0]);
179 Phi01 = Phi0;
180 ePhi01 = sqrt(fitresult->getCovariance5()[1][1]);
181 if (m_phi0InRad == false) { // unit of degrees.
182 ePhi01 *= 180 / M_PI;
183 }
184
185 Omega1 = fitresult->getOmega();
186 Mom1 = B2Vector3D(fitresult->getMomentum());
187 eOm1 = sqrt(fitresult->getCovariance5()[2][2]);
188 Z01 = fitresult->getZ0();
189 eZ01 = sqrt(fitresult->getCovariance5()[3][3]);
190 tanLambda1 = fitresult->getTanLambda();
191 etanL1 = sqrt(fitresult->getCovariance5()[4][4]);
192 Pt1 = fitresult->getTransverseMomentum();
193 ndf1 = ndf;
194 Pval1 = TrPval;
195 up = true;
196 }
197
198 if ((posSeed.Y() < 0 && !down) ||
199 (down && posSeed.Y()*posSeed2.Y() > 0 && ndf > ndf2)) {
200 charge2 = fitresult->getChargeSign();
201 posSeed2 = posSeed;
202 D02 = fitresult->getD0();
203 eD02 = sqrt(fitresult->getCovariance5()[0][0]);
204 Phi02 = Phi0;
205 ePhi02 = sqrt(fitresult->getCovariance5()[1][1]);
206 if (m_phi0InRad == false) { // unit of degrees.
207 ePhi02 *= 180 / M_PI;
208 }
209 Omega2 = fitresult->getOmega();
210 Mom2 = B2Vector3D(fitresult->getMomentum());
211 eOm2 = sqrt(fitresult->getCovariance5()[2][2]);
212 Z02 = fitresult->getZ0();
213 eZ02 = sqrt(fitresult->getCovariance5()[3][3]);
214 tanLambda2 = fitresult->getTanLambda();
215 etanL2 = sqrt(fitresult->getCovariance5()[4][4]);
216 Pt2 = fitresult->getTransverseMomentum();
217 ndf2 = ndf;
218 Pval2 = TrPval;
219 down = true;
220 }
221 }
222 if (m_bField && charge1 * charge2 == 0) return;
223 if (charge1 * charge2 >= 0 && up && down) {
224 charge = charge1;
225 tree->Fill();
226 }
227
228}
229
233
235{
236 tfile->cd();
237 tree->Write();
238 tfile->Close();
239}
240
241
DataType Y() const
access variable Y (= .at(1) without boundary check)
Definition B2Vector3.h:433
static ROOT::Math::XYZVector getFieldInTesla(const ROOT::Math::XYZVector &pos)
return the magnetic field at a given position in Tesla.
StoreObjPtr< EventT0 > m_eventTimeStoreObject
Event t0.
std::string m_recoTrackArrayName
Belle2::RecoTrack StoreArray name.e.
void initialize() override
Initializes the Module.
StoreArray< TrackFitResult > m_TrackFitResults
Track fit results.
void event() override
Event action (main routine).
double ndf1
degree of freedom of 1st track.
void terminate() override
Termination action.
bool m_phi0InRad
Unit of phi0, true: radian, false: degree.
double ndf2
degree of freedom of 2nd track.
std::string m_relRecoTrackTrackName
Relation between RecoTrack and Belle2:Track.
TVector3 posSeed1
seed position of the first track.
double etanL1
error on TanLambda of 1st track.
bool m_eventT0Extraction
run with event t0 extraction
void beginRun() override
Begin run action.
double eOm1
error on Omega of 1st track.
bool m_bField
Data are taken with B-field or not, if true, NDF=5 in cal P-value.
double eOm2
error on Omega of 2nd track.
double etanL2
error on TanLambda of 2nd track.
std::string m_trackArrayName
Belle2::Track StoreArray name.
TVector3 posSeed2
seed position of the second track.
TTree * tree
output tree, save info of each hit.
std::string m_trackFitResultArrayName
Belle2::TrackFitResult StoreArray name.
double ePhi02
error on Phi0 of 2nd track.
StoreArray< RecoTrack > m_RecoTracks
Tracks.
double ePhi01
error on Phi0 of 1st track.
StoreObjPtr< EventMetaData > m_EventMetaData
Event metadata.
std::string m_outputFileName
Output file name.
bool m_storeTrackParErrors
Store error of track parameters or not.
static const ChargedStable muon
muon particle
Definition Const.h:660
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition Module.cc:208
Module()
Constructor.
Definition Module.cc:30
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Definition Module.h:80
This is the Reconstruction Event-Data Model Track.
Definition RecoTrack.h:79
ROOT::Math::XYZVector getPositionSeed() const
Return the position seed stored in the reco track. ATTENTION: This is not the fitted position.
Definition RecoTrack.h:480
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:621
Low-level class to create/modify relations between StoreArrays.
TO * getRelatedTo(const std::string &name="", const std::string &namedRelation="") const
Get the object to which this object has a relation.
const std::string & getName() const
Return name under which the object is saved in the DataStore.
Values of the result of a track fit with a given particle hypothesis.
TMatrixDSym getCovariance5() const
Getter for covariance matrix of perigee parameters in matrix form.
short getChargeSign() const
Return track charge (1 or -1).
double getOmega() const
Getter for omega.
double getD0() const
Getter for d0.
double getTransverseMomentum() const
Getter for the absolute value of the transverse momentum at the perigee.
double getTanLambda() const
Getter for tanLambda.
double getZ0() const
Getter for z0.
ROOT::Math::XYZVector getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
double getPhi0() const
Getter for phi0.
Class that bundles various TrackFitResults.
Definition Track.h:25
const TrackFitResult * getTrackFitResult(const Const::ChargedStable &chargedStable) const
Default Access to TrackFitResults.
Definition Track.cc:30
const TrackFitResult * getTrackFitResultWithClosestMass(const Const::ChargedStable &requestedType) const
Return the track fit for the fit hypothesis with the closest mass.
Definition Track.cc:104
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition B2Vector3.h:516
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Abstract base class for different kinds of events.
STL namespace.