Belle II Software development
Public Member Functions | Static Public Member Functions | Private Attributes | List of all members
HadronSaturation Class Reference
Packages » cdc

Class to perform the hadron saturation calibration. More...

#include <HadronSaturation.h>

Public Member Functions

 HadronSaturation ()
 Constructor: Sets the description, the properties and the parameters of the algorithm.
 
 HadronSaturation (double alpha, double gamma, double delta, double power, double ratio, int cosbins)
 set the input variables
 
virtual ~HadronSaturation ()
 Destructor.
 
void fillSample (TString infilename)
 fill the vectors below
 
void printEvents (int firstevent, int nevents)
 print a sample of events
 
void fitSaturation ()
 perform the hadron saturation fit
 
void clear ()
 clear the vectors
 
void setCosBins (int nbins)
 set the number of cosine bins
 
double myFunction (double alpha, double gamma, double delta, double power, double ratio)
 some helper functions for the hadron saturation correction
 

Static Public Member Functions

static void minuitFunction (int &, double *, double &result, double *para, int)
 functions for the hadron saturation correction
 

Private Attributes

int m_cosbins
 the number of cosine bins
 
double m_alpha
 the alpha parameter for the hadron saturation correction
 
double m_gamma
 the gamma parameter for the hadron saturation correction
 
double m_delta
 the delta parameter for the hadron saturation correction
 
double m_power
 the power parameter for the hadron saturation correction
 
double m_ratio
 the ratio parameter for the hadron saturation correction
 
std::vector< double > m_dedx
 a vector to hold dE/dx measurements
 
std::vector< double > m_dedxerror
 a vector to hold dE/dx errors
 
std::vector< double > m_betagamma
 a vector to hold beta-gamma values
 
std::vector< double > m_costheta
 a vector to hold cos(theta) values
 

Detailed Description

Class to perform the hadron saturation calibration.

Definition at line 34 of file HadronSaturation.h.

Constructor & Destructor Documentation

◆ HadronSaturation() [1/2]

HadronSaturation ( )

Constructor: Sets the description, the properties and the parameters of the algorithm.

Definition at line 14 of file HadronSaturation.cc.

14 :
15 m_cosbins(8),
16 m_alpha(0.044286100),
17 m_gamma(0.0019308200),
18 m_delta(0.020),
19 m_power(1.3205300),
20 m_ratio(1.0000)
21{
22 m_dedx.clear();
23 m_dedxerror.clear();
24 m_betagamma.clear();
25 m_costheta.clear();
26 HC_obj = this;
27}
Belle2::HadronSaturation::m_costheta
std::vector< double > m_costheta
a vector to hold cos(theta) values
Definition HadronSaturation.h:101
Belle2::HadronSaturation::m_delta
double m_delta
the delta parameter for the hadron saturation correction
Definition HadronSaturation.h:94
Belle2::HadronSaturation::m_betagamma
std::vector< double > m_betagamma
a vector to hold beta-gamma values
Definition HadronSaturation.h:100
Belle2::HadronSaturation::m_ratio
double m_ratio
the ratio parameter for the hadron saturation correction
Definition HadronSaturation.h:96
Belle2::HadronSaturation::m_dedx
std::vector< double > m_dedx
a vector to hold dE/dx measurements
Definition HadronSaturation.h:98
Belle2::HadronSaturation::m_gamma
double m_gamma
the gamma parameter for the hadron saturation correction
Definition HadronSaturation.h:93
Belle2::HadronSaturation::m_dedxerror
std::vector< double > m_dedxerror
a vector to hold dE/dx errors
Definition HadronSaturation.h:99
Belle2::HadronSaturation::m_alpha
double m_alpha
the alpha parameter for the hadron saturation correction
Definition HadronSaturation.h:92
Belle2::HadronSaturation::m_power
double m_power
the power parameter for the hadron saturation correction
Definition HadronSaturation.h:95
Belle2::HadronSaturation::m_cosbins
int m_cosbins
the number of cosine bins
Definition HadronSaturation.h:90

◆ HadronSaturation() [2/2]

HadronSaturation ( double alpha,
double gamma,
double delta,
double power,
double ratio,
int cosbins )

set the input variables

Definition at line 29 of file HadronSaturation.cc.

29 :
30 m_cosbins(cosbins),
31 m_alpha(alpha),
32 m_gamma(gamma),
33 m_delta(delta),
34 m_power(power),
35 m_ratio(ratio)
36{
37 m_dedx.clear();
38 m_dedxerror.clear();
39 m_betagamma.clear();
40 m_costheta.clear();
41 HC_obj = this;
42}

◆ ~HadronSaturation()

virtual ~HadronSaturation ( )
inlinevirtual

Destructor.

Definition at line 51 of file HadronSaturation.h.

51{ clear(); };

Member Function Documentation

