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