Belle II Software  release-05-01-25
PXDValidationAlgorithm.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2021 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Qingyuan Liu *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <pxd/calibration/PXDValidationAlgorithm.h>
12 //#include <pxd/dbobjects/PXDGainMapPar.h>
13 //#include <calibration/dbobjects/TestCalibMean.h>
14 
15 #include <string>
16 #include <algorithm>
17 #include <set>
18 
19 #include <sstream>
20 #include <iostream>
21 
22 #include <boost/format.hpp>
23 
24 //ROOT
25 #include <TRandom.h>
26 #include <TH1.h>
27 #include <TF1.h>
28 
29 using namespace std;
30 using boost::format;
31 using namespace Belle2;
32 
33 
34 // Anonymous namespace for data objects used by PXDValidationAlgorithm class
35 namespace {
36 
38  void getNumberOfBins(const std::shared_ptr<TH1I>& histo_ptr, unsigned short& nBinsU, unsigned short& nBinsV)
39  {
40  set<unsigned short> uBinSet;
41  set<unsigned short> vBinSet;
42 
43  // Loop over all bins of input histo
44  for (auto histoBin = 1; histoBin <= histo_ptr->GetXaxis()->GetNbins(); histoBin++) {
45  // The bin label contains the vxdid, uBin and vBin
46  string label = histo_ptr->GetXaxis()->GetBinLabel(histoBin);
47 
48  // Parse label string format to read sensorID, uBin and vBin
49  istringstream stream(label);
50  string token;
51  getline(stream, token, '_');
52  getline(stream, token, '_');
53  unsigned short uBin = std::stoi(token);
54 
55  getline(stream, token, '_');
56  unsigned short vBin = std::stoi(token);
57 
58  uBinSet.insert(uBin);
59  vBinSet.insert(vBin);
60  }
61 
62  if (uBinSet.empty() || vBinSet.empty()) {
63  B2FATAL("Not able to determine the grid size. Something is wrong with collected data.");
64  } else {
65  nBinsU = *uBinSet.rbegin() + 1;
66  nBinsV = *vBinSet.rbegin() + 1;
67  }
68  }
69 
71  unsigned short getNumberOfSensors(const std::shared_ptr<TH1I>& histo_ptr)
72  {
73  set<unsigned short> sensorSet;
74 
75  // Loop over all bins of input histo
76  for (auto histoBin = 1; histoBin <= histo_ptr->GetXaxis()->GetNbins(); histoBin++) {
77  // The bin label contains the vxdid, uBin and vBin
78  string label = histo_ptr->GetXaxis()->GetBinLabel(histoBin);
79 
80  // Parse label string format to read sensorID, uBin and vBin
81  istringstream stream(label);
82  string token;
83  getline(stream, token, '_');
84  VxdID sensorID(token);
85  sensorSet.insert(sensorID.getID());
86  }
87  return sensorSet.size();
88  }
89 
90 }
91 
92 
93 PXDValidationAlgorithm::PXDValidationAlgorithm():
94  CalibrationAlgorithm("PXDCDSTCollector")
95  , minTrackPoints(1000), save2DHists(false)
96  //,m_file(nullptr), m_tree(nullptr)
97  //,minTrackPoints(1000), safetyFactor(2.0), forceContinue(false), strategy(0)
98 {
99  setDescription(
100  " -------------------------- PXDValidationAlgorithm ---------------------------------\n"
101  " \n"
102  " Algorithm for filling validation histograms. \n"
103  " ----------------------------------------------------------------------------------------\n"
104  );
105 
106 
107 }
108 
109 PXDValidationAlgorithm::~PXDValidationAlgorithm()
110 {
111 }
112 
113 CalibrationAlgorithm::EResult PXDValidationAlgorithm::calibrate()
114 {
115  // Turn off ROOT GUI
116  gROOT->SetBatch();
117 
118  // Current Exp No. and Run No.
119  auto expRuns = getRunList();
120  auto expRunsAll = getRunListFromAllData();
121  // size should be 1
122  if (!expRuns.size())
123  return c_Failure;
124 
125  m_exp = expRuns.back().first;
126  m_run = expRuns.back().second;
127  B2INFO("Current ExpRuns: [("
128  << expRuns.front().first << ", " << expRuns.front().second << "), ("
129  << m_exp << ", " << m_run << ")]");
130 
131  // Get counter histograms and set pointers
132  auto cluster_counter = getObjectPtr<TH1I>("PXDTrackClusterCounter");
133  auto point_counter = getObjectPtr<TH1I>("PXDTrackPointCounter");
134 
135  // Extract number of sensors from counter histograms
136  auto nSensors = getNumberOfSensors(cluster_counter);
137 
138  // Extract the number of grid bins from counter histograms
139  unsigned short nBinsU = 0;
140  unsigned short nBinsV = 0;
141  getNumberOfBins(cluster_counter, nBinsU, nBinsV);
142 
143 
144  B2INFO("Start info collection using a " << nBinsU << "x" << nBinsV << " grid per sensor.");
145  B2INFO("Number of collected track points is " << point_counter->GetEntries()
146  << " in " << nSensors << " sensors.");
147 
148  if (point_counter->GetEntries() < minTrackPoints) {
149  B2WARNING("Not enough Data: Only " << point_counter->GetEntries() << " hits were collected but " << minTrackPoints
150  << " needed!");
151  return c_NotEnoughData;
152  }
153 
154 
155  auto hPassedHitsLayer1 = getObjectPtr<TH2F>("hPassedHitsLayer1");
156  auto hTotalHitsLayer1 = getObjectPtr<TH2F>("hTotalHitsLayer1");
157  auto hPassedHitsLayer2 = getObjectPtr<TH2F>("hPassedHitsLayer2");
158  auto hTotalHitsLayer2 = getObjectPtr<TH2F>("hTotalHitsLayer2");
159  if (!cluster_counter) return c_Failure;
160  if (!point_counter) return c_Failure;
161  if (!hPassedHitsLayer1) return c_Failure;
162  if (!hTotalHitsLayer1) return c_Failure;
163  if (!hPassedHitsLayer2) return c_Failure;
164  if (!hTotalHitsLayer2) return c_Failure;
165 
166  // Save the current directory to change back later
167  TDirectory* currentDir = gDirectory;
168  //if (!currentDir) currentDir = gDirectory;
169 
170  // Create TFile if not exist
171  if (!m_file) {
172  std::string fileName = (this->getPrefix()) + "Validation.root";
173  B2INFO("Creating file " << fileName);
174  m_file = std::make_shared<TFile>(fileName.c_str(), "RECREATE");
175  }
176 
177  // Create TTree if not exist
178  if (!m_tree) {
179 
180  B2INFO("Creating TTree.");
181  m_tree = std::make_shared<TTree>("tree", "PXD validation data");
182  // Define histograms out of m_file
183  currentDir->cd();
184  m_hD0 = new TH1F("hD0", "Corrected d0;#Delta d0/#sqrt{2} [cm];Counts", 100, -0.03, 0.03);
185  m_hZ0 = new TH1F("hZ0", "Corrected z0;#Delta z0/#sqrt{2} [cm];Counts", 100, -0.03, 0.03);
186  m_file->cd();
187 
188  m_tree->Branch<int>("exp", &m_exp);
189  m_tree->Branch<int>("run", &m_run);
190  m_tree->Branch("pxdid", &m_pxdid);
191  m_tree->Branch("uBin", &m_uBin);
192  m_tree->Branch("vBin", &m_vBin);
193  m_tree->Branch("nTrackPoints", &m_nTrackPoints);
194  m_tree->Branch("nTrackClusters", &m_nTrackClusters);
195  m_tree->Branch<TH1F>("hD0", &m_hD0, 32000, 0);
196  m_tree->Branch<TH1F>("hZ0", &m_hZ0, 32000, 0);
197 
198  if (save2DHists) {
199  m_hTrackPointsLayer1 = (TH2F*)hTotalHitsLayer1->Clone();
200  m_hTrackClustersLayer1 = (TH2F*)hPassedHitsLayer1->Clone();
201  m_hTrackPointsLayer2 = (TH2F*)hTotalHitsLayer2->Clone();
202  m_hTrackClustersLayer2 = (TH2F*)hPassedHitsLayer2->Clone();
203 
204  // associated histograms
205  m_tree->Branch<TH2F>("hTrackPointsLayer1", &m_hTrackPointsLayer1, 32000, 0);
206  m_tree->Branch<TH2F>("hTrackClustersLayer1", &m_hTrackClustersLayer1, 32000, 0);
207  m_tree->Branch<TH2F>("hTrackPointsLayer2", &m_hTrackPointsLayer2, 32000, 0);
208  m_tree->Branch<TH2F>("hTrackClustersLayer2", &m_hTrackClustersLayer2, 32000, 0);
209  }
210  }
211 
212  m_file->cd();
213 
214  if (save2DHists) {
215  // Just update 2D histograms
216  *m_hTrackPointsLayer1 = *hTotalHitsLayer1;
217  *m_hTrackClustersLayer1 = *hPassedHitsLayer1;
218  *m_hTrackPointsLayer2 = *hTotalHitsLayer2;
219  *m_hTrackClustersLayer2 = *hPassedHitsLayer2;
220 
221  }
222 
223  // Clear
224  m_pxdid.clear();
225  m_uBin.clear();
226  m_vBin.clear();
227  m_nTrackClusters.clear();
228  m_nTrackPoints.clear();
229 
230  // Get resolution trees and create histograms
231  auto tree_d0z0 = getObjectPtr<TTree>("tree_d0z0");
232  if (!tree_d0z0) return c_Failure;
233  float d0, z0;
234  tree_d0z0->SetBranchAddress("d0", &d0);
235  tree_d0z0->SetBranchAddress("z0", &z0);
236  //TH1F hD0("hD0", "Corrected d0;#Delta d0/#sqrt{2} [cm];Counts", 100, -0.03, 0.03);
237  //TH1F hZ0("hZ0", "Corrected z0;#Delta z0/#sqrt{2} [cm];Counts", 100, -0.03, 0.03);
238  m_hD0->Reset();
239  m_hZ0->Reset();
240  string cuts = "abs(d0)<0.03&&abs(z0)<0.03";
241 
242  // Fill histograms from tree
243  for (int i = 0; i < tree_d0z0->GetEntries(); i++) {
244  tree_d0z0->GetEntry(i);
245  if (fabs(d0) > 0.03 || fabs(z0) > 0.03)
246  continue;
247  m_hD0->Fill(d0);
248  m_hZ0->Fill(z0);
249  }
250 
251 
252 
253 
254  // Loop over all bins of input histo
255  for (auto histoBin = 1; histoBin <= cluster_counter->GetXaxis()->GetNbins(); histoBin++) {
256  // The bin label contains the vxdid, uBin and vBin
257  string label = cluster_counter->GetXaxis()->GetBinLabel(histoBin);
258 
259  // Parse label string format to read sensorID, uBin and vBin
260  istringstream stream(label);
261  string token;
262  getline(stream, token, '_');
263  VxdID sensorID(token);
264 
265  getline(stream, token, '_');
266  unsigned short uBin = std::stoi(token);
267 
268  getline(stream, token, '_');
269  unsigned short vBin = std::stoi(token);
270 
271  // Read back the counters for number of collected clusters
272  int numberOfClusters = cluster_counter->GetBinContent(histoBin);
273  int numberOfPoints = point_counter->GetBinContent(histoBin);
274 
275  m_pxdid.emplace_back(sensorID.getLayerNumber() * 1000 + sensorID.getLadderNumber() * 10 + sensorID.getSensorNumber());
276  m_uBin.emplace_back(uBin);
277  m_vBin.emplace_back(vBin);
278  m_nTrackPoints.emplace_back(numberOfPoints);
279  m_nTrackClusters.emplace_back(numberOfClusters);
280  }
281 
282  m_tree->Fill();
283 
284  // At the last run of the data
285  if (m_exp == expRunsAll.back().first &&
286  m_run == expRunsAll.back().second) {
287  B2INFO("Reached Final ExpRun: (" << m_exp << ", " << m_run << ")");
288  m_file->cd();
289  m_tree->Write();
290  B2INFO("Writing Successful.");
291  //if (m_file) m_file->Close();// Don't do it! We are using shared_ptr!
292  if (m_hD0) delete m_hD0;
293  if (m_hZ0) delete m_hZ0;
294  if (save2DHists) {
295  if (m_hTrackPointsLayer1) delete m_hTrackPointsLayer1;
296  if (m_hTrackClustersLayer1) delete m_hTrackClustersLayer1;
297  if (m_hTrackPointsLayer2) delete m_hTrackPointsLayer2;
298  if (m_hTrackClustersLayer2) delete m_hTrackClustersLayer2;
299  }
300  m_file->Close();
301 
302  }
303  currentDir->cd();
304 
305  // Saving dummy DBObject. Required by CAF?
306  //TestCalibMean* correction = new TestCalibMean(m_run, gRandom->Gaus());
307  //B2INFO("Saving calibration results.");
308  //saveCalibration(correction, "DummyCalibMean");
309 
310  // Always return ok to have run-by-run info
311  return c_OK;
312 }
313 
314 
315 bool PXDValidationAlgorithm::isBoundaryRequired(const Calibration::ExpRun& /*currentRun*/)
316 {
317  // First run in data as we iterate, but our boundaries weren't set manually already?
318  // Just set the first run to be a boundary and we are done.
319  if (m_boundaries.empty()) {
320  B2INFO("This is the first run encountered, let's say it is a boundary.");
321  return true;
322  } else {
323  //return true;// -> simple run-by-run
324  return false;
325  }
326 }
Belle2::PXDValidationAlgorithm::m_hD0
TH1F * m_hD0
Histogram of corrected d0 for each 2-track event.
Definition: PXDValidationAlgorithm.h:105
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::PXDValidationAlgorithm::m_nTrackPoints
std::vector< unsigned long > m_nTrackPoints
Vecotr of number of track points.
Definition: PXDValidationAlgorithm.h:99
Belle2::PXDValidationAlgorithm::m_hTrackPointsLayer2
TH2F * m_hTrackPointsLayer2
Histogram of intersection points for layer 2.
Definition: PXDValidationAlgorithm.h:117
Belle2::PXDValidationAlgorithm::m_vBin
std::vector< unsigned short > m_vBin
Vecotr of vBin.
Definition: PXDValidationAlgorithm.h:96
Belle2::VxdID::getLadderNumber
baseType getLadderNumber() const
Get the ladder id.
Definition: VxdID.h:108
Belle2::PXDValidationAlgorithm::m_hTrackClustersLayer2
TH2F * m_hTrackClustersLayer2
Histogram of track matched clusters for layer 2.
Definition: PXDValidationAlgorithm.h:120
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfuly =0 in Python.
Definition: CalibrationAlgorithm.h:51
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition: CalibrationAlgorithm.h:331
Belle2::PXDValidationAlgorithm::m_exp
int m_exp
Experiment number.
Definition: PXDValidationAlgorithm.h:84
Belle2::PXD::getNumberOfBins
void getNumberOfBins(const std::shared_ptr< TH1I > &histo_ptr, unsigned short &nBinsU, unsigned short &nBinsV)
Helper function to extract number of bins along u side and v side from counter histogram labels.
Definition: PXDCalibrationUtilities.cc:37
Belle2::PXDValidationAlgorithm::m_pxdid
std::vector< unsigned short > m_pxdid
Vector of PXD module id (DHE id)
Definition: PXDValidationAlgorithm.h:90
Belle2::PXDValidationAlgorithm::m_file
std::shared_ptr< TFile > m_file
Pointer for TFile.
Definition: PXDValidationAlgorithm.h:78
Belle2::PXDValidationAlgorithm::~PXDValidationAlgorithm
virtual ~PXDValidationAlgorithm() override
Destructor.
Definition: PXDValidationAlgorithm.cc:109
Belle2::CalibrationAlgorithm::getRunListFromAllData
std::vector< Calibration::ExpRun > getRunListFromAllData() const
Get the complete list of runs from inspection of collected data.
Definition: CalibrationAlgorithm.cc:311
Belle2::PXDValidationAlgorithm::calibrate
virtual EResult calibrate() override
Run algo on data.
Definition: PXDValidationAlgorithm.cc:113
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::CalibrationAlgorithm::c_Failure
@ c_Failure
Failed =3 in Python.
Definition: CalibrationAlgorithm.h:54
Belle2::PXDValidationAlgorithm::save2DHists
bool save2DHists
Flag to save 2D histograms for efficiency;.
Definition: PXDValidationAlgorithm.h:49
Belle2::PXDValidationAlgorithm::m_hZ0
TH1F * m_hZ0
Histogram of corrected z0 for each 2-track event.
Definition: PXDValidationAlgorithm.h:108
Belle2::PXDValidationAlgorithm::m_tree
std::shared_ptr< TTree > m_tree
Pointer for TTree of the validation info.
Definition: PXDValidationAlgorithm.h:81
Belle2::VxdID::getSensorNumber
baseType getSensorNumber() const
Get the sensor id.
Definition: VxdID.h:110
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition: CalibrationAlgorithm.h:50
Belle2::CalibrationAlgorithm::getRunList
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
Definition: CalibrationAlgorithm.h:276
Belle2::PXDValidationAlgorithm::minTrackPoints
int minTrackPoints
Minimum number of track points per sensor.
Definition: PXDValidationAlgorithm.h:46
Belle2::PXDValidationAlgorithm::m_hTrackClustersLayer1
TH2F * m_hTrackClustersLayer1
Histogram of track matched clusters for layer 1.
Definition: PXDValidationAlgorithm.h:114
Belle2::CalibrationAlgorithm::c_NotEnoughData
@ c_NotEnoughData
Needs more data =2 in Python.
Definition: CalibrationAlgorithm.h:53
Belle2::PXDValidationAlgorithm::isBoundaryRequired
virtual bool isBoundaryRequired(const Calibration::ExpRun &) override
Decide if a run should be a payload boundary. Only used in certain Python Algorithm Starategies.
Definition: PXDValidationAlgorithm.cc:315
Belle2::PXDValidationAlgorithm::m_hTrackPointsLayer1
TH2F * m_hTrackPointsLayer1
Histogram of intersection points for layer 1.
Definition: PXDValidationAlgorithm.h:111
Belle2::CalibrationAlgorithm
Base class for calibration algorithms.
Definition: CalibrationAlgorithm.h:47
Belle2::VxdID::getLayerNumber
baseType getLayerNumber() const
Get the layer id.
Definition: VxdID.h:106
Belle2::PXDValidationAlgorithm::m_nTrackClusters
std::vector< unsigned long > m_nTrackClusters
Vecotr of number of track matched clusters.
Definition: PXDValidationAlgorithm.h:102
Belle2::PXDValidationAlgorithm::m_uBin
std::vector< unsigned short > m_uBin
Vecotr of uBin.
Definition: PXDValidationAlgorithm.h:93
Belle2::PXDValidationAlgorithm::m_run
int m_run
Run number.
Definition: PXDValidationAlgorithm.h:87
Belle2::PXD::getNumberOfSensors
unsigned short getNumberOfSensors(const std::shared_ptr< TH1I > &histo_ptr)
Helper function to extract number of sensors from counter histogram labels.
Definition: PXDCalibrationUtilities.cc:70
Belle2::CalibrationAlgorithm::getPrefix
std::string getPrefix() const
Get the prefix used for getting calibration data.
Definition: CalibrationAlgorithm.h:156