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