BOSS 7.1.1
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | List of all members
HadronSaturation Class Reference

#include <HadronSaturation.h>

Public Member Functions

 HadronSaturation ()
 
 HadronSaturation (double alpha, double gamma, double delta, double power, double ratio)
 
virtual ~HadronSaturation ()
 
void setParameters (double par[])
 
void setParameters (std::string parfile)
 
void fillSample (TString infilename)
 
void printEvents (int firstevent, int nevents)
 
void fitSaturation ()
 
void clear ()
 
void setCosBins (int nbins)
 
void setFlag (int flag)
 
double myFunction (double alpha, double gamma, double delta, double power, double ratio)
 
double D2I (double cosTheta, double D, double alpha, double gamma, double delta, double power, double ratio) const
 
double I2D (double cosTheta, double I, double alpha, double gamma, double delta, double power, double ratio) const
 
double D2I (double cosTheta, double D=1) const
 
double I2D (double cosTheta, double I=1) const
 

Static Public Member Functions

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

Detailed Description

Definition at line 34 of file HadronSaturation.h.

Constructor & Destructor Documentation

◆ HadronSaturation() [1/2]

HadronSaturation::HadronSaturation ( )

Definition at line 5 of file HadronSaturation.cc.

5 :
6 m_flag(0),
7 m_cosbins(18),
8 m_alpha(0.0266879),
9 m_gamma(0.00378704),
10 m_delta(0.0038601),
11 m_power(1.8957),
12 m_ratio(3.09834)
13{
14 m_dedx.clear();
15 m_dedxerror.clear();
16 m_betagamma.clear();
17 m_costheta.clear();
18 HC_obj = this;
19}

◆ HadronSaturation() [2/2]

HadronSaturation::HadronSaturation ( double alpha,
double gamma,
double delta,
double power,
double ratio )

Definition at line 21 of file HadronSaturation.cc.

21 :
22 m_flag(0),
23 m_cosbins(18),
24 m_alpha(alpha),
25 m_gamma(gamma),
26 m_delta(delta),
27 m_power(power),
28 m_ratio(ratio)
29{
30 m_dedx.clear();
31 m_dedxerror.clear();
32 m_betagamma.clear();
33 m_costheta.clear();
34}
delta
const double delta
Definition CheckRung/check/gainVsDoca.cxx:32
alpha
const double alpha
Definition FastVertexFit.cxx:4

◆ ~HadronSaturation()

virtual HadronSaturation::~HadronSaturation ( )
inlinevirtual

Definition at line 40 of file HadronSaturation.h.

40{ clear(); };

Member Function Documentation

◆ clear()

void HadronSaturation::clear ( )

Definition at line 299 of file HadronSaturation.cc.

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

Referenced by fillSample(), and ~HadronSaturation().

◆ D2I() [1/2]

double HadronSaturation::D2I ( double cosTheta,
double D,
double alpha,
double gamma,
double delta,
double power,
double ratio ) const

Definition at line 119 of file HadronSaturation.cc.

119 {
120
121 double absCosTheta = fabs(cosTheta);
122 double projection = pow(absCosTheta, power) + delta;
123 double chargeDensity = D / projection;
124 double numerator = 1 + alpha * chargeDensity;
125 double denominator = 1 + gamma * chargeDensity;
126
127 double I = D * ratio * numerator / denominator;
128
129 return I;
130}
I
const DifComplex I

Referenced by HadronPrep::bgCosThetaFits(), HadronPrep::bgFits(), HadronCalibration::fitSigmaVsNHit(), HadronCalibration::fitSigmaVsSin(), ElectronCorrection::HadronCorrection(), myFunction(), and HadronCalibration::plotEfficiency().

◆ D2I() [2/2]

double HadronSaturation::D2I ( double cosTheta,
double D = 1 ) const
inline

Definition at line 74 of file HadronSaturation.h.

74 {
75
76 double absCosTheta = fabs(cosTheta);
77 double projection = pow(absCosTheta, m_power) + m_delta;
78 double chargeDensity = D / projection;
79 double numerator = 1 + m_alpha * chargeDensity;
80 double denominator = 1 + m_gamma * chargeDensity;
81
82 double I = D * m_ratio * numerator / denominator;
83
84 return I;
85 }

◆ fillSample()

void HadronSaturation::fillSample ( TString infilename)

Definition at line 61 of file HadronSaturation.cc.

