Belle II Software  release-08-01-10
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 <framework/dataobjects/Helix.h>
12 #include <mdst/dataobjects/ECLCluster.h>
13 #include <analysis/dataobjects/Particle.h>
14 #include <analysis/dataobjects/ParticleList.h>
15 #include <Math/ProbFuncMathCore.h>
16 
17 #include <TH1F.h>
18 #include <TH2F.h>
19 #include <TMath.h>
20 
21 using namespace std;
22 using namespace Belle2;
23 using namespace CDC;
24 
25 REG_MODULE(CDCFudgeFactorCalibrationCollector);
26 
27 CDCFudgeFactorCalibrationCollectorModule::CDCFudgeFactorCalibrationCollectorModule() : CalibrationCollectorModule()
28 {
29  setDescription("Collector module for cdc fudege calibration");
30  setPropertyFlags(c_ParallelProcessingCertified); // specify this flag if you need parallel processing
31  addParam("StoreNtuple", m_StoreNtuple, "Store ntuple other studies", true);
32  addParam("MinColinearityTheta", m_minCollinearityTheta, "cut on colinear of the two track by theta", 10.);
33  addParam("MinColinearityPhi0", m_minCollinearityPhi0, "cut on colinear of the two track by phi0", 10.);
34  addParam("DiMuonListName", m_DiMuonListName, "name of the di-muon list", std::string("vpho:mumu"));
35  addParam("GammaListName", m_GammaListName, "name of the gamma list", std::string("gamma:HLT"));
36 
37 }
38 
40 {
41 }
42 
44 {
45  m_Tracks.isRequired(m_trackArrayName);
47  m_DiMuonList.isRequired(m_DiMuonListName);
48 
49  if (m_StoreNtuple) {
50  auto m_tree = new TTree(m_treeName.c_str(), "tree for cdc calibration");
51  m_tree->Branch<Float_t>("pt_pos", &Pt_pos);
52  m_tree->Branch<Int_t>("exp_run", &exp_run);
53  m_tree->Branch<Float_t>("pt_neg", &Pt_neg);
54  m_tree->Branch<Float_t>("pt_pos_cm", &Pt_pos_cm);
55  m_tree->Branch<Float_t>("pt_neg_cm", &Pt_neg_cm);
56  m_tree->Branch<Float_t>("theta_pos_cm", &Theta_pos_cm);
57  m_tree->Branch<Float_t>("theta_neg_cm", &Theta_neg_cm);
58  m_tree->Branch<Float_t>("phi0_pos_cm", &Phi0_pos_cm);
59  m_tree->Branch<Float_t>("theta_neg_cm", &Phi0_neg_cm);
60 
61  m_tree->Branch<Float_t>("Pval_pos", &Pval_pos);
62  m_tree->Branch<Float_t>("Pval_neg", &Pval_neg);
63  m_tree->Branch<Float_t>("ndf_pos", &ndf_pos);
64  m_tree->Branch<Float_t>("ndf_neg", &ndf_neg);
65  m_tree->Branch<Float_t>("ncdc_pos", &nCDC_pos);
66  m_tree->Branch<Float_t>("ncdc_neg", &nCDC_neg);
67  m_tree->Branch<Float_t>("d0_pos", &D0_pos);
68  m_tree->Branch<Float_t>("d0_neg", &D0_neg);
69  m_tree->Branch<Float_t>("z0_pos", &Z0_pos);
70  m_tree->Branch<Float_t>("z0_neg", &Z0_neg);
71  m_tree->Branch<Float_t>("d0ip_pos", &D0ip_pos);
72  m_tree->Branch<Float_t>("d0ip_neg", &D0ip_neg);
73  m_tree->Branch<Float_t>("z0ip_pos", &Z0ip_pos);
74  m_tree->Branch<Float_t>("z0ip_neg", &Z0ip_neg);
75 
76 
77  registerObject<TTree>(m_treeName.c_str(), m_tree);
78  }
79  auto m_hEventT0 = new TH1F("hEventT0", "Event T0", 1000, -100, 100);
80  auto m_hNDF_pos = new TH1F("hNDF_pos", "NDF of positive track;NDF;Tracks", 71, -1, 70);
81  auto m_hNDF_neg = new TH1F("hNDF_neg", "NDF of negative track;NDF;Tracks", 71, -1, 70);
82  auto m_hPval_pos = new TH1F("hPval_pos", "p-values of pos tracks;pVal;Tracks", 1000, 0, 1);
83  auto m_hPval_neg = new TH1F("hPval_neg", "p-values of neg tra cks;pVal;Tracks", 1000, 0, 1);
84  auto m_hnCDC_pos = new TH1F("hnCDC_pos", "nCDC hit of positive track ; nCDC ; Tracks", 71, -1, 70);
85  auto m_hnCDC_neg = new TH1F("hnCDC_neg", "nCDC hit of negative track ; nCDC ; Tracks", 71, -1, 70);
86 
87 
88 
89  auto m_hdPt = new TH1F("hdPt", "#DeltaP_{t} ; #DeltaP_{t} ; Events ", 200, -0.5, 0.5);
90  auto m_hdD0 = new TH1F("hdD0", "#DeltaD_{0} ; #DeltaD_{0} ; Events ", 200, -0.1, 0.1);
91  auto m_hdZ0 = new TH1F("hdZ0", "#DeltaZ_{0} ; #DeltaD_{0} ; Events ", 200, -1.5, 1.5);
92 
93  auto m_hdPt_cm = new TH1F("hdPt_cm", "#DeltaP_{t} in c.m frame ; #DeltaP_{t} (c.m) ; Events ", 200, -0.5, 0.5);
94  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,
95  0.5);
96 
97  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);
98  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);
99 
100  registerObject<TH1F>("hEventT0", m_hEventT0);
101  registerObject<TH1F>("hNDF_pos", m_hNDF_pos);
102  registerObject<TH1F>("hNDF_neg", m_hNDF_neg);
103  registerObject<TH1F>("hnCDC_pos", m_hnCDC_pos);
104  registerObject<TH1F>("hnCDC_neg", m_hnCDC_neg);
105 
106  registerObject<TH1F>("hPval_pos", m_hPval_pos);
107  registerObject<TH1F>("hPval_neg", m_hPval_neg);
108 
109 
110  registerObject<TH1F>("hdPt", m_hdPt);
111  registerObject<TH1F>("hdD0", m_hdD0);
112  registerObject<TH1F>("hdZ0", m_hdZ0);
113 
114  registerObject<TH1F>("hdPt_cm", m_hdPt_cm);
115  registerObject<TH2F>("hdPtPt_cm", m_hdPtPt_cm);
116  registerObject<TH1F>("hdPhi0_cm", m_hdPhi0_cm);
117  registerObject<TH1F>("hdTheta_cm", m_hdTheta_cm);
118 }
119 
121 {
122  int nCandidates = m_DiMuonList->getListSize();
123  B2DEBUG(199, "Number of muon canndiate:" << nCandidates);
124  if (nCandidates < 1) return;
125 
126  // event with is fail to extract t0 will be exclude from analysis
127  if (m_eventTimeStoreObject.isValid() && m_eventTimeStoreObject->hasEventT0()) {
128  getObjectPtr<TH1F>("hEventT0")->Fill(m_eventTimeStoreObject->getEventT0());
129  } else {
130  return;
131  }
132  //store run/exp info
133  StoreObjPtr<EventMetaData> m_EventMetaData;
134  int run = m_EventMetaData->getRun();
135  int exp = m_EventMetaData->getExperiment();
136  exp_run = exp * 1000000 + run;
137 
138  /************************************/
141  double Eecl_neutral = 0;
142  int nG = gamma_list->getListSize();
143  B2DEBUG(199, "Number of gamma: " << nG);
144  for (int i = 0; i < nG; ++i) {
145  Particle* gamma = gamma_list->getParticle(i);
146  Eecl_neutral += gamma->getEnergy();
147  }
148  B2DEBUG(199, "Sum of neutral ECL " << Eecl_neutral);
149  /************************************/
150 
151  // const int nTr = m_Tracks.getEntries();
153  //now start to collect dimuon parameters
154  double theta_pos(0), theta_neg(0);
155  int charge_sum = 0;
156  double Eecl_trk = 0;
157 
158  for (int i = 0; i < nCandidates; ++i) {
159  Particle* part = m_DiMuonList->getParticle(i);
160  //vertex from vertex fit
161  ROOT::Math::XYZVector v0Vertex = part->getVertex();
162 
163  //loop over its daugter
164  for (int j = 0; j < 2; ++j) {
165  //we have two daugters.
166  const Particle* d0 = part->getDaughter(j);
167  short chg = d0->getCharge();
168 
169  const Belle2::TrackFitResult* fitresult = d0->getTrackFitResult();
170  if (!fitresult) {
171  B2WARNING("No track fit result found.");
172  break;
173  }
174  //get Cluter Energy
175  Eecl_trk += d0->getECLClusterEnergy();
176  if (chg > 0) {
177  ndf_pos = fitresult->getNDF();
178  Pval_pos = fitresult->getPValue();
179  Pt_pos = fitresult->getTransverseMomentum();
180  D0_pos = fitresult->getD0();
181  Z0_pos = fitresult->getZ0();
182  theta_pos = fitresult->getMomentum().Theta() * 180 / M_PI;
183  nCDC_pos = fitresult->getHitPatternCDC().getNHits();
184  ROOT::Math::PxPyPzEVector P4_pos = T.rotateLabToCms() * fitresult->get4Momentum();
185  Pt_pos_cm = P4_pos.Pt();
186  Theta_pos_cm = P4_pos.Theta() * 180 / M_PI;
187  Phi0_pos_cm = P4_pos.Phi() * 180 / M_PI;
188  UncertainHelix helix_pos = fitresult->getUncertainHelix();
189  helix_pos.passiveMoveBy(v0Vertex);
190  D0ip_pos = helix_pos.getD0();
191  Z0ip_pos = helix_pos.getZ0();
192  } else if (chg < 0) {
193  ndf_neg = fitresult->getNDF();
194  Pval_neg = fitresult->getPValue();
195  Pt_neg = fitresult->getTransverseMomentum();
196  D0_neg = fitresult->getD0();
197  Z0_neg = fitresult->getZ0();
198  theta_neg = fitresult->getMomentum().Theta() * 180 / M_PI;
199  nCDC_neg = fitresult->getHitPatternCDC().getNHits();
200  ROOT::Math::PxPyPzEVector P4_neg = T.rotateLabToCms() * fitresult->get4Momentum();
201  Pt_neg_cm = P4_neg.Pt();
202  Theta_neg_cm = P4_neg.Theta() * 180 / M_PI;
203  Phi0_neg_cm = P4_neg.Phi() * 180 / M_PI;
204  UncertainHelix helix_neg = fitresult->getUncertainHelix();
205  helix_neg.passiveMoveBy(v0Vertex);
206  D0ip_neg = helix_neg.getD0();
207  Z0ip_neg = helix_neg.getZ0();
208  } else {
209  continue;
210  }
211  //count #good charged track and sum of charged
212  charge_sum += chg;
213  }
214  // cut good charged track
215  if (charge_sum != 0) {continue;}
216  //cut according to the Line400 in the plotMumu_4ff.C
217  if (Pt_pos < 2.5 || Pt_neg < 2.5
218  || theta_pos < 45 || theta_pos > 125
219  || theta_neg < 45 || theta_neg > 125) {return;}
220 
221  // cut on Energy deposite in ECL
222  double Eecl_tot = Eecl_neutral + Eecl_trk;
223  if (Eecl_tot > 2 || Eecl_trk > 2) return;
224 
225  // B2INFO("Total Eecl_track = "<< Eecl_trk);
226  double dPt = (Pt_pos - Pt_neg) / sqrt(2);
227  double dD0 = (D0_pos + D0_neg) / sqrt(2);
228  double dZ0 = (Z0_pos - Z0_neg) / sqrt(2);
229  double Pt_cm = (Pt_pos_cm + Pt_neg_cm) / 2;
230  double dPt_cm = (Pt_pos_cm - Pt_neg_cm) / sqrt(2);
231  double dPhi0_cm = (180 - fabs(Phi0_pos_cm - Phi0_neg_cm)) / sqrt(2);
232  double dTheta_cm = (180 - fabs(Theta_pos_cm + Theta_neg_cm)) / sqrt(2);
233 
234  B2DEBUG(199, "ECL neutral - trk:" << Eecl_neutral << " - " << Eecl_trk);
235  B2DEBUG(199, "Pos: theta phi0 :" << Theta_pos_cm << " - " << Phi0_pos_cm);
236  B2DEBUG(199, "Neg: theta phi0 :" << Theta_neg_cm << " - " << Phi0_neg_cm);
237  B2DEBUG(199, "DeltaPhi Theta :" << dPhi0_cm << " " << dTheta_cm);
238 
239  //require back to back
240  if (fabs(dPhi0_cm) > m_minCollinearityPhi0 || fabs(dTheta_cm) > m_minCollinearityTheta) return;
241 
242  getObjectPtr<TH1F>("hPval_pos")->Fill(Pval_pos);
243  getObjectPtr<TH1F>("hPval_neg")->Fill(Pval_neg);
244  getObjectPtr<TH1F>("hNDF_pos")->Fill(ndf_pos);
245  getObjectPtr<TH1F>("hNDF_neg")->Fill(ndf_neg);
246  getObjectPtr<TH1F>("hnCDC_pos")->Fill(nCDC_pos);
247  getObjectPtr<TH1F>("hnCDC_neg")->Fill(nCDC_neg);
248 
249  getObjectPtr<TH1F>("hdPt")->Fill(dPt);
250  getObjectPtr<TH1F>("hdD0")->Fill(dD0);
251  getObjectPtr<TH1F>("hdZ0")->Fill(dZ0);
252 
253  getObjectPtr<TH1F>("hdPt_cm")->Fill(dPt_cm);
254  getObjectPtr<TH2F>("hdPtPt_cm")->Fill(Pt_cm, dPt_cm);
255  getObjectPtr<TH1F>("hdPhi0_cm")->Fill(dPhi0_cm);
256  getObjectPtr<TH1F>("hdTheta_cm")->Fill(dTheta_cm);
257 
258  if (m_StoreNtuple) {
259  getObjectPtr<TTree>(m_treeName.c_str())->Fill();
260  }
261  }
262 }
263 
265 {
266 }
Float_t Phi0_neg_cm
phi0 of the negative track in c.m frame.
Double_t m_minCollinearityTheta
Minimum requirement for accolinear theta in c.m frame.
Float_t nCDC_pos
Number of CDC hit of the positive track.
Double_t m_minCollinearityPhi0
Minimum requirement for accolinear phi0 in c.m frame.
Int_t exp_run
Exp and run numbers, encoded by exp*10^6+run.
Float_t Pt_pos_cm
Transeverse momentum of the positive track in c.m frame.
StoreArray< TrackFitResult > m_TrackFitResults
Track fit results.
Float_t Phi0_pos_cm
phi0 of the positive 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 Theta_neg_cm
theta of the negative track in c.m frame.
bool m_StoreNtuple
Option to store ntuple, =true: tree with these variables will be stored.
Float_t Pt_neg_cm
Transeverse momentum of the negative track in c.m frame.
StoreObjPtr< ParticleList > m_DiMuonList
List of the reconstructed dimion.
Float_t Pt_pos
Transeverse momentum of the positive track
Float_t nCDC_neg
Number of CDC hit of the negative track.
std::string m_trackFitResultArrayName
Belle2::TrackFitResult StoreArray name.
Float_t Theta_pos_cm
theta of the positive track in c.m frame.
Float_t Pt_neg
Transeverse momentum of the negative track
std::string m_GammaListName
List name for the reconstruted dimuon.
Calibration collector module base class.
unsigned short getNHits() const
Get the total Number of CDC hits in the fit.
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.
const ROOT::Math::LorentzRotation rotateLabToCms() const
Returns Lorentz transformation from Lab to CMS.
Class to store reconstructed particles.
Definition: Particle.h:75
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).
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:560
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Abstract base class for different kinds of events.