9#include <cdc/calibration/CDCdEdx/CDCDedxWireGainAlgorithm.h>
11#include <cdc/geometry/CDCGeometryPar.h>
37 setDescription(
"A calibration algorithm for CDC dE/dx wire gains");
50 B2FATAL(
"There is no valid payload for BadWires and/or Wirgain");
56 vector<int>* wire = 0;
57 ttree->SetBranchAddress(
"wire", &wire);
59 vector<double>* dedxhit = 0;
60 ttree->SetBranchAddress(
"dedxhit", &dedxhit);
63 map<int, vector<double>> wirededx;
66 array<TH1D*, 2> hdedxL;
67 string label[2] = {
"IL",
"OL"};
68 for (
int il = 0; il < 2; il++)
71 for (
int i = 0; i < ttree->GetEntries(); ++i) {
73 for (
unsigned int j = 0; j < wire->size(); ++j) {
74 int jwire = wire->at(j);
75 double jhitdedx = dedxhit->at(j);
76 wirededx[jwire].push_back(jhitdedx);
78 if (jwire < 1280) hdedxL[0]->Fill(jhitdedx);
79 else hdedxL[1]->Fill(jhitdedx);
85 for (
unsigned int jw = 0; jw < c_nwireCDC; ++jw)
86 if (wirededx[jw].size() <= 100) minstat++;
91 array<unsigned int, 2> minbinL, maxbinL;
92 for (
int il = 0; il < 2; il++) {
94 hdedxL[il]->SetTitle(Form(
"%s(%s);%d;%d", label[il].data(),
m_suffix.data(), minbinL[il], maxbinL[il]));
97 vector<double> vrel_mean;
98 vector<double> vdedx_means;
100 vector<TH1D*> hdedxhit(c_nwireCDC);
102 B2INFO(
"Creating CDCGeometryPar object");
105 vector<double> layermean(c_maxNSenseLayers);
107 int activelayers = 0;
108 double layeravg = 0.0;
111 for (
unsigned int il = 0; il < c_maxNSenseLayers; ++il) {
116 for (
unsigned int iw = 0; iw < cdcgeo.
nWiresInLayer(il); ++iw) {
121 for (
unsigned int ih = 0; ih < wirededx[jwire].size(); ++ih) {
122 hdedxhit[jwire]->Fill(wirededx[jwire][ih]);
125 unsigned int minbin, maxbin;
137 hdedxhit[jwire]->SetTitle(Form(
"dedxhit-dist, wire: %d (%s);%d;%d", jwire,
m_suffix.data(), minbin, maxbin));
140 if (
m_DBBadWires->getBadWireStatus(jwire) == kTRUE) dedxmean = 0.0;
143 vrel_mean.push_back(dedxmean);
147 vdedx_means.push_back(dedxmean * prewg);
148 B2INFO(
"merged-wireGain: [" << jwire <<
"], prewgvious = " << prewg <<
", rel = " << dedxmean <<
", merged = " << vdedx_means.at(
150 }
else vdedx_means.push_back(dedxmean);
153 if (vdedx_means.at(jwire) > 0) {
154 layermean[il] += vdedx_means.at(jwire);
159 if (activewires > 0) layermean[il] /= activewires;
160 else layermean[il] = 1.0;
163 if (il >= 8 && layermean[il] > 0) {
164 layeravg += layermean[il];
171 if (activelayers > 0) layeravg /= activelayers;
173 for (
unsigned int iw = 0; iw < c_nwireCDC; ++iw) {
174 vrel_mean.at(iw) /= layeravg;
175 vdedx_means.at(iw) /= layeravg;
210 if (cruns == 0) B2INFO(
"CDCDedxWireGain: start exp " << expRun.first <<
" and run " << expRun.second <<
"");
215 int estart = erStart.first;
216 int rstart = erStart.second;
219 int rend = erEnd.second;
224 else m_suffix = Form(
"e%d_r%dr%d", estart, rstart, rend);
231 B2INFO(
"dE/dx Calibration done for " << vdedx_means.size() <<
" CDC wires");
241 if (hdedxhit->Integral() < 100)
return 1.0;
243 if (binlow <= 0 || binhigh > hdedxhit->GetNbinsX())
return 1.0;
245 double binweights = 0., sumofbc = 0.;
246 for (
int ibin = binlow; ibin <= binhigh; ibin++) {
247 double bcdedx = hdedxhit->GetBinContent(ibin);
249 binweights += (bcdedx * hdedxhit->GetBinCenter(ibin));
253 if (sumofbc > 0)
return binweights / sumofbc;
262 double sum = hdedxhit->Integral();
263 if (sum <= 0 || hdedxhit->GetNbinsX() <= 0) {
264 binlow = 1; binhigh = 1;
268 binlow = 1.0; binhigh = 1.0;
269 double sumPer5 = 0.0, sumPer75 = 0.0;
270 for (
int ibin = 1; ibin <= hdedxhit->GetNbinsX(); ibin++) {
271 double bcdedx = hdedxhit->GetBinContent(ibin);
288 TCanvas cldedx(
"cldedx",
"IL/OL dedxhit dist", 900, 400);
291 for (
int il = 0; il < 2; il++) {
294 int minbin = stoi(hdedxL[il]->GetXaxis()->GetTitle());
295 int maxbin = stoi(hdedxL[il]->GetYaxis()->GetTitle());
296 double lowedge = hdedxL[il]->GetXaxis()->GetBinLowEdge(minbin);
297 double upedge = hdedxL[il]->GetXaxis()->GetBinUpEdge(maxbin);
299 hdedxL[il]->SetFillColor(kYellow);
300 hdedxL[il]->SetTitle(Form(
"%s, trunc(%0.02f - %0.02f);dedxhit;entries", hdedxL[il]->GetTitle(), lowedge, upedge));
301 hdedxL[il]->Draw(
"histo");
303 TH1D* hdedxLC = (TH1D*)hdedxL[il]->Clone(Form(
"%s_c", hdedxL[il]->GetName()));
304 hdedxLC->GetXaxis()->SetRange(minbin, maxbin);
305 hdedxLC->SetFillColor(kAzure + 1);
306 hdedxLC->Draw(
"same histo");
309 cldedx.SaveAs(Form(
"cdcdedx_wgcal_layerdedx_%s.pdf",
m_suffix.data()));
316 TCanvas ctmp(Form(
"cdcdedx_%s",
m_suffix.data()),
"", 1200, 1200);
318 ctmp.SetBatch(kTRUE);
321 psname << Form(
"cdcdedx_wgcal_%s.pdf[",
m_suffix.data());
322 ctmp.Print(psname.str().c_str());
324 psname << Form(
"cdcdedx_wgcal_%s.pdf",
m_suffix.data());
326 for (
unsigned int iw = 0; iw < hist.size(); iw++) {
328 int minbin = stoi(hist[iw]->GetXaxis()->GetTitle());
329 int maxbin = stoi(hist[iw]->GetYaxis()->GetTitle());
331 hist[iw]->SetFillColor(kYellow - 9);
332 hist[iw]->SetTitle(Form(
"%s, rel. #mu_{trunc} %0.03f;dedxhit;entries", hist[iw]->GetTitle(), vrel_mean.at(iw)));
335 hist[iw]->SetLineColor(kRed);
336 hist[iw]->SetLineWidth(2);
338 ctmp.cd(iw % 16 + 1);
339 hist[iw]->DrawCopy();
341 TH1D* hdedxhitC = (TH1D*)hist[iw]->Clone(Form(
"%sC", hist[iw]->GetName()));
342 hdedxhitC->GetXaxis()->SetRange(minbin, maxbin);
343 hdedxhitC->SetFillColor(kAzure + 1);
344 hdedxhitC->DrawCopy(
"same histo");
346 if (((iw + 1) % 16 == 0) || iw == (hist.size() - 1)) {
347 ctmp.Print(psname.str().c_str());
355 psname << Form(
"cdcdedx_wgcal_%s.pdf]",
m_suffix.data());
356 ctmp.Print(psname.str().c_str());
364 TCanvas cwconst(
"cwconst",
"", 900, 500);
365 TCanvas cwconstvar(
"cwconstvar",
"", 500, 400);
367 array<TH1D*, 2> hconstpw, hconstpwvar;
369 for (
int i = 0; i < 2; i++) {
371 hconstpw[i] =
new TH1D(Form(
"hconstpw_%d_%s", i,
m_suffix.data()),
"", c_nwireCDC, -0.5, 14335.5);
372 hconstpw[i]->SetTitle(Form(
"wiregain const (%s); wire numbers;<dedxhit>",
m_suffix.data()));
374 && i == 0) hconstpw[i]->SetTitle(Form(
"merged wiregain rel-const (%s), avg = %0.03f; wire numbers;<dedxhit>",
m_suffix.data(),
377 hconstpwvar[i] =
new TH1D(Form(
"hconstpwvar_%s",
m_suffix.data()),
"", 400, -0.5, 2.5);
378 hconstpwvar[i]->SetTitle(Form(
"wiregain const (%s); wire gains; nentries",
m_suffix.data()));
380 && i == 0) hconstpwvar[i]->SetTitle(Form(
"merged wiregain rel-const (%s), avg = %0.03f; wire gains; nentries",
m_suffix.data(),
383 for (
unsigned int iw = 0; iw < c_nwireCDC; iw++) {
385 double gain = vdedx_means.at(iw);
386 if (
m_isMerge && i == 1) gain = vrel_mean.at(iw);
387 hconstpw[i]->SetBinContent(iw + 1, gain);
388 hconstpwvar[i]->Fill(gain);
390 if (iw % 500 == 0) hconstpw[i]->GetXaxis()->SetBinLabel(iw + 1, Form(
"w%d", iw + 1));
393 hconstpw[i]->SetLineColor(i * 2 + 2);
394 hconstpw[i]->LabelsOption(
"u",
"X");
395 hconstpw[i]->GetYaxis()->SetRangeUser(-0.1, 3.5);
396 hconstpw[i]->LabelsDeflate();
398 hconstpwvar[i]->SetFillColor(i * 2 + 2);
399 hconstpwvar[i]->Scale(1 / hconstpwvar[i]->GetMaximum());
404 hconstpw[0]->Draw(
"");
405 if (
m_isMerge) hconstpw[1]->Draw(
"same");
408 hconstpwvar[0]->Draw(
"hist");
409 if (
m_isMerge) hconstpwvar[1]->Draw(
"hist same");
411 cwconst.SaveAs(Form(
"cdcdedx_wgcal_wireconst_%s.pdf",
m_suffix.data()));
412 cwconstvar.SaveAs(Form(
"cdcdedx_wgcal_wireconstvar_%s.pdf",
m_suffix.data()));
419 TH1D hlayeravg(Form(
"hlayeravg_%s",
m_suffix.data()),
"", layermean.size(), -0.5, 55.5);
420 hlayeravg.SetTitle(Form(
"layer gain avg (%s); layer numbers;<dedxhit>",
m_suffix.data()));
422 for (
unsigned int il = 0; il < layermean.size(); il++) {
423 hlayeravg.SetBinContent(il + 1, layermean[il]);
424 if (il % 2 == 0 || il == layermean.size() - 1) hlayeravg.GetXaxis()->SetBinLabel(il + 1, Form(
"L%d", il));
427 TCanvas clayeravg(
"clayeravg",
"clayeravg", 800, 500);
428 clayeravg.SetGridy(1);
430 gStyle->SetOptStat(
"ne");
431 hlayeravg.LabelsOption(
"u",
"X");
432 hlayeravg.SetLineColor(kBlue);
433 hlayeravg.GetYaxis()->SetRangeUser(-0.1, 3.5);
434 hlayeravg.SetTitle(Form(
"%s, avg = %0.03f (abs)", hlayeravg.GetTitle(), layeravg));
435 hlayeravg.LabelsDeflate();
437 TLine* tl =
new TLine(-0.5, layeravg, 55.5, layeravg);
438 tl->SetLineColor(kRed);
439 tl->DrawClone(
"same");
440 clayeravg.SaveAs(Form(
"cdcdedx_wgcal_layeravg_%s.pdf",
m_suffix.data()));
450 TCanvas clconst(
"clconst",
"", 800, 500);
451 clconst.Divide(2, 2);
452 clconst.SetBatch(kTRUE);
454 stringstream psnameL;
455 psnameL << Form(
"cdcdedx_wgcal_layerconst_%s.pdf[",
m_suffix.data());
456 clconst.Print(psnameL.str().c_str());
457 psnameL.str(
""); psnameL << Form(
"cdcdedx_wgcal_layerconst_%s.pdf",
m_suffix.data());
461 for (
unsigned int il = 0; il < layermean.size(); ++il) {
464 TH1D hconstpl(Form(
"hconstpwvar_l%d_%s", il,
m_suffix.data()),
"", nwires, jwire, jwire + nwires);
465 hconstpl.SetTitle(Form(
"abs-const, layer: %d (%s); wire numbers;<dedxhit>", il,
m_suffix.data()));
467 for (
unsigned int iw = 0; iw < nwires; ++iw) {
468 hconstpl.SetBinContent(iw + 1, vdedx_means.at(jwire));
470 if (iw % 10 == 0) hconstpl.GetXaxis()->SetBinLabel(iw + 1, Form(
"w%d", jwire));
472 if (iw % 15 == 0) hconstpl.GetXaxis()->SetBinLabel(iw + 1, Form(
"w%d", jwire));
477 double lmean = layermean.at(il) / layeravg;
479 clconst.cd(il % 4 + 1);
480 gStyle->SetOptStat(
"ne");
482 hconstpl.SetTitle(Form(
"%s, avg = %0.03f", hconstpl.GetTitle(), lmean));
484 if (il < 8) hconstpl.GetYaxis()->SetRangeUser(-0.1, 4.0);
485 else hconstpl.GetYaxis()->SetRangeUser(-0.1, 2.0);
486 hconstpl.SetFillColor(kAzure - 1);
487 hconstpl.LabelsOption(
"u",
"X");
488 hconstpl.DrawCopy(
"hist");
490 TLine* tlc =
new TLine(jwire - nwires, lmean, jwire, lmean);
491 tlc->SetLineColor(kRed);
492 tlc->DrawClone(
"same");
494 if ((il + 1) % 4 == 0 || (il + 1) == layermean.size()) {
495 clconst.Print(psnameL.str().c_str());
503 psnameL << Form(
"cdcdedx_wgcal_layerconst_%s.pdf]",
m_suffix.data());
504 clconst.Print(psnameL.str().c_str());
511 TCanvas cstats(
"cstats",
"cstats", 800, 400);
512 cstats.SetBatch(kTRUE);
518 hestats->SetName(Form(
"hestats_%s",
m_suffix.data()));
519 hestats->SetStats(0);
520 hestats->DrawCopy(
"");
526 htstats->SetName(Form(
"htstats_%s",
m_suffix.data()));
527 htstats->SetStats(0);
528 htstats->DrawCopy(
"");
531 cstats.Print(Form(
"cdcdedx_wgcal_stats_%s.pdf",
m_suffix.data()));
void plotWGPerLayer(const std::vector< double > &vdedx_means, const std::vector< double > &layermean, double layeravg)
function to draw WG per layer
void plotLayerGain(const std::vector< double > &layermean, double layeravg)
function to draw layer gains
void plotLayerDist(std::array< TH1D *, 2 > hdedxL)
function to draw dE/dx for inner/outer layer
bool m_isMerge
merge payload at the time of calibration
double m_truncMax
max trunc range for mean
double m_truncMin
min trunc range for mean
void getTruncatedBins(TH1D *hdedxhit, unsigned int &binlow, unsigned int &binhigh)
function to get bins of truncation from histogram
void plotWireGain(const std::vector< double > &vdedx_means, const std::vector< double > &vrel_mean, double layeravg)
function to draw wire gains
void getExpRunInfo()
function to get extract calibration run/exp
void createPayload(const std::vector< double > &vdedx_tmeans)
function to finally store new payload after full calibration
void plotWireDist(const std::vector< TH1D * > &hist, const std::vector< double > &vrel_mean)
function to draw dE/dx histograms for each wire
bool m_isMakePlots
Save arithmetic and truncated mean for the 'dedx' values.
DBObjPtr< CDCDedxBadWires > m_DBBadWires
Bad wire DB object.
std::string m_suffix
suffix string to separate plots
double getTruncationMean(TH1D *hdedxhit, int binlow, int binhigh)
function to get mean of truncation from histogram
int m_dedxBins
number of bins for dedx histogram
DBObjPtr< CDCGeometry > m_cdcGeo
Geometry of CDC.
CDCDedxWireGainAlgorithm()
Constructor: Sets the description, the properties and the parameters of the algorithm.
DBObjPtr< CDCDedxWireGain > m_DBWireGains
Wire gain DB object.
void plotEventStats()
function to draw statstics
virtual EResult calibrate() override
Wire gain algorithm.
double m_dedxMax
max dedx range for wiregain cal
double m_dedxMin
min dedx range for wiregain cal
bool m_isWireTruc
method of trunc range for mean
dE/dx wire gain calibration constants
The Class for CDC Geometry Parameters.
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
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.