Belle II Software development
T0CalibrationAlgorithm.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#include <cdc/calibration/T0CalibrationAlgorithm.h>
9#include <calibration/CalibrationAlgorithm.h>
10#include <cdc/dbobjects/CDCTimeZeros.h>
11#include <cdc/dataobjects/WireID.h>
12#include <cdc/geometry/CDCGeometryPar.h>
13
14#include <TError.h>
15#include <TROOT.h>
16#include <TGraphErrors.h>
17#include <TF1.h>
18#include <TFile.h>
19#include <TTree.h>
20#include <TStopwatch.h>
21
22#include <framework/database/DBObjPtr.h>
23#include <framework/database/IntervalOfValidity.h>
24#include <framework/logging/Logger.h>
25
26using namespace std;
27using namespace Belle2;
28using namespace CDC;
29
31{
32
34 " -------------------------- T0 Calibration Algorithm -------------------------\n"
35 );
36}
37
39{
40
41 B2INFO("Creating histograms");
42 Float_t x;
43 Float_t t_mea;
44 Float_t t_fit;
45 Float_t ndf;
46 Float_t Pval;
47 UShort_t IWire;
48 UChar_t lay;
49
50 auto tree = getObjectPtr<TTree>("tree");
51 tree->SetBranchAddress("lay", &lay);
52 tree->SetBranchAddress("IWire", &IWire);
53 tree->SetBranchAddress("x_u", &x);
54 tree->SetBranchAddress("t", &t_mea);
55 tree->SetBranchAddress("t_fit", &t_fit);
56 tree->SetBranchAddress("ndf", &ndf);
57 tree->SetBranchAddress("Pval", &Pval);
58
59
60 /* Disable unused branch */
61 std::vector<TString> list_vars = {"lay", "IWire", "x_u", "t", "t_fit", "Pval", "ndf"};
62 tree->SetBranchStatus("*", 0);
63
64 for (TString brname : list_vars) {
65 tree->SetBranchStatus(brname, 1);
66 }
67
68
69 double halfCSize[56];
70
72
73 for (int i = 0; i < 56; ++i) {
74 double R = cdcgeo.senseWireR(i);
75 double nW = cdcgeo.nWiresInLayer(i);
76 halfCSize[i] = M_PI * R / nW;
77 }
78
79 m_hTotal = new TH1F("hTotal", "hTotal", 30, -15, 15);
80
81 //for each channel
82 for (int il = 0; il < 56; il++) {
83 const int nW = cdcgeo.nWiresInLayer(il);
84 for (int ic = 0; ic < nW; ++ic) {
85 m_h1[il][ic] = new TH1F(Form("hdT_L%d_W%d", il, ic), Form("L%d_cell%d", il, ic), 30, -15, 15);
86 }
87 }
88
89 //for each board
90 for (int ib = 0; ib < 300; ++ib) {
91 m_hT0b[ib] = new TH1F(Form("hdT_b%d", ib), Form("boardID_%d", ib), 100, -20, 20);
92 }
93
94 //read data
95 const Long64_t nEntries = tree->GetEntries();
96 B2INFO("Number of entries: " << nEntries);
97 TStopwatch timer;
98 timer.Start();
99 for (Long64_t i = 0; i < nEntries; ++i) {
100 tree->GetEntry(i);
101 double xmax = halfCSize[lay] - 0.1;
102 if ((fabs(x) < m_xmin) || (fabs(x) > xmax)
103 || (ndf < m_ndfmin)
104 || (Pval < m_Pvalmin)) continue; /*select good region*/
105 //each channel
106 m_hTotal->Fill(t_mea - t_fit);
107 m_h1[lay][IWire]->Fill(t_mea - t_fit);
108 //each board
109 int boardID = cdcgeo.getBoardID(WireID(lay, IWire));
110 m_hT0b[boardID]->Fill(t_mea - t_fit);
111 }
112 timer.Stop();
113 B2INFO("Finish making histogram for all channels");
114 B2INFO("Time to fill histograms: " << timer.RealTime() << "s");
115
116}
117
119{
120 B2INFO("Start calibration");
121
122 gROOT->SetBatch(1);
123 gErrorIgnoreLevel = 3001;
124
125 // We are potentially using data from several runs at once during execution
126 // (which may have different DBObject values). So in general you would need to
127 // average them, or aply them to the correct collector data.
128
129 // However since this is the geometry lets assume it is fixed for now.
130 const auto exprun = getRunList()[0];
131 B2INFO("ExpRun used for DB Geometry : " << exprun.first << " " << exprun.second);
132 updateDBObjPtrs(1, exprun.second, exprun.first);
133
134 // CDCGeometryPar basically constructs a ton of objects and other DB objects.
135 // Normally we'd call updateDBObjPtrs to set the values of the requested DB objects.
136 // But in CDCGeometryPar the DB objects get used during the constructor so they must
137 // be set before/during the constructor.
138
139 // Since we are avoiding using the DataStore EventMetaData, we need to pass in
140 // an EventMetaData object to be used when constructing the DB objects.
141
142 B2INFO("Creating CDCGeometryPar object");
144
145 auto hEvtT0 = getObjectPtr<TH1F>("hEventT0");
146 double dEventT0 = hEvtT0->GetMean();
147 createHisto();
148 TH1F* hm_All = new TH1F("hm_All", "mean of #DeltaT distribution for all chanels", 500, -10, 10);
149 TH1F* hs_All = new TH1F("hs_All", "#sigma of #DeltaT distribution for all chanels", 100, 0, 10);
151
152 TF1* g1 = new TF1("g1", "gaus", -100, 100);
153 g1->SetParLimits(1, -20, 20);
154 vector<double> b, db, Sb, dSb;
155 vector<double> c[56];
156 vector<double> dc[56];
157 vector<double> s[56];
158 vector<double> ds[56];
159
160 B2INFO("Gaus fitting for whole channel");
161 double par[3];
162 m_hTotal->SetDirectory(0);
163 double mean = m_hTotal->GetMean();
164 m_hTotal->Fit("g1", "Q", "", mean - 15, mean + 15);
165 g1->GetParameters(par);
166
167 B2INFO("Gaus fitting for each board");
168 for (int ib = 1; ib < 300; ++ib) {
169 // Set Delta_T0=0 again to make sure there is no strange case
170 // in which T0 might be initialed with a strange value
171 dtb[ib] = 0;
172 err_dtb[ib] = 0;
173
174 if (m_hT0b[ib]->Integral(1, m_hT0b[ib]->GetNbinsX()) < 50) {
175 //set error to large number as a flag of bad value
176 err_dtb[ib] = 50; continue;
177 }
178 mean = m_hT0b[ib]->GetMean();
179 m_hT0b[ib]->SetDirectory(0);
180 m_hT0b[ib]->Fit("g1", "Q", "", mean - 15, mean + 15);
181 g1->GetParameters(par);
182 b.push_back(ib);
183 db.push_back(0.0);
184 Sb.push_back(par[1]);
185 dSb.push_back(g1->GetParError(1));
186 dtb[ib] = par[1] + dEventT0;// add dEvtT0 here
187 err_dtb[ib] = g1->GetParError(1);
188 }
189 B2INFO("Gaus fitting for each cell");
190 for (int ilay = 0; ilay < 56; ++ilay) {
191 const unsigned int nW = cdcgeo.nWiresInLayer(ilay);
192 for (unsigned int iwire = 0; iwire < nW; ++iwire) {
193
194 const int n = m_h1[ilay][iwire]->Integral(1, m_h1[ilay][iwire]->GetNbinsX()) ;
195 B2DEBUG(21, "layer " << ilay << " wire " << iwire << " entries " << n);
196 // Set Delta_T0=0 again to make sure there is no strange case
197 // in which T0 might be initialed with a strange value
198 dt[ilay][iwire] = 0;
199 err_dt[ilay][iwire] = 0;
200
201 if (n < 30) {
202 //set error to large number as a flag of bad value
203 err_dt[ilay][iwire] = 50; continue;
204 }
205
206 mean = m_h1[ilay][iwire]->GetMean();
207 m_h1[ilay][iwire]->SetDirectory(0);
208 m_h1[ilay][iwire]->Fit("g1", "Q", "", mean - 15, mean + 15);
209 g1->GetParameters(par);
210 c[ilay].push_back(iwire);
211 dc[ilay].push_back(0.0);
212 s[ilay].push_back(par[1]);
213 ds[ilay].push_back(g1->GetParError(1));
214
215 dt[ilay][iwire] = par[1] + dEventT0; // add dEvtT0 here;
216 err_dt[ilay][iwire] = g1->GetParError(1);
217 hm_All->Fill(par[1]);// mean of gauss fitting.
218 hs_All->Fill(par[2]); // sigma of gauss fitting.
219 }
220 }
221
222 if (m_storeHisto) {
223 B2INFO("Storing histograms");
224 auto hNDF = getObjectPtr<TH1F>("hNDF");
225 auto hPval = getObjectPtr<TH1F>("hPval");
226 TFile* fout = new TFile(m_histName.c_str(), "RECREATE");
227 fout->cd();
228 TGraphErrors* gr[56];
229 TDirectory* top = gDirectory;
230
231 //store NDF, P-val. EventT0 histogram for monitoring during calibration
232 if (hNDF && hPval && hEvtT0) {
233 hEvtT0->Write();
234 hPval->Write();
235 hNDF->Write();
236 }
237 m_hTotal->Write();
238 hm_All->Write();
239 hs_All->Write();
240 TDirectory* subDir[56];
241 for (int ilay = 0; ilay < 56; ++ilay) {
242 subDir[ilay] = top ->mkdir(Form("lay_%d", ilay));
243 subDir[ilay]->cd();
244 const unsigned int nW = cdcgeo.nWiresInLayer(ilay);
245 for (unsigned int iwire = 0; iwire < nW; ++iwire) {
246 m_h1[ilay][iwire]->Write();
247 }
248 }
249
250 top->cd();
251 TDirectory* corrT0 = top->mkdir("DeltaT0");
252 corrT0->cd();
253
254 if (b.size() > 2) {
255 TGraphErrors* grb = new TGraphErrors(b.size(), &b.at(0), &Sb.at(0), &db.at(0), &dSb.at(0));
256 grb->SetMarkerColor(2);
257 grb->SetMarkerSize(1.0);
258 grb->SetTitle("#DeltaT0;BoardID;#DeltaT0[ns]");
259 grb->SetMaximum(10);
260 grb->SetMinimum(-10);
261 grb->SetName("Board");
262 grb->Write();
263 }
264 for (int sl = 0; sl < 56; ++sl) {
265 if (c[sl].size() < 2) continue;
266 gr[sl] = new TGraphErrors(c[sl].size(), &c[sl].at(0), &s[sl].at(0), &dc[sl].at(0), &ds[sl].at(0));
267 gr[sl]->SetMarkerColor(2);
268 gr[sl]->SetMarkerSize(1.0);
269 gr[sl]->SetTitle(Form("Layer_%d;IWire;#LT t_{mea}-t_{fit} #GT [ns]", sl));
270 gr[sl]->SetMaximum(10);
271 gr[sl]->SetMinimum(-10);
272 gr[sl]->SetName(Form("lay%d", sl));
273 gr[sl]->Write();
274 }
275 fout->Close();
276 }
277 B2INFO("Writing constants...");
278 int Ngood_T0 = write();
279 B2INFO("Checking conversion conditons...");
280 double rms_dT = hm_All->GetRMS();
281 double n_below = hm_All->Integral(0, hm_All->GetXaxis()->FindBin(m_offsetMeanDt - 5 * rms_dT));
282 double n_upper = hm_All->Integral(hm_All->GetXaxis()->FindBin(m_offsetMeanDt + 5 * rms_dT), hm_All->GetXaxis()->GetNbins() - 1);
283 int N_remaining = n_below + n_upper - (14112 - Ngood_T0);
284 if (N_remaining < 0) N_remaining = 0;
285 B2INFO("+ Number of channel which are properly calibrated : " << Ngood_T0);
286 B2INFO("+ Number of channel which still need to be calibrated are: " << N_remaining << "(requirement <" << m_maxBadChannel << ")");
287 B2INFO("+ Median of Delta_T - offset:" << hm_All->GetMean() - m_offsetMeanDt << "(requirement: <" << m_maxMeanDt << ")");
288 B2INFO("+ RMS of Delta_T dist. :" << rms_dT << " (Requirement: <" << m_maxRMSDt << ")");
289
290 if (fabs(hm_All->GetMean() - m_offsetMeanDt) < m_maxMeanDt
291 && fabs(hm_All->GetRMS()) < m_maxRMSDt
292 && N_remaining < m_maxBadChannel) {
293 B2INFO("T0 Calibration Finished:");
294 return c_OK;
295 } else {
296 B2INFO("Need more iteration ...");
297 return c_Iterate;
298 }
299}
300
302{
304 CDCTimeZeros* tz = new CDCTimeZeros();
305
306 TH1F* hT0B[300];
307 for (int ib = 0; ib < 300; ++ib) {
308 hT0B[ib] = new TH1F(Form("hT0B%d", ib), Form("boardID_%d", ib), 9000, 0, 9000);
309 }
310
311 TH1F* hT0_all = new TH1F("hT0_all", "T0 distribution", 9000, 0, 9000);
312 double T0;
313 // get old T0 to calculate T0 mean of each board
314 for (int ilay = 0; ilay < 56; ++ilay) {
315 const unsigned int nW = cdcgeo.nWiresInLayer(ilay);
316 for (unsigned int iwire = 0; iwire < nW; ++iwire) {
317 WireID wireid(ilay, iwire);
318 T0 = cdcgeo.getT0(wireid);
319 hT0_all->Fill(T0);
320 hT0B[cdcgeo.getBoardID(wireid)]->Fill(T0);
321 }
322 }
323
324 //get Nominal T0
325 double T0_average = getMeanT0(hT0_all);
326 if (fabs(T0_average - m_commonT0) > 50) {B2WARNING("Large difference between common T0 (" << m_commonT0 << ") and aveage value" << T0_average);}
327 // get average T0 for each board and also apply T0 correction for T0 of that board
328 double T0B[300];
329 for (int i = 0; i < 300; ++i) {T0B[i] = m_commonT0;}
330
331
332 for (int ib = 1; ib < 300; ++ib) {
333 T0B[ib] = getMeanT0(hT0B[ib]);
334 if (fabs(T0B[ib] - T0_average) > 25) {
335 B2WARNING("T0 of Board " << ib << " (= " << T0B[ib] << ") is too different with common T0: " << T0_average <<
336 "\n It will be replaced by common T0");
337 T0B[ib] = T0_average;
338 continue;
339 }
340 //correct T0 board
341 if (abs(err_dtb[ib]) < 2 && abs(dtb[ib]) < 20) {
342 T0B[ib] -= dtb[ib];
343 }
344 }
345
346 //correct T0 and write
347 double dT;
348 int N_good = 0;
349 for (int ilay = 0; ilay < 56; ++ilay) {
350 const unsigned int nW = cdcgeo.nWiresInLayer(ilay);
351 for (unsigned int iwire = 0; iwire < nW; ++iwire) {
352 WireID wireid(ilay, iwire);
353 int bID = cdcgeo.getBoardID(wireid);
354
355 //get old T0, replace with common T0 of board or average over all channel.
356 T0 = cdcgeo.getT0(wireid);
357 if (fabs(T0 - T0B[bID]) > 25 || fabs(T0 - T0_average) > 25) {
358 B2WARNING("T0 of wireID L-W: " << ilay << "--" << iwire << " (= " << T0
359 << ") is too different with common T0 of board: " << bID << " = " << T0B[bID] <<
360 "\n It will be replaced by common T0 of the Board");
361 T0 = T0B[bID];
362 }
363
364 //select DeltaT
365 dT = 0;
366 if (abs(err_dt[ilay][iwire]) < 2 && abs(dt[ilay][iwire]) < 20) {
367 dT = dt[ilay][iwire];
368 N_good += 1;
369 }
370
371 //export
372 tz->setT0(wireid, T0 - dT);
373 }
374 }
375
376 if (m_textOutput == true) {
378 }
379 saveCalibration(tz, "CDCTimeZeros");
380 return N_good;
381}
383{
384 double mean1 = h1->GetMean();
385 h1->GetXaxis()->SetRangeUser(mean1 - 50, 8000);
386 double mean2 = h1->GetMean();
387 while (fabs(mean1 - mean2) > 0.5) {
388 mean1 = mean2;
389 h1->GetXaxis()->SetRangeUser(mean1 - 50, 8000);
390 mean2 = h1->GetMean();
391 }
392 return mean2;
393
394}
double R
typedef autogenerated by FFTW
Database object for timing offset (t0).
Definition: CDCTimeZeros.h:26
void outputToFile(std::string fileName) const
Output the contents in text file format.
Definition: CDCTimeZeros.h:123
void setT0(unsigned short iCLayer, unsigned short iWire, float t0)
Set t0 in the list.
Definition: CDCTimeZeros.h:40
The Class for CDC Geometry Parameters.
unsigned short getBoardID(const WireID &wID) const
Returns frontend board id. corresponding to the wire id.
float getT0(const WireID &wireID) const
Returns t0 parameter of the specified sense wire.
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
double senseWireR(int layerId) const
Returns radius of sense wire in each layer.
std::string m_outputT0FileName
output t0 file name for text file
double m_offsetMeanDt
offset dT distribution caused by event timing extraction;
TH1F * m_hTotal
1D histogram of delta T whole channel
double dt[56][385]
dt of each channel
void createHisto()
create histo for each channel
double getMeanT0(TH1F *h1)
calculate mean of the T0 distribution
double m_maxRMSDt
RMS of dT distribution of all channels.
DBObjPtr< CDCGeometry > m_cdcGeo
Geometry of CDC.
double m_commonT0
A common T0 of all channels.
double m_maxMeanDt
Mean of dT distribution of all channels;.
TH1F * m_h1[56][385]
1D histogram for each channel
double err_dt[56][385]
error of dt of each channel
bool m_textOutput
output text file if true
EResult calibrate() override
Run algo on data.
TH1F * m_hT0b[300]
1D histogram for each board
double m_Pvalmin
minimum pvalue required
double err_dtb[300]
error of dt of board
double m_maxBadChannel
Number of channels which has DeltaT0 larger then 0.5.
Base class for calibration algorithms.
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
void updateDBObjPtrs(const unsigned int event, const int run, const int experiment)
Updates any DBObjPtrs by calling update(event) for DBStore.
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.
EResult
The result of calibration.
@ c_OK
Finished successfully =0 in Python.
@ c_Iterate
Needs iteration =1 in Python.
Class to identify a wire inside the CDC.
Definition: WireID.h:34
Abstract base class for different kinds of events.
STL namespace.