Belle II Software development
CDCFudgeFactorCalibrationCollector.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/cdcCalibrationCollector/CDCFudgeFactorCalibrationCollector.h"
10#include "analysis/utility/PCmsLabTransform.h"
11#include <mdst/dataobjects/HitPatternVXD.h>
12#include <mdst/dataobjects/EventLevelTrackingInfo.h>
13#include <analysis/dataobjects/Particle.h>
14#include <analysis/dataobjects/ParticleList.h>
15#include <analysis/variables/PIDVariables.h>
16
17#include <TH1F.h>
18#include <TH2F.h>
19
20using namespace std;
21using namespace Belle2;
22using namespace CDC;
23
24REG_MODULE(CDCFudgeFactorCalibrationCollector);
25
27{
28 setDescription("Collector module for cdc fudege calibration");
29 setPropertyFlags(c_ParallelProcessingCertified); // specify this flag if you need parallel processing
30 addParam("StoreNtuple", m_StoreNtuple, "Store ntuple other studies", true);
31 addParam("MinColinearityTheta", m_minCollinearityTheta, "cut on colinear of the two track by theta", 10.);
32 addParam("MinColinearityPhi0", m_minCollinearityPhi0, "cut on colinear of the two track by phi0", 10.);
33 addParam("DiMuonListName", m_DiMuonListName, "name of the di-muon list", std::string("vpho:mumu"));
34 addParam("GammaListName", m_GammaListName, "name of the gamma list", std::string("gamma:HLT"));
35
36}
37
41
43{
44 m_Tracks.isRequired(m_trackArrayName);
47
48 if (m_StoreNtuple) {
49 auto m_tree = new TTree(m_treeName.c_str(), "tree for cdc calibration");
50 m_tree->Branch<Int_t>("exp_run", &expRun);
51 m_tree->Branch<Float_t>("pt_pos", &ptPos);
52 m_tree->Branch<Float_t>("pt_neg", &ptNeg);
53 m_tree->Branch<Float_t>("pz_pos", &pzPos);
54 m_tree->Branch<Float_t>("pz_neg", &pzNeg);
55
56 m_tree->Branch<Float_t>("pt_pos_cm", &ptPosCm);
57 m_tree->Branch<Float_t>("pt_neg_cm", &ptNegCm);
58 m_tree->Branch<Float_t>("pz_pos_cm", &pzPosCm);
59 m_tree->Branch<Float_t>("pz_neg_cm", &pzNegCm);
60
61 m_tree->Branch<Float_t>("theta_pos_cm", &thetaPosCm);
62 m_tree->Branch<Float_t>("theta_neg_cm", &thetaNegCm);
63
64 m_tree->Branch<Float_t>("phi0_pos_cm", &phi0PosCm);
65 m_tree->Branch<Float_t>("theta_neg_cm", &phi0NegCm);
66
67 m_tree->Branch<Float_t>("Pval_pos", &pvalPos);
68 m_tree->Branch<Float_t>("Pval_neg", &pvalNeg);
69
70 m_tree->Branch<Float_t>("ndf_pos", &ndfPos);
71 m_tree->Branch<Float_t>("ndf_neg", &ndfNeg);
72
73 m_tree->Branch<Float_t>("ncdc_pos", &ncdcPos);
74 m_tree->Branch<Float_t>("ncdc_neg", &ncdcNeg);
75 m_tree->Branch<Float_t>("npxd_pos", &npxdPos);
76 m_tree->Branch<Float_t>("npxd_neg", &npxdNeg);
77 m_tree->Branch<Float_t>("nsvd_pos", &nsvdPos);
78 m_tree->Branch<Float_t>("nsvd_neg", &nsvdNeg);
79
80 m_tree->Branch<Float_t>("nextra_cdchit", &nExtraCDCHits);
81 m_tree->Branch<Float_t>("ecl_track", &eclTrack);
82 m_tree->Branch<Float_t>("ecl_neutral", &eclNeutral);
83
84 m_tree->Branch<Float_t>("d0_pos", &d0Pos);
85 m_tree->Branch<Float_t>("d0_neg", &d0Neg);
86 m_tree->Branch<Float_t>("z0_pos", &z0Pos);
87 m_tree->Branch<Float_t>("z0_neg", &z0Neg);
88 m_tree->Branch<Float_t>("d0ip_pos", &d0ipPos);
89 m_tree->Branch<Float_t>("d0ip_neg", &d0ipNeg);
90 m_tree->Branch<Float_t>("z0ip_pos", &z0ipPos);
91 m_tree->Branch<Float_t>("z0ip_neg", &z0ipNeg);
92 m_tree->Branch<Float_t>("muid_pos", &muidPos);
93 m_tree->Branch<Float_t>("muid_neg", &muidNeg);
94 m_tree->Branch<Float_t>("eid_pos", &eidPos);
95 m_tree->Branch<Float_t>("eid_neg", &eidNeg);
96
97 registerObject<TTree>(m_treeName.c_str(), m_tree);
98 }
99 auto m_hEventT0 = new TH1F("hEventT0", "Event T0", 200, -100, 100);
100 auto m_hExtraCDCHit = new TH1F("hExtraCDCHit", "Extra cdc hits", 500, 0, 5000);
101 auto m_hNDF_pos = new TH1F("hNDF_pos", "NDF of positive track;NDF;Tracks", 71, -1, 70);
102 auto m_hNDF_neg = new TH1F("hNDF_neg", "NDF of negative track;NDF;Tracks", 71, -1, 70);
103 auto m_hPval_pos = new TH1F("hPval_pos", "p-values of pos tracks;pVal;Tracks", 1000, 0, 1);
104 auto m_hPval_neg = new TH1F("hPval_neg", "p-values of neg tra cks;pVal;Tracks", 1000, 0, 1);
105 auto m_hnCDC_pos = new TH1F("hnCDC_pos", "nCDC hit of positive track ; nCDC ; Tracks", 71, -1, 70);
106 auto m_hnCDC_neg = new TH1F("hnCDC_neg", "nCDC hit of negative track ; nCDC ; Tracks", 71, -1, 70);
107
108 auto m_hdPt = new TH1F("hdPt", "#DeltaP_{t} ; #DeltaP_{t} ; Events ", 200, -0.5, 0.5);
109 auto m_hdD0 = new TH1F("hdD0", "#DeltaD_{0} ; #DeltaD_{0} ; Events ", 200, -0.1, 0.1);
110 auto m_hdZ0 = new TH1F("hdZ0", "#DeltaZ_{0} ; #DeltaD_{0} ; Events ", 200, -1.5, 1.5);
111
112 auto m_hdPt_cm = new TH1F("hdPt_cm", "#DeltaP_{t} in c.m frame ; #DeltaP_{t} (c.m) ; Events ", 200, -0.5, 0.5);
113 auto m_hdPtPt_cm = new TH2F("hdPtPt_cm", "#DeltaP_{t}:P_{t} in c.m frame ;P_{t} GeV/c ; #DeltaP_{t} (c.m) ", 50, 2., 7., 200, -0.5,
114 0.5);
115
116 auto m_hdPhi0_cm = new TH1F("hdPhi0_cm", "#DeltaPhi_{0} in c.m frame ; #Delta#Phi_{0} in c.m ; Events ", 200, -0.7, 0.7);
117 auto m_hdTheta_cm = new TH1F("hdTheta_cm", "#Delta#theta in c.m frame ; #Delta#theta in c.m ; Events ", 200, -3.0, 3.0);
118
119 registerObject<TH1F>("hEventT0", m_hEventT0);
120 registerObject<TH1F>("hExtraCDCHit", m_hExtraCDCHit);
121 registerObject<TH1F>("hNDF_pos", m_hNDF_pos);
122 registerObject<TH1F>("hNDF_neg", m_hNDF_neg);
123 registerObject<TH1F>("hnCDC_pos", m_hnCDC_pos);
124 registerObject<TH1F>("hnCDC_neg", m_hnCDC_neg);
125
126 registerObject<TH1F>("hPval_pos", m_hPval_pos);
127 registerObject<TH1F>("hPval_neg", m_hPval_neg);
128
129
130 registerObject<TH1F>("hdPt", m_hdPt);
131 registerObject<TH1F>("hdD0", m_hdD0);
132 registerObject<TH1F>("hdZ0", m_hdZ0);
133
134 registerObject<TH1F>("hdPt_cm", m_hdPt_cm);
135 registerObject<TH2F>("hdPtPt_cm", m_hdPtPt_cm);
136 registerObject<TH1F>("hdPhi0_cm", m_hdPhi0_cm);
137 registerObject<TH1F>("hdTheta_cm", m_hdTheta_cm);
138}
139
141{
142 int nCandidates = m_DiMuonList->getListSize();
143 B2DEBUG(29, "Number of muon canndiate:" << nCandidates);
144 if (nCandidates < 1) return;
145
146 // event with is fail to extract t0 will be exclude from analysis
147 if (m_eventTimeStoreObject.isValid() && m_eventTimeStoreObject->hasEventT0()) {
148 getObjectPtr<TH1F>("hEventT0")->Fill(m_eventTimeStoreObject->getEventT0());
149 } else {
150 return;
151 }
152 //store run/exp info
153 StoreObjPtr<EventMetaData> m_EventMetaData;
154 int run = m_EventMetaData->getRun();
155 int exp = m_EventMetaData->getExperiment();
156 expRun = exp * 1000000 + run;
157
158 /************************************/
161 int nG = gamma_list->getListSize();
162 B2DEBUG(29, "Number of gamma: " << nG);
163 eclNeutral = 0;
164 eclTrack = 0;
165 for (int i = 0; i < nG; ++i) {
166 Particle* gamma = gamma_list->getParticle(i);
167 eclNeutral += gamma->getEnergy();
168 }
169 B2DEBUG(29, "Sum of neutral ECL " << eclNeutral);
170 /************************************/
172 //now start to collect dimuon parameters
173 double thetaPos(0), thetaNeg(0);
174 int charge_sum = 0;
175
176 for (int i = 0; i < nCandidates; ++i) {
177 Particle* part = m_DiMuonList->getParticle(i);
178 //vertex from vertex fit
179 ROOT::Math::XYZVector v0Vertex = part->getVertex();
180
181 //loop over its daugter
182 for (int j = 0; j < 2; ++j) {
183 //we have two daugters.
184 const Particle* d0 = part->getDaughter(j);
185 short chg = d0->getCharge();
186
187 const Belle2::TrackFitResult* fitresult = d0->getTrackFitResult();
188 if (!fitresult) {
189 B2WARNING("No track fit result found.");
190 break;
191 }
192 //get Cluter Energy
193 eclTrack += d0->getECLClusterEnergy();
194 double muid = Variable::muonID(d0);
195 double eid = Variable::electronID(d0);
196
197 if (chg > 0) {
198 ndfPos = fitresult->getNDF();
199 pvalPos = fitresult->getPValue();
200 ptPos = fitresult->getTransverseMomentum();
201 pzPos = fitresult->getMomentum().Z();
202 d0Pos = fitresult->getD0();
203 z0Pos = fitresult->getZ0();
204 thetaPos = fitresult->getMomentum().Theta() * 180 / M_PI;
205 ncdcPos = fitresult->getHitPatternCDC().getNHits();
206 nsvdPos = fitresult->getHitPatternVXD().getNSVDHits();
207 npxdPos = fitresult->getHitPatternVXD().getNPXDHits();
208 ROOT::Math::PxPyPzEVector P4_pos = T.rotateLabToCms() * fitresult->get4Momentum();
209 ptPosCm = P4_pos.Pt();
210 thetaPosCm = P4_pos.Theta() * 180 / M_PI;
211 phi0PosCm = P4_pos.Phi() * 180 / M_PI;
212 UncertainHelix helix_pos = fitresult->getUncertainHelix();
213 helix_pos.passiveMoveBy(v0Vertex);
214 d0ipPos = helix_pos.getD0();
215 z0ipPos = helix_pos.getZ0();
216 muidPos = muid; eidPos = eid;
217 } else if (chg < 0) {
218 ndfNeg = fitresult->getNDF();
219 pvalNeg = fitresult->getPValue();
220 ptNeg = fitresult->getTransverseMomentum();
221 pzNeg = fitresult->getMomentum().Z();
222 d0Neg = fitresult->getD0();
223 z0Neg = fitresult->getZ0();
224 thetaNeg = fitresult->getMomentum().Theta() * 180 / M_PI;
225 ncdcNeg = fitresult->getHitPatternCDC().getNHits();
226 nsvdNeg = fitresult->getHitPatternVXD().getNSVDHits();
227 npxdNeg = fitresult->getHitPatternVXD().getNPXDHits();
228
229 ROOT::Math::PxPyPzEVector P4_neg = T.rotateLabToCms() * fitresult->get4Momentum();
230 ptNegCm = P4_neg.Pt();
231 thetaNegCm = P4_neg.Theta() * 180 / M_PI;
232 phi0NegCm = P4_neg.Phi() * 180 / M_PI;
233 UncertainHelix helix_neg = fitresult->getUncertainHelix();
234 helix_neg.passiveMoveBy(v0Vertex);
235 d0ipNeg = helix_neg.getD0();
236 z0ipNeg = helix_neg.getZ0();
237 muidNeg = muid;
238 eidNeg = eid;
239 } else {
240 continue;
241 }
242 //count #good charged track and sum of charged
243 charge_sum += chg;
244 }
245 //save extra cdc hit for background monitoring
247 if (!elti) nExtraCDCHits = -1;
248 else nExtraCDCHits = elti->getNCDCHitsNotAssigned();
249 getObjectPtr<TH1F>("hExtraCDCHit")->Fill(nExtraCDCHits);
250 // cut good charged track
251 if (charge_sum != 0) {continue;}
252 //cut according to the Line400 in the plotMumu_4ff.C
253 if (ptPos < 2.5 || ptNeg < 2.5
254 || thetaPos < 45 || thetaPos > 125
255 || thetaNeg < 45 || thetaNeg > 125) {return;}
256
257 // cut on Energy deposite in ECL
258 // Keep the same as in the dimuon study script,
259 // Although it is not necessary to keep this as cut on Etot is enough
260 double eclTot = eclNeutral + eclTrack;
261 if (eclTot > 2 || eclTrack > 2) return;
262
263 // B2INFO("Total Eecl_track = "<< Eecl_trk);
264 double dPt = (ptPos - ptNeg) / sqrt(2);
265 double dD0 = (d0Pos + d0Neg) / sqrt(2);
266 double dZ0 = (z0Pos - z0Neg) / sqrt(2);
267 double Pt_cm = (ptPosCm + ptNegCm) / 2;
268 double dPt_cm = (ptPosCm - ptNegCm) / sqrt(2);
269 double dPhi0_cm = (180 - fabs(phi0PosCm - phi0NegCm)) / sqrt(2);
270 double dTheta_cm = (180 - fabs(thetaPosCm + thetaNegCm)) / sqrt(2);
271
272 //require back to back
273 if (fabs(dPhi0_cm) > m_minCollinearityPhi0 || fabs(dTheta_cm) > m_minCollinearityTheta) return;
274
275 getObjectPtr<TH1F>("hPval_pos")->Fill(pvalPos);
276 getObjectPtr<TH1F>("hPval_neg")->Fill(pvalNeg);
277 getObjectPtr<TH1F>("hNDF_pos")->Fill(ndfPos);
278 getObjectPtr<TH1F>("hNDF_neg")->Fill(ndfNeg);
279 getObjectPtr<TH1F>("hnCDC_pos")->Fill(ncdcPos);
280 getObjectPtr<TH1F>("hnCDC_neg")->Fill(ncdcNeg);
281
282 getObjectPtr<TH1F>("hdPt")->Fill(dPt);
283 getObjectPtr<TH1F>("hdD0")->Fill(dD0);
284 getObjectPtr<TH1F>("hdZ0")->Fill(dZ0);
285
286 getObjectPtr<TH1F>("hdPt_cm")->Fill(dPt_cm);
287 getObjectPtr<TH2F>("hdPtPt_cm")->Fill(Pt_cm, dPt_cm);
288 getObjectPtr<TH1F>("hdPhi0_cm")->Fill(dPhi0_cm);
289 getObjectPtr<TH1F>("hdTheta_cm")->Fill(dTheta_cm);
290
291 if (m_StoreNtuple) {
292 getObjectPtr<TTree>(m_treeName.c_str())->Fill();
293 }
294 }
295}
296
Double_t m_minCollinearityTheta
Minimum requirement for accolinear theta in c.m frame.
Double_t m_minCollinearityPhi0
Minimum requirement for accolinear phi0 in c.m frame.
Float_t pzPos
Longitudinal momentum of the positive track.
Int_t expRun
Exp and run numbers, encoded by exp*10^6+run.
StoreArray< TrackFitResult > m_TrackFitResults
Track fit results.
Float_t nExtraCDCHits
Number of CDC hits not assigned to any tracks.
Float_t ptNeg
Transeverse momentum of the negative track.
Float_t thetaNegCm
theta of the negative track in c.m frame.
std::string m_DiMuonListName
List name for the reconstruted dimuon.
void collect() override
Event action, collect information for calibration.
Float_t ptPos
Transeverse momentum of the positive track.
bool m_StoreNtuple
Option to store ntuple, =true: tree with these variables will be stored.
StoreObjPtr< ParticleList > m_DiMuonList
List of the reconstructed dimion.
Float_t thetaPosCm
theta of the positive track in c.m frame.
Float_t pzPosCm
Longitudinal momentum of the positive track in c.m frame.
std::string m_trackFitResultArrayName
Belle2::TrackFitResult StoreArray name.
Float_t ptPosCm
Transeverse momentum of the positive track in c.m frame.
Float_t ptNegCm
Transeverse momentum of the negative track in c.m frame.
Float_t pzNegCm
Longitudinal momentum of the negative track in c.m frame.
std::string m_GammaListName
List name for the reconstruted dimuon.
Float_t pzNeg
Longitudinal momentum of the negative track.
void registerObject(std::string name, T *obj)
Register object with a name, takes ownership, do not access the pointer beyond prepare()
CalibrationCollectorModule()
Constructor. Sets the default prefix for calibration dataobjects.
T * getObjectPtr(std::string name)
Calls the CalibObjManager to get the requested stored collector data.
unsigned short getNHits() const
Get the total Number of CDC hits in the fit.
unsigned short getNPXDHits() const
Get total number of hits in the PXD.
unsigned short getNSVDHits() const
Get total number of hits in the SVD.
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
@ 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
Class to hold Lorentz transformations from/to CMS and boost vector.
Class to store reconstructed particles.
Definition Particle.h:76
double getEnergy() const
Returns total energy.
Definition Particle.h:555
ROOT::Math::XYZVector getVertex() const
Returns vertex position (POCA for charged, IP for neutral FS particles)
Definition Particle.h:651
const Particle * getDaughter(unsigned i) const
Returns a pointer to the i-th daughter particle.
Definition Particle.cc:662
Type-safe access to single objects in the data store.
Definition StoreObjPtr.h:96
Values of the result of a track fit with a given particle hypothesis.
float getNDF() const
Getter for number of degrees of freedom of the track fit.
double getPValue() const
Getter for Chi2 Probability of the track fit.
ROOT::Math::PxPyPzEVector get4Momentum() const
Getter for the 4Momentum at the closest approach of the track in the r/phi projection.
double getD0() const
Getter for d0.
double getTransverseMomentum() const
Getter for the absolute value of the transverse momentum at the perigee.
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.
HitPatternCDC getHitPatternCDC() const
Getter for the hit pattern in the CDC;.
UncertainHelix getUncertainHelix() const
Conversion to framework Uncertain Helix (i.e., with covariance).
HitPatternVXD getHitPatternVXD() const
Getter for the hit pattern in the VXD;.
This class represents an ideal helix in perigee parameterization including the covariance matrix of t...
double passiveMoveBy(const ROOT::Math::XYZVector &by)
Moves origin of the coordinate system (passive transformation) by the given vector.
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
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Abstract base class for different kinds of events.
STL namespace.