Belle II Software  release-05-01-25
CDCDedxRunGainAlgorithm.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/CDCDedxRunGainAlgorithm.h>
12 
13 #include <TF1.h>
14 #include <TCanvas.h>
15 #include <TLine.h>
16 #include <iostream>
17 #include <fstream>
18 #include <TH1I.h>
19 
20 using namespace Belle2;
21 
22 //-----------------------------------------------------------------
23 // Implementation
24 //-----------------------------------------------------------------
26  CalibrationAlgorithm("CDCDedxElectronCollector"),
27  isMakePlots(true),
28  isMergePayload(true),
29  fsrun(""),
30  fSigLim(2.5),
31  fdEdxBins(1000),
32  fdEdxMin(0.0),
33  fdEdxMax(5.0),
34  fAdjust(1.00)
35 {
36  // Set module properties
37  setDescription("A calibration algorithm for CDC dE/dx run gains");
38 }
39 
40 //-----------------------------------------------------------------
41 // Run the calibration
42 //-----------------------------------------------------------------
44 {
45  B2INFO("Preparing dE/dx calibration for CDC dE/dx run gain");
46 
47  const auto expRun = getRunList()[0];
48  updateDBObjPtrs(1, expRun.second, expRun.first);
49  double ExistingRG = m_DBRunGain->getRunGain();
50 
51  // Get data objects
52  auto ttree = getObjectPtr<TTree>("tree");
53  if (!ttree)return c_NotEnoughData;
54 
55  double dedx;
56  int run = 0, runtemp = 0;
57  TString status = "";
58  ttree->SetBranchAddress("dedx", &dedx);
59  ttree->SetBranchAddress("run", &run);
60 
61  TH1D* hDedx = new TH1D("hDedx", "hDedx;dEdx rel (no hadsat);entries", fdEdxBins, fdEdxMin, fdEdxMax);
62  TH1D* hDedxAbs = new TH1D("hDedxAbs", "hDedxAbs;dEdx abs(no hadsat);entries", fdEdxBins, fdEdxMin, fdEdxMax);
63 
64  for (int i = 0; i < ttree->GetEntries(); ++i) {
65  ttree->GetEvent(i);
66  if (dedx <= 0) continue;
67  hDedx->Fill(dedx);
68  hDedxAbs->Fill(dedx * ExistingRG);
69  runtemp = run;
70  }
71 
72  //min 500 good events otherwise merge to next run
73  if (hDedx->Integral() < 1000) {
74  B2INFO("Not enough data for run (" << runtemp << ") so far only " << hDedx->Integral() << " therefore adding next run");
75  fsrun += Form("%d_", runtemp);
76  delete hDedx;
77  delete hDedxAbs;
78  return c_NotEnoughData;
79  }
80 
81  B2INFO("Total track for this run: " << runtemp << " = " << hDedx->Integral());
82  fsrun += Form("%d", runtemp);
83  hDedx->SetName(Form("%s_%s", hDedx->GetName(), fsrun.Data()));
84  hDedxAbs->SetName(Form("%s_%s", hDedxAbs->GetName(), fsrun.Data()));
85 
86  //Fit absolute run for the reference only
87  double RGabs_Mean = 1.0;
88  FitGaussianWRange(hDedxAbs, status);
89  if (status != "FitOK") {
90  hDedxAbs->SetTitle(Form("dE/dx abs, Run %s (%s)", fsrun.Data(), status.Data()));
91  } else {
92  RGabs_Mean = hDedxAbs->GetFunction("gaus")->GetParameter(1);
93  double RGabs_MeanErr = hDedxAbs->GetFunction("gaus")->GetParError(1);
94  double RGabs_Sigma = hDedxAbs->GetFunction("gaus")->GetParameter(2);
95  hDedxAbs->SetTitle(Form("dE/dx abs, Run %s (%s), #mu_{fit}: %0.04f#pm%0.04f,, #sigma_{fit}: %0.04f", fsrun.Data(), status.Data(),
96  RGabs_Mean, RGabs_MeanErr, RGabs_Sigma));
97  }
98 
99  //Fit relative run for actual calibration
100  double RGrel_Mean = 1.0, RGrel_MeanErr = -99.0, RGrel_Sigma = -99.0;
101  FitGaussianWRange(hDedx, status);
102  if (status != "FitOK") {
103  RGrel_Mean = hDedx->GetMean(); //mean for hist for fit failure
104  hDedx->SetTitle(Form("dE/dx abs, Run %s (%s)", fsrun.Data(), status.Data()));
105  } else {
106  RGrel_Mean = hDedx->GetFunction("gaus")->GetParameter(1);
107  RGrel_MeanErr = hDedx->GetFunction("gaus")->GetParError(1);
108  RGrel_Sigma = hDedx->GetFunction("gaus")->GetParameter(2);
109  hDedx->SetTitle(Form("dE/dx rel, Run %s (%s), #mu_{fit}: %0.04f#pm%0.04f,, #sigma_{fit}: %0.04f", fsrun.Data(), status.Data(),
110  RGrel_Mean, RGrel_MeanErr, RGrel_Sigma));
111  }
112 
113  if (RGrel_Mean <= 0.0)RGrel_Mean = 1.0; //safe button
114 
115  //for validation purpose
116  if (isMakePlots) {
117 
118  std::ofstream fileRGout;
119  fileRGout.open(Form("fRG_Constants_%s.txt", fsrun.Data()));
120  fileRGout << expRun.second << " " << ExistingRG << " " << RGrel_Mean << " " << RGrel_MeanErr << " " << RGrel_Sigma << "\n";
121  fileRGout.close();
122 
123  TCanvas* ctmp = new TCanvas("tmp", "tmp", 1000, 500);
124  ctmp->SetBatch(kTRUE);
125  ctmp->Divide(2, 1);
126 
127  ctmp->cd(1);
128  hDedx->SetFillColorAlpha(kYellow, 0.30);
129  hDedx->DrawCopy("");
130 
131  TLine* tl = new TLine();
132  tl->SetLineColor(kRed);
133  tl->SetX1(RGrel_Mean); tl->SetX2(RGrel_Mean);
134  tl->SetY1(0); tl->SetY2(hDedx->GetMaximum());
135  tl->DrawClone("same");
136 
137  ctmp->cd(2);
138  hDedxAbs->SetFillColorAlpha(kBlue, 0.10);
139  hDedxAbs->DrawCopy("");
140 
141  TLine* tlAbs = new TLine();
142  tlAbs->SetLineColor(kRed);
143  tlAbs->SetX1(RGabs_Mean); tlAbs->SetX2(RGabs_Mean);
144  tlAbs->SetY1(0); tlAbs->SetY2(hDedxAbs->GetMaximum());
145  tlAbs->DrawClone("same");
146 
147  ctmp->Print(Form("cdcdedx_rungain_fits_%s_%s.pdf", fsrun.Data(), status.Data()));
148 
149  TCanvas* cstats = new TCanvas("cstats", "cstats", 1000, 500);
150  cstats->SetBatch(kTRUE);
151  cstats->Divide(2, 1);
152  cstats->cd(1);
153  auto hestats = getObjectPtr<TH1I>("hestats");
154  if (hestats)hestats->DrawCopy("");
155  cstats->cd(2);
156  auto htstats = getObjectPtr<TH1I>("htstats");
157  if (htstats)htstats->DrawCopy("");
158  cstats->Print(Form("cdcdedx_rungain_stats_%s.pdf", fsrun.Data()));
159 
160  delete ctmp;
161  delete tl;
162  delete tlAbs;
163  delete cstats;
164  }
165 
166  B2INFO("Saving new rung for (Exp, Run) : (" << expRun.first << "," << expRun.second << ")");
167  generateNewPayloads(RGrel_Mean, ExistingRG);
168 
169  fsrun.Resize(0);
170  delete hDedx;
171  delete hDedxAbs;
172  return c_OK;
173 }
174 
175 
176 void CDCDedxRunGainAlgorithm::generateNewPayloads(double RGrel_Mean, double ExistingRG)
177 {
178  if (isMergePayload) {
179  B2INFO("--> RunGain: Previous = " << ExistingRG << ", Relative = " << RGrel_Mean << ", Adjustment = " << fAdjust <<
180  ", Merged (saved) = " << RGrel_Mean *
181  ExistingRG * fAdjust);
182  RGrel_Mean *= ExistingRG;
183  RGrel_Mean *= fAdjust; //default is 1.0
184  } else {
185  B2INFO("--> RunGain: Previous = " << ExistingRG << ", Relative (saved) = " << RGrel_Mean);
186  }
187 
188  CDCDedxRunGain* gain = new CDCDedxRunGain(RGrel_Mean);
189  saveCalibration(gain, "CDCDedxRunGain");
190 }
191 
192 
193 void CDCDedxRunGainAlgorithm::FitGaussianWRange(TH1D*& temphist, TString& status)
194 {
195 
196  double histmean = temphist->GetMean();
197  double histrms = temphist->GetRMS();
198  temphist->GetXaxis()->SetRangeUser(histmean - 5.0 * histrms, histmean + 5.0 * histrms);
199 
200  int fs = temphist->Fit("gaus", "Q0");
201  if (fs != 0) {
202  B2INFO(Form("\tFit (round 1) for hist (%s) failed (status = %d)", temphist->GetName(), fs));
203  status = "FitFailed";
204  return;
205  } else {
206  double mean = temphist->GetFunction("gaus")->GetParameter(1);
207  double width = temphist->GetFunction("gaus")->GetParameter(2);
208  temphist->GetXaxis()->SetRangeUser(mean - 5.0 * width, mean + 5.0 * width);
209  fs = temphist->Fit("gaus", "QR", "", mean - fSigLim * width, mean + fSigLim * width);
210  if (fs != 0) {
211  B2INFO(Form("\tFit (round 2) for hist (%s) failed (status = %d)", temphist->GetName(), fs));
212  status = "FitFailed";
213  return;
214  } else {
215  temphist->GetXaxis()->SetRangeUser(mean - 5.0 * width, mean + 5.0 * width);
216  B2INFO(Form("\tFit for hist (%s) sucessfull (status = %d)", temphist->GetName(), fs));
217  status = "FitOK";
218  }
219  }
220 }
Belle2::CDCDedxRunGainAlgorithm::isMakePlots
bool isMakePlots
produce plots for status
Definition: CDCDedxRunGainAlgorithm.h:98
Belle2::CDCDedxRunGainAlgorithm::fdEdxMin
double fdEdxMin
min dedx range for gain cal
Definition: CDCDedxRunGainAlgorithm.h:105
Belle2::CDCDedxRunGainAlgorithm::fSigLim
double fSigLim
fit range limit based on sigma
Definition: CDCDedxRunGainAlgorithm.h:103
Belle2::CDCDedxRunGainAlgorithm::fAdjust
double fAdjust
factor to adjust dedx gain
Definition: CDCDedxRunGainAlgorithm.h:107
Belle2::CDCDedxRunGainAlgorithm::isMergePayload
bool isMergePayload
merge payload at the of calibration
Definition: CDCDedxRunGainAlgorithm.h:99
Belle2::CDCDedxRunGainAlgorithm::fsrun
TString fsrun
flag to indentify low stats runs
Definition: CDCDedxRunGainAlgorithm.h:102
Belle2::CDCDedxRunGain
dE/dx run gain calibration constants
Definition: CDCDedxRunGain.h:33
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::CDCDedxRunGainAlgorithm::CDCDedxRunGainAlgorithm
CDCDedxRunGainAlgorithm()
Constructor: Sets the description, the properties and the parameters of the algorithm.
Definition: CDCDedxRunGainAlgorithm.cc:25
Belle2::CDCDedxRunGainAlgorithm::calibrate
virtual EResult calibrate() override
Run algorithm.
Definition: CDCDedxRunGainAlgorithm.cc:43
Belle2::CDCDedxRunGainAlgorithm::fdEdxBins
int fdEdxBins
number of bins for dedx histogram
Definition: CDCDedxRunGainAlgorithm.h:104
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfuly =0 in Python.
Definition: CalibrationAlgorithm.h:51
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition: CalibrationAlgorithm.h:331
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::CDCDedxRunGainAlgorithm::generateNewPayloads
void generateNewPayloads(double RunGainConst, double ExistingRG)
function to store new payload after full calibration
Definition: CDCDedxRunGainAlgorithm.cc:176
Belle2::CDCDedxRunGainAlgorithm::FitGaussianWRange
void FitGaussianWRange(TH1D *&temphist, TString &status)
function to perform fit run by run
Definition: CDCDedxRunGainAlgorithm.cc:193
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
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::CalibrationAlgorithm::c_NotEnoughData
@ c_NotEnoughData
Needs more data =2 in Python.
Definition: CalibrationAlgorithm.h:53
Belle2::CDCDedxRunGainAlgorithm::fdEdxMax
double fdEdxMax
max dedx range for gain cal
Definition: CDCDedxRunGainAlgorithm.h:106
Belle2::CalibrationAlgorithm
Base class for calibration algorithms.
Definition: CalibrationAlgorithm.h:47
Belle2::CDCDedxRunGainAlgorithm::m_DBRunGain
DBObjPtr< CDCDedxRunGain > m_DBRunGain
Run gain DB object.
Definition: CDCDedxRunGainAlgorithm.h:100