Belle II Software development
PXDValidationAlgorithm.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
9#include <pxd/calibration/PXDValidationAlgorithm.h>
10#include <vxd/dataobjects/VxdID.h>
11
12#include <string>
13#include <algorithm>
14#include <set>
15
16#include <sstream>
17#include <iostream>
18
19#include <boost/format.hpp>
20
21//ROOT
22#include <TRandom.h>
23#include <TH1.h>
24#include <TF1.h>
25
26using namespace std;
27using boost::format;
28using namespace Belle2;
29
30
31// Anonymous namespace for data objects used by PXDValidationAlgorithm class
32namespace {
33
35 void getNumberOfBins(const std::shared_ptr<TH1I>& histo_ptr, unsigned short& nBinsU, unsigned short& nBinsV)
36 {
37 set<unsigned short> uBinSet;
38 set<unsigned short> vBinSet;
39
40 // Loop over all bins of input histo
41 for (auto histoBin = 1; histoBin <= histo_ptr->GetXaxis()->GetNbins(); histoBin++) {
42 // The bin label contains the vxdid, uBin and vBin
43 string label = histo_ptr->GetXaxis()->GetBinLabel(histoBin);
44
45 // Parse label string format to read sensorID, uBin and vBin
46 istringstream stream(label);
47 string token;
48 getline(stream, token, '_');
49 getline(stream, token, '_');
50 unsigned short uBin = std::stoi(token);
51
52 getline(stream, token, '_');
53 unsigned short vBin = std::stoi(token);
54
55 uBinSet.insert(uBin);
56 vBinSet.insert(vBin);
57 }
58
59 if (uBinSet.empty() || vBinSet.empty()) {
60 B2FATAL("Not able to determine the grid size. Something is wrong with collected data.");
61 } else {
62 nBinsU = *uBinSet.rbegin() + 1;
63 nBinsV = *vBinSet.rbegin() + 1;
64 }
65 }
66
68 unsigned short getNumberOfSensors(const std::shared_ptr<TH1I>& histo_ptr)
69 {
70 set<unsigned short> sensorSet;
71
72 // Loop over all bins of input histo
73 for (auto histoBin = 1; histoBin <= histo_ptr->GetXaxis()->GetNbins(); histoBin++) {
74 // The bin label contains the vxdid, uBin and vBin
75 string label = histo_ptr->GetXaxis()->GetBinLabel(histoBin);
76
77 // Parse label string format to read sensorID, uBin and vBin
78 istringstream stream(label);
79 string token;
80 getline(stream, token, '_');
81 VxdID sensorID(token);
82 sensorSet.insert(sensorID.getID());
83 }
84 return sensorSet.size();
85 }
86
87}
88
89
91 CalibrationAlgorithm("PXDCDSTCollector")
92 , minTrackPoints(1000), save2DHists(false), saveD0Z0(false), binsD0Z0(600)
93 , m_exp(-1), m_run(-1), m_hD0(nullptr), m_hZ0(nullptr)
96{
98 " -------------------------- PXDValidationAlgorithm ---------------------------------\n"
99 " \n"
100 " Algorithm for filling validation histograms. \n"
101 " ----------------------------------------------------------------------------------------\n"
102 );
103
104
105}
106
110
112{
113 // Turn off ROOT GUI
114 gROOT->SetBatch();
115
116 // Current Exp No. and Run No.
117 auto expRuns = getRunList();
118 auto expRunsAll = getRunListFromAllData();
119 // size should be 1
120 if (!expRuns.size())
121 return c_Failure;
122
123 m_exp = expRuns.back().first;
124 m_run = expRuns.back().second;
125 B2INFO("Current ExpRuns: [("
126 << expRuns.front().first << ", " << expRuns.front().second << "), ("
127 << m_exp << ", " << m_run << ")]");
128
129 // Get counter histograms and set pointers
130 auto cluster_counter = getObjectPtr<TH1I>("PXDTrackClusterCounter");
131 auto point_counter = getObjectPtr<TH1I>("PXDTrackPointCounter");
132 auto selected_point_counter = getObjectPtr<TH1I>("PXDSelTrackPointCounter"); // can be empty
133 auto selected_cluster_counter = getObjectPtr<TH1I>("PXDSelTrackClusterCounter"); // can be empty
134 if (!cluster_counter) return c_NotEnoughData;
135 if (!point_counter) return c_NotEnoughData;
136
137 // Extract number of sensors from counter histograms
138 auto nSensors = getNumberOfSensors(cluster_counter);
139
140 // Extract the number of grid bins from counter histograms
141 unsigned short nBinsU = 0;
142 unsigned short nBinsV = 0;
143 getNumberOfBins(cluster_counter, nBinsU, nBinsV);
144
145
146 B2INFO("Start info collection using a " << nBinsU << "x" << nBinsV << " grid per sensor.");
147 B2INFO("Number of collected track points is " << point_counter->GetEntries()
148 << " in " << nSensors << " sensors.");
149
150 if (point_counter->GetEntries() < minTrackPoints) {
151 B2WARNING("Not enough Data: Only " << point_counter->GetEntries() << " hits were collected but " << minTrackPoints
152 << " needed!");
153 return c_NotEnoughData;
154 }
155
156
157 auto hPassedHitsLayer1 = getObjectPtr<TH2F>("hPassedHitsLayer1");
158 auto hTotalHitsLayer1 = getObjectPtr<TH2F>("hTotalHitsLayer1");
159 auto hPassedHitsLayer2 = getObjectPtr<TH2F>("hPassedHitsLayer2");
160 auto hTotalHitsLayer2 = getObjectPtr<TH2F>("hTotalHitsLayer2");
161 if (!hPassedHitsLayer1) return c_NotEnoughData;
162 if (!hTotalHitsLayer1) return c_NotEnoughData;
163 if (!hPassedHitsLayer2) return c_NotEnoughData;
164 if (!hTotalHitsLayer2) return c_NotEnoughData;
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", binsD0Z0, -0.03, 0.03);
185 m_hZ0 = new TH1F("hZ0", "Corrected z0;#Delta z0/#sqrt{2} [cm];Counts", binsD0Z0, -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("nSelTrackPoints", &m_nSelTrackPoints);
196 m_tree->Branch("nSelTrackClusters", &m_nSelTrackClusters);
197 m_tree->Branch<TH1F>("hD0", &m_hD0, 32000, 0);
198 m_tree->Branch<TH1F>("hZ0", &m_hZ0, 32000, 0);
199
200 if (saveD0Z0) {
201 m_tree->Branch("d0", &m_d0);
202 m_tree->Branch("z0", &m_z0);
203 }
204
205 if (save2DHists) {
206 m_hTrackPointsLayer1 = (TH2F*)hTotalHitsLayer1->Clone();
207 m_hTrackClustersLayer1 = (TH2F*)hPassedHitsLayer1->Clone();
208 m_hTrackPointsLayer2 = (TH2F*)hTotalHitsLayer2->Clone();
209 m_hTrackClustersLayer2 = (TH2F*)hPassedHitsLayer2->Clone();
210
211 // associated histograms
212 m_tree->Branch<TH2F>("hTrackPointsLayer1", &m_hTrackPointsLayer1, 32000, 0);
213 m_tree->Branch<TH2F>("hTrackClustersLayer1", &m_hTrackClustersLayer1, 32000, 0);
214 m_tree->Branch<TH2F>("hTrackPointsLayer2", &m_hTrackPointsLayer2, 32000, 0);
215 m_tree->Branch<TH2F>("hTrackClustersLayer2", &m_hTrackClustersLayer2, 32000, 0);
216 }
217 }
218
219 m_file->cd();
220
221 if (save2DHists) {
222 // Just update 2D histograms
223 *m_hTrackPointsLayer1 = *hTotalHitsLayer1;
224 *m_hTrackClustersLayer1 = *hPassedHitsLayer1;
225 *m_hTrackPointsLayer2 = *hTotalHitsLayer2;
226 *m_hTrackClustersLayer2 = *hPassedHitsLayer2;
227
228 }
229
230 // Clear
231 m_pxdid.clear();
232 m_uBin.clear();
233 m_vBin.clear();
234 m_nTrackClusters.clear();
235 m_nTrackPoints.clear();
236 m_nSelTrackPoints.clear();
237 m_nSelTrackClusters.clear();
238 m_d0.clear();
239 m_z0.clear();
240
241 // Get resolution trees and create histograms
242 auto tree_d0z0 = getObjectPtr<TTree>("tree_d0z0");
243 if (!tree_d0z0) return c_NotEnoughData;
244 float d0, z0;
245 tree_d0z0->SetBranchAddress("d0", &d0);
246 tree_d0z0->SetBranchAddress("z0", &z0);
247 m_hD0->Reset();
248 m_hZ0->Reset();
249
250 // Fill histograms from tree
251 for (int i = 0; i < tree_d0z0->GetEntries(); i++) {
252 tree_d0z0->GetEntry(i);
253 if (saveD0Z0) {
254 m_d0.emplace_back(d0);
255 m_z0.emplace_back(z0);
256 }
257 if (fabs(d0) > 0.03 || fabs(z0) > 0.03)
258 continue;
259 m_hD0->Fill(d0);
260 m_hZ0->Fill(z0);
261 }
262
263
264 // Loop over all bins of input histo
265 for (auto histoBin = 1; histoBin <= cluster_counter->GetXaxis()->GetNbins(); histoBin++) {
266 // The bin label contains the vxdid, uBin and vBin
267 string label = cluster_counter->GetXaxis()->GetBinLabel(histoBin);
268
269 // Parse label string format to read sensorID, uBin and vBin
270 istringstream stream(label);
271 string token;
272 getline(stream, token, '_');
273 VxdID sensorID(token);
274
275 getline(stream, token, '_');
276 unsigned short uBin = std::stoi(token);
277
278 getline(stream, token, '_');
279 unsigned short vBin = std::stoi(token);
280
281 // Read back the counters for number of collected clusters
282 int numberOfClusters = cluster_counter->GetBinContent(histoBin);
283 int numberOfPoints = point_counter->GetBinContent(histoBin);
284
285 m_pxdid.emplace_back(sensorID.getLayerNumber() * 1000 + sensorID.getLadderNumber() * 10 + sensorID.getSensorNumber());
286 m_uBin.emplace_back(uBin);
287 m_vBin.emplace_back(vBin);
288 m_nTrackPoints.emplace_back(numberOfPoints);
289 m_nSelTrackPoints.push_back(selected_point_counter->GetBinContent(histoBin));
290 m_nSelTrackClusters.push_back(selected_cluster_counter->GetBinContent(histoBin));
291 m_nTrackClusters.emplace_back(numberOfClusters);
292 }
293
294 m_tree->Fill();
295
296 // At the last run of the data
297 if (m_exp == expRunsAll.back().first &&
298 m_run == expRunsAll.back().second) {
299 B2INFO("Reached Final ExpRun: (" << m_exp << ", " << m_run << ")");
300 m_file->cd();
301 m_tree->Write();
302 B2INFO("Writing Successful.");
303 //if (m_file) m_file->Close();// Don't do it! We are using shared_ptr!
304 if (m_hD0) delete m_hD0;
305 if (m_hZ0) delete m_hZ0;
306 if (save2DHists) {
311 }
312 m_file->Close();
313
314 }
315 currentDir->cd();
316
317 // Always return ok to have run-by-run info
318 return c_OK;
319}
320
321
322bool PXDValidationAlgorithm::isBoundaryRequired(const Calibration::ExpRun& /*currentRun*/)
323{
324 // First run in data as we iterate, but our boundaries weren't set manually already?
325 // Just set the first run to be a boundary and we are done.
326 if (m_boundaries.empty()) {
327 B2INFO("This is the first run encountered, let's say it is a boundary.");
328 return true;
329 } else {
330 //return true;// -> simple run-by-run
331 return false;
332 }
333}
std::vector< Calibration::ExpRun > m_boundaries
When using the boundaries functionality from isBoundaryRequired, this is used to store the boundaries...
std::vector< Calibration::ExpRun > getRunListFromAllData() const
Get the complete list of runs from inspection of collected data.
void setDescription(const std::string &description)
Set algorithm description (in constructor)
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
std::string getPrefix() const
Get the prefix used for getting calibration data.
EResult
The result of calibration.
@ c_OK
Finished successfully =0 in Python.
@ c_NotEnoughData
Needs more data =2 in Python.
CalibrationAlgorithm(const std::string &collectorModuleName)
Constructor - sets the prefix for collected objects (won't be accesses until execute(....
TH1F * m_hD0
Histogram of corrected d0 for each 2-track event.
TH2F * m_hTrackPointsLayer2
Histogram of intersection points for layer 2.
TH2F * m_hTrackClustersLayer1
Histogram of track matched clusters for layer 1.
TH2F * m_hTrackPointsLayer1
Histogram of intersection points for layer 1.
std::vector< unsigned long > m_nTrackClusters
Vector of number of track matched clusters.
virtual ~PXDValidationAlgorithm() override
Destructor.
TH2F * m_hTrackClustersLayer2
Histogram of track matched clusters for layer 2.
std::vector< unsigned long > m_nTrackPoints
Vector of number of track points.
std::vector< float > m_z0
Vector of delta z0 over sqrt(2)
bool saveD0Z0
Flag to save delta d0 (z0) over sqrt(2) to TTree;.
std::vector< unsigned long > m_nSelTrackPoints
Vector of number of track points outside of defective pixels.
PXDValidationAlgorithm()
Constructor set the prefix to PXDValidationAlgorithm.
int minTrackPoints
Minimum number of track points per sensor.
std::shared_ptr< TFile > m_file
Pointer for TFile.
bool save2DHists
Flag to save 2D histograms for efficiency;.
std::vector< unsigned short > m_uBin
Vector of uBin.
std::vector< float > m_d0
Vector of delta d0 over sqrt(2)
virtual EResult calibrate() override
Run algo on data.
int binsD0Z0
Bin size for d0/z0 histogram.
std::vector< unsigned long > m_nSelTrackClusters
Vector of number of track clusters outside of defective pixels.
std::vector< unsigned short > m_vBin
Vector of vBin.
std::shared_ptr< TTree > m_tree
Pointer for TTree of the validation info.
virtual bool isBoundaryRequired(const Calibration::ExpRun &) override
Decide if a run should be a payload boundary. Only used in certain Python Algorithm Starategies.
TH1F * m_hZ0
Histogram of corrected z0 for each 2-track event.
std::vector< unsigned short > m_pxdid
Vector of PXD module id (DHE id)
Class to uniquely identify a any structure of the PXD and SVD.
Definition VxdID.h:32
baseType getSensorNumber() const
Get the sensor id.
Definition VxdID.h:99
baseType getLadderNumber() const
Get the ladder id.
Definition VxdID.h:97
baseType getLayerNumber() const
Get the layer id.
Definition VxdID.h:95
std::shared_ptr< T > getObjectPtr(const std::string &name, const std::vector< Calibration::ExpRun > &requestedRuns)
Get calibration data object by name and list of runs, the Merge function will be called to generate t...
unsigned short getNumberOfSensors(const std::shared_ptr< TH1I > &histo_ptr)
Helper function to extract number of sensors from counter histogram labels.
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.
Abstract base class for different kinds of events.
STL namespace.