Belle II Software  release-05-02-19
BKLMAnaModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Yinghui GUAN *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 /* Own header. */
12 #include <klm/bklm/modules/bklmAna/BKLMAnaModule.h>
13 
14 /* Belle 2 headers. */
15 #include <framework/dataobjects/EventMetaData.h>
16 #include <framework/datastore/StoreObjPtr.h>
17 #include <mdst/dataobjects/Track.h>
18 #include <mdst/dataobjects/TrackFitResult.h>
19 #include <klm/dataobjects/KLMMuidLikelihood.h>
20 
21 /* CLHEP headers. */
22 #include <CLHEP/Units/SystemOfUnits.h>
23 
24 using namespace std;
25 using namespace Belle2;
26 using namespace CLHEP;
27 
28 //-----------------------------------------------------------------
29 // Register the Module
30 //-----------------------------------------------------------------
31 REG_MODULE(BKLMAna)
32 
33 //-----------------------------------------------------------------
34 // Implementation
35 //-----------------------------------------------------------------
36 
38  m_file(nullptr),
39  m_extTree(nullptr),
40  m_run(0),
41  m_nExtHit(0),
42  m_hdistance(nullptr),
43  m_totalMom(nullptr),
44  m_passMom(nullptr),
45  m_effiMom(nullptr),
46  m_totalYX(nullptr),
47  m_passYX(nullptr),
48  m_effiYX(nullptr),
49  m_totalYZ(nullptr),
50  m_passYZ(nullptr),
51  m_effiYZ(nullptr)
52 {
53  for (int i = 0; i < 200; ++i) {
54  m_extx[i] = 0.0;
55  m_exty[i] = 0.0;
56  m_extz[i] = 0.0;
57  }
58  for (int i = 0; i < 15; ++i) {
59  m_totalThephi[i] = nullptr;
60  m_passThephi[i] = nullptr;
61  m_effiThephi[i] = nullptr;
62  m_totalTrkThephi[i] = nullptr;
63  m_passTrkThephi[i] = nullptr;
64  m_effiTrkThephi[i] = nullptr;
65  }
66  setDescription("analyze bklm efficiency associated to CDC, check performance of bklm et al.");
67  addParam("filename", m_filename, "Output root filename", string("bklmana.root"));
68 }
69 
70 BKLMAnaModule::~BKLMAnaModule()
71 {
72 }
73 
74 void BKLMAnaModule::initialize()
75 {
76  hits2D.isRequired();
77  extHits.isRequired();
78  tracks.isRequired();
79 
80 
81  m_file = new TFile(m_filename.c_str(), "RECREATE");
82 
83  m_extTree = new TTree("exthit", "ext hit");
84  m_extTree->Branch("run", &m_run, "m_run/I");
85  m_extTree->Branch("nExtHit", &m_nExtHit, "m_nExtHit/I");
86  m_extTree->Branch("x", &m_extx, "m_extx[m_nExtHit]/F");
87  m_extTree->Branch("y", &m_exty, "m_exty[m_nExtHit]/F");
88  m_extTree->Branch("z", &m_extz, "m_extz[m_nExtHit]/F");
89 
90  float phiBins = 36;
91  float phiMin = 0.0;
92  float phiMax = 360.0;
93  float thetaBins = 19;
94  float thetaMin = 35.0;
95  float thetaMax = 130.0;
96  char hname[50];
97  for (int iL = 0; iL < 15 ; iL ++) {
98  //based on theta phi of ExtHit position
99  sprintf(hname, "denominator_Layer%i", iL + 1);
100  m_totalThephi[iL] = new TH2F(hname, hname, phiBins, phiMin, phiMax, thetaBins, thetaMin, thetaMax);
101  m_totalThephi[iL]->GetXaxis()->SetTitle("phi");
102  m_totalThephi[iL]->GetYaxis()->SetTitle("theta");
103  sprintf(hname, "numerator_Layer%i", iL + 1);
104  m_passThephi[iL] = new TH2F(hname, hname, phiBins, phiMin, phiMax, thetaBins, thetaMin, thetaMax);
105  sprintf(hname, "effi_Layer%i", iL + 1);
106  m_passThephi[iL]->GetXaxis()->SetTitle("phi");
107  m_passThephi[iL]->GetYaxis()->SetTitle("theta");
108  m_effiThephi[iL] = new TH2F(hname, hname, phiBins, phiMin, phiMax, thetaBins, thetaMin, thetaMax);
109  m_effiThephi[iL]->GetXaxis()->SetTitle("phi");
110  m_effiThephi[iL]->GetYaxis()->SetTitle("theta");
111  //based on theta phi of trk
112  sprintf(hname, "Denominator_Layer%i", iL + 1);
113  m_totalTrkThephi[iL] = new TH2F(hname, hname, phiBins, phiMin, phiMax, thetaBins, thetaMin, thetaMax);
114  m_totalTrkThephi[iL]->GetXaxis()->SetTitle("trk.phi");
115  m_totalTrkThephi[iL]->GetYaxis()->SetTitle("trk.theta");
116  sprintf(hname, "Numerator_Layer%i", iL + 1);
117  m_passTrkThephi[iL] = new TH2F(hname, hname, phiBins, phiMin, phiMax, thetaBins, thetaMin, thetaMax);
118  sprintf(hname, "Effi_Layer%i", iL + 1);
119  m_passTrkThephi[iL]->GetXaxis()->SetTitle("trk.phi");
120  m_passTrkThephi[iL]->GetYaxis()->SetTitle("trk.theta");
121  m_effiTrkThephi[iL] = new TH2F(hname, hname, phiBins, phiMin, phiMax, thetaBins, thetaMin, thetaMax);
122  m_effiTrkThephi[iL]->GetXaxis()->SetTitle("trk.phi");
123  m_effiTrkThephi[iL]->GetYaxis()->SetTitle("trk.theta");
124  }
125  float gmin = -350;
126  float gmax = 350;
127  int gNbin = 150;
128  float mommin = 0.;
129  float mommax = 15;
130  int mNbin = 30;
131 
132  m_hdistance = new TH1F("m_hdistance", " distance between mathced extHit and bklmHit2d ", 100, 0, 30);
133 
134  m_totalMom = new TH1F("m_totalMom", " denominator vs. p", mNbin, mommin, mommax);
135  m_passMom = new TH1F("m_passMom", " numerator vs. p", mNbin, mommin, mommax);
136  m_effiMom = new TH1F("m_effiMom", " effi. vs. p", mNbin, mommin, mommax);
137  m_effiMom->GetXaxis()->SetTitle("p (GeV)");
138  m_effiMom->GetYaxis()->SetTitle("Efficiency (GeV)");
139 
140  m_totalYX = new TH2F("m_totalYX", " denominator Y vs. X", gNbin, gmin, gmax, gNbin, gmin, gmax);
141  m_passYX = new TH2F("m_passYX", " numerator Y vs. X", gNbin, gmin, gmax, gNbin, gmin, gmax);
142  m_totalYZ = new TH2F("m_totalYZ", " denominator Y vs. Z", gNbin, gmin, gmax, gNbin, gmin, gmax);
143  m_passYZ = new TH2F("m_passYZ", " numerator Y vs. Z", gNbin, gmin, gmax, gNbin, gmin, gmax);
144  m_effiYX = new TH2F("m_effiYX", " effi. Y vs. X", gNbin, gmin, gmax, gNbin, gmin, gmax);
145  m_effiYZ = new TH2F("m_effiYZ", " effi. Y vs. Z", gNbin, gmin, gmax, gNbin, gmin, gmax);
146  m_effiYX->GetXaxis()->SetTitle("x (cm)");
147  m_effiYX->GetYaxis()->SetTitle("y (cm)");
148  m_effiYZ->GetXaxis()->SetTitle("z (cm)");
149  m_effiYZ->GetYaxis()->SetTitle("y (cm)");
150 
151 }
152 
153 void BKLMAnaModule::beginRun()
154 {
155 }
156 
157 void BKLMAnaModule::event()
158 {
159  StoreObjPtr<EventMetaData> eventMetaData("EventMetaData", DataStore::c_Event);
160  //unsigned long eventNumber = eventMetaData->getEvent();
161  unsigned long runNumber = eventMetaData->getRun();
162  //unsigned long expNumber = eventMetaData->getExperiment();
163  //StoreArray<RecoTrack> recoTracks;
164  //numRecoTrk = recoTracks.getEntries();
165  //StoreArray<BKLMTrack> bklmtracks;
166  //StoreArray<TrackFitResult> trackFitResults;
167  //StoreArray<MuidHit> muidHits;
168  //StoreArray<Muid> muids;
169  //StoreArray<MCParticle> mcParticles;
170 
171  //set<int> m_pointUsed;
172  //m_pointUsed.clear();
173  //check extHits
174  //all ExtHit in bklm scope in each event, should not be many
175  int nExtHit = 0;
176  for (int t = 0; t < extHits.getEntries(); t++) {
177  ExtHit* exthit = extHits[t];
178  if (exthit->getDetectorID() != Const::EDetector::BKLM) continue;
179  m_extx[nExtHit] = exthit->getPosition()[0];
180  m_exty[nExtHit] = exthit->getPosition()[1];
181  m_extz[nExtHit] = exthit->getPosition()[2];
182  nExtHit++;
183  if (nExtHit > 199) break;
184  }
185  m_run = runNumber;
186  m_nExtHit = nExtHit;
187 
188 //the second way, require muid
189  for (int k = 0; k < tracks.getEntries(); k++) {
190  Track* track = tracks[k];
191  // load the muon fit hypothesis or the hypothesis which is the clostes in mass to a muon
192  // the tracking will not always fit a muon hypothesis
193  const TrackFitResult* fitres = track->getTrackFitResultWithClosestMass(Belle2::Const::muon);
194  double mom = fitres->getMomentum().Mag();
195  //double pt = fitres->getTransverseMomentum();
196  TLorentzVector p4 = fitres->get4Momentum();
197  double trkphi = p4.Vect().Phi() * 180.0 / CLHEP::pi;
198  double trktheta = p4.Vect().Theta() * 180.0 / CLHEP::pi;
199  if (trkphi < 0) trkphi = trkphi + 360.0;
200  RelationVector<BKLMHit2d> relatedHit2D = track->getRelationsTo<BKLMHit2d>();
201  RelationVector<ExtHit> relatedExtHit = track->getRelationsTo<ExtHit>();
202  for (unsigned int t = 0; t < relatedExtHit.size(); t++) {
203  ExtHit* exthit = relatedExtHit[t];
204  if (exthit->getDetectorID() != Const::EDetector::BKLM) continue;
205  int module = exthit->getCopyID();
206  int section = BKLMElementNumbers::getSectionByModule(module);
207  int sector = BKLMElementNumbers::getSectorByModule(module);
208  int layer = BKLMElementNumbers::getLayerByModule(module);
209  bool crossed = false; // should be only once ?
210  KLMMuidLikelihood* muid = track->getRelatedTo<KLMMuidLikelihood>();
211  if (!muid) continue;
212  int extPattern = muid->getExtLayerPattern();
213  if ((extPattern & (1 << (layer - 1))) != 0) crossed = true;
214  if (!crossed) continue;
215 
216  TVector3 extMom = exthit->getMomentum();
217  TVector3 extVec = exthit->getPosition();
218  bool matched = false;
219  m_totalYX->Fill(extVec[0], extVec[1]);
220  m_totalYZ->Fill(extVec[2], extVec[1]);
221  float phi = extVec.Phi() * 180.0 / CLHEP::pi;
222  float theta = extVec.Theta() * 180.0 / CLHEP::pi;
223  if (phi < 0) phi = phi + 360.0;
224  m_totalThephi[layer - 1]->Fill(phi, theta);
225  m_totalTrkThephi[layer - 1]->Fill(trkphi, trktheta);
226  m_totalMom->Fill(mom);
227  //look for mateched BKLM2dHit
228  //for (unsigned int mHit = 0; mHit < relatedHit2D.size(); mHit++) {
229  // BKLMHit2d* hit = relatedHit2D[mHit];
230  for (int mHit = 0; mHit < hits2D.getEntries(); mHit++) {
231  BKLMHit2d* hit = hits2D[mHit];
232  //if(!hit->inRPC()) continue;
233  if (hit->getSection() != section) continue;
234  if (hit->getSector() != sector) continue;
235  if (hit->getLayer() != layer) continue;
236  TVector3 position = hit->getGlobalPosition();
237  TVector3 distance = extVec - position;
238  //on same track, same sector, same layer, we should believe extHit and BKLMHit2d are matched.
239  //let's record the distance to check, should be small
240  m_hdistance->Fill(distance.Mag());
241  if (distance.Mag() < 20) matched = true;
242  //m_pointUsed.insert(m);
243  if (matched) break;
244  }
245  if (matched) {
246  m_passYX->Fill(extVec[0], extVec[1]);
247  m_passYZ->Fill(extVec[2], extVec[1]);
248  m_passTrkThephi[layer - 1]->Fill(trkphi, trktheta);
249  m_passThephi[layer - 1]->Fill(phi, theta);
250  m_passMom->Fill(mom);
251  }
252  }//end of loop ext hit
253  }//end of loop tracks
254  m_extTree->Fill();
255 }
256 
257 void BKLMAnaModule::endRun()
258 {
259 }
260 
261 void BKLMAnaModule::terminate()
262 {
263  for (int iL = 0; iL < 15; iL ++) {
264  for (int i = 0; i < m_totalThephi[iL]->GetNbinsX(); i++) {
265  for (int j = 0; j < m_totalThephi[iL]->GetNbinsY(); j++) {
266  float num = m_passThephi[iL]->GetBinContent(i + 1, j + 1);
267  float denom = m_totalThephi[iL]->GetBinContent(i + 1, j + 1);
268  if (num > 0) {
269  m_effiThephi[iL]->SetBinContent(i + 1, j + 1, num / denom);
270  m_effiThephi[iL]->SetBinError(i + 1, j + 1, sqrt(num * (denom - num) / (denom * denom * denom)));
271  } else {
272  m_effiThephi[iL]->SetBinContent(i + 1, j + 1, 0);
273  m_effiThephi[iL]->SetBinError(i + 1, j + 1, 0);
274  }
275  }
276  }
277  }
278 
279  for (int iL = 0; iL < 15; iL ++) {
280  for (int i = 0; i < m_totalTrkThephi[iL]->GetNbinsX(); i++) {
281  for (int j = 0; j < m_totalTrkThephi[iL]->GetNbinsY(); j++) {
282  float num = m_passTrkThephi[iL]->GetBinContent(i + 1, j + 1);
283  float denom = m_totalTrkThephi[iL]->GetBinContent(i + 1, j + 1);
284  if (num > 0) {
285  m_effiTrkThephi[iL]->SetBinContent(i + 1, j + 1, num / denom);
286  m_effiTrkThephi[iL]->SetBinError(i + 1, j + 1, sqrt(num * (denom - num) / (denom * denom * denom)));
287  } else {
288  m_effiTrkThephi[iL]->SetBinContent(i + 1, j + 1, 0);
289  m_effiTrkThephi[iL]->SetBinError(i + 1, j + 1, 0);
290  }
291  }
292  }
293  }
294 
295 
296  for (int i = 0; i < m_totalYX->GetNbinsX(); i++) {
297  for (int j = 0; j < m_totalYX->GetNbinsY(); j++) {
298  float num = m_passYX->GetBinContent(i + 1, j + 1);
299  float denom = m_totalYX->GetBinContent(i + 1, j + 1);
300  if (num > 0) {
301  m_effiYX->SetBinContent(i + 1, j + 1, num / denom);
302  m_effiYX->SetBinError(i + 1, j + 1, sqrt(num * (denom - num) / (denom * denom * denom)));
303  } else {
304  m_effiYX->SetBinContent(i + 1, j + 1, 0);
305  m_effiYX->SetBinError(i + 1, j + 1, 0);
306  }
307 
308  num = m_passYZ->GetBinContent(i + 1, j + 1);
309  denom = m_totalYZ->GetBinContent(i + 1, j + 1);
310  if (num > 0) {
311  m_effiYZ->SetBinContent(i + 1, j + 1, num / denom);
312  m_effiYZ->SetBinError(i + 1, j + 1, sqrt(num * (denom - num) / (denom * denom * denom)));
313  } else {
314  m_effiYZ->SetBinContent(i + 1, j + 1, 0);
315  m_effiYZ->SetBinError(i + 1, j + 1, 0);
316  }
317  }
318  }
319 
320  for (int i = 0; i < m_totalMom->GetNbinsX(); i++) {
321  float num = m_passMom->GetBinContent(i + 1);
322  float denom = m_totalMom->GetBinContent(i + 1);
323  if (num > 0) {
324  m_effiMom->SetBinContent(i + 1, num / denom);
325  m_effiMom->SetBinError(i + 1, sqrt(num * (denom - num) / (denom * denom * denom)));
326  } else {
327  m_effiMom->SetBinContent(i + 1, 0);
328  m_effiMom->SetBinError(i + 1, 0);
329  }
330  }
331 
332  m_file->cd();
333  m_hdistance->Write();
334  m_totalYX->Write();
335  m_passYX->Write();
336  m_totalYZ->Write();
337  m_passYZ->Write();
338  m_effiYX->Write();
339  m_effiYZ->Write();
340  m_totalMom->Write();
341  m_passMom->Write();
342  m_effiMom->Write();
343  for (int i = 0; i < 15; i++) {
344  m_totalThephi[i]->Write();
345  m_passThephi[i]->Write();
346  m_effiThephi[i]->Write();
347  m_totalTrkThephi[i]->Write();
348  m_passTrkThephi[i]->Write();
349  m_effiTrkThephi[i]->Write();
350  }
351  m_extTree->Write();
352  m_file->Close();
353 }
354 
Belle2::TrackFitResult::get4Momentum
TLorentzVector get4Momentum() const
Getter for the 4Momentum at the closest approach of the track in the r/phi projection.
Definition: TrackFitResult.h:125
Belle2::ExtHit::getPosition
TVector3 getPosition() const
Get position of this extrapolation hit.
Definition: ExtHit.h:153
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::ExtHit::getDetectorID
Const::EDetector getDetectorID() const
Get detector ID of this extrapolation hit.
Definition: ExtHit.h:130
Belle2::TrackFitResult
Values of the result of a track fit with a given particle hypothesis.
Definition: TrackFitResult.h:59
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::ExtHit
Store one Ext hit as a ROOT object.
Definition: ExtHit.h:40
Belle2::BKLMAnaModule
A class defining a module that perform efficiencies studies on bklm.
Definition: BKLMAnaModule.h:44
Belle2::RelationVector
Class for type safe access to objects that are referred to in relations.
Definition: DataStore.h:38
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::ExtHit::getCopyID
int getCopyID() const
Get detector-element ID of sensitive element within detector.
Definition: ExtHit.h:134
Belle2::ExtHit::getMomentum
TVector3 getMomentum() const
Get momentum at this extrapolation hit.
Definition: ExtHit.h:157
Belle2::Const::muon
static const ChargedStable muon
muon particle
Definition: Const.h:534
Belle2::Track
Class that bundles various TrackFitResults.
Definition: Track.h:35
Belle2::KLMMuidLikelihood
Class to store the likelihoods from KLM with additional informations related to the extrapolation.
Definition: KLMMuidLikelihood.h:38
Belle2::BKLMHit2d
Store one BKLM strip hit as a ROOT object.
Definition: BKLMHit2d.h:42