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 15 of file HadronSaturation.cc.

15 :
16 m_cosbins(8),
17 m_alpha(0.044286100),
18 m_gamma(0.0019308200),
19 m_delta(0.020),
20 m_power(1.3205300),
21 m_ratio(1.0000)
22{
23 m_dedx.clear();
24 m_dedxerror.clear();
25 m_betagamma.clear();
26 m_costheta.clear();
27 HC_obj = this;
28}
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 30 of file HadronSaturation.cc.

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

◆ ~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 298 of file HadronSaturation.cc.

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

◆ fillSample()

void fillSample ( TString infilename)

fill the vectors below

Definition at line 45 of file HadronSaturation.cc.

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

◆ fitSaturation()

void fitSaturation ( )

perform the hadron saturation fit

Definition at line 176 of file HadronSaturation.cc.

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

◆ minuitFunction()

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

functions for the hadron saturation correction

Definition at line 170 of file HadronSaturation.cc.

171{
172 result = HC_obj->myFunction(par[0], par[1], par[2], par[3], par[4]);
173}
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:99

◆ myFunction()

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

some helper functions for the hadron saturation correction

Definition at line 99 of file HadronSaturation.cc.

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

◆ printEvents()

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

print a sample of events

Definition at line 85 of file HadronSaturation.cc.

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

◆ 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: