Belle II Software  release-08-01-10
SVDClusterEvaluationModule.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 #include <svd/modules/svdPerformance/SVDClusterEvaluationModule.h>
9 #include <tracking/dataobjects/RecoTrack.h>
10 
11 using namespace std;
12 using namespace Belle2;
13 
14 REG_MODULE(SVDClusterEvaluation);
15 
16 SVDClusterEvaluationModule::SVDClusterEvaluationModule(): Module()
17  , m_interCoor(nullptr)
18  , m_interSigma(nullptr)
19  , m_clsCoor(nullptr)
20  , m_clsResid(nullptr)
21  , m_clsMinResid(nullptr)
22  , m_clsResid2D(nullptr)
23 {
24 
25  setDescription("This module check performances of SVD reconstruction of VXD TB data");
26 
27  addParam("outputFileName", m_rootFileName, "Name of output root file.", std::string("SVDClusterEvaluation_output.root"));
28 
29  addParam("LayerUnderStudy", m_theLayer, "Number of the layer under study. If 0, then all layers are plotted", int(0));
30  addParam("InterceptSigmaMax", m_interSigmaMax,
31  "Max of the histogram that contains the intercept statistical error. Default is OK for Phase2.", double(0.35));
32  addParam("uFiducialLength", m_uFiducial,
33  "length to be subtracted from the U-edge to consider intercepts inside the sensor. Positive values reduce the area; negative values increase the area",
34  double(0));
35  addParam("vFiducialLength", m_vFiducial,
36  "length to be subtracted from the V-edge to consider intercepts inside the sensor. Positive values reduce the area; negative values increase the area",
37  double(0));
38  addParam("efficiency_nSigma", m_nSigma, " number of residual sigmas for the determination of the efficiency", float(5));
39  addParam("efficiency_halfWidth", m_halfWidth, " window half width for the determination of the efficiency", float(0.05));
40  addParam("ClustersName", m_ClusterName, "Name of DUTs Cluster Store Array.", std::string(""));
41  addParam("InterceptsName", m_InterceptName, "Name of Intercept Store Array.", std::string(""));
42  addParam("TracksName", m_TrackName, "Name of Track Store Array.", std::string(""));
43  addParam("UbinWidth", m_UbinWidth, "Histograms U-bin width (in um)", double(10));
44  addParam("VbinWidth", m_VbinWidth, "Histograms V-bin width (in um)", double(10));
45  addParam("groupNstrips", m_groupNstrips, "How many strips group together in the 2D residual VS position plot", int(128));
46 }
47 
49 {
50 
51  m_eventMetaData.isRequired();
53  m_svdIntercepts.isRequired(m_InterceptName);
54  m_tracks.isRequired(m_TrackName);
55 
56  m_rootFilePtr = new TFile(m_rootFileName.c_str(), "RECREATE");
57 
60 
65 
70 
71  // create new root file
72  m_rootFilePtr = new TFile(m_rootFileName.c_str(), "RECREATE");
73 
74  //tree initialization
75  m_tree = new TTree("tree", "RECREATE");
76  b_experiment = m_tree->Branch("exp", &m_experiment, "exp/i");
77  b_run = m_tree->Branch("run", &m_run, "run/i");
78  b_layer = m_tree->Branch("layer", &m_layer, "layer/i");
79  b_ladder = m_tree->Branch("ladder", &m_ladder, "ladder/i");
80  b_sensor = m_tree->Branch("sensor", &m_sensor, "sensor/i");
81  b_interU = m_tree->Branch("interU", &m_interU, "interU/F");
82  b_interV = m_tree->Branch("interV", &m_interV, "interV/F");
83  b_interErrU = m_tree->Branch("interErrU", &m_interErrU, "interErrU/F");
84  b_interErrV = m_tree->Branch("interErrV", &m_interErrV, "interErrV/F");
85  b_interUprime = m_tree->Branch("interUprime", &m_interUprime, "interUprime/F");
86  b_interVprime = m_tree->Branch("interVprime", &m_interVprime, "interVprime/F");
87  b_interErrUprime = m_tree->Branch("interErrUprime", &m_interErrUprime, "interErrUprime/F");
88  b_interErrVprime = m_tree->Branch("interErrVprime", &m_interErrVprime, "interErrVprime/F");
89  b_residU = m_tree->Branch("residU", &m_residU, "residU/F");
90  b_residV = m_tree->Branch("residV", &m_residV, "residV/F");
91  b_clUpos = m_tree->Branch("clUpos", &m_clUpos, "clUpos/F");
92  b_clVpos = m_tree->Branch("clVpos", &m_clVpos, "clVpos/F");
93  b_clUcharge = m_tree->Branch("clUcharge", &m_clUcharge, "clUcharge/F");
94  b_clVcharge = m_tree->Branch("clVcharge", &m_clVcharge, "clVcharge/F");
95  b_clUsnr = m_tree->Branch("clUsnr", &m_clUsnr, "clUsnr/F");
96  b_clVsnr = m_tree->Branch("clVsnr", &m_clVsnr, "clVsnr/F");
97  b_clUsize = m_tree->Branch("clUsize", &m_clUsize, "clUsize/i");
98  b_clVsize = m_tree->Branch("clVsize", &m_clVsize, "clVsize/i");
99  b_clUtime = m_tree->Branch("clUtime", &m_clUtime, "clUtime/F");
100  b_clVtime = m_tree->Branch("clVtime", &m_clVtime, "clVtime/F");
101 
102  //tree initialization
103  m_treeSummary = new TTree("summary", "RECREATE");
104  bs_experiment = m_treeSummary->Branch("exp", &m_experiment, "exp/i");
105  bs_run = m_treeSummary->Branch("run", &ms_run, "run/i");
106  bs_layer = m_treeSummary->Branch("layer", &ms_layer, "layer/i");
107  bs_ladder = m_treeSummary->Branch("ladder", &ms_ladder, "ladder/i");
108  bs_sensor = m_treeSummary->Branch("sensor", &ms_sensor, "sensor/i");
109  bs_effU = m_treeSummary->Branch("effU", &ms_effU, "effU/F");
110  bs_effV = m_treeSummary->Branch("effV", &ms_effV, "effU/F");
111  bs_effErrU = m_treeSummary->Branch("effErrU", &ms_effErrU, "effErrU/F");
112  bs_effErrV = m_treeSummary->Branch("effErrV", &ms_effErrV, "effErrU/F");
113  bs_nIntercepts = m_treeSummary->Branch("nIntercepts", &ms_nIntercepts, "nIntercepts/i");
114  bs_residU = m_treeSummary->Branch("residU", &ms_residU, "residU/F");
115  bs_residV = m_treeSummary->Branch("residV", &ms_residV, "residU/F");
116  bs_misU = m_treeSummary->Branch("misU", &ms_misU, "misU/F");
117  bs_misV = m_treeSummary->Branch("misV", &ms_misV, "misU/F");
118  bs_statU = m_treeSummary->Branch("statU", &ms_statU, "statU/F");
119  bs_statV = m_treeSummary->Branch("statV", &ms_statV, "statU/F");
120 
121 }
122 
123 
125 {
126 
127  if (m_interCoor == nullptr) {
128 
129  //INTERCEPTS
131 
133 
134  //CLUSTERS
136 
138 
139  B2DEBUG(10, "Empty histograms have beein created");
140  B2DEBUG(10, "Large sensors, U side: width = " << m_width_LargeS_U << " cm, bin width = " << m_UbinWidth << " cm -> Nbins = " <<
142  B2DEBUG(10, "Large sensors, V side: width = " << m_width_LargeS_V << " cm, bin width = " << m_UbinWidth << " cm -> Nbins = " <<
144 
145  B2DEBUG(10, "Small sensors, U side: width = " << m_width_SmallS_U << " cm, bin width = " << m_UbinWidth << " cm -> Nbins = " <<
147  B2DEBUG(10, "Small sensors, V side: width = " << m_width_SmallS_V << " cm, bin width = " << m_UbinWidth << " cm -> Nbins = " <<
149  }
150 
151 
152 
153 
154 
155 
156 }
157 
159 {
160 
161  //tree variables - event
163  m_run = meta->getRun();
164  m_experiment = meta->getExperiment();
165 
166  // int nEvent = m_eventMetaData->getEvent();
167  // B2DEBUG(10, "nEvent = " << nEvent << ": n intercepts = " << m_svdIntercepts.getEntries() << "n clusters DUT = " << m_svdClusters.getEntries());
168  bool isU = true;
169 
170  //intercepts
171  for (int inter = 0 ; inter < m_svdIntercepts.getEntries(); inter++) {
172 
173  if (!isRelatedToTrack(m_svdIntercepts[inter]))
174  continue;
175 
176  B2DEBUG(10, "this intercept is related to a good track");
177 
178  VxdID::baseType theVxdID = (VxdID::baseType)m_svdIntercepts[inter]->getSensorID();
179  double coorU = m_svdIntercepts[inter]->getCoorU();
180  double coorV = m_svdIntercepts[inter]->getCoorV();
181  double sigmaU = m_svdIntercepts[inter]->getSigmaU();
182  double sigmaV = m_svdIntercepts[inter]->getSigmaV();
183 
184  //tree variables - sensor
185  m_layer = VxdID(theVxdID).getLayerNumber();
186  m_ladder = VxdID(theVxdID).getLadderNumber();
187  m_sensor = VxdID(theVxdID).getSensorNumber();
188  //tree variables - intercept
189  m_interU = coorU;
190  m_interV = coorV;
191  m_interErrU = sigmaU;
192  m_interErrV = sigmaV;
193  m_interUprime = m_svdIntercepts[inter]->getUprime();
194  m_interVprime = m_svdIntercepts[inter]->getVprime();
195  m_interErrUprime = m_svdIntercepts[inter]->getSigmaUprime();
196  m_interErrVprime = m_svdIntercepts[inter]->getSigmaVprime();
197 
198  const VXD::SensorInfoBase& theSensorInfo = m_geoCache.getSensorInfo(theVxdID);
199  if (theSensorInfo.inside(coorU, coorV, -m_uFiducial, -m_vFiducial)) {
200  B2DEBUG(10, "intercept is inside fiducial area");
201 
202  m_interCoor->fill(theVxdID, isU, coorU, coorV);
203  m_interSigma->fill(theVxdID, isU, sigmaU);
204  m_interSigma->fill(theVxdID, !isU, sigmaV);
205 
206  double minresidU = 999;
207  bool minfoundU = false;
208  double minresidV = 999;
209  bool minfoundV = false;
210  int idU = -99;
211  int idV = -99;
212  m_residU = -99;
213  m_clUpos = -99;
214  m_clUcharge = -99;
215  m_clUsnr = -99;
216  m_clUsize = -99;
217  m_clUtime = -99;
218  m_residV = -99;
219  m_clVpos = -99;
220  m_clVcharge = -99;
221  m_clVsnr = -99;
222  m_clVsize = -99;
223  m_clVtime = -99;
224 
225  //loop on clusters
226  for (int cls = 0 ; cls < m_svdClusters.getEntries(); cls++) {
227 
228  VxdID::baseType clVxdID = (VxdID::baseType)m_svdClusters[cls]->getSensorID();
229  if (clVxdID != theVxdID)
230  continue;
231 
232  double interCoor = coorV;
233  double clPos = m_svdClusters[cls]->getPosition();
234  // double interSigma = sigmaV;
235  if (m_svdClusters[cls]->isUCluster()) {
236  interCoor = coorU;
237  //interSigma = sigmaU;
238  clPos = m_svdClusters[cls]->getPosition(coorV);
239  }
240  double resid = interCoor - clPos;
241  m_clsResid->fill(theVxdID, m_svdClusters[cls]->isUCluster(), resid);
242  m_clsResid2D->fill(theVxdID, m_svdClusters[cls]->isUCluster(), clPos, resid);
243 
244  //looking for the minimal residual
245  if (m_svdClusters[cls]->isUCluster()) {
246  if (fabs(resid) < fabs(minresidU)) {
247  minfoundU = true;
248  minresidU = resid;
249  idU = cls;
250  }
251  } else {
252  if (fabs(resid) < fabs(minresidV)) {
253  minfoundV = true;
254  minresidV = resid;
255  idV = cls;
256  }
257  }
258  }
259  if (minfoundU) {
260  m_clsMinResid->fill(theVxdID, true, minresidU);
261  m_residU = minresidU;
262  m_clUpos = m_svdClusters[idU]->getPosition(coorV);
263  m_clUcharge = m_svdClusters[idU]->getCharge();
264  m_clUsnr = m_svdClusters[idU]->getSNR();
265  m_clUsize = (int)m_svdClusters[idU]->getSize();
266  m_clUtime = m_svdClusters[idU]->getClsTime();
267  }
268  if (minfoundV) {
269  m_clsMinResid->fill(theVxdID, false, minresidV);
270  m_residV = minresidV;
271  m_clVpos = m_svdClusters[idV]->getPosition();
272  m_clVcharge = m_svdClusters[idV]->getCharge();
273  m_clVsnr = m_svdClusters[idV]->getSNR();
274  m_clVsize = (int)m_svdClusters[idV]->getSize();
275  m_clVtime = m_svdClusters[idV]->getClsTime();
276  }
277 
278  //fill only if inside fiducial area
279  m_tree->Fill();
280 
281  }
282 
283 
284  }
285 
286  //clusters
287  for (int cls = 0 ; cls < m_svdClusters.getEntries(); cls++) {
288 
289  VxdID::baseType theVxdID = (VxdID::baseType)m_svdClusters[cls]->getSensorID();
290  m_clsCoor->fill(theVxdID, m_svdClusters[cls]->isUCluster(), m_svdClusters[cls]->getPosition());
291 
292  }
293 
294 
295 }
296 
297 
299 {
300 
302  ms_run = meta->getRun();
303  ms_experiment = meta->getExperiment();
304 
305  if (m_rootFilePtr != nullptr) {
306  m_rootFilePtr->cd();
307 
308  m_tree->Write();
309 
310  const int Nsensors = 172;//L6
311  // float sensors[Nsensors]; //sensor identificator
312  // float sensorsErr[Nsensors]; //sensor identificator
313  float residU[Nsensors]; //U residuals
314  float residV[Nsensors]; //V residuals
315  float misU[Nsensors]; //U misalignment
316  float misV[Nsensors]; //V misalignment
317  // float resolU[Nsensors]; //U residuals
318  // float resolV[Nsensors]; //V residuals
319  float effU[Nsensors];
320  float effV[Nsensors];
321  float effUErr[Nsensors];
322  float effVErr[Nsensors];
323  TString sensorU[Nsensors];
324  TString sensorV[Nsensors];
325 
326  for (int i = 0; i < Nsensors; i++) {
327  // sensors[i] = i;
328  // sensorsErr[i] = 0;
329  residU[i] = 0;
330  residV[i] = 0;
331  misU[i] = 0;
332  misV[i] = 0;
333  // resolU[i] = 0;
334  // resolV[i] = 0;
335  effU[i] = -1;
336  effV[i] = -1;
337  effUErr[i] = 0;
338  effVErr[i] = 0;
339  sensorU[i] = "";
340  sensorV[i] = "";
341  }
342 
343  TH1F* h_residU = new TH1F("hResidU", "U Residuals", 1, 0, 1);
344  h_residU->SetCanExtend(TH1::kAllAxes);
345  h_residU->SetStats(0);
346  h_residU->GetXaxis()->SetTitle("sensor");
347  h_residU->GetYaxis()->SetTitle("U residuals (#mum)");
348  TH1F* h_residV = new TH1F("hResidV", "V Residuals", 1, 0, 1);
349  h_residV->SetCanExtend(TH1::kAllAxes);
350  h_residV->SetStats(0);
351  h_residV->GetXaxis()->SetTitle("sensor");
352  h_residV->GetYaxis()->SetTitle("V residuals (#mum)");
353 
354  TH1F* h_statU = new TH1F("hStatU", "U Intercept Statistical Error", 1, 0, 1);
355  h_statU->SetCanExtend(TH1::kAllAxes);
356  h_statU->SetStats(0);
357  h_statU->GetXaxis()->SetTitle("sensor");
358  h_statU->GetYaxis()->SetTitle("U extrap. error (#mum)");
359  TH1F* h_statV = new TH1F("hStatV", "V Intercept Statistical Error", 1, 0, 1);
360  h_statV->SetCanExtend(TH1::kAllAxes);
361  h_statV->SetStats(0);
362  h_statV->GetXaxis()->SetTitle("sensor");
363  h_statV->GetYaxis()->SetTitle("V extrap. error (#mum)");
364 
365  TH1F* h_misU = new TH1F("hMisU", "U Residual Misalignment", 1, 0, 1);
366  h_misU->SetCanExtend(TH1::kAllAxes);
367  h_misU->SetStats(0);
368  h_misU->GetXaxis()->SetTitle("sensor");
369  h_misU->GetYaxis()->SetTitle("U misalignment (#mum)");
370  TH1F* h_misV = new TH1F("hMisV", "V Residual Misalignment", 1, 0, 1);
371  h_misV->SetCanExtend(TH1::kAllAxes);
372  h_misV->SetStats(0);
373  h_misV->GetXaxis()->SetTitle("sensor");
374  h_misV->GetYaxis()->SetTitle("V misalignment (#mum)");
375 
376 
377  TH1F* h_effU = new TH1F("hEffU", Form("U-Side Summary, %.1f#sigma or #pm%.1f mm", m_nSigma, m_halfWidth * 10), 1, 0, 1);
378  h_effU->SetCanExtend(TH1::kAllAxes);
379  h_effU->SetStats(0);
380  h_effU->GetXaxis()->SetTitle("sensor");
381  h_effU->GetYaxis()->SetTitle("U efficiency");
382  TH1F* h_effV = new TH1F("hEffV", Form("V-Side Summary, %.1f#sigma or #pm%.1f mm", m_nSigma, m_halfWidth * 10), 1, 0, 1);
383  h_effV->SetCanExtend(TH1::kAllAxes);
384  h_effV->SetStats(0);
385  h_effV->GetXaxis()->SetTitle("sensor");
386  h_effV->GetYaxis()->SetTitle("V efficiency");
387 
388  TDirectory* oldDir = gDirectory;
389 
390  int s = 0; //sensor counter;
391 
392  for (auto layer : m_geoCache.getLayers(VXD::SensorInfoBase::SVD)) {
393  int currentLayer = layer.getLayerNumber();
394  ms_layer = currentLayer;
395 
396  if (m_theLayer != 0 && currentLayer != m_theLayer)
397  continue;
398 
399  TString interName = Form("interceptsL%d", layer.getLayerNumber());
400  TString clsName = Form("clustersL%d", layer.getLayerNumber());
401  TString residName = Form("residualsL%d", layer.getLayerNumber());
402  TDirectory* dir_inter = oldDir->mkdir(interName.Data());
403  TDirectory* dir_cls = oldDir->mkdir(clsName.Data());
404  TDirectory* dir_resid = oldDir->mkdir(residName.Data());
405  for (auto ladder : m_geoCache.getLadders(layer))
406  for (Belle2::VxdID sensor : m_geoCache.getSensors(ladder)) {
407  ms_ladder = (VxdID)sensor.getLadderNumber();
408  ms_sensor = (VxdID)sensor.getSensorNumber();
409 
410  dir_inter->cd();
411  (m_interCoor->getHistogram(sensor, 1))->Write();
412  for (int view = SVDHistograms<TH1F>::VIndex ; view < SVDHistograms<TH1F>::UIndex + 1; view++) {
413  dir_cls->cd();
414  (m_clsCoor->getHistogram(sensor, view))->Write();
415 
416  dir_inter->cd();
417  float stat = (m_interSigma->getHistogram(sensor, view))->GetMean();
418  ms_nIntercepts = (m_interSigma->getHistogram(sensor, view))->GetEntries();
419  (m_interSigma->getHistogram(sensor, view))->Write();
420 
421  dir_resid->cd();
422  TH1F* res = m_clsMinResid->getHistogram(sensor, view);
423  Double_t median, q;
424  q = 0.5; // 0.5 for "median"
425  {
426  if (view == SVDHistograms<TH1F>::UIndex) {
427  sensorU[s] = Form("%d.%d.%dU", currentLayer, ladder.getLadderNumber(), sensor.getSensorNumber());
428  B2DEBUG(10, "U-side efficiency for " << currentLayer << "." << ladder.getLadderNumber() << "." << sensor.getSensorNumber());
429 
430  if (res->GetEntries() > 0) {
431 
432  res->GetQuantiles(1, &median, &q);
433 
434  residU[s] = getOneSigma(res);
435  misU[s] = median;
436 
437  float halfWindow = m_nSigma * residU[s];
438  if (m_nSigma == 0)
439  halfWindow = m_halfWidth;
440 
441  int binMin = res->FindBin(misU[s] - halfWindow);
442  int binMax = res->FindBin(misU[s] + halfWindow);
443  B2DEBUG(10, "from " << misU[s] - halfWindow << " -> binMin = " << binMin);
444  B2DEBUG(10, "to " << misU[s] + halfWindow << " -> binMax = " << binMax);
445 
446  int num = 0;
447  for (int bin = binMin; bin < binMax + 1; bin++)
448  num = num + res->GetBinContent(bin);
449 
450  float bkg = 0;
451  for (int bin = 1; bin < binMin; bin++)
452  bkg = bkg + res->GetBinContent(bin);
453  for (int bin = binMax; bin < res->GetNbinsX() + 1; bin++)
454  bkg = bkg + res->GetBinContent(bin);
455  //remove background clusters estimated from sidebands
456  num = num - bkg * (binMax - binMin + 1.) / (binMin + res->GetNbinsX() - binMax - 1);
457 
458  if (ms_nIntercepts > 0) {
459  effU[s] = 1.*num / ms_nIntercepts;
460  //filling efficiency histogram
461  h_effU->Fill(sensorU[s], effU[s]);
462  if (effU[s] > 1)
463  B2WARNING("something is wrong! efficiency greater than 1: " << num << "/" << ms_nIntercepts);
464  effUErr[s] = sqrt(effU[s] * (1 - effU[s]) / ms_nIntercepts);
465  }
466  B2DEBUG(10, "num = " << num);
467  B2DEBUG(10, "den = " << ms_nIntercepts);
468  B2RESULT("U-side efficiency for " << currentLayer << "." << ladder.getLadderNumber() << "." << sensor.getSensorNumber() << " = " <<
469  effU[s] << " ± " << effUErr[s]);
470 
471  //filling summary Histograms for the U side
472  h_statU->Fill(sensorU[s], stat * m_cmTomicron);
473  h_residU->Fill(sensorU[s], residU[s]*m_cmTomicron);
474  h_misU->Fill(sensorU[s], misU[s]*m_cmTomicron);
475 
476  //finally set branch values
477  ms_residU = residU[s];
478  ms_misU = misU[s];
479  ms_effU = effU[s];
480  ms_effErrU = effUErr[s];
481  ms_statU = stat;
482  } else {
483  //set to some values if residual histogram is empty
484  ms_residU = -99;
485  ms_misU = -99;
486  ms_effU = -99;
487  ms_effErrU = -99;
488  ms_statU = -99;
489  }
490 
491 
492  } else { // V-side
493  sensorV[s] = Form("%d.%d.%dV", currentLayer, ladder.getLadderNumber(), sensor.getSensorNumber());
494  B2DEBUG(10, "V-side efficiency for " << currentLayer << "." << ladder.getLadderNumber() << "." << sensor.getSensorNumber());
495 
496  if (res->GetEntries() > 0) {
497 
498  res->GetQuantiles(1, &median, &q);
499 
500  residV[s] = getOneSigma(res);
501  misV[s] = median;
502  ms_misV = misV[s];
503 
504  float halfWindow = m_nSigma * residV[s];
505  if (m_nSigma == 0)
506  halfWindow = m_halfWidth;
507 
508  int binMin = res->FindBin(misV[s] - halfWindow);
509  int binMax = res->FindBin(misV[s] + halfWindow);
510  B2DEBUG(10, "from " << misV[s] - halfWindow << " -> binMin = " << binMin);
511  B2DEBUG(10, "to " << misV[s] + halfWindow << " -> binMax = " << binMax);
512 
513  //determine signal clusters
514  int num = 0;
515  for (int bin = binMin; bin < binMax + 1; bin++)
516  num = num + res->GetBinContent(bin);
517 
518  float bkg = 0;
519  for (int bin = 1; bin < binMin; bin++)
520  bkg = bkg + res->GetBinContent(bin);
521  for (int bin = binMax; bin < res->GetNbinsX() + 1; bin++)
522  bkg = bkg + res->GetBinContent(bin);
523  //remove background clusters estimated from sidebands
524  num = num - bkg * (binMax - binMin + 1.) / (binMin + res->GetNbinsX() - binMax - 1);
525 
526  if (ms_nIntercepts > 0) {
527  effV[s] = 1.*num / ms_nIntercepts;
528  //filling efficiency histogram
529  h_effV->Fill(sensorV[s], effV[s]);
530  if (effV[s] > 1)
531  B2WARNING("something is wrong! efficiency greater than 1: " << num << "/" << ms_nIntercepts);
532  effVErr[s] = sqrt(effV[s] * (1 - effV[s]) / ms_nIntercepts);
533  }
534  B2DEBUG(10, "num = " << num);
535  B2DEBUG(10, "den = " << ms_nIntercepts);
536  B2RESULT("V-side efficiency for " << currentLayer << "." << ladder.getLadderNumber() << "." << sensor.getSensorNumber() << " = " <<
537  effV[s] << " ± " << effVErr[s]);
538 
539  //filling summary Histograms for the V side
540  h_statV->Fill(sensorV[s], stat * m_cmTomicron);
541  h_residV->Fill(sensorV[s], residV[s]*m_cmTomicron);
542  h_misV->Fill(sensorV[s], misV[s]*m_cmTomicron);
543 
544  //finally set branch values
545  ms_residV = residV[s];
546  ms_misV = misV[s];
547  ms_effV = effV[s];
548  ms_effErrV = effVErr[s];
549  ms_statV = stat;
550  } else {
551  //set to some values if residual histogram is empty
552  ms_residV = -99;
553  ms_misV = -99;
554  ms_effV = -99;
555  ms_effErrV = -99;
556  ms_statV = -99;
557  }
558 
559  }
560 
561  }
562  B2DEBUG(50, "writing out resid histograms for " << sensor.getLayerNumber() << "." << sensor.getLadderNumber() << "." <<
563  sensor.getSensorNumber() << "." << view);
564  (m_clsResid->getHistogram(sensor, view))->Write();
565  (res)->Write();
566  (m_clsResid2D->getHistogram(sensor, view))->Write();
567 
568 
569  }
570  m_treeSummary->Fill();
571  s++;
572  }
573  }
574 
575 
576 
577  oldDir->cd();
578  m_treeSummary->Write();
579 
580  for (int bin = 0; bin < h_residU->GetNbinsX(); bin++)
581  h_residU->SetBinError(bin, 0.);
582  h_residU->Write();
583  for (int bin = 0; bin < h_residV->GetNbinsX(); bin++)
584  h_residV->SetBinError(bin, 0.);
585  h_residV->Write();
586  for (int bin = 0; bin < h_statU->GetNbinsX(); bin++)
587  h_statU->SetBinError(bin, 0.);
588  h_statU->Write();
589  for (int bin = 0; bin < h_statV->GetNbinsX(); bin++)
590  h_statV->SetBinError(bin, 0.);
591  h_statV->Write();
592  for (int bin = 0; bin < h_misU->GetNbinsX(); bin++)
593  h_misU->SetBinError(bin, 0.);
594  h_misU->Write();
595  for (int bin = 0; bin < h_misV->GetNbinsX(); bin++)
596  h_misV->SetBinError(bin, 0.);
597  h_misV->Write();
598  for (int bin = 0; bin < h_effU->GetNbinsX(); bin++)
599  h_effU->SetBinError(bin, 0.);
600  h_effU->Write();
601  for (int bin = 0; bin < h_effV->GetNbinsX(); bin++)
602  h_effV->SetBinError(bin, 0.);
603  h_effV->Write();
604  }
605 
606  m_rootFilePtr->Close();
607 }
608 
609 
610 
612 {
613 
614  RelationVector<RecoTrack> theRC = DataStore::getRelationsWithObj<RecoTrack>(inter);
615  if (theRC.size() == 0)
616  return false;
617 
618  RelationVector<Track> theTrack = theRC[0]->getRelationsWith<Track>(m_TrackName);
619 
620  if (theTrack.size() == 0)
621  return false;
622 
623  return true;
624 
625 }
626 
627 
629 {
630 
631  TH2F h_coorUV_LargeSensor("interCoor_Large_L@layerL@ladderS@sensor",
632  "Intercept 2D Coordinate (layer @layer, ladder @ladder, sensor @sensor)",
635  h_coorUV_LargeSensor.GetXaxis()->SetTitle("Intercept U coordinate (cm)");
636  h_coorUV_LargeSensor.GetYaxis()->SetTitle("Intercept V coordinate (cm)");
637 
638  TH2F h_coorUV_SmallSensor("interCoor_Small_L@layerL@ladderS@sensor",
639  "Intercept 2D Coordinate (layer @layer, ladder @ladder, sensor @sensor)",
642  h_coorUV_SmallSensor.GetXaxis()->SetTitle("Intercept U coordinate (cm)");
643  h_coorUV_SmallSensor.GetYaxis()->SetTitle("Intercept V coordinate (cm)");
644 
645 
646  m_interCoor = new SVDHistograms<TH2F>(h_coorUV_SmallSensor, h_coorUV_SmallSensor, h_coorUV_LargeSensor, h_coorUV_LargeSensor);
647 }
648 
649 
651 {
652 
653  TH1F h_sigmaU("interSigmaU_L@layerL@ladderS@sensor@view",
654  "U Intercept Sigma (layer @layer, ladder @ladder, sensor @sensor)",
655  100, 0, m_interSigmaMax);
656  h_sigmaU.GetXaxis()->SetTitle("Intercept U Error (cm)");
657 
658  TH1F h_sigmaV("interSigmaV_L@layerL@ladderS@sensor@view",
659  "V Intercept Sigma (layer @layer, ladder @ladder, sensor @sensor)",
660  100, 0, m_interSigmaMax);
661  h_sigmaV.GetXaxis()->SetTitle("Intercept V Error (cm)");
662 
663 
664  m_interSigma = new SVDHistograms<TH1F>(h_sigmaU, h_sigmaV, h_sigmaU, h_sigmaV);
665 }
666 
667 
669 {
670 
671 
672  TH1F h_clcoorU_LargeSensor("clsCoorU_LS_L@layerL@ladderS@sensor@view",
673  "Cluster U Coordinate (layer @layer, ladder @ladder, sensor @sensor, side@view/@side)",
675  h_clcoorU_LargeSensor.GetXaxis()->SetTitle("Cluster U coordinate (cm)");
676 
677  TH1F h_clcoorV_LargeSensor("clsCoorV_LS_L@layerL@ladderS@sensor@view",
678  "Cluster V Coordinate (layer @layer, ladder @ladder, sensor @sensor, side@view/@side)",
680  h_clcoorV_LargeSensor.GetXaxis()->SetTitle("Cluster V coordinate (cm)");
681 
682  TH1F h_clcoorU_SmallSensor("clsCoorU_SS_L@layerL@ladderS@sensor@view",
683  "Cluster U Coordinate (layer @layer, ladder @ladder, sensor @sensor, side@view/@side)",
685  h_clcoorU_SmallSensor.GetXaxis()->SetTitle("Cluster U coordinate (cm)");
686 
687  TH1F h_clcoorV_SmallSensor("clsCoorV_SS_L@layerL@ladderS@sensor@view",
688  "Cluster V Coordinate (layer @layer, ladder @ladder, sensor @sensor, side@view/@side)",
690  h_clcoorV_SmallSensor.GetXaxis()->SetTitle("Cluster V coordinate (cm)");
691 
692 
693  m_clsCoor = new SVDHistograms<TH1F>(h_clcoorU_SmallSensor, h_clcoorV_SmallSensor, h_clcoorU_LargeSensor, h_clcoorV_LargeSensor);
694 
695 }
696 
698 {
699 
700  float range = 0.5;
701  int NbinsU = 200;//range*0.0001*2/m_UbinWidth;
702  int NbinsV = 200;//range*0.0001*2/m_VbinWidth;
703 
704  //CLUSTER RESIDUALS
705  TH1F h_clresidU_LargeSensor("clsResidU_LS_L@layerL@ladderS@sensor@view",
706  "U Cluster Residuals (layer @layer, ladder @ladder, sensor @sensor, side@view/@side)",
707  // m_nBins_LargeS_U, -m_abs_LargeS_U, m_abs_LargeS_U);
708  NbinsU, -range, range);
709  h_clresidU_LargeSensor.GetXaxis()->SetTitle("residual (cm)");
710 
711  TH1F h_clresidV_LargeSensor("clsResidV_LS_L@layerL@ladderS@sensor@view",
712  "V Cluster Residuals (layer @layer, ladder @ladder, sensor @sensor, side@view/@side)",
713  // m_nBins_LargeS_V, -m_abs_LargeS_V, m_abs_LargeS_V);
714  NbinsV, -range, range);
715  h_clresidV_LargeSensor.GetXaxis()->SetTitle("residual (cm)");
716 
717  TH1F h_clresidU_SmallSensor("clsResidU_SS_L@layerL@ladderS@sensor@view",
718  "U Cluster Residuals (layer @layer, ladder @ladder, sensor @sensor, side@view/@side)",
719  // m_nBins_SmallS_U, -m_abs_SmallS_U, m_abs_SmallS_U);
720  NbinsU, -range, range);
721  h_clresidU_SmallSensor.GetXaxis()->SetTitle("residual (cm)");
722 
723  TH1F h_clresidV_SmallSensor("clsResidV_SS_L@layerL@ladderS@sensor@view",
724  "V Cluster Residuals (layer @layer, ladder @ladder, sensor @sensor, side@view/@side)",
725  // m_nBins_SmallS_V, -m_abs_SmallS_V, m_abs_SmallS_V);
726  NbinsU, -range, range);
727  h_clresidV_SmallSensor.GetXaxis()->SetTitle("residual (cm)");
728 
729 
730 
731  m_clsResid = new SVDHistograms<TH1F>(h_clresidU_SmallSensor, h_clresidV_SmallSensor, h_clresidU_LargeSensor,
732  h_clresidV_LargeSensor);
733 
734  //CLUSTER RESIDUALS VS CL POSITION
735  const int Nzones_768 = 768 / m_groupNstrips;
736  const int Nzones_512 = 512 / m_groupNstrips;
737 
738  TH2F h2_clresidU_LargeSensor("clsResid2DU_LS_L@layerL@ladderS@sensor@view",
739  "U Cluster Residuals VS U Cluster Position(layer @layer, ladder @ladder, sensor @sensor, side@view/@side)",
740  Nzones_768, -m_width_LargeS_U / 2, m_width_LargeS_U / 2, NbinsU, -range, range);
741  h2_clresidU_LargeSensor.GetYaxis()->SetTitle("residual (cm)");
742  h2_clresidU_LargeSensor.GetXaxis()->SetTitle("cluster position (cm)");
743 
744  TH2F h2_clresidV_LargeSensor("clsResid2DV_LS_L@layerL@ladderS@sensor@view",
745  "V Cluster Residuals (layer @layer, ladder @ladder, sensor @sensor, side@view/@side)",
746  Nzones_512, -m_width_LargeS_V / 2, m_width_LargeS_V / 2, NbinsV, -range, range);
747  h2_clresidV_LargeSensor.GetYaxis()->SetTitle("residual (cm)");
748  h2_clresidV_LargeSensor.GetXaxis()->SetTitle("cluster position (cm)");
749 
750  TH2F h2_clresidU_SmallSensor("clsResid2DU_SS_L@layerL@ladderS@sensor@view",
751  "U Cluster Residuals (layer @layer, ladder @ladder, sensor @sensor, side@view/@side)",
752  Nzones_768, -m_width_SmallS_U / 2, m_width_SmallS_U / 2, NbinsU, -range, range);
753  h2_clresidU_SmallSensor.GetYaxis()->SetTitle("residual (cm)");
754  h2_clresidU_SmallSensor.GetXaxis()->SetTitle("cluster position (cm)");
755 
756  TH2F h2_clresidV_SmallSensor("clsResid2DV_SS_L@layerL@ladderS@sensor@view",
757  "V Cluster Residuals (layer @layer, ladder @ladder, sensor @sensor, side@view/@side)",
758  Nzones_512, -m_width_SmallS_V / 2, m_width_SmallS_V / 2, NbinsU, -range, range);
759  h2_clresidV_SmallSensor.GetYaxis()->SetTitle("residual (cm)");
760  h2_clresidV_SmallSensor.GetXaxis()->SetTitle("cluster position (cm)");
761 
762  m_clsResid2D = new SVDHistograms<TH2F>(h2_clresidU_SmallSensor, h2_clresidV_SmallSensor, h2_clresidU_LargeSensor,
763  h2_clresidV_LargeSensor);
764 
765  //CLUSTER MINIMUM RESIDUAL
766  //CLUSTER RESIDUALS
767  TH1F h_clminresidU_LargeSensor("clsMinResidU_LS_L@layerL@ladderS@sensor@view",
768  "U Cluster Residuals (layer @layer, ladder @ladder, sensor @sensor, side@view/@side)",
769  // m_nBins_LargeS_U, -m_abs_LargeS_U, m_abs_LargeS_U);
770  NbinsU, -range, range);
771  h_clminresidU_LargeSensor.GetXaxis()->SetTitle("residual (cm)");
772 
773  TH1F h_clminresidV_LargeSensor("clsMinResidV_LS_L@layerL@ladderS@sensor@view",
774  "V Cluster Residuals (layer @layer, ladder @ladder, sensor @sensor, side@view/@side)",
775  // m_nBins_LargeS_V, -m_abs_LargeS_V, m_abs_LargeS_V);
776  NbinsV, -range, range);
777  h_clminresidV_LargeSensor.GetXaxis()->SetTitle("residual (cm)");
778 
779  TH1F h_clminresidU_SmallSensor("clsMinResidU_SS_L@layerL@ladderS@sensor@view",
780  "U Cluster Residuals (layer @layer, ladder @ladder, sensor @sensor, side@view/@side)",
781  // m_nBins_SmallS_U, -m_abs_SmallS_U, m_abs_SmallS_U);
782  NbinsU, -range, range);
783  h_clminresidU_SmallSensor.GetXaxis()->SetTitle("residual (cm)");
784 
785  TH1F h_clminresidV_SmallSensor("clsMinResidV_SS_L@layerL@ladderS@sensor@view",
786  "V Cluster Residuals (layer @layer, ladder @ladder, sensor @sensor, side@view/@side)",
787  // m_nBins_SmallS_V, -m_abs_SmallS_V, m_abs_SmallS_V);
788  NbinsU, -range, range);
789  h_clminresidV_SmallSensor.GetXaxis()->SetTitle("residual (cm)");
790 
791 
792 
793  m_clsMinResid = new SVDHistograms<TH1F>(h_clminresidU_SmallSensor, h_clminresidV_SmallSensor, h_clminresidU_LargeSensor,
794  h_clminresidV_LargeSensor);
795 
796 
797 }
798 
799 
801 {
802 
803  TH1F* h1_res = (TH1F*)h1->Clone("h1_res");
804  double probs[2] = {0.16, 1 - 0.16};
805  double quant[2] = {0, 0};
806  int nbinsHisto = h1_res->GetNbinsX();
807  h1_res->SetBinContent(1, h1_res->GetBinContent(0) + h1_res->GetBinContent(1));
808  h1_res->SetBinContent(nbinsHisto, h1_res->GetBinContent(nbinsHisto) + h1_res->GetBinContent(nbinsHisto + 1));
809  h1_res->SetBinContent(0, 0);
810  h1_res->SetBinContent(nbinsHisto + 1, 0);
811  h1_res->GetQuantiles(2, quant, probs);
812 
813  return (-quant[0] + quant[1]) / 2;
814 }
Base class for Modules.
Definition: Module.h:72
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
Class for type safe access to objects that are referred to in relations.
size_t size() const
Get number of relations.
TBranch * b_interErrV
intercept V position error
TTree * m_tree
pointer at tree containing the parameters
TBranch * bs_statV
intercept stat error V
int m_groupNstrips
number of strip in the group in 2D resid vs position
VXD::GeoCache & m_geoCache
the geo cache instance
float m_abs_LargeS_U
half width including safety margin, large sensor U
virtual void initialize() override
check StoreArrays & create rootfile
SVDHistograms< TH1F > * m_clsResid
cluster resid plots
double getOneSigma(TH1F *h)
returns one sigma using quantiles
double m_VbinWidth
histogram v-bin width in microns
void create_SVDHistograms_interCoor()
create intercept coordinates plots
StoreArray< SVDCluster > m_svdClusters
clusters
virtual void event() override
fill histograms
bool isRelatedToTrack(SVDIntercept *inter)
is the intercept related to a track
TBranch * b_interV
intercept V position
void create_SVDHistograms_interSigma()
create intercept error plots
float m_halfWidth
window half width for efficiency computation
void create_SVDHistograms_clsResid()
create slucter resid plots
virtual void endRun() override
write histogrmas
TBranch * b_interErrUprime
intercept U prime error
TBranch * b_interU
intercept U position
TBranch * b_interErrU
intercept U position error
float m_abs_SmallS_V
half width including safety margin, small sensor V
SVDHistograms< TH1F > * m_interSigma
intercept stat error plots
StoreObjPtr< EventMetaData > m_eventMetaData
event meta data
float m_nSigma
number of sigmas for efficiency computation
TBranch * bs_statU
intercept stat error U
SVDHistograms< TH2F > * m_clsResid2D
2D resid plots
std::string m_TrackName
Track StoreArray name.
TBranch * b_interErrVprime
intercept V prime error
SVDHistograms< TH1F > * m_clsMinResid
cluster minimum resid plots
virtual void beginRun() override
create histograms
double m_interSigmaMax
max of the histo of the intercept stat error
float m_interErrUprime
intercept U prime error
SVDHistograms< TH1F > * m_clsCoor
cluster coordinates plots
TTree * m_treeSummary
pointer at tree containing the summary parameters
SVDHistograms< TH2F > * m_interCoor
intercept coordinates plots
float m_interErrVprime
intercept V prime error
std::string m_ClusterName
SVDCluster StoreArray name.
int m_nBins_SmallS_U
number of bins for small sensor U
float m_interErrU
intercept U position error
float m_abs_LargeS_V
half width including safety margin, large sensor V
double m_UbinWidth
histogram u-bin width in microns
StoreArray< SVDIntercept > m_svdIntercepts
intercepts
int m_nBins_LargeS_U
number of bins for large sensor U
float m_interErrV
intercept V position error
TFile * m_rootFilePtr
pointer at root file used for storing histograms
int m_nBins_LargeS_V
number of bins for large sensor V
int m_nBins_SmallS_V
number of bins for small sensor V
float m_abs_SmallS_U
half width including safety margin, small sensor U
std::string m_InterceptName
SVDIntercept StoreArray name.
TBranch * bs_nIntercepts
number of intercepts
void create_SVDHistograms_clsCoor()
create cluster coordinates plots
template class for SVd histograms
Definition: SVDHistograms.h:24
void fill(const VxdID &vxdID, int view, Types ... args)
fill the histogram for
Definition: SVDHistograms.h:77
H * getHistogram(const VxdID &vxdID, int view)
get a reference to the histogram for
Definition: SVDHistograms.h:56
SVDIntercept stores the U,V coordinates and uncertainties of the intersection of a track with an SVD ...
Definition: SVDIntercept.h:22
bool isRequired(const std::string &name="")
Ensure this array/object has been registered previously.
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:216
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:96
Class that bundles various TrackFitResults.
Definition: Track.h:25
const std::set< Belle2::VxdID > getLayers(SensorInfoBase::SensorType sensortype=SensorInfoBase::VXD)
Return a set of all known Layers.
Definition: GeoCache.cc:176
const SensorInfoBase & getSensorInfo(Belle2::VxdID id) const
Return a referecne to the SensorInfo of a given SensorID.
Definition: GeoCache.cc:67
const std::set< Belle2::VxdID > & getSensors(Belle2::VxdID ladder) const
Return a set of all sensor IDs belonging to a given ladder.
Definition: GeoCache.cc:204
const std::set< Belle2::VxdID > & getLadders(Belle2::VxdID layer) const
Return a set of all ladder IDs belonging to a given layer.
Definition: GeoCache.cc:193
Base class to provide Sensor Information for PXD and SVD.
bool inside(double u, double v, double uTolerance=DBL_EPSILON, double vTolerance=DBL_EPSILON) const
Check wether a given point is inside the active area.
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:33
baseType getSensorNumber() const
Get the sensor id.
Definition: VxdID.h:100
unsigned short baseType
The base integer type for VxdID.
Definition: VxdID.h:36
baseType getLadderNumber() const
Get the ladder id.
Definition: VxdID.h:98
baseType getLayerNumber() const
Get the layer id.
Definition: VxdID.h:96
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition: Module.h:560
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Abstract base class for different kinds of events.