Belle II Software  release-05-02-19
SpacePointHelperFunctions.h
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Jakob Lettenbichler, Giulia Casarosa, James Webb *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #pragma once
12 
13 #include <vector>
14 
15 #include <svd/calibration/SVDClusterCalibrations.h>
16 #include <svd/dataobjects/SVDCluster.h>
17 
18 #include <framework/datastore/StoreArray.h>
19 
20 #include <vxd/dataobjects/VxdID.h>
21 
22 #include <unordered_map>
23 
24 #include <TH2.h>
25 #include <math.h>
26 #include <TFile.h>
27 
28 namespace Belle2 {
39  struct ClustersOnSensor {
40 
41  public:
42 
44  inline void addCluster(const SVDCluster* entry)
45  {
46  vxdID = entry->getSensorID();
47  if (entry->isUCluster() == true) { clustersU.push_back(entry); return; }
48  clustersV.push_back(entry);
49  }
50 
52  VxdID vxdID;
53 
58  std::vector<const SVDCluster*> clustersU;
59 
64  std::vector<const SVDCluster*> clustersV;
65 
66  };
67 
73  template <class SpacePointType> void provideSVDClusterSingles(const StoreArray<SVDCluster>& svdClusters,
74  StoreArray<SpacePointType>& spacePoints)
75  {
76  for (unsigned int i = 0; i < uint(svdClusters.getEntries()); ++i) {
77  const SVDCluster* currentCluster = svdClusters[i];
78  std::vector<const SVDCluster*> currentClusterCombi = { currentCluster };
79  SpacePointType* newSP = spacePoints.appendNew(currentClusterCombi);
80  newSP->addRelationTo(currentCluster);
81  }
82  }
83 
84 
85 
94  inline void findPossibleCombinations(const Belle2::ClustersOnSensor& aSensor,
95  std::vector< std::vector<const SVDCluster*> >& foundCombinations, SVDClusterCalibrations& clusterCal)
96  {
97 
98  for (const SVDCluster* uCluster : aSensor.clustersU) {
99  if (! clusterCal.isClusterInTime(uCluster->getSensorID(), 1, uCluster->getClsTime())) {
100  B2DEBUG(1, "Cluster rejected due to timing cut. Cluster time: " << uCluster->getClsTime());
101  continue;
102  }
103  for (const SVDCluster* vCluster : aSensor.clustersV) {
104  if (! clusterCal.isClusterInTime(vCluster->getSensorID(), 0, vCluster->getClsTime())) {
105  B2DEBUG(1, "Cluster rejected due to timing cut. Cluster time: " << vCluster->getClsTime());
106  continue;
107  }
108 
109  if (! clusterCal.areClusterTimesCompatible(vCluster->getSensorID(), uCluster->getClsTime(), vCluster->getClsTime())) {
110  B2DEBUG(1, "Cluster combination rejected due to timing cut. Cluster time U (" << uCluster->getClsTime() <<
111  ") is incompatible with Cluster time V (" << vCluster->getClsTime() << ")");
112  continue;
113  }
114 
115 
116  foundCombinations.push_back({uCluster, vCluster});
117 
118 
119  }
120  }
121 
122 
123 
124 
125  }
126 
132  inline void spPDFName(const VxdID& sensor, int uSize, int vSize, int maxClusterSize, std::string& PDFName,
133  std::string& errorPDFName, bool useLegacyNaming)
134  {
135  if (useLegacyNaming == true) {
136 
137  if (uSize > maxClusterSize) uSize = maxClusterSize;
138  if (vSize > maxClusterSize) vSize = maxClusterSize;
139 
140  std::string sensorName;
141 
142  if (sensor.getLayerNumber() == 3) sensorName = "l3";
143  if (sensor.getLayerNumber() > 3 && sensor.getSensorNumber() == 1) sensorName = "trap";
144  if (sensor.getLayerNumber() > 3 && sensor.getSensorNumber() > 1) sensorName = "large";
145 
146  PDFName = sensorName + std::to_string(uSize) + std::to_string(vSize);
147  errorPDFName = "error" + PDFName;
148  } else {
149 
150  if (uSize > maxClusterSize) uSize = maxClusterSize;
151  if (vSize > maxClusterSize) vSize = maxClusterSize;
152 
153  int layer = sensor.getLayerNumber();
154  int ladder = sensor.getLadderNumber();
155  int sens = sensor.getSensorNumber();
156 
157  PDFName = std::to_string(layer) + "." + std::to_string(ladder) + "." + std::to_string(sens) + "." + std::to_string(
158  uSize) + "." + std::to_string(vSize);
159  errorPDFName = PDFName + "_Error";
160  }
161 
162 
163 
164  }
165 
166 
174  inline void calculatePairingProb(TFile* pdfFile, std::vector<const SVDCluster*>& clusters, double& prob, double& error,
175  bool useLegacyNaming)
176  {
177 
178  int maxSize;
179  int pdfEntries = pdfFile->GetListOfKeys()->GetSize();
180  if (useLegacyNaming == true) {
181  maxSize = floor(sqrt((pdfEntries - 4) / 6)); //4(time+size)+3(sensors)*2(prob/error)*size^2(u/v combo.)
182  } else {
183  maxSize = floor(sqrt((pdfEntries - 4) / 344)); //4(time+size)+172(sensorType)*2(prob/error)*size^2(u/v combo.)
184  }
185  std::string chargeProbInput;
186  std::string chargeErrorInput;
187 
188  spPDFName(clusters[0]->getSensorID(), clusters[0]->getSize(), clusters[1]->getSize(), maxSize,
189  chargeProbInput, chargeErrorInput, useLegacyNaming);
190  std::string timeProbInput = "timeProb";
191  std::string timeErrorInput = "timeError";
192  std::string sizeProbInput = "sizeProb";
193  std::string sizeErrorInput = "sizeError";
194 
195 
196  TH2F* chargePDF = nullptr;
197  TH2F* chargeError = nullptr;
198  TH2F* timePDF = nullptr;
199  TH2F* timeError = nullptr;
200  TH2F* sizePDF = nullptr;
201  TH2F* sizeError = nullptr;
202 
203  pdfFile->GetObject(chargeProbInput.c_str(), chargePDF);
204  pdfFile->GetObject(chargeErrorInput.c_str(), chargeError);
205  pdfFile->GetObject(timeProbInput.c_str(), timePDF);
206  pdfFile->GetObject(timeErrorInput.c_str(), timeError);
207  pdfFile->GetObject(sizeProbInput.c_str(), sizePDF);
208  pdfFile->GetObject(sizeErrorInput.c_str(), sizeError);
209 
210  int xChargeBin = chargePDF->GetXaxis()->FindFixBin(clusters[0]->getCharge());
211  int yChargeBin = chargePDF->GetYaxis()->FindFixBin(clusters[1]->getCharge());
212 
213  int xTimeBin = timePDF->GetXaxis()->FindFixBin(clusters[0]->getClsTime());
214  int yTimeBin = timePDF->GetYaxis()->FindFixBin(clusters[1]->getClsTime());
215 
216 
217  int xSizeBin = sizePDF->GetXaxis()->FindFixBin(clusters[0]->getSize());
218  int ySizeBin = sizePDF->GetYaxis()->FindFixBin(clusters[1]->getSize());
219 
220  double chargeProb = chargePDF->GetBinContent(xChargeBin, yChargeBin);
221  double timeProb = timePDF->GetBinContent(xTimeBin, yTimeBin);
222  double sizeProb = sizePDF->GetBinContent(xSizeBin, ySizeBin);
223  double chargeProbError = chargePDF->GetBinContent(xChargeBin, yChargeBin);
224  double timeProbError = timePDF->GetBinContent(xTimeBin, yTimeBin);
225  double sizeProbError = sizePDF->GetBinContent(xSizeBin, ySizeBin);
226 
227 
228  if (chargeProbError == 0) {
229  B2DEBUG(1, "svdClusterProbabilityEstimator has not been run, spacePoint QI will return zero!");
230  }
231 
232  prob = chargeProb * timeProb * sizeProb * clusters[0]->getQuality() * clusters[1]->getQuality();
233  error = prob * sqrt(pow(timeProb * sizeProb * clusters[0]->getQuality() * clusters[1]->getQuality() * chargeProbError , 2) +
234  pow(chargeProb * sizeProb * clusters[0]->getQuality() * clusters[1]->getQuality() * timeProbError , 2) +
235  pow(chargeProb * timeProb * clusters[0]->getQuality() * clusters[1]->getQuality() * sizeProbError , 2) +
236  pow(chargeProb * timeProb * sizeProb * clusters[1]->getQuality() * clusters[0]->getQualityError() , 2) +
237  pow(chargeProb * timeProb * sizeProb * clusters[0]->getQuality() * clusters[1]->getQualityError() , 2));
238  }
239 
247  template <class SpacePointType> void provideSVDClusterCombinations(const StoreArray<SVDCluster>& svdClusters,
248  StoreArray<SpacePointType>& spacePoints, SVDClusterCalibrations& clusterCal, bool useQualityEstimator, TFile* pdfFile,
249  bool useLegacyNaming)
250  {
251  std::unordered_map<VxdID::baseType, ClustersOnSensor>
252  activatedSensors; // collects one entry per sensor, each entry will contain all Clusters on it TODO: better to use a sorted vector/list?
253  std::vector<std::vector<const SVDCluster*> >
254  foundCombinations; // collects all combinations of Clusters which were possible (condition: 1u+1v-Cluster on the same sensor)
255 
256  double probability;
257  double error;
258 
259  // sort Clusters by sensor. After the loop, each entry of activatedSensors contains all U and V-type clusters on that sensor
260  for (unsigned int i = 0; i < uint(svdClusters.getEntries()); ++i) {
261  SVDCluster* currentCluster = svdClusters[i];
262 
263  activatedSensors[currentCluster->getSensorID().getID()].addCluster(currentCluster);
264  }
265 
266 
267  for (auto& aSensor : activatedSensors)
268  findPossibleCombinations(aSensor.second, foundCombinations, clusterCal);
269 
270 
271  for (auto& clusterCombi : foundCombinations) {
272  SpacePointType* newSP = spacePoints.appendNew(clusterCombi);
273  if (useQualityEstimator == true) {
274  calculatePairingProb(pdfFile, clusterCombi, probability, error, useLegacyNaming);
275  newSP->setQualityEstimation(probability);
276  newSP->setQualityEstimationError(error);
277  }
278  for (auto* cluster : clusterCombi) {
279  newSP->addRelationTo(cluster, cluster->isUCluster() ? 1. : -1.);
280  }
281  }
282  }
283 
284 
286 } //Belle2 namespace
Belle2::StoreArray::appendNew
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:256
Belle2::provideSVDClusterCombinations
void provideSVDClusterCombinations(const StoreArray< SVDCluster > &svdClusters, StoreArray< SpacePointType > &spacePoints, SVDClusterCalibrations &clusterCal, bool useQualityEstimator, TFile *pdfFile, bool useLegacyNaming)
finds all possible combinations of U and V Clusters for SVDClusters.
Definition: SpacePointHelperFunctions.h:255
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
Belle2::findPossibleCombinations
void findPossibleCombinations(const Belle2::ClustersOnSensor &aSensor, std::vector< std::vector< const SVDCluster * > > &foundCombinations, SVDClusterCalibrations &clusterCal)
stores all possible 2-Cluster-combinations.
Definition: SpacePointHelperFunctions.h:102
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::ClustersOnSensor::clustersU
std::vector< const SVDCluster * > clustersU
stores all SVDclusters of U type.
Definition: SpacePointHelperFunctions.h:66
Belle2::ClustersOnSensor::addCluster
void addCluster(const SVDCluster *entry)
member function to automatically add the cluster to its corresponding entry
Definition: SpacePointHelperFunctions.h:52
Belle2::calculatePairingProb
void calculatePairingProb(TFile *pdfFile, std::vector< const SVDCluster * > &clusters, double &prob, double &error, bool useLegacyNaming)
Function to extract probability of correct (pair from signal hit) cluster pairing from preconfigured ...
Definition: SpacePointHelperFunctions.h:182
Belle2::VxdID::getID
baseType getID() const
Get the unique id.
Definition: VxdID.h:104
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::SVDCluster::getSensorID
VxdID getSensorID() const
Get the sensor ID.
Definition: SVDCluster.h:110
Belle2::provideSVDClusterSingles
void provideSVDClusterSingles(const StoreArray< SVDCluster > &svdClusters, StoreArray< SpacePointType > &spacePoints)
simply store one spacePoint for each existing SVDCluster.
Definition: SpacePointHelperFunctions.h:81
Belle2::SVDCluster
The SVD Cluster class This class stores all information about reconstructed SVD clusters.
Definition: SVDCluster.h:38
Belle2::StoreArray< SVDCluster >
Belle2::ClustersOnSensor::vxdID
VxdID vxdID
Id of sensor, TODO can be removed if struct is used in a map.
Definition: SpacePointHelperFunctions.h:60
Belle2::ClustersOnSensor
small struct for storing all clusters of the same sensor in one container.
Definition: SpacePointHelperFunctions.h:47
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226
Belle2::ClustersOnSensor::clustersV
std::vector< const SVDCluster * > clustersV
stores all SVDclusters of V type.
Definition: SpacePointHelperFunctions.h:72