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(); };
Belle2::HadronSaturation::clear
void clear()
clear the vectors
Definition: HadronSaturation.cc:257

Member Function Documentation

◆ clear()

void clear ( )

clear the vectors

Definition at line 257 of file HadronSaturation.cc.

258{
259
260 m_dedx.clear();
261 m_dedxerror.clear();
262 m_betagamma.clear();
263 m_costheta.clear();
264}

◆ 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->GetEvent(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}

◆ fitSaturation()

void fitSaturation ( )

perform the hadron saturation fit

Definition at line 148 of file HadronSaturation.cc.

149{
150
151 B2INFO("\t Performing the hadron saturation fit...");
152
153 // Construct the fitter
154 TFitter* minimizer = new TFitter(5);
155 minimizer->SetFCN(HadronSaturation::minuitFunction);
156
157 minimizer->SetParameter(0, "alpha", m_alpha, gRandom->Rndm() * 0.001, -5.0, 5.0);
158 minimizer->SetParameter(1, "gamma", m_gamma, gRandom->Rndm() * 0.001, -5.0, 5.0);
159 minimizer->SetParameter(2, "delta", m_delta, gRandom->Rndm() * 0.001, 0, 0.50);
160 minimizer->SetParameter(3, "power", m_power, gRandom->Rndm() * 0.01, 0.05, 2.5);
161 minimizer->SetParameter(4, "ratio", m_ratio, gRandom->Rndm() * 0.01, 0.5, 2);
162 minimizer->FixParameter(2);
163 minimizer->FixParameter(4);
164
165 // Set minuit fitting strategy
166 double arg[10];
167
168 double strategy(2.);
169 minimizer->ExecuteCommand("SET START", &strategy, 1);
170
171 double up(1.);
172 minimizer->ExecuteCommand("SET ERR", &up, 1);
173
174 // Suppress much of the output of MINUIT
175 arg[0] = -1;
176 minimizer->ExecuteCommand("SET PRINT", arg, 1);
177
178 double eps_machine(std::numeric_limits<double>::epsilon());
179 minimizer->ExecuteCommand("SET EPS", &eps_machine, 1);
180
181 double fitpar[5], fiterr[5];
182
183 for (int i = 0; i < 30; ++i) {
184
185 minimizer->SetParameter(0, "alpha", m_alpha, gRandom->Rndm() * 0.001, -5.0, 5.0);
186 minimizer->SetParameter(1, "gamma", m_gamma, gRandom->Rndm() * 0.001, -5.0, 5.0);
187 minimizer->SetParameter(2, "delta", m_delta, gRandom->Rndm() * 0.001, 0, 0.50);
188 minimizer->SetParameter(3, "power", m_power, gRandom->Rndm() * 0.01, 0.05, 2.5);
189 minimizer->SetParameter(4, "ratio", m_ratio, gRandom->Rndm() * 0.01, 0.5, 2);
190 minimizer->FixParameter(2);
191 minimizer->FixParameter(4);
192
193 double maxcalls(5000.), tolerance(0.1);
194 double arglist[] = {maxcalls, tolerance};
195 unsigned int nargs(2);
196 minimizer->ExecuteCommand("MIGRAD", arglist, nargs);
197 minimizer->ExecuteCommand("MIGRAD", arglist, nargs);
198 int status = minimizer->ExecuteCommand("HESSE", arglist, nargs);
199
200 minimizer->PrintResults(1, 0);
201 B2INFO("\t iter = " << i << ": Fit status: " << status);
202
203 int counter = 0;
204 while (status != 0 && counter < 5) {
205 minimizer->SetParameter(0, "alpha", m_alpha, gRandom->Rndm() * 0.001, -5.0, 5.0);
206 minimizer->SetParameter(1, "gamma", m_gamma, gRandom->Rndm() * 0.001, -5.0, 5.0);
207 minimizer->SetParameter(2, "delta", m_delta, gRandom->Rndm() * 0.001, 0, 0.50);
208 minimizer->SetParameter(3, "power", m_power, gRandom->Rndm() * 0.01, 0.05, 2.5);
209 minimizer->SetParameter(4, "ratio", m_ratio, gRandom->Rndm() * 0.01, 0.5, 2);
210 minimizer->FixParameter(2);
211 // minimizer->FixParameter(3);
212 minimizer->FixParameter(4);
213 counter++;
214 minimizer->ExecuteCommand("MIGRAD", arglist, nargs);
215 status = minimizer->ExecuteCommand("HESSE", arglist, nargs);
216 B2INFO("\t re-Fit iter: " << counter << ", status code: " << status);
217 }
218
219 if (status != 0) {
220 B2INFO("\t HadronSaturation::ERROR - BAD FIT!");
221 return;
222 }
223 }
224
225 for (int par = 0; par < 5; ++par) {
226 fitpar[par] = minimizer->GetParameter(par);
227 fiterr[par] = minimizer->GetParError(par);
228 }
229
230 B2INFO("\t Final Result: HadronSaturation: fit results");
231 B2INFO("\t alpha (" << m_alpha << "): " << fitpar[0] << " +- " << fiterr[0]);
232 B2INFO("\t gamma (" << m_gamma << "): " << fitpar[1] << " +- " << fiterr[1]);
233 B2INFO("\t delta (" << m_delta << "): " << fitpar[2] << " +- " << fiterr[2]);
234 B2INFO("\t power (" << m_power << "): " << fitpar[3] << " +- " << fiterr[3]);
235 B2INFO("\t ratio (" << m_ratio << "): " << fitpar[4] << " +- " << fiterr[4]);
236
237 // function value, estimated distance to minimum, errdef
238 // variable parameters, number of parameters
239 double chi2, edm, errdef;
240 int nvpar, nparx;
241 minimizer->GetStats(chi2, edm, errdef, nvpar, nparx);
242 B2INFO("\t\t Fit chi^2: " << chi2);
243
244 std::ofstream parfile;
245 parfile.open("sat-pars.fit.txt");
246
247 parfile << fitpar[0] << std::endl;
248 parfile << fitpar[1] << std::endl;
249 parfile << fitpar[2] << std::endl;
250 parfile << fitpar[3] << std::endl;
251 parfile << fitpar[4] << std::endl;
252 parfile.close();
253
254 delete minimizer;
255}
Belle2::HadronSaturation::minuitFunction
static void minuitFunction(int &, double *, double &result, double *para, int)
functions for the hadron saturation correction
Definition: HadronSaturation.cc:142

