BOSS 7.0.9
BESIII Offline Software System
Loading...
Searching...
No Matches
histgen/curve.h
Go to the documentation of this file.
1#ifndef CURVE_H
2#define CURVE_H
3
4/// Curve parameter vars
5int vflag[3] = { 0, 0, 0 };
6std::vector<double> cpar;
7std::vector<double> spar;
8#include "global.h"
9
10void ReadCurPara(int index) {
11
12 TFile* fcurve = new TFile(curve_files[index].c_str());
13
14 TTree* curve = (TTree*)fcurve->Get("curve");
15 TTree* sigma = (TTree*)fcurve->Get("sigma");
16 double c[100], s[100];
17 int csize, ssize;
18 curve->SetBranchAddress("curve", c);
19 curve->SetBranchAddress("CurveSize", &csize);
20 sigma->SetBranchAddress("sigma", s);
21 sigma->SetBranchAddress("SigmaSize", &ssize);
22
23 curve->GetEntry(0);
24 sigma->GetEntry(0);
25
26 vflag[0] = (int)c[0];
27 vflag[1] = (int)s[0];
28 if (c[0] < 0 && s[0] < 0)
29 vflag[2] = 1;
30 else
31 vflag[2] = 0;
32
33 cpar.clear();
34 spar.clear();
35
36 for (int i = 1; i < csize; i++) {
37 cpar.push_back(c[i]);
38 }
39 for (int i = 1; i < ssize; i++) {
40 spar.push_back(s[i]);
41 }
42 delete fcurve;
43}
44
45void dedx_pid_exp(
46 double mom,
47 double theta,
48 double Nohit,
49 double t0,
50 double dedx_exp[5],
51 double ex_sigma[5]) {
52
53 const int par_cand(5);
54 const float Charge_Mass[par_cand] = { 0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
55 double beta_G, beta, betterm, bethe_B;
56 const int trkalg(1);
57
58 int dedxflag = vflag[0];
59 int sigmaflag = vflag[1];
60 bool ifMC = false;
61 if(vflag[2] == 1) ifMC = true;
62
63 for( int it = 0; it < par_cand; it++ ) {
64 beta_G = mom/Charge_Mass[it];
65
66 if(dedxflag == 1){
67 beta = beta_G/sqrt(1+(beta_G)*(beta_G));
68 betterm = cpar[1]-log(cpar[2]+pow(1/beta_G,cpar[4]));
69 bethe_B = cpar[0]/pow(beta,cpar[3])*betterm-cpar[0];
70 }
71 else if(dedxflag == 2) {
72 double A=0, B=0,C=0;
73 double x = beta_G;
74 if(x<4.5)
75 A=1;
76 else if(x<10)
77 B=1;
78 else
79 C=1;
80 double partA = cpar[0]*pow(sqrt(x*x+1),cpar[2])/pow(x,cpar[2])*(cpar[1]-cpar[5]*log(pow(1/x,cpar[3])) ) -
81 cpar[4]+exp(cpar[6]+cpar[7]*x);
82 double partB = cpar[8]*pow(x,3)+cpar[9]*pow(x,2)+cpar[10]*x+cpar[11];
83 double partC = -cpar[12]*log(cpar[15]+pow(1/x,cpar[13]))+cpar[14];
84 bethe_B = 550*(A*partA+B*partB+C*partC);
85 }
86
87 if(ifMC){
88 double A=0, B=0,C=0;
89 double x = beta_G;
90 if(x<4.5)
91 A=1;
92 else if(x<10)
93 B=1;
94 else
95 C=1;
96 double partA = cpar[0]*pow(sqrt(x*x+1),cpar[2])/pow(x,cpar[2])*(cpar[1]-cpar[5]*log(pow(1/x,cpar[3])) ) -
97 cpar[4]+exp(cpar[6]+cpar[7]*x);
98 double partB = cpar[8]*pow(x,3)+cpar[9]*pow(x,2)+cpar[10]*x+cpar[11];
99 double partC = -cpar[12]*log(cpar[15]+pow(1/x,cpar[13]))+cpar[14];
100 bethe_B = 550*(A*partA+B*partB+C*partC);
101 //for fermi plateau ( the electron region) we just use 1.0
102 if(beta_G> 100) bethe_B=550*1.0;
103 }
104
105 if (Nohit > 0) {
106 dedx_exp[it] = bethe_B;
107 double sig_the = std::sin((double)theta);
108 double f_betagamma, g_sinth, h_nhit, i_t0;
109 if (ifMC){
110 if(sigmaflag == 6){
111 double nhit = (double)Nohit;
112 double sigma_bg=1.0;
113 double ndedx=bethe_B/550;
114 sigma_bg= spar[0]+spar[1]*ndedx+spar[2]*ndedx*ndedx;
115 double cor_nhit = 1.0;
116 if (nhit < 5) nhit = 5;
117 if (nhit <= 35)
118 cor_nhit = spar[3]*pow(nhit,4)+spar[4]*pow(nhit,3)+spar[5]*pow(nhit,2) +
119 spar[6]*nhit+spar[7];
120
121 double cor_sin= 1.0;
122 if(sig_the>0.99) sig_the=0.99;
123 cor_sin = spar[8]*pow(sig_the,4)+spar[9]*pow(sig_the,3) +
124 spar[10]*pow(sig_the,2)+spar[11]*sig_the+spar[12];
125
126 //sigma vs t0
127 double cor_t0 = 1.0;
128
129 //calculate sigma
130 if (trkalg == 1)
131 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
132 else
133 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * spar[13];
134 }// end of if sigmaversion == 6
135
136 else{
137 double x = beta_G;
138 double nhit = (double)Nohit;
139 double sigma_bg = 1.0;
140
141 if (x > 0.3)
142 sigma_bg = spar[0]*exp(spar[1]*x)+spar[2]*exp(spar[3]*pow(x,0.25))+spar[4];
143 else
144 sigma_bg = spar[5]*exp(spar[6]*x)+ spar[7];
145
146 double cor_nhit = 1.0;
147 if (nhit < 5) nhit = 5;
148 if (nhit <= 35)
149 cor_nhit = spar[8]*pow(nhit,4)+spar[9]*pow(nhit,3)+spar[10]*pow(nhit,2) +
150 spar[11]*nhit+spar[12];
151
152 //sigma vs sintheta
153 double cor_sin = 1.0;
154 if(sig_the>0.99) sig_the=0.99;
155 cor_sin = spar[13]*pow(sig_the,4)+spar[14]*pow(sig_the,3) +
156 spar[15]*pow(sig_the,2)+spar[16]*sig_the+spar[17];
157
158 //sigma vs t0
159 double cor_t0 = 1.0;
160
161 //calculate sigma
162 if (trkalg == 1)
163 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
164 else
165 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * spar[18];
166 } // end of else
167 }//end of MC
168 else {
169 if(sigmaflag == 1) {
170 f_betagamma=spar[0]*pow(beta_G,spar[1])+spar[2];
171 g_sinth=(spar[3]*sig_the*sig_the+spar[4])/(spar[3]*spar[5]*spar[5]+spar[4]);
172 h_nhit=(spar[6]*Nohit*Nohit+spar[7]*Nohit+spar[8]) /
173 (spar[6]*spar[9]*spar[9]+spar[7]*spar[9]+spar[8]);
174 if(spar[13] != 0)
175 i_t0 = (spar[10]*t0*t0+spar[11]*t0+spar[12]) /
176 (spar[10]*spar[13]*spar[13]+spar[11]*spar[13]+spar[12]);
177 else if(spar[13] == 0)
178 i_t0 =1;
179 //cout<<"f_betagamma : "<<f_betagamma<<" g_sinth : "<<g_sinth<<" h_nhit : "
180 //<<h_nhit<<" i_t0 : "<<i_t0<<endl;
181 ex_sigma[it]= f_betagamma* g_sinth * h_nhit * i_t0;
182 }
183 else if(sigmaflag == 2) {
184 double x = beta_G;
185 double nhit = (double)Nohit;
186 double sigma_bg=1.0;
187 if (x > 0.3)
188 sigma_bg = spar[0]*exp(spar[1]*x)+spar[2]*exp(spar[3]*pow(x,0.25))+spar[4];
189 else
190 sigma_bg = spar[5]*exp(spar[6]*x)+ spar[7];
191
192 double cor_nhit=1.0;
193 if(nhit<5) nhit=5;
194 if (nhit <= 35)
195 cor_nhit = spar[8]*pow(nhit,4)+spar[9]*pow(nhit,3)+spar[10]*pow(nhit,2) + spar[11]*nhit+spar[12];
196
197 double cor_sin= 1.0;
198 if(sig_the>0.99) sig_the = 0.99;
199 cor_sin = spar[13]*pow(sig_the,4)+spar[14]*pow(sig_the,3) +
200 spar[15]*pow(sig_the,2)+spar[16]*sig_the+spar[17];
201
202 double cor_t0 = 1;
203 if (t0 > 1200) t0 = 1200;
204 if (t0 > 800)
205 cor_t0 = spar[18]*pow(t0,2)+spar[19]*t0+spar[20];
206
207 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
208 }
209 else if(sigmaflag == 3) {
210 double x= beta_G;
211 double nhit = (double)Nohit;
212 double sigma_bg=1.0;
213 if (x > 0.3)
214 sigma_bg = spar[0]*exp(spar[1]*x)+spar[2]*exp(spar[3]*pow(x,0.25))+spar[4];
215 else
216 sigma_bg = spar[5]*exp(spar[6]*x)+ spar[7];
217
218 double cor_nhit = 1.0;
219 if (nhit < 5) nhit = 5;
220 if (nhit <= 35)
221 cor_nhit = spar[8]*pow(nhit,4)+spar[9]*pow(nhit,3)+spar[10]*pow(nhit,2) + spar[11]*nhit+spar[12];
222
223 double cor_sin= 1.0;
224 if(sig_the>0.99) sig_the=0.99;
225 cor_sin = spar[13]*pow(sig_the,4)+spar[14]*pow(sig_the,3) +
226 spar[15]*pow(sig_the,2)+spar[16]*sig_the+spar[17];
227
228 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
229 }
230 else if(sigmaflag == 4) {
231 double x= beta_G;
232 double nhit = (double)Nohit;
233 double sigma_bg=1.0;
234 if (x > 0.3)
235 sigma_bg = spar[0]*exp(spar[1]*x)+spar[2]*exp(spar[3]*pow(x,0.25))+spar[4];
236 else
237 sigma_bg = spar[5]*exp(spar[6]*x)+ spar[7];
238
239 double cor_nhit = 1.0;
240 if (nhit < 5) nhit = 5;
241 if (nhit <= 35)
242 cor_nhit = spar[8]*pow(nhit,4)+spar[9]*pow(nhit,3)+spar[10]*pow(nhit,2) + spar[11]*nhit+spar[12];
243
244 double cor_sin= 1.0;
245 if(sig_the>0.99) sig_the=0.99;
246 cor_sin = spar[13]*pow(sig_the,4)+spar[14]*pow(sig_the,3) +
247 spar[15]*pow(sig_the,2)+spar[16]*sig_the+spar[17];
248
249 if(trkalg==1)
250 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
251 else
252 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*spar[18];
253 }
254 else if(sigmaflag == 5) {
255 double x = beta_G;
256 double nhit = (double)Nohit;
257 double sigma_bg=1.0;
258 if (x > 0.3)
259 sigma_bg = spar[0]*exp(spar[1]*x)+spar[2]*exp(spar[3]*pow(x,0.25))+spar[4];
260 else
261 sigma_bg = spar[5]*exp(spar[6]*x)+ spar[7];
262
263 double cor_nhit=1.0;
264 if(nhit<5) nhit=5;
265 if (nhit <= 35)
266 cor_nhit = spar[8]*pow(nhit,4)+spar[9]*pow(nhit,3)+spar[10]*pow(nhit,2) +
267 spar[11]*nhit+spar[12];
268 double cor_sin= 1.0;
269 if(sig_the>0.99) sig_the = 0.99;
270 cor_sin = spar[13]*pow(sig_the,4)+spar[14]*pow(sig_the,3) +
271 spar[15]*pow(sig_the,2)+spar[16]*sig_the+spar[17];
272
273 double cor_t0 = 1;
274 if (t0 > 1200) t0 = 1200;
275 if (t0 > 800)
276 cor_t0 = spar[18]*pow(t0,2)+spar[19]*t0+spar[20];
277
278 if(trkalg==1)
279 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
280 else
281 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0*spar[21];
282 }
283 else if(sigmaflag == 6) {
284 double x= bethe_B/550;
285 double nhit = (double)Nohit;
286 double sigma_bg=1.0;
287
288 sigma_bg= spar[0]+spar[1]*x+spar[2]*x*x;
289
290 double cor_nhit = 1.0;
291 if (nhit < 5) nhit = 5;
292 if (nhit <= 35)
293 cor_nhit = spar[3]*pow(nhit,4)+spar[4]*pow(nhit,3)+spar[5]*pow(nhit,2) + spar[6]*nhit+spar[7];
294
295 double cor_sin= 1.0;
296 if(sig_the>0.99) sig_the=0.99;
297 cor_sin = spar[8]*pow(sig_the,4)+spar[9]*pow(sig_the,3) +
298 spar[10]*pow(sig_the,2)+spar[11]*sig_the+spar[12];
299
300 if(trkalg==1)
301 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
302 else
303 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*spar[13];
304 }
305 else if(sigmaflag == 7) {
306 double x= bethe_B/550;
307 double nhit = (double)Nohit;
308 double sigma_bg=1.0;
309
310 sigma_bg= spar[0]+spar[1]*x+spar[2]*x*x;
311
312
313 double cor_nhit=1.0;
314 if(nhit<5) nhit=5;
315 if (nhit <= 35)
316 cor_nhit = spar[3]*pow(nhit,4)+spar[4]*pow(nhit,3)+spar[5]*pow(nhit,2) +
317 spar[6]*nhit+spar[7];
318 double cor_sin= 1.0;
319 if(sig_the>0.99) sig_the = 0.99;
320 cor_sin = spar[8]*pow(sig_the,4)+spar[9]*pow(sig_the,3) +
321 spar[10]*pow(sig_the,2)+spar[11]*sig_the+spar[12];
322
323 double cor_t0 = 1;
324 if (t0 > 1200) t0 = 1200;
325 if (t0 > 800)
326 cor_t0 = spar[13]*pow(t0,2)+spar[14]*t0+spar[15];
327
328 if(trkalg==1)
329 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
330 else
331 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0*spar[16];
332 }
333 }
334 } // end of Nohit > 0
335 else{
336 dedx_exp[it] = 0.0;
337 ex_sigma[it] = 1000.0;
338 } //if Nohit > 0
339 } // end of loop in particle type
340}
341
342void dedx_pid_exp(int vflag[3], double dedx, int trkalg,
343 int Nohit, double mom, double theta, double t0,
344 double dedx_exp[5], double ex_sigma[5],
345 double pid_prob[5], double chi_dedx[5],
346 std::vector<double> & par, std::vector<double> & sig_par)
347{
348 const int par_cand( 5 );
349 const float Charge_Mass[par_cand] = {0.00051100, 0.10566, 0.13957, 0.4937, 0.93827 };
350 double beta_G, beta, betterm, bethe_B;
351
352 int dedxflag = vflag[0];
353 int sigmaflag = vflag[1];
354 bool ifMC = false;
355 if(vflag[2] == 1) ifMC = true;
356
357 int Nmax_prob(0);
358 float max_prob(-0.01);
359 int ndf;
360 float chi2;
361
362 for( int it = 0; it < par_cand; it++ ) {
363 beta_G = mom/Charge_Mass[it];
364
365 if(dedxflag == 1){
366 beta = beta_G/sqrt(1+(beta_G)*(beta_G));
367 betterm = par[1]-log(par[2]+pow(1/beta_G,par[4]));
368 bethe_B = par[0]/pow(beta,par[3])*betterm-par[0];
369 }
370 else if(dedxflag == 2) {
371 double A=0, B=0,C=0;
372 double x = beta_G;
373 if(x<4.5)
374 A=1;
375 else if(x<10)
376 B=1;
377 else
378 C=1;
379 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);
380 double partB = par[8]*pow(x,3)+par[9]*pow(x,2)+par[10]*x+par[11];
381 double partC = -par[12]*log(par[15]+pow(1/x,par[13]))+par[14];
382 bethe_B = 550*(A*partA+B*partB+C*partC);
383 //for fermi plateau ( the electron region) we just use 1.0
384 if(beta_G> 100) bethe_B=550*1.0;
385 }
386
387 if (ifMC) {
388 double A=0, B=0,C=0;
389 double x = beta_G;
390 if(x<4.5)
391 A=1;
392 else if(x<10)
393 B=1;
394 else
395 C=1;
396 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);
397 double partB = par[8]*pow(x,3)+par[9]*pow(x,2)+par[10]*x+par[11];
398 double partC = -par[12]*log(par[15]+pow(1/x,par[13]))+par[14];
399 bethe_B = 550*(A*partA+B*partB+C*partC);
400 //for fermi plateau ( the electron region) we just use 1.0
401 if(beta_G> 100) bethe_B=550*1.0;
402 }
403
404
405 if (Nohit > 0) {
406 dedx_exp[it] = bethe_B;
407 double sig_the = std::sin((double)theta);
408 double f_betagamma, g_sinth, h_nhit, i_t0;
409
410 if (ifMC) {
411
412 if(sigmaflag == 6) {
413
414 double nhit = (double)Nohit;
415 double sigma_bg=1.0;
416 double ndedx=bethe_B/550;
417 sigma_bg= sig_par[0]+sig_par[1]*ndedx+sig_par[2]*ndedx*ndedx;
418
419 double cor_nhit = 1.0;
420 if (nhit < 5) nhit = 5;
421 if (nhit <= 35)
422 cor_nhit = sig_par[3]*pow(nhit,4)+sig_par[4]*pow(nhit,3)+sig_par[5]*pow(nhit,2) +
423 sig_par[6]*nhit+sig_par[7];
424
425 double cor_sin= 1.0;
426 if(sig_the>0.99) sig_the=0.99;
427 cor_sin = sig_par[8]*pow(sig_the,4)+sig_par[9]*pow(sig_the,3) +
428 sig_par[10]*pow(sig_the,2)+sig_par[11]*sig_the+sig_par[12];
429
430 //sigma vs t0
431 double cor_t0 = 1.0;
432
433 //calculate sigma
434 if (trkalg == 1)
435 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
436 else
437 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[13];
438
439 }
440
441 else{
442 double x= beta_G;
443 double nhit = (double)Nohit;
444 double sigma_bg=1.0;
445
446 if (x > 0.3)
447 sigma_bg = sig_par[0]*exp(sig_par[1]*x)+sig_par[2]*exp(sig_par[3]*pow(x,0.25))+sig_par[4];
448 else
449 sigma_bg = sig_par[5]*exp(sig_par[6]*x)+ sig_par[7];
450
451
452 double cor_nhit = 1.0;
453 if (nhit < 5) nhit = 5;
454 if (nhit <= 35)
455 cor_nhit = sig_par[8]*pow(nhit,4)+sig_par[9]*pow(nhit,3)+sig_par[10]*pow(nhit,2) +
456 sig_par[11]*nhit+sig_par[12];
457
458 double cor_sin= 1.0;
459 if(sig_the>0.99) sig_the=0.99;
460 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
461 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
462
463 //sigma vs t0
464 double cor_t0 = 1.0;
465
466 //calculate sigma
467 if (trkalg == 1)
468 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0;
469 else
470 ex_sigma[it] = 550 * sigma_bg * cor_nhit * cor_sin * cor_t0 * sig_par[18];
471
472 }
473 }//end of MC
474 else {
475 if(sigmaflag == 1) {
476 f_betagamma=sig_par[0]*pow(beta_G,sig_par[1])+sig_par[2];
477 g_sinth=(sig_par[3]*sig_the*sig_the+sig_par[4])/(sig_par[3]*sig_par[5]*sig_par[5]+sig_par[4]);
478 h_nhit=(sig_par[6]*Nohit*Nohit+sig_par[7]*Nohit+sig_par[8]) /
479 (sig_par[6]*sig_par[9]*sig_par[9]+sig_par[7]*sig_par[9]+sig_par[8]);
480 if(sig_par[13] != 0)
481 i_t0 = (sig_par[10]*t0*t0+sig_par[11]*t0+sig_par[12]) /
482 (sig_par[10]*sig_par[13]*sig_par[13]+sig_par[11]*sig_par[13]+sig_par[12]);
483 else if(sig_par[13] == 0)
484 i_t0 =1;
485 //cout<<"f_betagamma : "<<f_betagamma<<" g_sinth : "<<g_sinth<<" h_nhit : "
486 //<<h_nhit<<" i_t0 : "<<i_t0<<endl;
487 ex_sigma[it]= f_betagamma* g_sinth * h_nhit * i_t0;
488 }
489 else if(sigmaflag == 2) {
490 double x = beta_G;
491 double nhit = (double)Nohit;
492 double sigma_bg=1.0;
493 if (x > 0.3)
494 sigma_bg = sig_par[0]*exp(sig_par[1]*x)+sig_par[2]*exp(sig_par[3]*pow(x,0.25))+sig_par[4];
495 else
496 sigma_bg = sig_par[5]*exp(sig_par[6]*x)+ sig_par[7];
497
498 double cor_nhit=1.0;
499 if(nhit<5) nhit=5;
500 if (nhit <= 35)
501 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];
502
503 double cor_sin= 1.0;
504 if(sig_the>0.99) sig_the = 0.99;
505 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
506 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
507
508 double cor_t0 = 1;
509 if (t0 > 1200) t0 = 1200;
510 if (t0 > 800)
511 cor_t0 = sig_par[18]*pow(t0,2)+sig_par[19]*t0+sig_par[20];
512
513 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
514 }
515 else if(sigmaflag == 3) {
516 double x= beta_G;
517 double nhit = (double)Nohit;
518 double sigma_bg=1.0;
519 if (x > 0.3)
520 sigma_bg = sig_par[0]*exp(sig_par[1]*x)+sig_par[2]*exp(sig_par[3]*pow(x,0.25))+sig_par[4];
521 else
522 sigma_bg = sig_par[5]*exp(sig_par[6]*x)+ sig_par[7];
523
524 double cor_nhit = 1.0;
525 if (nhit < 5) nhit = 5;
526 if (nhit <= 35)
527 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];
528
529 double cor_sin= 1.0;
530 if(sig_the>0.99) sig_the=0.99;
531 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
532 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
533
534 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
535 }
536 else if(sigmaflag == 4) {
537 double x= beta_G;
538 double nhit = (double)Nohit;
539 double sigma_bg=1.0;
540 if (x > 0.3)
541 sigma_bg = sig_par[0]*exp(sig_par[1]*x)+sig_par[2]*exp(sig_par[3]*pow(x,0.25))+sig_par[4];
542 else
543 sigma_bg = sig_par[5]*exp(sig_par[6]*x)+ sig_par[7];
544
545 double cor_nhit = 1.0;
546 if (nhit < 5) nhit = 5;
547 if (nhit <= 35)
548 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];
549
550 double cor_sin= 1.0;
551 if(sig_the>0.99) sig_the=0.99;
552 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
553 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
554
555 if(trkalg==1)
556 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
557 else
558 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*sig_par[18];
559 }
560 else if(sigmaflag == 5) {
561 double x = beta_G;
562 double nhit = (double)Nohit;
563 double sigma_bg=1.0;
564 if (x > 0.3)
565 sigma_bg = sig_par[0]*exp(sig_par[1]*x)+sig_par[2]*exp(sig_par[3]*pow(x,0.25))+sig_par[4];
566 else
567 sigma_bg = sig_par[5]*exp(sig_par[6]*x)+ sig_par[7];
568
569 double cor_nhit=1.0;
570 if(nhit<5) nhit=5;
571 if (nhit <= 35)
572 cor_nhit = sig_par[8]*pow(nhit,4)+sig_par[9]*pow(nhit,3)+sig_par[10]*pow(nhit,2) +
573 sig_par[11]*nhit+sig_par[12];
574 double cor_sin= 1.0;
575 if(sig_the>0.99) sig_the = 0.99;
576 cor_sin = sig_par[13]*pow(sig_the,4)+sig_par[14]*pow(sig_the,3) +
577 sig_par[15]*pow(sig_the,2)+sig_par[16]*sig_the+sig_par[17];
578
579 double cor_t0 = 1;
580 if (t0 > 1200) t0 = 1200;
581 if (t0 > 800)
582 cor_t0 = sig_par[18]*pow(t0,2)+sig_par[19]*t0+sig_par[20];
583
584 if(trkalg==1)
585 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
586 else
587 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0*sig_par[21];
588 }
589 else if(sigmaflag == 6) {
590 double x= bethe_B/550;
591 double nhit = (double)Nohit;
592 double sigma_bg=1.0;
593
594 sigma_bg= sig_par[0]+sig_par[1]*x+sig_par[2]*x*x;
595
596 double cor_nhit = 1.0;
597 if (nhit < 5) nhit = 5;
598 if (nhit <= 35)
599 cor_nhit = sig_par[3]*pow(nhit,4)+sig_par[4]*pow(nhit,3)+sig_par[5]*pow(nhit,2) + sig_par[6]*nhit+sig_par[7];
600
601 double cor_sin= 1.0;
602 if(sig_the>0.99) sig_the=0.99;
603 cor_sin = sig_par[8]*pow(sig_the,4)+sig_par[9]*pow(sig_the,3) +
604 sig_par[10]*pow(sig_the,2)+sig_par[11]*sig_the+sig_par[12];
605
606 if(trkalg==1)
607 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin;
608 else
609 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*sig_par[13];
610 }
611 else if(sigmaflag == 7) {
612 double x= bethe_B/550;
613 double nhit = (double)Nohit;
614 double sigma_bg=1.0;
615
616 sigma_bg= sig_par[0]+sig_par[1]*x+sig_par[2]*x*x;
617
618
619 double cor_nhit=1.0;
620 if(nhit<5) nhit=5;
621 if (nhit <= 35)
622 cor_nhit = sig_par[3]*pow(nhit,4)+sig_par[4]*pow(nhit,3)+sig_par[5]*pow(nhit,2) +
623 sig_par[6]*nhit+sig_par[7];
624 double cor_sin= 1.0;
625 if(sig_the>0.99) sig_the = 0.99;
626 cor_sin = sig_par[8]*pow(sig_the,4)+sig_par[9]*pow(sig_the,3) +
627 sig_par[10]*pow(sig_the,2)+sig_par[11]*sig_the+sig_par[12];
628
629 double cor_t0 = 1;
630 if (t0 > 1200) t0 = 1200;
631 if (t0 > 800)
632 cor_t0 = sig_par[13]*pow(t0,2)+sig_par[14]*t0+sig_par[15];
633
634 if(trkalg==1)
635 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0;
636 else
637 ex_sigma[it]=550*sigma_bg*cor_nhit*cor_sin*cor_t0*sig_par[16];
638 }
639 }
640
641 double dedx_correc = dedx;
642 chi_dedx[it] = (dedx_correc - dedx_exp[it])/ex_sigma[it];
643 chi2 = chi_dedx[it]*chi_dedx[it];
644 ndf=1;
645 pid_prob[it] = 0.0;//prob_(&chi2,&ndf);
646 //if(it ==0 ) cout<<"runflag: "<<runflag<<" dedx : "<<dedx<<" chi_dedx: "
647 //<<chi_dedx[it] <<" ptrk: "<<mom<<endl;
648 if (it == -999) { // here a debug flag
649 std::cout << " mom = " << mom <<"exp"<< dedx_exp[it]
650 << " sigma "<<ex_sigma[it] <<" prob "<<pid_prob[it]
651 << std::endl;
652 }
653 if( pid_prob[it] > max_prob ){
654 max_prob = pid_prob[it];
655 Nmax_prob = it;
656 }
657 }
658 else{
659 dedx_exp[it] = 0.0;
660 ex_sigma[it] = 1000.0;
661 pid_prob[it] = 0.0;
662 chi_dedx[it] = 999.0;
663 } //if Nohit > 0
664 }
665 // std::cout<<"MdcDedxCorrection::dedx_pid_exp(blum)!!!"<<std::endl;
666}
667
668#endif
sigma
TTree * sigma
Definition: CurConst_dEdXCal.cxx:6
curve
TTree * curve
Definition: CurConst_dEdXCal.cxx:5
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
s
XmlRpcServer s
Definition: HelloServer.cpp:11
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
curve_files
const std::string curve_files[2]
Definition: check/curve/global.h:28
vflag
int vflag[3]
Curve parameter vars.
Definition: histgen/curve.h:5
ReadCurPara
void ReadCurPara(int index)
Definition: histgen/curve.h:10
spar
std::vector< double > spar
Definition: histgen/curve.h:7
dedx_pid_exp
void dedx_pid_exp(double mom, double theta, double Nohit, double t0, double dedx_exp[5], double ex_sigma[5])
Definition: histgen/curve.h:45
cpar
std::vector< double > cpar
Definition: histgen/curve.h:6