12 #include <cdc/calibration/SpaceResolutionCalibrationAlgorithm.h>
13 #include <cdc/geometry/CDCGeometryPar.h>
14 #include <cdc/dbobjects/CDCSpaceResols.h>
16 #include <framework/database/DBObjPtr.h>
17 #include <framework/logging/Logger.h>
30 typedef std::array<float, 3> array3;
31 SpaceResolutionCalibrationAlgorithm::SpaceResolutionCalibrationAlgorithm() :
35 " -------------------------- Position Resolution Calibration Algorithm -------------------------\n"
40 B2INFO(
"Creating histograms");
45 for (
int i = 0; i < 50; ++i) {
46 yb.push_back(-0.07 + i * (0.14 / 50));
48 for (
int i = 0; i < 50; ++i) {
49 yu.push_back(-0.08 + i * (0.16 / 50));
55 for (
int i = 1; i < np; ++i) {
59 for (
int il = 0; il < 56; ++il) {
60 for (
int lr = 0; lr < 2; ++lr) {
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));
75 auto tree = getObjectPtr<TTree>(
"tree");
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);
99 std::vector<TString> list_vars = {
"lay",
"ndf",
"Pval",
"x_u",
"x_b",
"x_mea",
"weight",
"alpha",
"theta"};
100 tree->SetBranchStatus(
"*", 0);
102 for (TString brname : list_vars) {
103 tree->SetBranchStatus(brname, 1);
107 const Long64_t nEntries = tree->GetEntries();
108 B2INFO(
"Number of entries: " << nEntries);
111 for (Long64_t i = 0; i < nEntries; ++i) {
113 if (std::fabs(x_b) < 0.02 || std::fabs(x_u) < 0.02)
continue;
130 int ilr = x_u > 0 ? 1 : 0;
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);
137 absRes_u = fabs(x_mea) - fabs(x_u);
138 absRes_b = fabs(x_mea) - fabs(x_b);
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);
146 B2INFO(
"Start to obtain the biased and unbiased sigmas");
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);
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;
161 ofstream ofss(
"IntReso.dat");
165 for (
int il = 0; il < 56; ++il) {
166 for (
int lr = 0; lr < 2; ++lr) {
170 B2DEBUG(21,
"layer-lr-al-th " << il <<
" - " << lr <<
" - " << al <<
" - " << th);
171 if (
m_hBiased[il][lr][al][th]->GetEntries() < 5000) {
176 auto* proYb =
m_hBiased[il][lr][al][th]->ProjectionY();
177 auto* proYu =
m_hUnbiased[il][lr][al][th]->ProjectionY();
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);
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());
193 B2DEBUG(21,
"Nentries: " <<
m_hBiased[il][lr][al][th]->GetEntries());
194 m_hBiased[il][lr][al][th]->SetDirectory(0);
200 m_hBiased[il][lr][al][th]->FitSlicesY(g0b, firstbin, ib1, minEntry);
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,
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",
214 m_hBiased[il][lr][al][th]->FitSlicesY(gb, ib1 + 1, np, minEntry);
216 m_hMeanBiased[il][lr][al][th]->Add((TH1F*)gDirectory->Get(Form(
"hb_%d_%d_%d_%d_1", il, lr, al, th)));
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());
224 m_hUnbiased[il][lr][al][th]->FitSlicesY(g0u, firstbin, ib1, minEntry);
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",
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",
238 m_hUnbiased[il][lr][al][th]->FitSlicesY(gu, ib1 + 1, np, minEntry);
240 m_hMeanUnbiased[il][lr][al][th]->Add((TH1F*)gDirectory->Get(Form(
"hu_%d_%d_%d_%d_1", il, lr, al, th)));
242 m_hSigmaUnbiased[il][lr][al][th]->Add((TH1F*)gDirectory->Get(Form(
"hu_%d_%d_%d_%d_2", il, lr, al, th)));
244 B2WARNING(
"sliced histo not found");
260 if (
m_hSigmaBiased[il][lr][al][th]->GetBinContent(j) == 0)
continue;
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;
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 <<
" " <<
282 if (xl.size() < 7 || xl.size() >
Max_np) {
284 B2WARNING(
"number of element might out of range");
continue;
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));
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));
302 gDirectory->Delete(
"hu_%d_%d_%d_%d_0");
314 B2INFO(
"Start calibration");
315 gPrintViaErrorHandler =
true;
317 gErrorIgnoreLevel = 3001;
320 B2INFO(
"ExpRun used for DB Geometry : " << exprun.first <<
" " << exprun.second);
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);
333 for (
int i = 0; i < 56; ++i) {
334 for (
int lr = 0; lr < 2; ++lr) {
337 if (!
m_gFit[i][lr][al][th])
continue;
341 B2DEBUG(199,
"xmax for fitting: " << upFit);
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);
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]);
355 for (
int j = 0; j < 10; j++) {
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") {
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);
372 B2DEBUG(21,
"Prob of fit: " << func->GetProb());
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);
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]);
404 int nFitCompleted = 0;
405 for (
int l = 0; l < 56; ++l) {
406 for (
int lr = 0; lr < 2; ++lr) {
417 if (
static_cast<double>(nFitCompleted) / nTotal <
m_threshold) {
418 B2WARNING(
"Less than " <<
m_threshold * 100 <<
" % of Sigmas were fitted.");
427 B2INFO(
"saving histograms");
429 TFile* ff =
new TFile(
m_histName.c_str(),
"RECREATE");
430 TDirectory* top = gDirectory;
431 TDirectory* Direct[56];
433 auto hNDF = getObjectPtr<TH1F>(
"hNDF");
434 auto hPval = getObjectPtr<TH1F>(
"hPval");
435 auto hEvtT0 = getObjectPtr<TH1F>(
"hEventT0");
437 if (hNDF && hPval && hEvtT0) {
444 for (
int il = 0; il < 56; ++il) {
446 Direct[il] = gDirectory->mkdir(Form(
"lay_%d", il));
449 for (
int lr = 0; lr < 2; ++lr) {
452 if (!
m_graph[il][lr][al][th])
continue;
453 if (!
m_gFit[il][lr][al][th])
continue;
461 m_graph[il][lr][al][th]->Write();
462 m_gFit[il][lr][al][th]->Write();
469 B2INFO(
"Finish store histogram");
474 B2INFO(
"Writing calibrated sigma's");
480 const float deg2rad = M_PI / 180.0;
483 std::array<float, 3> alpha3 = {
m_lowerAlpha[i]* deg2rad,
492 std::array<float, 3> theta3 = {
m_lowerTheta[i]* deg2rad,
502 for (
int iCL = 0; iCL < 56; ++iCL) {
503 for (
int iLR = 1; iLR >= 0; --iLR) {
504 std::vector<float> sgbuff;
507 for (
int i = 0; i < 7; ++i) {
508 sgbuff.push_back(
m_sigma[iCL][iLR][ialpha][itheta][i]);
514 for (
int i = 0; i < 7; ++i) {
515 sgbuff.push_back(
m_sigmaPost[iCL][iLR][ialpha][itheta][i]);
531 B2RESULT(
"Histos succesfully fitted: " << nfitted);
532 B2RESULT(
"Histos fit failure: " << nfailure);
539 B2INFO(
"Prepare calibration of space resolution");
541 const double rad2deg = 180 / M_PI;
550 array3 alpha = dbSigma->getAlphaBin(i);
558 array3 theta = dbSigma->getThetaBin(i);
565 for (
unsigned short iCL = 0; iCL < 56; ++iCL) {
566 for (
unsigned short iLR = 0; iLR < 2; ++iLR) {
569 const std::vector<float> params = dbSigma->getSigmaParams(iCL, iLR, iA, iT);
570 unsigned short np = params.size();
572 for (
unsigned short i = 0; i < np; ++i) {