BOSS 7.0.9
BESIII Offline Software System
Loading...
Searching...
No Matches
DedxPID.cxx
Go to the documentation of this file.
1#include <fstream>
2#include <cmath>
3#include <cstdlib>
4#include "Math/ChebyshevPol.h"
5#include "ParticleID/DedxPID.h"
6#ifndef BEAN
7#include "EvtRecEvent/EvtRecTrack.h"
8#include "MdcRecEvent/RecMdcTrack.h"
9#include "MdcRecEvent/RecMdcDedx.h"
10#else
11#include "RootEventDataVer.h"
12#include "DstEvtRecTracks.h"
13#endif
14
15DedxPID * DedxPID::m_pointer = 0;
16
17DedxPID * DedxPID::instance() {
18 if(!m_pointer) m_pointer = new DedxPID();
19 return m_pointer;
20}
21
22DedxPID::DedxPID():ParticleIDBase() {
23 m_readstate=0;
24}
25
26void DedxPID::init() {
27 for(int i = 0; i < 5; i++) {
28 m_chi[i] = 99.0;
29 m_prob[i] = -1.0;
30 m_offset[i] = 99.0;
31 m_sigma[i] = 1.0;
32 }
33 m_chimin = 99.;
34 m_pdfmin =99;
35 m_ndof = 0;
36 m_goodHits = -99;
37 m_normPH = -99;
38 m_probPH = -99;
39 m_nhitcutdx=5;
40}
41
42void DedxPID::calculate() {
43 // int rundedx = getRunNo();
44 if(!m_readstate) {
45 inputpar();
46 m_readstate=1;
47 }
48 if(particleIDCalculation() == 0) m_ndof=1;
49}
50
51int DedxPID::particleIDCalculation() {
52 // int rundedx2 = getRunNo();
53 int runNo = getRunNo();
54 int nhitcutdedx=getNhitCutDx();
55 int irc = -1;
56 EvtRecTrack* recTrk = PidTrk();
57 if(!(recTrk->isMdcTrackValid())) return irc;
58 RecMdcTrack* mdcTrk = recTrk->mdcTrack();
59
60 double ptrk = mdcTrk->p();
61 int charge = mdcTrk->charge();
62 if(ptrk>5) return irc;
63 double cost = cos(mdcTrk->theta());
64 // double sig_the= sin(mdcTrk->theta());
65
66 if(!(recTrk->isMdcDedxValid())) return irc;
67 RecMdcDedx* dedxTrk = recTrk->mdcDedx();
68
69 //if((dedxTrk->normPH()>30)||(dedxTrk->normPH()<0)) return irc;
70 m_goodHits = dedxTrk->numGoodHits();
71 //if(dedxTrk->numGoodHits()<nhitcutdedx) return irc;
72 m_normPH = dedxTrk->normPH();
73 m_probPH = dedxTrk->probPH();
74 // calculate chi and min chi
75 double chitemp = 99.;
76 double pdftemp = 0;
77 // double testchi[5];
78 // double testptrk[5];
79 // double testcost[5];
80 for(int i = 0; i < 5; i++) {
81 //add by CHEN Tong
82 if (recTrk->isMdcKalTrackValid())
83 {
84 RecMdcKalTrack *kalTrk = recTrk->mdcKalTrack();
85 if (i == 0)
86 {
87 kalTrk->setPidType(RecMdcKalTrack::electron);
88 }
89 else if (i == 1)
90 {
91 kalTrk->setPidType(RecMdcKalTrack::muon);
92 }
93 else if (i == 2)
94 {
95 kalTrk->setPidType(RecMdcKalTrack::pion);
96 }
97 else if (i == 3)
98 {
99 kalTrk->setPidType(RecMdcKalTrack::kaon);
100 }
101 else if (i == 4)
102 {
103 kalTrk->setPidType(RecMdcKalTrack::proton);
104 }
105 ptrk = kalTrk->p();
106 cost = cos(kalTrk->theta());
107 }
108 //CT
109 double sep = dedxTrk->chi(i);
110
111#ifndef BEAN
112 string sftver = getenv("BES_RELEASE");
113 string sft;
114 sft.assign(sftver,0,5);
115 if(sft=="6.6.0"||sft=="6.5.5") {
116 m_chi[i] = CorrDedx(i,ptrk,cost,sep,charge);
117 }
118 else
119 m_chi[i]=sep;
120#else
121// This is BEAN part:
122#if (ROOTEVENTDATA_VERSION_NUMERIC == ROOTEVENTDATA_VER(6,5,5) ||\
123 ROOTEVENTDATA_VERSION_NUMERIC == ROOTEVENTDATA_VER(6,6,0) )
124 m_chi[i] = CorrDedx(i,ptrk,cost,sep,charge);
125#else
126 m_chi[i]=sep;
127#endif
128#endif
129
130 m_offset[i] = offsetDedx(i, ptrk, cost);
131 m_sigma[i] = sigmaDedx(i, ptrk, cost);
132 //add by CHEN Tong
133 if (runNo >= 9947 && runNo <=10878 || runNo>= 27255 && runNo <= 28236 || runNo>= 52940 && runNo <= 54976 || runNo>= 55861 && runNo <= 56546 || runNo>= 56788 && runNo <= 59015)
134 {
135 m_offsetCorr[i] = offsetCorr(i, charge, ptrk, cost);
136 m_sigmaCorr[i] = sigmaCorr(i, charge, ptrk, cost);
137 m_chi[i] = (sep - m_offsetCorr[i]) / m_sigmaCorr[i];
138 }
139 //CT
140 if(fabs(m_chi[i]) < chitemp) chitemp = fabs(m_chi[i]);
141 double ppp = pdfCalculate(m_chi[i],1);
142 if(fabs(ppp) > pdftemp) pdftemp = fabs(ppp);
143
144 }
145 m_chimin = chitemp;
146 m_pdfmin = pdftemp;
147 if(m_chimin > chiMinCut()) return irc;
148 if(pdftemp < pdfCalculate(pdfMinSigmaCut(),1)) return irc;
149
150
151 // calculate prob
152
153 for(int i = 0; i < 5; i++)
154 m_prob[i] = probCalculate(m_chi[i]*m_chi[i], 1);
155
156 m_ndof = 1;
157 irc = 0;
158 return irc;
159}
160
161//
162// dE/dx Correction routines
163//
164
165
166
167double DedxPID::offsetDedx(int n, double ptrk, double cost) {
168 return 0;
169}
170
171double DedxPID::CorrDedx(int n, double ptrk, double cost,double chi,int charge) {
172 int rundedx2 = getRunNo();
173 double offset = 0.0;
174 double offsetp = 0.0;
175 double offsetc = 0.0;
176 double sigcos = 1;
177 double sigp = 1;
178 double chicor=chi;
179 // double gb = ptrk/xmass(n);
180
181 switch(n) {
182 case 0: { // Electron
183 break;
184 }
185
186 case 1: {// Muon
187 break;
188 }
189
190 case 2: {// Pion
191 // double ptemp = ptrk;
192 double costm = cost;
193 if(ptrk<0.1||ptrk>1) break;
194 int index = int((ptrk-0.1)/0.05);
195 if(index<=0) index=1;
196 if(index>=17) index=16;
197
198 if(fabs(costm)>=0.8) break;
199 int index1 = int((costm+0.8)/0.1);
200 if(index1<=0) index1=1;
201 if(index1>=15) index1=14;
202
203 //psipp data
204 if(rundedx2>=11414&&rundedx2<=14604) {
205 offsetp = cal_par(index,m_psipp_pi_ptrk_offset,ptrk,0.125,0.05);
206 sigp = cal_par(index,m_psipp_pi_ptrk_sigma,ptrk,0.125,0.05);
207 offsetc = cal_par(index1,m_psipp_pi_theta_offset,costm,-0.75,0.1);
208 sigcos = cal_par(index1,m_psipp_pi_theta_sigma,costm,-0.75,0.1);
209 }
210 //psipp mc
211 if(rundedx2<=-11414&&rundedx2>=-14604) {
212 offsetp = cal_par(index,m_psipp_mc_pi_ptrk_offset,ptrk,0.125,0.05);
213 sigp = cal_par(index,m_psipp_mc_pi_ptrk_sigma,ptrk,0.125,0.05);
214 offsetc = cal_par(index1,m_psipp_mc_pi_theta_offset,costm,-0.75,0.1);
215 sigcos = cal_par(index1,m_psipp_mc_pi_theta_sigma,costm,-0.75,0.1);
216 }
217
218 offset=offsetp+sigp*offsetc;
219 chicor=(chicor-offset)/(sigcos*sigp);
220 break;
221 }
222
223 case 3: {// Kaon
224 // double ptemp = ptrk;
225 double costm = cost;
226 if(ptrk<0.3||ptrk>0.8) break;
227 offset=0;
228 int index = int((ptrk-0.3)/0.1);
229 if(index<=0) index=1;
230 if(index>=4) index=3;
231
232 int index1 = int((costm+0.9)/0.1);
233 if(index1<=0) index1=1;
234 if(index1>=17) index1=16;
235 //data Jpsi
236 if(rundedx2>=9947&&rundedx2<=10878) {
237 if(charge>0) {
238 offsetp = cal_par(index,m_jpsi_kap_ptrk_offset,ptrk,0.35,0.1);
239 sigp = cal_par(index,m_jpsi_kap_ptrk_sigma,ptrk,0.35,0.1);
240 if(fabs(costm)<=0.83) {
241 offsetc = cal_par(index1,m_jpsi_kap_theta_offset,costm,-0.85,0.1);
242 sigcos = cal_par(index1,m_jpsi_kap_theta_sigma,costm,-0.85,0.1);
243 }
244 }
245 if(charge<0) {
246 offsetp = cal_par(index,m_jpsi_kam_ptrk_offset,ptrk,0.35,0.1);
247 sigp = cal_par(index,m_jpsi_kam_ptrk_sigma,ptrk,0.35,0.1);
248 if(fabs(costm)<=0.83) {
249 offsetc = cal_par(index1,m_jpsi_kam_theta_offset,costm,-0.85,0.1);
250 sigcos = cal_par(index1,m_jpsi_kam_theta_sigma,costm,-0.85,0.1);
251 }
252 }
253 }
254
255 //mc Jpsi
256 if(rundedx2<=-9947&&rundedx2>=-10878) {
257 if(charge>0) {
258 offsetp = cal_par(index,m_jpsi_mc_kap_ptrk_offset,ptrk,0.35,0.1);
259 sigp = cal_par(index,m_jpsi_mc_kap_ptrk_sigma,ptrk,0.35,0.1);
260 if(fabs(costm)<=0.83) {
261 offsetc = cal_par(index1,m_jpsi_mc_kap_theta_offset,costm,-0.85,0.1);
262 sigcos = cal_par(index1,m_jpsi_mc_kap_theta_sigma,costm,-0.85,0.1);
263 }
264 }
265 if(charge<0) {
266 offsetp = cal_par(index,m_jpsi_mc_kam_ptrk_offset,ptrk,0.35,0.1);
267 sigp = cal_par(index,m_jpsi_mc_kam_ptrk_sigma,ptrk,0.35,0.1);
268 if(fabs(costm)<=0.83) {
269 offsetc = cal_par(index1,m_jpsi_mc_kam_theta_offset,costm,-0.85,0.1);
270 sigcos = cal_par(index1,m_jpsi_mc_kam_theta_sigma,costm,-0.85,0.1);
271 }
272 }
273 }
274
275 //data Psip
276 if(rundedx2>=8093&&rundedx2<=9025) {
277 if(ptrk<0.3||ptrk>1.2) break;
278 index = int((ptrk-0.3)/0.1);
279 if(index<=0) index=1;
280 if(index>=8) index=7;
281 if(charge>0) {
282 offsetp = cal_par(index,m_psip_kap_ptrk_offset,ptrk,0.35,0.1);
283 sigp = cal_par(index,m_psip_kap_ptrk_sigma,ptrk,0.35,0.1);
284 }
285 if(charge<0) {
286 offsetp = cal_par(index,m_psip_kam_ptrk_offset,ptrk,0.35,0.1);
287 sigp = cal_par(index,m_psip_kam_ptrk_sigma,ptrk,0.35,0.1);
288 }
289 }
290
291 //mc Psip
292 if(rundedx2<=-8093&&rundedx2>=-9025) {
293 // if(ptrk < 0.4) ptrk = 0.4;
294 if(ptrk<0.3||ptrk>1.2) break;
295 index = int((ptrk-0.3)/0.1);
296 if(index<=0) index=1;
297 if(index>=8) index=7;
298 if(charge>0) {
299 offsetp = cal_par(index,m_psip_mc_kap_ptrk_offset,ptrk,0.35,0.1);
300 sigp = cal_par(index,m_psip_mc_kap_ptrk_sigma,ptrk,0.35,0.1);
301 }
302 if(charge<0) {
303 offsetp = cal_par(index,m_psip_mc_kam_ptrk_offset,ptrk,0.35,0.1);
304 sigp = cal_par(index,m_psip_mc_kam_ptrk_sigma,ptrk,0.35,0.1);
305 }
306 }
307
308
309 //psipp kaon data
310 if(rundedx2>=11414&&rundedx2<=14604) {
311 if(ptrk<0.15||ptrk>1) break;
312 index = int((ptrk-0.15)/0.05);
313 if(index<=0) index=1;
314 if(index>=16) index=15;
315 if(fabs(costm)>=0.8) break;
316 index1 = int((costm+0.8)/0.1);
317 if(index1<=0) index1=1;
318 if(index1>=15) index1=14;
319
320 offsetp = cal_par(index,m_psipp_ka_ptrk_offset,ptrk,0.175,0.05);
321 sigp = cal_par(index,m_psipp_ka_ptrk_sigma,ptrk,0.175,0.05);
322 offsetc = cal_par(index1,m_psipp_ka_theta_offset,costm,-0.75,0.1);
323 sigcos = cal_par(index1,m_psipp_ka_theta_sigma,costm,-0.75,0.1);
324 }
325 //psipp kaon mc
326 if(rundedx2<=-11414&&rundedx2>=-14604) {
327 if(ptrk<0.15||ptrk>1) break;
328 index = int((ptrk-0.15)/0.05);
329 if(index<=0) index=1;
330 if(index>=16) index=15;
331 if(fabs(costm)>=0.8) break;
332 index1 = int((costm+0.8)/0.1);
333 if(index1<=0) index1=1;
334 if(index1>=15) index1=14;
335 offsetp = cal_par(index,m_psipp_mc_ka_ptrk_offset,ptrk,0.175,0.05);
336 sigp = cal_par(index,m_psipp_mc_ka_ptrk_sigma,ptrk,0.175,0.05);
337 offsetc = cal_par(index1,m_psipp_mc_ka_theta_offset,costm,-0.75,0.1);
338 sigcos = cal_par(index1,m_psipp_mc_ka_theta_sigma,costm,-0.75,0.1);
339 }
340
341 offset=offsetp+sigp*offsetc;
342 chicor=(chicor-offset)/(sigcos*sigp);
343 break;
344 }
345
346 case 4 : { // Proton
347 // double ptemp = ptrk;
348 double costm = cost;
349 if(ptrk<0.3||ptrk>1.1) break;
350 int index = int((ptrk-0.3)/0.1);
351 if(index<=0) index=1;
352 if(index>=7) index=6;
353
354 int index1 = int((costm+0.9)/0.1);
355 if(index1<=0) index1=1;
356 if(index1>=17) index1=16;
357
358 // double plog = log(ptemp);
359 offset=0;
360 if(rundedx2>=9947&&rundedx2<=10878) {
361 if(charge>0) {
362 offsetp = cal_par(index,m_jpsi_protonp_ptrk_offset,ptrk,0.35,0.1);
363 sigp = cal_par(index,m_jpsi_protonp_ptrk_sigma,ptrk,0.35,0.1);
364 if(fabs(costm)<=0.83) {
365 offsetc = cal_par(index1,m_jpsi_protonp_theta_offset,costm,-0.85,0.1);
366 sigcos = cal_par(index1,m_jpsi_protonp_theta_sigma,costm,-0.85,0.1);
367 }
368 }
369 if(charge<0) {
370 offsetp = cal_par(index,m_jpsi_protonm_ptrk_offset,ptrk,0.35,0.1);
371 sigp = cal_par(index,m_jpsi_protonm_ptrk_sigma,ptrk,0.35,0.1);
372 if(fabs(costm)<=0.83) {
373 offsetc = cal_par(index1,m_jpsi_protonm_theta_offset,costm,-0.85,0.1);
374 sigcos = cal_par(index1,m_jpsi_protonm_theta_sigma,costm,-0.85,0.1);
375 }
376 }
377 }
378
379 //mc JPsi
380 if(rundedx2<=-9947&&rundedx2>=-10878) {
381 if(charge>0) {
382 offsetp = cal_par(index,m_jpsi_mc_protonp_ptrk_offset,ptrk,0.35,0.1);
383 sigp = cal_par(index,m_jpsi_mc_protonp_ptrk_sigma,ptrk,0.35,0.1);
384 if(fabs(costm)<=0.83) {
385 offsetc = cal_par(index1,m_jpsi_mc_protonp_theta_offset,costm,-0.85,0.1);
386 sigcos = cal_par(index1,m_jpsi_mc_protonp_theta_sigma,costm,-0.85,0.1);
387 }
388 }
389 if(charge<0) {
390 offsetp = cal_par(index,m_jpsi_mc_protonm_ptrk_offset,ptrk,0.35,0.1);
391 sigp = cal_par(index,m_jpsi_mc_protonm_ptrk_sigma,ptrk,0.35,0.1);
392 if(fabs(costm)<=0.83) {
393 offsetc = cal_par(index1,m_jpsi_mc_protonm_theta_offset,costm,-0.85,0.1);
394 sigcos = cal_par(index1,m_jpsi_mc_protonm_theta_sigma,costm,-0.85,0.1);
395 }
396 }
397 }
398
399 //data Psip
400 if(rundedx2>=8093&&rundedx2<=9025) {
401 if(charge>0) {
402 offsetp = cal_par(index,m_psip_protonp_ptrk_offset,ptrk,0.35,0.1);
403 sigp = cal_par(index,m_psip_protonp_ptrk_sigma,ptrk,0.35,0.1);
404 }
405 if(charge<0) {
406 offsetp = cal_par(index,m_psip_protonm_ptrk_offset,ptrk,0.35,0.1);
407 sigp = cal_par(index,m_psip_protonm_ptrk_sigma,ptrk,0.35,0.1);
408 }
409 }
410
411 //mc Psip
412 if(rundedx2<=-8093&&rundedx2>=-9025) {
413 if(charge>0) {
414 offsetp = cal_par(index,m_psip_mc_protonp_ptrk_offset,ptrk,0.35,0.1);
415 sigp = cal_par(index,m_psip_mc_protonp_ptrk_sigma,ptrk,0.35,0.1);
416 }
417 if(charge<0) {
418 offsetp = cal_par(index,m_psip_mc_protonm_ptrk_offset,ptrk,0.35,0.1);
419 sigp = cal_par(index,m_psip_mc_protonm_ptrk_sigma,ptrk,0.35,0.1);
420 }
421 }
422
423 //psipp proton data
424 if(rundedx2>=11414&&rundedx2<=14604) {
425 if(ptrk<0.2||ptrk>1.1) break;
426 index = int((ptrk-0.2)/0.05);
427 if(index<=0) index=1;
428 if(index>=17) index=16;
429 if(fabs(costm)>=0.83) break;
430 index1 = int((costm+0.9)/0.1);
431 if(index1<=0) index1=1;
432 if(index1>=17) index1=16;
433
434 offsetp = cal_par(index,m_psipp_proton_ptrk_offset,ptrk,0.225,0.05);
435 sigp = cal_par(index,m_psipp_proton_ptrk_sigma,ptrk,0.225,0.05);
436 offsetc = cal_par(index1,m_psipp_proton_theta_offset,costm,-0.85,0.1);
437 sigcos = cal_par(index1,m_psipp_proton_theta_sigma,costm,-0.85,0.1);
438 }
439 //psipp proton mc
440 if(rundedx2<=-11414&&rundedx2>=-14604) {
441 if(ptrk<0.2||ptrk>1.1) break;
442 index = int((ptrk-0.2)/0.1);
443 if(index<=0) index=1;
444 if(index>=8) index=7;
445 if(fabs(costm)>=0.83) break;
446 index1 = int((costm+0.9)/0.1);
447 if(index1<=0) index1=1;
448 if(index1>=17) index1=16;
449 offsetp = cal_par(index,m_psipp_mc_proton_ptrk_offset,ptrk,0.25,0.1);
450 sigp = cal_par(index,m_psipp_mc_proton_ptrk_sigma,ptrk,0.25,0.1);
451 offsetc = cal_par(index1,m_psipp_mc_proton_theta_offset,costm,-0.85,0.1);
452 sigcos = cal_par(index1,m_psipp_mc_proton_theta_sigma,costm,-0.85,0.1);
453 }
454 offset=offsetp+sigp*offsetc;
455 chicor=(chicor-offset)/(sigcos*sigp);
456 break;
457 }
458
459 default:
460 offset = 0.0;
461 break;
462 }
463 // offset = 0.0;
464 return chicor;
465}
466
467double DedxPID::sigmaDedx(int n, double ptrk, double cost) {
468
469 /* int rundedx3 = getRunNo();
470 double sigma = 1.0;
471 double sigp = 1.0;
472 double sigmac = 1.0;
473 double gb = ptrk/xmass(n);
474 switch(n) {
475
476 case 0: {// Electron
477 double ptemp = ptrk;
478 double costm = cost;
479 break;
480 }
481
482 case 1: {// Muon
483 double ptemp = ptrk;
484 double costm = cost;
485 break;
486 }
487
488 case 2: {// Pion
489 double ptemp = ptrk;
490 double costm = cost;
491 break;
492 }
493
494 case 3: { // Kaon
495 double ptemp = ptrk;
496 double costm = cost;
497 break;
498 }
499
500
501 case 4: {// Proton
502 double ptemp = ptrk;
503 double costm = cost;
504 break;
505 }
506
507 default:
508 sigma = 1.0;
509 break;
510 }
511 */
512 // sigma = 1.2;
513 // sigma =1.0;
514 return 1;
515 // return sigma;
516}
517
518double DedxPID::mypol3(double x, double par0, double par1, double par2, double par3)
519{
520 double y = x;
521 return par0 + (par1 * y) +(par2 * y * y) + (par3 * y * y * y);
522
523}
524
525double DedxPID::mypol5(double x, double par0, double par1, double par2, double par3, double par4, double par5)
526{
527 double y = x;
528 return par0 + (par1 * y) +(par2 * y * y) + (par3 * y * y * y) + (par4 * y * y * y *y)+ (par5 * y * y * y * y * y);
529
530}
531
532void DedxPID::inputpar() {
533
534 //Jpsi ka+ momentum correction
535 std::string jpsi_kap_mom = path + "/share/JPsi/kaon/dedx_kap.txt";
536 std::string jpsi_kap_mom_mc = path + "/share/JPsi/kaon/dedx_kap_mc.txt";
537 ifstream inputmomdata6(jpsi_kap_mom.c_str(),std::ios_base::in);
538 if ( !inputmomdata6 ) {
539 cout << " can not open: " << jpsi_kap_mom << endl;
540 exit(1);
541 }
542 ifstream inputmomdata6mc(jpsi_kap_mom_mc.c_str(),std::ios_base::in);
543 if ( !inputmomdata6mc ) {
544 cout << " can not open: " << jpsi_kap_mom_mc << endl;
545 exit(1);
546 }
547 for(int i=0; i<12; i++) {
548 inputmomdata6>>m_jpsi_kap_ptrk_offset[i];
549 inputmomdata6>>m_jpsi_kap_ptrk_sigma[i];
550 inputmomdata6mc>>m_jpsi_mc_kap_ptrk_offset[i];
551 inputmomdata6mc>>m_jpsi_mc_kap_ptrk_sigma[i];
552 }
553
554 //Jpsi ka- momentum correction
555 std::string jpsi_kam_mom = path + "/share/JPsi/kaon/dedx_kam.txt";
556 std::string jpsi_kam_mom_mc = path + "/share/JPsi/kaon/dedx_kam_mc.txt";
557 ifstream inputmomdata7(jpsi_kam_mom.c_str(),std::ios_base::in);
558 if ( !inputmomdata7 ) {
559 cout << " can not open: " << jpsi_kam_mom << endl;
560 exit(1);
561 }
562 ifstream inputmomdata7mc(jpsi_kam_mom_mc.c_str(),std::ios_base::in);
563 if ( !inputmomdata7mc ) {
564 cout << " can not open: " << jpsi_kam_mom_mc << endl;
565 exit(1);
566 }
567 for(int i=0; i<12; i++) {
568 inputmomdata7>>m_jpsi_kam_ptrk_offset[i];
569 inputmomdata7>>m_jpsi_kam_ptrk_sigma[i];
570 inputmomdata7mc>>m_jpsi_mc_kam_ptrk_offset[i];
571 inputmomdata7mc>>m_jpsi_mc_kam_ptrk_sigma[i];
572
573 }
574
575
576 //Jpsi ka+ theta correction
577 std::string jpsi_kap_the = path + "/share/JPsi/kaon/dedx_kap_theta.txt";
578 std::string jpsi_kap_the_mc = path + "/share/JPsi/kaon/dedx_kap_theta_mc.txt";
579 ifstream inputmomdata8(jpsi_kap_the.c_str(),std::ios_base::in);
580 if ( !inputmomdata8 ) {
581 cout << " can not open: " << jpsi_kap_the << endl;
582 exit(1);
583 }
584 ifstream inputmomdata8mc(jpsi_kap_the_mc.c_str(),std::ios_base::in);
585 if ( !inputmomdata8mc ) {
586 cout << " can not open: " << jpsi_kap_the_mc << endl;
587 exit(1);
588 }
589 for(int i=0; i<18; i++) {
590 inputmomdata8>>m_jpsi_kap_theta_offset[i];
591 inputmomdata8>>m_jpsi_kap_theta_sigma[i];
592 inputmomdata8mc>>m_jpsi_mc_kap_theta_offset[i];
593 inputmomdata8mc>>m_jpsi_mc_kap_theta_sigma[i];
594 }
595
596 //Jpsi ka- theta correction
597 std::string jpsi_kam_the = path + "/share/JPsi/kaon/dedx_kam_theta.txt";
598 std::string jpsi_kam_the_mc = path + "/share/JPsi/kaon/dedx_kam_theta_mc.txt";
599 ifstream inputmomdata9(jpsi_kam_the.c_str(),std::ios_base::in);
600 if ( !inputmomdata9 ) {
601 cout << " can not open: " << jpsi_kam_the << endl;
602 exit(1);
603 }
604 ifstream inputmomdata9mc(jpsi_kam_the_mc.c_str(),std::ios_base::in);
605 if ( !inputmomdata9mc ) {
606 cout << " can not open: " << jpsi_kam_the_mc << endl;
607 exit(1);
608 }
609 for(int i=0; i<18; i++) {
610 inputmomdata9>>m_jpsi_kam_theta_offset[i];
611 inputmomdata9>>m_jpsi_kam_theta_sigma[i];
612 inputmomdata9mc>>m_jpsi_mc_kam_theta_offset[i];
613 inputmomdata9mc>>m_jpsi_mc_kam_theta_sigma[i];
614 }
615
616 //Jpsi proton+ momentum correction
617 std::string jpsi_protonp_mom = path + "/share/JPsi/proton/dedx_protonp.txt";
618 std::string jpsi_protonp_mom_mc = path + "/share/JPsi/proton/dedx_protonp_mc.txt";
619 ifstream inputmomdata12(jpsi_protonp_mom.c_str(),std::ios_base::in);
620 if ( !inputmomdata12 ) {
621 cout << " can not open: " << jpsi_protonp_mom << endl;
622 exit(1);
623 }
624 ifstream inputmomdata12mc(jpsi_protonp_mom_mc.c_str(),std::ios_base::in);
625 if ( !inputmomdata12mc ) {
626 cout << " can not open: " << jpsi_protonp_mom_mc << endl;
627 exit(1);
628 }
629 for(int i=0; i<8; i++) {
630 inputmomdata12>>m_jpsi_protonp_ptrk_offset[i];
631 inputmomdata12>>m_jpsi_protonp_ptrk_sigma[i];
632 inputmomdata12mc>>m_jpsi_mc_protonp_ptrk_offset[i];
633 inputmomdata12mc>>m_jpsi_mc_protonp_ptrk_sigma[i];
634 }
635
636 //Jpsi proton- momentum correction
637 std::string jpsi_protonm_mom = path + "/share/JPsi/proton/dedx_protonm.txt";
638 std::string jpsi_protonm_mom_mc = path + "/share/JPsi/proton/dedx_protonm_mc.txt";
639 ifstream inputmomdata13(jpsi_protonm_mom.c_str(),std::ios_base::in);
640 if ( !inputmomdata13 ) {
641 cout << " can not open: " << jpsi_protonm_mom << endl;
642 exit(1);
643 }
644 ifstream inputmomdata13mc(jpsi_protonm_mom_mc.c_str(),std::ios_base::in);
645 if ( !inputmomdata13mc ) {
646 cout << " can not open: " << jpsi_protonm_mom_mc << endl;
647 exit(1);
648 }
649 for(int i=0; i<8; i++) {
650 inputmomdata13>>m_jpsi_protonm_ptrk_offset[i];
651 inputmomdata13>>m_jpsi_protonm_ptrk_sigma[i];
652 inputmomdata13mc>>m_jpsi_mc_protonm_ptrk_offset[i];
653 inputmomdata13mc>>m_jpsi_mc_protonm_ptrk_sigma[i];
654 }
655
656 //Jpsi proton+ theta correction
657 std::string jpsi_protonp_the = path + "/share/JPsi/proton/dedx_protonp_theta.txt";
658 std::string jpsi_protonp_the_mc = path + "/share/JPsi/proton/dedx_protonp_theta_mc.txt";
659
660 ifstream inputmomdata14(jpsi_protonp_the.c_str(),std::ios_base::in);
661 if ( !inputmomdata14 ) {
662 cout << " can not open: " << jpsi_protonp_the << endl;
663 exit(1);
664 }
665 ifstream inputmomdata14mc(jpsi_protonp_the_mc.c_str(),std::ios_base::in);
666 if ( !inputmomdata14mc ) {
667 cout << " can not open: " << jpsi_protonp_the_mc << endl;
668 exit(1);
669 }
670 for(int i=0; i<18; i++) {
671 inputmomdata14>>m_jpsi_protonp_theta_offset[i];
672 inputmomdata14>>m_jpsi_protonp_theta_sigma[i];
673 inputmomdata14mc>>m_jpsi_mc_protonp_theta_offset[i];
674 inputmomdata14mc>>m_jpsi_mc_protonp_theta_sigma[i];
675 }
676
677 //Jpsi proton- theta correction
678 std::string jpsi_protonm_the = path + "/share/JPsi/proton/dedx_protonm_theta.txt";
679 std::string jpsi_protonm_the_mc = path + "/share/JPsi/proton/dedx_protonm_theta_mc.txt";
680 ifstream inputmomdata15(jpsi_protonm_the.c_str(),std::ios_base::in);
681 if ( !inputmomdata15 ) {
682 cout << " can not open: " << jpsi_protonm_the << endl;
683 exit(1);
684 }
685 ifstream inputmomdata15mc(jpsi_protonm_the_mc.c_str(),std::ios_base::in);
686 if ( !inputmomdata15mc ) {
687 cout << " can not open: " << jpsi_protonm_the_mc << endl;
688 exit(1);
689 }
690 for(int i=0; i<18; i++) {
691 inputmomdata15>>m_jpsi_protonm_theta_offset[i];
692 inputmomdata15>>m_jpsi_protonm_theta_sigma[i];
693 inputmomdata15mc>>m_jpsi_mc_protonm_theta_offset[i];
694 inputmomdata15mc>>m_jpsi_mc_protonm_theta_sigma[i];
695 }
696
697
698
699
700 // Psip ka+ momentum correction
701 std::string psip_kap_mom = path + "/share/Psip/kaon/dedx_kap.txt";
702 std::string psip_kap_mom_mc = path + "/share/Psip/kaon/dedx_kap_mc.txt";
703 ifstream inputmomdata24(psip_kap_mom.c_str(),std::ios_base::in);
704 if ( !inputmomdata24 ) {
705 cout << " can not open: " << psip_kap_mom << endl;
706 exit(1);
707 }
708 ifstream inputmomdata24mc(psip_kap_mom_mc.c_str(),std::ios_base::in);
709 if ( !inputmomdata24mc ) {
710 cout << " can not open: " << psip_kap_mom_mc << endl;
711 exit(1);
712 }
713 for(int i=0; i<9; i++) {
714 inputmomdata24>>m_psip_kap_ptrk_offset[i];
715 inputmomdata24>>m_psip_kap_ptrk_sigma[i];
716 inputmomdata24mc>>m_psip_mc_kap_ptrk_offset[i];
717 inputmomdata24mc>>m_psip_mc_kap_ptrk_sigma[i];
718 }
719
720 //Psip ka- momentum correction
721 std::string psip_kam_mom = path + "/share/Psip/kaon/dedx_kam.txt";
722 std::string psip_kam_mom_mc = path + "/share/Psip/kaon/dedx_kam_mc.txt";
723 ifstream inputmomdata25(psip_kam_mom.c_str(),std::ios_base::in);
724 if ( !inputmomdata25 ) {
725 cout << " can not open: " << psip_kam_mom << endl;
726 exit(1);
727 }
728 ifstream inputmomdata25mc(psip_kam_mom_mc.c_str(),std::ios_base::in);
729 if ( !inputmomdata25mc ) {
730 cout << " can not open: " << psip_kam_mom_mc << endl;
731 exit(1);
732 }
733 for(int i=0; i<9; i++) {
734 inputmomdata25>>m_psip_kam_ptrk_offset[i];
735 inputmomdata25>>m_psip_kam_ptrk_sigma[i];
736 inputmomdata25mc>>m_psip_mc_kam_ptrk_offset[i];
737 inputmomdata25mc>>m_psip_mc_kam_ptrk_sigma[i];
738 }
739
740
741 // Psip proton+ momentum correction
742 std::string psip_protonp_mom = path + "/share/Psip/proton/dedx_protonp.txt";
743 std::string psip_protonp_mom_mc = path + "/share/Psip/proton/dedx_protonp_mc.txt";
744 ifstream inputmomdata26(psip_protonp_mom.c_str(),std::ios_base::in);
745 if ( !inputmomdata26 ) {
746 cout << " can not open: " << psip_protonp_mom << endl;
747 exit(1);
748 }
749 ifstream inputmomdata26mc(psip_protonp_mom_mc.c_str(),std::ios_base::in);
750 if ( !inputmomdata26mc ) {
751 cout << " can not open: " << psip_protonp_mom_mc << endl;
752 exit(1);
753 }
754 for(int i=0; i<9; i++) {
755 inputmomdata26>>m_psip_protonp_ptrk_offset[i];
756 inputmomdata26>>m_psip_protonp_ptrk_sigma[i];
757 inputmomdata26mc>>m_psip_mc_protonp_ptrk_offset[i];
758 inputmomdata26mc>>m_psip_mc_protonp_ptrk_sigma[i];
759 }
760
761 //Psip proton- momentum correction
762 std::string psip_protonm_mom = path + "/share/Psip/proton/dedx_protonm.txt";
763 std::string psip_protonm_mom_mc = path + "/share/Psip/proton/dedx_protonm_mc.txt";
764 ifstream inputmomdata27(psip_protonm_mom.c_str(),std::ios_base::in);
765 if ( !inputmomdata27 ) {
766 cout << " can not open: " << psip_protonm_mom << endl;
767 exit(1);
768 }
769 ifstream inputmomdata27mc(psip_protonm_mom_mc.c_str(),std::ios_base::in);
770 if ( !inputmomdata27mc ) {
771 cout << " can not open: " << psip_protonm_mom_mc << endl;
772 exit(1);
773 }
774 for(int i=0; i<9; i++) {
775 inputmomdata27>>m_psip_protonm_ptrk_offset[i];
776 inputmomdata27>>m_psip_protonm_ptrk_sigma[i];
777 inputmomdata27mc>>m_psip_mc_protonm_ptrk_offset[i];
778 inputmomdata27mc>>m_psip_mc_protonm_ptrk_sigma[i];
779 }
780
781 //Psipp pi momentum correction
782 std::string psipp_pi_mom = path + "/share/Psipp/pion/dedx_pi.txt";
783 std::string psipp_pi_mom_mc = path + "/share/Psipp/pion/dedx_pi_mc.txt";
784 ifstream inputmomdata28(psipp_pi_mom.c_str(),std::ios_base::in);
785 if ( !inputmomdata28 ) {
786 cout << " can not open: " << psipp_pi_mom << endl;
787 exit(1);
788 }
789 ifstream inputmomdata28mc(psipp_pi_mom_mc.c_str(),std::ios_base::in);
790 if ( !inputmomdata28mc ) {
791 cout << " can not open: " << psipp_pi_mom_mc << endl;
792 exit(1);
793 }
794 for(int i=0; i<18; i++) {
795 inputmomdata28>>m_psipp_pi_ptrk_offset[i];
796 inputmomdata28>>m_psipp_pi_ptrk_sigma[i];
797 inputmomdata28mc>>m_psipp_mc_pi_ptrk_offset[i];
798 inputmomdata28mc>>m_psipp_mc_pi_ptrk_sigma[i];
799 }
800
801 //Psipp pi theta correction
802 std::string psipp_pi_the = path + "/share/Psipp/pion/dedx_pi_theta.txt";
803 std::string psipp_pi_the_mc = path + "/share/Psipp/pion/dedx_pi_theta_mc.txt";
804 ifstream inputmomdata29(psipp_pi_the.c_str(),std::ios_base::in);
805 if ( !inputmomdata29 ) {
806 cout << " can not open: " << psipp_pi_the << endl;
807 exit(1);
808 }
809 ifstream inputmomdata29mc(psipp_pi_the_mc.c_str(),std::ios_base::in);
810 if ( !inputmomdata29mc ) {
811 cout << " can not open: " << psipp_pi_the_mc << endl;
812 exit(1);
813 }
814 for(int i=0; i<16; i++) {
815 inputmomdata29>>m_psipp_pi_theta_offset[i];
816 inputmomdata29>>m_psipp_pi_theta_sigma[i];
817 inputmomdata29mc>>m_psipp_mc_pi_theta_offset[i];
818 inputmomdata29mc>>m_psipp_mc_pi_theta_sigma[i];
819 }
820
821 //Psipp ka momentum correction
822 std::string psipp_ka_mom = path + "/share/Psipp/kaon/dedx_ka.txt";
823 std::string psipp_ka_mom_mc = path + "/share/Psipp/kaon/dedx_ka_mc.txt";
824 ifstream inputmomdata30(psipp_ka_mom.c_str(),std::ios_base::in);
825 if ( !inputmomdata30 ) {
826 cout << " can not open: " << psipp_ka_mom << endl;
827 exit(1);
828 }
829 ifstream inputmomdata30mc(psipp_ka_mom_mc.c_str(),std::ios_base::in);
830 if ( !inputmomdata30mc ) {
831 cout << " can not open: " << psipp_ka_mom_mc << endl;
832 exit(1);
833 }
834 for(int i=0; i<17; i++) {
835 inputmomdata30>>m_psipp_ka_ptrk_offset[i];
836 inputmomdata30>>m_psipp_ka_ptrk_sigma[i];
837 inputmomdata30mc>>m_psipp_mc_ka_ptrk_offset[i];
838 inputmomdata30mc>>m_psipp_mc_ka_ptrk_sigma[i];
839 }
840
841 //Psipp ka theta correction
842 std::string psipp_ka_the = path + "/share/Psipp/kaon/dedx_ka_theta.txt";
843 std::string psipp_ka_the_mc = path + "/share/Psipp/kaon/dedx_ka_theta_mc.txt";
844 ifstream inputmomdata31(psipp_ka_the.c_str(),std::ios_base::in);
845 if ( !inputmomdata31 ) {
846 cout << " can not open: " << psipp_ka_the << endl;
847 exit(1);
848 }
849 ifstream inputmomdata31mc(psipp_ka_the_mc.c_str(),std::ios_base::in);
850 if ( !inputmomdata31mc ) {
851 cout << " can not open: " << psipp_ka_the_mc << endl;
852 exit(1);
853 }
854 for(int i=0; i<16; i++) {
855 inputmomdata31>>m_psipp_ka_theta_offset[i];
856 inputmomdata31>>m_psipp_ka_theta_sigma[i];
857 inputmomdata31mc>>m_psipp_mc_ka_theta_offset[i];
858 inputmomdata31mc>>m_psipp_mc_ka_theta_sigma[i];
859 }
860
861
862 //Psipp proton momentum correction
863 std::string psipp_proton_mom = path + "/share/Psipp/proton/dedx_proton.txt";
864 std::string psipp_proton_mom_mc = path + "/share/Psipp/proton/dedx_proton_mc.txt";
865 ifstream inputmomdata32(psipp_proton_mom.c_str(),std::ios_base::in);
866 if ( !inputmomdata32 ) {
867 cout << " can not open: " << psipp_proton_mom << endl;
868 exit(1);
869 }
870 ifstream inputmomdata32mc(psipp_proton_mom_mc.c_str(),std::ios_base::in);
871 if ( !inputmomdata32mc ) {
872 cout << " can not open: " << psipp_proton_mom_mc << endl;
873 exit(1);
874 }
875 for(int i=0; i<18; i++) {
876 inputmomdata32>>m_psipp_proton_ptrk_offset[i];
877 inputmomdata32>>m_psipp_proton_ptrk_sigma[i];
878 }
879 for(int i=0; i<9; i++) {
880 inputmomdata32mc>>m_psipp_mc_proton_ptrk_offset[i];
881 inputmomdata32mc>>m_psipp_mc_proton_ptrk_sigma[i];
882 }
883
884 //Psipp proton theta correction
885 std::string psipp_proton_the = path + "/share/Psipp/proton/dedx_proton_theta.txt";
886 std::string psipp_proton_the_mc = path + "/share/Psipp/proton/dedx_proton_theta_mc.txt";
887 ifstream inputmomdata33(psipp_proton_the.c_str(),std::ios_base::in);
888 if ( !inputmomdata33 ) {
889 cout << " can not open: " << psipp_proton_the << endl;
890 exit(1);
891 }
892 ifstream inputmomdata33mc(psipp_proton_the_mc.c_str(),std::ios_base::in);
893 if ( !inputmomdata33mc ) {
894 cout << " can not open: " << psipp_proton_the_mc << endl;
895 exit(1);
896 }
897 for(int i=0; i<18; i++) {
898 inputmomdata33>>m_psipp_proton_theta_offset[i];
899 inputmomdata33>>m_psipp_proton_theta_sigma[i];
900 inputmomdata33mc>>m_psipp_mc_proton_theta_offset[i];
901 inputmomdata33mc>>m_psipp_mc_proton_theta_sigma[i];
902 }
903
904}
905
906double DedxPID::iterate(double ptrk,double *mean,double *p) {
907 double p1,p2,p3;
908 p2=((mean[0]-mean[1])*(p[1]*p[1]-p[2]*p[2])-(mean[1]-mean[2])*(p[0]*p[0]-p[1]*p[1]))/((p[0]-p[1])*(p[1]*p[1]-p[2]*p[2])-(p[1]-p[2])*(p[0]*p[0]-p[1]*p[1]));
909 p3=((p[0]-p[1])*(mean[1]-mean[2])-(p[1]-p[2])*(mean[0]-mean[1]))/((p[0]-p[1])*(p[1]*p[1]-p[2]*p[2])-(p[1]-p[2])*(p[0]*p[0]-p[1]*p[1]));
910 p1=mean[0]-p2*p[0]-p3*p[0]*p[0];
911 double mean1 = p1+p2*ptrk+p3*ptrk*ptrk;
912 return mean1;
913}
914
915double DedxPID::cal_par(int index1,double *m_jpsi_pip_ptrk_offset,double ptrk,double begin,double bin) {
916 double mean1[3],p[3];
917 p[0]=begin+(index1-1)*bin;
918 p[1]=begin+index1*bin;
919 p[2]=begin+(index1+1)*bin;
920 mean1[0]=m_jpsi_pip_ptrk_offset[index1-1];
921 mean1[1]=m_jpsi_pip_ptrk_offset[index1];
922 mean1[2]=m_jpsi_pip_ptrk_offset[index1+1];
923 double res=iterate(ptrk,mean1,p);
924 return res;
925}
926
927//add by CHEN Tong
928double DedxPID::offsetCorr(int n, int charge, double ptrk, double cost)
929{
930 double cosbin[30] = {-0.865, -0.7, -0.5, -0.325, -0.225, -0.19, -0.17, -0.15, -0.13, -0.11, -0.09, -0.07, -0.05, -0.03, -0.01, 0.01, 0.03, 0.05, 0.07, 0.09, 0.11, 0.13, 0.15, 0.17, 0.19, 0.225, 0.325, 0.5, 0.7, 0.865};
931 double corr_offset = 0;
932 if (n == 0 || n == 1 || n == 2 && ptrk > 0.20 || n == 3 && (ptrk > 0.43 || ptrk < 0.2 || fabs(cost) > 0.2) || n == 4 && (ptrk > 0.43 || ptrk < 0.2))
933 return 0;
934 else if (n == 2)
935 {
936 int tmp = 0;
937 double par_pip_offset[30][7] = {
938 {-0.273711, 0.418852, -0.111059, 0.011915, 0.0595808, -0.064644, -0.00930414},
939 {-0.269974, 0.295548, -0.0805116, 0.112013, -0.0844919, -0.0101669, 0.0472545},
940 {-0.248076, 0.230898, -0.0954416, 0.00461465, -0.0513009, 0.0841044, -0.0320536},
941 {-0.20576, 0.15908, -0.0827095, -0.0534723, 0.011227, 0.0325395, -0.0518127},
942 {-0.168736, 0.129392, -0.05107, -0.126828, 0.126172, -0.0489699, -0.0373492},
943 {-0.176773, 0.163439, -0.100698, -0.111413, 0.115099, -0.0726195, -0.00490596},
944 {-0.0385498, -0.0807258, 0.0897135, -0.237934, 0.243725, -0.108172, -0.0188745},
945 {-0.109386, 0.119386, -0.154018, -0.119937, 0.110189, -0.0662777, -0.0217286},
946 {-0.0942173, 0.171144, -0.0806635, -0.182887, 0.22318, -0.137083, 0.0174856},
947 {-0.0209689, 0.0959425, -0.00407716, -0.186367, 0.309657, -0.158974, 0.0899474},
948 {-0.0742976, 0.194453, -0.245653, -0.119011, 0.169814, -0.0985074, -0.0145252},
949 {-0.00874818, 0.381593, -0.398284, 0.133183, 0.0303175, -0.024684, -0.0535397},
950 {0.264664, 0.128337, -0.322216, 0.00757683, 0.072456, -0.0374562, -0.0347229},
951 {0.756099, -0.426524, 0.0622856, -0.194491, 0.148391, -0.027508, -0.0244935},
952 {1.01453, -0.255592, -0.28767, 0.285283, -0.249035, 0.257705, -0.188833},
953 {0.992926, -0.470497, -0.193231, 0.07265, -0.14322, 0.124669, -0.147956},
954 {0.461868, 0.0132674, -0.211087, 0.057901, 0.0263037, 0.0836409, -0.108405},
955 {0.1606, 0.0849392, -0.189776, -0.0816648, 0.0954246, -0.00679478, -0.0341808},
956 {0.00423913, 0.0757551, 0.0444685, -0.281896, 0.379112, -0.212121, 0.0594902},
957 {-0.253426, 0.326755, -0.225738, -0.0691578, 0.109788, -0.0812613, -0.0482903},
958 {-0.344417, 0.464815, -0.268857, 0.0137453, 0.146423, -0.0908125, 0.0131566},
959 {0.180595, -0.471044, 0.600536, -0.680218, 0.666876, -0.345305, 0.215572},
960 {-0.202407, 0.290958, -0.0843309, -0.0730275, 0.157111, -0.0836953, 0.0374609},
961 {-0.222727, 0.258844, -0.0676259, -0.123459, 0.18615, -0.0793805, 0.0245286},
962 {-0.320709, 0.351355, -0.22756, 0.019593, 0.0261016, 0.0191958, -0.0449903},
963 {-0.275517, 0.305623, -0.147753, -0.0221624, 0.0685188, -0.00113907, -0.0110694},
964 {-0.237632, 0.222056, -0.0967139, -0.0467247, 0.0196372, 0.013321, -0.0537844},
965 {-0.254455, 0.247849, -0.0727208, -0.0359449, -0.0293587, 0.0650545, -0.0368777},
966 {-0.260665, 0.296247, -0.0244744, 0.0393139, -0.0729511, 0.0271136, 0.00654422},
967 {-0.312108, 0.443109, -0.0993222, 0.0326532, 4.76077e-05, -0.026312, 0.00211329}
968 };
969 while (cost >= cosbin[tmp] && tmp != 28)
970 {
971 tmp++;
972 }
973 if (tmp == 0)
974 tmp += 1;
975 double par_cos[7];
976 for (int j = 0; j < 7; j++)
977 {
978 double cosbin_tmp[3] = {cosbin[tmp - 1], cosbin[tmp], cosbin[tmp + 1]};
979 double par_pip_offset_tmp[3] = {par_pip_offset[tmp - 1][j], par_pip_offset[tmp][j], par_pip_offset[tmp + 1][j]};
980 par_cos[j] = interpolation(cost, cosbin_tmp, par_pip_offset_tmp);
981 }
982 double ptrk_tmp = (ptrk - 0.17) / 0.1;
983 corr_offset = ROOT::Math::ChebyshevN(6, ptrk_tmp, par_cos);
984 return corr_offset;
985 }
986 else if (n == 3)
987 {
988 double par_kp_offset[3][5] = {
989 {0.00, 0.00, 0.00, 0.00, 0.00},
990 {0.90, 0.90, 0.75, 0.15, 0.00},
991 {0.00, 0.00, 0.00, 0.00, 0.00}
992 };
993 double p_bin[6] = {0.175, 0.225, 0.275, 0.325, 0.375, 0.425};
994 int bin_p = (ptrk - 0.175) / 0.05;
995 double int_p1 = (par_kp_offset[0][bin_p] - par_kp_offset[1][bin_p]) * fabs(cost) / 0.1 + par_kp_offset[1][bin_p];
996 double int_p2 = (par_kp_offset[0][bin_p + 1] - par_kp_offset[1][bin_p + 1]) * fabs(cost) / 0.1 + par_kp_offset[1][bin_p + 1];
997 corr_offset = (int_p2 - int_p1) * (ptrk - p_bin[bin_p]) + int_p1;
998 return corr_offset;
999 }
1000 else if (n == 4)
1001 {
1002 int tmp = 0;
1003 double par_p_offset[30][6] = {
1004 {-0.826976, 1.26319, 0.0168621, -0.350471, 0.208162, -0.0422268},
1005 {-0.655279, 1.17488, -0.624155, 0.0140827, 0.16105, -0.119258},
1006 {-0.316389, 0.678636, -0.746556, 0.476035, -0.189224, 0.0435277},
1007 {0.151886, -0.0224519, -0.287477, 0.274142, -0.153749, 0.0734597},
1008 {0.448262, -0.362655, -0.0415588, 0.110434, -0.0833664, 0.0432422},
1009 {0.622144, -0.564482, 0.13891, -0.0124711, -0.00501423, -0.00480475},
1010 {0.70422, -0.680696, 0.15964, -0.0148526, -0.0430492, -0.0110849},
1011 {0.910677, -0.964731, 0.367566, -0.107582, 0.0527411, 0},
1012 {1.0095, -1.10729, 0.394171, -0.119731, 0.0157713, 0},
1013 {1.22613, -1.46063, 0.636896, -0.195013, 0.0624278, 0},
1014 {1.54499, -1.90653, 0.818578, -0.312619, 0.0763852, 0},
1015 {1.98669, -2.52962, 1.10066, -0.335902, 0.0635311, 0},
1016 {2.58488, -3.35617, 1.53409, -0.555987, 0.160924, 0},
1017 {3.2245, -4.24687, 1.85847, -0.643344, 0.1679, 0},
1018 {3.88496, -4.6193, 1.77549, -0.482142, 0.0735315, 0},
1019 {3.72702, -4.64294, 1.89108, -0.510733, 0.0914982, 0},
1020 {2.44759, -3.709, 1.62889, -0.530121, 0.122125, 0},
1021 {1.48439, -2.52461, 1.10631, -0.38799, 0.0801175, 0},
1022 {0.90893, -1.66597, 0.723416, -0.209917, 0.0630468, 0},
1023 {0.560033, -1.07386, 0.398224, -0.103849, 0.00455556, 0},
1024 {0.419805, -0.809094, 0.318519, -0.0969289, 0.0118289, 0},
1025 {0.287319, -0.436034, 0.0837204, 0.0581752, -0.0808111, 0.0359385},
1026 {0.258884, -0.357818, 0.0475784, 0.0384534, -0.0402014, -0.000481745},
1027 {0.21925, -0.246146, 0.00484353, 0.0531687, -0.031646, 0.0224279},
1028 {0.16258, -0.193137, -0.0671331, 0.0707033, -0.0645643, 0.00244481},
1029 {0.145212, -0.095416, -0.152871, 0.155015, -0.0852023, 0.0387718},
1030 {0.00824374, 0.111452, -0.341958, 0.296775, -0.163263, 0.0740253},
1031 {-0.348253, 0.699121, -0.755615, 0.479183, -0.188644, 0.043546},
1032 {-0.662414, 1.16992, -0.664637, 0.0413632, 0.145736, -0.110265},
1033 {-0.852959, 1.28192, -0.0518724, -0.297676, 0.164896, -0.0247413}
1034 };
1035 while (cost >= cosbin[tmp] && tmp != 28)
1036 {
1037 tmp++;
1038 }
1039 if (tmp == 0)
1040 tmp += 1;
1041 double par_cos[6];
1042 for (int j = 0; j < 6; j++)
1043 {
1044 double cosbin_tmp[3] = {cosbin[tmp - 1], cosbin[tmp], cosbin[tmp + 1]};
1045 double par_p_offset_tmp[3] = {par_p_offset[tmp - 1][j], par_p_offset[tmp][j], par_p_offset[tmp + 1][j]};
1046 par_cos[j] = interpolation(cost, cosbin_tmp, par_p_offset_tmp);
1047 }
1048 double ptrk_tmp = (ptrk - 0.33) / 0.1;
1049 corr_offset = ROOT::Math::Chebyshev5(ptrk_tmp, par_cos[0], par_cos[1], par_cos[2], par_cos[3], par_cos[4], par_cos[5]);
1050 if (cost > 0.83 && ptrk < 0.3)
1051 return 2 * corr_offset;
1052 return corr_offset;
1053 }
1054}
1055
1056double DedxPID::sigmaCorr(int n, int charge, double ptrk, double cost)
1057{
1058 double cosbin[30] = {-0.865, -0.7, -0.5, -0.325, -0.225, -0.19, -0.17, -0.15, -0.13, -0.11, -0.09, -0.07, -0.05, -0.03, -0.01, 0.01, 0.03, 0.05, 0.07, 0.09, 0.11, 0.13, 0.15, 0.17, 0.19, 0.225, 0.325, 0.5, 0.7, 0.865};
1059 double corr_sigma = 1;
1060 if (n == 0 || n == 1 || n == 2 && ptrk > 0.20 || n == 3 && (ptrk > 0.43 || ptrk < 0.2 || fabs(cost) > 0.2) || n == 4 && (ptrk > 0.43 || ptrk < 0.2))
1061 return 1;
1062 else if (n == 2)
1063 {
1064 int tmp = 0;
1065 double par_pip_sigma[30][8] = {
1066 {0.980745, -0.0320824, 0.148076, -0.0185231, -0.0287245, -0.0146609, 0.0291458, -0.0213545},
1067 {0.909279, 0.171988, 0.0389074, -0.0877263, 0.0104833, 0.0392911, -0.0306531, 0.00307295},
1068 {0.952309, 0.168354, 0.0100906, -0.0781473, 0.0874947, -0.0399515, -0.0040687, 0.0273449},
1069 {1.01305, 0.110082, 0.0356032, -0.0320006, 0.0398265, 0.00283998, -0.00353269, 0.0323046},
1070 {1.10214, -0.0464266, 0.0844326, -0.0133869, -0.0233443, 0.0470427, -0.0233875, 0.00796884},
1071 {1.13722, -0.0694168, 0.0810207, -0.00461111, -0.0657331, 0.0435981, -0.0305237, 0.000891762},
1072 {1.12494, -0.0402269, 0.039213, 0.0410197, -0.0849904, 0.0588127, -0.0110273, -0.0102289},
1073 {1.19926, -0.124029, 0.0915294, 0.0051942, -0.0546915, 0.0774849, -0.0186289, 0.00478784},
1074 {1.29678, -0.285694, 0.230571, -0.127148, 0.0295359, -0.0319002, 0.0336125, -0.0387693},
1075 {1.34107, -0.299083, 0.21921, -0.115439, 0.0098274, -0.0225279, 0.0177997, -0.0240218},
1076 {1.37443, -0.292003, 0.170567, -0.0571655, -0.0168999, 0.0336943, 0.00351935, -0.0213522},
1077 {1.53528, -0.469584, 0.286359, -0.136798, -0.0158706, -0.00297735, 0.00902674, -0.0365105},
1078 {1.76254, -0.697432, 0.422368, -0.120523, -0.070091, 0.094326, -0.030649, -0.0066399},
1079 {1.94569, -0.911353, 0.454219, 0.0156298, -0.319211, 0.379627, -0.238878, 0.0963599},
1080 {2.43193, -1.66075, 1.11878, -0.488709, 0.0340719, 0.160446, -0.0920272, 0.0514607},
1081 {2.0932, -1.05913, 0.580453, -0.0946158, -0.236143, 0.29226, -0.176364, 0.0320375},
1082 {1.92749, -0.929276, 0.553641, -0.170987, -0.098254, 0.135239, -0.075694, -0.0032753},
1083 {1.875, -0.988737, 0.679303, -0.365717, 0.0721803, 0.0117188, -0.00415287, -0.0189484},
1084 {1.58964, -0.578412, 0.43666, -0.221821, 0.051573, -0.00910371, 0.0148332, -0.036281},
1085 {1.37654, -0.303804, 0.192889, -0.0902226, -0.00981288, -0.0260383, 0.0428074, -0.0312006},
1086 {1.33936, -0.328366, 0.265274, -0.150611, 0.0536288, -0.0520476, 0.0561347, -0.0397721},
1087 {1.21193, -0.144586, 0.0934788, -0.0178466, -0.0531063, 0.0124997, 0.0153325, -0.0307769},
1088 {1.17235, -0.0870064, 0.0758076, 0.00109802, -0.0561468, 0.030234, 0.011407, -0.0287855},
1089 {1.16807, -0.127506, 0.118003, -0.0407013, -0.00810034, -0.00531038, 0.0271216, -0.0245021},
1090 {1.10302, -0.0128453, 0.0366757, 0.0430047, -0.0677478, 0.0658191, -0.0213085, 0.0231517},
1091 {1.09825, -0.0253014, 0.0807453, -0.0315216, -0.00550889, 0.0126133, 0.00452082, 0.00183966},
1092 {1.0199, 0.0937959, 0.0449512, -0.0539498, 0.0544941, -0.0364448, 0.0205853, 0.0103639},
1093 {0.964461, 0.159884, 0.0289211, -0.0837419, 0.0984221, -0.0551503, 0.00340403, 0.0213683},
1094 {0.92457, 0.159631, 0.0423011, -0.07931, 0.0114676, 0.0488719, -0.0348463, 0.00461647},
1095 {1.00597, -0.0291862, 0.131163, -0.03148, -0.024414, -0.039103, 0.0488831, -0.0346602}
1096 };
1097 while (cost >= cosbin[tmp] && tmp != 28)
1098 {
1099 tmp++;
1100 }
1101 if (tmp == 0)
1102 tmp += 1;
1103 double par_cos[8];
1104 for (int j = 0; j < 8; j++)
1105 {
1106 double cosbin_tmp[3] = {cosbin[tmp - 1], cosbin[tmp], cosbin[tmp + 1]};
1107 double par_pip_sigma_tmp[3] = {par_pip_sigma[tmp - 1][j], par_pip_sigma[tmp][j], par_pip_sigma[tmp + 1][j]};
1108 par_cos[j] = interpolation(cost, cosbin_tmp, par_pip_sigma_tmp);
1109 }
1110 double ptrk_tmp = (ptrk - 0.17) / 0.1;
1111 double corr_sigma = ROOT::Math::ChebyshevN(7, ptrk_tmp, par_cos);
1112 if (corr_sigma < 1)
1113 return 1;
1114 else
1115 return corr_sigma;
1116 }
1117 else if (n == 3)
1118 {
1119 double par_kp_sigma[3][5] = {
1120 {1.00, 1.00, 1.00, 1.00, 1.00},
1121 {1.80, 1.80, 1.51, 1.41, 1.00},
1122 {1.00, 1.00, 1.00, 1.00, 1.00}
1123 };
1124 double p_bin[6] = {0.175, 0.225, 0.275, 0.325, 0.375, 0.425};
1125 int bin_p = (ptrk - 0.175) / 0.05;
1126 double int_p1 = (par_kp_sigma[0][bin_p] - par_kp_sigma[1][bin_p]) * fabs(cost) / 0.1 + par_kp_sigma[1][bin_p];
1127 double int_p2 = (par_kp_sigma[0][bin_p + 1] - par_kp_sigma[1][bin_p + 1]) * fabs(cost) / 0.1 + par_kp_sigma[1][bin_p + 1];
1128 corr_sigma = (int_p2 - int_p1) * (ptrk - p_bin[bin_p]) + int_p1;
1129 return corr_sigma;
1130 }
1131 else if (n == 4)
1132 {
1133 int tmp = 0;
1134 double par_p_sigma[30][8] = {
1135 {0.794024, 0.0425693, 0.0236678, -0.0382406, 0.0695961, -0.0580967, 0.035697, -0.0135807},
1136 {0.832773, -0.00113245, -0.031817, 0.0606602, -0.0447306, -0.00903627, 0.025789, -0.0195913},
1137 {0.908858, -0.087108, 0.0549567, 0.00174534, -0.0270899, 0.0429156, -0.0280865, 0.0188789},
1138 {1.04046, -0.246353, 0.133491, -0.049544, 0.0180147, 0, 0, 0},
1139 {1.25697, -0.492783, 0.244496, -0.0930121, 0.0267921, 0, 0, 0},
1140 {1.40495, -0.656157, 0.341844, -0.13557, 0.0444445, 0, 0, 0},
1141 {1.48819, -0.722884, 0.375376, -0.133594, 0.0550627, 0, 0, 0},
1142 {1.73349, -1.02811, 0.545484, -0.22501, 0.0867905, 0, 0, 0},
1143 {1.86727, -1.11375, 0.566508, -0.209777, 0.0683113, 0, 0, 0},
1144 {2.17391, -1.50475, 0.78278, -0.317744, 0.0926452, 0, 0, 0},
1145 {2.4923, -1.78499, 0.944323, -0.412239, 0.144967, 0, 0, 0},
1146 {2.96861, -2.37577, 1.24553, -0.50482, 0.135875, 0, 0, 0},
1147 {3.31789, -2.67592, 1.3589, -0.552132, 0.210676, 0, 0, 0},
1148 {3.7896, -3.26956, 1.68685, -0.702016, 0.155152, 0, 0, 0},
1149 {3.86579, -3.22667, 1.54792, -0.607399, 0.174962, 0, 0, 0},
1150 {3.91034, -3.35332, 1.66152, -0.618642, 0.152329, 0, 0, 0},
1151 {3.33904, -2.66857, 1.27733, -0.460009, 0.0800364, 0, 0, 0},
1152 {2.97639, -2.27119, 1.1101, -0.410801, 0.120688, 0, 0, 0},
1153 {2.55881, -1.83093, 0.905367, -0.339415, 0.107886, 0, 0, 0},
1154 {2.34426, -1.713, 0.866848, -0.370726, 0.087723, 0, 0, 0},
1155 {1.98031, -1.24099, 0.640849, -0.232078, 0.0726222, 0, 0, 0},
1156 {1.74302, -0.984163, 0.490949, -0.172, 0.0443975, 0, 0, 0},
1157 {1.56317, -0.802742, 0.388115, -0.14842, 0.0359668, 0, 0, 0},
1158 {1.44037, -0.668254, 0.352312, -0.120142, 0.0549672, 0, 0, 0},
1159 {1.34493, -0.583195, 0.310501, -0.130395, 0.0447765, 0, 0, 0},
1160 {1.22836, -0.433327, 0.229097, -0.0728195, 0.022962, 0, 0, 0},
1161 {1.05117, -0.246895, 0.142671, -0.0529643, 0.016318, 0, 0, 0},
1162 {0.909469, -0.0691198, 0.0377954, 0.019234, -0.0322931, 0.0460066, -0.0270032, 0.02252},
1163 {0.843402, -0.0106399, -0.0217012, 0.0502854, -0.0341327, -0.0117776, 0.0292822, -0.0190088},
1164 {0.826268, -0.00178627, 0.0679738, -0.065918, 0.0696007, -0.0648257, 0.0328222, -0.00459817}
1165 };
1166 while (cost >= cosbin[tmp] && tmp != 28)
1167 {
1168 tmp++;
1169 }
1170 if (tmp == 0)
1171 tmp += 1;
1172 double par_cos[8];
1173 for (int j = 0; j < 8; j++)
1174 {
1175 double cosbin_tmp[3] = {cosbin[tmp - 1], cosbin[tmp], cosbin[tmp + 1]};
1176 double par_p_sigma_tmp[3] = {par_p_sigma[tmp - 1][j], par_p_sigma[tmp][j], par_p_sigma[tmp + 1][j]};
1177 par_cos[j] = interpolation(cost, cosbin_tmp, par_p_sigma_tmp);
1178 }
1179 double ptrk_tmp = (ptrk - 0.33) / 0.1;
1180 corr_sigma = ROOT::Math::ChebyshevN(7, ptrk_tmp, par_cos);
1181 if (corr_sigma < 1)
1182 return 1;
1183 else
1184 return corr_sigma;
1185 }
1186}
1187
1188double DedxPID::interpolation(double cost, double *costheta, double *par)
1189{
1190 double ux = (costheta[0] - costheta[1]) * (costheta[0] - costheta[2]);
1191 double uy = (costheta[1] - costheta[0]) * (costheta[1] - costheta[2]);
1192 double uz = (costheta[2] - costheta[0]) * (costheta[2] - costheta[1]);
1193 double bx = par[0] / ux;
1194 double by = par[1] / uy;
1195 double bz = par[2] / uz;
1196 double c1 = bx + by + bz;
1197 double c2 = bx * (costheta[1] + costheta[2]) + by * (costheta[0] + costheta[2]) + bz * (costheta[0] + costheta[1]);
1198 double c3 = bx * costheta[1] * costheta[2] + by * costheta[0] * costheta[2] + bz * costheta[0] * costheta[1];
1199 return c1 * cost * cost - c2 * cost + c3;
1200}
1201//CT
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
runNo
int runNo
Definition: DQA_TO_DB.cxx:12
n
const Int_t n
Definition: DataBase/tau_mode.c:65
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
bin
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per bin
Definition: FoamA.h:85
DedxPID::mypol3
double mypol3(double x, double par0, double par1, double par2, double par3)
Definition: DedxPID.cxx:518
DedxPID::mypol5
double mypol5(double x, double par0, double par1, double par2, double par3, double par4, double par5)
Definition: DedxPID.cxx:525
DedxPID::getNhitCutDx
int getNhitCutDx() const
Definition: DedxPID.h:34
DedxPID::offsetCorr
double offsetCorr(int n, int charge, double ptrk, double cost)
Definition: DedxPID.cxx:928
DedxPID::particleIDCalculation
int particleIDCalculation()
Definition: DedxPID.cxx:51
DedxPID::CorrDedx
double CorrDedx(int n, double ptrk, double cost, double chi, int charge)
Definition: DedxPID.cxx:171
DedxPID::sigmaDedx
double sigmaDedx(int n, double ptrk, double cost)
Definition: DedxPID.cxx:467
DedxPID::inputpar
void inputpar()
Definition: DedxPID.cxx:532
DedxPID::iterate
double iterate(double ptrk, double *mean, double *p)
Definition: DedxPID.cxx:906
DedxPID::init
void init()
Definition: DedxPID.cxx:26
DedxPID::offsetDedx
double offsetDedx(int n, double ptrk, double cost)
Definition: DedxPID.cxx:167
DedxPID::interpolation
double interpolation(double cost, double *costheta, double *par)
Definition: DedxPID.cxx:1188
DedxPID::offset
double offset(int n) const
Definition: DedxPID.h:28
DedxPID::cal_par
double cal_par(int index1, double *m_jpsi_pip_ptrk_offset, double ptrk, double begin, double bin)
Definition: DedxPID.cxx:915
DedxPID::calculate
void calculate()
Definition: DedxPID.cxx:42
DedxPID::chi
double chi(int n) const
Definition: DedxPID.h:26
DedxPID::instance
static DedxPID * instance()
Definition: DedxPID.cxx:17
DedxPID::sigmaCorr
double sigmaCorr(int n, int charge, double ptrk, double cost)
Definition: DedxPID.cxx:1056
DstMdcDedx::probPH
double probPH() const
Definition: DstMdcDedx.h:66
DstMdcDedx::numGoodHits
int numGoodHits() const
Definition: DstMdcDedx.h:64
DstMdcDedx::normPH
double normPH() const
Definition: DstMdcDedx.h:67
DstMdcDedx::chi
double chi(int i) const
Definition: DstMdcDedx.h:58
DstMdcKalTrack::theta
const double theta() const
Definition: DstMdcKalTrack.h:70
DstMdcKalTrack::setPidType
static void setPidType(PidType pidType)
Definition: DstMdcKalTrack.h:59
DstMdcKalTrack::p
const double p() const
Definition: DstMdcKalTrack.h:82
DstMdcTrack::theta
const double theta() const
Definition: DstMdcTrack.h:59
DstMdcTrack::charge
const int charge() const
Definition: DstMdcTrack.h:53
DstMdcTrack::p
const double p() const
Definition: DstMdcTrack.h:58
EvtRecTrack::isMdcDedxValid
bool isMdcDedxValid()
Definition: EvtRecTrack.h:45
EvtRecTrack::mdcDedx
RecMdcDedx * mdcDedx()
Definition: EvtRecTrack.h:55
EvtRecTrack::isMdcKalTrackValid
bool isMdcKalTrackValid()
Definition: EvtRecTrack.h:44
EvtRecTrack::isMdcTrackValid
bool isMdcTrackValid()
Definition: EvtRecTrack.h:43
EvtRecTrack::mdcTrack
RecMdcTrack * mdcTrack()
Definition: EvtRecTrack.h:53
EvtRecTrack::mdcKalTrack
RecMdcKalTrack * mdcKalTrack()
Definition: EvtRecTrack.h:54
ParticleIDBase::chiMinCut
double chiMinCut() const
Definition: ParticleIDBase.h:38
ParticleIDBase::PidTrk
EvtRecTrack * PidTrk() const
Definition: ParticleIDBase.h:36
ParticleIDBase::probCalculate
double probCalculate(double chi2, int n)
Definition: ParticleIDBase.cxx:83
ParticleIDBase::getRunNo
double getRunNo() const
Definition: ParticleIDBase.h:44
ParticleIDBase::pdfCalculate
double pdfCalculate(double offset, double sigma)
Definition: ParticleIDBase.cxx:91
ParticleIDBase::pdfMinSigmaCut
double pdfMinSigmaCut() const
Definition: ParticleIDBase.h:42
ParticleIDBase::path
static std::string path
Definition: ParticleIDBase.h:102
y
double y[1000]
Definition: draw_charge_space_coarse_II.cxx:8
costheta
float costheta
Definition: prepare_data.cxx:28
charge
float charge
Definition: prepare_data.cxx:28
ptrk
float ptrk
Definition: prepare_data.cxx:28