BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
DedxCalibParameters.cxx File Reference
#include "TGraphErrors.h"
#include <TMath.h>
#include <cmath>
#include <iostream>
#include <vector>
#include "DedxCalibAlg/DedxCalibParameters.h"

Go to the source code of this file.

Functions

float vavset_ (double *, double *, int *)
float vavden_ (double *)
float prob_ (float *, int *)
double mylan (double *x, double *par)
double landaun (double *x, double *par)
double Landau (double *x, double *par)
double Vavilov (double *x, double *par)
double AsymGauss (double *x, double *par)
void dedx_pid_exp_old (int landau, int runflag, float dedx, int Nohit, float mom, float theta, float t0, float lsamp, double dedx_exp[5], double ex_sigma[5], double pid_prob[5], double chi_dedx[5])
void dedx_pid_exp (int vflag[3], float dedx, int trkalg, int Nohit, float mom, float theta, float t0, float lsamp, double dedx_exp[5], double ex_sigma[5], double pid_prob[5], double chi_dedx[5], std::vector< double > &par, std::vector< double > &sig_par)
double SpaceChargeCorrec (double m_theta, double mom, int Particle, double dEdx)

Function Documentation

◆ AsymGauss()

double AsymGauss ( double * x,
double * par )

Definition at line 43 of file DedxCalibParameters.cxx.

43 {
44 double delta, a, b, c, d, fitval;
45 a = TMath::Sqrt( log( 4.0 ) );
46 b = ( x[0] - par[1] ) / par[2];
47 c = TMath::Power( par[3], 2.0 );
48 delta = ( 1 + TMath::SinH( par[3] * a ) / a ) * b;
49 d = TMath::Power( log( delta ), 2.0 ) / c + c;
50 fitval = par[0] * TMath::Exp( -0.5 * d );
51 return fitval;
52}

Referenced by DedxCalibDocaEAng::AnalyseHists(), DedxCalibEAng::AnalyseHists(), DedxCalibLayerGain::AnalyseHists(), and DedxCalibWireGain::AnalyseHists().

◆ dedx_pid_exp()

void dedx_pid_exp ( int vflag[3],
float dedx,
int trkalg,
int Nohit,
float mom,
float theta,
float t0,
float lsamp,
double dedx_exp[5],
double ex_sigma[5],
double pid_prob[5],
double chi_dedx[5],
std::vector< double > & par,
std::vector< double > & sig_par )

Definition at line 162 of file DedxCalibParameters.cxx.

