Belle II Software  release-05-02-19
SVDChannelMappingModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Giulia Casarosa, Eugenio Paoloni *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <svd/modules/svdCalibration/SVDChannelMappingModule.h>
12 #include <vxd/geometry/GeoCache.h>
13 
14 #include <TH2F.h>
15 
16 using namespace std;
17 using namespace Belle2;
18 
19 //-----------------------------------------------------------------
20 // Register the Module
21 //-----------------------------------------------------------------
22 REG_MODULE(SVDChannelMapping)
23 
24 //-----------------------------------------------------------------
25 // Implementation
26 //-----------------------------------------------------------------
27 
29  : Module()
30  , hInterDictionary(172, [](const Belle2::VxdID & vxdid) {return (size_t)vxdid.getID(); })
31 , n_events(0)
32 {
33  //Set module properties
34  setDescription("SVD Channel Mapping Verification Module");
35  setPropertyFlags(c_ParallelProcessingCertified);
36 
37  addParam("outputFileName", m_rootFileName, "Name of output root file.", std::string("SVDChannelMapping.root"));
38 
39  addParam("ShaperDigitsName", m_SVDShaperDigitsName,
40  "name of the list of SVDShaperDigits", std::string(""));
41 
42  addParam("ClustersName", m_SVDClustersName,
43  "name of the list of SVDClusters", std::string(""));
44 
45  addParam("InterceptsName", m_InterceptsName,
46  "name of the list of interceptions", std::string(""));
47 
48 }
49 
50 
51 
52 void SVDChannelMappingModule::initialize()
53 {
54 
55  m_histoList_digits = new TList;
56  m_histoList_clusters = new TList;
57 
58  createHistosDictionaries();
59 
60  m_shapers.isRequired(m_SVDShaperDigitsName);
61  m_clusters.isRequired(m_SVDClustersName);
62  m_Intercepts.isRequired(m_InterceptsName);
63 
64  n_events = 0;
65 
66  m_rootFilePtr = new TFile(m_rootFileName.c_str(), "RECREATE");
67 
68 }
69 
70 void SVDChannelMappingModule::event()
71 {
72 
73  n_events++;
74 
75 
76 
77  for (auto& it : m_Intercepts)
78  fillSensorInterHistos(&it);
79 }
80 
81 void SVDChannelMappingModule::createHistosDictionaries()
82 {
83 
84  // VXD::GeoCache& aGeometry = VXD::GeoCache::getInstance();
85 
86  string name; //name of the histogram
87  string title; //title of the histogram
88  TH2F* tmp2D; //temporary 2D histo used to set axis title
89  TH1F* tmp1D; //temporary 1D histo used to set axis title
90 
91 
92  std::set<Belle2::VxdID> svdLayers = m_aGeometry.getLayers(VXD::SensorInfoBase::SVD);
93  std::set<Belle2::VxdID>::iterator itSvdLayers = svdLayers.begin();
94 
95  while (itSvdLayers != svdLayers.end()) {
96 
97  std::set<Belle2::VxdID> svdLadders = m_aGeometry.getLadders(*itSvdLayers);
98  std::set<Belle2::VxdID>::iterator itSvdLadders = svdLadders.begin();
99 
100  while (itSvdLadders != svdLadders.end()) {
101 
102  std::set<Belle2::VxdID> svdSensors = m_aGeometry.getSensors(*itSvdLadders);
103  std::set<Belle2::VxdID>::iterator itSvdSensors = svdSensors.begin();
104 
105  while (itSvdSensors != svdSensors.end()) {
106 
107  string sensorid = std::to_string(itSvdSensors->getLayerNumber()) + "_" + std::to_string(itSvdSensors->getLadderNumber()) + "_" +
108  std::to_string(itSvdSensors->getSensorNumber());
109 
110 
111  // ------ HISTOGRAMS WITH A FILL PER INTERCEPT -------
112 
113  // coor U and V
114  name = "hCoorU_" + sensorid;
115  title = "U coordinate of the extrapolation in U for sensor " + sensorid;
116  tmp1D = new TH1F(name.c_str(), title.c_str(), 100, -2.9, 2.9);
117  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
118  (
119  (Belle2::VxdID)*itSvdSensors,
121  tmp1D,
122  [](TH1 * hPtr, const SVDIntercept * inter) { hPtr->Fill(inter->getCoorU()); }
123  )
124  )
125  );
126  m_histoList_digits->Add(tmp1D);
127 
128  name = "hCoorV_" + sensorid;
129  title = "V coordinate of the extrapolation in V for sensor " + sensorid;
130  tmp1D = new TH1F(name.c_str(), title.c_str(), 100, -6.15, 6.15);
131  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
132  (
133  (Belle2::VxdID)*itSvdSensors,
135  tmp1D,
136  [](TH1 * hPtr, const SVDIntercept * inter) { hPtr->Fill(inter->getCoorV()); }
137  )
138  )
139  );
140  m_histoList_digits->Add(tmp1D);
141 
142  // Intercept U vs V coordinate
143  name = "hCoorU_vs_CoorV_" + sensorid;
144  title = "U vs V intercept (cm) " + sensorid;
145  tmp2D = new TH2F(name.c_str(), title.c_str(), 100, -2.9, 2.9, 100, -6.15, 6.15);
146  tmp2D->GetXaxis()->SetTitle("intercept U coor (cm)");
147  tmp2D->GetYaxis()->SetTitle("intercept V coor (cm)");
148  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
149  (
150  (Belle2::VxdID)*itSvdSensors,
152  tmp2D,
153  [](TH1 * hPtr, const SVDIntercept * inter) { hPtr->Fill(inter->getCoorU(), inter->getCoorV()); }
154  )
155  )
156  );
157  m_histoList_digits->Add(tmp2D);
158 
159  // sigma U and V
160  name = "hStatErrU_" + sensorid;
161  title = "stat error of the extrapolation in U for sensor " + sensorid;
162  tmp1D = new TH1F(name.c_str(), title.c_str(), 100, 0, 0.35);
163  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
164  (
165  (Belle2::VxdID)*itSvdSensors,
167  tmp1D,
168  [](TH1 * hPtr, const SVDIntercept * inter) { hPtr->Fill(inter->getSigmaU()); }
169  )
170  )
171  );
172  m_histoList_digits->Add(tmp1D);
173 
174  name = "hStatErrV_" + sensorid;
175  title = "stat error of the extrapolation in V for sensor " + sensorid;
176  tmp1D = new TH1F(name.c_str(), title.c_str(), 100, 0, 0.35);
177  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
178  (
179  (Belle2::VxdID)*itSvdSensors,
181  tmp1D,
182  [](TH1 * hPtr, const SVDIntercept * inter) { hPtr->Fill(inter->getSigmaV()); }
183  )
184  )
185  );
186  m_histoList_digits->Add(tmp1D);
187 
188  //---------------
189  //digits
190  //1D residuals
191  name = "hDigitResidU_" + sensorid;
192  title = "U residuals = intercept - digit, for sensor " + sensorid;
193  tmp1D = new TH1F(name.c_str(), title.c_str(), 1000, -2.9, 2.9);
194  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
195  (
196  (Belle2::VxdID)*itSvdSensors,
198  tmp1D,
199  [this](TH1 * hPtr, const SVDIntercept * inter) {
200  StoreArray<SVDShaperDigit> SVDShaperDigits(this->m_SVDShaperDigitsName);
201 
202  for (auto& it : SVDShaperDigits)
203  if ((int)it.getSensorID() == (int)inter->getSensorID()) {
204  if (it.isUStrip()) {
205  const VXD::SensorInfoBase& aSensorInfo = m_aGeometry.getSensorInfo(it.getSensorID());
206  hPtr->Fill(inter->getCoorU() - aSensorInfo.getUCellPosition(it.getCellID()));
207  }
208  }
209  }
210  )
211  )
212  );
213  m_histoList_digits->Add(tmp1D);
214 
215 
216  name = "hDigitResidV_" + sensorid;
217  title = "V residuals = intercept - digit, for sensor " + sensorid;
218  tmp1D = new TH1F(name.c_str(), title.c_str(), 1000, -6.15, 6.15);
219  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
220  (
221  (Belle2::VxdID)*itSvdSensors,
223  tmp1D,
224  [this](TH1 * hPtr, const SVDIntercept * inter) {
225  StoreArray<SVDShaperDigit> SVDShaperDigits(this->m_SVDShaperDigitsName);
226 
227  for (auto& it : SVDShaperDigits)
228  if ((int)it.getSensorID() == (int)inter->getSensorID()) {
229  if (!it.isUStrip()) {
230  const VXD::SensorInfoBase& aSensorInfo = m_aGeometry.getSensorInfo(it.getSensorID());
231  hPtr->Fill(inter->getCoorV() - aSensorInfo.getVCellPosition(it.getCellID()));
232  }
233  }
234  }
235  )
236  )
237  );
238  m_histoList_digits->Add(tmp1D);
239 
240  //residual U,V vs coordinate U,V
241  name = "hDigitResidU_vs_DigitU_" + sensorid;
242  title = "U residual (cm) vs digit U (cm) " + sensorid;
243  tmp2D = new TH2F(name.c_str(), title.c_str(), 1000, -2.9, 2.9, 1000, -2.9, 2.9);
244  tmp2D->GetYaxis()->SetTitle("U resid (cm)");
245  tmp2D->GetXaxis()->SetTitle("U digit (cm)");
246  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
247  (
248  (Belle2::VxdID)*itSvdSensors,
250  tmp2D,
251  [this](TH1 * hPtr, const SVDIntercept * inter) {
252  StoreArray<SVDShaperDigit> SVDShaperDigits(this->m_SVDShaperDigitsName);
253 
254  for (auto& it : SVDShaperDigits)
255  if (((int)it.getSensorID() == (int)inter->getSensorID()) && it.isUStrip()) {
256  const VXD::SensorInfoBase& aSensorInfo = m_aGeometry.getSensorInfo(it.getSensorID());
257  double resid = inter->getCoorU() - aSensorInfo.getUCellPosition(it.getCellID());
258  hPtr->Fill(aSensorInfo.getUCellPosition(it.getCellID()), resid);
259  }
260  }
261  )
262  )
263  );
264  m_histoList_digits->Add(tmp2D);
265 
266  name = "hDigitResidV_vs_DigitV_" + sensorid;
267  title = "V residual (cm) vs digit V (cm) " + sensorid;
268  tmp2D = new TH2F(name.c_str(), title.c_str(), 1000, -6.15, 6.15, 1000, -6.15, 6.15);
269  tmp2D->GetYaxis()->SetTitle("V resid (cm)");
270  tmp2D->GetXaxis()->SetTitle("V digit (cm)");
271  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
272  (
273  (Belle2::VxdID)*itSvdSensors,
275  tmp2D,
276  [this](TH1 * hPtr, const SVDIntercept * inter) {
277  StoreArray<SVDShaperDigit> SVDShaperDigits(this->m_SVDShaperDigitsName);
278 
279  for (auto& it : SVDShaperDigits)
280  if (((int)it.getSensorID() == (int)inter->getSensorID()) && (!it.isUStrip())) {
281  const VXD::SensorInfoBase& aSensorInfo = m_aGeometry.getSensorInfo(it.getSensorID());
282  double resid = inter->getCoorV() - aSensorInfo.getVCellPosition(it.getCellID());
283  hPtr->Fill(aSensorInfo.getVCellPosition(it.getCellID()), resid);
284  }
285  }
286  )
287  )
288  );
289  m_histoList_digits->Add(tmp2D);
290 
291 
292  // scatter plot: U,V intercept in cm VS U,V cell position
293  name = "hCoorU_vs_DigitU_" + sensorid;
294  title = "U intercept (cm) vs U digit (cm) " + sensorid;
295  tmp2D = new TH2F(name.c_str(), title.c_str(), 1000, -2.9, 2.9, 1000, -2.9, 2.9);
296  tmp2D->GetXaxis()->SetTitle("intercept U coor (cm)");
297  tmp2D->GetYaxis()->SetTitle("digit U coor (cm)");
298  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
299  (
300  (Belle2::VxdID)*itSvdSensors,
302  tmp2D,
303  [this](TH1 * hPtr, const SVDIntercept * inter) {
304  StoreArray<SVDShaperDigit> SVDShaperDigits(this->m_SVDShaperDigitsName);
305 
306  for (auto& it : SVDShaperDigits)
307  if (((int)it.getSensorID() == (int)inter->getSensorID()) && (it.isUStrip())) {
308  const VXD::SensorInfoBase& aSensorInfo = m_aGeometry.getSensorInfo(it.getSensorID());
309  hPtr->Fill(inter->getCoorU(), aSensorInfo.getUCellPosition(it.getCellID()));
310  }
311  }
312  )
313  )
314  );
315  m_histoList_digits->Add(tmp2D);
316 
317  name = "hCoorV_vs_DigitV_" + sensorid;
318  title = "V intercept (cm) vs V digit (ID) " + sensorid;
319  tmp2D = new TH2F(name.c_str(), title.c_str(), 1000, -6.15, 6.15, 1000, -6.15, 6.15);
320  tmp2D->GetXaxis()->SetTitle("intercept V coor (cm)");
321  tmp2D->GetYaxis()->SetTitle("digi V coor (cm)");
322  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
323  (
324  (Belle2::VxdID)*itSvdSensors,
326  tmp2D,
327  [this](TH1 * hPtr, const SVDIntercept * inter) {
328  StoreArray<SVDShaperDigit> SVDShaperDigits(this->m_SVDShaperDigitsName);
329 
330  for (auto& it : SVDShaperDigits) {
331  if (((int)it.getSensorID() == (int)inter->getSensorID()) && (!it.isUStrip())) {
332  const VXD::SensorInfoBase& aSensorInfo = m_aGeometry.getSensorInfo(it.getSensorID());
333  hPtr->Fill(inter->getCoorV(), aSensorInfo.getVCellPosition(it.getCellID()));
334  }
335  }
336  }
337  )
338  )
339  );
340  m_histoList_digits->Add(tmp2D);
341 
342  //-----------------
343  //clusters
344  //1D residuals
345  name = "hClusterResidU_" + sensorid;
346  title = "U residuals = intercept - cluster, for sensor " + sensorid;
347  tmp1D = new TH1F(name.c_str(), title.c_str(), 1000, -2.9, 2.9);
348  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
349  (
350  (Belle2::VxdID)*itSvdSensors,
352  tmp1D,
353  [this](TH1 * hPtr, const SVDIntercept * inter) {
354  StoreArray<SVDCluster> SVDClusters(this->m_SVDClustersName);
355 
356  for (auto& it : SVDClusters)
357  if ((int)it.getSensorID() == (int)inter->getSensorID()) {
358  if (it.isUCluster()) {
359  hPtr->Fill(inter->getCoorU() - it.getPosition());
360  }
361  }
362  }
363  )
364  )
365  );
366  m_histoList_clusters->Add(tmp1D);
367 
368 
369  name = "hClusterResidV_" + sensorid;
370  title = "V residuals = intercept - cluster, for sensor " + sensorid;
371  tmp1D = new TH1F(name.c_str(), title.c_str(), 1000, -6.15, 6.15);
372  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
373  (
374  (Belle2::VxdID)*itSvdSensors,
376  tmp1D,
377  [this](TH1 * hPtr, const SVDIntercept * inter) {
378  StoreArray<SVDCluster> SVDClusters(this->m_SVDClustersName);
379 
380  for (auto& it : SVDClusters)
381  if ((int)it.getSensorID() == (int)inter->getSensorID()) {
382  if (!it.isUCluster()) {
383  hPtr->Fill(inter->getCoorV() - it.getPosition());
384  }
385  }
386  }
387  )
388  )
389  );
390  m_histoList_clusters->Add(tmp1D);
391 
392  //residual U,V vs coordinate U,V
393  name = "hClusterResidU_vs_ClusterU_" + sensorid;
394  title = "U residual (cm) vs cluster U (cm) " + sensorid;
395  tmp2D = new TH2F(name.c_str(), title.c_str(), 1000, -2.9, 2.9, 1000, -2.9, 2.9);
396  tmp2D->GetYaxis()->SetTitle("U resid (cm)");
397  tmp2D->GetXaxis()->SetTitle("U cluster (cm)");
398  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
399  (
400  (Belle2::VxdID)*itSvdSensors,
402  tmp2D,
403  [this](TH1 * hPtr, const SVDIntercept * inter) {
404  StoreArray<SVDCluster> SVDClusters(this->m_SVDClustersName);
405 
406  for (auto& it : SVDClusters)
407  if (((int)it.getSensorID() == (int)inter->getSensorID()) && it.isUCluster()) {
408  double resid = inter->getCoorU() - it.getPosition();
409  hPtr->Fill(it.getPosition(), resid);
410  }
411  }
412  )
413  )
414  );
415  m_histoList_clusters->Add(tmp2D);
416 
417  name = "hClusterResidV_vs_ClusterV_" + sensorid;
418  title = "V residual (cm) vs cluster V (cm) " + sensorid;
419  tmp2D = new TH2F(name.c_str(), title.c_str(), 1000, -6.15, 6.15, 1000, -6.15, 6.15);
420  tmp2D->GetYaxis()->SetTitle("V resid (cm)");
421  tmp2D->GetXaxis()->SetTitle("V cluster (cm)");
422  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
423  (
424  (Belle2::VxdID)*itSvdSensors,
426  tmp2D,
427  [this](TH1 * hPtr, const SVDIntercept * inter) {
428  StoreArray<SVDCluster> SVDClusters(this->m_SVDClustersName);
429 
430  for (auto& it : SVDClusters)
431  if (((int)it.getSensorID() == (int)inter->getSensorID()) && (!it.isUCluster())) {
432  double resid = inter->getCoorV() - it.getPosition();
433  hPtr->Fill(it.getPosition(), resid);
434  }
435  }
436  )
437  )
438  );
439  m_histoList_clusters->Add(tmp2D);
440 
441 
442  // scatter plot: U,V intercept in cm VS U,V cell position
443  name = "hCoorU_vs_ClusterU_" + sensorid;
444  title = "U intercept (cm) vs U cluster (cm) " + sensorid;
445  tmp2D = new TH2F(name.c_str(), title.c_str(), 1000, -2.9, 2.9, 1000, -2.9, 2.9);
446  tmp2D->GetXaxis()->SetTitle("intercept U coor (cm)");
447  tmp2D->GetYaxis()->SetTitle("cluster U coor (cm)");
448  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
449  (
450  (Belle2::VxdID)*itSvdSensors,
452  tmp2D,
453  [this](TH1 * hPtr, const SVDIntercept * inter) {
454  StoreArray<SVDCluster> SVDClusters(this->m_SVDClustersName);
455 
456  for (auto& it : SVDClusters)
457  if (((int)it.getSensorID() == (int)inter->getSensorID()) && (it.isUCluster())) {
458  hPtr->Fill(inter->getCoorU(), it.getPosition());
459  }
460  }
461  )
462  )
463  );
464  m_histoList_clusters->Add(tmp2D);
465 
466  name = "hCoorV_vs_ClusterV_" + sensorid;
467  title = "V intercept (cm) vs V cluster (ID) " + sensorid;
468  tmp2D = new TH2F(name.c_str(), title.c_str(), 1000, -6.15, 6.15, 1000, -6.15, 6.15);
469  tmp2D->GetXaxis()->SetTitle("intercept V coor (cm)");
470  tmp2D->GetYaxis()->SetTitle("cluster V coor (cm)");
471  hInterDictionary.insert(pair< Belle2::VxdID, InterHistoAndFill >
472  (
473  (Belle2::VxdID)*itSvdSensors,
475  tmp2D,
476  [this](TH1 * hPtr, const SVDIntercept * inter) {
477  StoreArray<SVDCluster> SVDClusters(this->m_SVDClustersName);
478 
479  for (auto& it : SVDClusters) {
480  if (((int)it.getSensorID() == (int)inter->getSensorID()) && (!it.isUCluster())) {
481  hPtr->Fill(inter->getCoorV(), it.getPosition());
482  }
483  }
484  }
485  )
486  )
487  );
488  m_histoList_clusters->Add(tmp2D);
489 
490  ++itSvdSensors;
491  }
492  ++itSvdLadders;
493  }
494  ++itSvdLayers;
495  }
496 
497 }
498 
499 void SVDChannelMappingModule::fillSensorInterHistos(const SVDIntercept* inter)
500 {
501 
502  auto its = hInterDictionary.equal_range(inter->getSensorID());
503 
504  for (auto it = its.first; it != its.second; ++it) {
505  InterHistoAndFill aInterHistoAndFill = it->second;
506  aInterHistoAndFill.second(aInterHistoAndFill.first, inter);
507  }
508 
509 }
510 
511 void SVDChannelMappingModule::terminate()
512 {
513 
514  if (m_rootFilePtr != NULL) {
515  m_rootFilePtr->cd();
516 
517 
518  TDirectory* oldDir = gDirectory;
519  TObject* obj;
520 
521  TDirectory* dir_digits = oldDir->mkdir("digits");
522  dir_digits->cd();
523  TIter nextH_digits(m_histoList_digits);
524  while ((obj = nextH_digits()))
525  obj->Write();
526 
527  TDirectory* dir_clusters = oldDir->mkdir("clusters");
528  dir_clusters->cd();
529  TIter nextH_clusters(m_histoList_clusters);
530  while ((obj = nextH_clusters()))
531  obj->Write();
532 
533  }
534 }
Belle2::SVDIntercept
SVDIntercept stores the U,V coordinates and uncertainties of the intersection of a track with an SVD ...
Definition: SVDIntercept.h:32
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::VXD::SensorInfoBase
Base class to provide Sensor Information for PXD and SVD.
Definition: SensorInfoBase.h:40
Belle2::SVDChannelMappingModule::InterHistoAndFill
std::pair< TH1 *, std::function< void(TH1 *, const SVDIntercept *) > > InterHistoAndFill
typedef: histograms to be filled once per intercept + filling function
Definition: SVDChannelMappingModule.h:84
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::VXD::SensorInfoBase::getVCellPosition
double getVCellPosition(int vID) const
Return the position of a specific strip/pixel in v direction.
Definition: SensorInfoBase.h:191
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::SVDChannelMappingModule
The Channel Mapping Check Module.
Definition: SVDChannelMappingModule.h:49
Belle2::VXD::SensorInfoBase::getUCellPosition
double getUCellPosition(int uID, int vID=-1) const
Return the position of a specific strip/pixel in u direction.
Definition: SensorInfoBase.h:179
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33