61 {
62
63 std::cout << "HadronSaturation: filling sample events..." << std::endl;
64
65 // clear the vectors to be filled
66 clear();
67
68 std::vector< TString > types;
69 types.push_back("pion");
70 types.push_back("kaon");
71 types.push_back("proton");
72 types.push_back("muon");
73 types.push_back("electron");
74
75 // fill the containers with a previously prepared sample
76 TFile* satFile = new TFile(infilename);
77
78 for( int i = 0; i < 4; ++i ){
79 TTree* satTree = (TTree*)satFile->Get(types[i]);
80
81 double satdedx, saterror, satbg, satcosth;
82 satTree->SetBranchAddress("dedx",&satdedx);
83 satTree->SetBranchAddress("error",&saterror);
84 satTree->SetBranchAddress("bg",&satbg);
85 satTree->SetBranchAddress("costh",&satcosth);
86
87 // fill the vectors
88 for( unsigned int i = 0; i < satTree->GetEntries(); ++i ){
89 satTree->GetEvent(i);
90
91 if( saterror == 0 ) continue;
92
93 m_dedx.push_back(satdedx);
94 m_dedxerror.push_back(saterror);
95 m_betagamma.push_back(satbg);
96 m_costheta.push_back(satcosth);
97 }
98 }
99
100 satFile->Close();
101}

Referenced by HadronInterface::SaturationCorrection().

◆ fitSaturation()

void HadronSaturation::fitSaturation ( )

Definition at line 205 of file HadronSaturation.cc.

205 {
206
207 std::cout << "Performing the hadron saturation fit..." << std::endl;
208
209 // Construct the fitter
210 TFitter* minimizer = new TFitter(5);
211 minimizer->SetFCN(HadronSaturation::minuitFunction);
212
213 TRandom* rand = new TRandom();
214
215 minimizer->SetParameter(0,"alpha",m_alpha,rand->Rndm()*0.001,0.0,0.5);
216 minimizer->SetParameter(1,"gamma",m_gamma,rand->Rndm()*0.001,-2,2);
217 minimizer->SetParameter(2,"delta",m_delta,rand->Rndm()*0.001,-2,2);
218 minimizer->SetParameter(3,"power",m_power,rand->Rndm()*0.01,-5,5);
219 minimizer->SetParameter(4,"ratio",m_ratio,rand->Rndm()*0.01,-100,100);
220
221 // Set minuit fitting strategy
222 double arg[10];
223 arg[0] = 2;
224 minimizer->ExecuteCommand("SET STR",arg,1);
225
226 // Suppress much of the output of MINUIT
227 arg[0] = -1;
228 minimizer->ExecuteCommand("SET PRINT",arg,1);
229
230 // Minimize with MIGRAD
231 arg[0] = 10000;
232
233 double fitpar[5], fiterr[5];
234 for( int i = 0; i < 20; ++i ){
235
236 minimizer->SetParameter(0,"alpha",m_alpha,rand->Rndm()*0.001,0.0,0.5);
237 minimizer->SetParameter(1,"gamma",m_gamma,rand->Rndm()*0.001,-2,2);
238 minimizer->SetParameter(2,"delta",m_delta,rand->Rndm()*0.001,-2,2);
239 minimizer->SetParameter(3,"power",m_power,rand->Rndm()*0.01,-5,5);
240 minimizer->SetParameter(4,"ratio",m_ratio,rand->Rndm()*0.01,-100,100);
241
242 int status = minimizer->ExecuteCommand("MIGRAD",arg,1);
243 status = minimizer->ExecuteCommand("MIGRAD",arg,1);
244 minimizer->PrintResults(1,0);
245 std::cout << "Fit status: " << status << std::endl;
246
247 int counter = 0;
248 while( status != 0 && counter < 10 ){
249 counter++;
250
251 minimizer->SetParameter(0,"alpha",m_alpha,rand->Rndm()*0.001,0.0,0.5);
252 minimizer->SetParameter(1,"gamma",m_gamma,rand->Rndm()*0.001,-2,2);
253 minimizer->SetParameter(2,"delta",m_delta,rand->Rndm()*0.001,-2,2);
254 minimizer->SetParameter(3,"power",m_power,rand->Rndm()*0.01,-5,5);
255 minimizer->SetParameter(4,"ratio",m_ratio,rand->Rndm()*0.01,-100,100);
256
257 status = minimizer->ExecuteCommand("MIGRAD",arg,1);
258 std::cout << "Fit status: " << status << std::endl;
259 }
260 if( status != 0 ){
261 std::cout << "HadronSaturation::ERROR - BAD FIT!" << std::endl;
262 return;
263 }
264 }
265
266 for( int par = 0; par < 5; ++par ){
267 fitpar[par] = minimizer->GetParameter(par);
268 fiterr[par] = minimizer->GetParError(par);
269 }
270
271 std::cout << "HadronSaturation: fit results" << std::endl;
272 std::cout << "\t alpha (" << m_alpha << "): " << fitpar[0] << " +- " << fiterr[0] << std::endl;
273 std::cout << "\t gamma (" << m_gamma << "): " << fitpar[1] << " +- " << fiterr[1] << std::endl;
274 std::cout << "\t delta (" << m_delta << "): " << fitpar[2] << " +- " << fiterr[2] << std::endl;
275 std::cout << "\t power (" << m_power << "): " << fitpar[3] << " +- " << fiterr[3] << std::endl;
276 std::cout << "\t ratio (" << m_ratio << "): " << fitpar[4] << " +- " << fiterr[4] << std::endl;
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 std::cout << "\t\t Fit chi^2: " << chi2 << std::endl;
284
285 std::ofstream parfile;
286 parfile.open("sat-pars.fit.txt");
287
288 parfile << fitpar[0] << std::endl;
289 parfile << fitpar[1] << std::endl;
290 parfile << fitpar[2] << std::endl;
291 parfile << fitpar[3] << std::endl;
292 parfile << fitpar[4] << std::endl;
293 parfile.close();
294
295 delete minimizer;
296}
arg
double arg(const EvtComplex &c)
Definition EvtComplex.hh:227
HadronSaturation::minuitFunction
static void minuitFunction(int &nDim, double *gout, double &result, double *para, int flg)
Definition HadronSaturation.cc:200

