Belle II Software development
SpaceResolutionCalibrationAlgorithm.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 <cdc/calibration/SpaceResolutionCalibrationAlgorithm.h>
10#include <cdc/dbobjects/CDCSpaceResols.h>
11
12#include <framework/database/DBObjPtr.h>
13#include <framework/logging/Logger.h>
14
15#include <TError.h>
16#include <TF1.h>
17#include <TH1F.h>
18#include <TH2F.h>
19#include <TMinuit.h>
20#include <TROOT.h>
21#include <TStopwatch.h>
22
23using namespace std;
24using namespace Belle2;
25using namespace CDC;
26
27typedef std::array<float, 3> array3;
28SpaceResolutionCalibrationAlgorithm::SpaceResolutionCalibrationAlgorithm() :
29 CalibrationAlgorithm("CDCCalibrationCollector")
30{
31 setDescription(
32 " -------------------------- Position Resolution Calibration Algorithm -------------------------\n"
33 );
34}
35void SpaceResolutionCalibrationAlgorithm::createHisto()
36{
37 B2INFO("Creating histograms");
38 const int np = floor(1 / m_binWidth);
39
40 vector<double> yu;
41 vector <double> yb;
42 for (int i = 0; i < 50; ++i) {
43 yb.push_back(-0.07 + i * (0.14 / 50));
44 }
45 for (int i = 0; i < 50; ++i) {
46 yu.push_back(-0.08 + i * (0.16 / 50));
47 }
48
49 vector<double> xbin;
50 xbin.push_back(0.);
51 xbin.push_back(0.02);
52 for (int i = 1; i < np; ++i) {
53 xbin.push_back(i * m_binWidth);
54 }
55
56 for (int il = 0; il < 56; ++il) {
57 for (int lr = 0; lr < 2; ++lr) {
58 for (int al = 0; al < m_nAlphaBins; ++al) {
59 for (int th = 0; th < m_nThetaBins; ++th) {
60 m_hBiased[il][lr][al][th] = new TH2F(Form("hb_%d_%d_%d_%d", il, lr, al, th),
61 Form("lay_%d_lr%d_al_%3.0f_th_%3.0f;Drift Length [cm];#DeltaX", il, lr, m_iAlpha[al], m_iTheta[th]),
62 xbin.size() - 1, &xbin.at(0), yb.size() - 1, &yb.at(0));
63 m_hUnbiased[il][lr][al][th] = new TH2F(Form("hu_%d_%d_%d_%d", il, lr, al, th),
64 Form("lay_%d_lr%d_al_%3.0f_th_%3.0f;Drift Length [cm];#DeltaX", il, lr, m_iAlpha[al], m_iTheta[th]),
65 xbin.size() - 1, &xbin.at(0), yu.size() - 1, &yu.at(0));
66 }
67 }
68 }
69 }
70
71
72 auto tree = getObjectPtr<TTree>("tree");
73
74 UChar_t lay;
75 Float_t w;
76 Float_t x_u;
77 Float_t x_b;
78 Float_t x_mea;
79 Float_t Pval;
80 Float_t alpha;
81 Float_t theta;
82 Float_t ndf;
83 Float_t absRes_u;
84 Float_t absRes_b;
85 tree->SetBranchAddress("lay", &lay);
86 tree->SetBranchAddress("ndf", &ndf);
87 tree->SetBranchAddress("Pval", &Pval);
88 tree->SetBranchAddress("x_u", &x_u);
89 tree->SetBranchAddress("x_b", &x_b);
90 tree->SetBranchAddress("x_mea", &x_mea);
91 tree->SetBranchAddress("weight", &w);
92 tree->SetBranchAddress("alpha", &alpha);
93 tree->SetBranchAddress("theta", &theta);
94
95 /* Disable unused branch */
96 std::vector<TString> list_vars = {"lay", "ndf", "Pval", "x_u", "x_b", "x_mea", "weight", "alpha", "theta"};
97 tree->SetBranchStatus("*", 0);
98
99 for (TString brname : list_vars) {
100 tree->SetBranchStatus(brname, 1);
101 }
102
103
104 const Long64_t nEntries = tree->GetEntries();
105 B2INFO("Number of entries: " << nEntries);
106 int ith = -99;
107 int ial = -99;
108 TStopwatch timer;
109 timer.Start();
110 for (Long64_t i = 0; i < nEntries; ++i) {
111 tree->GetEntry(i);
112 if (std::fabs(x_b) < 0.02 || std::fabs(x_u) < 0.02) continue;
113 if (Pval < m_minPval || ndf < m_minNdf) continue;
114
115 for (int k = 0; k < m_nAlphaBins; ++k) {
116 if (alpha < m_upperAlpha[k]) {
117 ial = k;
118 break;
119 }
120 }
121
122 for (int j = 0; j < m_nThetaBins; ++j) {
123 if (theta < m_upperTheta[j]) {
124 ith = j;
125 break;
126 }
127 }
128
129 int ilr = x_u > 0 ? 1 : 0;
130
131 if (ial == -99 || ith == -99) {
132 TString command = Form("Error in alpha=%3.2f and theta = %3.2f>> error", alpha, theta);
133 B2FATAL("ERROR" << command);
134 }
135
136 absRes_u = fabs(x_mea) - fabs(x_u);
137 absRes_b = fabs(x_mea) - fabs(x_b);
138
139 int ilay = static_cast<int>(lay);
140 m_hUnbiased[ilay][ilr][ial][ith]->Fill(fabs(x_u), absRes_u, w);
141 m_hBiased[ilay][ilr][ial][ith]->Fill(fabs(x_b), absRes_b, w);
142 }
143
144 timer.Stop();
145 B2INFO("Time to fill histograms: " << timer.RealTime() << "s");
146
147 B2INFO("Start to obtain the biased and unbiased sigmas...");
148 TF1* gb = new TF1("gb", "gaus", -0.05, 0.05);
149 TF1* gu = new TF1("gu", "gaus", -0.06, 0.06);
150 TF1* g0b = new TF1("g0b", "gaus", -0.015, 0.07);
151 TF1* g0u = new TF1("g0u", "gaus", -0.015, 0.08);
152
153 std::vector<double> sigma;
154 std::vector<double> dsigma;
155 std::vector<double> s2;
156 std::vector<double> ds2;
157 std::vector<double> xl;
158 std::vector<double> dxl;
159 std::vector<double> dxl0;
160
161 ofstream ofss("IntReso.dat");
162 const int ib1 = int(0.1 / m_binWidth) + 1;
163 int firstbin = 1;
164 int minEntry = 10;
165 for (int il = 0; il < 56; ++il) {
166 for (int lr = 0; lr < 2; ++lr) {
167 for (int al = 0; al < m_nAlphaBins; ++al) {
168 for (int th = 0; th < m_nThetaBins; ++th) {
169
170 B2DEBUG(21, "layer-lr-al-th " << il << " - " << lr << " - " << al << " - " << th);
171 if (m_hBiased[il][lr][al][th]->GetEntries() < 5000) {
172 m_fitStatus[il][lr][al][th] = -1;
173 continue;
174 }
175
176 auto* proYb = m_hBiased[il][lr][al][th]->ProjectionY();
177 auto* proYu = m_hUnbiased[il][lr][al][th]->ProjectionY();
178
179 g0b->SetParLimits(0, 0, m_hBiased[il][lr][al][th]->GetEntries() * 5);
180 g0u->SetParLimits(0, 0, m_hUnbiased[il][lr][al][th]->GetEntries() * 5);
181 g0b->SetParLimits(1, -0.01, 0.004);
182 g0u->SetParLimits(1, -0.01, 0.004);
183 g0b->SetParLimits(2, 0.0, proYb->GetRMS() * 5);
184 g0u->SetParLimits(2, 0.0, proYu->GetRMS() * 5);
185
186 g0b->SetParameter(0, m_hBiased[il][lr][al][th]->GetEntries());
187 g0u->SetParameter(0, m_hUnbiased[il][lr][al][th]->GetEntries());
188 g0b->SetParameter(1, 0);
189 g0u->SetParameter(1, 0);
190 g0b->SetParameter(2, proYb->GetRMS());
191 g0u->SetParameter(2, proYu->GetRMS());
192
193 B2DEBUG(21, "Nentries: " << m_hBiased[il][lr][al][th]->GetEntries());
194 m_hBiased[il][lr][al][th]->SetDirectory(0);
195 m_hUnbiased[il][lr][al][th]->SetDirectory(0);
196
197 // With biased track fit result
198
199 // Apply slice fit for the region near sense wire
200 m_hBiased[il][lr][al][th]->FitSlicesY(g0b, firstbin, ib1, minEntry);
201
202 // mean
203 m_hMeanBiased[il][lr][al][th] = static_cast<TH1F*>(gDirectory->Get(Form("hb_%d_%d_%d_%d_1", il, lr, al,
204 th))->Clone(Form("hb_%d_%d_%d_%d_m", il,
205 lr, al,
206 th)));
207 // sigma
208 m_hSigmaBiased[il][lr][al][th] = static_cast<TH1F*>(gDirectory->Get(Form("hb_%d_%d_%d_%d_2", il, lr, al,
209 th))->Clone(Form("hb_%d_%d_%d_%d_s",
210 il, lr, al,
211 th)));
212 m_hMeanBiased[il][lr][al][th]->SetDirectory(0);
213 m_hSigmaBiased[il][lr][al][th]->SetDirectory(0);
214
215 //Apply slice fit for other regions
216 m_hBiased[il][lr][al][th]->FitSlicesY(gb, ib1 + 1, np, minEntry);
217 // mean
218 m_hMeanBiased[il][lr][al][th]->Add(static_cast<TH1F*>(gDirectory->Get(Form("hb_%d_%d_%d_%d_1", il, lr, al, th))));
219 //sigma
220 m_hSigmaBiased[il][lr][al][th]->Add(static_cast<TH1F*>(gDirectory->Get(Form("hb_%d_%d_%d_%d_2", il, lr, al, th))));
221 B2DEBUG(21, "entries (2nd): " << m_hSigmaBiased[il][lr][al][th]->GetEntries());
222
223 // With unbiased track fit result
224
225 // Apply slice fit for the region near sense wire
226 m_hUnbiased[il][lr][al][th]->FitSlicesY(g0u, firstbin, ib1, minEntry);
227 // mean
228 m_hMeanUnbiased[il][lr][al][th] = static_cast<TH1F*>(gDirectory->Get(Form("hu_%d_%d_%d_%d_1", il, lr, al,
229 th))->Clone(Form("hu_%d_%d_%d_%d_m",
230 il, lr, al,
231 th)));
232 // sigma
233 m_hSigmaUnbiased[il][lr][al][th] = static_cast<TH1F*>(gDirectory->Get(Form("hu_%d_%d_%d_%d_2", il, lr, al,
234 th))->Clone(Form("hu_%d_%d_%d_%d_s",
235 il, lr, al,
236 th)));
237 m_hMeanUnbiased[il][lr][al][th]->SetDirectory(0);
238 m_hSigmaUnbiased[il][lr][al][th]->SetDirectory(0);
239
240
241 //Apply slice fit for other regions
242 m_hUnbiased[il][lr][al][th]->FitSlicesY(gu, ib1 + 1, np, minEntry);
243 //mean
244 m_hMeanUnbiased[il][lr][al][th]->Add(static_cast<TH1F*>(gDirectory->Get(Form("hu_%d_%d_%d_%d_1", il, lr, al, th))));
245 //sigma
246 m_hSigmaUnbiased[il][lr][al][th]->Add(static_cast<TH1F*>(gDirectory->Get(Form("hu_%d_%d_%d_%d_2", il, lr, al, th))));
247 if (!m_hSigmaUnbiased[il][lr][al][th] || !m_hSigmaBiased[il][lr][al][th]) {
248 B2WARNING("sliced histo not found");
249 m_fitStatus[il][lr][al][th] = -1;
250 continue;
251 }
252 //clean up container before adding new values.
253 xl.clear();
254 dxl.clear();
255 dxl0.clear();
256 sigma.clear();
257 dsigma.clear();
258 s2.clear();
259 ds2.clear();
260
261
262 for (int j = 1; j < m_hSigmaUnbiased[il][lr][al][th]->GetNbinsX(); j++) {
263 if (m_hSigmaUnbiased[il][lr][al][th]->GetBinContent(j) == 0) continue;
264 if (m_hSigmaBiased[il][lr][al][th]->GetBinContent(j) == 0) continue;
265 double sb = m_hSigmaBiased[il][lr][al][th]->GetBinContent(j);
266 double su = m_hSigmaUnbiased[il][lr][al][th]->GetBinContent(j);
267
268 double dsb = m_hSigmaBiased[il][lr][al][th]->GetBinError(j);
269 double dsu = m_hSigmaUnbiased[il][lr][al][th]->GetBinError(j);
270 double XL = m_hSigmaBiased[il][lr][al][th]->GetXaxis()->GetBinCenter(j);
271 double dXL = (m_hSigmaBiased[il][lr][al][th]->GetXaxis()->GetBinWidth(j)) / 2;
272 double s_int = std::sqrt(sb * su);
273 double ds_int = 0.5 * s_int * (dsb / sb + dsu / su);
274 if (ds_int > 0.02) continue;
275 xl.push_back(XL);
276 dxl.push_back(dXL);
277 dxl0.push_back(0.);
278 sigma.push_back(s_int);
279 dsigma.push_back(ds_int);
280 s2.push_back(s_int * s_int);
281 ds2.push_back(2 * s_int * ds_int);
282 ofss << il << " " << lr << " " << al << " " << th << " " << j << " " << XL << " " << dXL << " " << s_int << " " <<
283 ds_int << endl;
284 }
285
286 if (xl.size() < 7 || xl.size() > Max_np) {
287 m_fitStatus[il][lr][al][th] = -1;
288 B2WARNING("number of element might out of range"); continue;
289 }
290
291 //Intrinsic resolution
292 B2DEBUG(21, "Create Histo for layer-lr: " << il << " " << lr);
293 m_graph[il][lr][al][th] = new TGraphErrors(xl.size(), &xl.at(0), &sigma.at(0), &dxl.at(0), &dsigma.at(0));
294 m_graph[il][lr][al][th]->SetMarkerSize(0.5);
295 m_graph[il][lr][al][th]->SetMarkerStyle(8);
296 m_graph[il][lr][al][th]->SetTitle(Form("Layer_%d lr%d #alpha = %3.0f #theta = %3.0f", il, lr, m_iAlpha[al], m_iTheta[th]));
297 m_graph[il][lr][al][th]->SetName(Form("lay%d_lr%d_al%d_th%d", il, lr, al, th));
298
299 //square of sigma for fitting
300 m_gFit[il][lr][al][th] = new TGraphErrors(xl.size(), &xl.at(0), &s2.at(0), &dxl0.at(0), &ds2.at(0));
301 m_gFit[il][lr][al][th]->SetMarkerSize(0.5);
302 m_gFit[il][lr][al][th]->SetMarkerStyle(8);
303 m_gFit[il][lr][al][th]->SetTitle(Form("L%d lr%d #alpha = %3.0f #theta = %3.0f ", il, lr, m_iAlpha[al], m_iTheta[th]));
304 m_gFit[il][lr][al][th]->SetName(Form("sigma2_lay%d_lr%d_al%d_th%d", il, lr, al, th));
305
306 gDirectory->Delete("hu_%d_%d_%d_%d_0");
307 }
308 }
309 }
310 }
311 ofss.close();
312
313}
314
315CalibrationAlgorithm::EResult SpaceResolutionCalibrationAlgorithm::calibrate()
316{
317
318 B2INFO("Start calibration");
319 gPrintViaErrorHandler = true; // Suppress huge log output from TMinuit
320 gROOT->SetBatch(1);
321 gErrorIgnoreLevel = 3001;
322
323 const auto exprun = getRunList()[0];
324 B2INFO("ExpRun used for DB Geometry : " << exprun.first << " " << exprun.second);
325 updateDBObjPtrs(1, exprun.second, exprun.first);
326 // B2INFO("Creating CDCGeometryPar object");
327 // CDC::CDCGeometryPar::Instance(&(*m_cdcGeo));
328
329 prepare();
330 createHisto();
331
332 TF1* func = new TF1("func", "[0]/(x*x + [1])+[2]* x+[3]+[4]*exp([5]*(x-[6])*(x-[6]))", 0, 1.);
333 TH1F* hprob = new TH1F("h1", "", 20, 0, 1);
334 double upFit;
335 double intp6;
336
337 for (int i = 0; i < 56; ++i) {
338 for (int lr = 0; lr < 2; ++lr) {
339 for (int al = 0; al < m_nAlphaBins; ++al) {
340 for (int th = 0; th < m_nThetaBins; ++th) {
341 if (!m_gFit[i][lr][al][th]) continue;
342 if (m_fitStatus[i][lr][al][th] != -1) { /*if graph exist, do fitting*/
343 upFit = getUpperBoundaryForFit(m_gFit[i][lr][al][th]);
344 intp6 = upFit + 0.2;
345 B2DEBUG(199, "xmax for fitting: " << upFit);
346
347 func->SetParameters(5E-6, 0.007, 1E-4, 1E-5, 0.00008, -30, intp6);
348 func->SetParLimits(0, 1E-7, 1E-4);
349 func->SetParLimits(1, 0.0045, 0.02);
350 func->SetParLimits(2, 1E-6, 0.0005);
351 func->SetParLimits(3, 1E-8, 0.0005);
352 func->SetParLimits(4, 0., 0.001);
353 func->SetParLimits(5, -40, 0.);
354 func->SetParLimits(6, intp6 - 0.5, intp6 + 0.2);
355
356 B2DEBUG(21, "Fitting for layer: " << i << "lr: " << lr << " ial" << al << " ith:" << th);
357 B2DEBUG(21, "Fit status before fit:" << m_fitStatus[i][lr][al][th]);
358
359 for (int j = 0; j < 10; j++) {
360
361 B2DEBUG(21, "loop: " << j);
362 B2DEBUG(21, "Int p6: " << intp6);
363 B2DEBUG(21, "Number of Point: " << m_gFit[i][lr][al][th]->GetN());
364 Int_t stat = m_gFit[i][lr][al][th]->Fit("func", "MQE", "", 0.05, upFit);
365 B2DEBUG(21, "stat of fit" << stat);
366 std::string Fit_status = gMinuit->fCstatu.Data();
367 B2DEBUG(21, "FIT STATUS: " << Fit_status);
368 if (Fit_status == "OK" || Fit_status == "SUCCESSFUL" || Fit_status == "CALL LIMIT"
369 || Fit_status == "PROBLEMS") {//need to found better way
370 if (fabs(func->Eval(0.3)) > 0.00035 || func->Eval(0.3) < 0) {
371 func->SetParameters(5E-6, 0.007, 1E-4, 1E-7, 0.0007, -30, intp6 + 0.05 * j);
372 func->SetParLimits(6, intp6 + 0.05 * j - 0.5, intp6 + 0.05 * j + 0.2);
373 // func->SetParameters(defaultparsmall);
374 m_fitStatus[i][lr][al][th] = 0;
375 } else {
376 B2DEBUG(21, "Prob of fit: " << func->GetProb());
377 m_fitStatus[i][lr][al][th] = 1;
378 break;
379 }
380 } else {
381 m_fitStatus[i][lr][al][th] = 0;
382 func->SetParameters(5E-6, 0.007, 1E-4, 1E-7, 0.0007, -30, intp6 + 0.05 * j);
383 func->SetParLimits(6, intp6 + 0.05 * j - 0.5, intp6 + 0.05 * j + 0.2);
384 upFit += 0.025;
385 if (j == 9) {
386 // TCanvas* c1 = new TCanvas("c1", "", 600, 600);
387 // m_gFit[i][lr][al][th]->Draw();
388 // c1->SaveAs(Form("Sigma_Fit_Error_%s_%d_%d_%d_%d.png", Fit_status.c_str(), i, lr, al, th));
389 // B2WARNING("Fit error: " << i << " " << lr << " " << al << " " << th);
390 }
391 }
392 }
393 if (m_fitStatus[i][lr][al][th] == 1) {
394 B2DEBUG(21, "ProbFit: Lay_lr_al_th: " << i << " " << lr << " " << al << " " << th << func->GetProb());
395 hprob->Fill(func->GetProb());
396 func->GetParameters(m_sigma[i][lr][al][th]);
397 }
398 }
399 }
400 }
401 }
402 }
403
404 write();
405 storeHisto();
406
407 const int nTotal = 56 * 2 * m_nAlphaBins * m_nThetaBins;
408 int nFitCompleted = 0;
409 for (int l = 0; l < 56; ++l) {
410 for (int lr = 0; lr < 2; ++lr) {
411 for (int al = 0; al < m_nAlphaBins; ++al) {
412 for (int th = 0; th < m_nThetaBins; ++th) {
413 if (m_fitStatus[l][lr][al][th] == 1) {
414 nFitCompleted += 1;
415 }
416 }
417 }
418 }
419 }
420
421 if (static_cast<double>(nFitCompleted) / nTotal < m_threshold) {
422 B2WARNING("Less than " << m_threshold * 100 << " % of Sigmas were fitted.");
423 return c_NotEnoughData;
424 }
425
426 return c_OK;
427}
428
429void SpaceResolutionCalibrationAlgorithm::storeHisto()
430{
431 B2INFO("saving histograms");
432
433 TFile* ff = new TFile(m_histName.c_str(), "RECREATE");
434 TDirectory* top = gDirectory;
435 TDirectory* Direct[56];
436
437 auto hNDF = getObjectPtr<TH1F>("hNDF");
438 auto hPval = getObjectPtr<TH1F>("hPval");
439 auto hEvtT0 = getObjectPtr<TH1F>("hEventT0");
440 //store NDF, P-val. EventT0 histogram for monitoring during calibration
441 if (hNDF && hPval && hEvtT0) {
442 hEvtT0->Write();
443 hPval->Write();
444 hNDF->Write();
445 }
446
447
448 for (int il = 0; il < 56; ++il) {
449 top->cd();
450 Direct[il] = gDirectory->mkdir(Form("lay_%d", il));
451 Direct[il]->cd();
452
453 for (int lr = 0; lr < 2; ++lr) {
454 for (int al = 0; al < m_nAlphaBins; ++al) {
455 for (int th = 0; th < m_nThetaBins; ++th) {
456 if (!m_graph[il][lr][al][th]) continue;
457 if (!m_gFit[il][lr][al][th]) continue;
458 if (m_fitStatus[il][lr][al][th] == 1) {
459 m_hBiased[il][lr][al][th]->Write();
460 m_hUnbiased[il][lr][al][th]->Write();
461 m_hMeanBiased[il][lr][al][th]->Write();
462 m_hSigmaBiased[il][lr][al][th]->Write();
463 m_hMeanUnbiased[il][lr][al][th]->Write();
464 m_hSigmaUnbiased[il][lr][al][th]->Write();
465 m_graph[il][lr][al][th]->Write();
466 m_gFit[il][lr][al][th]->Write();
467 }
468 }
469 }
470 }
471 }
472 ff->Close();
473 B2INFO("Finish store histogram");
474
475}
476void SpaceResolutionCalibrationAlgorithm::write()
477{
478 B2INFO("Writing calibrated sigma's");
479 int nfitted = 0;
480 int nfailure = 0;
481
482 CDCSpaceResols* dbSigma = new CDCSpaceResols();
483
484 const float deg2rad = M_PI / 180.0;
485
486 for (unsigned short i = 0; i < m_nAlphaBins; ++i) {
487 std::array<float, 3> alpha3 = {m_lowerAlpha[i]* deg2rad,
488 m_upperAlpha[i]* deg2rad,
489 m_iAlpha[i]* deg2rad
490 };
491 dbSigma->setAlphaBin(alpha3);
492 }
493
494
495 for (unsigned short i = 0; i < m_nThetaBins; ++i) {
496 std::array<float, 3> theta3 = {m_lowerTheta[i]* deg2rad,
497 m_upperTheta[i]* deg2rad,
498 m_iTheta[i]* deg2rad
499 };
500 dbSigma->setThetaBin(theta3);
501 }
502
503 dbSigma->setSigmaParamMode(m_sigmaParamMode);
504 for (int ialpha = 0; ialpha < m_nAlphaBins; ++ialpha) {
505 for (int itheta = 0; itheta < m_nThetaBins; ++itheta) {
506 for (int iCL = 0; iCL < 56; ++iCL) {
507 for (int iLR = 1; iLR >= 0; --iLR) {
508 std::vector<float> sgbuff;
509 if (m_fitStatus[iCL][iLR][ialpha][itheta] == 1) {
510 nfitted += 1; // inclement number of successfully fitted sigma's
511 for (int i = 0; i < 7; ++i) {
512 sgbuff.push_back(m_sigma[iCL][iLR][ialpha][itheta][i]);
513 }
514 } else {
515 //B2WARNING("Fitting error and old sigma will be used. (Layer " << iCL << ") (lr = " << iLR <<
516 // ") (al = " << ialpha << ") (th = " << itheta << ")");
517 nfailure += 1; // inclement number of fit failed sigma's
518 for (int i = 0; i < 7; ++i) {
519 sgbuff.push_back(m_sigmaPost[iCL][iLR][ialpha][itheta][i]);
520 }
521 }
522 dbSigma->setSigmaParams(iCL, iLR, ialpha, itheta, sgbuff);
523 }
524 }
525 }
526 }
527
528 if (m_textOutput == true) {
529 dbSigma->outputToFile(m_outputFileName);
530 }
531
532 saveCalibration(dbSigma, "CDCSpaceResols");
533
534 B2RESULT("Number of histogram: " << 56 * 2 * m_nAlphaBins * m_nThetaBins);
535 B2RESULT("Histos successfully fitted: " << nfitted);
536 B2RESULT("Histos fit failure: " << nfailure);
537
538
539}
540
541void SpaceResolutionCalibrationAlgorithm::prepare()
542{
543 B2INFO("Prepare calibration of space resolution");
544
545 const double rad2deg = 180 / M_PI;
546
547 DBObjPtr<CDCSpaceResols> dbSigma;
548
549 m_nAlphaBins = dbSigma->getNoOfAlphaBins();
550 m_nThetaBins = dbSigma->getNoOfThetaBins();
551
552 B2INFO("Number of alpha bins: " << m_nAlphaBins);
553 for (int i = 0; i < m_nAlphaBins; ++i) {
554 array3 alpha = dbSigma->getAlphaBin(i);
555 m_lowerAlpha[i] = alpha[0] * rad2deg;
556 m_upperAlpha[i] = alpha[1] * rad2deg;
557 m_iAlpha[i] = alpha[2] * rad2deg;
558 }
559
560 B2INFO("Number of theta bins: " << m_nThetaBins);
561 for (int i = 0; i < m_nThetaBins; ++i) {
562 array3 theta = dbSigma->getThetaBin(i);
563 m_lowerTheta[i] = theta[0] * rad2deg;
564 m_upperTheta[i] = theta[1] * rad2deg;
565 m_iTheta[i] = theta[2] * rad2deg;
566 }
567 m_sigmaParamModePost = dbSigma->getSigmaParamMode();
568
569 for (unsigned short iCL = 0; iCL < 56; ++iCL) {
570 for (unsigned short iLR = 0; iLR < 2; ++iLR) {
571 for (unsigned short iA = 0; iA < m_nAlphaBins; ++iA) {
572 for (unsigned short iT = 0; iT < m_nThetaBins; ++iT) {
573 const std::vector<float> params = dbSigma->getSigmaParams(iCL, iLR, iA, iT);
574 unsigned short np = params.size();
575 // std::cout <<"np4sigma= " << np << std::endl;
576 for (unsigned short i = 0; i < np; ++i) {
577 m_sigmaPost[iCL][iLR][iA][iT][i] = params[i];
578 }
579 }
580 }
581 }
582 }
583}
E
R E
internal precision of FFTW codelets
Definition DiscreteCosineTransform_31points.cc:36
Belle2::CDCSpaceResols
Database object for space resolutions.
Definition CDCSpaceResols.h:26
Belle2::CDCSpaceResols::setSigmaParams
void setSigmaParams(const SigmaID sigmaID, const std::vector< float > &params)
Set sigma parameters for the specified id.
Definition CDCSpaceResols.h:101
Belle2::CDCSpaceResols::setThetaBin
void setThetaBin(const array3 &theta)
Set theta-angle bin (rad)
Definition CDCSpaceResols.h:64
Belle2::CDCSpaceResols::outputToFile
void outputToFile(std::string fileName) const
Output the contents in text file format.
Definition CDCSpaceResols.h:351
Belle2::CDCSpaceResols::setAlphaBin
void setAlphaBin(const array3 &alpha)
Set alpha-angle bin (rad)
Definition CDCSpaceResols.h:50
Belle2::CDCSpaceResols::setSigmaParamMode
void setSigmaParamMode(unsigned short mode)
Set sigma parameterization mode.
Definition CDCSpaceResols.h:85
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_hSigmaBiased
TH1F * m_hSigmaBiased[56][2][Max_nalpha][Max_ntheta]
sigma histogram of biased residual
Definition SpaceResolutionCalibrationAlgorithm.h:124
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::prepare
void prepare()
Prepare the calibration of space resolution.
Definition SpaceResolutionCalibrationAlgorithm.cc:541
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_threshold
double m_threshold
minimal requirement for the fraction of fitted results
Definition SpaceResolutionCalibrationAlgorithm.h:117
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_sigmaParamMode
unsigned short m_sigmaParamMode
sigma mode for this calibration.
Definition SpaceResolutionCalibrationAlgorithm.h:137
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_sigmaParamModePost
unsigned short m_sigmaParamModePost
sigma mode before this calibration.
Definition SpaceResolutionCalibrationAlgorithm.h:140
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_hMeanBiased
TH1F * m_hMeanBiased[56][2][Max_nalpha][Max_ntheta]
mean histogram biased residual
Definition SpaceResolutionCalibrationAlgorithm.h:123
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_sigmaPost
double m_sigmaPost[56][2][18][7][8]
sigma parameters before calibration
Definition SpaceResolutionCalibrationAlgorithm.h:139
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::getUpperBoundaryForFit
double getUpperBoundaryForFit(TGraphErrors *graph)
search max point at boundary region
Definition SpaceResolutionCalibrationAlgorithm.h:81
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_hBiased
TH2F * m_hBiased[56][2][Max_nalpha][Max_ntheta]
2D histogram of biased residual
Definition SpaceResolutionCalibrationAlgorithm.h:121
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::Max_np
static const unsigned short Max_np
Maximum number of point =1/binwidth.
Definition SpaceResolutionCalibrationAlgorithm.h:109
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_hSigmaUnbiased
TH1F * m_hSigmaUnbiased[56][2][Max_nalpha][Max_ntheta]
sigma histogram of ubiased residual
Definition SpaceResolutionCalibrationAlgorithm.h:126
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_hMeanUnbiased
TH1F * m_hMeanUnbiased[56][2][Max_nalpha][Max_ntheta]
mean histogram of unbiased residual
Definition SpaceResolutionCalibrationAlgorithm.h:125
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_gFit
TGraphErrors * m_gFit[56][2][18][7]
sigma*sigma graph for fit
Definition SpaceResolutionCalibrationAlgorithm.h:119
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_fitStatus
int m_fitStatus[56][2][Max_nalpha][Max_ntheta]
Fit flag; 1:OK ; 0:error.
Definition SpaceResolutionCalibrationAlgorithm.h:127
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_hUnbiased
TH2F * m_hUnbiased[56][2][Max_nalpha][Max_ntheta]
2D histogram of unbiased residual
Definition SpaceResolutionCalibrationAlgorithm.h:122
Belle2::CalibrationAlgorithm::saveCalibration
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
Definition CalibrationAlgorithm.cc:297
Belle2::CalibrationAlgorithm::updateDBObjPtrs
static void updateDBObjPtrs(const unsigned int event, const int run, const int experiment)
Updates any DBObjPtrs by calling update(event) for DBStore.
Definition CalibrationAlgorithm.cc:405
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition CalibrationAlgorithm.h:321
Belle2::CalibrationAlgorithm::getRunList
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
Definition CalibrationAlgorithm.h:266
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition CalibrationAlgorithm.h:40
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfully =0 in Python.
Definition CalibrationAlgorithm.h:41
Belle2::CalibrationAlgorithm::c_NotEnoughData
@ c_NotEnoughData
Needs more data =2 in Python.
Definition CalibrationAlgorithm.h:43
Belle2::CalibrationAlgorithm::CalibrationAlgorithm
CalibrationAlgorithm(const std::string &collectorModuleName)
Constructor - sets the prefix for collected objects (won't be accesses until execute(....
Definition CalibrationAlgorithm.h:140
Belle2::DBObjPtr
Class for accessing objects in the database.
Definition DBObjPtr.h:21
Belle2::CalibrationAlgorithm::getObjectPtr
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...
Definition CalibrationAlgorithm.h:414
Belle2::TOP::func
commonly used functions
Definition func.h:22
Belle2
Abstract base class for different kinds of events.
Definition MillepedeAlgorithm.h:17
std
STL namespace.