11 #include <reconstruction/calibration/CDCDedxRunGainAlgorithm.h>
45 B2INFO(
"Preparing dE/dx calibration for CDC dE/dx run gain");
52 auto ttree = getObjectPtr<TTree>(
"tree");
56 int run = 0, runtemp = 0;
58 ttree->SetBranchAddress(
"dedx", &dedx);
59 ttree->SetBranchAddress(
"run", &run);
64 for (
int i = 0; i < ttree->GetEntries(); ++i) {
66 if (dedx <= 0)
continue;
68 hDedxAbs->Fill(dedx * ExistingRG);
73 if (hDedx->Integral() < 1000) {
74 B2INFO(
"Not enough data for run (" << runtemp <<
") so far only " << hDedx->Integral() <<
" therefore adding next run");
75 fsrun += Form(
"%d_", runtemp);
81 B2INFO(
"Total track for this run: " << runtemp <<
" = " << hDedx->Integral());
82 fsrun += Form(
"%d", runtemp);
83 hDedx->SetName(Form(
"%s_%s", hDedx->GetName(),
fsrun.Data()));
84 hDedxAbs->SetName(Form(
"%s_%s", hDedxAbs->GetName(),
fsrun.Data()));
87 double RGabs_Mean = 1.0;
89 if (status !=
"FitOK") {
90 hDedxAbs->SetTitle(Form(
"dE/dx abs, Run %s (%s)",
fsrun.Data(), status.Data()));
92 RGabs_Mean = hDedxAbs->GetFunction(
"gaus")->GetParameter(1);
93 double RGabs_MeanErr = hDedxAbs->GetFunction(
"gaus")->GetParError(1);
94 double RGabs_Sigma = hDedxAbs->GetFunction(
"gaus")->GetParameter(2);
95 hDedxAbs->SetTitle(Form(
"dE/dx abs, Run %s (%s), #mu_{fit}: %0.04f#pm%0.04f,, #sigma_{fit}: %0.04f",
fsrun.Data(), status.Data(),
96 RGabs_Mean, RGabs_MeanErr, RGabs_Sigma));
100 double RGrel_Mean = 1.0, RGrel_MeanErr = -99.0, RGrel_Sigma = -99.0;
102 if (status !=
"FitOK") {
103 RGrel_Mean = hDedx->GetMean();
104 hDedx->SetTitle(Form(
"dE/dx abs, Run %s (%s)",
fsrun.Data(), status.Data()));
106 RGrel_Mean = hDedx->GetFunction(
"gaus")->GetParameter(1);
107 RGrel_MeanErr = hDedx->GetFunction(
"gaus")->GetParError(1);
108 RGrel_Sigma = hDedx->GetFunction(
"gaus")->GetParameter(2);
109 hDedx->SetTitle(Form(
"dE/dx rel, Run %s (%s), #mu_{fit}: %0.04f#pm%0.04f,, #sigma_{fit}: %0.04f",
fsrun.Data(), status.Data(),
110 RGrel_Mean, RGrel_MeanErr, RGrel_Sigma));
113 if (RGrel_Mean <= 0.0)RGrel_Mean = 1.0;
118 std::ofstream fileRGout;
119 fileRGout.open(Form(
"fRG_Constants_%s.txt",
fsrun.Data()));
120 fileRGout << expRun.second <<
" " << ExistingRG <<
" " << RGrel_Mean <<
" " << RGrel_MeanErr <<
" " << RGrel_Sigma <<
"\n";
123 TCanvas* ctmp =
new TCanvas(
"tmp",
"tmp", 1000, 500);
124 ctmp->SetBatch(kTRUE);
128 hDedx->SetFillColorAlpha(kYellow, 0.30);
131 TLine* tl =
new TLine();
132 tl->SetLineColor(kRed);
133 tl->SetX1(RGrel_Mean); tl->SetX2(RGrel_Mean);
134 tl->SetY1(0); tl->SetY2(hDedx->GetMaximum());
135 tl->DrawClone(
"same");
138 hDedxAbs->SetFillColorAlpha(kBlue, 0.10);
139 hDedxAbs->DrawCopy(
"");
141 TLine* tlAbs =
new TLine();
142 tlAbs->SetLineColor(kRed);
143 tlAbs->SetX1(RGabs_Mean); tlAbs->SetX2(RGabs_Mean);
144 tlAbs->SetY1(0); tlAbs->SetY2(hDedxAbs->GetMaximum());
145 tlAbs->DrawClone(
"same");
147 ctmp->Print(Form(
"cdcdedx_rungain_fits_%s_%s.pdf",
fsrun.Data(), status.Data()));
149 TCanvas* cstats =
new TCanvas(
"cstats",
"cstats", 1000, 500);
150 cstats->SetBatch(kTRUE);
151 cstats->Divide(2, 1);
153 auto hestats = getObjectPtr<TH1I>(
"hestats");
154 if (hestats)hestats->DrawCopy(
"");
156 auto htstats = getObjectPtr<TH1I>(
"htstats");
157 if (htstats)htstats->DrawCopy(
"");
158 cstats->Print(Form(
"cdcdedx_rungain_stats_%s.pdf",
fsrun.Data()));
166 B2INFO(
"Saving new rung for (Exp, Run) : (" << expRun.first <<
"," << expRun.second <<
")");
179 B2INFO(
"--> RunGain: Previous = " << ExistingRG <<
", Relative = " << RGrel_Mean <<
", Adjustment = " <<
fAdjust <<
180 ", Merged (saved) = " << RGrel_Mean *
182 RGrel_Mean *= ExistingRG;
185 B2INFO(
"--> RunGain: Previous = " << ExistingRG <<
", Relative (saved) = " << RGrel_Mean);
196 double histmean = temphist->GetMean();
197 double histrms = temphist->GetRMS();
198 temphist->GetXaxis()->SetRangeUser(histmean - 5.0 * histrms, histmean + 5.0 * histrms);
200 int fs = temphist->Fit(
"gaus",
"Q0");
202 B2INFO(Form(
"\tFit (round 1) for hist (%s) failed (status = %d)", temphist->GetName(), fs));
203 status =
"FitFailed";
206 double mean = temphist->GetFunction(
"gaus")->GetParameter(1);
207 double width = temphist->GetFunction(
"gaus")->GetParameter(2);
208 temphist->GetXaxis()->SetRangeUser(mean - 5.0 * width, mean + 5.0 * width);
209 fs = temphist->Fit(
"gaus",
"QR",
"", mean -
fSigLim * width, mean +
fSigLim * width);
211 B2INFO(Form(
"\tFit (round 2) for hist (%s) failed (status = %d)", temphist->GetName(), fs));
212 status =
"FitFailed";
215 temphist->GetXaxis()->SetRangeUser(mean - 5.0 * width, mean + 5.0 * width);
216 B2INFO(Form(
"\tFit for hist (%s) sucessfull (status = %d)", temphist->GetName(), fs));