165 {
166 const int par_cand( 5 );
167 const float Charge_Mass[par_cand] = { 0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
168 double beta_G, beta, betterm, bethe_B;
169
170 int dedxflag = vflag[0];
171 int sigmaflag = vflag[1];
172 bool ifMC = false;
173 if ( vflag[2] == 1 ) ifMC = true;
174
175 int Nmax_prob( 0 );
176 float max_prob( -0.01 );
177 int ndf;
178 float chi2;
179
180 for ( int it = 0; it < par_cand; it++ )
181 {
182 beta_G = mom / Charge_Mass[it];
183
184 if ( dedxflag == 1 )
185 {
186 beta = beta_G / sqrt( 1 + ( beta_G ) * ( beta_G ) );
187 betterm = par[1] - log( par[2] + pow( 1 / beta_G, par[4] ) );
188 bethe_B = par[0] / pow( beta, par[3] ) * betterm - par[0];
189 }
190 else if ( dedxflag == 2 )
191 {
192 double A = 0, B = 0, C = 0;
193 double x = beta_G;
194 if ( x < 4.5 ) A = 1;
195 else if ( x < 10 ) B = 1;
196 else C = 1;
197 double partA = par[0] * pow( sqrt( x * x + 1 ), par[2] ) / pow( x, par[2] ) *
198 ( par[1] - par[5] * log( pow( 1 / x, par[3] ) ) ) -
199 par[4] + exp( par[6] + par[7] * x );
200 double partB = par[8] * pow( x, 3 ) + par[9] * pow( x, 2 ) + par[10] * x + par[11];
201 double partC = -par[12] * log( par[15] + pow( 1 / x, par[13] ) ) + par[14];
202 bethe_B = 550 * ( A * partA + B * partB + C * partC );
203 // for fermi plateau ( the electron region) we just use 1.0
204 if ( beta_G > 100 ) bethe_B = 550 * 1.0;
205 }
206
207 if ( ifMC )
208 {
209 double A = 0, B = 0, C = 0;
210 double x = beta_G;
211 if ( x < 4.5 ) A = 1;
212 else if ( x < 10 ) B = 1;
213 else C = 1;
214 double partA = par[0] * pow( sqrt( x * x + 1 ), par[2] ) / pow( x, par[2] ) *
215 ( par[1] - par[5] * log( pow( 1 / x, par[3] ) ) ) -
216 par[4] + exp( par[6] + par[7] * x );
217 double partB = par[8] * pow( x, 3 ) + par[9] * pow( x, 2 ) + par[10] * x + par[11];
218 double partC = -par[12] * log( par[15] + pow( 1 / x, par[13] ) ) + par[14];
219 bethe_B = 550 * ( A * partA + B * partB + C * partC );
220 // for fermi plateau ( the electron region) we just use 1.0
221 if ( beta_G > 100 ) bethe_B = 550 * 1.0;
222 }
223
224 if ( Nohit > 0 )
225 {
226 dedx_exp[it] = bethe_B;
227 double sig_the = std::sin( (double)theta );
228 double f_betagamma, g_sinth, h_nhit, i_t0;
229
230 if ( ifMC )
231 {
232
233 double x = beta_G;
234 double nhit = (double)Nohit;
235 double sigma_bg = 1.0;
236 if ( x > 0.3 )
237 sigma_bg = sig_par[0] * exp( sig_par[1] * x ) +
238 sig_par[2] * exp( sig_par[3] * pow( x, 0.25 ) ) + sig_par[4];
239 else sigma_bg = sig_par[5] * exp( sig_par[6] * x ) + sig_par[7];
240
241 double cor_nhit = 1.0;
242 if ( nhit < 5 ) nhit = 5;
243 if ( nhit <= 35 )
244 cor_nhit = sig_par[8] * pow( nhit, 4 ) + sig_par[9] * pow( nhit, 3 ) +
245 sig_par[10] * pow( nhit, 2 ) + sig_par[11] * nhit + sig_par[12];
246
247 double cor_sin = 1.0;
248 if ( sig_the > 0.99 ) sig_the = 0.99;
249 cor_sin = sig_par[13] * pow( sig_the, 4 ) + sig_par[14] * pow( sig_the, 3 ) +
250 sig_par[15] * pow( sig_the, 2 ) + sig_par[16] * sig_the + sig_par[17];
251
252 // sigma vs t0
253 double cor_t0 = 1.0;
254
255 // calculate sigma
256 if ( trkalg == 1 ) ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
257 else ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[18];
258 }
259 else
260 {
261 if ( sigmaflag == 1 )
262 {
263 f_betagamma = sig_par[0] * pow( beta_G, sig_par[1] ) + sig_par[2];
264 g_sinth = ( sig_par[3] * sig_the * sig_the + sig_par[4] ) /
265 ( sig_par[3] * sig_par[5] * sig_par[5] + sig_par[4] );
266 h_nhit =
267 ( sig_par[6] * Nohit * Nohit + sig_par[7] * Nohit + sig_par[8] ) /
268 ( sig_par[6] * sig_par[9] * sig_par[9] + sig_par[7] * sig_par[9] + sig_par[8] );
269 if ( sig_par[13] != 0 )
270 i_t0 = ( sig_par[10] * t0 * t0 + sig_par[11] * t0 + sig_par[12] ) /
271 ( sig_par[10] * sig_par[13] * sig_par[13] + sig_par[11] * sig_par[13] +
272 sig_par[12] );
273 else if ( sig_par[13] == 0 ) i_t0 = 1;
274 // cout<<"f_betagamma : "<<f_betagamma<<" g_sinth : "<<g_sinth<<" h_nhit
275 // : "
276 //<<h_nhit<<" i_t0 : "<<i_t0<<endl;
277 ex_sigma[it] = f_betagamma * g_sinth * h_nhit * i_t0;
278 }
279 else if ( sigmaflag == 2 )
280 {
281 double x = beta_G;
282 double nhit = (double)Nohit;
283 double sigma_bg = 1.0;
284 if ( x > 0.3 )
285 sigma_bg = sig_par[0] * exp( sig_par[1] * x ) +
286 sig_par[2] * exp( sig_par[3] * pow( x, 0.25 ) ) + sig_par[4];
287 else sigma_bg = sig_par[5] * exp( sig_par[6] * x ) + sig_par[7];
288
289 double cor_nhit = 1.0;
290 if ( nhit < 5 ) nhit = 5;
291 if ( nhit <= 35 )
292 cor_nhit = sig_par[8] * pow( nhit, 4 ) + sig_par[9] * pow( nhit, 3 ) +
293 sig_par[10] * pow( nhit, 2 ) + sig_par[11] * nhit + sig_par[12];
294
295 double cor_sin = 1.0;
296 if ( sig_the > 0.99 ) sig_the = 0.99;
297 cor_sin = sig_par[13] * pow( sig_the, 4 ) + sig_par[14] * pow( sig_the, 3 ) +
298 sig_par[15] * pow( sig_the, 2 ) + sig_par[16] * sig_the + sig_par[17];
299
300 double cor_t0 = 1;
301 if ( t0 > 1200 ) t0 = 1200;
302 if ( t0 > 800 ) cor_t0 = sig_par[18] * pow( t0, 2 ) + sig_par[19] * t0 + sig_par[20];
303
304 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
305 }
306 else if ( sigmaflag == 3 )
307 {
308 double x = beta_G;
309 double nhit = (double)Nohit;
310 double sigma_bg = 1.0;
311 if ( x > 0.3 )
312 sigma_bg = sig_par[0] * exp( sig_par[1] * x ) +
313 sig_par[2] * exp( sig_par[3] * pow( x, 0.25 ) ) + sig_par[4];
314 else sigma_bg = sig_par[5] * exp( sig_par[6] * x ) + sig_par[7];
315
316 double cor_nhit = 1.0;
317 if ( nhit < 5 ) nhit = 5;
318 if ( nhit <= 35 )
319 cor_nhit = sig_par[8] * pow( nhit, 4 ) + sig_par[9] * pow( nhit, 3 ) +
320 sig_par[10] * pow( nhit, 2 ) + sig_par[11] * nhit + sig_par[12];
321
322 double cor_sin = 1.0;
323 if ( sig_the > 0.99 ) sig_the = 0.99;
324 cor_sin = sig_par[13] * pow( sig_the, 4 ) + sig_par[14] * pow( sig_the, 3 ) +
325 sig_par[15] * pow( sig_the, 2 ) + sig_par[16] * sig_the + sig_par[17];
326
327 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin;
328 }
329 else if ( sigmaflag == 4 )
330 {
331 double x = beta_G;
332 double nhit = (double)Nohit;
333 double sigma_bg = 1.0;
334 if ( x > 0.3 )
335 sigma_bg = sig_par[0] * exp( sig_par[1] * x ) +
336 sig_par[2] * exp( sig_par[3] * pow( x, 0.25 ) ) + sig_par[4];
337 else sigma_bg = sig_par[5] * exp( sig_par[6] * x ) + sig_par[7];
338
339 double cor_nhit = 1.0;
340 if ( nhit < 5 ) nhit = 5;
341 if ( nhit <= 35 )
342 cor_nhit = sig_par[8] * pow( nhit, 4 ) + sig_par[9] * pow( nhit, 3 ) +
343 sig_par[10] * pow( nhit, 2 ) + sig_par[11] * nhit + sig_par[12];
344
345 double cor_sin = 1.0;
346 if ( sig_the > 0.99 ) sig_the = 0.99;
347 cor_sin = sig_par[13] * pow( sig_the, 4 ) + sig_par[14] * pow( sig_the, 3 ) +
348 sig_par[15] * pow( sig_the, 2 ) + sig_par[16] * sig_the + sig_par[17];
349
350 if ( trkalg == 1 ) ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin;
351 else ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * sig_par[18];
352 }
353 else if ( sigmaflag == 5 )
354 {
355 double x = beta_G;
356 double nhit = (double)Nohit;
357 double sigma_bg = 1.0;
358 if ( x > 0.3 )
359 sigma_bg = sig_par[0] * exp( sig_par[1] * x ) +
360 sig_par[2] * exp( sig_par[3] * pow( x, 0.25 ) ) + sig_par[4];
361 else sigma_bg = sig_par[5] * exp( sig_par[6] * x ) + sig_par[7];
362
363 double cor_nhit = 1.0;
364 if ( nhit < 5 ) nhit = 5;
365 if ( nhit <= 35 )
366 cor_nhit = sig_par[8] * pow( nhit, 4 ) + sig_par[9] * pow( nhit, 3 ) +
367 sig_par[10] * pow( nhit, 2 ) + sig_par[11] * nhit + sig_par[12];
368 double cor_sin = 1.0;
369 if ( sig_the > 0.99 ) sig_the = 0.99;
370 cor_sin = sig_par[13] * pow( sig_the, 4 ) + sig_par[14] * pow( sig_the, 3 ) +
371 sig_par[15] * pow( sig_the, 2 ) + sig_par[16] * sig_the + sig_par[17];
372
373 double cor_t0 = 1;
374 if ( t0 > 1200 ) t0 = 1200;
375 if ( t0 > 800 ) cor_t0 = sig_par[18] * pow( t0, 2 ) + sig_par[19] * t0 + sig_par[20];
376
377 if ( trkalg == 1 ) ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
378 else ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[21];
379 }
380 }
381
382 double dedx_correc = dedx;
383 chi_dedx[it] = ( dedx_correc - dedx_exp[it] ) / ex_sigma[it];
384 chi2 = chi_dedx[it] * chi_dedx[it];
385 ndf = 1;
386 pid_prob[it] = prob_( &chi2, &ndf );
387 // if(it ==0 ) cout<<"runflag: "<<runflag<<" dedx : "<<dedx<<" chi_dedx: "
388 //<<chi_dedx[it] <<" ptrk: "<<mom<<endl;
389 if ( it == -999 )
390 { // here a debug flag
391 std::cout << " mom = " << mom << "exp" << dedx_exp[it] << " sigma " << ex_sigma[it]
392 << " prob " << pid_prob[it] << std::endl;
393 }
394 if ( pid_prob[it] > max_prob )
395 {
396 max_prob = pid_prob[it];
397 Nmax_prob = it;
398 }
399 }
400 else
401 {
402 dedx_exp[it] = 0.0;
403 ex_sigma[it] = 1000.0;
404 pid_prob[it] = 0.0;
405 chi_dedx[it] = 999.0;
406 } // if Nohit > 0
407 }
408}
prob_
float prob_(float *, int *)
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:209
C
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations C
Definition RRes.h:29

Referenced by DedxCalib::set_dEdx().

◆ dedx_pid_exp_old()

void dedx_pid_exp_old ( int landau,
int runflag,
float dedx,
int Nohit,
float mom,
float theta,
float t0,
float lsamp,
double dedx_exp[5],
double ex_sigma[5],
double pid_prob[5],
double chi_dedx[5] )

Definition at line 55 of file DedxCalibParameters.cxx.

57 {
58 double par[5], sigma_par[4], sigma_index_nhit, sigma_index_sin;
59
60 if ( runflag == 1 )
61 {
62 par[0] = HV1_curvep0;
63 par[1] = HV1_curvep1;
64 par[2] = HV1_curvep2;
65 par[3] = HV1_curvep3;
66 par[4] = HV1_curvep4;
67 sigma_par[0] = HV1_sigmap0;
68 sigma_par[1] = HV1_sigmap1;
69 sigma_par[2] = HV1_sigmap2;
70 sigma_par[3] = HV1_sigmap3;
71 sigma_index_nhit = HV1_index_nhit;
72 sigma_index_sin = HV1_index_sin;
73 }
74 else if ( runflag == 2 )
75 {
76 par[0] = HV2_curvep0;
77 par[1] = HV2_curvep1;
78 par[2] = HV2_curvep2;
79 par[3] = HV2_curvep3;
80 par[4] = HV2_curvep4;
81 sigma_par[0] = HV2_sigmap0;
82 sigma_par[1] = HV2_sigmap1;
83 sigma_par[2] = HV2_sigmap2;
84 sigma_par[3] = HV2_sigmap3;
85 sigma_index_nhit = HV2_index_nhit;
86 sigma_index_sin = HV2_index_sin;
87 }
88
89 const int par_cand( 5 );
90 const float Charge_Mass[par_cand] = { 0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
91 double beta_G, beta, betterm, bethe_B, sig_param;
92
93 int Nmax_prob( 0 );
94 float max_prob( -0.01 );
95 int ndf;
96 float chi2;
97
98 for ( int it = 0; it < par_cand; it++ )
99 {
100 beta_G = mom / Charge_Mass[it];
101 beta = beta_G / sqrt( 1 + ( beta_G ) * ( beta_G ) );
102 betterm = par[1] - log( par[2] + pow( 1 / beta_G, par[4] ) );
103 bethe_B = par[0] / pow( beta, par[3] ) * betterm - par[0];
104
105 if ( Nohit > 0 )
106 {
107 dedx_exp[it] = bethe_B;
108 double sig_the = std::sin( (double)theta );
109
110 if ( runflag < 3 && runflag > 0 )
111 {
112 if ( landau == 0 )
113 {
114 sig_param = 1.6 * std::sin( (double)theta ) / ( lsamp * double( Nohit ) );
115 ex_sigma[it] = 0.05 * bethe_B * sqrt( 50.0 * sig_param );
116 }
117 else
118 {
119 // currently use one sigmap0
120 if ( beta_G < 4 )
121 { sig_param = sigma_par[1] + sigma_par[2] * std::pow( beta_G, sigma_par[3] ); }
122 else { sig_param = sigma_par[0]; }
123 // double sig_the=std::sin( (double)theta );
124 sig_the = std::pow( sig_the, sigma_index_sin );
125 double sig_n;
126 sig_n = 35.0 / double( Nohit );
127 sig_n = std::pow( sig_n, sigma_index_nhit );
128 ex_sigma[it] = sig_param * sig_the * sig_n;
129 }
130 }
131
132 double dedx_correc;
133 if ( runflag == 2 ) dedx_correc = SpaceChargeCorrec( theta, mom, it, dedx );
134 else dedx_correc = dedx;
135 chi_dedx[it] = ( dedx_correc - dedx_exp[it] ) / ex_sigma[it];
136 chi2 = chi_dedx[it] * chi_dedx[it];
137 ndf = 1;
138 pid_prob[it] = prob_( &chi2, &ndf );
139 // if(it ==0 ) cout<<"runflag: "<<runflag<<" dedx : "<<dedx<<" chi_dedx: "
140 //<<chi_dedx[it] <<" ptrk: "<<mom<<endl;
141 if ( it == -999 )
142 { // here a debug flag
143 std::cout << " mom = " << mom << "exp" << dedx_exp[it] << " sigma " << ex_sigma[it]
144 << " prob " << pid_prob[it] << std::endl;
145 }
146 if ( pid_prob[it] > max_prob )
147 {
148 max_prob = pid_prob[it];
149 Nmax_prob = it;
150 }
151 }
152 else
153 {
154 dedx_exp[it] = 0.0;
155 ex_sigma[it] = 1000.0;
156 pid_prob[it] = 0.0;
157 chi_dedx[it] = 999.0;
158 }
159 }
160}
SpaceChargeCorrec
double SpaceChargeCorrec(double m_theta, double mom, int Particle, double dEdx)
Definition DedxCalibParameters.cxx:410
HV2_curvep2
const double HV2_curvep2
Definition DedxCalibParameters.h:28
HV2_index_nhit
const double HV2_index_nhit
Definition DedxCalibParameters.h:32
HV2_curvep0
const double HV2_curvep0
Definition DedxCalibParameters.h:26
HV2_sigmap3
const double HV2_sigmap3
Definition DedxCalibParameters.h:36
HV2_sigmap0
const double HV2_sigmap0
Definition DedxCalibParameters.h:33
HV1_curvep1
const double HV1_curvep1
Definition DedxCalibParameters.h:14
HV1_sigmap0
const double HV1_sigmap0
Definition DedxCalibParameters.h:20
HV2_curvep3
const double HV2_curvep3
Definition DedxCalibParameters.h:29
HV2_sigmap1
const double HV2_sigmap1
Definition DedxCalibParameters.h:34
HV1_sigmap3
const double HV1_sigmap3
Definition DedxCalibParameters.h:23
HV1_index_nhit
const double HV1_index_nhit
Definition DedxCalibParameters.h:19
HV1_sigmap1
const double HV1_sigmap1
Definition DedxCalibParameters.h:21
HV2_curvep1
const double HV2_curvep1
Definition DedxCalibParameters.h:27
HV2_curvep4
const double HV2_curvep4
Definition DedxCalibParameters.h:30
HV1_sigmap2
const double HV1_sigmap2
Definition DedxCalibParameters.h:22
HV1_curvep3
const double HV1_curvep3
Definition DedxCalibParameters.h:16
HV2_sigmap2
const double HV2_sigmap2
Definition DedxCalibParameters.h:35
HV1_curvep0
const double HV1_curvep0
Definition DedxCalibParameters.h:13
HV1_index_sin
const double HV1_index_sin
Definition DedxCalibParameters.h:18
HV1_curvep4
const double HV1_curvep4
Definition DedxCalibParameters.h:17
HV2_index_sin
const double HV2_index_sin
Definition DedxCalibParameters.h:31
HV1_curvep2
const double HV1_curvep2
Definition DedxCalibParameters.h:15

Referenced by DedxCalib::set_dEdx().

◆ Landau()

double Landau ( double * x,
double * par )

Definition at line 30 of file DedxCalibParameters.cxx.

30 {
31 double fitval = TMath::Landau( x[0], par[0], par[1], kFALSE );
32 return fitval;
33}

Referenced by DedxCalibDocaEAng::AnalyseHists(), DedxCalibEAng::AnalyseHists(), DedxCalibLayerGain::AnalyseHists(), and DedxCalibWireGain::AnalyseHists().

◆ landaun()

double landaun ( double * x,
double * par )

Definition at line 22 of file DedxCalibParameters.cxx.

22 {
23
24 double kk = ( x[0] - par[1] ) / par[2];
25 double fterm = kk + exp( -1 * kk );
26 double fitval = par[0] * exp( -0.5 * fterm );
27 return fitval;
28}

Referenced by DedxCalibDocaEAng::AnalyseHists(), DedxCalibEAng::AnalyseHists(), DedxCalibLayerGain::AnalyseHists(), and DedxCalibWireGain::AnalyseHists().

◆ mylan()

double mylan ( double * x,
double * par )

Definition at line 15 of file DedxCalibParameters.cxx.

15 {
16 double kk = ( x[0] - par[1] ) / par[2];
17 double fterm = kk + exp( -1 * kk );
18 double fitval = par[0] * exp( par[3] * fterm );
19 return fitval;
20}

Referenced by DedxCalibDocaEAng::AnalyseHists(), DedxCalibEAng::AnalyseHists(), DedxCalibLayerGain::AnalyseHists(), and DedxCalibWireGain::AnalyseHists().

◆ prob_()

float prob_ ( float * ,
int *  )

Referenced by dedx_pid_exp(), MdcDedxCorrection::dedx_pid_exp(), dedx_pid_exp_old(), MdcDedxCorrection::dedx_pid_exp_old(), and Lpav::prob().

◆ SpaceChargeCorrec()

double SpaceChargeCorrec ( double m_theta,
double mom,
int Particle,
double dEdx )

Definition at line 410 of file DedxCalibParameters.cxx.

410 {
411 const int par_cand( 5 );
412 const float Charge_Mass[par_cand] = { 0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
413 double beta_G;
414 double e_Par[5] = { 143.349, 1.7315, 0.192616, 2.90437, 1.08248 };
415 double Beta_Gamma[22] = { 0.373026, 0.479605, 0.586184, 0.692763, 0.799342, 782.779,
416 1565.56, 2348.34, 17.2727, 18.1245, 1.43297, 2.14946,
417 12.1803, 13.6132, 6.62515, 10.4109, 14.1967, 17.9825,
418 21.7683, 26.0274, 30.7596, 35.4919 };
419 double K_par[22] = { 4.64411e-05, 5.86544e-05, 8.05289e-05, 8.46981e-05, 8.92014e-05,
420 4.74517e-05, 4.51684e-05, 5.32732e-05, 6.12803e-05, 6.14592e-05,
421 8.08608e-05, 6.73184e-05, 5.46448e-05, 6.1377e-05, 6.57385e-05,
422 7.03053e-05, 6.61171e-05, 6.86824e-05, 6.246e-05, 7.25988e-05,
423 7.11034e-05, 6.24924e-05 };
424 double D_par[22] = { 0.0871504, 0.0956379, 0.117193, 0.118647, 0.127203, 0.0566449,
425 0.0529198, 0.0642525, 0.0764562, 0.081341, 0.0952263, 0.0987536,
426 0.0639901, 0.0845994, 0.0777062, 0.0823206, 0.0783874, 0.079537,
427 0.0815792, 0.0849875, 0.0824751, 0.0776296 };
428 double DSqr_par[22] = { 0.00759519, 0.0091466, 0.0137341, 0.0140772, 0.0161807,
429 0.00320864, 0.00280051, 0.00412839, 0.00584555, 0.00661636,
430 0.00906805, 0.00975227, 0.00409473, 0.00715706, 0.00603826,
431 0.00677668, 0.00614458, 0.00632613, 0.00665516, 0.00722288,
432 0.00680214, 0.00602635 };
433
434 beta_G = mom / Charge_Mass[Particle];
435 if ( beta_G < 0.3 ) beta_G = 0.3;
436 double bet = beta_G / TMath::Sqrt( beta_G * beta_G + 1 );
437 double fterm = TMath::Log( e_Par[2] + 1 / pow( beta_G, e_Par[4] ) );
438 double fitval =
439 e_Par[0] / pow( bet, e_Par[3] ) * ( e_Par[1] - pow( bet, e_Par[3] ) - fterm );
440 TGraphErrors* gr1 = new TGraphErrors( 22, Beta_Gamma, K_par, 0, 0 );
441 TGraphErrors* gr2 = new TGraphErrors( 22, Beta_Gamma, DSqr_par, 0, 0 );
442
443 double par[3];
444 par[0] = fitval;
445 par[1] = gr1->Eval( m_theta );
446 par[2] = gr2->Eval( m_theta );
447 Double_t y = fabs( cos( m_theta ) );
448 double electron_par[3] = { 334.032, 6.20658e-05, 0.00525673 };
449 double arg = TMath::Sqrt( y * y + par[2] );
450 // double cal_factor =par[0]*TMath::Exp(-(par[1]* par[0])/arg);
451 double cal_factor = TMath::Exp( -( par[1] * par[0] ) / arg );
452 double arg_electron = TMath::Sqrt( y * y + electron_par[2] );
453 // double electron_factor = electron_par[0]*TMath::Exp(-(electron_par[1]*
454 // electron_par[0])/arg);
455 double electron_factor = TMath::Exp( -( electron_par[1] * electron_par[0] ) / arg_electron );
456 // cout<<"cal_factor = "<<cal_factor<<" electron_factor = "<<electron_factor<<endl;
457 double dedx_cal = dEdx / ( cal_factor / electron_factor );
458 // double dedx_cal = dEdx/cal_factor;
459 // cout<<"m_theta= "<<m_theta<<" y ="<<y<<" beta_G = "<<beta_G <<" dEdx =
460 // "<<dEdx<<" cal dedx = "<<dedx_cal<<endl;
461 delete gr1;
462 delete gr2;
463 return dedx_cal;
464}
arg
double arg(const EvtComplex &c)
Definition EvtComplex.hh:195
Particle
HepMC::GenParticle Particle
Definition GenIMCselector.h:16
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187

Referenced by dedx_pid_exp_old().

◆ vavden_()

float vavden_ ( double * )

Referenced by Vavilov().

◆ Vavilov()

double Vavilov ( double * x,
double * par )

Definition at line 35 of file DedxCalibParameters.cxx.

35 {
36 double kappa, beta2;
37 int mode = 0;
38 vavset_( &par[0], &par[1], &mode );
39 double fitval = vavden_( &x[0] );
40 return fitval;
41}
vavden_
float vavden_(double *)
vavset_
float vavset_(double *, double *, int *)

Referenced by DedxCalibDocaEAng::AnalyseHists(), DedxCalibEAng::AnalyseHists(), DedxCalibLayerGain::AnalyseHists(), and DedxCalibWireGain::AnalyseHists().

◆ vavset_()

float vavset_ ( double * ,
double * ,
int *  )

Referenced by Vavilov().