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