Belle II Software development
MCTrackCandClassifierModule.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 <tracking/modules/mcTrackCandClassifier/MCTrackCandClassifierModule.h>
10
11#include <pxd/dataobjects/PXDTrueHit.h>
12#include <svd/dataobjects/SVDTrueHit.h>
13#include <vxd/geometry/GeoCache.h>
14
15#include <framework/geometry/B2Vector3.h>
16#include <framework/geometry/BFieldManager.h>
17
18#include <TH2F.h>
19
20#include <cmath>
21
22using namespace Belle2;
23
25REG_MODULE(MCTrackCandClassifier);
26
29 , m_rootFilePtr(nullptr)
30{
31 //Set module properties
32 setDescription("This module is meant to classify the MCTrackCands as either ideal, fine and nasty");
34
35 addParam("MCParticlesName", m_mcParticlesName,
36 "Name of MC Particle collection.",
37 std::string(""));
38
39 addParam("MCTrackCandCollName", m_mcTrackCandsColName,
40 "Name of the input collection of MC track candidates",
41 std::string(""));
42
43 addParam("rootFileName", m_rootFileName,
44 "Name of the root file",
45 std::string("MCTrackCandClassifier.root"));
46
47 addParam("isInAnnulusCriterium", m_applyAnnulus,
48 "Require that the hit is in the expected annulus",
49 bool(true));
50 addParam("isInSemiplaneCriterium", m_applySemiplane,
51 "Require that the hit is in the expected semiplane",
52 bool(true));
53 addParam("isInFirstLapCriterium", m_applyLap,
54 "Require that the hit belong to the first lap in the transverse plane",
55 bool(true));
56 addParam("isInWedgePartCriterium", m_applyWedge,
57 "Require that the hit belong to the barrel part of the SVD",
58 bool(true));
59 addParam("removeBadHits", m_removeBadHits,
60 "Remove the clusters that do not satisfy the criteria from the idealMCTrackCands",
61 bool(true));
62
63 addParam("minNhits", m_minHit,
64 "Minimum number of 1D Clusters to classify the MCTrackCand as ideal",
65 int(5));
66
67 addParam("nSigma_dR", m_nSigma, "n sigma dR", int(3));
68
69 addParam("lapFraction", m_fraction, "Fraction of lap", double(1));
70
71 addParam("usePXD", m_usePXD, "Use the PXD or not", bool(true));
72}
73
74
76{
77 // MCParticles, MCTrackCands, MCTracks needed for this module
78 m_PXDClusters.isRequired();
79 m_SVDClusters.isRequired();
81
83 m_GenfitIdealMCTrackCands.registerInDataStore("idealMCTrackCands", DataStore::c_ErrorIfAlreadyRegistered);
84
85 //create list of histograms to be saved in the rootfile
86 m_histoList = new TList;
87
88 //set the ROOT File
89 m_rootFilePtr = new TFile(m_rootFileName.c_str(), "RECREATE");
90
91 //histograms to produce efficiency plots
92 Double_t bins_pt[9 + 1] = {0, 0.05, 0.1, 0.15, 0.2, 0.3, 0.5, 1, 2, 3.5}; //GeV/c
93 const Double_t bins_theta[10 + 1] = {0, 0.25, 0.5, 0.75, 0.75 + 0.32, 0.75 + 2 * 0.32, 0.75 + 3 * 0.32, 0.75 + 4 * 0.32, 0.75 + 5 * 0.32, 2.65, TMath::Pi()};
94 Double_t bins_lambda[10 + 1];
95 const Double_t width_lambda = TMath::Pi() / 10;
96 Double_t bins_phi[14 + 1];
97 const Double_t width_phi = 2 * TMath::Pi() / 14;
98 for (int bin = 0; bin < 14 + 1; bin++)
99 bins_phi[bin] = - TMath::Pi() + bin * width_phi;
100
101 for (int bin = 0; bin < 10 + 1; bin++) {
102 bins_lambda[bin] = - TMath::Pi() / 2 + bin * width_lambda;
103 B2DEBUG(21, bins_lambda[bin] << " " << bins_theta[bin]);
104 }
105
106 m_h3_MCParticle = createHistogram3D("h3MCParticle", "entry per MCParticle",
107 9, bins_pt, "p_{t} (GeV/c)",
108 10, bins_lambda, "#lambda",
109 14, bins_phi, "#phi" /*, m_histoList*/);
110
111 m_h3_idealMCTrackCand = (TH3F*)duplicateHistogram("h3idealMCTrackCand",
112 "entry per idealMCTrackCand",
113 m_h3_MCParticle /*, m_histoList*/);
114
115 m_h3_MCTrackCand = (TH3F*)duplicateHistogram("h3MCTrackCand",
116 "entry per MCTrackCand",
117 m_h3_MCParticle /*, m_histoList*/);
118
119
120 m_h1_thetaMS_SVD = new TH1F("h1thetaMS_SVD", "Multiple Scattering Angle (SVD)", 500, 0, 500);
122 m_h1_thetaMS_SVD->GetXaxis()->SetTitle("#theta_{MS} (mrad)");
123
124 m_h1_thetaMS_PXD = (TH1F*) duplicateHistogram("h1thetaMS_PXD", "Multiple Scattering Angle (PXD)", m_h1_thetaMS_SVD, m_histoList);
125
126 m_h1_dR = new TH1F("h1dR", "dR, annulus half width", 1000, 0, 5);
127 m_histoList->Add(m_h1_dR);
128 m_h1_dR->GetXaxis()->SetTitle("dR (cm)");
129
130 m_h1_dRoverR = new TH1F("h1dRoverR", "dR over helix radius", 1000, 0, 0.1);
132 m_h1_dRoverR->GetXaxis()->SetTitle("dR/R");
133
134 m_h1_distOVERdR = new TH1F("h1distOVERdR", "(hit radius - helix radius)/dR", 100, -5, 5);
136 m_h1_distOVERdR->GetXaxis()->SetTitle("(hit R - helix R)/dR");
137
138 m_h1_hitRadius_accepted = new TH1F("h1hitRadAccep", "hit radius for accepted hits", 100, 0, 500);
140 m_h1_hitRadius_accepted->GetXaxis()->SetTitle("hit radius (cm)");
141
142 m_h1_hitRadius_rejected = new TH1F("h1hitRadRejec", "hit radius for rejected hits", 100, 0, 500);
144 m_h1_hitRadius_rejected->GetXaxis()->SetTitle("hit radius (cm)");
145
146 m_h1_hitDistance_accepted = new TH1F("h1hitDistCentAccep", "hit distance from 0,0 for accepted hits", 100, 0, 15);
148 m_h1_hitDistance_accepted->GetXaxis()->SetTitle("hit distance (cm)");
149
150 m_h1_hitDistance_rejected = new TH1F("h1hitDistCentRejec", "hit distance from 0,0 for rejected hits", 100, 0, 15);
152 m_h1_hitDistance_rejected->GetXaxis()->SetTitle("hit distance (cm)");
153
154 m_h1_lapTime = new TH1F("h1LapTime", "lap time", 200, 0, 100);
156 m_h1_lapTime->GetXaxis()->SetTitle("time (ns)");
157
158 m_h1_timeDifference = (TH1F*)duplicateHistogram("h1TimeDiff", "Hit Time Difference",
160 m_h1_diffOVERlap = new TH1F("h1HitDiffOVERlap", "Hit Time Difference over Lap Time",
161 100, 0, 1.5);
163 m_h1_diffOVERlap->GetXaxis()->SetTitle("hit difference/lap");
164
165 m_h1_nGoodTrueHits = new TH1F("h1nTrueHitsGoods", "Number of True Hits for Accepted Tracks", 20, 0, 20);
167 m_h1_nGoodTrueHits->GetXaxis()->SetTitle("number of hits");
168
169 m_h1_nBadTrueHits = new TH1F("h1nTrueHitsBads", "Number of True Hits for Rejected Tracks", 10, 0, 10);
171 m_h1_nBadTrueHits->GetXaxis()->SetTitle("number of hits");
172
173 m_h1_nGood1dInfo = new TH1F("h1nGood1Dinfo", "Number of 1D Info for Accepted Tracks", 20, 0, 20);
175 m_h1_nGood1dInfo->GetXaxis()->SetTitle("number of hits");
176
177 m_h1_nBad1dInfo = new TH1F("h1nBad1Dinfo", "Number of 1D Info for Rejected Tracks", 20, 0, 20);
179 m_h1_nBad1dInfo->GetXaxis()->SetTitle("number of hits");
180
181 m_h1_firstRejectedHit = new TH1F("h1idealMCTCnHit", "idealMCTrackCands number of hits", 40, 0, 40);
183 m_h1_firstRejectedHit->GetXaxis()->SetTitle("# idealMCTrackCands hits");
184
185 m_h1_firstRejectedOVERMCHit = new TH1F("h1FirstRejOVERmc", "# idealMCTrackCands hits / # MCTrackCands hits", 100, 0, 1);
187 m_h1_firstRejectedOVERMCHit->GetXaxis()->SetTitle("# idealMCTrackCands hits / # MCTrackCands hits");
188
189 m_h1_MCTrackCandNhits = (TH1F*)duplicateHistogram("h1MCTrackCandNhits", "number of MCTrackCands hits", m_h1_firstRejectedHit,
191}
192
193
195{
196 nWedge = 0;
197 nBarrel = 0;
198
200}
201
202
204{
205 B2DEBUG(21, "+++++ 1. loop on MCTrackCands");
206
207 const VXD::GeoCache& aGeometry = VXD::GeoCache::getInstance();
208
209 //1.a retrieve the MCTrackCands
210 for (const genfit::TrackCand& mcTrackCand : m_GenfitMCTrackCands) {
211
212 int nGoodTrueHits = 0;
213 int nGood1Dinfo = 0;
214
215 B2DEBUG(21, " a NEW MCTrackCand ");
216
217 //1.b retrieve the MCParticle
218 RelationVector<MCParticle> MCParticles_fromMCTrackCand = DataStore::getRelationsWithObj<MCParticle>(&mcTrackCand);
219
220 B2DEBUG(21, "~~~ " << MCParticles_fromMCTrackCand.size() << " MCParticles related to this MCTrackCand");
221 for (int mcp = 0; mcp < (int)MCParticles_fromMCTrackCand.size(); mcp++) { //should be ONE
222
223 MCParticle mcParticle = *MCParticles_fromMCTrackCand[mcp];
224
225 B2DEBUG(21, " a NEW charged MC Particle, " << mcParticle.getIndex() << ", " << mcParticle.getPDG());
226
227 MCParticleInfo mcParticleInfo(mcParticle, m_magField);
228
229 B2Vector3D decayVertex = mcParticle.getProductionVertex();
230 ROOT::Math::XYZVector mom = mcParticle.getMomentum();
231 double charge = mcParticle.getCharge();
232 double omega = mcParticleInfo.getOmega();
233 double px = mom.x();
234 double py = mom.y();
235 double pt = mom.Rho();
236 double x = decayVertex.X();
237 double y = decayVertex.Y();
238 double R = 1 / std::abs(omega); //cm
239
240 m_h3_MCTrackCand->Fill(mcParticleInfo.getPt(), mcParticleInfo.getLambda(), mcParticleInfo.getPphi());
241
242 double alpha = R / pt * charge; //cm/GeV
243 double Cx = x + alpha * py; //cm
244 double Cy = y - alpha * px; //cm
245
246 ROOT::Math::XYZVector center(Cx, Cy, 0);
247
248 //recover Clusters and loop on them
249 int Nhits = mcTrackCand.getNHits();
250 m_h1_MCTrackCandNhits->Fill(Nhits);
251
252 int cluster = 0;
253 bool hasTrueHit = true;
254 bool isAccepted = true;
255 int firstRejectedHit = Nhits + 1;
256 double prevHitRadius = std::abs(1 / omega);
257
258 double lapTime = 2 * M_PI * mcParticle.getEnergy() / 0.299792 / m_magField.Z();
259 double FirstHitTime = -1;
260 double HitTime = -1;
261
262 bool isFirstSVDhit = true;
263
264 while (cluster < Nhits && isAccepted && hasTrueHit) {
265 int detId, hitId;
266 mcTrackCand.getHit(cluster, detId, hitId);
267
268 bool hasPXDCluster = false;
269 bool hasSVDuCluster = false;
270 bool hasSVDvCluster = false;
271
272 double uCoor = 0;
273 double vCoor = 0;
274 VxdID sensor = 0;
275
276 double thetaMS = 0;
277
278 if (detId == Const::PXD && m_usePXD) {
279
280 PXDCluster* aPXDCluster = m_PXDClusters[hitId];
281 RelationVector<PXDTrueHit> PXDTrueHit_fromPXDCluster = aPXDCluster->getRelationsWith<PXDTrueHit>();
282 if (PXDTrueHit_fromPXDCluster.size() == 0) {
283 B2WARNING("What's happening?!? no True Hit associated to the PXD Cluster");
284 hasTrueHit = false;
285 isAccepted = false;
286 continue;
287 }
288
289 PXDTrueHit* aPXDTrueHit = PXDTrueHit_fromPXDCluster[0];
290 thetaMS = compute_thetaMS(mcParticleInfo, aPXDTrueHit);
291 m_h1_thetaMS_PXD->Fill(thetaMS / 2 * 1000); //PXD
292
293 uCoor = aPXDTrueHit->getU();
294 vCoor = aPXDTrueHit->getV();
295 sensor = aPXDTrueHit->getSensorID();
296 if (cluster == 0) {
297 FirstHitTime = aPXDTrueHit->getGlobalTime();
298 HitTime = FirstHitTime;
299 } else
300 HitTime = aPXDTrueHit->getGlobalTime();
301
302 hasPXDCluster = true;
303 } else if (detId == Const::SVD) {
304 SVDCluster* aSVDCluster = m_SVDClusters[hitId];
305 RelationVector<SVDTrueHit> SVDTrueHit_fromSVDCluster = aSVDCluster->getRelationsWith<SVDTrueHit>();
306 if (SVDTrueHit_fromSVDCluster.size() == 0) {
307 B2WARNING("What's happening?!? no True Hit associated to the SVD Cluster");
308 hasTrueHit = false;
309 isAccepted = false;
310 continue;
311 }
312
313 SVDTrueHit* aSVDTrueHit = SVDTrueHit_fromSVDCluster[0];
314
315 thetaMS = compute_thetaMS(mcParticleInfo, aSVDTrueHit);
316 m_h1_thetaMS_SVD->Fill(thetaMS * 1000); //SVD
317
318 uCoor = aSVDTrueHit->getU();
319 vCoor = aSVDTrueHit->getV();
320 sensor = aSVDTrueHit->getSensorID();
321 if (isFirstSVDhit) {
322 FirstHitTime = aSVDTrueHit->getGlobalTime();
323 HitTime = FirstHitTime;
324 isFirstSVDhit = false;
325
326 } else
327 HitTime = aSVDTrueHit->getGlobalTime();
328 if (aSVDCluster->isUCluster())
329 hasSVDuCluster = true;
330 else
331 hasSVDvCluster = true;
332 } else {
333 cluster++;
334 continue;
335 }
336
337 const VXD::SensorInfoBase& aSensorInfo = aGeometry.getSensorInfo(sensor);
338 bool accepted4 = true;
339 if (m_applyWedge) {
340 if (aSensorInfo.getForwardWidth() != aSensorInfo.getBackwardWidth()) {
341 nWedge++;
342 accepted4 = false;
343 } else
344 nBarrel++;
345 }
346
347 ROOT::Math::XYZVector globalHit = aSensorInfo.pointToGlobal(ROOT::Math::XYZVector(uCoor, vCoor, 0), true);
348 double hitRadius = theDistance(center, globalHit);
349
350 bool accepted1 = true;
352 accepted1 = isInSemiPlane(semiPlane(decayVertex, center, globalHit), omega);
353
354 if (accepted1) {
355 B2DEBUG(21, " semiplane: ACCEPTED");
356 } else {
357 B2DEBUG(21, " semiplane: REJECTED, next track");
358 }
359
360 double dR = compute_dR(thetaMS, theDistance(ROOT::Math::XYZVector(0, 0, 0), globalHit));
361 m_h1_dR->Fill(dR);
362 m_h1_dRoverR->Fill(dR * std::abs(omega));
363 m_h1_distOVERdR->Fill((hitRadius - std::abs(1 / omega)) / dR);
364
365
366 bool accepted2 = true;
367 if (m_applyAnnulus)
368 accepted2 = isInAnnulus(hitRadius, prevHitRadius, dR);
369
370 prevHitRadius = hitRadius;
371
372 if (accepted2) {
373 B2DEBUG(21, " annulus: ACCEPTED");
374 } else {
375 B2DEBUG(21, " annulus: REJECTED, next track");
376 }
377
378 bool accepted3 = true;
379 if (m_applyLap)
380 accepted3 = isFirstLap(FirstHitTime, HitTime, lapTime);
381
382 if (accepted3) {
383 B2DEBUG(21, " lapTime: ACCEPTED");
384 } else {
385 B2DEBUG(21, " lapTime: REJECTED, next track");
386 }
387
388 if (accepted2 && accepted1 && accepted3 && accepted4) {
389 nGoodTrueHits ++;
390 m_h1_hitDistance_accepted->Fill(theDistance(ROOT::Math::XYZVector(0, 0, 0), globalHit));
391 m_h1_hitRadius_accepted->Fill(hitRadius);
392 } else {
393 m_h1_hitDistance_rejected->Fill(theDistance(ROOT::Math::XYZVector(0, 0, 0), globalHit));
394 m_h1_hitRadius_rejected->Fill(hitRadius);
395 if (m_removeBadHits)
396 firstRejectedHit = cluster;
397 isAccepted = false;
398 continue;
399 }
400
401 if (hasPXDCluster)
402 nGood1Dinfo = +2;
403 else {
404 if (hasSVDuCluster)
405 nGood1Dinfo++;
406 if (hasSVDvCluster)
407 nGood1Dinfo++;
408 }
409 if (hasPXDCluster || hasSVDuCluster || hasSVDvCluster)
410 B2DEBUG(21, "cluster: ACCEPTED (" << nGood1Dinfo << ")");
411
412 cluster++;
413 }//close loop on clusters
414
415 if (nGood1Dinfo >= m_minHit) {
416 B2DEBUG(21, " idealMCTrackCand FOUND!! " << nGood1Dinfo << " 1D infos (" << nGoodTrueHits << " good true hits)");
417 m_h3_idealMCTrackCand->Fill(mcParticleInfo.getPt(), mcParticleInfo.getLambda(), mcParticleInfo.getPphi());
418 m_h1_nGoodTrueHits->Fill(nGoodTrueHits);
419 m_h1_nGood1dInfo->Fill(nGood1Dinfo);
420
421 genfit::TrackCand* tmpTrackCand = new genfit::TrackCand(mcTrackCand);
422
423 if ((int)firstRejectedHit <= (int)mcTrackCand.getNHits()) {
424 tmpTrackCand->reset();
425 for (int hit = 0; hit < firstRejectedHit; hit++)
426 if (mcTrackCand.getHit(hit))
427 tmpTrackCand->addHit(mcTrackCand.getHit(hit));
428 tmpTrackCand->sortHits();
429 }
430 m_GenfitIdealMCTrackCands.appendNew(*tmpTrackCand);
431
432 m_h1_firstRejectedHit->Fill(tmpTrackCand->getNHits());
433 m_h1_firstRejectedOVERMCHit->Fill((float)tmpTrackCand->getNHits() / mcTrackCand.getNHits());
434 } else {
435 B2DEBUG(21, " too few good hits (" << nGood1Dinfo << ") to track this one ( vs " << nGoodTrueHits << " true hits)");
436 m_h1_nBadTrueHits->Fill(nGoodTrueHits);
437 m_h1_nBad1dInfo->Fill(nGood1Dinfo);
438 }
439
440 B2DEBUG(21, "");
441 }//close loop on MCParticles
442 }//close loop on MCTrackCands
443}
444
445
447{
448 B2INFO("** MCTrackCandClassifier parameters **");
449 B2INFO("rootfilename = " << m_rootFileName);
450 B2INFO("use PXD information = " << m_usePXD);
451 B2INFO("--> classification criteria:");
452 if (m_applyAnnulus)
453 B2INFO(" -) |d - R| < " << m_nSigma << " dL thetaMS");
455 B2INFO(" -) hit in the expected semiplane");
456 if (m_applyLap)
457 B2INFO(" -) HitTime < " << m_fraction << " lap time");
458 if (m_applyWedge)
459 B2INFO(" -) hit must be in the barrel part of the VXD");
460 B2INFO("");
461
462 double num = 0;
463 double den = 0;
464
465 num = m_h3_idealMCTrackCand->GetEntries();
466 den = m_h3_MCTrackCand->GetEntries();
467 double efficiency = num / den ;
468 double efficiencyErr = sqrt(efficiency * (1 - efficiency)) / sqrt(den);
469
470 B2INFO("");
471 B2INFO("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");
472 B2INFO("~ MCTrackCandClassifier ~ SHORT SUMMARY ~");
473 B2INFO("");
474 B2INFO(" + overall:");
475 B2INFO(" fraction of ideal MCTrackCands = (" << efficiency * 100 << " +/- " << efficiencyErr * 100 << ")% ");
476 B2INFO("");
477 B2INFO(" # idealMCTrackCand = " << num);
478 B2INFO(" # MCTrackCand = " << den);
479 B2INFO("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");
480 B2INFO("");
481 B2DEBUG(21, " nWedge = " << nWedge);
482 B2DEBUG(21, " nBarrel = " << nBarrel);
483}
484
485
487{
490
491 if (m_rootFilePtr != nullptr) {
492 m_rootFilePtr->cd();
493
494 TIter nextH(m_histoList);
495 TObject* obj;
496 while ((obj = nextH()))
497 obj->Write();
498
499 m_rootFilePtr->Close();
500 }
501}
502
503
504double MCTrackCandClassifierModule::semiPlane(ROOT::Math::XYZVector vertex, ROOT::Math::XYZVector center, ROOT::Math::XYZVector hit)
505{
506 ROOT::Math::XYZVector err = center - vertex;
507
508 double semiPlane = err.Y() / err.X() * hit.X() + err.Y() / err.X() * vertex.x() - vertex.Y();
509
510 B2DEBUG(21, "");
511 B2DEBUG(21, " SEMI-PLANE defined by: y + " << err.Y() / err.X() << " x + " << err.Y() / err.X()*vertex.x() - vertex.Y() << " = 0");
512 B2DEBUG(21, " with: center(" << center.X() << "," << center.Y() << ")");
513 B2DEBUG(21, " decayV(" << vertex.X() << "," << vertex.Y() << ")");
514 B2DEBUG(21, " vector(" << err.X() << "," << err.Y() << ")");
515 B2DEBUG(21, " y SLOPE = " << semiPlane << " VS y HIT = " << hit.Y());
516 B2DEBUG(21, " HIT - SLOPE = " << - semiPlane + hit.Y());
517
518 if (vertex.X() < center.X())
519 return hit.Y() - semiPlane;
520 else
521 return semiPlane - hit.Y();
522}
523
524
526{
527 if (semiPlane * omega > 0)
528 return true;
529 else
530 return false;
531}
532
533
534double MCTrackCandClassifierModule::theDistance(ROOT::Math::XYZVector center, ROOT::Math::XYZVector hit)
535{
536 double xSquared = TMath::Power(center.X() - hit.X(), 2);
537 double ySquared = TMath::Power(center.Y() - hit.Y(), 2);
538
539 return TMath::Sqrt(xSquared + ySquared);
540}
541
542
543bool MCTrackCandClassifierModule::isInAnnulus(double hitDistance, double R, double dR)
544{
545 bool accepted = false;
546
547 B2DEBUG(21, "");
548 B2DEBUG(21, " ANNULUS defined between radii: " << R - dR << " and " << R + dR);
549 B2DEBUG(21, " hit distance = " << hitDistance);
550 B2DEBUG(21, " helix radius = " << R);
551 B2DEBUG(21, " dR = " << dR);
552
553 if ((hitDistance > R - dR) && (hitDistance < R + dR))
554 accepted = true;
555
556 return accepted;
557}
558
559
560bool MCTrackCandClassifierModule::isFirstLap(double FirstHitTime, double HitTime, double LapTime)
561{
562 bool accepted = false;
563
564 B2DEBUG(21, "");
565 B2DEBUG(21, " lapTime: " << LapTime);
566 B2DEBUG(21, " FirstHitTime = " << FirstHitTime);
567 B2DEBUG(21, " HitTime = " << HitTime);
568 B2DEBUG(21, " difference = " << HitTime - FirstHitTime);
569
570 m_h1_lapTime->Fill(LapTime);
571 m_h1_timeDifference->Fill(HitTime - FirstHitTime);
572 m_h1_diffOVERlap->Fill((HitTime - FirstHitTime) / LapTime);
573
574 if (HitTime - FirstHitTime < m_fraction * LapTime)
575 accepted = true;
576
577 return accepted;
578}
579
580
581TH3F* MCTrackCandClassifierModule::createHistogram3D(const char* name, const char* title,
582 Int_t nbinsX, Double_t minX, Double_t maxX,
583 const char* titleX,
584 Int_t nbinsY, Double_t minY, Double_t maxY,
585 const char* titleY,
586 Int_t nbinsZ, Double_t minZ, Double_t maxZ,
587 const char* titleZ,
588 TList* histoList)
589{
590 TH3F* h = new TH3F(name, title, nbinsX, minX, maxX, nbinsY, minY, maxY, nbinsZ, minZ, maxZ);
591
592 h->GetXaxis()->SetTitle(titleX);
593 h->GetYaxis()->SetTitle(titleY);
594 h->GetZaxis()->SetTitle(titleZ);
595
596 if (histoList)
597 histoList->Add(h);
598
599 return h;
600}
601
602
603TH3F* MCTrackCandClassifierModule::createHistogram3D(const char* name, const char* title,
604 Int_t nbinsX, Double_t* binsX,
605 const char* titleX,
606 Int_t nbinsY, Double_t* binsY,
607 const char* titleY,
608 Int_t nbinsZ, Double_t* binsZ,
609 const char* titleZ,
610 TList* histoList)
611{
612 TH3F* h = new TH3F(name, title, nbinsX, binsX, nbinsY, binsY, nbinsZ, binsZ);
613
614 h->GetXaxis()->SetTitle(titleX);
615 h->GetYaxis()->SetTitle(titleY);
616 h->GetZaxis()->SetTitle(titleZ);
617
618 if (histoList)
619 histoList->Add(h);
620
621 return h;
622}
623
624
625TH1* MCTrackCandClassifierModule::duplicateHistogram(const char* newname, const char* newtitle,
626 TH1* h, TList* histoList)
627{
628 TH1F* h1 = dynamic_cast<TH1F*>(h);
629 TH2F* h2 = dynamic_cast<TH2F*>(h);
630 TH3F* h3 = dynamic_cast<TH3F*>(h);
631
632 TH1* newh = nullptr;
633
634 if (h1)
635 newh = new TH1F(*h1);
636 if (h2)
637 newh = new TH2F(*h2);
638 if (h3)
639 newh = new TH3F(*h3);
640
641 if (newh == nullptr) {
642 B2ERROR("In function duplicateHistogram: newh is a nullptr. This shouldn't happen."\
643 "Don't continue creation of duplicate histogram in this case and return nullptr.");
644 return nullptr;
645 }
646
647 newh->SetName(newname);
648 newh->SetTitle(newtitle);
649
650 if (histoList)
651 histoList->Add(newh);
652
653 return newh;
654}
655
656
658{
659 //normalized to MCTrackCands
660 TH1F* h_effMCTC_pt = createHistogramsRatio("heffMCTCpt", "fraction of idealMCTrackCand VS pt", m_h3_idealMCTrackCand,
661 m_h3_MCTrackCand, true, 0);
662 histoList->Add(h_effMCTC_pt);
663
664 TH1F* h_effMCTC_theta = createHistogramsRatio("heffMCTCtheta", "fraction of idealMCTrackCandVS #lambda", m_h3_idealMCTrackCand,
665 m_h3_MCTrackCand, true, 1);
666 histoList->Add(h_effMCTC_theta);
667
668 TH1F* h_effMCTC_phi = createHistogramsRatio("heffMCTCphi", "fraction of idealMCTrackCand VS #phi", m_h3_idealMCTrackCand,
669 m_h3_MCTrackCand, true, 2);
670 histoList->Add(h_effMCTC_phi);
671}
672
673
675{
676 //normalized to MCTrackCands
677 TH1F* h_ineffMCTC_pt = createHistogramsRatio("hineffMCTCpt", "1 - fraction of idealMCTrackCand VS pt", m_h3_idealMCTrackCand,
678 m_h3_MCTrackCand, false, 0);
679 histoList->Add(h_ineffMCTC_pt);
680
681 TH1F* h_ineffMCTC_theta = createHistogramsRatio("hineffMCTCtheta", "1 - fraction of idealMCTrackCandVS #lambda",
683 histoList->Add(h_ineffMCTC_theta);
684
685 TH1F* h_ineffMCTC_phi = createHistogramsRatio("hineffMCTCphi", "1 - fraction of idealMCTrackCand VS #phi", m_h3_idealMCTrackCand,
686 m_h3_MCTrackCand, false, 2);
687 histoList->Add(h_ineffMCTC_phi);
688}
689
690
691TH1F* MCTrackCandClassifierModule::createHistogramsRatio(const char* name, const char* title,
692 TH1* hNum, TH1* hDen, bool isEffPlot,
693 int axisRef)
694{
695 TH1F* h1den = dynamic_cast<TH1F*>(hDen);
696 TH1F* h1num = dynamic_cast<TH1F*>(hNum);
697 TH2F* h2den = dynamic_cast<TH2F*>(hDen);
698 TH2F* h2num = dynamic_cast<TH2F*>(hNum);
699 TH3F* h3den = dynamic_cast<TH3F*>(hDen);
700 TH3F* h3num = dynamic_cast<TH3F*>(hNum);
701
702 TH1* hden = nullptr;
703 TH1* hnum = nullptr;
704
705 if (h1den) {
706 hden = new TH1F(*h1den);
707 hnum = new TH1F(*h1num);
708 }
709 if (h2den) {
710 hden = new TH2F(*h2den);
711 hnum = new TH2F(*h2num);
712 }
713 if (h3den) {
714 hden = new TH3F(*h3den);
715 hnum = new TH3F(*h3num);
716 }
717
718 if (hden == nullptr or hnum == nullptr) {
719 B2ERROR("In function createHistogramsRatio: either hden or hnum are a nullptr. This shouldn't happen."\
720 "Don't continue creatio of histogram ratios in this case and return nullptr.");
721 return nullptr;
722 }
723
724 TAxis* the_axis;
725 TAxis* the_other1;
726 TAxis* the_other2;
727
728 if (axisRef == 0) {
729 the_axis = hden->GetXaxis();
730 the_other1 = hden->GetYaxis();
731 the_other2 = hden->GetZaxis();
732 } else if (axisRef == 1) {
733 the_axis = hden->GetYaxis();
734 the_other1 = hden->GetXaxis();
735 the_other2 = hden->GetZaxis();
736 } else if (axisRef == 2) {
737 the_axis = hden->GetZaxis();
738 the_other1 = hden->GetXaxis();
739 the_other2 = hden->GetYaxis();
740 } else
741 return nullptr;
742
743
744 TH1F* h;
745 if (the_axis->GetXbins()->GetSize())
746 h = new TH1F(name, title, the_axis->GetNbins(), (the_axis->GetXbins())->GetArray());
747 else
748 h = new TH1F(name, title, the_axis->GetNbins(), the_axis->GetXmin(), the_axis->GetXmax());
749 h->GetXaxis()->SetTitle(the_axis->GetTitle());
750
751 h->GetYaxis()->SetRangeUser(0.00001, 1);
752
753 Int_t bin = 0;
754
755 for (int the_bin = 1; the_bin < the_axis->GetNbins() + 1; the_bin++) {
756
757 double num = 0;
758 double den = 0 ;
759
760 for (int other1_bin = 1; other1_bin < the_other1->GetNbins() + 1; other1_bin++)
761 for (int other2_bin = 1; other2_bin < the_other2->GetNbins() + 1; other2_bin++) {
762
763 if (axisRef == 0) bin = hden->GetBin(the_bin, other1_bin, other2_bin);
764 else if (axisRef == 1) bin = hden->GetBin(other1_bin, the_bin, other2_bin);
765 else if (axisRef == 2) bin = hden->GetBin(other1_bin, other2_bin, the_bin);
766
767 if (hden->IsBinUnderflow(bin))
768 B2DEBUG(21, " bin = " << bin << "(" << the_bin << "," << other1_bin << "," << other2_bin << "), UNDERFLOW");
769 if (hden->IsBinOverflow(bin))
770 B2DEBUG(21, " bin = " << bin << "(" << the_bin << "," << other1_bin << "," << other2_bin << "), OVERFLOW");
771
772 num += hnum->GetBinContent(bin);
773 den += hden->GetBinContent(bin);
774 }
775 double eff = 0;
776 double err = 0;
777
778 if (den > 0) {
779 eff = (double)num / den;
780 err = sqrt(eff * (1 - eff)) / sqrt(den);
781 }
782
783 if (isEffPlot) {
784 h->SetBinContent(the_bin, eff);
785 h->SetBinError(the_bin, err);
786 } else {
787 h->SetBinContent(the_bin, 1 - eff);
788 h->SetBinError(the_bin, err);
789 }
790
791 }
792
793 return h;
794}
795
796
797float MCTrackCandClassifierModule::compute_dR(double thetaMS, double hitDistance)
798{
799 double dL;
800 if (hitDistance < 1.8) //L1
801 dL = 0.4;
802 else if (hitDistance < 3) //L2
803 dL = 0.8;
804 else if (hitDistance < 5) //L3
805 dL = 1.6;
806 else if (hitDistance < 9) //L4
807 dL = 4.2;
808 else if (hitDistance < 12) //L5
809 dL = 2.4;
810 else dL = 3.1;
811
812 if ((hitDistance < 3) && (hitDistance > 1.2))
813 thetaMS = thetaMS / 2;
814
815 double dR = m_nSigma * dL * thetaMS;
816
817 return dR;
818};
819
820
822{
823 // double thetaMS = 0.0136 * 14 * sqrt(0.008); //SVD, PXD is half of it
824 double thetaMS = 0.0136 * 14; //SVD, PXD is half of it
825
826 double p = mcParticleInfo.getP();
827 // double pt = mcParticleInfo.getPt();
828 double E = mcParticleInfo.getEnergy();
829
830 double X = sqrt(pow(aTrueHit->getEntryU() - aTrueHit->getExitU(), 2) +
831 pow(aTrueHit->getEntryV() - aTrueHit->getExitV(), 2) +
832 pow(aTrueHit->getEntryW() - aTrueHit->getExitW(), 2));
833
834 double X0 = 21.82; // g cm-2
835 double rho = 2.329; // g cm-3
836 thetaMS = thetaMS / (p * p / E) * sqrt(X / X0 * rho);
837
838 return thetaMS;
839};
R E
internal precision of FFTW codelets
double R
typedef autogenerated by FFTW
DataType X() const
access variable X (= .at(0) without boundary check)
Definition B2Vector3.h:431
DataType Y() const
access variable Y (= .at(1) without boundary check)
Definition B2Vector3.h:433
@ c_ErrorIfAlreadyRegistered
If the object/array was already registered, produce an error (aborting initialisation).
Definition DataStore.h:72
static RelationVector< T > getRelationsWithObj(const TObject *object, const std::string &name="", const std::string &namedRelation="")
Get the relations between an object and other objects in a store array.
Definition DataStore.h:412
This struct is used by the TrackingPerformanceEvaluation Module to save information of reconstructed ...
double getPt()
Getter for transverse momentum.
double getLambda()
Getter for Lambda.
double getEnergy()
Getter for energy.
double getOmega()
Getter for Omega.
double getPphi()
Getter for phi of momentum vector.
double getP()
Getter for magnitut of momentum.
A Class to store the Monte Carlo particle information.
Definition MCParticle.h:32
float getEnergy() const
Return particle energy in GeV.
Definition MCParticle.h:136
int getIndex() const
Get 1-based index of the particle in the corresponding MCParticle list.
Definition MCParticle.h:219
ROOT::Math::XYZVector getProductionVertex() const
Return production vertex position.
Definition MCParticle.h:178
float getCharge() const
Return the particle charge defined in TDatabasePDG.
Definition MCParticle.cc:36
int getPDG() const
Return PDG code of particle.
Definition MCParticle.h:101
ROOT::Math::XYZVector getMomentum() const
Return momentum.
Definition MCParticle.h:187
bool m_applySemiplane
Whether to require that the hit is in the expected semiplane.
bool m_applyWedge
Whether to require that the hit belong to the barrel part of the SVD.
void addInefficiencyPlots(TList *graphList=nullptr)
Function to create inefficiency plots and add them to list.
ROOT::Math::XYZVector m_magField
magnetic field needed set particle info
StoreArray< genfit::TrackCand > m_GenfitMCTrackCands
MC Genfit TrackCands StoreArray.
int nWedge
Counter for hits on wedged sensors.
StoreArray< SVDCluster > m_SVDClusters
SVDClusters StoreArray.
bool isInAnnulus(double hitDistance, double R, double dR)
Function to check if hitDistance is within a given annulus.
void initialize() override
Initializes the Module.
std::string m_mcParticlesName
MCParticle list name.
int nBarrel
Counter for hits on barrel sensors.
TH1 * duplicateHistogram(const char *newname, const char *newtitle, TH1 *h, TList *histoList=nullptr)
Function to clone a histogram.
float compute_thetaMS(MCParticleInfo &mcParticleInfo, VXDTrueHit *aTrueHit)
Calculate thetaMS.
void terminate() override
Termination action.
bool m_applyLap
Whether to require that the hit belongs to the first lap in the transverse plane.
MCTrackCandClassifierModule()
Constructor of the module.
StoreArray< genfit::TrackCand > m_GenfitIdealMCTrackCands
Ideal Genfit TrackCands StoreArray.
bool isFirstLap(double FirstHitTime, double HitTime, double LapTime)
Function to check if a hitTime is within a given lapTime, under consideration of m_fraction and with ...
void addEfficiencyPlots(TList *graphList=nullptr)
Function to create efficiency plots and add them to list.
StoreArray< PXDCluster > m_PXDClusters
PXDClusters StoreArray.
double semiPlane(ROOT::Math::XYZVector vertex, ROOT::Math::XYZVector center, ROOT::Math::XYZVector hit)
Function to get semiplane.
double theDistance(ROOT::Math::XYZVector center, ROOT::Math::XYZVector hit)
Get distance between two points.
TH3F * createHistogram3D(const char *name, const char *title, Int_t nbinsX, Double_t minX, Double_t maxX, const char *titleX, Int_t nbinsY, Double_t minY, Double_t maxY, const char *titleY, Int_t nbinsZ, Double_t minZ, Double_t maxZ, const char *titleZ, TList *histoList=nullptr)
Create a 3D ROOT Histogram.
TH1F * createHistogramsRatio(const char *name, const char *title, TH1 *hNum, TH1 *hDen, bool isEffPlot, int axisRef)
Function to create a ratio histogram from two histograms.
TFile * m_rootFilePtr
Pointer to root file used for storing histograms.
StoreArray< MCParticle > m_MCParticles
MCParticles StoreArray.
bool m_applyAnnulus
Whether to require that the hit is in the expected annulus.
std::string m_mcTrackCandsColName
TrackCand list name.
bool isInSemiPlane(double semiPlane, double omega)
Function to check if a omega value is in a given semiPlane.
float compute_dR(double thetaMS, double omega)
Calculate dR.
bool m_removeBadHits
Whether to remove the clusters that do not satisfy the criteria from the idealMCTrackCands.
int m_minHit
Minimum number of 1D Clusters to classify the MCTrackCand as ideal.
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
The PXD Cluster class This class stores all information about reconstructed PXD clusters The position...
Definition PXDCluster.h:30
Class PXDTrueHit - Records of tracks that either enter or leave the sensitive volume.
Definition PXDTrueHit.h:31
Class for type safe access to objects that are referred to in relations.
size_t size() const
Get number of relations.
RelationVector< T > getRelationsWith(const std::string &name="", const std::string &namedRelation="") const
Get the relations between this object and another store array.
The SVD Cluster class This class stores all information about reconstructed SVD clusters.
Definition SVDCluster.h:29
bool isUCluster() const
Get the direction of strips.
Definition SVDCluster.h:110
Class SVDTrueHit - Records of tracks that either enter or leave the sensitive volume.
Definition SVDTrueHit.h:33
static const double T
[tesla]
Definition Unit.h:120
Class VXDTrueHit - Records of tracks that either enter or leave the sensitive volume.
Definition VXDTrueHit.h:36
float getV() const
Return local v coordinate of hit.
Definition VXDTrueHit.h:76
float getGlobalTime() const
Return the time when the track reached its midpoint.
Definition VXDTrueHit.h:94
float getEntryU() const
Return local u coordinate of hit when entering silicon.
Definition VXDTrueHit.h:80
float getExitW() const
Return local w coordinate of hit at the endpoint of the track.
Definition VXDTrueHit.h:90
float getEntryW() const
Return local w coordinate of the start point of the track.
Definition VXDTrueHit.h:84
VxdID getSensorID() const
Return the Sensor ID.
Definition VXDTrueHit.h:70
float getU() const
Return local u coordinate of hit.
Definition VXDTrueHit.h:74
float getExitU() const
Return local u coordinate of hit at the endpoint of the track.
Definition VXDTrueHit.h:86
float getExitV() const
Return local v coordinate of hit at the endpoint of the track.
Definition VXDTrueHit.h:88
float getEntryV() const
Return local v coordinate of the start point of the track.
Definition VXDTrueHit.h:82
Class to facilitate easy access to sensor information of the VXD like coordinate transformations or p...
Definition GeoCache.h:39
const SensorInfoBase & getSensorInfo(Belle2::VxdID id) const
Return a reference to the SensorInfo of a given SensorID.
Definition GeoCache.cc:67
static GeoCache & getInstance()
Return a reference to the singleton instance.
Definition GeoCache.cc:214
Base class to provide Sensor Information for PXD and SVD.
ROOT::Math::XYZVector pointToGlobal(const ROOT::Math::XYZVector &local, bool reco=false) const
Convert a point from local to global coordinates.
double getForwardWidth() const
Convenience Wrapper to return width at forward side.
double getBackwardWidth() const
Convenience Wrapper to return width at backward side.
Class to uniquely identify a any structure of the PXD and SVD.
Definition VxdID.h:33
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
static void getField(const double *pos, double *field)
return the magnetic field at a given position.
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Abstract base class for different kinds of events.