Belle II Software  release-05-01-25
SVDROIFinderAnalysisModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2011 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Giulia Casarosa, Eugenio Paoloni *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <tracking/modules/svdROIFinder/SVDROIFinderAnalysisModule.h>
12 #include <framework/datastore/StoreArray.h>
13 #include <framework/datastore/RelationIndex.h>
14 #include <framework/logging/Logger.h>
15 
16 #include <svd/dataobjects/SVDShaperDigit.h>
17 #include <svd/dataobjects/SVDTrueHit.h>
18 #include <iostream>
19 #include <TVector3.h>
20 
21 #include <vxd/geometry/GeoCache.h>
22 
23 
24 using namespace std;
25 using namespace Belle2;
26 
27 //-----------------------------------------------------------------
28 // Register the Module
29 //-----------------------------------------------------------------
30 REG_MODULE(SVDROIFinderAnalysis)
31 
32 //-----------------------------------------------------------------
33 // Implementation
34 //-----------------------------------------------------------------
35 
36 SVDROIFinderAnalysisModule::SVDROIFinderAnalysisModule() : Module()
37  , m_recoTrackListName()
38  , m_SVDInterceptListName()
39  , m_ROIListName()
40  , m_rootFilePtr(nullptr)
41  , m_rootFileName("")
42  , m_writeToRoot(false)
43  , m_rootEvent(-1)
44  //efficiency graphs
45  , m_gEff2(nullptr)
46  , m_gEff(nullptr)
47  , m_h1DigitsPerParticle(nullptr)
48  , m_h1RecoTracksPerParticle(nullptr)
49  , m_h1digiIn(nullptr)
50  , m_h1digiOut2(nullptr)
51  , m_h1digiOut3(nullptr)
52  , m_h1digiOut4(nullptr)
53  , m_h1digiOut5(nullptr)
54  //tracks with no digit in ROI
55  , m_h1TrackOneDigiIn(nullptr)
56  , m_h1nnotINtrack2(nullptr)
57  , m_h1nnotINtrack3(nullptr)
58  , m_h1nnotINtrack4(nullptr)
59  , m_h1nnotINtrack5(nullptr)
60  //all tracks
61  , m_h1Track(nullptr)
62  , m_h1Track_pt(nullptr)
63  , m_h1Track_phi(nullptr)
64  , m_h1Track_lambda(nullptr)
65  , m_h1Track_cosTheta(nullptr)
66  , m_h1Track_pVal(nullptr)
67  , m_h1Track_nSVDhits(nullptr)
68  , m_h1Track_nCDChits(nullptr)
69  //tracks with at least one digit in ROI
70  , m_h1INtrack1(nullptr)
71  , m_h1INtrack1_pt(nullptr)
72  , m_h1INtrack1_phi(nullptr)
73  , m_h1INtrack1_lambda(nullptr)
74  , m_h1INtrack1_cosTheta(nullptr)
75  , m_h1INtrack1_pVal(nullptr)
76  , m_h1INtrack1_nSVDhits(nullptr)
77  , m_h1INtrack1_nCDChits(nullptr)
78  //tracks with no intercept
79  , m_h1notINtrack5(nullptr)
80  , m_h1notINtrack5_pt(nullptr)
81  , m_h1notINtrack5_phi(nullptr)
82  , m_h1notINtrack5_lambda(nullptr)
83  , m_h1notINtrack5_cosTheta(nullptr)
84  , m_h1notINtrack5_pVal(nullptr)
85  , m_h1notINtrack5_nSVDhits(nullptr)
86  , m_h1notINtrack5_nCDChits(nullptr)
87  //digits inside ROI
88  , m_h1PullU(nullptr)
89  , m_h1PullV(nullptr)
90  , m_h2sigmaUphi(nullptr)
91  , m_h2sigmaVphi(nullptr)
92  , m_h1ResidU(nullptr)
93  , m_h1ResidV(nullptr)
94  , m_h1SigmaU(nullptr)
95  , m_h1SigmaV(nullptr)
96  , m_h1GlobalTime(nullptr)
97  //digits outside2 ROI
98  , m_h2sigmaUphi_out2(nullptr)
99  , m_h2sigmaVphi_out2(nullptr)
100  , m_h1ResidU_out2(nullptr)
101  , m_h1ResidV_out2(nullptr)
102  , m_h1SigmaU_out2(nullptr)
103  , m_h1SigmaV_out2(nullptr)
104  , m_h1GlobalTime_out2(nullptr)
105  //digits outside3 ROI
106  , m_h2sigmaUphi_out3(nullptr)
107  , m_h2sigmaVphi_out3(nullptr)
108  , m_h1ResidU_out3(nullptr)
109  , m_h1ResidV_out3(nullptr)
110  , m_h1SigmaU_out3(nullptr)
111  , m_h1SigmaV_out3(nullptr)
112  , m_h1GlobalTime_out3(nullptr)
113  //digits outside4 ROI
114  , m_h2sigmaUphi_out4(nullptr)
115  , m_h2sigmaVphi_out4(nullptr)
116  , m_h1SigmaU_out4(nullptr)
117  , m_h1SigmaV_out4(nullptr)
118  , m_h1GlobalTime_out4(nullptr)
119  //digits outside5 ROI
120  , m_h1GlobalTime_out5(nullptr)
121 
122  //ROI stuff
123  , m_h2ROIbottomLeft(nullptr)
124  , m_h2ROItopRight(nullptr)
125  , m_h2ROIuMinMax(nullptr)
126  , m_h2ROIvMinMax(nullptr)
127  , m_h1totROIs(nullptr)
128  , m_h1okROIs(nullptr)
129  , m_h1totUstrips(nullptr)
130  , m_h1totVstrips(nullptr)
131 
132  , m_h1effPerTrack(nullptr)
133 
134  //variables
135  , m_globalTime(0.)
136  , m_coorU(0.)
137  , m_coorV(0.)
138  , m_sigmaU(0.)
139  , m_sigmaV(0.)
140  , m_vxdID(-1)
141 
142  , m_coormc(0.)
143  , m_idmc(-1)
144  , m_vxdIDmc(-1)
145  , m_pTmc(0.)
146  , m_momXmc(0.)
147  , m_momYmc(0.)
148  , m_momZmc(0.)
149  , m_thetamc(0.)
150  , m_costhetamc(0.)
151  , m_phimc(0.)
152  , m_lambdamc(0.)
153 
154  , Ntrack(0)
155  , NtrackHit(0)
156  , n_notINtrack2(0)
157  , n_notINtrack3(0)
158  , n_notINtrack4(0)
159  , n_notINtrack5(0)
160 
161  , n_rois(0)
162  , n_OKrois(0)
163  , m_nGoodROIs(0)
164  , n_intercepts(0)
165  , n_tracks(0)
166  , n_tracksWithDigits(0)
167  , n_tracksWithDigitsInROI(0)
168 
169 
170  , n_svdDigit(0)
171  , n_svdDigitInROI(0)
172  , n_notINdigit2(0)
173  , n_notINdigit3(0)
174  , n_notINdigit4(0)
175  , n_notINdigit5(0)
176 
177  //vectors
178  , nsvdDigit{0}
179  , nsvdDigitInROI{0}
180  , nnotINdigit2{0}
181  , nnotINdigit3{0}
182  , nnotINdigit4{0}
183  , nnotINdigit5{0}
184  , TrackOneDigiIn{0}
185  , nnotINtrack2{0}
186  , nnotINtrack3{0}
187  , nnotINtrack4{0}
188  , nnotINtrack5{0}
189 
190 {
191  //Set module properties
192  setDescription("This module performs the analysis of the SVDROIFinder module output");
193 
194  addParam("isSimulation", m_isSimulation,
195  "set true if you want to evaluate efficiency on simulation", bool(true));
196  addParam("writeToRoot", m_writeToRoot,
197  "set true if you want to save the informations in a root file named by parameter 'rootFileName'", bool(true));
198 
199  addParam("rootFileName", m_rootFileName,
200  "fileName used for . Will be ignored if parameter 'writeToRoot' is false (standard)",
201  string("svdDataRedAnalysis"));
202 
203  addParam("recoTrackListName", m_recoTrackListName,
204  "name of the input collection of RecoTracks", std::string(""));
205 
206  addParam("shapers", m_shapersName,
207  "name of the input collection of SVDShaperDigits", std::string(""));
208 
209  addParam("SVDInterceptListName", m_SVDInterceptListName,
210  "name of the list of interceptions", std::string(""));
211 
212  addParam("ROIListName", m_ROIListName,
213  "name of the list of ROIs", std::string(""));
214 
215  m_rootEvent = 0;
216 }
217 
218 SVDROIFinderAnalysisModule::~SVDROIFinderAnalysisModule()
219 {
220 }
221 
222 
223 void SVDROIFinderAnalysisModule::initialize()
224 {
225 
226  m_shapers.isRequired(m_shapersName);
227  m_trackList.isRequired(m_recoTrackListName);
228  m_ROIs.isRequired(m_ROIListName);
229  m_SVDIntercepts.isRequired(m_SVDInterceptListName);
230  m_mcParticles.isRequired();
231 
232  n_rois = 0;
233  n_OKrois = 0; //simulation
234  m_nGoodROIs = 0; //data
235  n_intercepts = 0;
236  n_tracks = 0;
237  n_svdDigit = 0;
238  n_svdDigitInROI = 0;
239 
240  n_tracksWithDigits = 0;
241  n_tracksWithDigitsInROI = 0;
242  NtrackHit = 0;
243  n_notINtrack2 = 0;
244  n_notINtrack3 = 0;
245  n_notINtrack4 = 0;
246  n_notINtrack5 = 0;
247 
248  for (int i = 0; i < 6; i++) {
249  nsvdDigit[i] = 0;
250  nsvdDigitInROI[i] = 0;
251  }
252 
253  if (m_writeToRoot == true) {
254  m_rootFileName += ".root";
255  m_rootFilePtr = new TFile(m_rootFileName.c_str(), "RECREATE");
256  } else
257  m_rootFilePtr = nullptr;
258 
259 
260  m_h1GlobalTime = new TH1F("hGlobalTime", "global time for SVDShaperDigits contained in ROI", 200, -100, 100);
261  m_h1PullU = new TH1F("hPullU", "U pulls for SVDShaperDigits contained in ROI", 100, -6, 6);
262  m_h1PullV = new TH1F("hPullV", "V pulls for SVDShaperDigits contained in ROI", 100, -6, 6);
263  m_h2sigmaUphi = new TH2F("hsigmaUvsPhi", "sigmaU vs phi digits in ROI", 100, -180, 180, 100, 0, 0.35);
264  m_h2sigmaVphi = new TH2F("hsigmaVvsPhi", "sigmaU vs phi digits in ROI", 100, -180, 180, 100, 0, 0.4);
265  m_h1ResidU = new TH1F("hResidU", "U resid for SVDShaperDigits contained in ROI", 100, -0.5, 0.5);
266  m_h1ResidV = new TH1F("hResidV", "V resid for SVDShaperDigits contained in ROI", 100, -0.5, 0.5);
267  m_h1SigmaU = new TH1F("hSigmaU", "sigmaU for SVDShaperDigits contained in ROI", 100, 0, 0.35);
268  m_h1SigmaV = new TH1F("hSigmaV", "sigmaV for SVDShaperDigits contained in ROI", 100, 0, 0.35);
269 
270 
271  m_h1GlobalTime_out2 = new TH1F("hGlobalTime_out2", "global time for SVDShaperDigits not contained in ROI", 200, -100, 100);
272  m_h2sigmaUphi_out2 = new TH2F("hsigmaUvsPhi_out2", "sigmaU vs phi digits not contained in ROI", 100, -180, 180, 100, 0, 0.35);
273  m_h2sigmaVphi_out2 = new TH2F("hsigmaVvsPhi_out2", "sigmaU vs phi digits not contained in ROI", 100, -180, 180, 100, 0, 0.4);
274  m_h1ResidU_out2 = new TH1F("hResidU_out2", "U resid for SVDShaperDigits not contained in ROI", 100, -2.5, 2.5);
275  m_h1ResidV_out2 = new TH1F("hResidV_out2", "V resid for SVDShaperDigits not contained in ROI", 100, -2.5, 2.5);
276  m_h1SigmaU_out2 = new TH1F("hSigmaU_out2", "sigmaU for SVDShaperDigits not contained in ROI", 100, 0, 0.35);
277  m_h1SigmaV_out2 = new TH1F("hSigmaV_out2", "sigmaV for SVDShaperDigits not contained in ROI", 100, 0, 0.35);
278 
279  m_h1GlobalTime_out3 = new TH1F("hGlobalTime_out3", "global time for SVDShaperDigits not contained in ROI", 200, -100, 100);
280  m_h2sigmaUphi_out3 = new TH2F("hsigmaUvsPhi_out3", "sigmaU vs phi digits not contained in ROI", 100, -180, 180, 100, 0, 0.35);
281  m_h2sigmaVphi_out3 = new TH2F("hsigmaVvsPhi_out3", "sigmaU vs phi digits not contained in ROI", 100, -180, 180, 100, 0, 0.4);
282  m_h1ResidU_out3 = new TH1F("hResidU_out3", "U resid for SVDShaperDigits not contained in ROI", 100, -2.5, 2.5);
283  m_h1ResidV_out3 = new TH1F("hResidV_out3", "V resid for SVDShaperDigits not contained in ROI", 100, -2.5, 2.5);
284  m_h1SigmaU_out3 = new TH1F("hSigmaU_out3", "sigmaU for SVDShaperDigits not contained in ROI", 100, 0, 0.35);
285  m_h1SigmaV_out3 = new TH1F("hSigmaV_out3", "sigmaV for SVDShaperDigits not contained in ROI", 100, 0, 0.35);
286 
287 
288  m_h1GlobalTime_out4 = new TH1F("hGlobalTime_out4", "global time for SVDShaperDigits not contained in ROI", 200, -100, 100);
289  m_h2sigmaUphi_out4 = new TH2F("hsigmaUvsPhi_out4", "sigmaU vs phi digits not contained in ROI", 100, -180, 180, 100, 0, 0.35);
290  m_h2sigmaVphi_out4 = new TH2F("hsigmaVvsPhi_out4", "sigmaU vs phi digits not contained in ROI", 100, -180, 180, 100, 0, 0.4);
291  m_h1SigmaU_out4 = new TH1F("hSigmaU_out4", "sigmaU for SVDShaperDigits not contained in ROI", 100, 0, 0.35);
292  m_h1SigmaV_out4 = new TH1F("hSigmaV_out4", "sigmaV for SVDShaperDigits not contained in ROI", 100, 0, 0.35);
293 
294 
295  m_h1GlobalTime_out5 = new TH1F("hGlobalTime_out5", "global time for SVDShaperDigits not contained in ROI", 200, -100, 100);
296 
297 
298 
299  m_h1totROIs = new TH1F("h1TotNROIs", "number of all ROIs", 110, 0, 110);
300  m_h1okROIs = new TH1F("h1OkNROIs", "number of all ROIs containing a SVDShaperDigit", 110, 0, 110);
301 
302  m_h1totUstrips = new TH1F("h1TotUstrips", "number of U strips in ROIs", 100, 0, 250000);
303  m_h1totVstrips = new TH1F("h1TotVstrips", "number of V strips in ROIs", 100, 0, 250000);
304 
305  m_h1effPerTrack = new TH1F("heffPerTrack", "fraction of digits in ROI per track", 100, -0.02, 1.02);
306 
307 
308 
309  m_h2ROIbottomLeft = new TH2F("h2ROIbottomLeft", "u,v ID of the bottom left pixel", 650, -200, 450, 1300, -300, 1000);
310  m_h2ROItopRight = new TH2F("h2ROItopRight", "u,v ID of the top right pixel", 650, -200, 450, 1300, -300, 1000);
311 
312  m_h2ROIuMinMax = new TH2F("h2ROIuMinMax", "u Min vs Max", 650, -200, 450, 650, -200, 450);
313  m_h2ROIvMinMax = new TH2F("h2ROIvMinMax", "v Min vs Max", 1300, -300, 1000, 1300, -300, 1000);
314 
315 
316  m_h1DigitsPerParticle = new TH1F("h1DigitsPerPart", "Number of SVDShaperDigits per Particle", 50, 0, 50);
317  m_h1RecoTracksPerParticle = new TH1F("h1RecoTracksPerPart", "Number of RecoTracks per Particle", 10, 0, 10);
318 
319 
320  //analysis
321  Double_t lowBin[6 + 1];
322  for (int i = 0; i < 6; i++)
323  lowBin[i] = pt[i] - ptErr[i];
324  lowBin[6] = pt[5] + ptErr[5];
325 
326  m_h1TrackOneDigiIn = new TH1F("hTracksDigiIn", "Tracks with at least one digit contained in a ROI", 6, lowBin);
327  m_h1nnotINtrack2 = new TH1F("h1outROITrack", "Tracks with ROI with correct VxdID but no digits inside ROI", 6, lowBin);
328  m_h1nnotINtrack3 = new TH1F("h1noROITrack", "Tracks with ROI with wrong VxdID but no digits inside ROI", 6, lowBin);
329  m_h1nnotINtrack4 = new TH1F("h1wrongVxdIDTrack", "Tracks with no ROI, Intercept with correct VxdID", 6, lowBin);
330  m_h1nnotINtrack5 = new TH1F("h1noInterTrack", "Tracks with no Intercept matching a VxdID of digits", 6, lowBin);
331 
332  m_h1notINtrack5 = new TH1F("hNoInterTrack", "track with no intercepts", 20, 0, 20);
333  m_h1notINtrack5_pt = new TH1F("hNoInterTrack_pT", "track with no intercepts", 100, 0, 6);
334  m_h1notINtrack5_phi = new TH1F("h1NoInterTrack_phi", "hNoInterTrack_phi", 100, -180, 180);
335  m_h1notINtrack5_lambda = new TH1F("h1NoInterTrack_lambda", "hNoInterTrack_lambda", 100, -180, 180);
336  m_h1notINtrack5_cosTheta = new TH1F("h1NoInterTrack_cosTheta", "hNoInterTrack_cosTheta", 100, -1, 1);
337  m_h1notINtrack5_pVal = new TH1F("hNoInterTrack_pVal", "track with no intercepts", 100, 0, 1);
338  m_h1notINtrack5_nSVDhits = new TH1F("hNoInterTrack_nSVDhits", "track with no intercepts", 50, 0, 50);
339  m_h1notINtrack5_nCDChits = new TH1F("hNoInterTrack_nCDChits", "track with no intercepts", 100, 0, 100);
340 
341  m_h1INtrack1 = new TH1F("hINTrack", "track with at least one digit inside ROI", 20, 0, 20);
342  m_h1INtrack1_pt = new TH1F("hINTrack_pT", "track with at least one digit inside ROI", 100, 0, 6);
343  m_h1INtrack1_phi = new TH1F("h1INTrack_phi", "hINTrack_phi", 100, -180, 180);
344  m_h1INtrack1_lambda = new TH1F("h1INTrack_lambda", "hINTrack_lambda", 100, -180, 180);
345  m_h1INtrack1_cosTheta = new TH1F("h1INTrack_cosTheta", "hINTrack_cosTheta", 100, -1, 1);
346  m_h1INtrack1_pVal = new TH1F("h1INTrack_pVal", "track with no intercepts", 100, 0, 1);
347  m_h1INtrack1_nSVDhits = new TH1F("h1INTrack_nSVDhits", "track with no intercepts", 50, 0, 50);
348  m_h1INtrack1_nCDChits = new TH1F("h1INTrack_nCDChits", "track with no intercepts", 100, 0, 100);
349 
350  m_h1Track = new TH1F("hTrack", "all tracks", 20, 0, 20);
351  m_h1Track_pt = new TH1F("hTrack_pT", "all tracks with digits", 100, 0, 6);
352  m_h1Track_lambda = new TH1F("h1Track_lambda", "hTrack_lambda", 100, -180, 180);
353  m_h1Track_phi = new TH1F("h1Track_phi", "hTrack_phi", 100, -180, 180);
354  m_h1Track_cosTheta = new TH1F("h1Track_cosTheta", "hTrack_cos theta", 100, -1, 1);
355  m_h1Track_pVal = new TH1F("h1Track_pVal", "track with no intercepts", 100, 0, 1);
356  m_h1Track_nSVDhits = new TH1F("h1Track_nSVDhits", "track with no intercepts", 50, 0, 50);
357  m_h1Track_nCDChits = new TH1F("h1Track_nCDChits", "track with no intercepts", 100, 0, 100);
358 
359  m_h1digiIn = new TH1F("hdigiIn", "digits inside ROI", 6, lowBin);
360  m_h1digiOut2 = new TH1F("hdigiOut2", "ROI exists with with correct VxdID but no digits inside ROI", 6, lowBin);
361  m_h1digiOut3 = new TH1F("hdigiOut3", "ROI exists with with wrong VxdID", 6, lowBin);
362  m_h1digiOut4 = new TH1F("hdigiOut4", "ROI does not exist, but intercept has correct VxdID", 6, lowBin);
363  m_h1digiOut5 = new TH1F("hdigiOut5", "no ROI, no Intercpets with correct VXDid", 6, lowBin);
364 
365  // m_h2_VXDhitsPR_xy = new TH2F("hNoInteTrack_SVDhitsXY", "SVD Hits Missed by the VXDTF", 200, -15, 15, 200, -15, 15);
366 
367  // m_h2_VXDhitsPR_rz = new TH2F("hNoInteTrack_SVDhitsRZ", "SVD Hits Missed by the VXDTF, r_{T} z", 200, -30, 40, 200, 0, 15);
368 
369 
370 }
371 
372 void SVDROIFinderAnalysisModule::beginRun()
373 {
374  m_rootEvent = 0;
375 }
376 
377 
378 void SVDROIFinderAnalysisModule::event()
379 {
380 
381  typedef RelationIndex < RecoTrack, SVDIntercept>::range_from SVDInterceptsFromRecoTracks;
382  typedef RelationIndex < RecoTrack, SVDIntercept>::iterator_from SVDInterceptIteratorType;
383  typedef RelationIndex < SVDShaperDigit, SVDTrueHit>::range_from SVDTrueHitFromSVDShaperDigit;
384  typedef RelationIndex < SVDShaperDigit, SVDTrueHit>::iterator_from SVDTrueHitIteratorType;
385  RelationIndex < SVDShaperDigit, SVDTrueHit >
386  relDigitTrueHit(DataStore::relationName(DataStore::arrayName<SVDShaperDigit>(""),
387  DataStore::arrayName<SVDTrueHit>("")));
388  RelationIndex < RecoTrack, SVDIntercept >
389  recoTrackToSVDIntercept(DataStore::relationName(m_recoTrackListName, m_SVDInterceptListName));
390 
391  double tmpGlobalTime;
392  int tmpNGlobalTime;
393 
394  NtrackHit = 0;
395  Ntrack = 0;
396 
397  B2DEBUG(1, " ++++++++++++++ SVDROIFinderAnalysisModule");
398 
399  int nROIs = 0;
400 
401 
402 
403  //ROIs general
404  for (int i = 0; i < (int)m_ROIs.getEntries(); i++) { //loop on ROIlist
405 
406  m_h2ROIbottomLeft->Fill(m_ROIs[i]->getMinUid(), m_ROIs[i]->getMinVid());
407  m_h2ROItopRight->Fill(m_ROIs[i]->getMaxUid(), m_ROIs[i]->getMaxVid());
408  m_h2ROIuMinMax->Fill(m_ROIs[i]->getMinUid(), m_ROIs[i]->getMaxUid());
409  m_h2ROIvMinMax->Fill(m_ROIs[i]->getMinVid(), m_ROIs[i]->getMaxVid());
410 
411  for (int s = 0; s < m_shapers.getEntries(); s++) {
412  if (m_ROIs[i]->Contains(*(m_shapers[s]))) {
413  m_nGoodROIs++;
414  break;
415  }
416 
417  }
418 
419 
420  bool isOK = false;
421 
422  for (int j = 0; j < (int)m_mcParticles.getEntries(); j++) {
423  MCParticle* aMcParticle = m_mcParticles[j];
424 
425  // continue only if MCParticle has a related SVDShaperDigit and RecoTrack
426  RelationVector<SVDShaperDigit> svdDigits_MCParticle = aMcParticle->getRelationsFrom<SVDShaperDigit>();
427 
428  if (!isOK)
429  //loop on SVDShaperDigits
430  for (unsigned int iSVDShaperDigit = 0; iSVDShaperDigit < svdDigits_MCParticle.size(); iSVDShaperDigit++)
431  if (m_ROIs[i]->Contains(*(svdDigits_MCParticle[iSVDShaperDigit]))) {
432  nROIs++;
433  isOK = true;
434  break;
435  }
436  }
437  }
438 
439  m_h1totROIs->Fill(m_ROIs.getEntries());
440  n_rois += m_ROIs.getEntries();
441 
442  //RecoTrack general
443  n_tracks += m_trackList.getEntries();
444 
445  //SVDIntercepts general
446  n_intercepts += m_SVDIntercepts.getEntries();
447 
448  Int_t n_NoInterceptTracks = 0;
449 
450  //loop on MCParticles
451  for (int j = 0; j < (int)m_mcParticles.getEntries(); j++) {
452 
453  MCParticle* aMcParticle = m_mcParticles[j];
454 
455  // continue only if MCParticle has a related SVDShaperDigit and RecoTrack
456  RelationVector<SVDShaperDigit> svdDigits_MCParticle = aMcParticle->getRelationsFrom<SVDShaperDigit>();
457  RelationVector<RecoTrack> recoTracks_MCParticle = aMcParticle->getRelationsWith<RecoTrack>();
458 
459  m_h1DigitsPerParticle->Fill(svdDigits_MCParticle.size());
460  if (svdDigits_MCParticle.size() == 0)
461  continue;
462 
463  // m_h1RecoTracksPerParticle->Fill(recoTracks_MCParticle.size());
464  // if (recoTracks_MCParticle.size() == 0)
465  // continue;
466 
467  Ntrack++;
468 
469  B2DEBUG(1, "Number of RecoTracks = " << recoTracks_MCParticle.size() << " and SVDShaperDigits = " << svdDigits_MCParticle.size() <<
470  " related to this MCParticle");
471 
472  //retrieve general informations of MCParticle
473  m_momXmc = (aMcParticle->getMomentum()).X();
474  m_momYmc = (aMcParticle->getMomentum()).Y();
475  m_momZmc = (aMcParticle->getMomentum()).Z();
476  m_phimc = (aMcParticle->getMomentum()).Phi() * 180 / 3.1415;
477  m_thetamc = (aMcParticle->getMomentum()).Theta() * 180 / 3.1415;
478  m_costhetamc = (aMcParticle->getMomentum()).CosTheta();
479  m_lambdamc = 90 - m_thetamc;
480  m_pTmc = (aMcParticle->getMomentum()).Perp();
481 
482 
483  bool part_outsideROI = false;
484  bool part_noROI = false;
485  bool part_wrongVxdID = false;
486  bool part_noInter = false;
487  bool hasOneDigitInROI = false;
488 
489  Int_t nDigitsInRoiPerTrack = 0;
490  Int_t nDigitsPerTrack = 0;
491 
492  //loop on SVDShaperDigits
493  for (unsigned int iSVDShaperDigit = 0; iSVDShaperDigit < svdDigits_MCParticle.size(); iSVDShaperDigit++) {
494 
495  bool isU = svdDigits_MCParticle[iSVDShaperDigit]->isUStrip();
496 
497  bool hasIntercept = false;
498  bool hasROI = false;
499  bool interceptRightVxdID = false;
500  bool MissingHit = true;
501 
502  n_svdDigit ++ ;
503  nDigitsPerTrack++;
504 
505  SVDTrueHitFromSVDShaperDigit SVDTrueHits = relDigitTrueHit.getElementsFrom(*svdDigits_MCParticle[iSVDShaperDigit]);
506  SVDTrueHitIteratorType theSVDTrueHitIterator = SVDTrueHits.begin();
507  SVDTrueHitIteratorType theSVDTrueHitIteratorEnd = SVDTrueHits.end();
508  tmpGlobalTime = 0;
509  tmpNGlobalTime = 0;
510 
511  for (; theSVDTrueHitIterator != theSVDTrueHitIteratorEnd; theSVDTrueHitIterator++) {
512  tmpGlobalTime = tmpGlobalTime + theSVDTrueHitIterator->to->getGlobalTime();
513  tmpNGlobalTime++;
514  }
515 
516  m_globalTime = tmpGlobalTime / tmpNGlobalTime;
517  m_idmc = svdDigits_MCParticle[iSVDShaperDigit]->getCellID();
518  m_vxdIDmc = svdDigits_MCParticle[iSVDShaperDigit]->getSensorID();
519 
520  VXD::GeoCache& aGeometry = VXD::GeoCache::getInstance();
521  const VXD::SensorInfoBase& aSensorInfo = aGeometry.getSensorInfo(m_vxdIDmc);
522 
523  if (isU)
524  m_coormc = aSensorInfo.getUCellPosition(m_idmc);
525  else
526  m_coormc = aSensorInfo.getVCellPosition(m_idmc);
527 
528  if (m_pTmc > 1) nsvdDigit[5]++;
529  if (m_pTmc <= 1 && m_pTmc > 0.5) nsvdDigit[4]++;
530  if (m_pTmc <= 0.5 && m_pTmc > 0.3) nsvdDigit[3]++;
531  if (m_pTmc <= 0.3 && m_pTmc > 0.2) nsvdDigit[2]++;
532  if (m_pTmc <= 0.2 && m_pTmc > 0.1) nsvdDigit[1]++;
533  if (m_pTmc <= 0.1) nsvdDigit[0]++;
534 
535 
536  // for (int i = 0; i < (int)recoTracks_MCParticle.size(); i++) { //loop on related RecoTracks
537  for (int i = 0; i < (int)m_trackList.getEntries(); i++) { //loop on all RecoTracks
538 
539  // SVDInterceptsFromRecoTracks SVDIntercepts = recoTrackToSVDIntercept.getElementsFrom(recoTracks_MCParticle[i]);
540  SVDInterceptsFromRecoTracks SVDIntercepts = recoTrackToSVDIntercept.getElementsFrom(m_trackList[i]);
541 
542  SVDInterceptIteratorType theSVDInterceptIterator = SVDIntercepts.begin();
543  SVDInterceptIteratorType theSVDInterceptIteratorEnd = SVDIntercepts.end();
544 
545 
546  for (; theSVDInterceptIterator != theSVDInterceptIteratorEnd; theSVDInterceptIterator++) {
547 
548  const SVDIntercept* theIntercept = theSVDInterceptIterator->to;
549 
550  if (theIntercept) {
551 
552  hasIntercept = true;
553 
554  m_coorU = theIntercept->getCoorU();
555  m_coorV = theIntercept->getCoorV();
556  m_sigmaU = theIntercept->getSigmaU();
557  m_sigmaV = theIntercept->getSigmaV();
558  m_vxdID = theIntercept->getSensorID();
559 
560  if (m_vxdID == m_vxdIDmc)
561  interceptRightVxdID = true;
562  else
563  continue;
564 
565  const ROIid* theROIid = theIntercept->getRelatedTo<ROIid>(m_ROIListName);
566 
567  if (theROIid) {
568 
569  hasROI = true;
570 
571  if (theROIid->Contains(*(svdDigits_MCParticle[iSVDShaperDigit]))) { //CASO1
572 
573  if (MissingHit) {
574  nDigitsInRoiPerTrack++;
575 
576  m_h1GlobalTime->Fill(m_globalTime);
577 
578  if (isU) {
579  m_h1PullU->Fill((m_coorU - m_coormc) / m_sigmaU);
580  m_h1ResidU->Fill(m_coorU - m_coormc);
581  m_h2sigmaUphi->Fill(m_phimc, m_sigmaU);
582  m_h1SigmaU->Fill(m_sigmaU);
583  } else {
584  m_h1PullV->Fill((m_coorV - m_coormc) / m_sigmaV);
585  m_h1ResidV->Fill(m_coorV - m_coormc);
586  m_h2sigmaVphi->Fill(m_phimc, m_sigmaV);
587  m_h1SigmaV->Fill(m_sigmaV);
588  }
589 
590  hasOneDigitInROI = true;
591  n_svdDigitInROI++;
592 
593  if (m_pTmc > 1) nsvdDigitInROI[5]++;
594  if (m_pTmc <= 1 && m_pTmc > 0.5) nsvdDigitInROI[4]++;
595  if (m_pTmc <= 0.5 && m_pTmc > 0.3) nsvdDigitInROI[3]++;
596  if (m_pTmc <= 0.3 && m_pTmc > 0.2) nsvdDigitInROI[2]++;
597  if (m_pTmc <= 0.2 && m_pTmc > 0.1) nsvdDigitInROI[1]++;
598  if (m_pTmc <= 0.1) nsvdDigitInROI[0]++;
599 
600  MissingHit = false;
601  }
602 
603 
604  break; // To avoid double counting (intercepts)
605  } //if theROIid contains
606  } //if (theROIid)
607  } //if (theintercept)
608  } //(end loop on intercept list)
609 
610  if (!MissingHit)
611  break;// To avoid double counting (recoTracks)
612 
613  } //(end loop on recoTracks)
614 
615 
616  if (MissingHit) {
617 
618  if (hasROI && hasIntercept && interceptRightVxdID) {
619  part_outsideROI = true;
620 
621  n_notINdigit2 ++;
622 
623  m_h1GlobalTime_out2->Fill(m_globalTime);
624  if (isU) {
625  m_h1ResidU_out2->Fill(m_coorU - m_coormc);
626  m_h2sigmaUphi_out2->Fill(m_phimc, m_sigmaU);
627  m_h1SigmaU_out2->Fill(m_sigmaU);
628  } else {
629  m_h1ResidV_out2->Fill(m_coorV - m_coormc);
630  m_h2sigmaVphi_out2->Fill(m_phimc, m_sigmaV);
631  m_h1SigmaV_out2->Fill(m_sigmaV);
632  }
633 
634  if (m_pTmc > 1) nnotINdigit2[5]++;
635  if (m_pTmc <= 1 && m_pTmc > 0.5) nnotINdigit2[4]++;
636  if (m_pTmc <= 0.5 && m_pTmc > 0.3) nnotINdigit2[3]++;
637  if (m_pTmc <= 0.3 && m_pTmc > 0.2) nnotINdigit2[2]++;
638  if (m_pTmc <= 0.2 && m_pTmc > 0.1) nnotINdigit2[1]++;
639  if (m_pTmc <= 0.1) nnotINdigit2[0]++;
640 
641  } else if (!hasROI && hasIntercept && interceptRightVxdID) {
642  part_noROI = true;
643 
644  n_notINdigit3 ++;
645 
646  m_h1GlobalTime_out3->Fill(m_globalTime);
647  if (isU) {
648  m_h1ResidU_out3->Fill(m_coorU - m_coormc);
649  m_h2sigmaUphi_out3->Fill(m_phimc, m_sigmaU);
650  m_h1SigmaU_out3->Fill(m_sigmaU);
651  } else {
652  m_h1ResidV_out3->Fill(m_coorV - m_coormc);
653  m_h2sigmaVphi_out3->Fill(m_phimc, m_sigmaV);
654  m_h1SigmaV_out3->Fill(m_sigmaV);
655  }
656 
657  if (m_pTmc > 1) nnotINdigit3[5]++;
658  if (m_pTmc <= 1 && m_pTmc > 0.5) nnotINdigit3[4]++;
659  if (m_pTmc <= 0.5 && m_pTmc > 0.3) nnotINdigit3[3]++;
660  if (m_pTmc <= 0.3 && m_pTmc > 0.2) nnotINdigit3[2]++;
661  if (m_pTmc <= 0.2 && m_pTmc > 0.1) nnotINdigit3[1]++;
662  if (m_pTmc <= 0.1) nnotINdigit3[0]++;
663 
664  } else if (hasIntercept && !interceptRightVxdID) {
665  part_wrongVxdID = true;
666 
667  n_notINdigit4 ++;
668 
669  m_h1GlobalTime_out4->Fill(m_globalTime);
670  if (isU) {
671  m_h2sigmaUphi_out4->Fill(m_phimc, m_sigmaU);
672  m_h1SigmaU_out4->Fill(m_sigmaU);
673  } else {
674  m_h2sigmaVphi_out4->Fill(m_phimc, m_sigmaV);
675  m_h1SigmaV_out4->Fill(m_sigmaV);
676  }
677 
678  if (m_pTmc > 1) nnotINdigit4[5]++;
679  if (m_pTmc <= 1 && m_pTmc > 0.5) nnotINdigit4[4]++;
680  if (m_pTmc <= 0.5 && m_pTmc > 0.3) nnotINdigit4[3]++;
681  if (m_pTmc <= 0.3 && m_pTmc > 0.2) nnotINdigit4[2]++;
682  if (m_pTmc <= 0.2 && m_pTmc > 0.1) nnotINdigit4[1]++;
683  if (m_pTmc <= 0.1) nnotINdigit4[0]++;
684 
685  } else if (!hasIntercept) {
686  part_noInter = true;
687 
688  n_notINdigit5 ++;
689 
690  m_h1GlobalTime_out5->Fill(m_globalTime);
691 
692  if (m_pTmc > 1) nnotINdigit5[5]++;
693  if (m_pTmc <= 1 && m_pTmc > 0.5) nnotINdigit5[4]++;
694  if (m_pTmc <= 0.5 && m_pTmc > 0.3) nnotINdigit5[3]++;
695  if (m_pTmc <= 0.3 && m_pTmc > 0.2) nnotINdigit5[2]++;
696  if (m_pTmc <= 0.2 && m_pTmc > 0.1) nnotINdigit5[1]++;
697  if (m_pTmc <= 0.1) nnotINdigit5[0]++;
698  }
699  }
700  } //end loop on digits
701 
702  m_h1effPerTrack->Fill((float) nDigitsInRoiPerTrack / nDigitsPerTrack);
703  m_h1Track_pt->Fill(m_pTmc);
704  m_h1Track_phi->Fill(m_phimc);
705  m_h1Track_lambda->Fill(m_lambdamc);
706  m_h1Track_cosTheta->Fill(m_costhetamc);
707 
708  if (hasOneDigitInROI) {
709  NtrackHit++;
710  if (m_pTmc > 1) TrackOneDigiIn[5]++;
711  if (m_pTmc <= 1 && m_pTmc > 0.5) TrackOneDigiIn[4]++;
712  if (m_pTmc <= 0.5 && m_pTmc > 0.3) TrackOneDigiIn[3]++;
713  if (m_pTmc <= 0.3 && m_pTmc > 0.2) TrackOneDigiIn[2]++;
714  if (m_pTmc <= 0.2 && m_pTmc > 0.1) TrackOneDigiIn[1]++;
715  if (m_pTmc <= 0.1) TrackOneDigiIn[0]++;
716 
717  m_h1INtrack1_pt->Fill(m_pTmc);
718  m_h1INtrack1_phi->Fill(m_phimc);
719  m_h1INtrack1_lambda->Fill(m_lambdamc);
720  m_h1INtrack1_cosTheta->Fill(m_costhetamc);
721  } else if (part_outsideROI) { //CASO2
722  n_notINtrack2++;
723  if (m_pTmc > 1) nnotINtrack2[5]++;
724  if (m_pTmc <= 1 && m_pTmc > 0.5) nnotINtrack2[4]++;
725  if (m_pTmc <= 0.5 && m_pTmc > 0.3) nnotINtrack2[3]++;
726  if (m_pTmc <= 0.3 && m_pTmc > 0.2) nnotINtrack2[2]++;
727  if (m_pTmc <= 0.2 && m_pTmc > 0.1) nnotINtrack2[1]++;
728  if (m_pTmc <= 0.1) nnotINtrack2[0]++;
729  } else if (part_noROI) { //CASO3
730  n_notINtrack3++;
731  if (m_pTmc > 1) nnotINtrack3[5]++;
732  if (m_pTmc <= 1 && m_pTmc > 0.5) nnotINtrack3[4]++;
733  if (m_pTmc <= 0.5 && m_pTmc > 0.3) nnotINtrack3[3]++;
734  if (m_pTmc <= 0.3 && m_pTmc > 0.2) nnotINtrack3[2]++;
735  if (m_pTmc <= 0.2 && m_pTmc > 0.1) nnotINtrack3[1]++;
736  if (m_pTmc <= 0.1) nnotINtrack3[0]++;
737  } else if (part_wrongVxdID) { //CASO4
738  n_notINtrack4++;
739  if (m_pTmc > 1) nnotINtrack4[5]++;
740  if (m_pTmc <= 1 && m_pTmc > 0.5) nnotINtrack4[4]++;
741  if (m_pTmc <= 0.5 && m_pTmc > 0.3) nnotINtrack4[3]++;
742  if (m_pTmc <= 0.3 && m_pTmc > 0.2) nnotINtrack4[2]++;
743  if (m_pTmc <= 0.2 && m_pTmc > 0.1) nnotINtrack4[1]++;
744  if (m_pTmc <= 0.1) nnotINtrack4[0]++;
745  } else if (part_noInter) { //CASO5
746  n_notINtrack5++;
747  if (m_pTmc > 1) nnotINtrack5[5]++;
748  if (m_pTmc <= 1 && m_pTmc > 0.5) nnotINtrack5[4]++;
749  if (m_pTmc <= 0.5 && m_pTmc > 0.3) nnotINtrack5[3]++;
750  if (m_pTmc <= 0.3 && m_pTmc > 0.2) nnotINtrack5[2]++;
751  if (m_pTmc <= 0.2 && m_pTmc > 0.1) nnotINtrack5[1]++;
752  if (m_pTmc <= 0.1) nnotINtrack5[0]++;
753 
754  n_NoInterceptTracks++;
755  m_h1notINtrack5_phi->Fill(m_phimc);
756  m_h1notINtrack5_lambda->Fill(m_lambdamc);
757  m_h1notINtrack5_cosTheta->Fill(m_costhetamc);
758  m_h1notINtrack5_pt->Fill(m_pTmc);
759  }
760 
761  } // close loop on MCParticlet
762 
763 
764  m_h1notINtrack5->Fill(n_NoInterceptTracks);
765  m_h1Track->Fill(Ntrack);
766  m_h1INtrack1->Fill(NtrackHit);
767 
768  n_tracksWithDigits += Ntrack;
769  n_tracksWithDigitsInROI += NtrackHit;
770 
771  m_rootEvent++;
772  if ((m_ROIs.getEntries() > 0) && m_isSimulation) {
773  B2RESULT(" o SVDROIFinder ANALYSIS: tot ROIs = " << m_ROIs.getEntries() << ", ok ROIs = " << nROIs);
774  B2RESULT(" o : NtrackHit/Ntrack = " << NtrackHit << "/ " << Ntrack << " = " <<
775  (double)NtrackHit / Ntrack);
776  for (int i = 0; i < m_ROIs.getEntries(); i++) {
777  VxdID sensor = m_ROIs[i]->getSensorID();
778  B2RESULT(i << " ROI " << sensor.getLadderNumber() << "." << sensor.getLayerNumber() << "." << sensor.getSensorNumber() <<
779  ": " << m_ROIs[i]->getMinUid() << ", " << m_ROIs[i]->getMinVid() << ", " << m_ROIs[i]->getMaxUid() << ", " <<
780  m_ROIs[i]->getMaxVid());
781  }
782  }
783  if (nROIs > m_ROIs.getEntries()) B2RESULT(" HOUSTON WE HAVE A PROBLEM!");
784 
785  m_h1okROIs->Fill(nROIs);
786  n_OKrois += nROIs;
787 }
788 
789 
790 void SVDROIFinderAnalysisModule::endRun()
791 {
792 }
793 
794 
795 void SVDROIFinderAnalysisModule::terminate()
796 {
797 
798  Double_t epsilon[6];
799  Double_t epsilonErr[6];
800  double epsilonTot = (double)n_svdDigitInROI / (double) n_svdDigit;
801  Double_t epsilon2[6];
802  Double_t epsilon2Err[6];
803  double epsilon2Tot = (double)n_tracksWithDigitsInROI / (double) n_tracksWithDigits;
804 
805  for (int i = 0; i < 6; i++) {
806  m_h1digiOut2->SetBinContent(i + 1, nnotINdigit2[i]);
807  m_h1digiOut3->SetBinContent(i + 1, nnotINdigit3[i]);
808  m_h1digiOut4->SetBinContent(i + 1, nnotINdigit4[i]);
809  m_h1digiOut5->SetBinContent(i + 1, nnotINdigit5[i]);
810  m_h1digiIn->SetBinContent(i + 1, nsvdDigitInROI[i]);
811  }
812 
813  for (int i = 0; i < 6; i++) {
814  m_h1nnotINtrack2->SetBinContent(i + 1, nnotINtrack2[i]);
815  m_h1nnotINtrack3->SetBinContent(i + 1, nnotINtrack3[i]);
816  m_h1nnotINtrack4->SetBinContent(i + 1, nnotINtrack4[i]);
817  m_h1nnotINtrack5->SetBinContent(i + 1, nnotINtrack5[i]);
818  m_h1TrackOneDigiIn->SetBinContent(i + 1, TrackOneDigiIn[i]);
819  }
820 
821  B2RESULT(" ROI Analysis Summary ");
822  B2RESULT("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");
823  B2RESULT("");
824  B2RESULT(" number of tracks = " << n_tracks);
825  B2RESULT(" number of Intercepts = " << n_intercepts);
826  B2RESULT(" number of ROIs = " << n_rois);
827  if (m_isSimulation)
828  B2RESULT(" number of GOOD ROIs = " << n_OKrois);
829  // B2RESULT("");
830  B2RESULT(" average number of ROIs = " << m_h1totROIs->GetMean());
831  if (m_isSimulation) {
832  B2RESULT(" average number of ROIs w digits = " << m_h1okROIs->GetMean());
833  B2RESULT("");
834  B2RESULT("");
835  B2RESULT("tracks w digits: " << n_tracksWithDigits);
836  B2RESULT("tracks w digits in ROI: " << n_tracksWithDigitsInROI);
837 
838  B2RESULT("efficiency PTD : " << epsilon2Tot << " +/- " << sqrt(epsilon2Tot * (1 - epsilon2Tot) / n_tracksWithDigits));
839 
840  Int_t totTrackOneDigiIn = 0; //not used for the moment, added to double check
841  Int_t totnnotINtrack2 = 0;
842  Int_t totnnotINtrack3 = 0;
843  Int_t totnnotINtrack4 = 0;
844  Int_t totnnotINtrack5 = 0;
845 
846  Int_t totTrack[6];
847 
848  for (int j = 0; j < m_h1TrackOneDigiIn->GetNbinsX(); j++) {
849  totTrackOneDigiIn = totTrackOneDigiIn + m_h1TrackOneDigiIn->GetBinContent(j + 1);
850  totnnotINtrack2 = totnnotINtrack2 + m_h1nnotINtrack2->GetBinContent(j + 1);
851  totnnotINtrack3 = totnnotINtrack3 + m_h1nnotINtrack3->GetBinContent(j + 1);
852  totnnotINtrack4 = totnnotINtrack4 + m_h1nnotINtrack4->GetBinContent(j + 1);
853  totnnotINtrack5 = totnnotINtrack5 + m_h1nnotINtrack5->GetBinContent(j + 1);
854 
855  totTrack[j] = m_h1nnotINtrack5->GetBinContent(j + 1) + m_h1nnotINtrack4->GetBinContent(j + 1) + m_h1nnotINtrack3->GetBinContent(
856  j + 1) + m_h1nnotINtrack2->GetBinContent(j + 1) + m_h1TrackOneDigiIn->GetBinContent(j + 1);
857  }
858 
859  B2RESULT(" out ROI = " << totnnotINtrack2);
860  B2RESULT(" no ROI = " << totnnotINtrack3);
861  B2RESULT(" wrongVxdID = " << totnnotINtrack4);
862  B2RESULT(" no Inter = " << totnnotINtrack5);
863  B2RESULT("");
864 
865  B2RESULT(" svdDigit : " << n_svdDigit);
866  B2RESULT(" svdDigitIn : " << n_svdDigitInROI);
867 
868  B2RESULT(" eff DGT: " << epsilonTot << " +/- " << sqrt(epsilonTot * (1 - epsilonTot) / n_svdDigit));
869  B2RESULT(" inefficiency (SVDShaperDigits): ");
870  B2RESULT(" out ROI: " << n_notINdigit2);
871  B2RESULT(" no ROI: " << n_notINdigit3);
872  B2RESULT(" wrongVxdID: " << n_notINdigit4);
873  B2RESULT(" noInter: " << n_notINdigit5);
874  B2RESULT("");
875 
876 
877  B2RESULT(" pT > 1 : " << pt[5]);
878  B2RESULT(" out ROI: " << nnotINdigit2[5]);
879  B2RESULT(" no ROI: " << nnotINdigit3[5]);
880  B2RESULT(" wrongVxdID: " << nnotINdigit4[5]);
881  B2RESULT(" noInter: " << nnotINdigit5[5]);
882  B2RESULT(" svdDigit : " << nsvdDigit[5]);
883  B2RESULT(" svdDigitIn : " << nsvdDigitInROI[5]);
884  if ((nsvdDigit[5] - nsvdDigitInROI[5]) != (nnotINdigit2[5] + nnotINdigit3[5] + nnotINdigit4[5] + nnotINdigit5[5]))
885  B2RESULT(" svdDigitOut : " << nsvdDigit[5] - nsvdDigitInROI[5] << " != " << nnotINdigit2[5] + nnotINdigit3[5] + nnotINdigit4[5] +
886  nnotINdigit5[5]);
887  epsilon[5] = (double)nsvdDigitInROI[5] / (double) nsvdDigit[5];
888  epsilonErr[5] = sqrt(epsilon[5] * (1 - epsilon[5]) / nsvdDigit[5]);
889  B2RESULT(" efficiency : " << epsilon[5] << " +/- " << epsilonErr[5]);
890  epsilon2[5] = (double)TrackOneDigiIn[5] / (double) totTrack[5] ;
891  epsilon2Err[5] = sqrt(epsilon2[5] * (1 - epsilon2[5]) / totTrack[5]);
892  B2RESULT(" efficiency2 : " << epsilon2[5] << " +/- " << epsilon2Err[5]);
893  B2RESULT("");
894 
895  B2RESULT(" 0.5 < pT < 1 : " << pt[4]);
896  B2RESULT(" out ROI: " << nnotINdigit2[4]);
897  B2RESULT(" no ROI: " << nnotINdigit3[4]);
898  B2RESULT(" wrongVxdID: " << nnotINdigit4[4]);
899  B2RESULT(" noInter: " << nnotINdigit5[4]);
900  B2RESULT(" svdDigit : " << nsvdDigit[4]);
901  B2RESULT(" svdDigitIn : " << nsvdDigitInROI[4]);
902  if ((nsvdDigit[4] - nsvdDigitInROI[4]) != (nnotINdigit2[4] + nnotINdigit3[4] + nnotINdigit4[4] + nnotINdigit5[4]))
903  B2RESULT(" svdDigitOut : " << nsvdDigit[4] - nsvdDigitInROI[4] << " != " << nnotINdigit2[4] + nnotINdigit3[4] + nnotINdigit4[4] +
904  nnotINdigit5[4]);
905  epsilon[4] = (double)nsvdDigitInROI[4] / (double) nsvdDigit[4];
906  epsilonErr[4] = sqrt(epsilon[4] * (1 - epsilon[4]) / nsvdDigit[4]);
907  B2RESULT(" efficiency : " << epsilon[4] << " +/- " << epsilonErr[4]);
908  epsilon2[4] = (double)TrackOneDigiIn[4] / (double) totTrack[4] ;
909  epsilon2Err[4] = sqrt(epsilon2[4] * (1 - epsilon2[4]) / totTrack[4]);
910  B2RESULT(" efficiency2 : " << epsilon2[4] << " +/- " << epsilon2Err[4]);
911 
912  B2RESULT("");
913  B2RESULT(" 0.3 < pT < 0.5 : " << pt[3]);
914  B2RESULT(" out ROI: " << nnotINdigit2[3]);
915  B2RESULT(" no ROI: " << nnotINdigit3[3]);
916  B2RESULT(" wrongVxdID: " << nnotINdigit4[3]);
917  B2RESULT(" noInter: " << nnotINdigit5[3]);
918  B2RESULT(" svdDigit : " << nsvdDigit[3]);
919  B2RESULT(" svdDigitIn : " << nsvdDigitInROI[3]);
920  if ((nsvdDigit[3] - nsvdDigitInROI[3]) != (nnotINdigit2[3] + nnotINdigit3[3] + nnotINdigit4[3] + nnotINdigit5[3]))
921  B2RESULT(" svdDigitOut : " << nsvdDigit[3] - nsvdDigitInROI[3] << " != " << nnotINdigit2[3] + nnotINdigit3[3] + nnotINdigit4[3] +
922  nnotINdigit5[3]);
923  epsilon[3] = (double)nsvdDigitInROI[3] / (double) nsvdDigit[3];
924  epsilonErr[3] = sqrt(epsilon[3] * (1 - epsilon[3]) / nsvdDigit[3]);
925  B2RESULT(" efficiency : " << epsilon[3] << " +/- " << epsilonErr[3]);
926  epsilon2[3] = (double)TrackOneDigiIn[3] / (double) totTrack[3];
927  epsilon2Err[3] = sqrt(epsilon2[3] * (1 - epsilon2[3]) / totTrack[3]);
928  B2RESULT(" efficiency2 : " << epsilon2[3] << " +/- " << epsilon2Err[3]);
929 
930  B2RESULT("");
931  B2RESULT(" 0.2 < pT < 0.3 : " << pt[2]);
932  B2RESULT(" out ROI: " << nnotINdigit2[2]);
933  B2RESULT(" no ROI: " << nnotINdigit3[2]);
934  B2RESULT(" wrongVxdID: " << nnotINdigit4[2]);
935  B2RESULT(" noInter: " << nnotINdigit5[2]);
936  B2RESULT(" svdDigit : " << nsvdDigit[2]);
937  B2RESULT(" svdDigitIn : " << nsvdDigitInROI[2]);
938  if ((nsvdDigit[2] - nsvdDigitInROI[2]) != (nnotINdigit2[2] + nnotINdigit3[2] + nnotINdigit4[2] + nnotINdigit5[2]))
939  B2RESULT(" svdDigitOut : " << nsvdDigit[2] - nsvdDigitInROI[2] << " != " << nnotINdigit2[2] + nnotINdigit3[2] + nnotINdigit4[2] +
940  nnotINdigit5[2]);
941  epsilon[2] = (double)nsvdDigitInROI[2] / (double) nsvdDigit[2];
942  epsilonErr[2] = sqrt(epsilon[2] * (1 - epsilon[2]) / nsvdDigit[2]);
943  B2RESULT(" efficiency : " << epsilon[2] << " +/- " << epsilonErr[2]);
944  epsilon2[2] = (double)TrackOneDigiIn[2] / (double) totTrack[2] ;
945  epsilon2Err[2] = sqrt(epsilon2[2] * (1 - epsilon2[2]) / totTrack[2]);
946  B2RESULT(" efficiency2 : " << epsilon2[2] << " +/- " << epsilon2Err[2]);
947 
948  B2RESULT("");
949  B2RESULT(" 0.1 < pT < 0.2 : " << pt[1]);
950  B2RESULT(" out ROI: " << nnotINdigit2[1]);
951  B2RESULT(" no ROI: " << nnotINdigit3[1]);
952  B2RESULT(" wrongVxdID: " << nnotINdigit4[1]);
953  B2RESULT(" noInter: " << nnotINdigit5[1]);
954  B2RESULT(" svdDigit : " << nsvdDigit[1]);
955  B2RESULT(" svdDigitIn : " << nsvdDigitInROI[1]);
956  if ((nsvdDigit[1] - nsvdDigitInROI[1]) != (nnotINdigit2[1] + nnotINdigit3[1] + nnotINdigit4[1] + nnotINdigit5[1]))
957  B2RESULT(" svdDigitOut : " << nsvdDigit[1] - nsvdDigitInROI[1] << " ?=? " << nnotINdigit2[1] + nnotINdigit3[1] + nnotINdigit4[1] +
958  nnotINdigit5[1]);
959  epsilon[1] = (double)nsvdDigitInROI[1] / (double) nsvdDigit[1];
960  epsilonErr[1] = sqrt(epsilon[1] * (1 - epsilon[1]) / nsvdDigit[1]);
961  B2RESULT(" efficiency : " << epsilon[1] << " +/- " << epsilonErr[1]);
962  epsilon2[1] = (double)TrackOneDigiIn[1] / (double) totTrack[1] ;
963  epsilon2Err[1] = sqrt(epsilon2[1] * (1 - epsilon2[1]) / totTrack[1]);
964  B2RESULT(" efficiency2 : " << epsilon2[1] << " +/- " << epsilon2Err[1]);
965 
966  B2RESULT("");
967  B2RESULT(" pT < 0.1 : " << pt[0]);
968  B2RESULT(" out ROI: " << nnotINdigit2[0]);
969  B2RESULT(" no ROI: " << nnotINdigit3[0]);
970  B2RESULT(" wrongVxdID: " << nnotINdigit4[0]);
971  B2RESULT(" noInter: " << nnotINdigit5[0]);
972  B2RESULT(" svdDigit : " << nsvdDigit[0]);
973  B2RESULT(" svdDigitIn : " << nsvdDigitInROI[0]);
974  if ((nsvdDigit[0] - nsvdDigitInROI[0]) != (nnotINdigit2[0] + nnotINdigit3[0] + nnotINdigit4[0] + nnotINdigit5[0]))
975  B2RESULT(" svdDigitOut : " << nsvdDigit[0] - nsvdDigitInROI[0] << " ?=? " << nnotINdigit2[0] + nnotINdigit3[0] + nnotINdigit4[0] +
976  nnotINdigit5[0]);
977  epsilon[0] = (double)nsvdDigitInROI[0] / (double) nsvdDigit[0];
978  epsilonErr[0] = sqrt(epsilon[0] * (1 - epsilon[0]) / nsvdDigit[0]);
979  B2RESULT(" efficiency : " << epsilon[0] << " +/- " << epsilonErr[0]);
980  epsilon2[0] = (double)TrackOneDigiIn[0] / (double) totTrack[0] ;
981  epsilon2Err[0] = sqrt(epsilon2[0] * (1 - epsilon2[0]) / totTrack[0]);
982  B2RESULT(" efficiency2 : " << epsilon2[0] << " +/- " << epsilon2Err[0]);
983 
984  B2RESULT("legend:");
985  B2RESULT(" CASO2: if (ROI exists but no SVDShaperDigit inside)");
986  B2RESULT(" CASO3: if (ROI does not exist, intercept with correct VxdID)");
987  B2RESULT(" CASO4: if (intercept with wrong VxdID)");
988  B2RESULT(" CASO5: if (intercept does not exist)");
989 
990  B2RESULT("==========================");
991  }
992  B2RESULT(" tot ROIs = " << n_rois);
993  B2RESULT(" ROIs with digits = " << m_nGoodROIs);
994  B2RESULT(" SVD INefficiency = " << 1 - (float)m_nGoodROIs / n_rois);
995 
996  m_gEff2 = new TGraphErrors(6, pt, epsilon2, ptErr, epsilon2Err);
997  m_gEff2->SetName("g_eff2");
998  m_gEff2->SetTitle("Normalized to MCParticles with digits and related track");
999  m_gEff = new TGraphErrors(6, pt, epsilon, ptErr, epsilonErr);
1000  m_gEff->SetName("g_eff");
1001  m_gEff->SetTitle("Normalized to digits of MCParticles with digits and related track");
1002 
1003 
1004 
1005  if (m_rootFilePtr != nullptr) {
1006  m_rootFilePtr->cd(); //important! without this the famework root I/O (SimpleOutput etc) could mix with the root I/O of this module
1007 
1008  TDirectory* oldDir = gDirectory;
1009  TDirectory* m_digiDir = oldDir->mkdir("digits");
1010  TDirectory* m_tracksDir = oldDir->mkdir("tracks");
1011  TDirectory* m_notINtrack5 = oldDir->mkdir("notINtrack5");
1012  TDirectory* m_INtrack1 = oldDir->mkdir("INtrack1");
1013  TDirectory* m_alltracks = oldDir->mkdir("alltracks");
1014  TDirectory* m_in = oldDir->mkdir("digi_in");
1015  TDirectory* m_out2 = oldDir->mkdir("digi_out2");
1016  TDirectory* m_out3 = oldDir->mkdir("digi_out3");
1017  TDirectory* m_out4 = oldDir->mkdir("digi_out4");
1018  TDirectory* m_out5 = oldDir->mkdir("digi_out5");
1019  TDirectory* m_ROIDir = oldDir->mkdir("roi");
1020 
1021  m_h1DigitsPerParticle->Write();
1022  m_h1RecoTracksPerParticle->Write();
1023 
1024  m_gEff->Write();
1025  m_gEff2->Write();
1026  m_h1effPerTrack->Write();
1027 
1028 
1029  m_digiDir->cd();
1030  m_h1digiIn->Write();
1031  m_h1digiOut2->Write();
1032  m_h1digiOut3->Write();
1033  m_h1digiOut4->Write();
1034  m_h1digiOut5->Write();
1035 
1036  m_tracksDir->cd();
1037  m_h1TrackOneDigiIn->Write();
1038  m_h1nnotINtrack2->Write();
1039  m_h1nnotINtrack3->Write();
1040  m_h1nnotINtrack4->Write();
1041  m_h1nnotINtrack5->Write();
1042 
1043  m_notINtrack5->cd();
1044  m_h1notINtrack5->Write();
1045  m_h1notINtrack5_pt->Write();
1046  m_h1notINtrack5_phi->Write();
1047  m_h1notINtrack5_lambda->Write();
1048  m_h1notINtrack5_cosTheta->Write();
1049 
1050  m_INtrack1->cd();
1051  m_h1INtrack1->Write();
1052  m_h1INtrack1_pt->Write();
1053  m_h1INtrack1_phi->Write();
1054  m_h1INtrack1_lambda->Write();
1055  m_h1INtrack1_cosTheta->Write();
1056 
1057  m_alltracks->cd();
1058  m_h1Track->Write();
1059  m_h1Track_pt->Write();
1060  m_h1Track_phi->Write();
1061  m_h1Track_lambda->Write();
1062  m_h1Track_cosTheta->Write();
1063 
1064  m_in->cd();
1065  m_h1GlobalTime->Write();
1066  m_h1PullU->Write();
1067  m_h1PullV->Write();
1068  m_h1ResidU->Write();
1069  m_h1ResidV->Write();
1070  m_h1SigmaU->Write();
1071  m_h1SigmaV->Write();
1072  m_h2sigmaUphi->Write();
1073  m_h2sigmaVphi->Write();
1074 
1075  m_out2->cd();
1076  m_h1GlobalTime_out2->Write();
1077  m_h1ResidU_out2->Write();
1078  m_h1ResidV_out2->Write();
1079  m_h1SigmaU_out2->Write();
1080  m_h1SigmaV_out2->Write();
1081  m_h2sigmaUphi_out2->Write();
1082  m_h2sigmaVphi_out2->Write();
1083 
1084  m_out3->cd();
1085  m_h1GlobalTime_out3->Write();
1086  m_h1ResidU_out3->Write();
1087  m_h1ResidV_out3->Write();
1088  m_h1SigmaU_out3->Write();
1089  m_h1SigmaV_out3->Write();
1090  m_h2sigmaUphi_out3->Write();
1091  m_h2sigmaVphi_out3->Write();
1092 
1093  m_out4->cd();
1094  m_h1GlobalTime_out4->Write();
1095  m_h1SigmaU_out4->Write();
1096  m_h1SigmaV_out4->Write();
1097  m_h2sigmaUphi_out4->Write();
1098  m_h2sigmaVphi_out4->Write();
1099 
1100  m_out5->cd();
1101  m_h1GlobalTime_out5->Write();
1102 
1103  m_ROIDir->cd();
1104  m_h1totUstrips->Write();
1105  m_h1totVstrips->Write();
1106  m_h1okROIs->Write();
1107  m_h1totROIs->Write();
1108 
1109 
1110  m_h2ROIbottomLeft->Write();
1111  m_h2ROItopRight->Write();
1112  m_h2ROIuMinMax->Write();
1113  m_h2ROIvMinMax->Write();
1114  m_rootFilePtr->Close();
1115 
1116  }
1117 
1118 }
1119 
Belle2::RelationVector::size
size_t size() const
Get number of relations.
Definition: RelationVector.h:98
Belle2::SVDIntercept
SVDIntercept stores the U,V coordinates and uncertainties of the intersection of a track with an SVD ...
Definition: SVDIntercept.h:32
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::RelationsInterface::getRelationsWith
RelationVector< T > getRelationsWith(const std::string &name="", const std::string &namedRelation="") const
Get the relations between this object and another store array.
Definition: RelationsObject.h:232
Belle2::RelationIndex::iterator_from
index_from::const_iterator iterator_from
Element iterator of the from side index.
Definition: RelationIndex.h:103
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::RelationIndex
Provides access to fast ( O(log n) ) bi-directional lookups on a specified relation.
Definition: RelationIndex.h:84
Belle2::VXD::SensorInfoBase
Base class to provide Sensor Information for PXD and SVD.
Definition: SensorInfoBase.h:40
Belle2::RelationsInterface::getRelatedTo
TO * getRelatedTo(const std::string &name="", const std::string &namedRelation="") const
Get the object to which this object has a relation.
Definition: RelationsObject.h:250
Belle2::SVDROIFinderAnalysisModule
This module performs the analysis of the SVD data reduction module performances
Definition: SVDROIFinderAnalysisModule.h:44
Belle2::SVDShaperDigit
The SVD ShaperDigit class.
Definition: SVDShaperDigit.h:46
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::RecoTrack
This is the Reconstruction Event-Data Model Track.
Definition: RecoTrack.h:78
Belle2::RelationIndex::getElementsFrom
range_from getElementsFrom(const FROM *from) const
Return a range of all elements pointing from the given object.
Definition: RelationIndex.h:157
Belle2::RelationVector
Class for type safe access to objects that are referred to in relations.
Definition: DataStore.h:38
Belle2::RelationIndex::range_from
boost::iterator_range< iterator_from > range_from
Iterator range [first,second) of the from side.
Definition: RelationIndex.h:109
Belle2::VXD::SensorInfoBase::getVCellPosition
double getVCellPosition(int vID) const
Return the position of a specific strip/pixel in v direction.
Definition: SensorInfoBase.h:191
Belle2::ROIid
ROIid stores the U and V ids and the sensor id of the Region Of Interest.
Definition: ROIid.h:35
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::RelationsInterface::getRelationsFrom
RelationVector< FROM > getRelationsFrom(const std::string &name="", const std::string &namedRelation="") const
Get the relations that point from another store array to this object.
Definition: RelationsObject.h:214
Belle2::ROIid::Contains
bool Contains(const Belle2::PXDRawHit &thePXDRawHit) const
true if the ROI contains the thePXDRawHit
Definition: ROIid.cc:27
Belle2::VXD::GeoCache
Class to faciliate easy access to sensor information of the VXD like coordinate transformations or pi...
Definition: GeoCache.h:41
Belle2::MCParticle::getMomentum
TVector3 getMomentum() const
Return momentum.
Definition: MCParticle.h:209
Belle2::VXD::SensorInfoBase::getUCellPosition
double getUCellPosition(int uID, int vID=-1) const
Return the position of a specific strip/pixel in u direction.
Definition: SensorInfoBase.h:179
Belle2::VXD::GeoCache::getSensorInfo
const SensorInfoBase & getSensorInfo(Belle2::VxdID id) const
Return a referecne to the SensorInfo of a given SensorID.
Definition: GeoCache.cc:68
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:43