Belle II Software development
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
11using namespace std;
12using namespace Belle2;
13
14REG_MODULE(SVDClusterEvaluation);
15
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();
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
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 m_rootFilePtr->Close();
606 }
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.
STL namespace.