20 #include <reconstruction/calibration/CDCDedx1DCellAlgorithm.h>
33 feaLE(-TMath::Pi() / 2),
34 feaUE(+TMath::Pi() / 2),
41 setDescription(
"A calibration algorithm for the CDC dE/dx entrance angle cleanup correction");
52 auto ttree = getObjectPtr<TTree>(
"tree");
55 std::vector<double>* dedxhit = 0, *enta = 0;
56 std::vector<int>* layer = 0;
58 ttree->SetBranchAddress(
"dedxhit", &dedxhit);
59 ttree->SetBranchAddress(
"layer", &layer);
60 ttree->SetBranchAddress(
"entaRS", &enta);
65 std::vector<int> globalbins;
66 for (
int ibin = 0; ibin <
fnEntaBinG; ibin++)globalbins.push_back(ibin);
77 std::vector<TH1D*> hILdEdxhitInEntaBin(
fnEntaBinL, 0);
78 std::vector<TH1D*> hOLdEdxhitInEntaBin(
fnEntaBinL, 0);
80 Double_t ifeaLE = 0, ifeaUE = 0;
89 ifeaUE = ifeaLE +
feaBS;
92 hILdEdxhitInEntaBin[iea] =
new TH1D(Form(
"hILdEdxhitInEntaBin%d", iea),
"bla-bla", 250, 0., 5.);
93 hILdEdxhitInEntaBin[iea]->SetTitle(Form(
"IL: dedxhit in EntA = (%0.03f to %0.03f)", ifeaLE, ifeaUE));
94 hILdEdxhitInEntaBin[iea]->GetXaxis()->SetTitle(
"dedxhits in Inner Layer");
95 hILdEdxhitInEntaBin[iea]->GetYaxis()->SetTitle(
"Entries");
97 hOLdEdxhitInEntaBin[iea] =
new TH1D(Form(
"hOLdEdxhitInEntaBin%d", iea),
"bla-bla", 250, 0., 5.);
98 hOLdEdxhitInEntaBin[iea]->SetTitle(Form(
"OL: dedxhit in EntA = (%0.03f to %0.03f)", ifeaLE, ifeaUE));
99 hOLdEdxhitInEntaBin[iea]->GetXaxis()->SetTitle(
"dedxhits in Outer Layer");
100 hOLdEdxhitInEntaBin[iea]->GetYaxis()->SetTitle(
"Entries");
105 hILEntaG->GetXaxis()->SetTitle(
"Entrance angle (#theta)");
106 hILEntaG->GetYaxis()->SetTitle(
"Entries");
109 hOLEntaG->GetXaxis()->SetTitle(
"Entrance angle (#theta)");
110 hOLEntaG->GetYaxis()->SetTitle(
"Entries");
115 TH1D* hILEntaL =
new TH1D(
"hILEntaL",
"EntA: Inner Layer (assym bin)",
fnEntaBinL, RmapEntaValue);
116 hILEntaL->GetXaxis()->SetTitle(
"Entrance angle (#theta)");
117 hILEntaL->GetYaxis()->SetTitle(
"Entries");
119 TH1D* hOLEntaL =
new TH1D(
"hOLEntaL",
"EntA: Outer Layer (assym bin)",
fnEntaBinL, RmapEntaValue);
120 hOLEntaL->GetYaxis()->SetTitle(
"Entries");
121 hOLEntaL->GetXaxis()->SetTitle(
"Entrance angle (#theta)");
123 TH1D* hILdEdx_all =
new TH1D(
"hILdEdx_all",
"", 250, 0., 5.);
124 TH1D* hOLdEdx_all =
new TH1D(
"hOLdEdx_all",
"", 250, 0., 5.);
127 for (
int i = 0; i < ttree->GetEntries(); ++i) {
131 for (
unsigned int j = 0; j < dedxhit->size(); ++j) {
133 if (dedxhit->at(j) == 0)
continue;
135 Double_t ieaHit = enta->at(j);
136 if (ieaHit < -TMath::Pi() / 2.0) ieaHit += TMath::Pi() / 2.0;
137 else if (ieaHit > TMath::Pi() / 2.0) ieaHit -= TMath::Pi() / 2.0;
138 if (abs(ieaHit) > TMath::Pi() / 2.0)
continue;
146 if (layer->at(j) < 8) {
147 hILEntaG->Fill(ieaHit);
149 hILdEdx_all->Fill(dedxhit->at(j));
150 hILdEdxhitInEntaBin[ibinEA]->Fill(dedxhit->at(j));
152 hOLEntaG->Fill(ieaHit);
154 hOLdEdx_all->Fill(dedxhit->at(j));
155 hOLdEdxhitInEntaBin[ibinEA]->Fill(dedxhit->at(j));
161 TCanvas* ctmpde =
new TCanvas(
"hEntaDist",
"Enta distributions", 400, 400);
163 ctmpde->SetCanvasSize(800, 400);
164 ctmpde->Divide(2, 1);
165 ctmpde->cd(1); gPad->SetLogy(); hOLEntaG->SetMarkerColor(kBlue); hOLEntaG->Draw(
""); hILEntaG->Draw(
"same");
166 ctmpde->cd(2); gPad->SetLogy(); hOLEntaL->SetMarkerColor(kBlue); hOLEntaL->Draw(
""); hILEntaL->Draw(
"same");
168 hOLEntaG->Draw(
""); hILEntaG->Draw(
"same");
170 ctmpde->SaveAs(Form(
"hEntaDistributions_%s.pdf",
fSetPrefix.data()));
174 double InsumPer5 = 0.0, InsumPer75 = 0.0;
175 double OutsumPer5 = 0.0, OutsumPer75 = 0.0;
176 double InLayInt = hILdEdx_all->Integral();
177 double OutLayInt = hOLdEdx_all->Integral();
179 Int_t lBinInLayer = 1, hBinInLayer = 1;
180 Int_t lBinOutLayer = 1, hBinOutLayer = 1;
182 for (
int ibin = 1; ibin <= hILdEdx_all->GetNbinsX(); ibin++) {
184 if (InsumPer5 <= 0.05 * InLayInt) {
185 InsumPer5 += hILdEdx_all->GetBinContent(ibin);
189 if (InsumPer75 <= 0.75 * InLayInt) {
190 InsumPer75 += hILdEdx_all->GetBinContent(ibin);
194 if (OutsumPer5 <= 0.05 * OutLayInt) {
195 OutsumPer5 += hOLdEdx_all->GetBinContent(ibin);
199 if (OutsumPer75 <= 0.75 * OutLayInt) {
200 OutsumPer75 += hOLdEdx_all->GetBinContent(ibin);
207 TCanvas* ctem =
new TCanvas(
"Layerhisto",
"Inner and Outer Histo", 1200, 600);
210 hILdEdx_all->SetMarkerColor(kBlue);
211 hILdEdx_all->Draw(
"histo");
213 TLine* tlInF =
new TLine();
214 tlInF->SetLineColor(kBlue);
215 tlInF->SetX1(InsumPer5); tlInF->SetX2(InsumPer5);
216 tlInF->SetY1(0); tlInF->SetY2(hILdEdx_all->GetMaximum());
217 tlInF->DrawClone(
"same");
219 TLine* tlInL =
new TLine();
220 tlInL->SetLineColor(kBlue);
221 tlInL->SetX1(InsumPer75); tlInL->SetX2(InsumPer75);
222 tlInL->SetY1(0); tlInL->SetY2(hILdEdx_all->GetMaximum());
223 tlInL->DrawClone(
"same");
226 hOLdEdx_all->Draw(
"histo");
227 hOLdEdx_all->SetMarkerColor(kRed);
229 TLine* tlOutF =
new TLine();
230 tlOutF->SetLineColor(kBlue);
231 tlOutF->SetX1(OutsumPer5); tlOutF->SetX2(OutsumPer5);
232 tlOutF->SetY1(0); tlOutF->SetY2(hOLdEdx_all->GetMaximum());
233 tlOutF->DrawClone(
"same");
235 TLine* tlOutL =
new TLine();
236 tlOutL->SetLineColor(kBlue);
237 tlOutL->SetX1(OutsumPer75); tlOutL->SetX2(OutsumPer75);
238 tlOutL->SetY1(0); tlOutL->SetY2(hOLdEdx_all->GetMaximum());
239 tlOutL->DrawClone(
"same");
240 ctem->SaveAs(Form(
"Layerhistodedxhit_OneDCorr_%s.pdf",
fSetPrefix.data()));
244 TCanvas* ctmp =
new TCanvas(
"tmp",
"tmp", 1200, 1200);
246 std::stringstream psname; psname << Form(
"dedx_1dcell_%s.pdf[",
fSetPrefix.data());
247 TLine* tl =
new TLine();
248 tl->SetLineColor(kRed);
250 TBox* tb =
new TBox();
253 ctmp->Print(psname.str().c_str());
254 psname.str(
""); psname << Form(
"dedx_1dcell_%s.pdf",
fSetPrefix.data());
258 std::vector<std::vector<double>> onedcors;
259 std::vector<double> onedcorIorOL, onedcorIorOLtemp;
262 hILEntaConst->GetXaxis()->SetTitle(
"Entrance angle (#theta)");
263 hILEntaConst->GetYaxis()->SetTitle(
"Constant");
266 hOLEntaConst->GetXaxis()->SetTitle(
"Entrance angle (#theta)");
267 hOLEntaConst->GetYaxis()->SetTitle(
"Constant");
269 for (
int iIOLayer = 0; iIOLayer <= 1; iIOLayer++) {
271 int startfrom = 1, endat = 1;
274 startfrom = lBinInLayer; endat = hBinInLayer;
275 }
else if (iIOLayer == 1) {
276 startfrom = lBinOutLayer; endat = hBinOutLayer;
281 Int_t ieaprime = iea;
285 if (iIOLayer == 0)htemp = (TH1D*)hILdEdxhitInEntaBin[ieaprime - 1]->Clone(Form(
"hL%d_Ea%d", iIOLayer, iea));
286 else if (iIOLayer == 1)htemp = (TH1D*)hOLdEdxhitInEntaBin[ieaprime - 1]->Clone(Form(
"hL%d_Ea%d", iIOLayer, iea));
289 double truncMean = 1.0;
290 if (htemp->Integral() < 250) truncMean = 1.0;
292 double binweights = 0.0;
294 for (
int ibin = startfrom; ibin <= endat; ibin++) {
296 if (htemp->GetBinContent(ibin) > 0) {
297 binweights += (htemp->GetBinContent(ibin) * htemp->GetBinCenter(ibin));
298 sumofbc += htemp->GetBinContent(ibin);
301 if (sumofbc > 0)truncMean = (double)(binweights / sumofbc);
302 else truncMean = 1.0;
305 if (truncMean <= 0)truncMean = 1.0;
306 onedcorIorOLtemp.push_back(truncMean);
309 ctmp->cd(((iea - 1) % 16) + 1);
310 htemp->SetFillColorAlpha(kGreen, 0.30);
311 htemp->SetTitle(Form(
"%s, #mean = %0.2f", htemp->GetTitle(), truncMean));
312 if (truncMean >= 1.02 || truncMean <= 0.98)htemp->SetFillColor(kYellow);
313 if (truncMean >= 1.05 || truncMean <= 0.95)htemp->SetFillColor(kOrange);
314 if (truncMean >= 1.10 || truncMean <= 0.90)htemp->SetFillColor(kRed);
315 htemp->DrawClone(
"hist");
316 tl->SetLineColor(kRed);
317 tl->SetX1(truncMean); tl->SetX2(truncMean);
318 tl->SetY1(0); tl->SetY2(htemp->GetMaximum());
319 tl->DrawClone(
"same");
321 tb->SetLineColor(kPink);
323 tb->SetFillColorAlpha(kPink, 0.15);
324 tb->SetX1(htemp->GetBinLowEdge(startfrom));
326 tb->SetX2(htemp->GetBinLowEdge(endat));
327 tb->SetY2(htemp->GetMaximum());
328 tb->DrawClone(
"same");
329 if (iea % 16 == 0)ctmp->Print(psname.str().c_str());
339 onedcorIorOL.push_back(onedcorIorOLtemp.at(ibinEA));
340 if (iIOLayer == 0)hILEntaConst->SetBinContent(iea + 1, onedcorIorOLtemp.at(ibinEA));
341 else if (iIOLayer == 1)hOLEntaConst->SetBinContent(iea + 1, onedcorIorOLtemp.at(ibinEA));
344 onedcors.push_back(onedcorIorOL);
345 onedcorIorOL.clear();
346 onedcorIorOLtemp.clear();
351 psname.str(
""); psname << Form(
"dedx_1dcell_%s.pdf]",
fSetPrefix.data());
352 ctmp->Print(psname.str().c_str());
355 TCanvas* cconst =
new TCanvas(
"FinalConstantHistoMap",
"Inner and Outer Histo", 800, 400);
356 cconst->Divide(2, 1);
357 cconst->cd(1); hILEntaConst->Draw();
358 cconst->cd(2); hOLEntaConst->Draw();
359 cconst->SaveAs(Form(
"FinalOneDConstantMap_%s.pdf",
fSetPrefix.data()));
362 B2INFO(
"dE/dx Calibration done for 1D cleanup correction");