BOSS 7.0.9
BESIII Offline Software System
Loading...
Searching...
No Matches
Functions
DedxCalibParameters.cxx File Reference
#include <TMath.h>
#include <cmath>
#include <vector>
#include <iostream>
#include "TGraphErrors.h"
#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 50 of file DedxCalibParameters.cxx.

51{
52 double delta, a, b, c, d, fitval;
53 a = TMath::Sqrt(log(4.0));
54 b = (x[0]-par[1])/par[2];
55 c = TMath::Power(par[3],2.0);
56 delta = (1+ TMath::SinH(par[3]*a)/a)*b;
57 d = TMath::Power(log(delta),2.0)/c+c;
58 fitval = par[0]*TMath::Exp(-0.5*d);
59 return fitval;
60}
delta
const double delta
Definition: CheckRung/check/gainVsDoca.cxx:32
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
b
const double b
Definition: slope.cxx:9

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 165 of file DedxCalibParameters.cxx.

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

Referenced by draw_check_curve(), get_hist_range(), and 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 63 of file DedxCalibParameters.cxx.

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

35{
36 double fitval = TMath::Landau(x[0], par[0], par[1], kFALSE);
37 return fitval;
38}

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

◆ landaun()

double landaun ( double *  x,
double *  par 
)

Definition at line 25 of file DedxCalibParameters.cxx.

26{
27
28 double kk=(x[0]-par[1])/par[2];
29 double fterm=kk+exp(-1*kk);
30 double fitval=par[0]*exp(-0.5*fterm);
31 return fitval;
32}

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

◆ mylan()

double mylan ( double *  x,
double *  par 
)

Definition at line 17 of file DedxCalibParameters.cxx.

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

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

◆ prob_()

float prob_ ( float *  ,
int *   
)

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

◆ SpaceChargeCorrec()

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

Definition at line 406 of file DedxCalibParameters.cxx.

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

Referenced by dedx_pid_exp_old().

◆ vavden_()

float vavden_ ( double *  )

Referenced by Vavilov().

◆ Vavilov()

double Vavilov ( double *  x,
double *  par 
)

Definition at line 41 of file DedxCalibParameters.cxx.

42{
43 double kappa, beta2;
44 int mode = 0;
45 vavset_(&par[0], &par[1], &mode);
46 double fitval = vavden_(&x[0]);
47 return fitval;
48}
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().