9#include <cdc/calibration/CDCdEdx/CDCDedx1DCellAlgorithm.h>
49 setDescription(
"A calibration algorithm for the CDC dE/dx entrance angle cleanup correction");
61 B2FATAL(
"There is no valid previous payload for CDCDedx1DCell");
67 std::vector<double>* dedxhit = 0, *enta = 0;
68 std::vector<int>* layer = 0;
69 double pt = 0, costh = 0;
71 ttree->SetBranchAddress(
"dedxhit", &dedxhit);
72 ttree->SetBranchAddress(
"entaRS", &enta);
73 ttree->SetBranchAddress(
"layer", &layer);
74 ttree->SetBranchAddress(
"pt", &pt);
75 ttree->SetBranchAddress(
"costh", &costh);
90 std::vector<TH1D*> hdedxhit[2];
94 TH2D* hptcosth =
new TH2D(
"hptcosth",
"pt vs costh dist;pt;costh", 1000, -8.0, 8.0, 1000, -1.0, 1.0);
99 for (
int i = 0; i < ttree->GetEntries(); ++i) {
103 if (std::abs(costh) >
m_cosMax)
continue;
109 if (std::abs(pt) >
m_ptMax)
continue;
112 int rand = gRandom->Integer(100);
113 if (rand < 10) hptcosth->Fill(pt, costh);
115 for (
unsigned int j = 0; j < dedxhit->size(); ++j) {
117 if (dedxhit->at(j) == 0)
continue;
119 double entaval = enta->at(j);
122 if (ibin < 0 || ibin >
m_eaBin)
continue;
125 if (layer->at(j) < 8)mL = 0;
128 hdedxlay[mL]->Fill(dedxhit->at(j));
129 if (rand < 10) hentalay[mL]->Fill(entaval);
133 hdedxhit[mL][jbinea]->Fill(dedxhit->at(j));
137 for (
int il = 0; il < 2; il++) {
139 int minlay = 0, maxlay = 0;
143 hdedxlay[il]->SetTitle(Form(
"%s;%d;%d", hdedxlay[il]->GetTitle(), minlay, maxlay));
146 std::vector<double>tempconst;
156 TH1D* htemp = (TH1D*)hdedxhit[il][jea]->Clone(Form(
"h_%s_b%d_c",
m_label[il].data(), jea));
158 int minbin = 1, maxbin = 1;
169 tempconst.push_back(dedxmean);
171 hdedxhit[il][iea]->SetTitle(Form(
"%s, #mu_{truc} = %0.5f;%d;%d", hdedxhit[il][iea]->GetTitle(), dedxmean, minbin, maxbin));
175 std::vector<double>layerconst;
178 for (
int iea = 0; iea <
m_eaBin; iea++) {
181 layerconst.push_back(tempconst.at(jea));
213 for (
int il = 0; il < 2; il++) {
217 delete hdedxhit[il][iea];
222 for (
int il = 0; il < 2; il++) {
235 if (cruns == 0) B2INFO(
"CDCDedxBadWires: start exp " << expRun.first <<
" and run " << expRun.second <<
"");
240 int estart = erStart.first;
241 int rstart = erStart.second;
244 int eend = erEnd.first;
245 int rend = erEnd.second;
249 m_runExp = Form(
"Range (%d:%d,%d:%d)", estart, rstart, eend, rend);
251 else m_suffix = Form(
"e%d_r%dr%d", estart, rstart, rend);
258 std::map<int, std::vector<double>> bounds;
259 std::map<int, std::vector<int>> steps;
261 const std::array<int, 2> nDev{8, 4};
262 bounds[0] = {0, 108, 123, 133, 158, 183, 193, 208, 316};
263 steps[0] = {9, 3, 2, 1, 1, 2, 3, 9};
264 bounds[1] = {0, 38, 158, 278, 316};
265 steps[1] = {2, 1, 1, 2};
267 for (
int il = 0; il < 2; il++) {
269 for (
int ibin = 0; ibin <= nDev[il]; ibin++) bounds[il][ibin] = bounds[il][ibin] *
m_binSplit;
271 int ieaprime = -1, temp = -99, ibin = 0;
276 for (
int iea = 0; iea <
m_eaBin; iea++) {
279 if (iea %
int(bounds[il][ibin + 1]) == 0 && iea > 0) ibin++;
280 int diff = iea - int(bounds[il][ibin]);
281 if (diff % steps[il][ibin] == 0) ieaprime++;
282 }
else ieaprime = iea;
286 if (ieaprime != temp) {
289 double binvalue = pastbin + binwidth;
291 if (std::abs(binvalue) < 1e-5)binvalue = 0;
304 for (
int il = 0; il < 2; il++) {
306 std::string title = Form(
"dedxhit dist (%s): %s ; dedxhit;entries",
m_label[il].data(),
m_runExp.data());
308 hdedxlay[il]->SetTitle(Form(
"%s", title.data()));
313 if (
isVarBins) title = Form(
"entaRS dist (variable bins): %s: (%s); entaRS (#alpha);entries",
m_label[il].data(),
m_runExp.data());
314 else title = Form(
"entaRS dist (sym. bins): %s: (%s); entaRS (#alpha);entries",
m_label[il].data(),
m_runExp.data());
316 hentalay[il] =
new TH1D(Form(
"hentalay%s",
m_label[il].data()),
"",
m_eaBinLocal[il], nvarBins);
317 hentalay[il]->SetTitle(Form(
"%s", title.data()));
323 double width = max - min;
325 if (
isPrintLog) B2INFO(
"bin: " << iea <<
" ], min:" << min <<
" , max: " << max <<
" , width: " << width);
327 title = Form(
"%s: entaRS = (%0.03f to %0.03f)",
m_label[il].data(), min, max);
329 hdedxhit[il][iea]->SetTitle(Form(
"%s", title.data()));
339 double sum = hist->Integral();
340 if (sum <= 0 || hist->GetNbinsX() <= 0) {
341 binlow = 1; binhigh = 1;
345 binlow = 1.0; binhigh = 1.0;
346 double sumPer5 = 0.0, sumPer75 = 0.0;
347 for (
int ibin = 1; ibin <= hist->GetNbinsX(); ibin++) {
348 double bcdedx = hist->GetBinContent(ibin);
366 if (hist->Integral() < 100)
return 1.0;
368 if (binlow <= 0 || binhigh > hist->GetNbinsX())
return 1.0;
370 double binweights = 0., sumofbc = 0.;
371 for (
int ibin = binlow; ibin <= binhigh; ibin++) {
372 double bcdedx = hist->GetBinContent(ibin);
374 binweights += (bcdedx * hist->GetBinCenter(ibin));
378 if (sumofbc > 0)
return binweights / sumofbc;
386 B2INFO(
"dE/dx one cell calibration: Generating payloads");
388 for (
unsigned int il = 0; il < 2; il++) {
394 B2ERROR(
"merging failed because of unmatch bins (old " <<
m_eaBin <<
" new " << nbins <<
")");
396 for (
unsigned int iea = 0; iea < nbins; iea++) {
402 double binsize = TMath::Pi() /
m_onedcors[il].size();
404 auto computeAverages = [&](
double angLow,
double angHigh,
const std::string & label) {
405 unsigned int binLow = std::floor((angLow + TMath::Pi() / 2.0) / binsize);
406 unsigned int binHigh = std::floor((angHigh + TMath::Pi() / 2.0) / binsize);
407 double sum_new = 0.0, sum_prev = 0.0;
410 for (
unsigned int iea = binLow; iea < binHigh; ++iea) {
416 double avg_new = (count > 0) ? sum_new / count : 1.0;
417 double avg_prev = (count > 0) ? sum_prev / count : 1.0;
418 return std::make_pair(avg_new, avg_prev);
421 double negLow = -0.75, negHigh = -0.25;
422 double posLow = 0.25, posHigh = 0.75;
425 negLow = -0.5; negHigh = -0.2;
426 posLow = 0.2; posHigh = 0.5;
429 auto [avgNewNeg, avgPrevNeg] = computeAverages(negLow, negHigh,
"Negative");
430 auto [avgNewPos, avgPrevPos] = computeAverages(posLow, posHigh,
"Positive");
432 double avgNew = (avgNewNeg + avgNewPos) / 2.0;
433 double avgPrev = (avgPrevNeg + avgPrevPos) / 2.0;
434 double scaleFactor = avgPrev / avgNew;
436 for (
unsigned int iea = 0; iea <
m_eaBin; iea++) {
447 B2INFO(
"dE/dx Calibration done for CDCDedx1DCell");
454 std::map<
int, std::vector<int>> steps)
457 TCanvas cmfactor(
"cmfactor",
"Merging factors", 800, 400);
458 cmfactor.Divide(2, 1);
460 for (
int il = 0; il < 2; il++) {
462 nvarBins = &bounds[il][0];
464 TH1I* hist =
new TH1I(Form(
"hist_%s",
m_label[il].data()),
"", nDev[il], nvarBins);
465 hist->SetTitle(Form(
"Merging factor for %s bins;binindex;merge-factors",
m_label[il].data()));
467 for (
int ibin = 0; ibin < nDev[il]; ibin++) hist->SetBinContent(ibin + 1, steps[il][ibin]);
470 hist->SetFillColor(kYellow);
475 cmfactor.SaveAs(Form(
"cdcdedx_1dcell_mergefactor%s.pdf",
m_suffix.data()));
476 cmfactor.SaveAs(Form(
"cdcdedx_1dcell_mergefactor%s.root",
m_suffix.data()));
483 TCanvas ctmp(
"tmp",
"tmp", 1200, 1200);
485 std::stringstream psname;
487 psname << Form(
"cdcdedx_1dcell_dedxhit%s.pdf[",
m_suffix.data());
488 ctmp.Print(psname.str().c_str());
490 psname << Form(
"cdcdedx_1dcell_dedxhit%s.pdf",
m_suffix.data());
492 for (
int il = 0; il < 2; il++) {
496 int minbin = std::stoi(hdedxhit[il][jea]->GetXaxis()->GetTitle());
497 int maxbin = std::stoi(hdedxhit[il][jea]->GetYaxis()->GetTitle());
499 ctmp.cd(jea % 16 + 1);
500 hdedxhit[il][jea]->SetFillColor(4 + il);
502 hdedxhit[il][jea]->SetTitle(Form(
"%s;dedxhit;entries", hdedxhit[il][jea]->GetTitle()));
503 hdedxhit[il][jea]->DrawClone(
"hist");
504 TH1D* htempC = (TH1D*)hdedxhit[il][jea]->Clone(Form(
"%sc2", hdedxhit[il][jea]->GetName()));
505 htempC->GetXaxis()->SetRange(minbin, maxbin);
506 htempC->SetFillColor(kGray);
507 htempC->DrawClone(
"same hist");
510 ctmp.Print(psname.str().c_str());
518 psname << Form(
"cdcdedx_1dcell_dedxhit%s.pdf]",
m_suffix.data());
519 ctmp.Print(psname.str().c_str());
526 TCanvas cdedxlayer(
"layerdedxhit",
"Inner and Outer Layer dedxhit dist", 900, 400);
527 cdedxlayer.Divide(2, 1);
529 for (
int il = 0; il < 2; il++) {
530 int minlay = 0, maxlay = 0;
532 minlay = std::stoi(hdedxlay[il]->GetXaxis()->GetTitle());
533 maxlay = std::stoi(hdedxlay[il]->GetYaxis()->GetTitle());
534 double lowedge = hdedxlay[il]->GetXaxis()->GetBinLowEdge(minlay);
535 double upedge = hdedxlay[il]->GetXaxis()->GetBinUpEdge(maxlay);
536 hdedxlay[il]->SetTitle(Form(
"%s, trunc #rightarrow: %0.02f - %0.02f;dedxhit;entries", hdedxlay[il]->GetTitle(), lowedge, upedge));
539 cdedxlayer.cd(il + 1);
540 hdedxlay[il]->SetFillColor(kYellow);
541 hdedxlay[il]->Draw(
"histo");
544 TH1D* hdedxlayC = (TH1D*)hdedxlay[il]->Clone(Form(
"hdedxlayC%d", il));
545 hdedxlayC->GetXaxis()->SetRange(minlay, maxlay);
546 hdedxlayC->SetFillColor(kAzure + 1);
547 hdedxlayC->Draw(
"same histo");
551 cdedxlayer.SaveAs(Form(
"cdcdedx_1dcell_dedxlay%s.pdf",
m_suffix.data()));
552 cdedxlayer.SaveAs(Form(
"cdcdedx_1dcell_dedxlay%s.root",
m_suffix.data()));
559 TCanvas ceadist(
"ceadist",
"Enta distributions", 800, 400);
560 ceadist.Divide(2, 1);
562 for (
int il = 0; il < 2; il++) {
566 hentalay[il]->SetFillColor(kYellow);
567 hentalay[il]->Draw(
"hist");
570 TCanvas cptcos(
"cptcos",
"pt vs costh dist.", 400, 400);
572 hptcosth->Draw(
"colz");
574 cptcos.SaveAs(Form(
"cdcdedx_ptcosth_%s.pdf",
m_suffix.data()));
575 ceadist.SaveAs(Form(
"cdcdedx_1dcell_enta%s.pdf",
m_suffix.data()));
576 ceadist.SaveAs(Form(
"cdcdedx_1dcell_enta%s.root",
m_suffix.data()));
583 TH1D* hconst, *hconstvar;
589 std::string title = Form(
"calibration const dist: %s: (%s); entaRS (#alpha); entries",
m_label[il].data(),
m_runExp.data());
590 hconst->SetTitle(Form(
"%s", title.data()));
592 hconstvar =
new TH1D(Form(
"hconstvar%s",
m_label[il].data()),
"",
m_eaBinLocal[il], nvarBins);
593 title = Form(
"calibration const dist (var bins): %s: (%s); entaRS (#alpha);entries",
m_label[il].data(),
m_runExp.data());
594 hconstvar->SetTitle(Form(
"%s", title.data()));
598 hconstvar->SetBinContent(iea + 1, tempconst.at(iea));
601 for (
int jea = 0; jea <
m_eaBin; jea++) hconst->SetBinContent(jea + 1, layerconst.at(jea));
603 gStyle->SetOptStat(
"ne");
604 TCanvas cconst(
"cconst",
"Calirbation Constants", 800, 400);
607 cconst.SetWindowSize(1000, 800);
611 hconst->SetFillColor(kYellow);
612 hconst->Draw(
"histo");
615 hconstvar->SetFillColor(kBlue);
616 hconstvar->Draw(
"hist");
618 cconst.SaveAs(Form(
"cdcdedx_1dcell_relconst%s_%s.pdf",
m_label[il].data(),
m_suffix.data()));
619 cconst.SaveAs(Form(
"cdcdedx_1dcell_relconst%s_%s.root",
m_label[il].data(),
m_suffix.data()));
630 TH1D* hnewconst[2], *holdconst[2];
631 double min[2], max[2];
633 for (
unsigned int il = 0; il < 2; il++) {
636 std::string title = Form(
"final calibration const dist (%s): %s; entaRS (#alpha); entries",
m_label[il].data(),
m_runExp.data());
638 hnewconst[il]->SetTitle(Form(
"%s", title.data()));
640 title = Form(
"old calibration const dist (%s): %s; entaRS (#alpha); entries",
m_label[il].data(),
m_runExp.data());
642 holdconst[il]->SetTitle(Form(
"%s", title.data()));
644 for (
unsigned int iea = 0; iea < nbins; iea++) {
646 holdconst[il]->SetBinContent(iea + 1, prev);
647 hnewconst[il]->SetBinContent(iea + 1,
m_onedcors[il][iea]);
649 min[il] = hnewconst[il]->GetMinimum();
650 max[il] = hnewconst[il]->GetMaximum();
654 if (max[1] < max[0])max[1] = max[0];
655 if (min[1] > min[0])min[1] = min[0];
657 gStyle->SetOptStat(
"ne");
658 TCanvas cfconst(
"cfconst",
"Final calirbation constants", 800, 400);
659 cfconst.Divide(2, 1);
661 for (
int il = 0; il < 2; il++) {
663 hnewconst[il]->GetYaxis()->SetRangeUser(min[1] * 0.95, max[1] * 1.05);
664 hnewconst[il]->SetLineColor(kBlack);
665 hnewconst[il]->Draw(
"histo");
666 holdconst[il]->SetLineColor(kRed);
667 holdconst[il]->Draw(
"histo same");
669 auto legend =
new TLegend(0.4, 0.75, 0.56, 0.85);
670 legend->AddEntry(holdconst[il],
"Old",
"lep");
671 legend->AddEntry(hnewconst[il],
"New",
"lep");
675 cfconst.SaveAs(Form(
"cdcdedx_1dcell_fconsts%s.pdf",
m_suffix.data()));
676 cfconst.SaveAs(Form(
"cdcdedx_1dcell_fconsts%s.root",
m_suffix.data()));
678 for (
int il = 0; il < 2; il++) {
679 delete hnewconst[il];
680 delete holdconst[il];
688 TCanvas cstats(
"cstats",
"cstats", 1000, 500);
689 cstats.SetBatch(kTRUE);
695 hestats->SetName(Form(
"hestats_%s",
m_suffix.data()));
696 hestats->SetStats(0);
697 hestats->DrawCopy(
"");
703 hestats->SetName(Form(
"htstats_%s",
m_suffix.data()));
704 htstats->DrawCopy(
"");
705 hestats->SetStats(0);
707 cstats.Print(Form(
"cdcdedx_1dcell_stats_%s.pdf",
m_suffix.data()));
CDCDedx1DCellAlgorithm()
Constructor: Sets the description, the properties and the parameters of the algorithm.
std::string m_label[2]
add inner/outer layer label
double m_eaMax
upper edge of entrance angle
void plotMergeFactor(std::map< int, std::vector< double > > bounds, const std::array< int, 2 > nDev, std::map< int, std::vector< int > > steps)
function to plot merging factor
double m_truncMax
upper threshold on truncation
int m_binSplit
multiply nbins by this factor in full range
std::array< std::vector< int >, 2 > m_binIndex
symm/Var bin numbers
double m_truncMin
lower threshold on truncation
double m_adjustFac
factor with that one what to adjust baseline
void getTruncatedBins(TH1D *hist, int &binlow, int &binhigh)
function to get bins of truncation from histogram
void CreateBinMapping()
class function to create vectors for bin mapping (Var->symm)
double m_chargeType
charge type for baseline adj
void getExpRunInfo()
function to get extract calibration run/exp
DBObjPtr< CDCDedx1DCell > m_DBOneDCell
One cell correction DB object.
double m_cosMax
a limit on cos theta
bool isPrintLog
print more debug information
std::array< std::vector< double >, 2 > m_binValue
enta Var bin values
std::string m_suffix
add suffix to all plot name
int m_eaB
reset # of bins for entrance angle for each experiment
double getTruncationMean(TH1D *hist, int binlow, int binhigh)
function to get truncated mean
double m_ptMax
a limit on transverse momentum
void plotConstants()
function to draw the old/new final constants
std::vector< int > m_eaBinLocal
bool isFixTrunc
true = fix window for all out/inner layers
bool isVarBins
true: if variable bin size is requested
void plotdedxHist(std::vector< TH1D * > hdedxhit[2])
function to draw the dE/dx histogram in enta bins
void defineHisto(std::vector< TH1D * > hdedxhit[2], TH1D *hdedxlay[2], TH1D *hentalay[2])
function to define histograms
double m_eaBW
binwdith of entrance angle bin
bool isRotSymm
if rotation symmetry requested
std::string m_runExp
add suffix to all plot name
void plotEventStats()
function to draw the stats plots
int rotationalBin(int nbin, int ibin)
class function to set rotation symmetry
virtual EResult calibrate() override
1D cell algorithm
void plotQaPars(TH1D *hentalay[2], TH2D *hptcosth)
function to draw pt vs costh and entrance angle distribution for Inner/Outer layer
double m_dedxMax
upper edge of dedxhit
void createPayload()
function to generate final constants
bool isMakePlots
produce plots for status
void plotRelConst(std::vector< double >tempconst, std::vector< double >layerconst, int il)
function to draw symm/Var layer constant
std::vector< std::vector< double > > m_onedcors
final vectors of calibration
bool isMerge
print more debug information
double m_dedxMin
lower edge of dedxhit
void plotLayerDist(TH1D *hdedxL[2])
function to draw dedx dist.
double m_eaMin
lower edge of entrance angle
dE/dx wire gain calibration constants
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
void updateDBObjPtrs(const unsigned int event, const int run, const int experiment)
Updates any DBObjPtrs by calling update(event) for DBStore.
void setDescription(const std::string &description)
Set algorithm description (in constructor)
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
EResult
The result of calibration.
@ c_OK
Finished successfully =0 in Python.
@ c_NotEnoughData
Needs more data =2 in Python.
CalibrationAlgorithm(const std::string &collectorModuleName)
Constructor - sets the prefix for collected objects (won't be accesses until execute(....
std::shared_ptr< T > getObjectPtr(const std::string &name, const std::vector< Calibration::ExpRun > &requestedRuns)
Get calibration data object by name and list of runs, the Merge function will be called to generate t...
Abstract base class for different kinds of events.