Belle II Software  release-05-02-19
CDCDedxWireGainAlgorithm.cc
1 /**************************************************************************
2 * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2016 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: jikumar, jvbennett *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <reconstruction/calibration/CDCDedxWireGainAlgorithm.h>
12 #include <TCanvas.h>
13 #include <TMath.h>
14 #include <TLine.h>
15 #include <TStyle.h>
16 #include <algorithm>
17 #include <iostream>
18 #include <fstream>
19 #include <TH1I.h>
20 #include <TPaveText.h>
21 #include <TText.h>
22 #include <TLegend.h>
23 using namespace Belle2;
24 
25 //-----------------------------------------------------------------
26 // Implementation
27 //-----------------------------------------------------------------
28 CDCDedxWireGainAlgorithm::CDCDedxWireGainAlgorithm() :
29  CalibrationAlgorithm("CDCDedxElectronCollector"),
30  m_badWireFPath(""),
31  m_badWireFName(""),
32  isMakePlots(true),
33  isMergePayload(true),
34  isLTruc(false),
35  isLayerScale(true),
36  fdEdxBins(250),
37  fdEdxMin(0.0),
38  fdEdxMax(5.0),
39  fTrucMin(0.05),
40  fTrucMax(0.75),
41  fStartRun(0)
42 {
43  // Set module properties
44  setDescription("A calibration algorithm for CDC dE/dx wire gains");
45 }
46 
47 
48 //-----------------------------------------------------------------
49 // Run the calibration
50 //-----------------------------------------------------------------
51 CalibrationAlgorithm::EResult CDCDedxWireGainAlgorithm::calibrate()
52 {
53 
54  // Get data objects
55  auto ttree = getObjectPtr<TTree>("tree");
56  if (ttree->GetEntries() < 100)return c_NotEnoughData;
57 
58  const auto expRun = getRunList()[0];
59  updateDBObjPtrs(1, expRun.second, expRun.first);
60  fStartRun = expRun.second;
61 
62  //checking existing payload pointer
63  if (!m_DBBadWires.isValid())
64  B2FATAL("There is no valid payload for CDCDedxBadWires");
65 
66  if (!m_DBWireGains.isValid())
67  B2FATAL("There is no valid payload for CDCDedxWireGain");
68 
69  std::vector<int>* wire = 0;
70  std::vector<double>* dedxhit = 0;
71  ttree->SetBranchAddress("wire", &wire);
72  ttree->SetBranchAddress("dedxhit", &dedxhit);
73 
74  // dedxhit vector to store dE/dx values for each wire
75  std::vector<std::vector<double>> wirededx(14336, std::vector<double>());
76  for (int i = 0; i < ttree->GetEntries(); ++i) {
77  ttree->GetEvent(i);
78  for (unsigned int j = 0; j < wire->size(); ++j) {
79  wirededx[wire->at(j)].push_back(dedxhit->at(j));
80  }
81  }
82 
83  //25-75 is estimated values for 250 bins
84  Int_t lBinILayer = 25, hBinILayer = 75;
85  Int_t lBinOLayer = 25, hBinOLayer = 75;
86 
87  if (!isLTruc) {
88 
89  TH1D* hILayer = new TH1D(Form("hILayer_%d", fStartRun), Form("Inner Layer (start run: %d);dedxhit;entries", fStartRun), fdEdxBins,
90  fdEdxMin, fdEdxMax);
91  TH1D* hOLayer = new TH1D(Form("hOLayer_%d", fStartRun), Form("Outer Layer (start run: %d);dedxhit;entries", fStartRun), fdEdxBins,
92  fdEdxMin, fdEdxMax);
93 
94  for (unsigned int jwire = 0; jwire < 14336; ++jwire) {
95  for (unsigned int jdedxhit = 0; jdedxhit < wirededx[jwire].size(); ++jdedxhit) {
96  double ihitdedx = wirededx[jwire][jdedxhit];
97  if (jwire > 0 && jwire < 160 * 8)hILayer->Fill(ihitdedx);
98  else hOLayer->Fill(ihitdedx);
99  }
100  }
101  getTrucationBins(hILayer, lBinILayer, hBinILayer);
102  getTrucationBins(hOLayer, lBinOLayer, hBinOLayer);
103 
104  if (isMakePlots) {
105  TCanvas* ctem = new TCanvas("Layerhisto", "Inner and Outer Layer dedxhit dist", 900, 400);
106  ctem->Divide(2, 1);
107  ctem->cd(1);
108  hILayer->SetFillColor(kYellow);
109  double lowedge = hILayer->GetXaxis()->GetBinLowEdge(lBinILayer);
110  double upedge = hILayer->GetXaxis()->GetBinUpEdge(hBinILayer);
111  hILayer->SetTitle(Form("%s, trunc range: %0.02f - %0.02f", hILayer->GetTitle(), lowedge, upedge));
112  hILayer->Draw("histo");
113  TH1D* hILayerClone = (TH1D*)hILayer->Clone(Form("hILClone_frun%d", fStartRun));
114  hILayerClone->GetXaxis()->SetRange(lBinILayer, hBinILayer);
115  hILayerClone->SetFillColor(kAzure + 1);
116  hILayerClone->Draw("same histo");
117  ctem->cd(2);
118  lowedge = hOLayer->GetXaxis()->GetBinLowEdge(lBinOLayer);
119  upedge = hOLayer->GetXaxis()->GetBinUpEdge(hBinOLayer);
120  hOLayer->SetTitle(Form("%s trunc range: %0.02f - %0.02f", hOLayer->GetTitle(), lowedge, upedge));
121  hOLayer->SetFillColor(kYellow);
122  hOLayer->Draw("histo");
123  TH1D* hOLayerClone = (TH1D*)hOLayer->Clone(Form("hOLClone_frun%d", fStartRun));
124  hOLayerClone->GetXaxis()->SetRange(lBinOLayer, hBinOLayer);
125  hOLayerClone->SetFillColor(kAzure + 1);
126  hOLayerClone->Draw("same histo");
127  ctem->SaveAs(Form("cdcdedx_wiregain_layerdists_frun%d.pdf", fStartRun));
128  delete ctem;
129  }
130  }
131 
132 
133  TCanvas* ctmp = new TCanvas("tmp", "tmp", 1200, 1200);
134  std::stringstream psname; psname << Form("cdcdedx_wiregain_wiredists_frun%d.pdf[", fStartRun);
135  if (isMakePlots) {
136  ctmp->Divide(4, 4);
137  ctmp->SetBatch(kTRUE);
138  ctmp->Print(psname.str().c_str());
139  psname.str(""); psname << Form("cdcdedx_wiregain_wiredists_frun%d.pdf", fStartRun);
140  }
141 
142  //initialisation of wire gains
143  double iWireTruncMean[14336];
144  for (Int_t jwire = 0; jwire < 14336; jwire++) {
145  iWireTruncMean[jwire] = 1.0;
146  }
147 
148  //Calculations part
149  TH1D* htempPerWire = new TH1D(Form("htempPerWire_frun%d", fStartRun), "blah-blah", fdEdxBins, fdEdxMin, fdEdxMax);
150 
151  for (unsigned int jwire = 0; jwire < 14336; ++jwire) {
152 
153  htempPerWire->SetName(Form("htempPerWire_%d_run%d", jwire, fStartRun));
154  htempPerWire->SetTitle(Form("dedxhit-dist, wire # = %d (start run: %d);dedxhit;entries", jwire, fStartRun));
155  for (unsigned int jdedxhit = 0; jdedxhit < wirededx[jwire].size(); ++jdedxhit) {
156  htempPerWire->Fill(wirededx[jwire][jdedxhit]);
157  }
158 
159  double truncMean = 1.0;
160  int startfrom = 1, endat = 1;
161 
162  if (htempPerWire->Integral() <= 50 || m_DBBadWires->getBadWireStatus(jwire)) {
163  truncMean = 0.0; //dead or bad wire
164  } else if (htempPerWire->Integral() > 50 && htempPerWire->Integral() < 1000) {
165  truncMean = 1.0; //partial dead wire
166  } else {
167  if (!isLTruc) {
168  if (jwire < 160 * 8) {
169  startfrom = lBinILayer; endat = hBinILayer;
170  } else {
171  startfrom = lBinOLayer; endat = hBinOLayer;
172  }
173  } else {
174  getTrucationBins(htempPerWire, startfrom, endat);
175  }
176 
177  double binweights = 0.0, sumofbc = 0.0;
178  for (int ibin = startfrom; ibin <= endat; ibin++) {
179  if (htempPerWire->GetBinContent(ibin) > 0) {
180  binweights += (htempPerWire->GetBinContent(ibin) * htempPerWire->GetBinCenter(ibin));
181  sumofbc += htempPerWire->GetBinContent(ibin);
182  }
183  }
184  if (sumofbc > 0)truncMean = binweights / sumofbc;
185  else truncMean = 1.0;
186  }
187 
188  if (truncMean < 0)truncMean = 0.0; // not <=0 as 0 is reserved for dead wires
189  iWireTruncMean[jwire] = truncMean;
190 
191  if (isMakePlots) {
192  ctmp->cd(jwire % 16 + 1);
193  htempPerWire->SetFillColor(kYellow);
194  htempPerWire->SetTitle(Form("%s, rel. #mu_{truc} %0.04f", htempPerWire->GetTitle(), iWireTruncMean[jwire]));
195  htempPerWire->DrawCopy("hist");
196  TH1D* htempPerWireClone = (TH1D*)htempPerWire->Clone(Form("htempPerWireClone_%d_frun%d", jwire, fStartRun));
197  htempPerWireClone->GetXaxis()->SetRange(startfrom, endat);
198  htempPerWireClone->SetFillColor(kAzure + 1);
199  htempPerWireClone->Draw("same histo");
200  // printf("--> wire %u, rel gain = %0.04f \n", jwire, iWireTruncMean[jwire]);
201  if ((jwire + 1) % 16 == 0)ctmp->Print(psname.str().c_str());
202  }
203  htempPerWire->Reset();
204  }//end of relative gains
205 
206  delete htempPerWire;
207  if (isMakePlots) {
208  psname.str(""); psname << Form("cdcdedx_wiregain_wiredists_frun%d.pdf]", fStartRun);
209  ctmp->Print(psname.str().c_str());
210  delete ctmp;
211 
212  TCanvas* cstats = new TCanvas("cstats", "cstats", 1000, 500);
213  cstats->SetBatch(kTRUE);
214  cstats->Divide(2, 1);
215  cstats->cd(1);
216  auto hestats = getObjectPtr<TH1I>("hestats");
217  if (hestats) {
218  hestats->SetName(Form("hestats_frun%d", fStartRun));
219  hestats->SetStats(0);
220  hestats->DrawCopy("");
221  }
222  cstats->cd(2);
223  auto htstats = getObjectPtr<TH1I>("htstats");
224  if (htstats) {
225  hestats->SetName(Form("htstats_frun%d", fStartRun));
226  htstats->DrawCopy("");
227  hestats->SetStats(0);
228  }
229  cstats->Print(Form("cdcdedx_wiregain_stats_frun%d.pdf", fStartRun));
230  delete cstats;
231 
232  }
233 
234  //changing to vector array
235  std::vector<double> dedxTruncmean;
236  for (Int_t jwire = 0; jwire < 14336; jwire++) {
237  dedxTruncmean.push_back(iWireTruncMean[jwire]);
238  }
239 
240  generateNewPayloads(dedxTruncmean);
241  return c_OK;
242 }
243 
244 
245 void CDCDedxWireGainAlgorithm::generateNewPayloads(std::vector<double> dedxTruncmean)
246 {
247 
248  const auto expRun = getRunList()[0];
249  updateDBObjPtrs(1, expRun.second, expRun.first);
250 
251  if (isMergePayload) {
252  //previous dead but active now ==> abs/merge gain = rel (this is abs by default)
253  //previous dead and dead now ==> abs/merge gain = 0.0
254  //previous active and active now ==> abs/merge gain = pre*rel
255  //previous active but dead now ==> abs/merge gain = 0.0
256  //bool refchange = m_DBWireGains.hasChanged(); //for future or manjor processing
257  B2INFO("Saving merged wiregains for (Exp, Run) : (" << expRun.first << "," << expRun.second << ")");
258  for (unsigned int iwire = 0; iwire < 14336; iwire++) {
259  double pre = m_DBWireGains->getWireGain(iwire);
260  double rel = dedxTruncmean.at(iwire);
261  if (pre != 0.0)dedxTruncmean.at(iwire) *= (double)m_DBWireGains->getWireGain(iwire); //merged or 0.0
262  B2INFO("WG for wire [" << iwire << "], previous = " << pre << ", relative = " << rel << ", merged = " << dedxTruncmean.at(iwire));
263  }
264  } else {
265  B2INFO("Saving relative wiregains for (Exp, Run) : (" << expRun.first << "," << expRun.second << ")");
266  }
267 
268  double layeravg = 1.0;
269  if (isLayerScale) {
270  layeravg = getLayerAverage(dedxTruncmean);
271  for (unsigned int iwire = 0; iwire < 14336; iwire++) {
272  //merged or 0.0 for all bad/dead wire
273  if (layeravg != 0.0) {
274  dedxTruncmean.at(iwire) /= layeravg;
275  }
276  }
277  }
278 
279  //saving final constants in a histograms for validation
280  if (isMakePlots) {
281 
282  TCanvas* cLConst = new TCanvas("cLConst", "cLConst", 1600, 1000);
283  std::stringstream psnameL; psnameL << Form("cdcdedx_wiregain_layerconst_frun%d.pdf[", fStartRun);
284  cLConst->Divide(2, 2);
285  cLConst->SetBatch(kTRUE);
286  cLConst->Print(psnameL.str().c_str());
287  psnameL.str(""); psnameL << Form("cdcdedx_wiregain_layerconst_frun%d.pdf", fStartRun);
288 
289  TH1D* hWGConst = new TH1D(Form("hWGConst_frun%d", fStartRun), Form("wiregain constant (start run: %d); wire numbers;<dedxhit>",
290  fStartRun), 14336, -0.5, 14335.5);
291  if (isMergePayload)hWGConst->SetTitle(Form("abs-const: %s", hWGConst->GetTitle()));
292  else hWGConst->SetTitle(Form("rel-const: %s", hWGConst->GetTitle()));
293 
294  TH1D* hWGConstVar = new TH1D(Form("hWGConstVar_frun%d", fStartRun), Form("wiregain variation (start run: %d); wire gains; nentries",
295  fStartRun), 400, -0.5, 3.5);
296  if (isMergePayload)hWGConstVar->SetTitle(Form("abs-const: %s", hWGConstVar->GetTitle()));
297  else hWGConstVar->SetTitle(Form("rel-const: %s", hWGConstVar->GetTitle()));
298 
299  std::ofstream fBadWG;
300  fBadWG.open(Form("cdcdedx_wiregain_badwire_frun%d.txt", fStartRun));
301 
302  std::ofstream fDeadWG_New;
303  fDeadWG_New.open(Form("cdcdedx_wiregain_deadwire_frun%d_new.txt", fStartRun));
304 
305  std::ofstream fDeadWG_Old;
306  fDeadWG_Old.open(Form("cdcdedx_wiregain_deadwire_frun%d_old.txt", fStartRun));
307 
308  int toWire = 0, countwire = 0;
309  int cnewdeadwires = 0, colddeadwires = 0, cbadwires = 0;
310  for (int iLayer = 0; iLayer < 56; iLayer++) {
311 
312  int iSuperLayer = (iLayer - 2) / 6;
313  if (iSuperLayer <= 0)iSuperLayer = 1; //hack for wire#
314  int nWireiLayer = 160 + (iSuperLayer - 1) * 32;
315 
316  int fromWire = countwire;
317  toWire = toWire + nWireiLayer;
318 
319  TH1D* hLayerConst = new TH1D(Form("hWireConst_L%d_frun%d", iLayer, fStartRun), "blah-blah", nWireiLayer, fromWire * 1.0,
320  toWire * 1.0);
321  if (isMergePayload)hLayerConst->SetTitle(Form("abs-const: Layer = %d (start run: %d); wire numbers;<dedxhit>", iLayer, fStartRun));
322  else hLayerConst->SetTitle(Form("rel-const: Layer = %d (start run: %d); wire numbers;<dedxhit>", iLayer, fStartRun));
323 
324  int iwire = 0, cnewdeadLwires = 0, colddeadLwires = 0, cbadLwires = 0;
325  for (int jwire = fromWire; jwire < toWire; jwire++) {
326 
327  countwire++;
328  iwire++;
329 
330  hLayerConst->SetBinContent(iwire, dedxTruncmean.at(jwire));
331  if (iLayer < 32 && (iwire % 10 == 0))hLayerConst->GetXaxis()->SetBinLabel(iwire, Form("w%d", jwire));
332  else if (iLayer >= 32 && (iwire % 15 == 0))hLayerConst->GetXaxis()->SetBinLabel(iwire, Form("w%d", jwire));
333 
334  hWGConstVar->Fill(dedxTruncmean.at(jwire));
335  hWGConst->SetBinContent(countwire, dedxTruncmean.at(jwire));
336 
337  //see new dead wires
338  if (m_DBBadWires->getBadWireStatus(jwire)) {
339  cbadwires++;
340  cbadLwires++;
341  fBadWG << jwire << "\n";
342  }
343 
344  //see new dead wires
345  if (dedxTruncmean.at(jwire) == 0 && !m_DBBadWires->getBadWireStatus(jwire)) {
346  cnewdeadwires++;
347  cnewdeadLwires++;
348  fDeadWG_New << jwire << "\n";
349  }
350 
351  //see old dead wires
352  if (m_DBWireGains->getWireGain(jwire) == 0 && !m_DBBadWires->getBadWireStatus(jwire)) {
353  colddeadwires++;
354  colddeadLwires++;
355  fDeadWG_Old << jwire << "\n";
356  }
357 
358  if ((countwire) % 500 == 0)hWGConst->GetXaxis()->SetBinLabel(countwire, Form("w%d", countwire));
359  }
360 
361  cLConst->cd(iLayer % 4 + 1);
362  gStyle->SetOptStat("ne");
363  double fraclbad = (100.0 * (cnewdeadLwires + cbadLwires)) / nWireiLayer;
364  if (isLayerScale)hLayerConst->SetTitle(Form("%s, avg = %0.04f, nDeadwires = %d->%d(%0.02f%%)", hLayerConst->GetTitle(),
365  flayerAvg.at(iLayer) / layeravg, colddeadLwires + cbadLwires, cnewdeadLwires + cbadLwires, fraclbad));
366  else hLayerConst->SetTitle(Form("%s, nDeadwires = %d->%d(%0.02f%%)", hLayerConst->GetTitle(), colddeadLwires + cbadLwires,
367  cnewdeadLwires + cbadLwires,
368  fraclbad));
369  if (iLayer < 8)hLayerConst->GetYaxis()->SetRangeUser(-0.1, 4.0);
370  else hLayerConst->GetYaxis()->SetRangeUser(-0.1, 2.0);
371  hLayerConst->SetFillColorAlpha(kAzure, 0.10);
372  hLayerConst->LabelsOption("u", "X");
373  hLayerConst->DrawCopy("hist");
374  if (isLayerScale) {
375  TLine* tlc = new TLine();
376  tlc->SetLineColor(kRed);
377  tlc->SetX1(fromWire); tlc->SetX2(toWire);
378  tlc->SetY1(flayerAvg.at(iLayer) / layeravg); tlc->SetY2(flayerAvg.at(iLayer) / layeravg);
379  tlc->DrawClone("same");
380  delete tlc;
381  }
382  if ((iLayer + 1) % 4 == 0)cLConst->Print(psnameL.str().c_str());
383  hLayerConst->Reset();
384  delete hLayerConst;
385  }
386 
387  psnameL.str(""); psnameL << Form("cdcdedx_wiregain_layerconst_frun%d.pdf]", fStartRun);
388  cLConst->Print(psnameL.str().c_str());
389  delete cLConst;
390 
391  TCanvas* cConst = new TCanvas("cConst", "cConst", 900, 500);
392  cConst->cd();
393  cConst->SetGridy(1);
394  hWGConst->LabelsOption("u", "X");
395  hWGConst->GetYaxis()->SetRangeUser(-0.1, hWGConst->GetMaximum() * 1.05);
396  double fracabad = (100.0 * (cnewdeadwires + cbadwires)) / 14336.0;
397  if (isMergePayload)hWGConst->SetTitle(Form("merged %s, nDeadwires = %d->%d (%0.02f%%)", hWGConst->GetTitle(),
398  colddeadwires + cbadwires,
399  cnewdeadwires + cbadwires, fracabad));
400  else hWGConst->SetTitle(Form("relative %s, nDeadwires = %d->%d (%0.02f%%)", hWGConst->GetTitle(), colddeadwires + cbadwires,
401  cnewdeadwires + cbadwires,
402  fracabad));
403  // hWGConst->SetStats(0);
404  hWGConst->LabelsDeflate();
405  hWGConst->Draw("");
406  cConst->SaveAs(Form("cdcdedx_wiregain_allconstants_frun%d.pdf", fStartRun));
407  delete cConst;
408 
409  TCanvas* cConstvar = new TCanvas("cConstvar", "cConstvar", 500, 400);
410  cConstvar->cd();
411  if (isMergePayload)hWGConstVar->SetTitle(Form("merged %s, nDeadwires = %d->%d (%0.02f%%)", hWGConstVar->GetTitle(),
412  colddeadwires + cbadwires,
413  cnewdeadwires + cbadwires,
414  fracabad));
415  else hWGConstVar->SetTitle(Form("relative %s, bad wire = %d->%d (%0.02f)", hWGConstVar->GetTitle(), colddeadwires + cbadwires,
416  cnewdeadwires + cbadwires,
417  fracabad));
418  //hWGConst->SetStats(0);
419  hWGConstVar->SetFillColorAlpha(kAzure, 0.10);
420  hWGConstVar->Draw("");
421  cConstvar->SaveAs(Form("cdcdedx_wiregain_constantsvar_frun%d.pdf", fStartRun));
422  delete cConstvar;
423 
424  fBadWG.close();
425  fDeadWG_New.close();
426  fDeadWG_Old.close();
427 
428  //create nice plots for bad wire status
429  plotBadWires(cnewdeadwires, colddeadwires, cbadwires);
430  }
431 
432  B2INFO("dE/dx Calibration done for " << dedxTruncmean.size() << " CDC wires");
433  CDCDedxWireGain* gains = new CDCDedxWireGain(dedxTruncmean);
434  saveCalibration(gains, "CDCDedxWireGain");
435 }
436 
437 
438 void CDCDedxWireGainAlgorithm::getTrucationBins(TH1D* htemp, int& binlow, int& binhigh)
439 {
440  //calculating truncation average
441  double TotalInt = htemp->Integral();
442  if (TotalInt <= 0 || htemp->GetNbinsX() <= 0) {
443  binlow = 1.0; binhigh = 1.0;
444  return;
445  } else {
446  binlow = 1.0; binhigh = 1.0;
447  double sumPer5 = 0.0, sumPer75 = 0.0;
448  for (int ibin = 1; ibin <= htemp->GetNbinsX(); ibin++) {
449 
450  if (sumPer5 <= fTrucMin * TotalInt) {
451  sumPer5 += htemp->GetBinContent(ibin);
452  binlow = ibin;
453  }
454 
455  if (sumPer75 <= fTrucMax * TotalInt) {
456  sumPer75 += htemp->GetBinContent(ibin);
457  binhigh = ibin;
458  }
459  }
460  }
461 }
462 
463 double CDCDedxWireGainAlgorithm::getLayerAverage(std::vector<double> tempWire)
464 {
465  //calculating layer average
466  double jLayerMeanAvg = 0.0, OutLayerMeanSum = 0.0;
467  int toWire = 0;
468  int countwire = 0, OutActLayer = 0;
469 
470  TH1D* hLayerAvg = new TH1D(Form("hLayerAvg_frun%d", fStartRun),
471  Form("Layer vs trunc mean avg (star run: %d); layer numbers;<dedxhit>", fStartRun), 56, -0.5, 55.5);
472 
473  for (int iLayer = 0; iLayer < 56; iLayer++) {
474 
475  int iSuperLayer = (iLayer - 2) / 6;
476  if (iSuperLayer <= 0)iSuperLayer = 1; //hack for wire#
477  int nWireiLayer = 160 + (iSuperLayer - 1) * 32;
478 
479  int fromWire = countwire; // or towire+countwire
480  toWire = toWire + nWireiLayer;
481 
482  double jLayerMeanSum = 0.0;
483  int jLayerActWires = 0;
484  std::cout << "iLayer = " << iLayer << ", from wire = " << fromWire << ", to wire = " << toWire << std::endl;
485  for (int jwire = fromWire; jwire < toWire; jwire++) {
486  countwire++;
487  if (tempWire.at(jwire) > 0.) { //active wire only
488  jLayerMeanSum += tempWire.at(jwire);
489  jLayerActWires++;
490  } else continue;
491  }
492 
493  if (jLayerActWires > 0)jLayerMeanAvg = jLayerMeanSum / jLayerActWires;
494  else jLayerMeanAvg = 0.0;
495 
496  flayerAvg.push_back(jLayerMeanAvg);
497  hLayerAvg->SetBinContent(iLayer + 1, jLayerMeanAvg);
498  if ((iLayer + 1) % 2 == 0)hLayerAvg->GetXaxis()->SetBinLabel(iLayer + 1, Form("L%d", iLayer));
499  std::cout << " \t --> sum = " << jLayerMeanSum << ", active wires = " << jLayerActWires
500  << ", average = " << jLayerMeanAvg << std::endl;
501 
502  if (iLayer >= 8 && jLayerMeanAvg > 0.0) {
503  OutLayerMeanSum += jLayerMeanAvg;
504  OutActLayer++; //if any layer is completely dead
505  }
506  }
507 
508  double OutLayerMeanAvg = 1.0;
509  if (OutActLayer > 0) OutLayerMeanAvg = OutLayerMeanSum / OutActLayer;
510 
511  if (isMakePlots) {
512  TCanvas* cLAvg = new TCanvas("clayerAvg", "clayerAvg", 800, 500);
513  cLAvg->SetGridy(1);
514  cLAvg->cd();
515  gStyle->SetOptStat("ne");
516  hLayerAvg->LabelsOption("u", "X");
517  hLayerAvg->SetLineColor(kBlue);
518  hLayerAvg->GetYaxis()->SetRangeUser(-0.1, hLayerAvg->GetMaximum() * 1.20);
519  if (isMergePayload)hLayerAvg->SetTitle(Form("%s, avg = %0.04f (abs)", hLayerAvg->GetTitle(), OutLayerMeanAvg));
520  else hLayerAvg->SetTitle(Form("%s, avg = %0.04f (rel)", hLayerAvg->GetTitle(), OutLayerMeanAvg));
521  hLayerAvg->LabelsDeflate();
522  hLayerAvg->Draw("");
523  TLine* tl = new TLine();
524  tl->SetLineColor(kRed);
525  tl->SetX1(-0.5); tl->SetX2(55.5);
526  tl->SetY1(OutLayerMeanAvg); tl->SetY2(OutLayerMeanAvg);
527  tl->DrawClone("same");
528  cLAvg->SaveAs(Form("cdcdedx_wiregain_layeravg_frun%d.pdf", fStartRun));
529  delete cLAvg;
530  }
531  return OutLayerMeanAvg;
532 }
533 
534 void CDCDedxWireGainAlgorithm::plotBadWires(int nDeadwires, int oDeadwires, int Badwires)
535 {
536 
537  TCanvas* cCDCWires = new TCanvas(Form("cCDCWires_frun%d", fStartRun), "CDC dE/dx bad wire status", 800, 800);
538  TH2F* hxyAll = getHistoPattern("all", "allwires");
539  hxyAll->SetMarkerStyle(20);
540  hxyAll->SetMarkerSize(0.2);
541  hxyAll->SetMarkerColor(kGray);
542  hxyAll->SetStats(0);
543  hxyAll->Draw();
544 
545 
546  TH2F* hxyBad = getHistoPattern(Form("cdcdedx_wiregain_badwire_frun%d.txt", fStartRun), "badwire");
547  if (hxyBad) {
548  hxyBad->SetTitle("");
549  hxyBad->SetMarkerStyle(20);
550  hxyBad->SetMarkerSize(0.3);
551  hxyBad->SetMarkerColor(kBlue);
552  hxyBad->SetStats(0);
553  hxyBad->Draw("same");
554  }
555 
556  TH2F* hxyOldDead = getHistoPattern(Form("cdcdedx_wiregain_deadwire_frun%d_old.txt", fStartRun), "olddead");
557  if (hxyOldDead) {
558  hxyOldDead->SetTitle("");
559  hxyOldDead->SetMarkerStyle(20);
560  hxyOldDead->SetMarkerSize(0.3);
561  hxyOldDead->SetMarkerColor(kBlack);
562  hxyOldDead->SetStats(0);
563  hxyOldDead->Draw("same");
564  }
565 
566  TH2F* hxyNewDead = getHistoPattern(Form("cdcdedx_wiregain_deadwire_frun%d_new.txt", fStartRun), "newdead");
567  if (hxyNewDead) {
568  hxyNewDead->SetTitle("");
569  hxyNewDead->SetMarkerStyle(20);
570  hxyNewDead->SetMarkerSize(0.25);
571  hxyNewDead->SetMarkerColor(kRed);
572  hxyNewDead->SetStats(0);
573  hxyNewDead->Draw("same");
574  }
575 
576  auto legend = new TLegend(0.72, 0.80, 0.90, 0.92);
577  legend->SetBorderSize(0);
578  legend->SetLineWidth(3);
579  legend->SetHeader(Form("Total Bad: %d (~%0.02f%%)", nDeadwires + Badwires, 100.*(nDeadwires + Badwires) / 14336.0));
580  legend->AddEntry(hxyBad, Form("badadc %d" , Badwires), "p");
581  legend->AddEntry(hxyOldDead, Form("olddead %d" , oDeadwires), "p");
582  legend->AddEntry(hxyNewDead, Form("newdead %d" , nDeadwires), "p");
583  legend->Draw();
584 
585  gStyle->SetLegendTextSize(0.025);
586  TPaveText* pt = new TPaveText(-0.30993, -1.470968, -0.3102707, -1.304516, "br");
587  pt->SetFillColor(0);
588  pt->SetFillStyle(3001);
589  pt->SetLineColor(2);
590  pt->SetTextFont(82);
591  pt->SetTextSize(0.02258064);
592  TText* t1 = pt->AddText("CDC-wire map: counter-clockwise and start from +x");
593  t1->SetTextColor(kGray + 1);
594  pt->Draw("same");
595 
596  cCDCWires->SaveAs(Form("cdcdedx_wiregain_wirestatus_frun%d.pdf", fStartRun));
597  delete cCDCWires;
598 }
599 
600 TH2F* CDCDedxWireGainAlgorithm::getHistoPattern(TString badFileName, TString suffix = "")
601 {
602 
603  int wire, nwire, twire;
604  double radius, phi, x, y;
605 
606  std::ifstream infile;
607  infile.open(Form("%s", badFileName.Data()));
608 
609  TH2F* temp = new TH2F(Form("temp_%s_frun%d", suffix.Data(), fStartRun), Form("wire status (start run: %d)", fStartRun), 2400, -1.2,
610  1.2, 2400, -1.2, 1.2);
611  if (!infile.fail()) {
612  int bwires = 0;
613  while (infile >> bwires) {
614  nwire = getIndexVal(bwires, "nwirelayer");
615  twire = getIndexVal(bwires, "twire");
616  radius = getIndexVal(bwires, "rwire");
617  wire = bwires - twire ;
618  phi = 2.*TMath::Pi() * (float(wire) / float(nwire));
619  x = radius * cos(phi);
620  y = radius * sin(phi);
621  temp->Fill(x, y);
622  }
623  } else {
624  for (int iwires = 0; iwires < 14336; iwires++) {
625  nwire = getIndexVal(iwires, "nwirelayer");
626  twire = getIndexVal(iwires, "twire");
627  radius = getIndexVal(iwires, "rwire");
628  wire = iwires - twire ;
629  phi = 2.*TMath::Pi() * (float(wire) / float(nwire));
630  x = radius * cos(phi);
631  y = radius * sin(phi);
632  temp->Fill(x, y);
633  }
634  }
635  return temp;
636 }
637 
638 double CDCDedxWireGainAlgorithm::getIndexVal(int iWire, TString what)
639 {
640  //radius of each CDC layer
641  double r[56] = {
642  16.80, 17.80, 18.80, 19.80, 20.80, 21.80, 22.80, 23.80,
643  25.70, 27.52, 29.34, 31.16, 32.98, 34.80,
644  36.52, 38.34, 40.16, 41.98, 43.80, 45.57,
645  47.69, 49.46, 51.28, 53.10, 54.92, 56.69,
646  58.41, 60.18, 62.00, 63.82, 65.64, 67.41,
647  69.53, 71.30, 73.12, 74.94, 76.76, 78.53,
648  80.25, 82.02, 83.84, 85.66, 87.48, 89.25,
649  91.37, 93.14, 94.96, 96.78, 98.60, 100.37,
650  102.09, 103.86, 105.68, 107.50, 109.32, 111.14
651  };
652 
653  Int_t totalWireiLayer = 0 ;
654  double myreturn = 0;
655  for (Int_t iLayer = 0; iLayer < 56; iLayer++) {
656  int iSuperLayer = (iLayer - 2) / 6;
657  if (iSuperLayer <= 0)iSuperLayer = 1;
658  int nWireiLayer = 160 + (iSuperLayer - 1) * 32;
659  totalWireiLayer += nWireiLayer;
660 
661  if (iWire < totalWireiLayer) {
662  if (what == "layer")myreturn = iLayer;
663  else if (what == "nwirelayer") myreturn = nWireiLayer;
664  else if (what == "twire") myreturn = totalWireiLayer - nWireiLayer;
665  else if (what == "rwire") myreturn = r[iLayer] / 100.;
666  else std::cout << "Invalid return :0 " << std::endl;
667  break;
668  }
669  }
670  return myreturn;
671 }
Belle2::CDCDedxWireGainAlgorithm::getTrucationBins
void getTrucationBins(TH1D *htemp, int &binlow, int &binhigh)
function to get bins of trunction from histogram
Definition: CDCDedxWireGainAlgorithm.cc:438
Belle2::CDCDedxWireGainAlgorithm::fTrucMin
double fTrucMin
min trunc range for mean
Definition: CDCDedxWireGainAlgorithm.h:128
Belle2::CDCDedxWireGainAlgorithm::plotBadWires
void plotBadWires(int nDeadwires, int oDeadwires, int Badwires)
function to plot bad wire status (then and now)
Definition: CDCDedxWireGainAlgorithm.cc:534
Belle2::CalibrationAlgorithm::saveCalibration
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
Definition: CalibrationAlgorithm.cc:290
Belle2::CDCDedxWireGainAlgorithm::isMergePayload
bool isMergePayload
merge payload at the of calibration
Definition: CDCDedxWireGainAlgorithm.h:122
Belle2::CDCDedxWireGainAlgorithm::m_DBBadWires
DBObjPtr< CDCDedxBadWires > m_DBBadWires
Bad wire DB object.
Definition: CDCDedxWireGainAlgorithm.h:134
Belle2::CDCDedxWireGainAlgorithm::getIndexVal
double getIndexVal(int iWire, TString what)
function to return various CDC indexing for given wire
Definition: CDCDedxWireGainAlgorithm.cc:638
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfuly =0 in Python.
Definition: CalibrationAlgorithm.h:51
Belle2::CDCDedxWireGainAlgorithm::isLTruc
bool isLTruc
method of trunc range for mean
Definition: CDCDedxWireGainAlgorithm.h:123
Belle2::CDCDedxWireGainAlgorithm::fdEdxMax
double fdEdxMax
max dedx range for wiregain cal
Definition: CDCDedxWireGainAlgorithm.h:127
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition: CalibrationAlgorithm.h:331
Belle2::CDCDedxWireGainAlgorithm::CDCDedxWireGainAlgorithm
CDCDedxWireGainAlgorithm()
Constructor: Sets the description, the properties and the parameters of the algorithm.
Definition: CDCDedxWireGainAlgorithm.cc:28
Belle2::CalibrationAlgorithm::updateDBObjPtrs
void updateDBObjPtrs(const unsigned int event, const int run, const int experiment)
Updates any DBObjPtrs by calling update(event) for DBStore.
Definition: CalibrationAlgorithm.cc:398
Belle2::CDCDedxWireGainAlgorithm::calibrate
virtual EResult calibrate() override
Wire gain algorithm.
Definition: CDCDedxWireGainAlgorithm.cc:51
Belle2::CDCDedxWireGainAlgorithm::fTrucMax
double fTrucMax
max trunc range for mean
Definition: CDCDedxWireGainAlgorithm.h:129
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::CDCDedxWireGainAlgorithm::generateNewPayloads
void generateNewPayloads(std::vector< double > dedxTruncmean)
function to finally store new payload after full calibration
Definition: CDCDedxWireGainAlgorithm.cc:245
Belle2::CDCDedxWireGainAlgorithm::fdEdxMin
double fdEdxMin
min dedx range for wiregain cal
Definition: CDCDedxWireGainAlgorithm.h:126
Belle2::CDCDedxWireGainAlgorithm::getLayerAverage
double getLayerAverage(std::vector< double > tempWire)
function to get layer avg from outer layers
Definition: CDCDedxWireGainAlgorithm.cc:463
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition: CalibrationAlgorithm.h:50
Belle2::CalibrationAlgorithm::getRunList
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
Definition: CalibrationAlgorithm.h:276
Belle2::CDCDedxWireGainAlgorithm::fdEdxBins
int fdEdxBins
number of bins for dedx histogram
Definition: CDCDedxWireGainAlgorithm.h:125
Belle2::CalibrationAlgorithm::c_NotEnoughData
@ c_NotEnoughData
Needs more data =2 in Python.
Definition: CalibrationAlgorithm.h:53
Belle2::CalibrationAlgorithm
Base class for calibration algorithms.
Definition: CalibrationAlgorithm.h:47
Belle2::CDCDedxWireGainAlgorithm::fStartRun
int fStartRun
boundary start at this run
Definition: CDCDedxWireGainAlgorithm.h:130
Belle2::CDCDedxWireGainAlgorithm::getHistoPattern
TH2F * getHistoPattern(TString badFileName, TString suffix)
function to plot wires in hist with input file
Definition: CDCDedxWireGainAlgorithm.cc:600
Belle2::CDCDedxWireGainAlgorithm::flayerAvg
std::vector< double > flayerAvg
layer wire avg of trun mean
Definition: CDCDedxWireGainAlgorithm.h:131
Belle2::CDCDedxWireGainAlgorithm::m_DBWireGains
DBObjPtr< CDCDedxWireGain > m_DBWireGains
Wire gain DB object.
Definition: CDCDedxWireGainAlgorithm.h:133
Belle2::CDCDedxWireGain
dE/dx wire gain calibration constants
Definition: CDCDedxWireGain.h:35
Belle2::CDCDedxWireGainAlgorithm::isMakePlots
bool isMakePlots
produce plots for status
Definition: CDCDedxWireGainAlgorithm.h:121
Belle2::CDCDedxWireGainAlgorithm::isLayerScale
bool isLayerScale
method of scaling layer avg
Definition: CDCDedxWireGainAlgorithm.h:124