Belle II Software  release-05-01-25
SVDSpacePointQICalibrationModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2018 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: James Webb *
7  * *
8  **************************************************************************/
9 
10 #include <svd/modules/svdSpacePointCreator/SpacePointHelperFunctions.h>
11 #include <svd/modules/svdSpacePointCreator/SVDSpacePointQICalibrationModule.h>
12 //#include <svd/modules/svdSpacePointCreator/SVDSpacePointCreatorModule.h>
13 
14 #include <math.h>
15 
16 using namespace std;
17 using namespace Belle2;
18 
19 
20 REG_MODULE(SVDSpacePointQICalibration)
21 
23  Module()
24 {
25 
26  setDescription("Generate PDFs used in assigning quality to SpacePoints.");
27  setPropertyFlags(c_ParallelProcessingCertified);
28 
29  // 1. Collections.
30  addParam("SVDClusters", m_svdClustersName,
31  "SVDCluster collection name", string(""));
32 
33  addParam("RecoTracks", m_recoTracksName,
34  "RecoTracks collection name", string(""));
35 
36  // 2.Modification parameters:
37  addParam("NameOfInstance", m_nameOfInstance,
38  "allows the user to set an identifier for this module. Usefull if one wants to use several instances of that module", string(""));
39  addParam("binSize", m_binSize, "Number of bins in charge distribution.",
40  int(50));
41 
42  addParam("maxClusterSize", m_maxClusterSize, "Max number of strips the PDF are separated into.",
43  int(5));
44 
45  addParam("useLegacyNaming", m_useLegacyNaming,
46  "Use legacy pdf naming", bool(true));
47 
48  addParam("outputFileName", m_outputFileName,
49  "Name of output file containing pdfs", std::string("spacePointQICalibration.root"));
50 }
51 
52 
53 
54 void SVDSpacePointQICalibrationModule::initialize()
55 {
56  // prepare all store- and relationArrays:
57  m_svdClusters.isRequired(m_svdClustersName);
58  m_recoTracks.isRequired(m_recoTracksName);
59 
60 
61 
62  VXD::GeoCache& geo = VXD::GeoCache::getInstance();
63  for (auto& layers : geo.getLayers(VXD::SensorInfoBase::SVD)) {
64  for (auto& ladders : geo.getLadders(layers)) {
65  for (auto& sensors : geo.getSensors(ladders)) {
66 
67  for (int uSize = 1; uSize <= m_maxClusterSize; uSize++) {
68  for (int vSize = 1; vSize <= m_maxClusterSize; vSize++) {
69  std::string sensorName = "";
70  std::string errorStringName = "";
71 
72 
73  spPDFName(sensors, uSize, vSize, m_maxClusterSize, sensorName, errorStringName, m_useLegacyNaming);
74 
75  std::string signalHistName = sensorName + "signal";
76  std::string backgroundHistName = sensorName + "background";
77 
78  if (signalHistMap.count(sensorName) == 0) {
79  TH2F* sigHist = new TH2F(signalHistName.c_str(), "", m_binSize, 0, 200000, m_binSize, 0, 200000);
80  TH2F* bkgHist = new TH2F(backgroundHistName.c_str(), "", m_binSize, 0, 200000, m_binSize, 0, 200000);
81 
82  signalHistMap[sensorName] = sigHist;
83  backgroundHistMap[sensorName] = bkgHist;
84  }
85  }
86  }
87  }
88  }
89  }
90  TH2F* timeSignal = new TH2F("timeSignal", "", 40, -100, 100, 40, -100, 100);
91  TH2F* timeBackground = new TH2F("timeBackground", "", 40, -100, 100, 40, -100, 100);
92  TH2F* sizeSignal = new TH2F("sizeSignal", "", 20, 0, 20, 20, 0, 20);
93  TH2F* sizeBackground = new TH2F("sizeBackground", "", 20, 0, 20, 20, 0, 20);
94 
95  signalHistMap["timeSignal"] = timeSignal;
96  signalHistMap["sizeSignal"] = sizeSignal;
97  backgroundHistMap["timeBackground"] = timeBackground;
98  backgroundHistMap["sizeBackground"] = sizeBackground;
99 
100 }
101 
102 
103 
104 void SVDSpacePointQICalibrationModule::event()
105 {
106  for (auto& track : m_recoTracks) {
107  if (track.wasFitSuccessful() == 1) {
108  if (track.hasSVDHits() == 1) {
109  for (auto& svdHit : track.getSVDHitList()) {
110  if (svdHit->isUCluster()) {
111  for (auto& svdHit2 : track.getSVDHitList()) {
112  if (svdHit2->isUCluster() == 0) {
113  if (svdHit->getSensorID() == svdHit2->getSensorID()) {
114  std::string pdfName;
115  std::string errorStringName;
116  spPDFName(svdHit->getSensorID(), svdHit->getSize(), svdHit2->getSize(), m_maxClusterSize, pdfName, errorStringName,
117  m_useLegacyNaming);
118  auto sigHist = signalHistMap.at(pdfName);
119  sigHist->Fill(svdHit->getCharge(), svdHit2->getCharge());
120  auto timeSigHist = signalHistMap.at("timeSignal");
121  timeSigHist->Fill(svdHit->getClsTime(), svdHit2->getClsTime());
122  auto sizeSigHist = signalHistMap.at("sizeSignal");
123  sizeSigHist->Fill(svdHit->getSize(), svdHit2->getSize());
124  }
125  }
126  }
127  }
128  }
129  }
130  }
131  }
132 
133 
134  for (auto& uCluster : m_svdClusters) {
135  if (uCluster.isUCluster() == 1) {
136  for (auto& vCluster : m_svdClusters) {
137  if (vCluster.isUCluster() == 0) {
138  if (uCluster.getSensorID() == vCluster.getSensorID()) {
139  std::string pdfName;
140  std::string errorStringName;
141  spPDFName(uCluster.getSensorID(), uCluster.getSize(), vCluster.getSize(), m_maxClusterSize, pdfName, errorStringName,
142  m_useLegacyNaming);
143  auto bkgHist = backgroundHistMap.at(pdfName);
144  bkgHist->Fill(uCluster.getCharge(), vCluster.getCharge());
145  auto timeBkgHist = backgroundHistMap.at("timeBackground");
146  timeBkgHist->Fill(uCluster.getClsTime(), vCluster.getClsTime());
147  auto sizeBkgHist = backgroundHistMap.at("sizeBackground");
148  sizeBkgHist->Fill(uCluster.getSize(), vCluster.getSize());
149  }
150  }
151  }
152  }
153  }
154 }
155 
156 
157 
158 void SVDSpacePointQICalibrationModule::terminate()
159 {
160  //Open rootfile
161  TFile* f = new TFile(m_outputFileName.c_str(), "RECREATE");
162  f->cd();
163  std::vector<std::string> usedSensors;
164  VXD::GeoCache& geo = VXD::GeoCache::getInstance();
165  for (auto& layers : geo.getLayers(VXD::SensorInfoBase::SVD)) {
166  for (auto& ladders : geo.getLadders(layers)) {
167  for (auto& sensors : geo.getSensors(ladders)) {
168  for (int uSize = 1; uSize <= m_maxClusterSize; uSize++) {
169  for (int vSize = 1; vSize <= m_maxClusterSize; vSize++) {
170  std::string sensorName;
171  std::string errorName;
172 
173 
174  spPDFName(sensors, uSize, vSize, m_maxClusterSize, sensorName, errorName, m_useLegacyNaming);
175 
176  if (std::find(usedSensors.begin(), usedSensors.end(), sensorName.c_str()) == usedSensors.end()) {
177  usedSensors.push_back(sensorName);
178  TH2F* probHist = new TH2F(sensorName.c_str(), "", m_binSize, 0, 200000, m_binSize, 0, 200000);
179  TH2F* errorHist = new TH2F(errorName.c_str(), "", m_binSize, 0, 200000, m_binSize, 0, 200000);
180  calculateProb(signalHistMap[sensorName], backgroundHistMap[sensorName], probHist);
181  calculateError(signalHistMap[sensorName], backgroundHistMap[sensorName], errorHist);
182 
183  probHist->Write();
184  errorHist->Write();
185  }
186  }
187  }
188  }
189  }
190  }
191  TH2F* timeProb = new TH2F("timeProb", "", 40, -100, 100, 40, -100, 100);
192  TH2F* timeError = new TH2F("timeError", "", 40, -100, 100, 40, -100, 100);
193  TH2F* sizeProb = new TH2F("sizeProb", "", 20, 0, 20, 20, 0, 20);
194  TH2F* sizeError = new TH2F("sizeError", "", 20, 0, 20, 20, 0, 20);
195  calculateProb(signalHistMap["timeSignal"], backgroundHistMap["timeBackground"], timeProb);
196  calculateError(signalHistMap["timeSignal"], backgroundHistMap["timeBackground"], timeError);
197  calculateProb(signalHistMap["sizeSignal"], backgroundHistMap["sizeBackground"], sizeProb);
198  calculateError(signalHistMap["sizeSignal"], backgroundHistMap["sizeBackground"], sizeError);
199  timeProb->Write();
200  timeError->Write();
201  sizeProb->Write();
202  sizeError->Write();
203 
204  f->Close();
205 }
206 
207 
208 void SVDSpacePointQICalibrationModule::calculateProb(TH2F* signal, TH2F* background, TH2F* probability)
209 {
210  signal->Divide(background);
211  probability->Add(signal);
212 
213 }
214 void SVDSpacePointQICalibrationModule::calculateError(TH2F* signal, TH2F* background, TH2F* error)
215 {
216  int imax = signal->GetXaxis()->GetNbins();
217  int jmax = signal->GetYaxis()->GetNbins();
218  for (int i = 1; i <= imax; i++) {
219  for (int j = 1; j <= jmax; j++) {
220  int bkg = background->GetBinContent(i, j);
221  int sig = signal->GetBinContent(i, j);
222  double var = ((sig + 1) * (sig + 2)) / ((bkg + 2) * (bkg + 3)) -
223  ((sig + 1) * (sig + 1)) / ((bkg + 2) * (bkg + 2));
224  double err = sqrt(var);
225  error->SetBinContent(i, j, err);
226  }
227  }
228 }
Belle2::spPDFName
void spPDFName(const VxdID &sensor, int uSize, int vSize, int maxClusterSize, std::string &PDFName, std::string &errorPDFName, bool useLegacyNaming)
Function to set name of PDF for spacePoint quality estimation.
Definition: SpacePointHelperFunctions.h:140
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::VXD::GeoCache::getLayers
const std::set< Belle2::VxdID > getLayers(SensorInfoBase::SensorType sensortype=SensorInfoBase::VXD)
Return a set of all known Layers.
Definition: GeoCache.cc:177
Belle2::VXD::GeoCache::getSensors
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:205
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::SVDSpacePointQICalibrationModule
Imports Clusters of the SVD detector and converts them to spacePoints.
Definition: SVDSpacePointQICalibrationModule.h:41
Belle2::VXD::GeoCache
Class to faciliate easy access to sensor information of the VXD like coordinate transformations or pi...
Definition: GeoCache.h:41
Belle2::VXD::GeoCache::getLadders
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:194