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