◆ minuitFunction()

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

functions for the hadron saturation correction

Definition at line 142 of file HadronSaturation.cc.

143{
144 result = HC_obj->myFunction(par[0], par[1], par[2], par[3], par[4]);
145}
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.

99{
100
101 unsigned int nevt = m_dedx.size();
102 double chisq = 0, dedxsum = 0;
103 std::vector< double > vdedxavg;
104
105 // Compute the average value (across cos(theta)) for each bin of beta-gamma.
106 // NOTE: the correction is not constrained to a certain value (1), it
107 // changes as a function of beta gamma...
108 CDCDedxHadSat hadsat;
109
110 double dedxcor;
111 for (unsigned int i = 0; i < nevt; i++) {
112
113 dedxcor = hadsat.D2I(m_costheta[i], hadsat.I2D(m_costheta[i], 1.0, alpha, gamma, delta, power, ratio) * m_dedx[i], alpha, gamma,
114 delta, power, ratio);
115
116 dedxsum += dedxcor;
117
118 if ((i + 1) % m_cosbins == 0) {
119 vdedxavg.push_back(dedxsum / m_cosbins);
120 B2INFO("\t --> average: = " << dedxsum / m_cosbins << ", " << m_cosbins << "");
121 dedxsum = 0;
122 }
123 }
124
125 // Construct a chi^2 value for the difference between the average in a cosine bin
126 // to the actual values
127 for (unsigned int i = 0; i < nevt; i++) {
128 dedxcor = hadsat.D2I(m_costheta[i], hadsat.I2D(m_costheta[i], 1.0, alpha, gamma, delta, power, ratio) * m_dedx[i], alpha, gamma,
129 delta, power, ratio);
130
131 int j = (int)i / m_cosbins;
132 chisq += pow((dedxcor - vdedxavg[j]) / m_dedxerror[i], 2);
133 B2INFO("\t " << i << ") " << dedxcor << "/" << vdedxavg[j] << ", error was "
134 << m_dedxerror[i] << " De = " << hadsat.I2D(m_costheta[i], 1.0, alpha, gamma, delta, power, ratio) <<
135 ": Final " << chisq);
136 }
137
138 return chisq;
139}
Belle2::CDCDedxHadSat
Class to hold the hadron saturation functions.
Definition: CDCDedxHadSat.h:31
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: