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

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

◆ 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 = static_cast<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:299

◆ fitSaturation()

void fitSaturation ( )

perform the hadron saturation fit

Definition at line 177 of file HadronSaturation.cc.

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

◆ minuitFunction()

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

functions for the hadron saturation correction

Definition at line 171 of file HadronSaturation.cc.

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