◆ clear()

void clear ( )

clear the vectors

Definition at line 297 of file HadronSaturation.cc.

298{
299
300 m_dedx.clear();
301 m_dedxerror.clear();
302 m_betagamma.clear();
303 m_costheta.clear();
304}

◆ fillSample()

void fillSample ( TString infilename)

fill the vectors below

Definition at line 44 of file HadronSaturation.cc.

45{
46
47 B2INFO("\n\t HadronSaturation: filling sample events...\n");
48
49 // clear the vectors to be filled
50 clear();
51
52 std::vector< TString > types;
53 types.push_back("pion");
54 types.push_back("kaon");
55 types.push_back("proton");
56 types.push_back("muon");
57
58 // fill the containers with a previously prepared sample
59 TFile* satFile = new TFile(infilename);
60
61 for (int i = 0; i < int(types.size()); ++i) {
62
63 TTree* satTree = (TTree*)satFile->Get(types[i]);
64 double satbg, satcosth, satdedx, satdedxerr;
65 satTree->SetBranchAddress("bg", &satbg);
66 satTree->SetBranchAddress("costh", &satcosth);
67 satTree->SetBranchAddress("dedx", &satdedx);
68 satTree->SetBranchAddress("dedxerr", &satdedxerr);
69
70 // fill the vectors
71 for (unsigned int j = 0; j < satTree->GetEntries(); ++j) {
72 satTree->GetEntry(j);
73 if (satdedxerr == 0) continue;
74 m_dedx.push_back(satdedx);
75 m_dedxerror.push_back(satdedxerr);
76 m_betagamma.push_back(satbg);
77 m_costheta.push_back(satcosth);
78 }
79 }
80 satFile->Close();
81}
Belle2::HadronSaturation::clear
void clear()
clear the vectors
Definition HadronSaturation.cc:297

◆ fitSaturation()

void fitSaturation ( )

perform the hadron saturation fit

Definition at line 175 of file HadronSaturation.cc.

176{
177
178 if (m_dedx.empty()) {
179 B2ERROR("HadronSaturation::fitSaturation - no data to fit!");
180 return;
181 }
182
183 B2INFO("\t Performing the hadron saturation fit...");
184
185 // Construct the fitter
186 TFitter* minimizer = new TFitter(5);
187 minimizer->SetFCN(HadronSaturation::minuitFunction);
188
189 auto setPar = [&](int idx, const char* name, double val, double step,
190 double low, double high, bool fix = false) {
191 minimizer->SetParameter(idx, name, val, step, low, high);
192 if (fix) minimizer->FixParameter(idx);
193 };
194
195 // Set initial parameters with reasonable ranges
196 setPar(0, "alpha", m_alpha, gRandom->Rndm() * 0.001, -5.0, 5.0);
197 setPar(1, "gamma", m_gamma, gRandom->Rndm() * 0.001, -5.0, 5.0);
198 setPar(2, "delta", m_delta, gRandom->Rndm() * 0.001, 0, 0.50, true);
199 setPar(3, "power", m_power, gRandom->Rndm() * 0.01, 0.05, 2.5);
200 setPar(4, "ratio", m_ratio, gRandom->Rndm() * 0.01, 0.5, 2, true);
201
202 // Set minuit fitting strategy
203 double arg[10];
204
205 double strategy(2.);
206 minimizer->ExecuteCommand("SET START", &strategy, 1);
207
208 double up(1.);
209 minimizer->ExecuteCommand("SET ERR", &up, 1);
210
211 // Suppress much of the output of MINUIT
212 arg[0] = -1;
213 minimizer->ExecuteCommand("SET PRINT", arg, 1);
214
215 double eps_machine(std::numeric_limits<double>::epsilon());
216 minimizer->ExecuteCommand("SET EPS", &eps_machine, 1);
217
218
219 double maxcalls(5000.), tolerance(0.1);
220 double arglist[] = {maxcalls, tolerance};
221 unsigned int nargs(2);
222
223 for (int i = 0; i < 30; ++i) {
224
225 minimizer->SetParameter(0, "alpha", m_alpha, gRandom->Rndm() * 0.001, -5.0, 5.0);
226 minimizer->SetParameter(1, "gamma", m_gamma, gRandom->Rndm() * 0.001, -5.0, 5.0);
227 minimizer->SetParameter(2, "delta", m_delta, gRandom->Rndm() * 0.001, 0, 0.50);
228 minimizer->SetParameter(3, "power", m_power, gRandom->Rndm() * 0.01, 0.05, 2.5);
229 minimizer->SetParameter(4, "ratio", m_ratio, gRandom->Rndm() * 0.01, 0.5, 2);
230 minimizer->FixParameter(2);
231 minimizer->FixParameter(4);
232
233
234 minimizer->ExecuteCommand("MIGRAD", arglist, nargs);
235 minimizer->ExecuteCommand("MIGRAD", arglist, nargs);
236 int status = minimizer->ExecuteCommand("HESSE", arglist, nargs);
237
238 minimizer->PrintResults(1, 0);
239 B2INFO("\t iter = " << i << ": Fit status: " << status);
240
241 int counter = 0;
242 while (status != 0 && counter < 5) {
243 minimizer->SetParameter(0, "alpha", m_alpha, gRandom->Rndm() * 0.001, -5.0, 5.0);
244 minimizer->SetParameter(1, "gamma", m_gamma, gRandom->Rndm() * 0.001, -5.0, 5.0);
245 minimizer->SetParameter(2, "delta", m_delta, gRandom->Rndm() * 0.001, 0, 0.50);
246 minimizer->SetParameter(3, "power", m_power, gRandom->Rndm() * 0.01, 0.05, 2.5);
247 minimizer->SetParameter(4, "ratio", m_ratio, gRandom->Rndm() * 0.01, 0.5, 2);
248 minimizer->FixParameter(2);
249 // minimizer->FixParameter(3);
250 minimizer->FixParameter(4);
251 counter++;
252 minimizer->ExecuteCommand("MIGRAD", arglist, nargs);
253 status = minimizer->ExecuteCommand("HESSE", arglist, nargs);
254 B2INFO("\t re-Fit iter: " << counter << ", status code: " << status);
255 }
256
257 if (status != 0) {
258 B2INFO("\t HadronSaturation::ERROR - BAD FIT!");
259 delete minimizer;
260 return;
261 }
262 }
263
264 double fitpar[5], fiterr[5];
265
266 for (int par = 0; par < 5; ++par) {
267 fitpar[par] = minimizer->GetParameter(par);
268 fiterr[par] = minimizer->GetParError(par);
269 }
270
271 B2INFO("\t Final Result: HadronSaturation: fit results");
272 B2INFO("\t alpha (" << m_alpha << "): " << fitpar[0] << " +- " << fiterr[0]);
273 B2INFO("\t gamma (" << m_gamma << "): " << fitpar[1] << " +- " << fiterr[1]);
274 B2INFO("\t delta (" << m_delta << "): " << fitpar[2] << " +- " << fiterr[2]);
275 B2INFO("\t power (" << m_power << "): " << fitpar[3] << " +- " << fiterr[3]);
276 B2INFO("\t ratio (" << m_ratio << "): " << fitpar[4] << " +- " << fiterr[4]);
277
278 // function value, estimated distance to minimum, errdef
279 // variable parameters, number of parameters
280 double chi2, edm, errdef;
281 int nvpar, nparx;
282 minimizer->GetStats(chi2, edm, errdef, nvpar, nparx);
283 B2INFO("\t\t Fit chi^2: " << chi2);
284
285 // write result file
286 std::ofstream parfile("sat-pars.fit.txt");
287 if (parfile.is_open()) {
288 for (int i = 0; i < 5; ++i) parfile << fitpar[i] << std::endl;
289 parfile.close();
290 } else {
291 B2WARNING("Unable to open sat-pars.fit.txt for writing fit parameters");
292 }
293
294 delete minimizer;
295}
Belle2::HadronSaturation::minuitFunction
static void minuitFunction(int &, double *, double &result, double *para, int)
functions for the hadron saturation correction
Definition HadronSaturation.cc:169

◆ minuitFunction()

void minuitFunction ( int & ,
double * ,
double & result,
double * para,
int  )
static

functions for the hadron saturation correction

Definition at line 169 of file HadronSaturation.cc.

170{
171 result = HC_obj->myFunction(par[0], par[1], par[2], par[3], par[4]);
172}
Belle2::HadronSaturation::myFunction
double myFunction(double alpha, double gamma, double delta, double power, double ratio)
some helper functions for the hadron saturation correction
Definition HadronSaturation.cc:98

◆ myFunction()

double myFunction ( double alpha,
double gamma,
double delta,
double power,
double ratio )

some helper functions for the hadron saturation correction

Definition at line 98 of file HadronSaturation.cc.

