Belle II Software  release-05-01-25
SpaceResolutionCalibration.cc
1 #include <iostream>
2 #include <iomanip>
3 #include <cdc/calibration/SpaceResolutionCalibration.h>
4 #include <cdc/geometry/CDCGeometryPar.h>
5 #include <cdc/calibration/CDCDatabaseImporter.h>
6 
7 #include <framework/database/DBObjPtr.h>
8 #include <framework/logging/Logger.h>
9 #include <framework/utilities/FileSystem.h>
10 
11 #include "TFile.h"
12 #include "TH1F.h"
13 #include "TH2F.h"
14 #include "TCanvas.h"
15 #include "TSystem.h"
16 #include "TChain.h"
17 #include "TROOT.h"
18 #include "TError.h"
19 #include "TMinuit.h"
20 
21 using namespace std;
22 using namespace Belle2;
23 using namespace CDC;
24 
25 SpaceResolutionCalibration::SpaceResolutionCalibration():
26  m_firstExperiment(0), m_firstRun(0),
27  m_lastExperiment(-1), m_lastRun(-1)
28 {
29  /*Space resolution calibration*/
30 }
31 void SpaceResolutionCalibration::createHisto()
32 {
33 
34  B2INFO("createHisto");
35  readSigma();
36  readProfile();
37  const int m_np = floor(1 / m_binWidth);
38 
39  TChain* tree = new TChain("tree");
40  tree->Add(m_inputRootFileNames.c_str());
41  B2INFO(" Open file name: " << m_inputRootFileNames.c_str());
42  if (!tree->GetBranch("ndf")) {
43  B2FATAL("input data do not exits, please check!");
44  gSystem->Exec("echo rootfile do not exits or something wrong >> error");
45  return;
46  }
47 
48  int lay;
49  double w;
50  double x_u;
51  double x_b;
52  double x_mea;
53  double Pval;
54  double alpha;
55  double theta;
56  double ndf;
57  double absRes_u;
58  double absRes_b;
59  tree->SetBranchAddress("lay", &lay);
60  tree->SetBranchAddress("ndf", &ndf);
61  tree->SetBranchAddress("Pval", &Pval);
62  tree->SetBranchAddress("x_u", &x_u);
63  tree->SetBranchAddress("x_b", &x_b);
64  tree->SetBranchAddress("x_mea", &x_mea);
65  tree->SetBranchAddress("weight", &w);
66  tree->SetBranchAddress("alpha", &alpha);
67  tree->SetBranchAddress("theta", &theta);
68 
69  /* Disable unused branch */
70  std::vector<TString> list_vars = {"lay", "ndf", "Pval", "x_u", "x_b", "x_mea", "weight", "alpha", "theta"};
71  tree->SetBranchStatus("*", 0);
72 
73  for (TString brname : list_vars) {
74  tree->SetBranchStatus(brname, 1);
75  }
76 
77 
78  vector<double> yu;
79  vector <double> yb;
80  for (int i = 0; i < 50; ++i) {
81  yb.push_back(-0.07 + i * (0.14 / 50));
82  }
83  for (int i = 0; i < 50; ++i) {
84  yu.push_back(-0.08 + i * (0.16 / 50));
85  }
86 
87  vector<double> xbin;
88  xbin.push_back(0.);
89  xbin.push_back(0.02);
90  for (int i = 1; i < m_np; ++i) {
91  xbin.push_back(i * m_binWidth);
92  }
93 
94  for (int il = 0; il < 56; ++il) {
95  for (int lr = 0; lr < 2; ++lr) {
96  for (int al = 0; al < m_nalpha; ++al) {
97  for (int th = 0; th < m_ntheta; ++th) {
98  hist_b[il][lr][al][th] = new TH2F(Form("hb_%d_%d_%d_%d", il, lr, al, th),
99  Form("lay_%d_lr%d_al_%3.0f_th_%3.0f;Drift Length [cm];#DeltaX", il, lr, ialpha[al], itheta[th]),
100  xbin.size() - 1, &xbin.at(0), yb.size() - 1, &yb.at(0));
101  hist_u[il][lr][al][th] = new TH2F(Form("hu_%d_%d_%d_%d", il, lr, al, th),
102  Form("lay_%d_lr%d_al_%3.0f_th_%3.0f;Drift Length [cm];#DeltaX", il, lr, ialpha[al], itheta[th]),
103  xbin.size() - 1, &xbin.at(0), yu.size() - 1, &yu.at(0));
104  }
105  }
106  }
107  }
108 
109 
110  const int nEntries = tree->GetEntries();
111  B2INFO("Number of entries: " << nEntries);
112  int ith = -99;
113  int ial = -99;
114  int ilr = -99;
115  for (int i = 0; i < nEntries; ++i) {
116  tree->GetEntry(i);
117  //cut
118  if (std::fabs(x_b) < 0.02 || std::fabs(x_u) < 0.02) continue;
119  if (Pval < m_Pvalmin) continue;
120  if (ndf < m_ndfmin) continue;
121  for (int k = 0; k < m_nalpha; ++k) {
122  if (alpha < u_alpha[k]) {
123  ial = k;
124  break;
125  }
126  }
127 
128  for (int j = 0; j < m_ntheta; ++j) {
129  if (theta < u_theta[j]) {
130  ith = j;
131  break;
132  }
133  }
134 
135  ilr = x_u > 0 ? 1 : 0;
136 
137  if (ial == -99 || ith == -99 || ilr == -99) {
138  TString command = Form("Error in alpha=%3.2f and theta = %3.2f>> error", alpha, theta);
139  // gSystem->Exec(command);
140  B2FATAL("ERROR" << command);
141  }
142  absRes_u = fabs(x_mea) - fabs(x_u);
143  absRes_b = fabs(x_mea) - fabs(x_b);
144 
145  hist_u[lay][ilr][ial][ith]->Fill(fabs(x_u), absRes_u, w);
146  hist_b[lay][ilr][ial][ith]->Fill(fabs(x_b), absRes_b, w);
147 
148  }
149 
150  B2INFO("Finish reading data");
151 
152  TF1* gb = new TF1("gb", "gaus", -0.05, 0.05);
153  TF1* gu = new TF1("gu", "gaus", -0.06, 0.06);
154  TF1* g0b = new TF1("g0b", "gaus", -0.015, 0.07);
155  TF1* g0u = new TF1("g0u", "gaus", -0.015, 0.08);
156  g0b->SetParLimits(1, -0.01, 0.004);
157  g0u->SetParLimits(1, -0.01, 0.004);
158 
159  std::vector<double> sigma;
160  std::vector<double> dsigma;
161  std::vector<double> s2;
162  std::vector<double> ds2;
163  std::vector<double> xl;
164  std::vector<double> dxl;
165  std::vector<double> dxl0;
166 
167  const int ib1 = int(0.1 / m_binWidth) + 1;
168  int firstbin = 1;
169  int minEntry = 10;
170  for (int il = 0; il < 56; ++il) {
171  for (int lr = 0; lr < 2; ++lr) {
172  for (int al = 0; al < m_nalpha; ++al) {
173  for (int th = 0; th < m_ntheta; ++th) {
174  //fit half gaus for first range near sense wire
175  B2DEBUG(199, "layer-lr-al-th " << il << " - " << lr << " - " << al << " - " << th);
176  if (hist_b[il][lr][al][th]->GetEntries() < 20000) {
177  m_fitflag[il][lr][al][th] = -1;
178  continue;
179  }
180  B2DEBUG(199, "Nentries: " << hist_b[il][lr][al][th]->GetEntries());
181  hist_b[il][lr][al][th]->SetDirectory(0);
182  hist_u[il][lr][al][th]->SetDirectory(0);
183 
184  hist_b[il][lr][al][th]->FitSlicesY(g0b, firstbin, ib1, minEntry);
185 
186  TH1F* hm1 = (TH1F*)gDirectory->Get(Form("hb_%d_%d_%d_%d_1", il, lr, al, th));
187  TH1F* hs1 = (TH1F*)gDirectory->Get(Form("hb_%d_%d_%d_%d_2", il, lr, al, th));
188  if (!hm1 || !hs1) {
189  m_fitflag[il][lr][al][th] = -1;
190  continue;
191  }
192 
193  hb_m[il][lr][al][th] = (TH1F*)hm1->Clone(Form("hb_%d_%d_%d_%d_m", il, lr, al, th));
194  hb_s[il][lr][al][th] = (TH1F*)hs1->Clone(Form("hb_%d_%d_%d_%d_s", il, lr, al, th));//sigma
195 
196  hb_m[il][lr][al][th]->SetDirectory(0);
197  hb_s[il][lr][al][th]->SetDirectory(0);
198 
199 
200  //slice other bin with full gaus func
201  hist_b[il][lr][al][th]->FitSlicesY(gb, ib1 + 1, m_np, minEntry);
202  hb_m[il][lr][al][th]->Add((TH1F*)gDirectory->Get(Form("hb_%d_%d_%d_%d_1", il, lr, al, th))); //mean
203  hb_s[il][lr][al][th]->Add((TH1F*)gDirectory->Get(Form("hb_%d_%d_%d_%d_2", il, lr, al, th))); //sigma
204  B2DEBUG(199, "entries (2nd): " << hb_s[il][lr][al][th]->GetEntries());
205 
206  //fit half gaus for first range near sense wire
207  hist_u[il][lr][al][th]->FitSlicesY(g0u, firstbin, ib1, minEntry);
208  hu_m[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", il, lr, al,
209  th));//mean
210  hu_s[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", il, lr, al,
211  th));//sigma
212  hu_m[il][lr][al][th]->SetDirectory(0);
213  hu_s[il][lr][al][th]->SetDirectory(0);
214 
215  //slice other bin with full gaus func
216  hist_u[il][lr][al][th]->FitSlicesY(gu, ib1 + 1, m_np, minEntry);
217  hu_m[il][lr][al][th]->Add((TH1F*)gDirectory->Get(Form("hu_%d_%d_%d_%d_1", il, lr, al, th))); //mean
218  hu_s[il][lr][al][th]->Add((TH1F*)gDirectory->Get(Form("hu_%d_%d_%d_%d_2", il, lr, al, th))); //sigma
219  if (!hu_s[il][lr][al][th] || !hb_s[il][lr][al][th]) {
220  B2WARNING("slice histo do not found");
221  m_fitflag[il][lr][al][th] = -1;
222  continue;
223  }
224  for (Int_t j = 1; j < hu_s[il][lr][al][th]->GetNbinsX(); j++) {
225  if (hu_s[il][lr][al][th]->GetBinContent(j) == 0) continue;
226  if (hb_s[il][lr][al][th]->GetBinContent(j) == 0) continue;
227  double sb = hb_s[il][lr][al][th]->GetBinContent(j);
228  double su = hu_s[il][lr][al][th]->GetBinContent(j);
229 
230  double dsb = hb_s[il][lr][al][th]->GetBinError(j);
231  double dsu = hu_s[il][lr][al][th]->GetBinError(j);
232  double XL = hb_s[il][lr][al][th]->GetXaxis()->GetBinCenter(j);
233  double dXL = (hb_s[il][lr][al][th]->GetXaxis()->GetBinWidth(j)) / 2;
234  double s_int = std::sqrt(sb * su);
235  double ds_int = 0.5 * s_int * (dsb / sb + dsu / su);
236  if (ds_int > 0.02) continue;
237  xl.push_back(XL);
238  dxl.push_back(dXL);
239  dxl0.push_back(0.);
240  sigma.push_back(s_int);
241  dsigma.push_back(ds_int);
242  s2.push_back(s_int * s_int);
243  ds2.push_back(2 * s_int * ds_int);
244  }
245 
246  if (xl.size() < 8 || xl.size() > Max_np) {
247  m_fitflag[il][lr][al][th] = -1;
248  B2WARNING("number of element might out of range"); continue;
249  }
250 
251  //Intrinsic resolution
252  B2DEBUG(199, "Create Histo for layer-lr: " << il << " " << lr);
253  gr[il][lr][al][th] = new TGraphErrors(xl.size(), &xl.at(0), &sigma.at(0), &dxl.at(0), &dsigma.at(0));
254  gr[il][lr][al][th]->SetMarkerSize(0.5);
255  gr[il][lr][al][th]->SetMarkerStyle(8);
256  // gr[il][lr][al][th]->SetMarkerColor(kBlack);
257  //gr[il][lr][al][th]->SetLineColor(1 + lr + al * 2 + th * 3);
258  gr[il][lr][al][th]->SetTitle(Form("Layer_%d_lr%d | #alpha = %3.0f | #theta = %3.0f", il, lr, ialpha[al], itheta[th]));
259  gr[il][lr][al][th]->SetName(Form("lay%d_lr%d_al%d_th%d", il, lr, al, th));
260 
261  //s2 for fitting
262  gfit[il][lr][al][th] = new TGraphErrors(xl.size(), &xl.at(0), &s2.at(0), &dxl0.at(0), &ds2.at(0));
263  gfit[il][lr][al][th]->SetMarkerSize(0.5);
264  gfit[il][lr][al][th]->SetMarkerStyle(8);
265  // gfit[il][lr][al][th]->SetMarkerColor(1 + lr + al * 2 + th * 3);
266  //gfit[il][lr][al][th]->SetLineColor(1 + lr + al * 2 + th * 3);
267  gfit[il][lr][al][th]->SetTitle(Form("L%d-lr%d | #alpha = %3.0f | #theta = %3.0f ", il, lr, ialpha[al], itheta[th]));
268  gfit[il][lr][al][th]->SetName(Form("sigma2_lay%d_lr%d_al%d_th%d", il, lr, al, th));
269 
270  xl.clear(); dxl.clear(); dxl0.clear(); sigma.clear(); dsigma.clear(); s2.clear(); ds2.clear();
271  gDirectory->Delete("hu_%d_%d_%d_%d_0");
272  }
273  }
274  }
275  }
276 }
277 
278 bool SpaceResolutionCalibration::calibrate()
279 {
280  gROOT->SetBatch(1);
281  gErrorIgnoreLevel = 3001;
282  createHisto();
283  B2INFO("Start to calibrate");
284  gSystem->Exec("mkdir -p Sigma_Fit_Err"); //create a folder to store error histo
285 
286  TF1* func = new TF1("func", "[0]/(x*x + [1])+[2]* x+[3]+[4]*exp([5]*(x-[6])*(x-[6]))", 0, 1.);
287  TH1F* hprob = new TH1F("h1", "", 20, 0, 1);
288  double upFit;
289  double intp6;
290 
291  for (int i = 0; i < 56; ++i) {
292  for (int lr = 0; lr < 2; ++lr) {
293  for (int al = 0; al < m_nalpha; ++al) {
294  for (int th = 0; th < m_ntheta; ++th) {
295  if (m_fitflag[i][lr][al][th] != -1) { /*if graph exist, do fitting*/
296 
297 
298  upFit = getUpperBoundaryForFit(gfit[i][lr][al][th]);
299  intp6 = upFit + 0.2;
300  B2DEBUG(199, "xmax for fitting: " << upFit);
301  // func->SetLineColor(1 + lr + al * 2 + th * 3);
302  func->SetParameters(5E-6, 0.007, 1E-4, 1E-5, 0.00008, -30, intp6);
303  func->SetParLimits(0, 1E-7, 1E-4);
304  func->SetParLimits(1, 0.00001, 0.02);
305  func->SetParLimits(2, 1E-6, 0.0005);
306  func->SetParLimits(3, 1E-8, 0.0005);
307  func->SetParLimits(4, 0., 0.001);
308  func->SetParLimits(5, -40, 0.);
309  func->SetParLimits(6, intp6 - 0.5, intp6 + 0.3);
310  B2DEBUG(199, "FITTING for layer: " << i << "lr: " << lr << " ial" << al << " ith:" << th);
311  B2DEBUG(199, "Fit flag before fit:" << m_fitflag[i][lr][al][th]);
312  // if(!(gfit[i][lr][al][th]->isValid())) continue;
313  // m_fitflag[i][lr][al][th] = 0;
314  for (int j = 0; j < 10; j++) {
315 
316  B2DEBUG(199, "loop: " << j);
317  B2DEBUG(199, "Int p6: " << intp6);
318  B2DEBUG(199, "Number of Point: " << gfit[i][lr][al][th]->GetN());
319  Int_t stat = gfit[i][lr][al][th]->Fit("func", "MQE", "", 0.05, upFit);
320  B2DEBUG(199, "stat of fit" << stat);
321  std::string Fit_status = gMinuit->fCstatu.Data();
322  B2DEBUG(199, "FIT STATUS: " << Fit_status);
323  // stat=gfit[i]->Fit(Form("ffit[%d]",i),"M "+Q,"",0.0,cellsize(i)+0.05+j*0.005);
324  if (Fit_status == "OK" ||
325  Fit_status == "SUCCESSFUL" ||
326  Fit_status == "CALL LIMIT" ||
327  Fit_status == "PROBLEMS") {
328  if (fabs(func->Eval(0.3)) > 0.00035 || func->Eval(0.3) < 0) {
329  func->SetParameters(5E-6, 0.007, 1E-4, 1E-7, 0.0007, -30, intp6 + 0.05 * j);
330  // func->SetParameters(defaultparsmall);
331  m_fitflag[i][lr][al][th] = 0;
332  } else {
333  B2DEBUG(199, "Prob of fit: " << func->GetProb());
334  m_fitflag[i][lr][al][th] = 1;
335  break;
336  }
337  } else {
338  m_fitflag[i][lr][al][th] = 0;
339  func->SetParameters(5E-6, 0.007, 1E-4, 1E-7, 0.0007, -30, intp6 + 0.05 * j);
340  // upFit += 0.025;
341  if (j == 9) {
342  TCanvas* c1 = new TCanvas("c1", "", 600, 600);
343  gfit[i][lr][al][th]->Draw();
344  c1->SaveAs(Form("Sigma_Fit_Err/%d_%d_al%d_th%d_%s.png", i, lr, al, th, Fit_status.c_str()));
345  B2WARNING("Fit error: " << i << " " << lr << " " << al << " " << th);
346  }
347  }
348  }
349  if (m_fitflag[i][lr][al][th] == 1) {
350  B2DEBUG(199, "ProbFit: Lay_lr_al_th: " << i << " " << lr << " " << al << " " << th << func->GetProb());
351  hprob->Fill(func->GetProb());
352  func->GetParameters(sigma_new[i][lr][al][th]);
353  }
354  }
355  }
356  }
357  }
358  }
359  storeHisto();
360  write();
361 
362  return true;
363 }
364 void SpaceResolutionCalibration::storeHisto()
365 {
366  B2INFO("storeHisto");
367  TFile* ff = new TFile("sigma_histo.root", "RECREATE");
368  TDirectory* top = gDirectory;
369  TDirectory* Direct[56];
370 
371  for (int il = 0; il < 56; ++il) {
372  top->cd();
373  Direct[il] = gDirectory->mkdir(Form("lay_%d", il));
374  Direct[il]->cd();
375  for (int lr = 0; lr < 2; ++lr) {
376  for (int al = 0; al < m_nalpha; ++al) {
377  for (int th = 0; th < m_ntheta; ++th) {
378  if (!gr[il][lr][al][th]) continue;
379  if (!gfit[il][lr][al][th]) continue;
380  if (m_fitflag[il][lr][al][th] == 1) {
381  hist_b[il][lr][al][th]->Write();
382  hist_u[il][lr][al][th]->Write();
383  hb_m[il][lr][al][th]->Write();
384  hb_s[il][lr][al][th]->Write();
385  hu_m[il][lr][al][th]->Write();
386  hu_s[il][lr][al][th]->Write();
387  gr[il][lr][al][th]->Write();
388  gfit[il][lr][al][th]->Write();
389  }
390  }
391  }
392  }
393  }
394  ff->Close();
395  B2INFO("Finish store histogram");
396 }
397 void SpaceResolutionCalibration::write()
398 {
399 
400  B2INFO("Exporting parameters...");
401  int nfitted = 0;
402  int nfailure = 0;
403  /* Write the fit params*/
404 
405  ofstream ofs(m_outputSigmaFileName.c_str());
406  ofs << m_nalpha << endl;
407  for (int i = 0; i < m_nalpha; ++i) {
408  ofs << std::setprecision(4) << l_alpha[i] << " " << std::setprecision(4) << u_alpha[i] << " " << std::setprecision(
409  4) << ialpha[i] << endl;
410  }
411 
412  ofs << m_ntheta << endl;
413  for (int i = 0; i < m_ntheta; ++i) {
414  ofs << std::setprecision(4) << l_theta[i] << " " << std::setprecision(4) << u_theta[i] << " " << std::setprecision(
415  4) << itheta[i] << endl;
416  }
417 
418  ofs << 0 << " " << 7 << endl; //mode and number of params;
419  for (int al = 0; al < m_nalpha; ++al) {
420  for (int th = 0; th < m_ntheta; ++th) {
421  for (int i = 0; i < 56; ++i) {
422  for (int lr = 1; lr >= 0; --lr) {
423  // ffit[i][lr][al][th]->GetParameters(par);
424  ofs << i << std::setw(4) << itheta[th] << std::setw(4) << ialpha[al] << std::setw(4) << lr << std::setw(15);
425  if (m_fitflag[i][lr][al][th] == 1) {
426  nfitted += 1;
427  for (int p = 0; p < 7; ++p) {
428  if (p != 6) { ofs << std::setprecision(7) << sigma_new[i][lr][al][th][p] << std::setw(15);}
429  if (p == 6) { ofs << std::setprecision(7) << sigma_new[i][lr][al][th][p] << std::endl;}
430  }
431  } else {
432  B2WARNING("Fitting error and old sigma will be used. (Layer " << i << ") (lr = " << lr << ") (al = " << al << ") (th = " << th <<
433  ")");
434  nfailure += 1;
435  int ial_old = 0;
436  int ith_old = 0;
437  for (int k = 0; k < nalpha_old; ++k) {
438  if (ialpha[al] < u_alpha_old[k]) {ial_old = k; break;}
439  }
440  for (int j = 0; j < ntheta_old; ++j) {
441  if (itheta[th] < u_theta_old[j]) {ith_old = j; break;}
442  }
443  for (int p = 0; p < 7; ++p) {
444  if (p != 6) { ofs << std::setprecision(7) << sigma_old[i][lr][ial_old][ith_old][p] << std::setw(15);}
445  if (p == 6) { ofs << std::setprecision(7) << sigma_old[i][lr][ial_old][ith_old][p] << std::endl;}
446  }
447  }
448  }
449  }
450  }
451  }
452  ofs.close();
453  B2RESULT("Number of histogram: " << 56 * 2 * m_nalpha * m_ntheta);
454  B2RESULT("Histos succesfully fitted: " << nfitted);
455  B2RESULT("Histos fit failure: " << nfailure);
456  if (m_useDB) {
457  CDCDatabaseImporter import(0, 0, -1, -1);
458  import.importSigma(m_outputSigmaFileName.c_str());
459  }
460 
461 
462 }
463 void SpaceResolutionCalibration::readSigma()
464 {
465  B2INFO("readSigma");
466  if (m_useDB) {
467  B2INFO("reading sigma from DB");
468  m_sResolFromDB = new DBObjPtr<CDCSpaceResols>;
469  readSigmaFromDB();
470  B2INFO("Number of theta bins from input sigma: " << ntheta_old);
471  } else {
472  B2INFO("Read sigma from text");
473  readSigmaFromText();
474  B2INFO("number of alpha bins from input sigma: " << nalpha_old);
475  }
476 }
477 
478 void SpaceResolutionCalibration::readSigmaFromText()
479 {
480  ifstream ifs;
481  std::string fileName1 = "/data/cdc" + m_sigmafile;
482  std::string fileName = FileSystem::findFile(fileName1);
483  if (fileName == "") {
484  fileName = FileSystem::findFile(m_sigmafile);
485  }
486  if (fileName == "") {
487  cout << "CDCGeometryPar: " << fileName << " not exist!" << endl;
488  } else {
489  cout << "CDCGeometryPar: " << fileName << " exists." << endl;
490  ifs.open(fileName);
491  if (!ifs) cout << "CDCGeometryPar: cannot open " << fileName << " !" << endl;
492  }
493 
494  //read alpha bin info.
495  if (ifs >> nalpha_old) {
496  if (nalpha_old == 0 || nalpha_old > 18) cout << "Fail to read alpha bins !" << endl;
497  } else {
498  cout << "Fail to read alpha bins !" << endl; return;
499  }
500  double alpha0, alpha1, alpha2;
501  for (unsigned short i = 0; i < nalpha_old; ++i) {
502  ifs >> alpha0 >> alpha1 >> alpha2;
503  l_alpha_old[i] = alpha0;
504  u_alpha_old[i] = alpha1;
505  ialpha_old[i] = alpha2;
506  }
507 
508  //read theta bin info.
509  if (ifs >> ntheta_old) {
510  if (ntheta_old == 0 || ntheta_old > 7) cout << "CDCGeometryPar: fail to read theta bins !" << endl;
511  } else {
512  cout << "CDCGeometryPar: fail to read theta bins !" << endl;
513  }
514  double theta0, theta1, theta2;
515  for (unsigned short i = 0; i < ntheta_old; ++i) {
516  ifs >> theta0 >> theta1 >> theta2;
517  l_theta_old[i] = theta0;
518  u_theta_old[i] = theta1;
519  itheta_old[i] = theta2;
520  }
521 
522  unsigned short np = 0;
523  unsigned short iCL, iLR;
524  double theta, alpha;
525  unsigned nRead = 0;
526 
527  ifs >> m_sigmaParamMode_old >> np;
528  double sigma[8];
529  // if (m_sigmaParamMode < 0 || m_sigmaParamMode > 1) cout<<"CDCGeometryPar: invalid sigma-parameterization mode read !"<<endl;
530  //if (m_sigmaParamMode == 1) cout<<"CDCGeometryPar: sigma-parameterization mode=1 not ready yet"<<endl;
531  // if (np <= 0 || np > nSigmaParams) cout<<"CDCGeometryPar: no. of sigma-params. outside limits !"<<endl;
532  const double epsi = 0.1;
533  while (ifs >> iCL) {
534  ifs >> theta >> alpha >> iLR;
535  for (int i = 0; i < np; ++i) {
536  ifs >> sigma[i];
537  }
538  ++nRead;
539 
540  int ith = -99;
541  for (unsigned short i = 0; i < ntheta_old; ++i) {
542  if (fabs(theta - itheta_old[i]) < epsi) {
543  ith = i;
544  break;
545  }
546  }
547  if (ith < 0) cout << "CDCGeometryPar: thetas in sigma.dat are inconsistent !" << endl;
548 
549  int ial = -99;
550  for (unsigned short i = 0; i < nalpha_old; ++i) {
551  if (fabs(alpha - ialpha_old[i]) < epsi) {
552  ial = i;
553  break;
554  }
555  }
556  if (ial < 0) cout << "CDCGeometryPar: alphas in sigma.dat are inconsistent !" << endl;
557 
558  for (int i = 0; i < np; ++i) {
559  sigma_old[iCL][iLR][ial][ith][i] = sigma[i];
560  }
561  } //end of while loop
562  ifs.close();
563 }
564 void SpaceResolutionCalibration::readSigmaFromDB()
565 {
566  typedef std::array<float, 3> array3;
567  // std::cout <<"setSResol called" << std::endl;
568  nalpha_old = (*m_sResolFromDB)->getNoOfAlphaBins();
569  double rad2deg = 180 / M_PI;
570  for (unsigned short i = 0; i < nalpha_old; ++i) {
571  // m_alphaPoints[i] = (*dbXT_old)->getAlphaPoint(i);
572  array3 alpha = (*m_sResolFromDB)->getAlphaBin(i);
573  l_alpha_old[i] = alpha[0] * rad2deg;
574  u_alpha_old[i] = alpha[1] * rad2deg;
575  ialpha_old[i] = alpha[2] * rad2deg;
576  // std::cout << m_alphaPoints[i]*180./M_PI << std::endl;
577  }
578 
579  ntheta_old = (*m_sResolFromDB)->getNoOfThetaBins();
580  B2INFO("Ntheta: " << ntheta_old);
581  for (unsigned short i = 0; i < ntheta_old; ++i) {
582  // m_thetaPoints[i] = (*dbXT_old).getThetaPoint(i);
583  array3 theta = (*m_sResolFromDB)->getThetaBin(i);
584  l_theta_old[i] = theta[0] * rad2deg;
585  u_theta_old[i] = theta[1] * rad2deg;
586  itheta_old[i] = theta[2] * rad2deg;
587  }
588  m_sigmaParamMode_old = (*m_sResolFromDB)->getSigmaParamMode();
589 
590  for (unsigned short iCL = 0; iCL < MAX_N_SLAYERS; ++iCL) {
591  for (unsigned short iLR = 0; iLR < 2; ++iLR) {
592  for (unsigned short iA = 0; iA < nalpha_old; ++iA) {
593  for (unsigned short iT = 0; iT < ntheta_old; ++iT) {
594  const std::vector<float> params = (*m_sResolFromDB)->getSigmaParams(iCL, iLR, iA, iT);
595  unsigned short np = params.size();
596  // std::cout <<"np4sigma= " << np << std::endl;
597  for (unsigned short i = 0; i < np; ++i) {
598  sigma_old[iCL][iLR][iA][iT][i] = params[i];
599  }
600  }
601  }
602  }
603  }
604 
605 }
606 void SpaceResolutionCalibration::readProfile()
607 {
608  B2INFO("readProfile");
609  /*Read profile for xt*/
610  if (m_useProfileFromInputSigma) {
611  B2INFO("use Sigma bining from input Sigma");
612  m_nalpha = nalpha_old;
613  m_ntheta = ntheta_old;
614  B2INFO("Number of alpha bins: " << m_nalpha);
615  for (int i = 0; i < m_nalpha; ++i) {
616  l_alpha[i] = l_alpha_old[i]; u_alpha[i] = u_alpha_old[i]; ialpha[i] = ialpha_old[i];
617  B2INFO("" << i << " | " << l_alpha[i] << " " << u_alpha[i] << " " << ialpha[i]);
618  }
619  B2INFO("Number of theta bins: " << m_ntheta);
620  for (int i = 0; i < m_ntheta; ++i) {
621  l_theta[i] = l_theta_old[i]; u_theta[i] = u_theta_old[i]; itheta[i] = itheta_old[i];
622  B2INFO("" << i << " |" << l_theta[i] << " " << u_theta[i] << " " << itheta[i]);
623  }
624  } else {
625  B2INFO("use Sigma bining from profile file");
626  ifstream proxt(m_ProfileFileName.c_str());
627  if (!proxt) {
628  B2FATAL("file not found: " << m_ProfileFileName);
629  }
630  double dumy1, dumy2, dumy3;
631  proxt >> m_nalpha;
632  B2INFO("Number of alpha bins: " << m_nalpha);
633  if (m_nalpha > Max_nalpha) {B2FATAL("number of alpha bin excess limit; please increse uplimit: " << m_nalpha << " > " << Max_nalpha);}
634  for (int i = 0; i < m_nalpha; ++i) {
635  proxt >> dumy1 >> dumy2 >> dumy3;
636  l_alpha[i] = dumy1; u_alpha[i] = dumy2; ialpha[i] = dumy3;
637  B2INFO("" << i << " | " << l_alpha[i] << " " << u_alpha[i] << " " << ialpha[i]);
638  }
639  proxt >> m_ntheta;
640  B2INFO("Number of theta bins: " << m_ntheta);
641  if (m_ntheta > Max_ntheta) {B2FATAL("number of theta bin excess limit; please increse uplimit: " << m_ntheta << " > " << Max_ntheta);}
642  for (int i = 0; i < m_ntheta; ++i) {
643  proxt >> dumy1 >> dumy2 >> dumy3;
644  l_theta[i] = dumy1; u_theta[i] = dumy2; itheta[i] = dumy3;
645  B2INFO("" << i << " |" << l_theta[i] << " " << u_theta[i] << " " << itheta[i]);
646  }
647  }
648  B2INFO("Finish asssign sigma bining");
649 }
Belle2::CDC::SpaceResolutionCalibration::hist_u
TH2F * hist_u[56][2][Max_nalpha][Max_ntheta]
2D histogram of unbiased residual
Definition: SpaceResolutionCalibration.h:102
Belle2::CDC::SpaceResolutionCalibration::Max_np
static const unsigned short Max_np
Maximum number of point =1/binwidth.
Definition: SpaceResolutionCalibration.h:84
Belle2::CDC::SpaceResolutionCalibration::m_useDB
bool m_useDB
use db or text mode
Definition: SpaceResolutionCalibration.h:92
Belle2::CDC::SpaceResolutionCalibration::createHisto
virtual void createHisto()
create histogram
Definition: SpaceResolutionCalibration.cc:31
Belle2::CDC::SpaceResolutionCalibration::itheta
double itheta[7]
represented alphas of theta bins.
Definition: SpaceResolutionCalibration.h:126
Belle2::CDC::SpaceResolutionCalibration::m_useProfileFromInputSigma
bool m_useProfileFromInputSigma
Use binning from old sigma or new one form input.
Definition: SpaceResolutionCalibration.h:93
Belle2::CDC::SpaceResolutionCalibration::hb_m
TH1F * hb_m[56][2][Max_nalpha][Max_ntheta]
mean histogram of unbiased residual
Definition: SpaceResolutionCalibration.h:105
Belle2::CDC::SpaceResolutionCalibration::hu_s
TH1F * hu_s[56][2][Max_nalpha][Max_ntheta]
sigma histogram of biased residual
Definition: SpaceResolutionCalibration.h:104
Belle2::CDC::SpaceResolutionCalibration::write
virtual void write()
save calibration, in text file or db
Definition: SpaceResolutionCalibration.cc:397
Belle2::CDC::SpaceResolutionCalibration::hist_b
TH2F * hist_b[56][2][Max_nalpha][Max_ntheta]
2D histogram of biased residual
Definition: SpaceResolutionCalibration.h:101
Belle2::CDC::SpaceResolutionCalibration::l_theta_old
double l_theta_old[7]
Lower boundays of theta bins from input.
Definition: SpaceResolutionCalibration.h:134
Belle2::CDC::SpaceResolutionCalibration::m_ndfmin
double m_ndfmin
Minimum NDF
Definition: SpaceResolutionCalibration.h:86
Belle2::CDC::SpaceResolutionCalibration::Max_nalpha
static const int Max_nalpha
Maximum alpha bin.
Definition: SpaceResolutionCalibration.h:82
Belle2::CDC::SpaceResolutionCalibration::sigma_new
double sigma_new[56][2][18][7][8]
new sigma prameters.
Definition: SpaceResolutionCalibration.h:97
Belle2::CDC::SpaceResolutionCalibration::ntheta_old
int ntheta_old
number of theta bins from input
Definition: SpaceResolutionCalibration.h:130
Belle2::CDC::SpaceResolutionCalibration::m_fitflag
int m_fitflag[56][2][Max_nalpha][Max_ntheta]
Fit flag; 1:OK ; 0:error.
Definition: SpaceResolutionCalibration.h:107
Belle2::CDC::SpaceResolutionCalibration::nalpha_old
int nalpha_old
number of alpha bins from input
Definition: SpaceResolutionCalibration.h:129
Belle2::CDC::SpaceResolutionCalibration::u_alpha_old
double u_alpha_old[18]
Upper boundays of alpha bins from input.
Definition: SpaceResolutionCalibration.h:132
Belle2::CDC::SpaceResolutionCalibration::readProfile
virtual void readProfile()
read sigma bining (alpha, theta bining)
Definition: SpaceResolutionCalibration.cc:606
Belle2::CDC::SpaceResolutionCalibration::l_alpha_old
double l_alpha_old[18]
Lower boundays of alpha bins from input.
Definition: SpaceResolutionCalibration.h:131
Belle2::CDC::SpaceResolutionCalibration::m_outputSigmaFileName
std::string m_outputSigmaFileName
Output sigma file name.
Definition: SpaceResolutionCalibration.h:109
Belle2::CDC::SpaceResolutionCalibration::hu_m
TH1F * hu_m[56][2][Max_nalpha][Max_ntheta]
mean histogram biased residual
Definition: SpaceResolutionCalibration.h:103
Belle2::CDC::SpaceResolutionCalibration::l_theta
double l_theta[7]
Lower boundays of theta bins.
Definition: SpaceResolutionCalibration.h:124
Belle2::CDC::SpaceResolutionCalibration::l_alpha
double l_alpha[18]
Lower boundays of alpha bins.
Definition: SpaceResolutionCalibration.h:121
Belle2::DBObjPtr
Class for accessing objects in the database.
Definition: DBObjPtr.h:31
Belle2::CDC::SpaceResolutionCalibration::m_sResolFromDB
DBObjPtr< CDCSpaceResols > * m_sResolFromDB
Database for sigma.
Definition: SpaceResolutionCalibration.h:112
Belle2::CDC::SpaceResolutionCalibration::storeHisto
virtual void storeHisto()
store histogram
Definition: SpaceResolutionCalibration.cc:364
Belle2::CDC::SpaceResolutionCalibration::gr
TGraphErrors * gr[56][2][18][7]
sigma graph.
Definition: SpaceResolutionCalibration.h:100
Belle2::CDC::SpaceResolutionCalibration::itheta_old
double itheta_old[7]
represented alphas of theta bins from input.
Definition: SpaceResolutionCalibration.h:136
Belle2::CDC::SpaceResolutionCalibration::u_theta_old
double u_theta_old[7]
Upper boundays of theta bins from input.
Definition: SpaceResolutionCalibration.h:135
Belle2::CDC::SpaceResolutionCalibration::sigma_old
double sigma_old[56][2][18][7][8]
old sigma prameters.
Definition: SpaceResolutionCalibration.h:96
Belle2::CDC::SpaceResolutionCalibration::m_inputRootFileNames
std::string m_inputRootFileNames
Input root file names.
Definition: SpaceResolutionCalibration.h:110
Belle2::CDC::SpaceResolutionCalibration::gfit
TGraphErrors * gfit[56][2][18][7]
sigma*sigma graph for fit
Definition: SpaceResolutionCalibration.h:99
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::CDC::SpaceResolutionCalibration::u_theta
double u_theta[7]
Upper boundays of theta bins.
Definition: SpaceResolutionCalibration.h:125
Belle2::CDC::SpaceResolutionCalibration::ialpha_old
double ialpha_old[18]
represented alphas of alpha bins from input.
Definition: SpaceResolutionCalibration.h:133
Belle2::CDC::SpaceResolutionCalibration::m_ntheta
int m_ntheta
number of theta bins
Definition: SpaceResolutionCalibration.h:120
Belle2::CDC::SpaceResolutionCalibration::m_binWidth
double m_binWidth
width of each bin, unit cm
Definition: SpaceResolutionCalibration.h:88
Belle2::CDC::SpaceResolutionCalibration::m_ProfileFileName
std::string m_ProfileFileName
Profile file name.
Definition: SpaceResolutionCalibration.h:111
Belle2::CDC::SpaceResolutionCalibration::Max_ntheta
static const int Max_ntheta
maximum theta bin
Definition: SpaceResolutionCalibration.h:83
Belle2::CDC::SpaceResolutionCalibration::readSigma
virtual void readSigma()
read sigma from previous calibration, (input sigma)
Definition: SpaceResolutionCalibration.cc:463
Belle2::CDC::SpaceResolutionCalibration::m_sigmafile
std::string m_sigmafile
Sigma file name, for text mode.
Definition: SpaceResolutionCalibration.h:113
Belle2::CDC::SpaceResolutionCalibration::m_Pvalmin
double m_Pvalmin
Minimum Prob(chi2) of track.
Definition: SpaceResolutionCalibration.h:87
Belle2::CDC::SpaceResolutionCalibration::ialpha
double ialpha[18]
represented alphas of alpha bins.
Definition: SpaceResolutionCalibration.h:123
Belle2::CDC::SpaceResolutionCalibration::calibrate
virtual bool calibrate()
Run algo on data.
Definition: SpaceResolutionCalibration.cc:278
Belle2::CDC::SpaceResolutionCalibration::m_sigmaParamMode_old
unsigned short m_sigmaParamMode_old
sigma mode from input.
Definition: SpaceResolutionCalibration.h:138
Belle2::CDC::SpaceResolutionCalibration::readSigmaFromDB
virtual void readSigmaFromDB()
read sigma from DB
Definition: SpaceResolutionCalibration.cc:564
Belle2::CDC::SpaceResolutionCalibration::hb_s
TH1F * hb_s[56][2][Max_nalpha][Max_ntheta]
sigma histogram of ubiased residual
Definition: SpaceResolutionCalibration.h:106
Belle2::FileSystem::findFile
static std::string findFile(const std::string &path, bool silent=false)
Search for given file or directory in local or central release directory, and return absolute path if...
Definition: FileSystem.cc:147
Belle2::CDCDatabaseImporter::importSigma
void importSigma(std::string fileName)
Import sigma table to the database.
Belle2::CDC::SpaceResolutionCalibration::m_nalpha
int m_nalpha
number of alpha bins
Definition: SpaceResolutionCalibration.h:119
Belle2::CDC::SpaceResolutionCalibration::readSigmaFromText
virtual void readSigmaFromText()
read sigma from text file
Definition: SpaceResolutionCalibration.cc:478
Belle2::CDCDatabaseImporter
CDC database importer.
Definition: CDCDatabaseImporter.h:32
Belle2::CDC::SpaceResolutionCalibration::getUpperBoundaryForFit
double getUpperBoundaryForFit(TGraphErrors *graph)
search max point at boundary region
Definition: SpaceResolutionCalibration.h:144
Belle2::CDC::SpaceResolutionCalibration::u_alpha
double u_alpha[18]
Upper boundays of alpha bins.
Definition: SpaceResolutionCalibration.h:122