Belle II Software development
SVDPositionErrorScaleFactorImporterModule.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 <svd/modules/svdCalibration/SVDPositionErrorScaleFactorImporterModule.h>
10#include <svd/calibration/SVDCoGOnlyErrorScaleFactors.h>
11#include <svd/calibration/SVDOldDefaultErrorScaleFactors.h>
12#include <vxd/geometry/GeoCache.h>
13#include <framework/datastore/StoreObjPtr.h>
14#include <framework/datastore/StoreArray.h>
15#include <framework/datastore/RelationVector.h>
16#include <framework/datastore/RelationArray.h>
17#include <framework/dataobjects/EventMetaData.h>
18
19using namespace Belle2;
20
21//-----------------------------------------------------------------
22// Register the Module
23//-----------------------------------------------------------------
24REG_MODULE(SVDPositionErrorScaleFactorImporter);
25
26//-----------------------------------------------------------------
27// Implementation
28//-----------------------------------------------------------------
29
31{
32 // Set module properties
33 setDescription("Module to produce a list of histograms showing the uploaded calibration constants");
34
35 // Parameter definitions
36 addParam("outputFileName", m_rootFileName, "Name of output root file.", std::string("SVDPositionErrorScaleFactors.root"));
37 addParam("posAlgorithm", m_posAlgorithm, "Position Algorithm.", std::string(m_posAlgorithm));
38 addParam("uniqueID", m_uniqueID, "Payload uniqueID.", std::string(m_uniqueID));
39 addParam("minPulls", m_min, "min of the pulls histograms.", float(m_min));
40 addParam("maxPulls", m_max, "max of the pulls histograms.", float(m_max));
41 addParam("nBinsPulls", m_nBins, "number of bins of the pulls histograms.", float(m_nBins));
42 addParam("noOutliers", m_noOutliers, "If True, removes outliers from cluster position error scale factor computation.",
43 bool(false));
44}
45
47{
49 m_truehits.isRequired();
50}
51
53{
54 //to avoid publicAllocationError:
55 delete m_rootFilePtr;
56
57 m_rootFilePtr = new TFile(m_rootFileName.c_str(), "RECREATE");
58 //tree initialization
59 m_tree = new TTree("scf", "RECREATE");
60 m_tree->Branch("exp", &m_exp, "exp/i");
61 m_tree->Branch("run", &m_run, "run/i");
62 m_tree->Branch("layer", &m_layer, "layer/i");
63 m_tree->Branch("ladder", &m_ladder, "ladder/i");
64 m_tree->Branch("sensor", &m_sensor, "sensor/i");
65 m_tree->Branch("side", &m_side, "side/i");
66 m_tree->Branch("size", &m_size, "size/i");
67 m_tree->Branch("clsCharge", &m_clsCharge, "clsCharge/F");
68 m_tree->Branch("clsTime", &m_clsTime, "clsTime/F");
69 m_tree->Branch("clsPos", &m_clsPos, "clsPos/F");
70 m_tree->Branch("clsErr", &m_clsErr, "clsErr/F");
71 m_tree->Branch("clsResid", &m_clsResid, "clsResid/F");
72 m_tree->Branch("clsPull", &m_clsPull, "clsPull/F");
73
74 //CLUSTER POSITION PULLS
75 TH1F hClsPullSize1("clusterPulls1_L@layerL@ladderS@sensor@view",
76 Form("Cluster %s Pulls for Size 1 in @layer.@ladder.@sensor @view/@side", m_posAlgorithm.c_str()),
78 hClsPullSize1.GetXaxis()->SetTitle("cluster pull");
79 m_hClsPullSize1 = new SVDHistograms<TH1F>(hClsPullSize1);
80
81 TH1F hClsPullSize2("clusterPulls2_L@layerL@ladderS@sensor@view",
82 Form("Cluster %s Pulls for Size 2 in @layer.@ladder.@sensor @view/@side", m_posAlgorithm.c_str()),
84 hClsPullSize2.GetXaxis()->SetTitle("cluster pull");
85 m_hClsPullSize2 = new SVDHistograms<TH1F>(hClsPullSize2);
86
87 TH1F hClsPullSize3("clusterPulls3_L@layerL@ladderS@sensor@view",
88 Form("Cluster %s Pulls for Size 3 in @layer.@ladder.@sensor @view/@side", m_posAlgorithm.c_str()),
90 hClsPullSize3.GetXaxis()->SetTitle("cluster pull");
91 m_hClsPullSize3 = new SVDHistograms<TH1F>(hClsPullSize3);
92
93 TH1F hClsPullSize4("clusterPulls4_L@layerL@ladderS@sensor@view",
94 Form("Cluster %s Pulls for Size 4 in @layer.@ladder.@sensor @view/@side", m_posAlgorithm.c_str()),
96 hClsPullSize4.GetXaxis()->SetTitle("cluster pull");
97 m_hClsPullSize4 = new SVDHistograms<TH1F>(hClsPullSize4);
98
99 TH1F hClsPullSize5("clusterPulls5_L@layerL@ladderS@sensor@view",
100 Form("Cluster %s Pulls for Size > 4 in @layer.@ladder.@sensor @view/@side", m_posAlgorithm.c_str()),
102 hClsPullSize5.GetXaxis()->SetTitle("cluster pull");
103 m_hClsPullSize5 = new SVDHistograms<TH1F>(hClsPullSize5);
104
105 for (int i = 0; i < 2; i++)
106 for (int s = 0; s < maxSize; s++) {
107 TString sside = "u";
108 if (i == 0) sside = "v";
109 m_hL3Pulls[s][i] = new TH1F(Form("l3_%s_size%d", sside.Data(), s + 1), Form("Size %d Cluster %s Pulls for L3 %s sensors",
110 s + 1, m_posAlgorithm.c_str(), sside.Data()), m_nBins, m_min, m_max);
111 m_hL3Pulls[s][i]->GetXaxis()->SetTitle("cluster pull");
112 m_hBWPulls[s][i] = new TH1F(Form("bw_%s_size%d", sside.Data(), s + 1), Form("Size %d Cluster %s Pulls for BW %s sensors",
113 s + 1, m_posAlgorithm.c_str(), sside.Data()), m_nBins, m_min, m_max);
114 m_hBWPulls[s][i]->GetXaxis()->SetTitle("cluster pull");
115 m_hFWPulls[s][i] = new TH1F(Form("fw_%s_size%d", sside.Data(), s + 1), Form("Size %d Cluster %s Pulls for FW %s sensors",
116 s + 1, m_posAlgorithm.c_str(), sside.Data()), m_nBins, m_min, m_max);
117 m_hFWPulls[s][i]->GetXaxis()->SetTitle("cluster pull");
118 m_hORPulls[s][i] = new TH1F(Form("or_%s_size%d", sside.Data(), s + 1), Form("Size %d Cluster %s Pulls for ORIGAMI %s sensors",
119 s + 1, m_posAlgorithm.c_str(), sside.Data()), m_nBins, m_min, m_max);
120 m_hORPulls[s][i]->GetXaxis()->SetTitle("cluster pull");
121 }
122
123
124}
125
127{
129 m_run = meta->getRun();
130 m_exp = meta->getExperiment();
131
132 for (const SVDCluster& cluster : m_clusters) {
133
134 RelationVector<SVDTrueHit> trueHit = cluster.getRelationsTo<SVDTrueHit>();
135
136 if (trueHit.size() == 0) continue;
137
138 m_size = cluster.getSize();
139 m_side = cluster.isUCluster();
140 m_clsCharge = cluster.getCharge();
141 m_clsTime = cluster.getClsTime();
142 m_clsPos = m_side ? cluster.getPosition(trueHit[0]->getV()) : cluster.getPosition();
143 m_clsErr = cluster.getPositionSigma();
144 m_clsResid = m_side ? m_clsPos - trueHit[0]->getU() : m_clsPos - trueHit[0]->getV();
146
147 m_layer = cluster.getSensorID().getLayerNumber();
148 m_ladder = cluster.getSensorID().getLadderNumber();
149 m_sensor = cluster.getSensorID().getSensorNumber();
150
151 m_tree->Fill();
152
153 if (m_size == 1) m_hClsPullSize1->fill(cluster.getSensorID(), m_side, m_clsPull);
154 else if (m_size == 2) m_hClsPullSize2->fill(cluster.getSensorID(), m_side, m_clsPull);
155 else if (m_size == 3) m_hClsPullSize3->fill(cluster.getSensorID(), m_side, m_clsPull);
156 else if (m_size == 4) m_hClsPullSize4->fill(cluster.getSensorID(), m_side, m_clsPull);
157 else m_hClsPullSize5->fill(cluster.getSensorID(), m_side, m_clsPull);
158
159 }
160}
161
163{
164
165 // 1. compute scale factors for each sensor (not used in the payload!)
166 // and add them to histogram title
167
168 //call for a geometry instance
170 std::set<Belle2::VxdID> svdLayers = aGeometry.getLayers(VXD::SensorInfoBase::SVD);
171 std::set<Belle2::VxdID>::iterator itSvdLayers = svdLayers.begin();
172
173 while ((itSvdLayers != svdLayers.end()) && (itSvdLayers->getLayerNumber() != 7)) { //loop on Layers
174
175 std::set<Belle2::VxdID> svdLadders = aGeometry.getLadders(*itSvdLayers);
176 std::set<Belle2::VxdID>::iterator itSvdLadders = svdLadders.begin();
177
178 while (itSvdLadders != svdLadders.end()) { //loop on Ladders
179
180 std::set<Belle2::VxdID> svdSensors = aGeometry.getSensors(*itSvdLadders);
181 std::set<Belle2::VxdID>::iterator itSvdSensors = svdSensors.begin();
182 B2DEBUG(1, " svd sensor info " << * (svdSensors.begin()));
183
184 while (itSvdSensors != svdSensors.end()) { //loop on sensors
185 B2DEBUG(1, " svd sensor info " << *itSvdSensors);
186
187 m_layer = (int)itSvdSensors->getLayerNumber();
188 m_ladder = (int)itSvdSensors->getLadderNumber();
189 m_sensor = (int)itSvdSensors->getSensorNumber();
191
192 for (m_side = 0; m_side < 2; m_side++) {
193
194 TH1F* h = m_hClsPullSize1->getHistogram(theVxdID, m_side);
195 h->SetTitle(Form("%s, scf = %1.2f", h->GetTitle(), oneSigma(h)));
196 h = m_hClsPullSize2->getHistogram(theVxdID, m_side);
197 h->SetTitle(Form("%s, scf = %1.2f", h->GetTitle(), oneSigma(h)));
198 h = m_hClsPullSize3->getHistogram(theVxdID, m_side);
199 h->SetTitle(Form("%s, scf = %1.2f", h->GetTitle(), oneSigma(h)));
200 h = m_hClsPullSize4->getHistogram(theVxdID, m_side);
201 h->SetTitle(Form("%s, scf = %1.2f", h->GetTitle(), oneSigma(h)));
202 h = m_hClsPullSize5->getHistogram(theVxdID, m_side);
203 h->SetTitle(Form("%s, scf = %1.2f", h->GetTitle(), oneSigma(h)));
204
205 //layer 3
206 if (theVxdID.getLayerNumber() == 3) {
212 }
213
214 //forward
215 else if (theVxdID.getSensorNumber() == 1) {
221 }
222 //backward
223 else if (theVxdID.getSensorNumber() == theVxdID.getLayerNumber() - 1) {
229 }
230 // origami
231 else {
237 }
238 }
239 ++itSvdSensors;
240 }
241 ++itSvdLadders;
242 }
243 ++itSvdLayers;
244 }
245
246
247 IntervalOfValidity iov(0, 0, -1, -1);
248
249 auto scfs = new Belle2::SVDPosErrScaleFactors;
250
251 //CoGOnly payload
252 auto payload_cogOnly = new Belle2::SVDCoGOnlyErrorScaleFactors::t_payload(*scfs, m_uniqueID);
253
254 //OldDefault payload
255 auto payload_oldDefault = new Belle2::SVDOldDefaultErrorScaleFactors::t_payload(*scfs, m_uniqueID);
256
257 // 2. compute scale factors using cumulative histograms
259
260 for (auto layer : geoCache.getLayers(VXD::SensorInfoBase::SVD))
261 for (auto ladder : geoCache.getLadders(layer))
262 for (Belle2::VxdID theVxdID : geoCache.getSensors(ladder))
263
264 for (m_side = 0; m_side < 2; m_side++) {
265
266 //layer 3
267 if (theVxdID.getLayerNumber() == 3) {
268 scfs->scaleError_clSize1 = oneSigma(m_hL3Pulls[0][m_side]);
269 scfs->scaleError_clSize2 = oneSigma(m_hL3Pulls[1][m_side]);
270 scfs->scaleError_clSize3 = oneSigma(m_hL3Pulls[2][m_side]);
271 scfs->scaleError_clSize4 = oneSigma(m_hL3Pulls[3][m_side]);
272 scfs->scaleError_clSize5 = oneSigma(m_hL3Pulls[4][m_side]);
273 }
274 //forward
275 else if (theVxdID.getSensorNumber() == 1) {
276 scfs->scaleError_clSize1 = oneSigma(m_hFWPulls[0][m_side]);
277 scfs->scaleError_clSize2 = oneSigma(m_hFWPulls[1][m_side]);
278 scfs->scaleError_clSize3 = oneSigma(m_hFWPulls[2][m_side]);
279 scfs->scaleError_clSize4 = oneSigma(m_hFWPulls[3][m_side]);
280 scfs->scaleError_clSize5 = oneSigma(m_hFWPulls[4][m_side]);
281 }
282 //backward
283 else if (theVxdID.getSensorNumber() == theVxdID.getLayerNumber() - 1) {
284 scfs->scaleError_clSize1 = oneSigma(m_hBWPulls[0][m_side]);
285 scfs->scaleError_clSize2 = oneSigma(m_hBWPulls[1][m_side]);
286 scfs->scaleError_clSize3 = oneSigma(m_hBWPulls[2][m_side]);
287 scfs->scaleError_clSize4 = oneSigma(m_hBWPulls[3][m_side]);
288 scfs->scaleError_clSize5 = oneSigma(m_hBWPulls[4][m_side]);
289 }
290 // origami
291 else {
292 scfs->scaleError_clSize1 = oneSigma(m_hORPulls[0][m_side]);
293 scfs->scaleError_clSize2 = oneSigma(m_hORPulls[1][m_side]);
294 scfs->scaleError_clSize3 = oneSigma(m_hORPulls[2][m_side]);
295 scfs->scaleError_clSize4 = oneSigma(m_hORPulls[3][m_side]);
296 scfs->scaleError_clSize5 = oneSigma(m_hORPulls[4][m_side]);
297 }
298
299 if (TString(m_posAlgorithm).Contains("CoGOnly"))
300 payload_cogOnly->set(theVxdID.getLayerNumber(), theVxdID.getLadderNumber(), theVxdID.getSensorNumber(), m_side, 1, *scfs);
301 if (TString(m_posAlgorithm).Contains("OldDefault"))
302 payload_oldDefault->set(theVxdID.getLayerNumber(), theVxdID.getLadderNumber(), theVxdID.getSensorNumber(), m_side, 1, *scfs);
303 }
304
305 if (TString(m_posAlgorithm).Contains("CoGOnly")) {
307
308 B2RESULT("SVDCoGOnlyErrorScaleFactors imported to database.");
309 }
310 if (TString(m_posAlgorithm).Contains("OldDefault")) {
312
313 B2RESULT("SVDOldDefaultErrorScaleFactors imported to database.");
314 }
315
316
317 //now write the rootfile
318 if (m_rootFilePtr != nullptr) {
319
320 m_rootFilePtr->cd();
321
322 //write the tree
323 m_tree->Write();
324
325 m_rootFilePtr->mkdir("sensor_pulls");
326 m_rootFilePtr->mkdir("pulls");
327
328 for (auto layer : geoCache.getLayers(VXD::SensorInfoBase::SVD))
329 for (auto ladder : geoCache.getLadders(layer))
330 for (Belle2::VxdID sensor : geoCache.getSensors(ladder))
331 for (int view = SVDHistograms<TH1F>::VIndex ; view < SVDHistograms<TH1F>::UIndex + 1; view++) {
332
333 //writing the histograms to root:
334
335 m_rootFilePtr->cd("sensor_pulls");
336 (m_hClsPullSize1->getHistogram(sensor, view))->Write();
337 (m_hClsPullSize2->getHistogram(sensor, view))->Write();
338 (m_hClsPullSize3->getHistogram(sensor, view))->Write();
339 (m_hClsPullSize4->getHistogram(sensor, view))->Write();
340 (m_hClsPullSize5->getHistogram(sensor, view))->Write();
341
342 }
343
344 for (int s = 0; s < maxSize; s++) {
345 m_rootFilePtr->cd("pulls");
346 m_hL3Pulls[s][1]->Write();
347 m_hFWPulls[s][1]->Write();
348 m_hBWPulls[s][1]->Write();
349 m_hORPulls[s][1]->Write();
350 m_hL3Pulls[s][0]->Write();
351 m_hFWPulls[s][0]->Write();
352 m_hBWPulls[s][0]->Write();
353 m_hORPulls[s][0]->Write();
354 }
355 m_rootFilePtr->Close();
356 B2RESULT("The rootfile containing the list of histograms has been filled and closed.");
357
358 }
359}
360
362{
363 TH1F* h1_res = (TH1F*)h1->Clone("h1_res");
364
365 double probs[2] = {0.16, 1 - 0.16};
366 double quant[2] = {0, 0};
367
368 int nbinsHisto = h1_res->GetNbinsX();
369
370 if (m_noOutliers) {
371 h1_res->SetBinContent(1, h1_res->GetBinContent(0) + h1_res->GetBinContent(1));
372 h1_res->SetBinContent(nbinsHisto, h1_res->GetBinContent(nbinsHisto) + h1_res->GetBinContent(nbinsHisto + 1));
373 h1_res->SetBinContent(0, 0);
374 h1_res->SetBinContent(nbinsHisto + 1, 0);
375 }
376
377 h1_res->GetQuantiles(2, quant, probs);
378
379 return (-quant[0] + quant[1]) / 2;
380
381}
A class that describes the interval of experiments/runs for which an object in the database is valid.
Base class for Modules.
Definition: Module.h:72
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
Class for type safe access to objects that are referred to in relations.
size_t size() const
Get number of relations.
The SVD Cluster class This class stores all information about reconstructed SVD clusters.
Definition: SVDCluster.h:29
static std::string name
name of SVDPosErrScaleFactors payload
SVDCalibrationsBase< SVDCalibrationsScalar< SVDPosErrScaleFactors > > t_payload
typedef for the of SVDPosErrScaleFactors payload of all SVD sensors
template class for SVd histograms
Definition: SVDHistograms.h:24
void fill(const VxdID &vxdID, int view, Types ... args)
fill the histogram for
Definition: SVDHistograms.h:77
H * getHistogram(const VxdID &vxdID, int view)
get a reference to the histogram for
Definition: SVDHistograms.h:56
static std::string name
name of SVDPosErrScaleFactors payload
SVDCalibrationsBase< SVDCalibrationsScalar< SVDPosErrScaleFactors > > t_payload
typedef for the of SVDPosErrScaleFactors payload of all SVD sensors
SVDHistograms< TH1F > * m_hClsPullSize1
cluster size 1, position error scale factor histo
TTree * m_tree
pointer at tree containing the mean and RMS of calibration constants
TH1F * m_hFWPulls[5][2]
cumulative FW pulls histograms [size][side]
TH1F * m_hORPulls[5][2]
cumulative ORIGAMI pulls histograms [size][side]
virtual void initialize() override
check presence of data objects
bool m_noOutliers
if True removes outliers from scale factor computation
SVDHistograms< TH1F > * m_hClsPullSize5
cluster size 5, position error scale factor histo
virtual void endRun() override
create cumulative histograms, compute scale factors, create payloads and import it to a localdb,...
SVDPositionErrorScaleFactorImporterModule()
Constructor: Sets the description, the properties and the parameters of the module.
SVDHistograms< TH1F > * m_hClsPullSize3
cluster size 3, position error scale factor histo
virtual void beginRun() override
initialize the TTrees and create histograms
TH1F * m_hL3Pulls[5][2]
cumulative L3 pulls histograms [size][side]
SVDHistograms< TH1F > * m_hClsPullSize2
cluster size 2, position error scale factor histo
SVDHistograms< TH1F > * m_hClsPullSize4
cluster size 4, position error scale factor histo
TH1F * m_hBWPulls[5][2]
cumulative BW pulls histograms [size][side]
double oneSigma(TH1F *)
computes the scale factor by requiring that 68% of the pull distribution is between ±1
TFile * m_rootFilePtr
pointer at root file used for storing histograms
Class SVDTrueHit - Records of tracks that either enter or leave the sensitive volume.
Definition: SVDTrueHit.h:33
bool isRequired(const std::string &name="")
Ensure this array/object has been registered previously.
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:96
Class to faciliate easy access to sensor information of the VXD like coordinate transformations or pi...
Definition: GeoCache.h:39
const std::set< Belle2::VxdID > getLayers(SensorInfoBase::SensorType sensortype=SensorInfoBase::VXD)
Return a set of all known Layers.
Definition: GeoCache.cc:176
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:204
static GeoCache & getInstance()
Return a reference to the singleton instance.
Definition: GeoCache.cc:214
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:193
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:33
baseType getSensorNumber() const
Get the sensor id.
Definition: VxdID.h:100
baseType getLayerNumber() const
Get the layer id.
Definition: VxdID.h:96
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition: Module.h:560
static Database & Instance()
Instance of a singleton Database.
Definition: Database.cc:42
bool storeData(const std::string &name, TObject *object, const IntervalOfValidity &iov)
Store an object in the database.
Definition: Database.cc:141
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Abstract base class for different kinds of events.
contains the scaling factors for the cluster position error