9#include <cdc/calibration/CDCdEdx/CDCDedxBadWireAlgorithm.h>
11#include <cdc/geometry/CDCGeometryPar.h>
54 B2FATAL(
"There is no valid payload for BadWire and/or Wirgain");
60 vector<int>* wire = 0;
61 ttree->SetBranchAddress(
"wire", &wire);
63 vector<double>* hitvar = 0;
64 if (
m_isADC) ttree->SetBranchAddress(
"adccorr", &hitvar);
65 else ttree->SetBranchAddress(
"dedxhit", &hitvar);
71 hvarall.SetTitle(Form(
"dist %s; %s; %s",
m_suffix.data(),
m_varName.data(),
"entries"));
73 map<int, vector<double>> vhitvar;
75 for (
int i = 0; i < ttree->GetEntries(); ++i) {
77 for (
unsigned int ih = 0; ih < wire->size(); ++ih) {
78 double ivalue = hitvar->at(ih);
79 vhitvar[wire->at(ih)].push_back(ivalue);
98 map<int, vector<double>> qapars;
99 vector<double> vdefectwires, vbadwires, vdeadwires;
101 for (
unsigned int jw = 0; jw <
c_nwireCDC; ++jw) {
102 int ncount = 0, tcount = 0;
104 for (
unsigned int jh = 0; jh < vhitvar[jw].size(); ++jh) {
105 double jvalue = vhitvar[jw][jh];
112 bool badwire =
false;
114 qapars[0].push_back(0);
115 qapars[1].push_back(0);
116 qapars[2].push_back(0);
119 nmean = nmean / ncount;
123 for (
unsigned int kh = 0; kh < vhitvar[jw].size(); ++kh) {
124 double kvalue = vhitvar[jw][kh];
125 if (kvalue <
m_varMax) nrms += pow(kvalue - nmean, 2);
128 nrms =
sqrt(nrms / ncount);
131 double badfrac = 0.0;
132 if (tcount > 0) badfrac = (1.0 * tcount) / (tcount + ncount);
135 qapars[0].push_back(nmean);
136 qapars[1].push_back(nrms);
137 qapars[2].push_back(badfrac);
141 vdefectwires.push_back(0.0);
142 if (ncount == 0) vdeadwires.push_back(jw);
143 else vbadwires.push_back(jw);
144 }
else vdefectwires.push_back(1.0);
163 B2INFO(
"dE/dx Badwire Calibration done: " << vdefectwires.size() <<
" wires");
178 if (cruns == 0) B2INFO(
"CDCDedxBadWires: start exp " << expRun.first <<
" and run " << expRun.second <<
"");
183 int estart = erStart.first;
184 int rstart = erStart.second;
187 int rend = erEnd.second;
192 else m_suffix = Form(
"e%d_r%dr%d", estart, rstart, rend);
198 map<
int, vector<double>>& vhitvar)
201 TList* bdlist =
new TList();
202 bdlist->SetName(
"badwires");
204 TList* goodlist =
new TList();
205 goodlist->SetName(
"goodwires");
207 TList* hflist =
new TList();
208 hflist->SetName(
"highfracwires");
210 for (
unsigned int jw = 0; jw <
c_nwireCDC; ++jw) {
215 hvarhf->SetTitle(Form(
"%s, wire = %d; %s; entries",
m_suffix.data(), jw,
m_varName.data()));
217 int ncount = 0, tcount = 0;
219 for (
unsigned int jh = 0; jh < vhitvar[jw].size(); ++jh) {
220 double jvalue = vhitvar[jw][jh];
226 if (jvalue <
m_varMax * 10.) hvarhf->Fill(jvalue / 10.);
230 double badfrac = 0.0;
231 if (tcount > 0) badfrac = (1.0 * tcount) / (tcount + ncount);
232 hvar->SetTitle(Form(
"%s, wire = %d; %s; %0.01f",
m_suffix.data(), jw,
m_varName.data(), badfrac * 100));
235 if (count(inwires.begin(), inwires.end(), jw)) isbad =
true;
239 hvar->SetLineWidth(2);
240 hvar->SetLineColor(kRed);
247 if (hvar->Integral() > 100) goodlist->Add(hvar);
263 string listname = list->GetName();
264 string sfx = Form(
"%s_%s", listname.data(),
m_suffix.data());
266 TCanvas ctmp(Form(
"cdcdedx_%s", sfx.data()),
"", 1200, 1200);
268 ctmp.SetBatch(kTRUE);
271 psname << Form(
"cdcdedx_bdcal_%s.pdf[", sfx.data());
272 ctmp.Print(psname.str().c_str());
274 psname << Form(
"cdcdedx_bdcal_%s.pdf", sfx.data());
276 for (
int ih = 0; ih < list->GetSize(); ih++) {
278 TH1D* hist = (TH1D*)list->At(ih);
280 double frac = stod(hist->GetYaxis()->GetTitle());
282 TPaveText* pinfo =
new TPaveText(0.40, 0.63, 0.89, 0.89,
"NBNDC");
285 pinfo->AddText(Form(
"N: %0.00f", hist->Integral()));
286 pinfo->AddText(Form(
"hf: %0.00f%%(%0.00f%%)", frac,
m_fracThres * 100));
289 ctmp.cd(ih % 16 + 1);
290 hist->GetYaxis()->SetTitle(
"entries");
291 hist->SetFillColor(color);
296 if (listname ==
"badwires") {
297 TH1D* histhf = (TH1D*)hflist->At(ih);
298 if (hist->GetMaximum() < histhf->GetMaximum()) hist->SetMaximum(histhf->GetMaximum() * 1.05);
299 histhf->SetFillColor(kGray);
301 histhf->Draw(
"same");
304 if (((ih + 1) % 16 == 0) || ih == (list->GetSize() - 1)) {
305 ctmp.Print(psname.str().c_str());
311 psname << Form(
"cdcdedx_bdcal_%s.pdf]", sfx.data());
312 ctmp.Print(psname.str().c_str());
319 TCanvas cmap(Form(
"cmap_%s",
m_suffix.data()),
"", 800, 800);
320 cmap.SetTitle(
"CDC dE/dx bad wire status");
324 hxyAll->SetTitle(Form(
"wire status map (%s)",
m_suffix.data()));
332 hxyBad->Draw(
"same");
339 hxyDead->Draw(
"same");
342 int ndefect = nbad + ndead;
343 auto leg =
new TLegend(0.68, 0.80, 0.90, 0.92);
344 leg->SetBorderSize(0);
345 leg->SetLineWidth(3);
346 leg->SetHeader(Form(
"total defective: %d (~%0.02f%%)", ndefect, 100.*(ndefect) /
c_nwireCDC));
347 leg->AddEntry(hxyBad, Form(
"bad #rightarrow %d", nbad),
"p");
348 leg->AddEntry(hxyDead, Form(
"dead #rightarrow %d", ndead),
"p");
351 gStyle->SetLegendTextSize(0.025);
352 TPaveText* pt =
new TPaveText(-0.30, -1.47, -0.31, -1.30,
"br");
355 TText* text = pt->AddText(
"CDC-wire map: counter-clockwise and start from +x");
356 text->SetTextColor(kGray + 1);
359 cmap.SaveAs(Form(
"cdcdedx_bdcal_wiremap_%s.pdf",
m_suffix.data()));
370 B2INFO(
"Creating CDCGeometryPar object");
373 TH2F* temp =
new TH2F(Form(
"temp_%s_%s",
m_suffix.data(), suffix.data()),
"", 2400, -1.2, 1.2, 2400, -1.2, 1.2);
377 for (
unsigned int ilay = 0; ilay < c_maxNSenseLayers; ++ilay) {
378 for (
unsigned int iwire = 0; iwire < cdcgeo.
nWiresInLayer(ilay); ++iwire) {
380 double phi = 2.*TMath::Pi() * (float(iwire) / float(cdcgeo.
nWiresInLayer(ilay)));
381 double radius = cdcgeo.
senseWireR(ilay) / 100.;
382 double x = radius * cos(phi);
383 double y = radius * sin(phi);
384 if (suffix ==
"all") {
388 if (count(inwires.begin(), inwires.end(), jwire)) {
402 string qaname[3] = {
"mean",
"rms",
"high_fraction"};
407 for (
int iqa = 0; iqa < 3; iqa++) {
409 TH1D histqa(Form(
"%s_%s", qaname[iqa].data(),
m_suffix.data()),
"",
c_nwireCDC, -0.5, 14335.5);
411 for (
unsigned int jw = 0; jw <
c_nwireCDC; jw++) {
412 if (iqa == 2) histqa.SetBinContent(jw + 1, qapars[iqa][jw] * 100);
413 else histqa.SetBinContent(jw + 1, qapars[iqa][jw]);
416 TCanvas c_pars(Form(
"c_pars_%d", iqa),
"", 800, 600);
420 histqa.SetTitle(Form(
"%s vs wires (%s); wire ; %s", qaname[iqa].data(),
m_suffix.data(), qaname[iqa].data()));
424 TLine* lmin =
new TLine(-0.5, linemin[iqa], 14335.5, linemin[iqa]);
425 lmin->SetLineColor(kRed);
427 TLine* lmax =
new TLine(-0.5, linemax[iqa], 14335.5, linemax[iqa]);
428 lmax->SetLineColor(kRed);
431 c_pars.Print(Form(
"cdcdedx_bdcal_%s_%s.root", qaname[iqa].data(),
m_suffix.data()));
432 c_pars.Print(Form(
"cdcdedx_bdcal_%s_%s.pdf", qaname[iqa].data(),
m_suffix.data()));
443 TCanvas cstats(
"cstats",
"cstats", 1000, 500);
444 cstats.SetBatch(kTRUE);
450 hestats->SetName(Form(
"htstats_%s",
m_suffix.data()));
451 hestats->SetStats(0);
452 hestats->DrawCopy(
"");
458 htstats->SetName(Form(
"htstats_%s",
m_suffix.data()));
459 htstats->SetStats(0);
460 htstats->DrawCopy(
"");
463 cstats.Print(Form(
"cdcdedx_bdcal_qastats_%s.pdf",
m_suffix.data()));
void setHistCosmetics(TH2F *hist, Color_t color)
function to change histogram styles
void plotBadWireMap(const std::vector< double > &vbadwires, const std::vector< double > &vdeadwires)
function to plot wire status map (all, bad and dead)
double m_varMax
max range for input variable
double m_rmsThres
rms Threshold accepted for good wire
double m_varMin
min range for input variable
std::string m_varName
std::string to set var name (adc or dedx)
void getExpRunInfo()
function to get extract calibration run/exp
unsigned int c_nwireCDC
number of wires in CDC
double m_arms
average rms of dedx for all wires
bool m_isMakePlots
produce plots for status
DBObjPtr< CDCDedxBadWires > m_DBBadWires
Badwire DB object.
std::string m_suffix
suffix std::string for naming plots
void setTextCosmetics(TPaveText *pt, double size)
function to change text styles
CDCDedxBadWireAlgorithm()
Constructor: Sets the description, the properties and the parameters of the algorithm.
DBObjPtr< CDCGeometry > m_cdcGeo
Geometry of CDC.
double m_meanThres
mean Threshold accepted for good wire
DBObjPtr< CDCDedxWireGain > m_DBWireGains
Wiregain DB object.
double m_fracThres
high-frac Threshold accepted for good wire
int m_varBins
number of bins for input variable
void printCanvas(TList *list, TList *hflist, Color_t color)
function to print canvas
void plotQaPars(std::map< int, std::vector< double > > &qapars)
function to plot the QA (decision) parameters
void plotEventStats()
function to draw the stats
virtual EResult calibrate() override
cdcdedx badwire algorithm
void plotWireDist(const std::vector< double > &inwires, std::map< int, std::vector< double > > &vhitvar)
function to draw per wire plots
double m_amean
average mean of dedx for all wires
bool m_isADC
Use adc if(true) else dedx for calibration.
TH2F * getHistoPattern(const std::vector< double > &inwires, const std::string &suffix, int &total)
function to get wire map with input file (all, bad and dead)
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.
double senseWireR(int layerId) const
Returns radius of sense wire in each layer.
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...
double sqrt(double a)
sqrt for double
Abstract base class for different kinds of events.