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] = (TH1F*)gDirectory->Get(Form("hb_%d_%d_%d_%d_1", il, lr, al, th))->Clone(Form("hb_%d_%d_%d_%d_m", il,
204 lr, al,
205 th));
206 // sigma
207 m_hSigmaBiased[il][lr][al][th] = (TH1F*)gDirectory->Get(Form("hb_%d_%d_%d_%d_2", il, lr, al, th))->Clone(Form("hb_%d_%d_%d_%d_s",
208 il, lr, al,
209 th));
210 m_hMeanBiased[il][lr][al][th]->SetDirectory(0);
211 m_hSigmaBiased[il][lr][al][th]->SetDirectory(0);
212
213 //Apply slice fit for other regions
214 m_hBiased[il][lr][al][th]->FitSlicesY(gb, ib1 + 1, np, minEntry);
215 // mean
216 m_hMeanBiased[il][lr][al][th]->Add((TH1F*)gDirectory->Get(Form("hb_%d_%d_%d_%d_1", il, lr, al, th)));
217 //sigma
218 m_hSigmaBiased[il][lr][al][th]->Add((TH1F*)gDirectory->Get(Form("hb_%d_%d_%d_%d_2", il, lr, al, th)));
219 B2DEBUG(21, "entries (2nd): " << m_hSigmaBiased[il][lr][al][th]->GetEntries());
220
221 // With unbiased track fit result
222
223 // Apply slice fit for the region near sense wire
224 m_hUnbiased[il][lr][al][th]->FitSlicesY(g0u, firstbin, ib1, minEntry);
225 // mean
226 m_hMeanUnbiased[il][lr][al][th] = (TH1F*)gDirectory->Get(Form("hu_%d_%d_%d_%d_1", il, lr, al, th))->Clone(Form("hu_%d_%d_%d_%d_m",
227 il, lr, al,
228 th));
229 // sigma
230 m_hSigmaUnbiased[il][lr][al][th] = (TH1F*)gDirectory->Get(Form("hu_%d_%d_%d_%d_2", il, lr, al, th))->Clone(Form("hu_%d_%d_%d_%d_s",
231 il, lr, al,
232 th));
233 m_hMeanUnbiased[il][lr][al][th]->SetDirectory(0);
234 m_hSigmaUnbiased[il][lr][al][th]->SetDirectory(0);
235
236
237 //Apply slice fit for other regions
238 m_hUnbiased[il][lr][al][th]->FitSlicesY(gu, ib1 + 1, np, minEntry);
239 //mean
240 m_hMeanUnbiased[il][lr][al][th]->Add((TH1F*)gDirectory->Get(Form("hu_%d_%d_%d_%d_1", il, lr, al, th)));
241 //sigma
242 m_hSigmaUnbiased[il][lr][al][th]->Add((TH1F*)gDirectory->Get(Form("hu_%d_%d_%d_%d_2", il, lr, al, th)));
243 if (!m_hSigmaUnbiased[il][lr][al][th] || !m_hSigmaBiased[il][lr][al][th]) {
244 B2WARNING("sliced histo not found");
245 m_fitStatus[il][lr][al][th] = -1;
246 continue;
247 }
248 //clean up container before adding new values.
249 xl.clear();
250 dxl.clear();
251 dxl0.clear();
252 sigma.clear();
253 dsigma.clear();
254 s2.clear();
255 ds2.clear();
256
257
258 for (int j = 1; j < m_hSigmaUnbiased[il][lr][al][th]->GetNbinsX(); j++) {
259 if (m_hSigmaUnbiased[il][lr][al][th]->GetBinContent(j) == 0) continue;
260 if (m_hSigmaBiased[il][lr][al][th]->GetBinContent(j) == 0) continue;
261 double sb = m_hSigmaBiased[il][lr][al][th]->GetBinContent(j);
262 double su = m_hSigmaUnbiased[il][lr][al][th]->GetBinContent(j);
263
264 double dsb = m_hSigmaBiased[il][lr][al][th]->GetBinError(j);
265 double dsu = m_hSigmaUnbiased[il][lr][al][th]->GetBinError(j);
266 double XL = m_hSigmaBiased[il][lr][al][th]->GetXaxis()->GetBinCenter(j);
267 double dXL = (m_hSigmaBiased[il][lr][al][th]->GetXaxis()->GetBinWidth(j)) / 2;
268 double s_int = std::sqrt(sb * su);
269 double ds_int = 0.5 * s_int * (dsb / sb + dsu / su);
270 if (ds_int > 0.02) continue;
271 xl.push_back(XL);
272 dxl.push_back(dXL);
273 dxl0.push_back(0.);
274 sigma.push_back(s_int);
275 dsigma.push_back(ds_int);
276 s2.push_back(s_int * s_int);
277 ds2.push_back(2 * s_int * ds_int);
278 ofss << il << " " << lr << " " << al << " " << th << " " << j << " " << XL << " " << dXL << " " << s_int << " " <<
279 ds_int << endl;
280 }
281
282 if (xl.size() < 7 || xl.size() > Max_np) {
283 m_fitStatus[il][lr][al][th] = -1;
284 B2WARNING("number of element might out of range"); continue;
285 }
286
287 //Intrinsic resolution
288 B2DEBUG(21, "Create Histo for layer-lr: " << il << " " << lr);
289 m_graph[il][lr][al][th] = new TGraphErrors(xl.size(), &xl.at(0), &sigma.at(0), &dxl.at(0), &dsigma.at(0));
290 m_graph[il][lr][al][th]->SetMarkerSize(0.5);
291 m_graph[il][lr][al][th]->SetMarkerStyle(8);
292 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]));
293 m_graph[il][lr][al][th]->SetName(Form("lay%d_lr%d_al%d_th%d", il, lr, al, th));
294
295 //square of sigma for fitting
296 m_gFit[il][lr][al][th] = new TGraphErrors(xl.size(), &xl.at(0), &s2.at(0), &dxl0.at(0), &ds2.at(0));
297 m_gFit[il][lr][al][th]->SetMarkerSize(0.5);
298 m_gFit[il][lr][al][th]->SetMarkerStyle(8);
299 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]));
300 m_gFit[il][lr][al][th]->SetName(Form("sigma2_lay%d_lr%d_al%d_th%d", il, lr, al, th));
301
302 gDirectory->Delete("hu_%d_%d_%d_%d_0");
303 }
304 }
305 }
306 }
307 ofss.close();
308
309}
310
311CalibrationAlgorithm::EResult SpaceResolutionCalibrationAlgorithm::calibrate()
312{
313
314 B2INFO("Start calibration");
315 gPrintViaErrorHandler = true; // Suppress huge log output from TMinuit
316 gROOT->SetBatch(1);
317 gErrorIgnoreLevel = 3001;
318
319 const auto exprun = getRunList()[0];
320 B2INFO("ExpRun used for DB Geometry : " << exprun.first << " " << exprun.second);
321 updateDBObjPtrs(1, exprun.second, exprun.first);
322 // B2INFO("Creating CDCGeometryPar object");
323 // CDC::CDCGeometryPar::Instance(&(*m_cdcGeo));
324
325 prepare();
326 createHisto();
327
328 TF1* func = new TF1("func", "[0]/(x*x + [1])+[2]* x+[3]+[4]*exp([5]*(x-[6])*(x-[6]))", 0, 1.);
329 TH1F* hprob = new TH1F("h1", "", 20, 0, 1);
330 double upFit;
331 double intp6;
332
333 for (int i = 0; i < 56; ++i) {
334 for (int lr = 0; lr < 2; ++lr) {
335 for (int al = 0; al < m_nAlphaBins; ++al) {
336 for (int th = 0; th < m_nThetaBins; ++th) {
337 if (!m_gFit[i][lr][al][th]) continue;
338 if (m_fitStatus[i][lr][al][th] != -1) { /*if graph exist, do fitting*/
339 upFit = getUpperBoundaryForFit(m_gFit[i][lr][al][th]);
340 intp6 = upFit + 0.2;
341 B2DEBUG(199, "xmax for fitting: " << upFit);
342
343 func->SetParameters(5E-6, 0.007, 1E-4, 1E-5, 0.00008, -30, intp6);
344 func->SetParLimits(0, 1E-7, 1E-4);
345 func->SetParLimits(1, 0.0045, 0.02);
346 func->SetParLimits(2, 1E-6, 0.0005);
347 func->SetParLimits(3, 1E-8, 0.0005);
348 func->SetParLimits(4, 0., 0.001);
349 func->SetParLimits(5, -40, 0.);
350 func->SetParLimits(6, intp6 - 0.5, intp6 + 0.2);
351
352 B2DEBUG(21, "Fitting for layer: " << i << "lr: " << lr << " ial" << al << " ith:" << th);
353 B2DEBUG(21, "Fit status before fit:" << m_fitStatus[i][lr][al][th]);
354
355 for (int j = 0; j < 10; j++) {
356
357 B2DEBUG(21, "loop: " << j);
358 B2DEBUG(21, "Int p6: " << intp6);
359 B2DEBUG(21, "Number of Point: " << m_gFit[i][lr][al][th]->GetN());
360 Int_t stat = m_gFit[i][lr][al][th]->Fit("func", "MQE", "", 0.05, upFit);
361 B2DEBUG(21, "stat of fit" << stat);
362 std::string Fit_status = gMinuit->fCstatu.Data();
363 B2DEBUG(21, "FIT STATUS: " << Fit_status);
364 if (Fit_status == "OK" || Fit_status == "SUCCESSFUL" || Fit_status == "CALL LIMIT"
365 || Fit_status == "PROBLEMS") {//need to found better way
366 if (fabs(func->Eval(0.3)) > 0.00035 || func->Eval(0.3) < 0) {
367 func->SetParameters(5E-6, 0.007, 1E-4, 1E-7, 0.0007, -30, intp6 + 0.05 * j);
368 func->SetParLimits(6, intp6 + 0.05 * j - 0.5, intp6 + 0.05 * j + 0.2);
369 // func->SetParameters(defaultparsmall);
370 m_fitStatus[i][lr][al][th] = 0;
371 } else {
372 B2DEBUG(21, "Prob of fit: " << func->GetProb());
373 m_fitStatus[i][lr][al][th] = 1;
374 break;
375 }
376 } else {
377 m_fitStatus[i][lr][al][th] = 0;
378 func->SetParameters(5E-6, 0.007, 1E-4, 1E-7, 0.0007, -30, intp6 + 0.05 * j);
379 func->SetParLimits(6, intp6 + 0.05 * j - 0.5, intp6 + 0.05 * j + 0.2);
380 upFit += 0.025;
381 if (j == 9) {
382 // TCanvas* c1 = new TCanvas("c1", "", 600, 600);
383 // m_gFit[i][lr][al][th]->Draw();
384 // c1->SaveAs(Form("Sigma_Fit_Error_%s_%d_%d_%d_%d.png", Fit_status.c_str(), i, lr, al, th));
385 // B2WARNING("Fit error: " << i << " " << lr << " " << al << " " << th);
386 }
387 }
388 }
389 if (m_fitStatus[i][lr][al][th] == 1) {
390 B2DEBUG(21, "ProbFit: Lay_lr_al_th: " << i << " " << lr << " " << al << " " << th << func->GetProb());
391 hprob->Fill(func->GetProb());
392 func->GetParameters(m_sigma[i][lr][al][th]);
393 }
394 }
395 }
396 }
397 }
398 }
399
400 write();
401 storeHisto();
402
403 const int nTotal = 56 * 2 * m_nAlphaBins * m_nThetaBins;
404 int nFitCompleted = 0;
405 for (int l = 0; l < 56; ++l) {
406 for (int lr = 0; lr < 2; ++lr) {
407 for (int al = 0; al < m_nAlphaBins; ++al) {
408 for (int th = 0; th < m_nThetaBins; ++th) {
409 if (m_fitStatus[l][lr][al][th] == 1) {
410 nFitCompleted += 1;
411 }
412 }
413 }
414 }
415 }
416
417 if (static_cast<double>(nFitCompleted) / nTotal < m_threshold) {
418 B2WARNING("Less than " << m_threshold * 100 << " % of Sigmas were fitted.");
419 return c_NotEnoughData;
420 }
421
422 return c_OK;
423}
424
425void SpaceResolutionCalibrationAlgorithm::storeHisto()
426{
427 B2INFO("saving histograms");
428
429 TFile* ff = new TFile(m_histName.c_str(), "RECREATE");
430 TDirectory* top = gDirectory;
431 TDirectory* Direct[56];
432
433 auto hNDF = getObjectPtr<TH1F>("hNDF");
434 auto hPval = getObjectPtr<TH1F>("hPval");
435 auto hEvtT0 = getObjectPtr<TH1F>("hEventT0");
436 //store NDF, P-val. EventT0 histogram for monitoring during calibration
437 if (hNDF && hPval && hEvtT0) {
438 hEvtT0->Write();
439 hPval->Write();
440 hNDF->Write();
441 }
442
443
444 for (int il = 0; il < 56; ++il) {
445 top->cd();
446 Direct[il] = gDirectory->mkdir(Form("lay_%d", il));
447 Direct[il]->cd();
448
449 for (int lr = 0; lr < 2; ++lr) {
450 for (int al = 0; al < m_nAlphaBins; ++al) {
451 for (int th = 0; th < m_nThetaBins; ++th) {
452 if (!m_graph[il][lr][al][th]) continue;
453 if (!m_gFit[il][lr][al][th]) continue;
454 if (m_fitStatus[il][lr][al][th] == 1) {
455 m_hBiased[il][lr][al][th]->Write();
456 m_hUnbiased[il][lr][al][th]->Write();
457 m_hMeanBiased[il][lr][al][th]->Write();
458 m_hSigmaBiased[il][lr][al][th]->Write();
459 m_hMeanUnbiased[il][lr][al][th]->Write();
460 m_hSigmaUnbiased[il][lr][al][th]->Write();
461 m_graph[il][lr][al][th]->Write();
462 m_gFit[il][lr][al][th]->Write();
463 }
464 }
465 }
466 }
467 }
468 ff->Close();
469 B2INFO("Finish store histogram");
470
471}
472void SpaceResolutionCalibrationAlgorithm::write()
473{
474 B2INFO("Writing calibrated sigma's");
475 int nfitted = 0;
476 int nfailure = 0;
477
478 CDCSpaceResols* dbSigma = new CDCSpaceResols();
479
480 const float deg2rad = M_PI / 180.0;
481
482 for (unsigned short i = 0; i < m_nAlphaBins; ++i) {
483 std::array<float, 3> alpha3 = {m_lowerAlpha[i]* deg2rad,
484 m_upperAlpha[i]* deg2rad,
485 m_iAlpha[i]* deg2rad
486 };
487 dbSigma->setAlphaBin(alpha3);
488 }
489
490
491 for (unsigned short i = 0; i < m_nThetaBins; ++i) {
492 std::array<float, 3> theta3 = {m_lowerTheta[i]* deg2rad,
493 m_upperTheta[i]* deg2rad,
494 m_iTheta[i]* deg2rad
495 };
496 dbSigma->setThetaBin(theta3);
497 }
498
499 dbSigma->setSigmaParamMode(m_sigmaParamMode);
500 for (int ialpha = 0; ialpha < m_nAlphaBins; ++ialpha) {
501 for (int itheta = 0; itheta < m_nThetaBins; ++itheta) {
502 for (int iCL = 0; iCL < 56; ++iCL) {
503 for (int iLR = 1; iLR >= 0; --iLR) {
504 std::vector<float> sgbuff;
505 if (m_fitStatus[iCL][iLR][ialpha][itheta] == 1) {
506 nfitted += 1; // inclement number of successfully fitted sigma's
507 for (int i = 0; i < 7; ++i) {
508 sgbuff.push_back(m_sigma[iCL][iLR][ialpha][itheta][i]);
509 }
510 } else {
511 //B2WARNING("Fitting error and old sigma will be used. (Layer " << iCL << ") (lr = " << iLR <<
512 // ") (al = " << ialpha << ") (th = " << itheta << ")");
513 nfailure += 1; // inclement number of fit failed sigma's
514 for (int i = 0; i < 7; ++i) {
515 sgbuff.push_back(m_sigmaPost[iCL][iLR][ialpha][itheta][i]);
516 }
517 }
518 dbSigma->setSigmaParams(iCL, iLR, ialpha, itheta, sgbuff);
519 }
520 }
521 }
522 }
523
524 if (m_textOutput == true) {
525 dbSigma->outputToFile(m_outputFileName);
526 }
527
528 saveCalibration(dbSigma, "CDCSpaceResols");
529
530 B2RESULT("Number of histogram: " << 56 * 2 * m_nAlphaBins * m_nThetaBins);
531 B2RESULT("Histos successfully fitted: " << nfitted);
532 B2RESULT("Histos fit failure: " << nfailure);
533
534
535}
536
537void SpaceResolutionCalibrationAlgorithm::prepare()
538{
539 B2INFO("Prepare calibration of space resolution");
540
541 const double rad2deg = 180 / M_PI;
542
543 DBObjPtr<CDCSpaceResols> dbSigma;
544
545 m_nAlphaBins = dbSigma->getNoOfAlphaBins();
546 m_nThetaBins = dbSigma->getNoOfThetaBins();
547
548 B2INFO("Number of alpha bins: " << m_nAlphaBins);
549 for (int i = 0; i < m_nAlphaBins; ++i) {
550 array3 alpha = dbSigma->getAlphaBin(i);
551 m_lowerAlpha[i] = alpha[0] * rad2deg;
552 m_upperAlpha[i] = alpha[1] * rad2deg;
553 m_iAlpha[i] = alpha[2] * rad2deg;
554 }
555
556 B2INFO("Number of theta bins: " << m_nThetaBins);
557 for (int i = 0; i < m_nThetaBins; ++i) {
558 array3 theta = dbSigma->getThetaBin(i);
559 m_lowerTheta[i] = theta[0] * rad2deg;
560 m_upperTheta[i] = theta[1] * rad2deg;
561 m_iTheta[i] = theta[2] * rad2deg;
562 }
563 m_sigmaParamModePost = dbSigma->getSigmaParamMode();
564
565 for (unsigned short iCL = 0; iCL < 56; ++iCL) {
566 for (unsigned short iLR = 0; iLR < 2; ++iLR) {
567 for (unsigned short iA = 0; iA < m_nAlphaBins; ++iA) {
568 for (unsigned short iT = 0; iT < m_nThetaBins; ++iT) {
569 const std::vector<float> params = dbSigma->getSigmaParams(iCL, iLR, iA, iT);
570 unsigned short np = params.size();
571 // std::cout <<"np4sigma= " << np << std::endl;
572 for (unsigned short i = 0; i < np; ++i) {
573 m_sigmaPost[iCL][iLR][iA][iT][i] = params[i];
574 }
575 }
576 }
577 }
578 }
579}
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:537
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
void updateDBObjPtrs(const unsigned int event, const int run, const int experiment)
Updates any DBObjPtrs by calling update(event) for DBStore.
Definition CalibrationAlgorithm.cc:404
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.