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