Referenced by HadronInterface::SaturationCorrection().

◆ I2D() [1/2]

double HadronSaturation::I2D ( double cosTheta,
double I,
double alpha,
double gamma,
double delta,
double power,
double ratio ) const

Definition at line 133 of file HadronSaturation.cc.

133 {
134
135 double absCosTheta = fabs(cosTheta);
136 double projection = pow(absCosTheta, power) + delta;
137
138 double a = alpha/projection;
139 double b = 1 - gamma/projection * (I/ratio);
140 double c = -1.0 * I/ratio;
141
142 if( b == 0 && a == 0 ){
143 std::cout << "both a and b coefficiants for hadron correction are 0" << std::endl;
144 return I;
145 }
146
147 double D = (a != 0) ? (-b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a) : -c / b;
148 if( D < 0 ){
149 // std::cout << "D is less 0! will try another solution" << std::endl;
150 D = (a != 0) ? (-b - sqrt(b * b + 4.0 * a * c)) / (2.0 * a) : -c / b;
151 if( D < 0 ){
152 // std::cout << "D is still less 0! just return uncorrectecd value" << std::endl;
153 return I;
154 }
155 }
156
157 return D;
158}
b
const double b
Definition slope.cxx:9

Referenced by HadronPrep::bgCosThetaFits(), HadronPrep::bgFits(), HadronCalibration::fitSigmaVsNHit(), HadronCalibration::fitSigmaVsSin(), WidgetGenerator::generateEvents(), ElectronCorrection::HadronCorrection(), myFunction(), HadronCalibration::plotEfficiency(), and WidgetGenerator::simulateDedx().

◆ I2D() [2/2]

double HadronSaturation::I2D ( double cosTheta,
double I = 1 ) const
inline

Definition at line 87 of file HadronSaturation.h.

87 {
88
89 double absCosTheta = fabs(cosTheta);
90 double projection = pow(absCosTheta, m_power) + m_delta;
91
92 double a = m_alpha/projection;
93 double b = 1 - m_gamma/projection * (I/m_ratio);
94 double c = -1.0 * I/m_ratio;
95
96 if( b == 0 && a == 0 ){
97 std::cout << "both a and b coefficiants for hadron correction are 0" << std::endl;
98 return I;
99 }
100
101 double D = (a != 0) ? (-b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a) : -c / b;
102 if( D < 0 ){
103 std::cout << "D is less 0! will try another solution" << std::endl;
104 D = (a != 0) ? (-b - sqrt(b * b + 4.0 * a * c)) / (2.0 * a) : -c / b;
105 if( D < 0 ){
106 std::cout << "D is still less 0! just return uncorrectecd value" << std::endl;
107 return I;
108 }
109 }
110
111 return D;
112 }

◆ minuitFunction()

void HadronSaturation::minuitFunction ( int & nDim,
double * gout,
double & result,
double * para,
int flg )
static

Definition at line 200 of file HadronSaturation.cc.

200 {
201 result = HC_obj->myFunction(par[0],par[1],par[2],par[3],par[4]);
202}
HadronSaturation::myFunction
double myFunction(double alpha, double gamma, double delta, double power, double ratio)
Definition HadronSaturation.cc:161

Referenced by fitSaturation().

◆ myFunction()

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

Definition at line 161 of file HadronSaturation.cc.

