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