CGEM BOSS 6.6.5.h
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
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 nhitcutdedx=getNhitCutDx();
54 int irc = -1;
55 EvtRecTrack* recTrk = PidTrk();
56 if(!(recTrk->isMdcTrackValid())) return irc;
57 RecMdcTrack* mdcTrk = recTrk->mdcTrack();
58
59 double ptrk = mdcTrk->p();
60 int charge = mdcTrk->charge();
61 if(ptrk>5) return irc;
62 double cost = cos(mdcTrk->theta());
63 // double sig_the= sin(mdcTrk->theta());
64
65 if(!(recTrk->isMdcDedxValid())) return irc;
66 RecMdcDedx* dedxTrk = recTrk->mdcDedx();
67
68 if((dedxTrk->normPH()>30)||(dedxTrk->normPH()<0)) return irc;
69 m_goodHits = dedxTrk->numGoodHits();
70 if(dedxTrk->numGoodHits()<nhitcutdedx) return irc;
71 m_normPH = dedxTrk->normPH();
72 m_probPH = dedxTrk->probPH();
73 // calculate chi and min chi
74 double chitemp = 99.;
75 double pdftemp = 0;
76 // double testchi[5];
77 // double testptrk[5];
78 // double testcost[5];
79 for(int i = 0; i < 5; i++) {
80 double sep = dedxTrk->chi(i);
81
82#ifndef BEAN
83 string sftver = getenv("BES_RELEASE");
84 string sft;
85 sft.assign(sftver,0,5);
86 if(sft=="6.6.0"||sft=="6.5.5") {
87 m_chi[i] = CorrDedx(i,ptrk,cost,sep,charge);
88 }
89 else
90 m_chi[i]=sep;
91#else
92// This is BEAN part:
93#if (ROOTEVENTDATA_VERSION_NUMERIC == ROOTEVENTDATA_VER(6,5,5) ||\
94 ROOTEVENTDATA_VERSION_NUMERIC == ROOTEVENTDATA_VER(6,6,0) )
95 m_chi[i] = CorrDedx(i,ptrk,cost,sep,charge);
96#else
97 m_chi[i]=sep;
98#endif
99#endif
100
101 m_offset[i] = offsetDedx(i, ptrk, cost);
102 m_sigma[i] = sigmaDedx(i, ptrk, cost);
103 if(fabs(m_chi[i]) < chitemp) chitemp = fabs(m_chi[i]);
104 double ppp = pdfCalculate(m_chi[i],1);
105 if(fabs(ppp) > pdftemp) pdftemp = fabs(ppp);
106
107 }
108 m_chimin = chitemp;
109 m_pdfmin = pdftemp;
110 if(m_chimin > chiMinCut()) return irc;
111 if(pdftemp < pdfCalculate(pdfMinSigmaCut(),1)) return irc;
112
113
114 // calculate prob
115
116 for(int i = 0; i < 5; i++)
117 m_prob[i] = probCalculate(m_chi[i]*m_chi[i], 1);
118
119 m_ndof = 1;
120 irc = 0;
121 return irc;
122}
123
124//
125// dE/dx Correction routines
126//
127
128
129
130double DedxPID::offsetDedx(int n, double ptrk, double cost) {
131 return 0;
132}
133
134double DedxPID::CorrDedx(int n, double ptrk, double cost,double chi,int charge) {
135 int rundedx2 = getRunNo();
136 double offset = 0.0;
137 double offsetp = 0.0;
138 double offsetc = 0.0;
139 double sigcos = 1;
140 double sigp = 1;
141 double chicor=chi;
142 // double gb = ptrk/xmass(n);
143
144 switch(n) {
145 case 0: { // Electron
146 break;
147 }
148
149 case 1: {// Muon
150 break;
151 }
152
153 case 2: {// Pion
154 // double ptemp = ptrk;
155 double costm = cost;
156 if(ptrk<0.1||ptrk>1) break;
157 int index = int((ptrk-0.1)/0.05);
158 if(index<=0) index=1;
159 if(index>=17) index=16;
160
161 if(fabs(costm)>=0.8) break;
162 int index1 = int((costm+0.8)/0.1);
163 if(index1<=0) index1=1;
164 if(index1>=15) index1=14;
165
166 //psipp data
167 if(rundedx2>=11414&&rundedx2<=14604) {
168 offsetp = cal_par(index,m_psipp_pi_ptrk_offset,ptrk,0.125,0.05);
169 sigp = cal_par(index,m_psipp_pi_ptrk_sigma,ptrk,0.125,0.05);
170 offsetc = cal_par(index1,m_psipp_pi_theta_offset,costm,-0.75,0.1);
171 sigcos = cal_par(index1,m_psipp_pi_theta_sigma,costm,-0.75,0.1);
172 }
173 //psipp mc
174 if(rundedx2<=-11414&&rundedx2>=-14604) {
175 offsetp = cal_par(index,m_psipp_mc_pi_ptrk_offset,ptrk,0.125,0.05);
176 sigp = cal_par(index,m_psipp_mc_pi_ptrk_sigma,ptrk,0.125,0.05);
177 offsetc = cal_par(index1,m_psipp_mc_pi_theta_offset,costm,-0.75,0.1);
178 sigcos = cal_par(index1,m_psipp_mc_pi_theta_sigma,costm,-0.75,0.1);
179 }
180
181 offset=offsetp+sigp*offsetc;
182 chicor=(chicor-offset)/(sigcos*sigp);
183 break;
184 }
185
186 case 3: {// Kaon
187 // double ptemp = ptrk;
188 double costm = cost;
189 if(ptrk<0.3||ptrk>0.8) break;
190 offset=0;
191 int index = int((ptrk-0.3)/0.1);
192 if(index<=0) index=1;
193 if(index>=4) index=3;
194
195 int index1 = int((costm+0.9)/0.1);
196 if(index1<=0) index1=1;
197 if(index1>=17) index1=16;
198 //data Jpsi
199 if(rundedx2>=9947&&rundedx2<=10878) {
200 if(charge>0) {
201 offsetp = cal_par(index,m_jpsi_kap_ptrk_offset,ptrk,0.35,0.1);
202 sigp = cal_par(index,m_jpsi_kap_ptrk_sigma,ptrk,0.35,0.1);
203 if(fabs(costm)<=0.83) {
204 offsetc = cal_par(index1,m_jpsi_kap_theta_offset,costm,-0.85,0.1);
205 sigcos = cal_par(index1,m_jpsi_kap_theta_sigma,costm,-0.85,0.1);
206 }
207 }
208 if(charge<0) {
209 offsetp = cal_par(index,m_jpsi_kam_ptrk_offset,ptrk,0.35,0.1);
210 sigp = cal_par(index,m_jpsi_kam_ptrk_sigma,ptrk,0.35,0.1);
211 if(fabs(costm)<=0.83) {
212 offsetc = cal_par(index1,m_jpsi_kam_theta_offset,costm,-0.85,0.1);
213 sigcos = cal_par(index1,m_jpsi_kam_theta_sigma,costm,-0.85,0.1);
214 }
215 }
216 }
217
218 //mc Jpsi
219 if(rundedx2<=-9947&&rundedx2>=-10878) {
220 if(charge>0) {
221 offsetp = cal_par(index,m_jpsi_mc_kap_ptrk_offset,ptrk,0.35,0.1);
222 sigp = cal_par(index,m_jpsi_mc_kap_ptrk_sigma,ptrk,0.35,0.1);
223 if(fabs(costm)<=0.83) {
224 offsetc = cal_par(index1,m_jpsi_mc_kap_theta_offset,costm,-0.85,0.1);
225 sigcos = cal_par(index1,m_jpsi_mc_kap_theta_sigma,costm,-0.85,0.1);
226 }
227 }
228 if(charge<0) {
229 offsetp = cal_par(index,m_jpsi_mc_kam_ptrk_offset,ptrk,0.35,0.1);
230 sigp = cal_par(index,m_jpsi_mc_kam_ptrk_sigma,ptrk,0.35,0.1);
231 if(fabs(costm)<=0.83) {
232 offsetc = cal_par(index1,m_jpsi_mc_kam_theta_offset,costm,-0.85,0.1);
233 sigcos = cal_par(index1,m_jpsi_mc_kam_theta_sigma,costm,-0.85,0.1);
234 }
235 }
236 }
237
238 //data Psip
239 if(rundedx2>=8093&&rundedx2<=9025) {
240 if(ptrk<0.3||ptrk>1.2) break;
241 index = int((ptrk-0.3)/0.1);
242 if(index<=0) index=1;
243 if(index>=8) index=7;
244 if(charge>0) {
245 offsetp = cal_par(index,m_psip_kap_ptrk_offset,ptrk,0.35,0.1);
246 sigp = cal_par(index,m_psip_kap_ptrk_sigma,ptrk,0.35,0.1);
247 }
248 if(charge<0) {
249 offsetp = cal_par(index,m_psip_kam_ptrk_offset,ptrk,0.35,0.1);
250 sigp = cal_par(index,m_psip_kam_ptrk_sigma,ptrk,0.35,0.1);
251 }
252 }
253
254 //mc Psip
255 if(rundedx2<=-8093&&rundedx2>=-9025) {
256 // if(ptrk < 0.4) ptrk = 0.4;
257 if(ptrk<0.3||ptrk>1.2) break;
258 index = int((ptrk-0.3)/0.1);
259 if(index<=0) index=1;
260 if(index>=8) index=7;
261 if(charge>0) {
262 offsetp = cal_par(index,m_psip_mc_kap_ptrk_offset,ptrk,0.35,0.1);
263 sigp = cal_par(index,m_psip_mc_kap_ptrk_sigma,ptrk,0.35,0.1);
264 }
265 if(charge<0) {
266 offsetp = cal_par(index,m_psip_mc_kam_ptrk_offset,ptrk,0.35,0.1);
267 sigp = cal_par(index,m_psip_mc_kam_ptrk_sigma,ptrk,0.35,0.1);
268 }
269 }
270
271
272 //psipp kaon data
273 if(rundedx2>=11414&&rundedx2<=14604) {
274 if(ptrk<0.15||ptrk>1) break;
275 index = int((ptrk-0.15)/0.05);
276 if(index<=0) index=1;
277 if(index>=16) index=15;
278 if(fabs(costm)>=0.8) break;
279 index1 = int((costm+0.8)/0.1);
280 if(index1<=0) index1=1;
281 if(index1>=15) index1=14;
282
283 offsetp = cal_par(index,m_psipp_ka_ptrk_offset,ptrk,0.175,0.05);
284 sigp = cal_par(index,m_psipp_ka_ptrk_sigma,ptrk,0.175,0.05);
285 offsetc = cal_par(index1,m_psipp_ka_theta_offset,costm,-0.75,0.1);
286 sigcos = cal_par(index1,m_psipp_ka_theta_sigma,costm,-0.75,0.1);
287 }
288 //psipp kaon mc
289 if(rundedx2<=-11414&&rundedx2>=-14604) {
290 if(ptrk<0.15||ptrk>1) break;
291 index = int((ptrk-0.15)/0.05);
292 if(index<=0) index=1;
293 if(index>=16) index=15;
294 if(fabs(costm)>=0.8) break;
295 index1 = int((costm+0.8)/0.1);
296 if(index1<=0) index1=1;
297 if(index1>=15) index1=14;
298 offsetp = cal_par(index,m_psipp_mc_ka_ptrk_offset,ptrk,0.175,0.05);
299 sigp = cal_par(index,m_psipp_mc_ka_ptrk_sigma,ptrk,0.175,0.05);
300 offsetc = cal_par(index1,m_psipp_mc_ka_theta_offset,costm,-0.75,0.1);
301 sigcos = cal_par(index1,m_psipp_mc_ka_theta_sigma,costm,-0.75,0.1);
302 }
303
304 offset=offsetp+sigp*offsetc;
305 chicor=(chicor-offset)/(sigcos*sigp);
306 break;
307 }
308
309 case 4 : { // Proton
310 // double ptemp = ptrk;
311 double costm = cost;
312 if(ptrk<0.3||ptrk>1.1) break;
313 int index = int((ptrk-0.3)/0.1);
314 if(index<=0) index=1;
315 if(index>=7) index=6;
316
317 int index1 = int((costm+0.9)/0.1);
318 if(index1<=0) index1=1;
319 if(index1>=17) index1=16;
320
321 // double plog = log(ptemp);
322 offset=0;
323 if(rundedx2>=9947&&rundedx2<=10878) {
324 if(charge>0) {
325 offsetp = cal_par(index,m_jpsi_protonp_ptrk_offset,ptrk,0.35,0.1);
326 sigp = cal_par(index,m_jpsi_protonp_ptrk_sigma,ptrk,0.35,0.1);
327 if(fabs(costm)<=0.83) {
328 offsetc = cal_par(index1,m_jpsi_protonp_theta_offset,costm,-0.85,0.1);
329 sigcos = cal_par(index1,m_jpsi_protonp_theta_sigma,costm,-0.85,0.1);
330 }
331 }
332 if(charge<0) {
333 offsetp = cal_par(index,m_jpsi_protonm_ptrk_offset,ptrk,0.35,0.1);
334 sigp = cal_par(index,m_jpsi_protonm_ptrk_sigma,ptrk,0.35,0.1);
335 if(fabs(costm)<=0.83) {
336 offsetc = cal_par(index1,m_jpsi_protonm_theta_offset,costm,-0.85,0.1);
337 sigcos = cal_par(index1,m_jpsi_protonm_theta_sigma,costm,-0.85,0.1);
338 }
339 }
340 }
341
342 //mc JPsi
343 if(rundedx2<=-9947&&rundedx2>=-10878) {
344 if(charge>0) {
345 offsetp = cal_par(index,m_jpsi_mc_protonp_ptrk_offset,ptrk,0.35,0.1);
346 sigp = cal_par(index,m_jpsi_mc_protonp_ptrk_sigma,ptrk,0.35,0.1);
347 if(fabs(costm)<=0.83) {
348 offsetc = cal_par(index1,m_jpsi_mc_protonp_theta_offset,costm,-0.85,0.1);
349 sigcos = cal_par(index1,m_jpsi_mc_protonp_theta_sigma,costm,-0.85,0.1);
350 }
351 }
352 if(charge<0) {
353 offsetp = cal_par(index,m_jpsi_mc_protonm_ptrk_offset,ptrk,0.35,0.1);
354 sigp = cal_par(index,m_jpsi_mc_protonm_ptrk_sigma,ptrk,0.35,0.1);
355 if(fabs(costm)<=0.83) {
356 offsetc = cal_par(index1,m_jpsi_mc_protonm_theta_offset,costm,-0.85,0.1);
357 sigcos = cal_par(index1,m_jpsi_mc_protonm_theta_sigma,costm,-0.85,0.1);
358 }
359 }
360 }
361
362 //data Psip
363 if(rundedx2>=8093&&rundedx2<=9025) {
364 if(charge>0) {
365 offsetp = cal_par(index,m_psip_protonp_ptrk_offset,ptrk,0.35,0.1);
366 sigp = cal_par(index,m_psip_protonp_ptrk_sigma,ptrk,0.35,0.1);
367 }
368 if(charge<0) {
369 offsetp = cal_par(index,m_psip_protonm_ptrk_offset,ptrk,0.35,0.1);
370 sigp = cal_par(index,m_psip_protonm_ptrk_sigma,ptrk,0.35,0.1);
371 }
372 }
373
374 //mc Psip
375 if(rundedx2<=-8093&&rundedx2>=-9025) {
376 if(charge>0) {
377 offsetp = cal_par(index,m_psip_mc_protonp_ptrk_offset,ptrk,0.35,0.1);
378 sigp = cal_par(index,m_psip_mc_protonp_ptrk_sigma,ptrk,0.35,0.1);
379 }
380 if(charge<0) {
381 offsetp = cal_par(index,m_psip_mc_protonm_ptrk_offset,ptrk,0.35,0.1);
382 sigp = cal_par(index,m_psip_mc_protonm_ptrk_sigma,ptrk,0.35,0.1);
383 }
384 }
385
386 //psipp proton data
387 if(rundedx2>=11414&&rundedx2<=14604) {
388 if(ptrk<0.2||ptrk>1.1) break;
389 index = int((ptrk-0.2)/0.05);
390 if(index<=0) index=1;
391 if(index>=17) index=16;
392 if(fabs(costm)>=0.83) break;
393 index1 = int((costm+0.9)/0.1);
394 if(index1<=0) index1=1;
395 if(index1>=17) index1=16;
396
397 offsetp = cal_par(index,m_psipp_proton_ptrk_offset,ptrk,0.225,0.05);
398 sigp = cal_par(index,m_psipp_proton_ptrk_sigma,ptrk,0.225,0.05);
399 offsetc = cal_par(index1,m_psipp_proton_theta_offset,costm,-0.85,0.1);
400 sigcos = cal_par(index1,m_psipp_proton_theta_sigma,costm,-0.85,0.1);
401 }
402 //psipp proton mc
403 if(rundedx2<=-11414&&rundedx2>=-14604) {
404 if(ptrk<0.2||ptrk>1.1) break;
405 index = int((ptrk-0.2)/0.1);
406 if(index<=0) index=1;
407 if(index>=8) index=7;
408 if(fabs(costm)>=0.83) break;
409 index1 = int((costm+0.9)/0.1);
410 if(index1<=0) index1=1;
411 if(index1>=17) index1=16;
412 offsetp = cal_par(index,m_psipp_mc_proton_ptrk_offset,ptrk,0.25,0.1);
413 sigp = cal_par(index,m_psipp_mc_proton_ptrk_sigma,ptrk,0.25,0.1);
414 offsetc = cal_par(index1,m_psipp_mc_proton_theta_offset,costm,-0.85,0.1);
415 sigcos = cal_par(index1,m_psipp_mc_proton_theta_sigma,costm,-0.85,0.1);
416 }
417 offset=offsetp+sigp*offsetc;
418 chicor=(chicor-offset)/(sigcos*sigp);
419 break;
420 }
421
422 default:
423 offset = 0.0;
424 break;
425 }
426 // offset = 0.0;
427 return chicor;
428}
429
430double DedxPID::sigmaDedx(int n, double ptrk, double cost) {
431
432 /* int rundedx3 = getRunNo();
433 double sigma = 1.0;
434 double sigp = 1.0;
435 double sigmac = 1.0;
436 double gb = ptrk/xmass(n);
437 switch(n) {
438
439 case 0: {// Electron
440 double ptemp = ptrk;
441 double costm = cost;
442 break;
443 }
444
445 case 1: {// Muon
446 double ptemp = ptrk;
447 double costm = cost;
448 break;
449 }
450
451 case 2: {// Pion
452 double ptemp = ptrk;
453 double costm = cost;
454 break;
455 }
456
457 case 3: { // Kaon
458 double ptemp = ptrk;
459 double costm = cost;
460 break;
461 }
462
463
464 case 4: {// Proton
465 double ptemp = ptrk;
466 double costm = cost;
467 break;
468 }
469
470 default:
471 sigma = 1.0;
472 break;
473 }
474 */
475 // sigma = 1.2;
476 // sigma =1.0;
477 return 1;
478 // return sigma;
479}
480
481double DedxPID::mypol3(double x, double par0, double par1, double par2, double par3)
482{
483 double y = x;
484 return par0 + (par1 * y) +(par2 * y * y) + (par3 * y * y * y);
485
486}
487
488double DedxPID::mypol5(double x, double par0, double par1, double par2, double par3, double par4, double par5)
489{
490 double y = x;
491 return par0 + (par1 * y) +(par2 * y * y) + (par3 * y * y * y) + (par4 * y * y * y *y)+ (par5 * y * y * y * y * y);
492
493}
494
495void DedxPID::inputpar() {
496
497 //Jpsi ka+ momentum correction
498 std::string jpsi_kap_mom = path + "/share/JPsi/kaon/dedx_kap.txt";
499 std::string jpsi_kap_mom_mc = path + "/share/JPsi/kaon/dedx_kap_mc.txt";
500 ifstream inputmomdata6(jpsi_kap_mom.c_str(),std::ios_base::in);
501 if ( !inputmomdata6 ) {
502 cout << " can not open: " << jpsi_kap_mom << endl;
503 exit(1);
504 }
505 ifstream inputmomdata6mc(jpsi_kap_mom_mc.c_str(),std::ios_base::in);
506 if ( !inputmomdata6mc ) {
507 cout << " can not open: " << jpsi_kap_mom_mc << endl;
508 exit(1);
509 }
510 for(int i=0; i<12; i++) {
511 inputmomdata6>>m_jpsi_kap_ptrk_offset[i];
512 inputmomdata6>>m_jpsi_kap_ptrk_sigma[i];
513 inputmomdata6mc>>m_jpsi_mc_kap_ptrk_offset[i];
514 inputmomdata6mc>>m_jpsi_mc_kap_ptrk_sigma[i];
515 }
516
517 //Jpsi ka- momentum correction
518 std::string jpsi_kam_mom = path + "/share/JPsi/kaon/dedx_kam.txt";
519 std::string jpsi_kam_mom_mc = path + "/share/JPsi/kaon/dedx_kam_mc.txt";
520 ifstream inputmomdata7(jpsi_kam_mom.c_str(),std::ios_base::in);
521 if ( !inputmomdata7 ) {
522 cout << " can not open: " << jpsi_kam_mom << endl;
523 exit(1);
524 }
525 ifstream inputmomdata7mc(jpsi_kam_mom_mc.c_str(),std::ios_base::in);
526 if ( !inputmomdata7mc ) {
527 cout << " can not open: " << jpsi_kam_mom_mc << endl;
528 exit(1);
529 }
530 for(int i=0; i<12; i++) {
531 inputmomdata7>>m_jpsi_kam_ptrk_offset[i];
532 inputmomdata7>>m_jpsi_kam_ptrk_sigma[i];
533 inputmomdata7mc>>m_jpsi_mc_kam_ptrk_offset[i];
534 inputmomdata7mc>>m_jpsi_mc_kam_ptrk_sigma[i];
535
536 }
537
538
539 //Jpsi ka+ theta correction
540 std::string jpsi_kap_the = path + "/share/JPsi/kaon/dedx_kap_theta.txt";
541 std::string jpsi_kap_the_mc = path + "/share/JPsi/kaon/dedx_kap_theta_mc.txt";
542 ifstream inputmomdata8(jpsi_kap_the.c_str(),std::ios_base::in);
543 if ( !inputmomdata8 ) {
544 cout << " can not open: " << jpsi_kap_the << endl;
545 exit(1);
546 }
547 ifstream inputmomdata8mc(jpsi_kap_the_mc.c_str(),std::ios_base::in);
548 if ( !inputmomdata8mc ) {
549 cout << " can not open: " << jpsi_kap_the_mc << endl;
550 exit(1);
551 }
552 for(int i=0; i<18; i++) {
553 inputmomdata8>>m_jpsi_kap_theta_offset[i];
554 inputmomdata8>>m_jpsi_kap_theta_sigma[i];
555 inputmomdata8mc>>m_jpsi_mc_kap_theta_offset[i];
556 inputmomdata8mc>>m_jpsi_mc_kap_theta_sigma[i];
557 }
558
559 //Jpsi ka- theta correction
560 std::string jpsi_kam_the = path + "/share/JPsi/kaon/dedx_kam_theta.txt";
561 std::string jpsi_kam_the_mc = path + "/share/JPsi/kaon/dedx_kam_theta_mc.txt";
562 ifstream inputmomdata9(jpsi_kam_the.c_str(),std::ios_base::in);
563 if ( !inputmomdata9 ) {
564 cout << " can not open: " << jpsi_kam_the << endl;
565 exit(1);
566 }
567 ifstream inputmomdata9mc(jpsi_kam_the_mc.c_str(),std::ios_base::in);
568 if ( !inputmomdata9mc ) {
569 cout << " can not open: " << jpsi_kam_the_mc << endl;
570 exit(1);
571 }
572 for(int i=0; i<18; i++) {
573 inputmomdata9>>m_jpsi_kam_theta_offset[i];
574 inputmomdata9>>m_jpsi_kam_theta_sigma[i];
575 inputmomdata9mc>>m_jpsi_mc_kam_theta_offset[i];
576 inputmomdata9mc>>m_jpsi_mc_kam_theta_sigma[i];
577 }
578
579 //Jpsi proton+ momentum correction
580 std::string jpsi_protonp_mom = path + "/share/JPsi/proton/dedx_protonp.txt";
581 std::string jpsi_protonp_mom_mc = path + "/share/JPsi/proton/dedx_protonp_mc.txt";
582 ifstream inputmomdata12(jpsi_protonp_mom.c_str(),std::ios_base::in);
583 if ( !inputmomdata12 ) {
584 cout << " can not open: " << jpsi_protonp_mom << endl;
585 exit(1);
586 }
587 ifstream inputmomdata12mc(jpsi_protonp_mom_mc.c_str(),std::ios_base::in);
588 if ( !inputmomdata12mc ) {
589 cout << " can not open: " << jpsi_protonp_mom_mc << endl;
590 exit(1);
591 }
592 for(int i=0; i<8; i++) {
593 inputmomdata12>>m_jpsi_protonp_ptrk_offset[i];
594 inputmomdata12>>m_jpsi_protonp_ptrk_sigma[i];
595 inputmomdata12mc>>m_jpsi_mc_protonp_ptrk_offset[i];
596 inputmomdata12mc>>m_jpsi_mc_protonp_ptrk_sigma[i];
597 }
598
599 //Jpsi proton- momentum correction
600 std::string jpsi_protonm_mom = path + "/share/JPsi/proton/dedx_protonm.txt";
601 std::string jpsi_protonm_mom_mc = path + "/share/JPsi/proton/dedx_protonm_mc.txt";
602 ifstream inputmomdata13(jpsi_protonm_mom.c_str(),std::ios_base::in);
603 if ( !inputmomdata13 ) {
604 cout << " can not open: " << jpsi_protonm_mom << endl;
605 exit(1);
606 }
607 ifstream inputmomdata13mc(jpsi_protonm_mom_mc.c_str(),std::ios_base::in);
608 if ( !inputmomdata13mc ) {
609 cout << " can not open: " << jpsi_protonm_mom_mc << endl;
610 exit(1);
611 }
612 for(int i=0; i<8; i++) {
613 inputmomdata13>>m_jpsi_protonm_ptrk_offset[i];
614 inputmomdata13>>m_jpsi_protonm_ptrk_sigma[i];
615 inputmomdata13mc>>m_jpsi_mc_protonm_ptrk_offset[i];
616 inputmomdata13mc>>m_jpsi_mc_protonm_ptrk_sigma[i];
617 }
618
619 //Jpsi proton+ theta correction
620 std::string jpsi_protonp_the = path + "/share/JPsi/proton/dedx_protonp_theta.txt";
621 std::string jpsi_protonp_the_mc = path + "/share/JPsi/proton/dedx_protonp_theta_mc.txt";
622
623 ifstream inputmomdata14(jpsi_protonp_the.c_str(),std::ios_base::in);
624 if ( !inputmomdata14 ) {
625 cout << " can not open: " << jpsi_protonp_the << endl;
626 exit(1);
627 }
628 ifstream inputmomdata14mc(jpsi_protonp_the_mc.c_str(),std::ios_base::in);
629 if ( !inputmomdata14mc ) {
630 cout << " can not open: " << jpsi_protonp_the_mc << endl;
631 exit(1);
632 }
633 for(int i=0; i<18; i++) {
634 inputmomdata14>>m_jpsi_protonp_theta_offset[i];
635 inputmomdata14>>m_jpsi_protonp_theta_sigma[i];
636 inputmomdata14mc>>m_jpsi_mc_protonp_theta_offset[i];
637 inputmomdata14mc>>m_jpsi_mc_protonp_theta_sigma[i];
638 }
639
640 //Jpsi proton- theta correction
641 std::string jpsi_protonm_the = path + "/share/JPsi/proton/dedx_protonm_theta.txt";
642 std::string jpsi_protonm_the_mc = path + "/share/JPsi/proton/dedx_protonm_theta_mc.txt";
643 ifstream inputmomdata15(jpsi_protonm_the.c_str(),std::ios_base::in);
644 if ( !inputmomdata15 ) {
645 cout << " can not open: " << jpsi_protonm_the << endl;
646 exit(1);
647 }
648 ifstream inputmomdata15mc(jpsi_protonm_the_mc.c_str(),std::ios_base::in);
649 if ( !inputmomdata15mc ) {
650 cout << " can not open: " << jpsi_protonm_the_mc << endl;
651 exit(1);
652 }
653 for(int i=0; i<18; i++) {
654 inputmomdata15>>m_jpsi_protonm_theta_offset[i];
655 inputmomdata15>>m_jpsi_protonm_theta_sigma[i];
656 inputmomdata15mc>>m_jpsi_mc_protonm_theta_offset[i];
657 inputmomdata15mc>>m_jpsi_mc_protonm_theta_sigma[i];
658 }
659
660
661
662
663 // Psip ka+ momentum correction
664 std::string psip_kap_mom = path + "/share/Psip/kaon/dedx_kap.txt";
665 std::string psip_kap_mom_mc = path + "/share/Psip/kaon/dedx_kap_mc.txt";
666 ifstream inputmomdata24(psip_kap_mom.c_str(),std::ios_base::in);
667 if ( !inputmomdata24 ) {
668 cout << " can not open: " << psip_kap_mom << endl;
669 exit(1);
670 }
671 ifstream inputmomdata24mc(psip_kap_mom_mc.c_str(),std::ios_base::in);
672 if ( !inputmomdata24mc ) {
673 cout << " can not open: " << psip_kap_mom_mc << endl;
674 exit(1);
675 }
676 for(int i=0; i<9; i++) {
677 inputmomdata24>>m_psip_kap_ptrk_offset[i];
678 inputmomdata24>>m_psip_kap_ptrk_sigma[i];
679 inputmomdata24mc>>m_psip_mc_kap_ptrk_offset[i];
680 inputmomdata24mc>>m_psip_mc_kap_ptrk_sigma[i];
681 }
682
683 //Psip ka- momentum correction
684 std::string psip_kam_mom = path + "/share/Psip/kaon/dedx_kam.txt";
685 std::string psip_kam_mom_mc = path + "/share/Psip/kaon/dedx_kam_mc.txt";
686 ifstream inputmomdata25(psip_kam_mom.c_str(),std::ios_base::in);
687 if ( !inputmomdata25 ) {
688 cout << " can not open: " << psip_kam_mom << endl;
689 exit(1);
690 }
691 ifstream inputmomdata25mc(psip_kam_mom_mc.c_str(),std::ios_base::in);
692 if ( !inputmomdata25mc ) {
693 cout << " can not open: " << psip_kam_mom_mc << endl;
694 exit(1);
695 }
696 for(int i=0; i<9; i++) {
697 inputmomdata25>>m_psip_kam_ptrk_offset[i];
698 inputmomdata25>>m_psip_kam_ptrk_sigma[i];
699 inputmomdata25mc>>m_psip_mc_kam_ptrk_offset[i];
700 inputmomdata25mc>>m_psip_mc_kam_ptrk_sigma[i];
701 }
702
703
704 // Psip proton+ momentum correction
705 std::string psip_protonp_mom = path + "/share/Psip/proton/dedx_protonp.txt";
706 std::string psip_protonp_mom_mc = path + "/share/Psip/proton/dedx_protonp_mc.txt";
707 ifstream inputmomdata26(psip_protonp_mom.c_str(),std::ios_base::in);
708 if ( !inputmomdata26 ) {
709 cout << " can not open: " << psip_protonp_mom << endl;
710 exit(1);
711 }
712 ifstream inputmomdata26mc(psip_protonp_mom_mc.c_str(),std::ios_base::in);
713 if ( !inputmomdata26mc ) {
714 cout << " can not open: " << psip_protonp_mom_mc << endl;
715 exit(1);
716 }
717 for(int i=0; i<9; i++) {
718 inputmomdata26>>m_psip_protonp_ptrk_offset[i];
719 inputmomdata26>>m_psip_protonp_ptrk_sigma[i];
720 inputmomdata26mc>>m_psip_mc_protonp_ptrk_offset[i];
721 inputmomdata26mc>>m_psip_mc_protonp_ptrk_sigma[i];
722 }
723
724 //Psip proton- momentum correction
725 std::string psip_protonm_mom = path + "/share/Psip/proton/dedx_protonm.txt";
726 std::string psip_protonm_mom_mc = path + "/share/Psip/proton/dedx_protonm_mc.txt";
727 ifstream inputmomdata27(psip_protonm_mom.c_str(),std::ios_base::in);
728 if ( !inputmomdata27 ) {
729 cout << " can not open: " << psip_protonm_mom << endl;
730 exit(1);
731 }
732 ifstream inputmomdata27mc(psip_protonm_mom_mc.c_str(),std::ios_base::in);
733 if ( !inputmomdata27mc ) {
734 cout << " can not open: " << psip_protonm_mom_mc << endl;
735 exit(1);
736 }
737 for(int i=0; i<9; i++) {
738 inputmomdata27>>m_psip_protonm_ptrk_offset[i];
739 inputmomdata27>>m_psip_protonm_ptrk_sigma[i];
740 inputmomdata27mc>>m_psip_mc_protonm_ptrk_offset[i];
741 inputmomdata27mc>>m_psip_mc_protonm_ptrk_sigma[i];
742 }
743
744 //Psipp pi momentum correction
745 std::string psipp_pi_mom = path + "/share/Psipp/pion/dedx_pi.txt";
746 std::string psipp_pi_mom_mc = path + "/share/Psipp/pion/dedx_pi_mc.txt";
747 ifstream inputmomdata28(psipp_pi_mom.c_str(),std::ios_base::in);
748 if ( !inputmomdata28 ) {
749 cout << " can not open: " << psipp_pi_mom << endl;
750 exit(1);
751 }
752 ifstream inputmomdata28mc(psipp_pi_mom_mc.c_str(),std::ios_base::in);
753 if ( !inputmomdata28mc ) {
754 cout << " can not open: " << psipp_pi_mom_mc << endl;
755 exit(1);
756 }
757 for(int i=0; i<18; i++) {
758 inputmomdata28>>m_psipp_pi_ptrk_offset[i];
759 inputmomdata28>>m_psipp_pi_ptrk_sigma[i];
760 inputmomdata28mc>>m_psipp_mc_pi_ptrk_offset[i];
761 inputmomdata28mc>>m_psipp_mc_pi_ptrk_sigma[i];
762 }
763
764 //Psipp pi theta correction
765 std::string psipp_pi_the = path + "/share/Psipp/pion/dedx_pi_theta.txt";
766 std::string psipp_pi_the_mc = path + "/share/Psipp/pion/dedx_pi_theta_mc.txt";
767 ifstream inputmomdata29(psipp_pi_the.c_str(),std::ios_base::in);
768 if ( !inputmomdata29 ) {
769 cout << " can not open: " << psipp_pi_the << endl;
770 exit(1);
771 }
772 ifstream inputmomdata29mc(psipp_pi_the_mc.c_str(),std::ios_base::in);
773 if ( !inputmomdata29mc ) {
774 cout << " can not open: " << psipp_pi_the_mc << endl;
775 exit(1);
776 }
777 for(int i=0; i<16; i++) {
778 inputmomdata29>>m_psipp_pi_theta_offset[i];
779 inputmomdata29>>m_psipp_pi_theta_sigma[i];
780 inputmomdata29mc>>m_psipp_mc_pi_theta_offset[i];
781 inputmomdata29mc>>m_psipp_mc_pi_theta_sigma[i];
782 }
783
784 //Psipp ka momentum correction
785 std::string psipp_ka_mom = path + "/share/Psipp/kaon/dedx_ka.txt";
786 std::string psipp_ka_mom_mc = path + "/share/Psipp/kaon/dedx_ka_mc.txt";
787 ifstream inputmomdata30(psipp_ka_mom.c_str(),std::ios_base::in);
788 if ( !inputmomdata30 ) {
789 cout << " can not open: " << psipp_ka_mom << endl;
790 exit(1);
791 }
792 ifstream inputmomdata30mc(psipp_ka_mom_mc.c_str(),std::ios_base::in);
793 if ( !inputmomdata30mc ) {
794 cout << " can not open: " << psipp_ka_mom_mc << endl;
795 exit(1);
796 }
797 for(int i=0; i<17; i++) {
798 inputmomdata30>>m_psipp_ka_ptrk_offset[i];
799 inputmomdata30>>m_psipp_ka_ptrk_sigma[i];
800 inputmomdata30mc>>m_psipp_mc_ka_ptrk_offset[i];
801 inputmomdata30mc>>m_psipp_mc_ka_ptrk_sigma[i];
802 }
803
804 //Psipp ka theta correction
805 std::string psipp_ka_the = path + "/share/Psipp/kaon/dedx_ka_theta.txt";
806 std::string psipp_ka_the_mc = path + "/share/Psipp/kaon/dedx_ka_theta_mc.txt";
807 ifstream inputmomdata31(psipp_ka_the.c_str(),std::ios_base::in);
808 if ( !inputmomdata31 ) {
809 cout << " can not open: " << psipp_ka_the << endl;
810 exit(1);
811 }
812 ifstream inputmomdata31mc(psipp_ka_the_mc.c_str(),std::ios_base::in);
813 if ( !inputmomdata31mc ) {
814 cout << " can not open: " << psipp_ka_the_mc << endl;
815 exit(1);
816 }
817 for(int i=0; i<16; i++) {
818 inputmomdata31>>m_psipp_ka_theta_offset[i];
819 inputmomdata31>>m_psipp_ka_theta_sigma[i];
820 inputmomdata31mc>>m_psipp_mc_ka_theta_offset[i];
821 inputmomdata31mc>>m_psipp_mc_ka_theta_sigma[i];
822 }
823
824
825 //Psipp proton momentum correction
826 std::string psipp_proton_mom = path + "/share/Psipp/proton/dedx_proton.txt";
827 std::string psipp_proton_mom_mc = path + "/share/Psipp/proton/dedx_proton_mc.txt";
828 ifstream inputmomdata32(psipp_proton_mom.c_str(),std::ios_base::in);
829 if ( !inputmomdata32 ) {
830 cout << " can not open: " << psipp_proton_mom << endl;
831 exit(1);
832 }
833 ifstream inputmomdata32mc(psipp_proton_mom_mc.c_str(),std::ios_base::in);
834 if ( !inputmomdata32mc ) {
835 cout << " can not open: " << psipp_proton_mom_mc << endl;
836 exit(1);
837 }
838 for(int i=0; i<18; i++) {
839 inputmomdata32>>m_psipp_proton_ptrk_offset[i];
840 inputmomdata32>>m_psipp_proton_ptrk_sigma[i];
841 }
842 for(int i=0; i<9; i++) {
843 inputmomdata32mc>>m_psipp_mc_proton_ptrk_offset[i];
844 inputmomdata32mc>>m_psipp_mc_proton_ptrk_sigma[i];
845 }
846
847 //Psipp proton theta correction
848 std::string psipp_proton_the = path + "/share/Psipp/proton/dedx_proton_theta.txt";
849 std::string psipp_proton_the_mc = path + "/share/Psipp/proton/dedx_proton_theta_mc.txt";
850 ifstream inputmomdata33(psipp_proton_the.c_str(),std::ios_base::in);
851 if ( !inputmomdata33 ) {
852 cout << " can not open: " << psipp_proton_the << endl;
853 exit(1);
854 }
855 ifstream inputmomdata33mc(psipp_proton_the_mc.c_str(),std::ios_base::in);
856 if ( !inputmomdata33mc ) {
857 cout << " can not open: " << psipp_proton_the_mc << endl;
858 exit(1);
859 }
860 for(int i=0; i<18; i++) {
861 inputmomdata33>>m_psipp_proton_theta_offset[i];
862 inputmomdata33>>m_psipp_proton_theta_sigma[i];
863 inputmomdata33mc>>m_psipp_mc_proton_theta_offset[i];
864 inputmomdata33mc>>m_psipp_mc_proton_theta_sigma[i];
865 }
866
867}
868
869double DedxPID::iterate(double ptrk,double *mean,double *p) {
870 double p1,p2,p3;
871 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]));
872 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]));
873 p1=mean[0]-p2*p[0]-p3*p[0]*p[0];
874 double mean1 = p1+p2*ptrk+p3*ptrk*ptrk;
875 return mean1;
876}
877
878double DedxPID::cal_par(int index1,double *m_jpsi_pip_ptrk_offset,double ptrk,double begin,double bin) {
879 double mean1[3],p[3];
880 p[0]=begin+(index1-1)*bin;
881 p[1]=begin+index1*bin;
882 p[2]=begin+(index1+1)*bin;
883 mean1[0]=m_jpsi_pip_ptrk_offset[index1-1];
884 mean1[1]=m_jpsi_pip_ptrk_offset[index1];
885 mean1[2]=m_jpsi_pip_ptrk_offset[index1+1];
886 double res=iterate(ptrk,mean1,p);
887 return res;
888}
889
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
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:481
DedxPID::mypol5
double mypol5(double x, double par0, double par1, double par2, double par3, double par4, double par5)
Definition DedxPID.cxx:488
DedxPID::getNhitCutDx
int getNhitCutDx() const
Definition DedxPID.h:34
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:134
DedxPID::sigmaDedx
double sigmaDedx(int n, double ptrk, double cost)
Definition DedxPID.cxx:430
DedxPID::inputpar
void inputpar()
Definition DedxPID.cxx:495
DedxPID::iterate
double iterate(double ptrk, double *mean, double *p)
Definition DedxPID.cxx:869
DedxPID::init
void init()
Definition DedxPID.cxx:26
DedxPID::offsetDedx
double offsetDedx(int n, double ptrk, double cost)
Definition DedxPID.cxx:130
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:878
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
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
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:49
EvtRecTrack::mdcDedx
RecMdcDedx * mdcDedx()
Definition EvtRecTrack.h:63
EvtRecTrack::isMdcTrackValid
bool isMdcTrackValid()
Definition EvtRecTrack.h:47
EvtRecTrack::mdcTrack
RecMdcTrack * mdcTrack()
Definition EvtRecTrack.h:61
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:84
ParticleIDBase::getRunNo
double getRunNo() const
Definition ParticleIDBase.h:42
ParticleIDBase::pdfCalculate
double pdfCalculate(double offset, double sigma)
Definition ParticleIDBase.cxx:92
ParticleIDBase::pdfMinSigmaCut
double pdfMinSigmaCut() const
Definition ParticleIDBase.h:40
ParticleIDBase::path
static std::string path
Definition ParticleIDBase.h:102