161 {
162
163 unsigned int nevt = m_dedx.size();
164 double chisq = 0, dedxsum = 0;
165 std::vector< double > vdedxavg;
166
167 // Compute the average value (across cos(theta)) for each bin of beta-gamma.
168 // NOTE: the correction is not constrained to a certain value (1), it
169 // changes as a function of beta gamma...
170 double dedxcor = 0;
171 for( unsigned int i = 0; i < nevt; i++ ){
172 if( m_flag == 0 )
173 dedxcor = D2I(m_costheta[i],I2D(m_costheta[i],1.0,alpha,gamma,delta,power,ratio)*m_dedx[i],alpha,gamma,delta,power,ratio);
174 else
175 dedxcor = D2I(m_costheta[i],m_dedx[i],alpha,gamma,delta,power,ratio);
176 dedxsum += dedxcor;
177
178 if( (i+1)%m_cosbins == 0 ){
179 vdedxavg.push_back(dedxsum/m_cosbins);
180 dedxsum = 0;
181 }
182 }
183
184 // Construct a chi^2 value for the difference between the average in a cosine bin
185 // to the actual values
186 for( unsigned int i = 0; i < nevt; i++ ){
187 if( m_flag == 0 )
188 dedxcor = D2I(m_costheta[i],I2D(m_costheta[i],1.0,alpha,gamma,delta,power,ratio)*m_dedx[i],alpha,gamma,delta,power,ratio);
189 else
190 dedxcor = D2I(m_costheta[i],m_dedx[i],alpha,gamma,delta,power,ratio);
191
192 int j = (int)i/m_cosbins;
193 chisq += pow((dedxcor-vdedxavg[j])/m_dedxerror[i],2);
194 }
195
196 return chisq;
197}
nevt
Int_t nevt[10]
Definition DataBase/tau_mode.c:56
HadronSaturation::I2D
double I2D(double cosTheta, double I, double alpha, double gamma, double delta, double power, double ratio) const
Definition HadronSaturation.cc:133
HadronSaturation::D2I
double D2I(double cosTheta, double D, double alpha, double gamma, double delta, double power, double ratio) const
Definition HadronSaturation.cc:119

Referenced by minuitFunction().

◆ printEvents()

void HadronSaturation::printEvents ( int firstevent = 0,
int nevents = 10 )

Definition at line 104 of file HadronSaturation.cc.

104 {
105
106 if( firstevent < 0 || (nevents+firstevent) > m_dedx.size() ){
107 std::cout << "HadronSaturation: trying to print events out of range (" << m_dedx.size() << ")" << std::endl;
108 exit(1);
109 }
110
111 std::cout << std::endl << std::endl;
112 for( int i = firstevent; i < nevents; ++i ){
113 std::cout << "Event " << i << ":\t bg = " << m_betagamma[i] << "\t cos(theta) = " << m_costheta[i] << "\t dE/dx mean = " << m_dedx[i] << "\t dE/dx error = " << m_dedxerror[i] << std::endl;
114 }
115 std::cout << std::endl;
116}

Referenced by HadronInterface::SaturationCorrection().

◆ setCosBins()

void HadronSaturation::setCosBins ( int nbins)
inline

Definition at line 61 of file HadronSaturation.h.

61{ m_cosbins = nbins; }
nbins
const int nbins
Definition CheckRung/check/dEdxVsCos.cxx:18

◆ setFlag()

void HadronSaturation::setFlag ( int flag)
inline

Definition at line 64 of file HadronSaturation.h.

64{ m_flag = flag; }

◆ setParameters() [1/2]

void HadronSaturation::setParameters ( double par[])

Definition at line 37 of file HadronSaturation.cc.

37 {
38 m_alpha = par[0];
39 m_gamma = par[1];
40 m_delta = par[2];
41 m_power = par[3];
42 m_ratio = par[4];
43}

Referenced by HadronPrep::bgCosThetaFits(), HadronPrep::bgFits(), HadronCalibration::fitSigmaVsNHit(), HadronCalibration::fitSigmaVsSin(), ElectronCorrection::initializeParameters(), HadronCalibration::plotEfficiency(), HadronInterface::SaturationCorrection(), and setParameters().

◆ setParameters() [2/2]

void HadronSaturation::setParameters ( std::string parfile)

Definition at line 46 of file HadronSaturation.cc.

46 {
47 double hadpar[5];
48 std::ifstream parfile(infile.c_str());
49 if( !parfile.fail() ){
50 for( int i = 0; i < 5; ++i ){
51 parfile >> hadpar[i];
52 }
53 parfile.close();
54 setParameters(hadpar);
55 }
56 else
57 std::cout << "HadronSaturation: WARNING - parameter file " << infile << " does not exist, using defaults..." << std::endl;
58}
HadronSaturation::setParameters
void setParameters(double par[])
Definition HadronSaturation.cc:37
The documentation for this class was generated from the following files: