Belle II Software  release-05-02-19
SpaceResolutionCalibrationAlgorithm.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Makoto Uchida *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <iostream>
12 #include <cdc/calibration/SpaceResolutionCalibrationAlgorithm.h>
13 #include <cdc/geometry/CDCGeometryPar.h>
14 #include <cdc/dbobjects/CDCSpaceResols.h>
15 
16 #include <framework/database/DBObjPtr.h>
17 #include <framework/logging/Logger.h>
18 
19 #include "TF1.h"
20 #include "TH1F.h"
21 #include "TH2F.h"
22 #include "TROOT.h"
23 #include "TError.h"
24 #include "TMinuit.h"
25 
26 using namespace std;
27 using namespace Belle2;
28 using namespace CDC;
29 
30 typedef std::array<float, 3> array3;
31 SpaceResolutionCalibrationAlgorithm::SpaceResolutionCalibrationAlgorithm() :
32  CalibrationAlgorithm("CDCCalibrationCollector")
33 {
34  setDescription(
35  " -------------------------- Position Resolution Calibration Algorithm -------------------------\n"
36  );
37 }
38 void SpaceResolutionCalibrationAlgorithm::createHisto()
39 {
40  B2INFO("Creating histograms");
41  const int np = floor(1 / m_binWidth);
42 
43  vector<double> yu;
44  vector <double> yb;
45  for (int i = 0; i < 50; ++i) {
46  yb.push_back(-0.07 + i * (0.14 / 50));
47  }
48  for (int i = 0; i < 50; ++i) {
49  yu.push_back(-0.08 + i * (0.16 / 50));
50  }
51 
52  vector<double> xbin;
53  xbin.push_back(0.);
54  xbin.push_back(0.02);
55  for (int i = 1; i < np; ++i) {
56  xbin.push_back(i * m_binWidth);
57  }
58 
59  for (int il = 0; il < 56; ++il) {
60  for (int lr = 0; lr < 2; ++lr) {
61  for (int al = 0; al < m_nAlphaBins; ++al) {
62  for (int th = 0; th < m_nThetaBins; ++th) {
63  m_hBiased[il][lr][al][th] = new TH2F(Form("hb_%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), yb.size() - 1, &yb.at(0));
66  m_hUnbiased[il][lr][al][th] = new TH2F(Form("hu_%d_%d_%d_%d", il, lr, al, th),
67  Form("lay_%d_lr%d_al_%3.0f_th_%3.0f;Drift Length [cm];#DeltaX", il, lr, m_iAlpha[al], m_iTheta[th]),
68  xbin.size() - 1, &xbin.at(0), yu.size() - 1, &yu.at(0));
69  }
70  }
71  }
72  }
73 
74 
75  auto tree = getObjectPtr<TTree>("tree");
76 
77  UChar_t lay;
78  Float_t w;
79  Float_t x_u;
80  Float_t x_b;
81  Float_t x_mea;
82  Float_t Pval;
83  Float_t alpha;
84  Float_t theta;
85  Float_t ndf;
86  Float_t absRes_u;
87  Float_t absRes_b;
88  tree->SetBranchAddress("lay", &lay);
89  tree->SetBranchAddress("ndf", &ndf);
90  tree->SetBranchAddress("Pval", &Pval);
91  tree->SetBranchAddress("x_u", &x_u);
92  tree->SetBranchAddress("x_b", &x_b);
93  tree->SetBranchAddress("x_mea", &x_mea);
94  tree->SetBranchAddress("weight", &w);
95  tree->SetBranchAddress("alpha", &alpha);
96  tree->SetBranchAddress("theta", &theta);
97 
98  /* Disable unused branch */
99  std::vector<TString> list_vars = {"lay", "ndf", "Pval", "x_u", "x_b", "x_mea", "weight", "alpha", "theta"};
100  tree->SetBranchStatus("*", 0);
101 
102  for (TString brname : list_vars) {
103  tree->SetBranchStatus(brname, 1);
104  }
105 
106 
107  const Long64_t nEntries = tree->GetEntries();
108  B2INFO("Number of entries: " << nEntries);
109  int ith = -99;
110  int ial = -99;
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 
146  B2INFO("Start to obtain the biased and unbiased sigmas");
147 
148  TF1* gb = new TF1("gb", "gausn", -0.05, 0.05);
149  TF1* gu = new TF1("gu", "gausn", -0.06, 0.06);
150  TF1* g0b = new TF1("g0b", "gausn", -0.015, 0.07);
151  TF1* g0u = new TF1("g0u", "gausn", -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 
311 CalibrationAlgorithm::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 
425 void 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 }
472 void 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 succesfully fitted: " << nfitted);
532  B2RESULT("Histos fit failure: " << nfailure);
533 
534 
535 }
536 
537 void 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 }
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_sigmaPost
double m_sigmaPost[56][2][18][7][8]
sigma prameters before calibration
Definition: SpaceResolutionCalibrationAlgorithm.h:149
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::getUpperBoundaryForFit
double getUpperBoundaryForFit(TGraphErrors *graph)
search max point at boundary region
Definition: SpaceResolutionCalibrationAlgorithm.h:91
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_iAlpha
float m_iAlpha[18]
represented alphas of alpha bins.
Definition: SpaceResolutionCalibrationAlgorithm.h:143
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_threshold
double m_threshold
minimal requirement for the fraction of fitted results
Definition: SpaceResolutionCalibrationAlgorithm.h:127
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::createHisto
void createHisto()
create histogram
Definition: SpaceResolutionCalibrationAlgorithm.cc:38
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_sigmaParamMode
unsigned short m_sigmaParamMode
sigma mode for this calibration.
Definition: SpaceResolutionCalibrationAlgorithm.h:147
Belle2::CalibrationAlgorithm::saveCalibration
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
Definition: CalibrationAlgorithm.cc:290
Belle2::CDCSpaceResols::setAlphaBin
void setAlphaBin(const array3 &alpha)
Set alpha-angle bin (rad)
Definition: CDCSpaceResols.h:60
Belle2::CDCSpaceResols::setSigmaParamMode
void setSigmaParamMode(unsigned short mode)
Set sigma parameterization mode.
Definition: CDCSpaceResols.h:95
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_hUnbiased
TH2F * m_hUnbiased[56][2][Max_nalpha][Max_ntheta]
2D histogram of unbiased residual
Definition: SpaceResolutionCalibrationAlgorithm.h:132
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_sigma
double m_sigma[56][2][18][7][8]
new sigma prameters.
Definition: SpaceResolutionCalibrationAlgorithm.h:128
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_minNdf
double m_minNdf
Minimum NDF
Definition: SpaceResolutionCalibrationAlgorithm.h:121
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfuly =0 in Python.
Definition: CalibrationAlgorithm.h:51
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition: CalibrationAlgorithm.h:331
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:398
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_hSigmaBiased
TH1F * m_hSigmaBiased[56][2][Max_nalpha][Max_ntheta]
sigma histogram of biased residual
Definition: SpaceResolutionCalibrationAlgorithm.h:134
Belle2::DBObjPtr
Class for accessing objects in the database.
Definition: DBObjPtr.h:31
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_nThetaBins
int m_nThetaBins
number of theta bins
Definition: SpaceResolutionCalibrationAlgorithm.h:140
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::calibrate
EResult calibrate() override
Run algo on data.
Definition: SpaceResolutionCalibrationAlgorithm.cc:311
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_lowerAlpha
float m_lowerAlpha[18]
Lower boundays of alpha bins.
Definition: SpaceResolutionCalibrationAlgorithm.h:141
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_gFit
TGraphErrors * m_gFit[56][2][18][7]
sigma*sigma graph for fit
Definition: SpaceResolutionCalibrationAlgorithm.h:129
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_minPval
double m_minPval
Minimum Prob(chi2) of track.
Definition: SpaceResolutionCalibrationAlgorithm.h:122
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_histName
std::string m_histName
root file name
Definition: SpaceResolutionCalibrationAlgorithm.h:154
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_nAlphaBins
int m_nAlphaBins
number of alpha bins
Definition: SpaceResolutionCalibrationAlgorithm.h:139
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_hMeanUnbiased
TH1F * m_hMeanUnbiased[56][2][Max_nalpha][Max_ntheta]
mean histogram of unbiased residual
Definition: SpaceResolutionCalibrationAlgorithm.h:135
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_upperAlpha
float m_upperAlpha[18]
Upper boundays of alpha bins.
Definition: SpaceResolutionCalibrationAlgorithm.h:142
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_hSigmaUnbiased
TH1F * m_hSigmaUnbiased[56][2][Max_nalpha][Max_ntheta]
sigma histogram of ubiased residual
Definition: SpaceResolutionCalibrationAlgorithm.h:136
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::storeHisto
void storeHisto()
store histogram
Definition: SpaceResolutionCalibrationAlgorithm.cc:425
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_textOutput
bool m_textOutput
output text file if true
Definition: SpaceResolutionCalibrationAlgorithm.h:152
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition: CalibrationAlgorithm.h:50
Belle2::CalibrationAlgorithm::getRunList
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
Definition: CalibrationAlgorithm.h:276
Belle2::CDCSpaceResols
Database object for space resolutions.
Definition: CDCSpaceResols.h:36
Belle2::CDCSpaceResols::setSigmaParams
void setSigmaParams(const SigmaID sigmaID, const std::vector< float > &params)
Set sigma parameters for the specified id.
Definition: CDCSpaceResols.h:103
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_hBiased
TH2F * m_hBiased[56][2][Max_nalpha][Max_ntheta]
2D histogram of biased residual
Definition: SpaceResolutionCalibrationAlgorithm.h:131
Belle2::CalibrationAlgorithm::c_NotEnoughData
@ c_NotEnoughData
Needs more data =2 in Python.
Definition: CalibrationAlgorithm.h:53
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::write
void write()
save calibration, in text file or db
Definition: SpaceResolutionCalibrationAlgorithm.cc:472
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_iTheta
float m_iTheta[7]
represented alphas of theta bins.
Definition: SpaceResolutionCalibrationAlgorithm.h:146
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::prepare
void prepare()
Prepare the calibration of space resolution.
Definition: SpaceResolutionCalibrationAlgorithm.cc:537
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::Max_np
static const unsigned short Max_np
Maximum number of point =1/binwidth.
Definition: SpaceResolutionCalibrationAlgorithm.h:119
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_lowerTheta
float m_lowerTheta[7]
Lower boundays of theta bins.
Definition: SpaceResolutionCalibrationAlgorithm.h:144
Belle2::CalibrationAlgorithm
Base class for calibration algorithms.
Definition: CalibrationAlgorithm.h:47
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_graph
TGraphErrors * m_graph[56][2][18][7]
sigma graph.
Definition: SpaceResolutionCalibrationAlgorithm.h:130
Belle2::CDCSpaceResols::setThetaBin
void setThetaBin(const array3 &theta)
Set theta-angle bin (rad)
Definition: CDCSpaceResols.h:74
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_hMeanBiased
TH1F * m_hMeanBiased[56][2][Max_nalpha][Max_ntheta]
mean histogram biased residual
Definition: SpaceResolutionCalibrationAlgorithm.h:133
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_fitStatus
int m_fitStatus[56][2][Max_nalpha][Max_ntheta]
Fit flag; 1:OK ; 0:error.
Definition: SpaceResolutionCalibrationAlgorithm.h:137
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_outputFileName
std::string m_outputFileName
Output sigma filename.
Definition: SpaceResolutionCalibrationAlgorithm.h:153
Belle2::CDCSpaceResols::outputToFile
void outputToFile(std::string fileName) const
Output the contents in text file format.
Definition: CDCSpaceResols.h:345
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_sigmaParamModePost
unsigned short m_sigmaParamModePost
sigma mode before this calibration.
Definition: SpaceResolutionCalibrationAlgorithm.h:150
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_upperTheta
float m_upperTheta[7]
Upper boundays of theta bins.
Definition: SpaceResolutionCalibrationAlgorithm.h:145
Belle2::CDC::SpaceResolutionCalibrationAlgorithm::m_binWidth
double m_binWidth
width of each bin, unit cm
Definition: SpaceResolutionCalibrationAlgorithm.h:123