100{
101 unsigned int nevt = m_dedx.size();
102 if (nevt == 0) {
103 B2ERROR("HadronSaturation::myFunction called with no events!");
104 return 1e9; // large chi2 to signal bad fit
105 }
106
107 if (m_cosbins <= 0) {
108 B2ERROR("HadronSaturation::myFunction invalid cosbins = " << m_cosbins);
109 return 1e9;
110 }
111
112 double chisq = 0;
113 double dedxsum = 0;
114 std::vector<double> vdedxavg;
115
116 CDCDedxHadSat hadsat;
117
118 // --- compute averages ---
119 for (unsigned int i = 0; i < nevt; i++) {
120 double dedxcor = hadsat.D2I(
121 m_costheta[i],
122 hadsat.I2D(m_costheta[i], 1.0, alpha, gamma, delta, power, ratio) * m_dedx[i],
123 alpha, gamma, delta, power, ratio);
124
125 if (!std::isfinite(dedxcor)) {
126 B2WARNING("Non-finite dedxcor at event " << i);
127 continue;
128 }
129
130 dedxsum += dedxcor;
131
132 if ((i + 1) % m_cosbins == 0) {
133 vdedxavg.push_back(dedxsum / m_cosbins);
134 B2INFO("\t --> average: = " << dedxsum / m_cosbins << ", " << m_cosbins << "");
135 dedxsum = 0;
136 }
137 }
138
139 if (vdedxavg.empty()) {
140 B2ERROR("HadronSaturation::myFunction failed to compute averages!");
141 return 1e9;
142 }
143
144 // --- compute chi2 ---
145 for (unsigned int i = 0; i < nevt; i++) {
146 double dedxcor = hadsat.D2I(
147 m_costheta[i],
148 hadsat.I2D(m_costheta[i], 1.0, alpha, gamma, delta, power, ratio) * m_dedx[i],
149 alpha, gamma, delta, power, ratio);
150
151 if (!std::isfinite(dedxcor) || m_dedxerror[i] <= 0) continue;
152
153 int j = i / m_cosbins;
154 if (j >= int(vdedxavg.size())) {
155 B2WARNING("Index " << j << " out of range for vdedxavg!");
156 continue;
157 }
158
159 chisq += pow((dedxcor - vdedxavg[j]) / m_dedxerror[i], 2);
160 B2INFO("\t " << i << ") " << dedxcor << "/" << vdedxavg[j] << ", error was "
161 << m_dedxerror[i] << " De = " << hadsat.I2D(m_costheta[i], 1.0, alpha, gamma, delta, power, ratio) <<
162 ": Final " << chisq);
163 }
164
165 return chisq;
166}
Belle2::CDCDedxHadSat::I2D
double I2D(double cosTheta, double I) const
hadron saturation parameterization part 2
Definition CDCDedxHadSat.cc:71
Belle2::CDCDedxHadSat::D2I
double D2I(double cosTheta, double D) const
hadron saturation parameterization part 1
Definition CDCDedxHadSat.cc:48

◆ printEvents()

void printEvents ( int firstevent = 0,
int nevents = 50 )

print a sample of events

Definition at line 84 of file HadronSaturation.cc.

85{
86
87 if ((nevents + firstevent) > int(m_dedx.size())) nevents = m_dedx.size();
88
89 if (firstevent < 0 || (nevents + firstevent) > int(m_dedx.size()))
90 B2FATAL("HadronSaturation: trying to print events out of range (" << m_dedx.size() << ")");
91
92 for (int i = firstevent; i < nevents; ++i)
93 B2INFO("\t Event " << i << ":\t bg = " << m_betagamma[i] << "\t cos(theta) = " << m_costheta[i] << "\t dE/dx mean = " <<
94 m_dedx[i] << "\t dE/dx error = " << m_dedxerror[i]);
95
96}

◆ setCosBins()

void setCosBins ( int nbins)
inline

set the number of cosine bins

Definition at line 76 of file HadronSaturation.h.

76{ m_cosbins = nbins; }

Member Data Documentation

◆ m_alpha

double m_alpha
private

the alpha parameter for the hadron saturation correction

Definition at line 92 of file HadronSaturation.h.

◆ m_betagamma

std::vector< double > m_betagamma
private

a vector to hold beta-gamma values

Definition at line 100 of file HadronSaturation.h.

◆ m_cosbins

int m_cosbins
private

the number of cosine bins

Definition at line 90 of file HadronSaturation.h.

◆ m_costheta

std::vector< double > m_costheta
private

a vector to hold cos(theta) values

Definition at line 101 of file HadronSaturation.h.

◆ m_dedx

std::vector< double > m_dedx
private

a vector to hold dE/dx measurements

Definition at line 98 of file HadronSaturation.h.

◆ m_dedxerror

std::vector< double > m_dedxerror
private

a vector to hold dE/dx errors

Definition at line 99 of file HadronSaturation.h.

◆ m_delta

double m_delta
private

the delta parameter for the hadron saturation correction

Definition at line 94 of file HadronSaturation.h.

◆ m_gamma

double m_gamma
private

the gamma parameter for the hadron saturation correction

Definition at line 93 of file HadronSaturation.h.

◆ m_power

double m_power
private

the power parameter for the hadron saturation correction

Definition at line 95 of file HadronSaturation.h.

◆ m_ratio

double m_ratio
private

the ratio parameter for the hadron saturation correction

Definition at line 96 of file HadronSaturation.h.

The documentation for this class was generated from the following files: