Belle II Software development
CDCDedxBadWireAlgorithm.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <cdc/calibration/CDCdEdx/CDCDedxBadWireAlgorithm.h>
10
11#include <cdc/geometry/CDCGeometryPar.h>
12
13#include <TCanvas.h>
14#include <TLegend.h>
15#include <TLine.h>
16#include <TMath.h>
17#include <TStyle.h>
18
19using namespace Belle2;
20using namespace CDC;
21using namespace std;
22
23//-----------------------------------------------------------------
24// Implementation
25//-----------------------------------------------------------------
27 CalibrationAlgorithm("CDCDedxElectronCollector"),
28 c_nwireCDC(c_nSenseWires),
29 m_isMakePlots(true),
30 m_isADC(false),
31 m_varBins(100),
32 m_varMin(0.0),
33 m_varMax(7.0),
34 m_meanThres(1.0),
35 m_rmsThres(1.0),
36 m_fracThres(1.0),
37 m_varName("hitdedx"),
38 m_suffix("")
39{
40 // Set module properties
41 setDescription("A calibration algorithm for CDC dE/dx bad wires");
42}
43
44//-----------------------------------------------------------------
45// Run the calibration
46//-----------------------------------------------------------------
48{
49
51
52 //old wg for book-keeping previous bad+dead
53 if (!m_DBBadWires.isValid() || !m_DBWireGains.isValid())
54 B2FATAL("There is no valid payload for BadWire and/or Wirgain");
55
56 // Get data objects
57 auto ttree = getObjectPtr<TTree>("tree");
58 if (ttree->GetEntries() < 1000) return c_NotEnoughData;
59
60 vector<int>* wire = 0;
61 ttree->SetBranchAddress("wire", &wire);
62
63 vector<double>* hitvar = 0;
64 if (m_isADC) ttree->SetBranchAddress("adccorr", &hitvar);
65 else ttree->SetBranchAddress("dedxhit", &hitvar);
66
67 if (m_isADC) m_varName = "hitadc";
68 m_suffix = Form("%s_%s", m_varName.data(), m_suffix.data());
69
70 TH1D hvarall(Form("hvarall_%s", m_suffix.data()), "", m_varBins, m_varMin, m_varMax);
71 hvarall.SetTitle(Form("dist %s; %s; %s", m_suffix.data(), m_varName.data(), "entries"));
72
73 map<int, vector<double>> vhitvar;
74
75 for (int i = 0; i < ttree->GetEntries(); ++i) {
76 ttree->GetEvent(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);
80 hvarall.Fill(ivalue);
81 }
82 }
83
84 m_amean = hvarall.GetMean();
85 m_arms = hvarall.GetRMS();
86
87 // Commenting minstat cut on March 2024,
88 //it is not worth to skip the calibration if there is more than 5% bad/dead wires
89
90 //return if >5% bad wire or null histogram
91 // int minstat = 0;
92 // for (unsigned int jw = 0; jw < c_nwireCDC; ++jw)
93 // if (vhitvar[jw].size() <= 100) minstat++;
94 // if (minstat > 0.05 * c_nwireCDC) return c_NotEnoughData;
95
96 if (m_amean == 0 || m_arms == 0) return c_NotEnoughData;
97
98 map<int, vector<double>> qapars;
99 vector<double> vdefectwires, vbadwires, vdeadwires;
100
101 for (unsigned int jw = 0; jw < c_nwireCDC; ++jw) {
102 int ncount = 0, tcount = 0;
103 double nmean = 0.;
104 for (unsigned int jh = 0; jh < vhitvar[jw].size(); ++jh) {
105 double jvalue = vhitvar[jw][jh];
106 if (jvalue < m_varMax) {
107 ncount++;
108 nmean += jvalue;
109 } else tcount++;
110 }
111
112 bool badwire = false;
113 if (ncount < 100) {
114 qapars[0].push_back(0);
115 qapars[1].push_back(0);
116 qapars[2].push_back(0);
117 badwire = true; //partial dead
118 } else {
119 nmean = nmean / ncount;
120 if (abs(nmean - m_amean) / m_amean > m_meanThres) badwire = true;
121
122 double nrms = 0.;
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);
126 }
127
128 nrms = sqrt(nrms / ncount);
129 if (abs(nrms - m_arms) / m_arms > m_rmsThres) badwire = true;
130
131 double badfrac = 0.0;
132 if (tcount > 0) badfrac = (1.0 * tcount) / (tcount + ncount);
133 if (badfrac > m_fracThres) badwire = true;
134
135 qapars[0].push_back(nmean);
136 qapars[1].push_back(nrms);
137 qapars[2].push_back(badfrac);
138 }
139
140 if (badwire) {
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);
145 }
146
147
148 if (m_isMakePlots) {
149 //1. plot bad and good wire plots.
150 plotWireDist(vbadwires, vhitvar);
151
152 //2. plots wire status map
153 plotBadWireMap(vbadwires, vdeadwires);
154
155 //3. plot control parameters histograms
156 plotQaPars(qapars);
157
158 //4. plot statistics related histograms
160 }
161
162 // Save payloads
163 B2INFO("dE/dx Badwire Calibration done: " << vdefectwires.size() << " wires");
164 CDCDedxBadWires* c_badwires = new CDCDedxBadWires(vdefectwires);
165 saveCalibration(c_badwires, "CDCDedxBadWires");
166
167 m_suffix.clear();
168
169 return c_OK;
170}
171
172//------------------------------------
174{
175
176 int cruns = 0;
177 for (auto expRun : getRunList()) {
178 if (cruns == 0) B2INFO("CDCDedxBadWires: start exp " << expRun.first << " and run " << expRun.second << "");
179 cruns++;
180 }
181
182 const auto erStart = getRunList()[0];
183 int estart = erStart.first;
184 int rstart = erStart.second;
185
186 const auto erEnd = getRunList()[cruns - 1];
187 int rend = erEnd.second;
188
189 updateDBObjPtrs(1, rstart, estart);
190
191 if (m_suffix.length() > 0) m_suffix = Form("%s_e%d_r%dr%d", m_suffix.data(), estart, rstart, rend);
192 else m_suffix = Form("e%d_r%dr%d", estart, rstart, rend);
193}
194
195
196//------------------------------------
197void CDCDedxBadWireAlgorithm::plotWireDist(const vector<double>& inwires,
198 map<int, vector<double>>& vhitvar)
199{
200
201 TList* bdlist = new TList();
202 bdlist->SetName("badwires");
203
204 TList* goodlist = new TList();
205 goodlist->SetName("goodwires");
206
207 TList* hflist = new TList();
208 hflist->SetName("highfracwires");
209
210 for (unsigned int jw = 0; jw < c_nwireCDC; ++jw) {
211
212 TH1D* hvar = new TH1D(Form("%s_wire%d", m_suffix.data(), jw), "", m_varBins, m_varMin, m_varMax);
213
214 TH1D* hvarhf = new TH1D(Form("hf%s_wire%d", m_suffix.data(), jw), "", m_varBins, m_varMin, m_varMax);
215 hvarhf->SetTitle(Form("%s, wire = %d; %s; entries", m_suffix.data(), jw, m_varName.data()));
216
217 int ncount = 0, tcount = 0;
218
219 for (unsigned int jh = 0; jh < vhitvar[jw].size(); ++jh) {
220 double jvalue = vhitvar[jw][jh];
221 if (jvalue < m_varMax) {
222 ncount++;
223 hvar->Fill(jvalue);
224 } else {
225 tcount++;
226 if (jvalue < m_varMax * 10.) hvarhf->Fill(jvalue / 10.);
227 }
228 }
229
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));
233
234 bool isbad = false;
235 if (count(inwires.begin(), inwires.end(), jw)) isbad = true;
236
237 double oldwg = m_DBWireGains->getWireGain(jw);
238 if (oldwg == 0) {
239 hvar->SetLineWidth(2);
240 hvar->SetLineColor(kRed);
241 }
242
243 if (isbad) {
244 bdlist->Add(hvar);
245 hflist->Add(hvarhf);
246 } else {
247 if (hvar->Integral() > 100) goodlist->Add(hvar);
248 }
249 }
250
251 printCanvas(bdlist, hflist, kYellow - 9);
252 printCanvas(goodlist, hflist, kGreen);
253
254 delete bdlist;
255 delete goodlist;
256 delete hflist;
257}
258
259//------------------------------------
260void CDCDedxBadWireAlgorithm::printCanvas(TList* list, TList* hflist, Color_t color)
261{
262
263 string listname = list->GetName();
264 string sfx = Form("%s_%s", listname.data(), m_suffix.data());
265
266 TCanvas ctmp(Form("cdcdedx_%s", sfx.data()), "", 1200, 1200);
267 ctmp.Divide(4, 4);
268 ctmp.SetBatch(kTRUE);
269
270 stringstream psname;
271 psname << Form("cdcdedx_bdcal_%s.pdf[", sfx.data());
272 ctmp.Print(psname.str().c_str());
273 psname.str("");
274 psname << Form("cdcdedx_bdcal_%s.pdf", sfx.data());
275
276 for (int ih = 0; ih < list->GetSize(); ih++) {
277
278 TH1D* hist = (TH1D*)list->At(ih);
279
280 double frac = stod(hist->GetYaxis()->GetTitle());
281
282 TPaveText* pinfo = new TPaveText(0.40, 0.63, 0.89, 0.89, "NBNDC");
283 pinfo->AddText(Form("#mu: %0.2f(%0.2f#pm%0.2f)", hist->GetMean(), m_amean, m_meanThres * m_amean));
284 pinfo->AddText(Form("#sigma: %0.2f(%0.2f#pm%0.2f)", hist->GetRMS(), m_arms, m_rmsThres * m_arms));
285 pinfo->AddText(Form("N: %0.00f", hist->Integral()));
286 pinfo->AddText(Form("hf: %0.00f%%(%0.00f%%)", frac, m_fracThres * 100));
287 setTextCosmetics(pinfo, 0.04258064);
288
289 ctmp.cd(ih % 16 + 1);
290 hist->GetYaxis()->SetTitle("entries");
291 hist->SetFillColor(color);
292 hist->SetStats(0);
293 hist->Draw();
294 pinfo->Draw("same");
295
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);
300 histhf->SetStats(0);
301 histhf->Draw("same");
302 }
303
304 if (((ih + 1) % 16 == 0) || ih == (list->GetSize() - 1)) {
305 ctmp.Print(psname.str().c_str());
306 ctmp.Clear("D");
307 }
308 }
309
310 psname.str("");
311 psname << Form("cdcdedx_bdcal_%s.pdf]", sfx.data());
312 ctmp.Print(psname.str().c_str());
313}
314
315//------------------------------------
316void CDCDedxBadWireAlgorithm::plotBadWireMap(const vector<double>& vbadwires, const vector<double>& vdeadwires)
317{
318
319 TCanvas cmap(Form("cmap_%s", m_suffix.data()), "", 800, 800);
320 cmap.SetTitle("CDC dE/dx bad wire status");
321
322 int total = 0;
323 TH2F* hxyAll = getHistoPattern(vbadwires, "all", total);
324 hxyAll->SetTitle(Form("wire status map (%s)", m_suffix.data()));
325 setHistCosmetics(hxyAll, kGray);
326 hxyAll->Draw();
327
328 int nbad = 0.0;
329 TH2F* hxyBad = getHistoPattern(vbadwires, "bad", nbad);
330 if (hxyBad) {
331 setHistCosmetics(hxyBad, kRed);
332 hxyBad->Draw("same");
333 }
334
335 int ndead = 0.0;
336 TH2F* hxyDead = getHistoPattern(vdeadwires, "dead", ndead);
337 if (hxyDead) {
338 setHistCosmetics(hxyDead, kBlack);
339 hxyDead->Draw("same");
340 }
341
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");
349 leg->Draw();
350
351 gStyle->SetLegendTextSize(0.025);
352 TPaveText* pt = new TPaveText(-0.30, -1.47, -0.31, -1.30, "br");
353 setTextCosmetics(pt, 0.02258064);
354
355 TText* text = pt->AddText("CDC-wire map: counter-clockwise and start from +x");
356 text->SetTextColor(kGray + 1);
357 pt->Draw("same");
358
359 cmap.SaveAs(Form("cdcdedx_bdcal_wiremap_%s.pdf", m_suffix.data()));
360
361 delete hxyAll;
362 delete hxyBad;
363 delete hxyDead;
364}
365
366//------------------------------------
367TH2F* CDCDedxBadWireAlgorithm::getHistoPattern(const vector<double>& inwires, const string& suffix, int& total)
368{
369
370 B2INFO("Creating CDCGeometryPar object");
372
373 TH2F* temp = new TH2F(Form("temp_%s_%s", m_suffix.data(), suffix.data()), "", 2400, -1.2, 1.2, 2400, -1.2, 1.2);
374
375 int jwire = -1;
376 total = 0;
377 for (unsigned int ilay = 0; ilay < c_maxNSenseLayers; ++ilay) {
378 for (unsigned int iwire = 0; iwire < cdcgeo.nWiresInLayer(ilay); ++iwire) {
379 jwire++;
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") {
385 total++;
386 temp->Fill(x, y);
387 } else {
388 if (count(inwires.begin(), inwires.end(), jwire)) {
389 temp->Fill(x, y);
390 total++;
391 }
392 }
393 }
394 }
395 return temp;
396}
397
398//------------------------------------
399void CDCDedxBadWireAlgorithm::plotQaPars(map<int, vector<double>>& qapars)
400{
401
402 string qaname[3] = {"mean", "rms", "high_fraction"};
403
404 double linemin[3] = {m_amean* (1 - m_meanThres), m_arms* (1 - m_rmsThres), m_fracThres * 100};
405 double linemax[3] = {m_amean* (1 + m_meanThres), m_arms* (1 + m_rmsThres), m_fracThres * 100};
406
407 for (int iqa = 0; iqa < 3; iqa++) {
408
409 TH1D histqa(Form("%s_%s", qaname[iqa].data(), m_suffix.data()), "", c_nwireCDC, -0.5, 14335.5);
410
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]);
414 }
415
416 TCanvas c_pars(Form("c_pars_%d", iqa), "", 800, 600);
417 c_pars.cd();
418 gPad->SetGridy();
419
420 histqa.SetTitle(Form("%s vs wires (%s); wire ; %s", qaname[iqa].data(), m_suffix.data(), qaname[iqa].data()));
421 histqa.SetStats(0);
422 histqa.Draw();
423
424 TLine* lmin = new TLine(-0.5, linemin[iqa], 14335.5, linemin[iqa]);
425 lmin->SetLineColor(kRed);
426 lmin->Draw("same");
427 TLine* lmax = new TLine(-0.5, linemax[iqa], 14335.5, linemax[iqa]);
428 lmax->SetLineColor(kRed);
429 lmax->Draw("same");
430
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()));
433
434 delete lmax;
435 delete lmin;
436 }
437}
438
439//------------------------------------
441{
442
443 TCanvas cstats("cstats", "cstats", 1000, 500);
444 cstats.SetBatch(kTRUE);
445 cstats.Divide(2, 1);
446
447 cstats.cd(1);
448 auto hestats = getObjectPtr<TH1I>("hestats");
449 if (hestats) {
450 hestats->SetName(Form("htstats_%s", m_suffix.data()));
451 hestats->SetStats(0);
452 hestats->DrawCopy("");
453 }
454
455 cstats.cd(2);
456 auto htstats = getObjectPtr<TH1I>("htstats");
457 if (htstats) {
458 htstats->SetName(Form("htstats_%s", m_suffix.data()));
459 htstats->SetStats(0);
460 htstats->DrawCopy("");
461 }
462
463 cstats.Print(Form("cdcdedx_bdcal_qastats_%s.pdf", m_suffix.data()));
464}
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
Definition beamHelpers.h:28
Abstract base class for different kinds of events.
STL namespace.