CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
TofCPID.cxx
Go to the documentation of this file.
1#include <fstream>
2#include <cmath>
3#include <cstdlib>
4
5#include "ParticleID/TofCPID.h"
6
7#ifndef BEAN
8#include "MdcRecEvent/RecMdcTrack.h"
9#include "TofRecEvent/RecTofTrack.h"
10#include "EvtRecEvent/EvtRecTrack.h"
11#include "DstEvent/TofHitStatus.h"
12#else
13#include "TofHitStatus.h"
14#endif
15
16
17TofCPID * TofCPID::m_pointer = 0;
18
19TofCPID * TofCPID::instance() {
20 if(!m_pointer) m_pointer = new TofCPID();
21 return m_pointer;
22}
23
24TofCPID::TofCPID():ParticleIDBase() {
25 m_pars[0]= -0.208289;
26 m_pars[1]= 1.63092;
27 m_pars[2]= -0.0733139;
28 m_pars[3]= 1.02385;
29 m_pars[4]= 0.00145052;
30 m_pars[5]= -1.72471e-06;
31 m_pars[6]= 5.92726e-10;
32 m_pars[7]= -0.0645428;
33 m_pars[8]= -0.00143637;
34 m_pars[9]= -0.133817;
35 m_pars[10]= 0.0101188;
36 m_pars[11]= -5.07622;
37 m_pars[12]= 5.31472;
38 m_pars[13]= -0.443142;
39 m_pars[14]= -12.1083;
40 m_readstate=0;
41}
42
43void TofCPID::init() {
44 for(int i = 0; i < 5; i++) {
45 m_chi[i] = 99.0;
46 m_prob[i] = -1.0;
47 m_offset[i] = 99.0;
48 m_sigma[i] = 1.0;
49 }
50 m_chimin = 99.;
51 m_pdfmin =99;
52 m_ndof = 0;
53}
54
55void TofCPID::calculate() {
56 int runtof = getRunNo();
57 if(!m_readstate) {
58 std::cout<<"read tofC"<<std::endl;
59 std::string tofdata_mom_file = path + "/share/pidparatof/tofpdata.txt";
60 ifstream inputmomdata(tofdata_mom_file.c_str(),std::ios_base::in);
61 if ( !inputmomdata ) {
62 cout << " can not open: " << tofdata_mom_file << endl;
63 exit(1);
64 }
65
66 std::string tofdata_theta_file = path + "/share/pidparatof/tofthetadata.txt";
67 ifstream inputthetadata(tofdata_theta_file.c_str(),std::ios_base::in);
68 if ( !inputthetadata ) {
69 cout << " can not open: " << tofdata_theta_file << endl;
70 exit(1);
71 }
72
73 std::string tofdata_endcap_file = path + "/share/pidparatof/tofendcapdata.txt";
74 ifstream inputendcapdata(tofdata_endcap_file.c_str(),std::ios_base::in);
75 if ( !inputendcapdata ) {
76 cout << " can not open: " << tofdata_endcap_file << endl;
77 exit(1);
78 }
79
80
81 std::string tofmc_mom_file = path + "/share/pidparatof/tofpmc.txt";
82 ifstream inputmommc(tofmc_mom_file.c_str(),std::ios_base::in);
83 if ( !inputmommc ) {
84 cout << " can not open: " << tofmc_mom_file << endl;
85 exit(1);
86 }
87
88 std::string tofmc_theta_file = path + "/share/pidparatof/tofthetamc.txt";
89 ifstream inputthetamc(tofmc_theta_file.c_str(),std::ios_base::in);
90 if ( !inputthetamc ) {
91 cout << " can not open: " << tofmc_theta_file << endl;
92 exit(1);
93 }
94
95 std::string tofmc_endcap_file = path + "/share/pidparatof/tofendcapmc.txt";
96 ifstream inputendcapmc(tofmc_endcap_file.c_str(),std::ios_base::in);
97 if ( !inputendcapmc ) {
98 cout << " can not open: " << tofmc_endcap_file << endl;
99 exit(1);
100 }
101
102
103 if(runtof>0)
104 {
105 for(int i=0; i<5; i++)
106 {
107 for(int j=0; j<8; j++)
108 {
109 inputthetadata>>m_thetapara[i][j];
110 }
111 }
112
113 for(int i=0; i<5; i++)
114 {
115 for(int j=0; j<12; j++)
116 {
117 inputmomdata>>m_momentpara[i][j];
118 }
119 }
120
121 for(int i=0; i<5; i++)
122 {
123 for(int j=0; j<4; j++)
124 {
125 inputendcapdata>>m_endcappara[i][j];
126 }
127 }
128
129 } else
130 {
131 for(int i=0; i<5; i++)
132 {
133 for(int j=0; j<8; j++)
134 {
135 inputthetamc>>m_thetapara[i][j];
136 }
137 }
138
139 for(int i=0; i<5; i++)
140 {
141 for(int j=0; j<12; j++)
142 {
143 inputmommc>>m_momentpara[i][j];
144 }
145 }
146
147 for(int i=0; i<5; i++)
148 {
149 for(int j=0; j<4; j++)
150 {
151 inputendcapmc>>m_endcappara[i][j];
152 }
153 }
154
155 }
156 m_readstate=1;
157 }
158 if(particleIDCalculation() == 0) m_ndof=1;
159}
160
161int TofCPID::particleIDCalculation() {
162 /*
163 cout<<"m_momentpara[2][2]="<<m_momentpara[2][2]<<endl;
164 cout<<"m_momentpara[2][3]="<<m_momentpara[2][3]<<endl;
165 cout<<"m_momentpara[3][2]="<<m_momentpara[3][2]<<endl;
166 cout<<"m_momentpara[3][3]="<<m_momentpara[3][3]<<endl;
167 cout<<"m_thetapara[2][2]="<<m_thetapara[2][2]<<endl;
168 cout<<"m_thetapara[2][3]="<<m_thetapara[2][3]<<endl;
169 cout<<"m_thetapara[3][2]="<<m_thetapara[3][2]<<endl;
170 cout<<"m_thetapara[3][3]="<<m_thetapara[3][3]<<endl;
171 cout<<"m_endcappara[2][2]="<<m_endcappara[2][2]<<endl;
172 cout<<"m_endcappara[2][3]="<<m_endcappara[2][3]<<endl;
173 cout<<"m_endcappara[3][2]="<<m_endcappara[3][2]<<endl;
174 cout<<"m_endcappara[3][3]="<<m_endcappara[3][3]<<endl;
175 */
176 int irc = -1;
177
178 EvtRecTrack* recTrk = PidTrk();
179 if(!(recTrk->isMdcTrackValid())) return irc;
180 RecMdcTrack* mdcTrk = recTrk->mdcTrack();
181 double ptrk = mdcTrk->p();
182 // double charge = mdcTrk->charge();
183 double cost = cos(mdcTrk->theta());
184 if(!(recTrk->isTofTrackValid())) return irc;
185
186#ifndef BEAN
187 SmartRefVector<RecTofTrack> tofTrk = recTrk->tofTrack();
188 SmartRefVector<RecTofTrack>::iterator it;//=tofTrk.begin();
189#else
190 const std::vector<TTofTrack* >& tofTrk = recTrk->tofTrack();
191 std::vector<TTofTrack* >::const_iterator it;//=tofTrk.begin();
192#endif
193
194 TofHitStatus *hitst = new TofHitStatus;
195 std::vector<int> tofccount;
196 int goodtofctrk=0;
197 for(it = tofTrk.begin(); it!=tofTrk.end(); it++,goodtofctrk++) {
198 unsigned int st = (*it)->status();
199 hitst->setStatus(st);
200 // if( !(hitst->is_barrel()) ) continue;
201 // if( !(hitst->is_counter()) ) continue;
202 // if( hitst->layer()==1 ) tofccount.push_back(goodtofctrk);
203 if(hitst->is_cluster()) tofccount.push_back(goodtofctrk);
204 }
205 delete hitst;
206 if(tofccount.size()!=1) return irc;//not tof2 track or more than 1 tracks
207 it = tofTrk.begin()+tofccount[0];
208 double tof = (*it)->tof();
209 m_tofc = tof;
210 // int qual = (*it)->quality();
211 // int cntr = (*it)->tofID();
212 double path = ((*it)->path())*10.0;//the unit from mm to cm
213 m_pathc = path;
214 m_phc = (*it)->ph(); //no change
215 m_zhitc = ((*it)->zrhit())*10;//the unit from mm to cm
216 double beta2 = path*path/velc()/velc()/tof/tof;
217 m_mass2 = ptrk * ptrk * (1/beta2 -1);
218 if ((m_mass2>20)||(m_mass2<-1)) return irc;
219 if(tof <=0 ) return irc;
220 double chitemp = 99.;
221 double pdftemp = 0;
222 // double sigma_tmp= (*it)->sigma(0);
223 double testchi[5];
224 double testpdf[5];
225 for(int i = 0; i < 5; i++) {
226 /*
227 m_offset[i] = tof - (*it)->texp(i);//- offsetTofC(i, cntr, ptrk, m_zhit1, m_ph1,charge);
228 if(sigma_tmp!=0) m_sigma[i] = 1.1*sigma_tmp/1000.;
229 else
230 m_sigma[i]=sigmaTofC(i, cntr,ptrk,m_zhitc, m_phc,charge);
231 m_chi[i] = m_offset[i]/m_sigma[i];
232 */
233 double sep = tof - (*it)->texp(i)-(*it)->toffset(i);
234 m_chi[i] = (sep - offsetTofC(i, ptrk, cost))/sigmaTofC(i, ptrk, cost);
235 m_offset[i] = offsetTofC(i, ptrk, cost);
236 m_sigma[i] = sigmaTofC(i, ptrk, cost);
237 testchi[i]=sep;
238 if(fabs(m_chi[i]) < chitemp) chitemp = fabs(m_chi[i]);
239 double ppp = pdfCalculate(m_chi[i],1);
240 testpdf[i]=ppp;
241 if(fabs(ppp) > pdftemp) pdftemp = fabs(ppp);
242 }
243 m_chimin = chitemp;
244 m_pdfmin = pdftemp;
245 if(pdftemp < pdfCalculate(pdfMinSigmaCut(),1)) return irc;
246 if(fabs(m_chimin) > chiMinCut()) return irc;
247 for(int i = 0; i < 5; i++) {
248 m_prob[i] = probCalculate(m_chi[i]*m_chi[i], 1);
249 }
250
251 irc = 0;
252 return irc;
253
254
255
256}
257
258//
259// dE/dx Correction routines
260//
261
262double TofCPID::offsetTofC(int n, double ptrk, double cost) {
263 int rundedx2 = getRunNo();
264 double offset = 0.0;
265 double offsetp = 0.0;
266 double offsetc = 0.0;
267 double sigcos = 0.0;
268 double sigp = 0.0;
269 // double gb = ptrk/xmass(n);
270
271 switch(n) {
272 case 0: { // Electron
273 double ptemp = ptrk;
274 double costm = cost;
275
276 if(rundedx2>0)
277 { if(ptrk < 0.3) ptemp = 0.3;
278 if(ptrk > 1.3) ptemp = 1.3;
279 }
280 else
281 { if(ptrk < 0.3) ptemp = 0.3;
282 if(ptrk > 1.3) ptemp = 1.3;
283 }
284
285 double plog = log(ptemp);
286 double costcos = cos(costm);
287 offsetp= mypol5(plog,m_momentpara[0][0],m_momentpara[0][1],m_momentpara[0][2],m_momentpara[0][3],m_momentpara[0][4],m_momentpara[0][5]);
288 sigp=mypol5(plog,m_momentpara[0][6],m_momentpara[0][7],m_momentpara[0][8],m_momentpara[0][9],m_momentpara[0][10],m_momentpara[0][11]);
289
290 if(costm<-0.83) {
291 offsetc=m_endcappara[0][0];
292 sigcos=m_endcappara[0][2];
293 }
294 if(costm>0.83) {
295 offsetc=m_endcappara[0][1];
296 sigcos=m_endcappara[0][3];
297 }
298 if(fabs(costm)<=0.83)
299 {
300 offsetc=mypol3(costcos,m_thetapara[0][0],m_thetapara[0][1],m_thetapara[0][2],m_thetapara[0][3]);
301 sigcos=mypol3(costcos,m_thetapara[0][4],m_thetapara[0][5],m_thetapara[0][6],m_thetapara[0][7]);
302 }
303
304
305 offset=offsetc+sigcos*offsetp;
306 //offset=offsetc;
307 offset=offsetp+sigp*offsetc;
308 break;
309 }
310
311 case 1: {// Muon
312 double ptemp = ptrk;
313 double costm = cost;
314 if(rundedx2>0)
315 { if(ptrk < 0.3) ptemp = 0.3;
316 if(ptrk > 1.3) ptemp = 1.3;
317 }
318 else
319 { if(ptrk < 0.3) ptemp = 0.3;
320 if(ptrk > 1.3) ptemp = 1.3;
321 }
322
323 double plog = log(ptemp);
324 double costcos = cos(costm);
325 offsetp= mypol5(plog,m_momentpara[1][0],m_momentpara[1][1],m_momentpara[1][2],m_momentpara[1][3],m_momentpara[1][4],m_momentpara[1][5]);
326 sigp=mypol5(plog,m_momentpara[1][6],m_momentpara[1][7],m_momentpara[1][8],m_momentpara[1][9],m_momentpara[1][10],m_momentpara[1][11]);
327
328 if(costm<-0.83) {
329 offsetc=m_endcappara[1][0];
330 sigcos=m_endcappara[1][2];
331 }
332 if(costm>0.83) {
333 offsetc=m_endcappara[1][1];
334 sigcos=m_endcappara[1][3];
335 }
336 if(fabs(costm)<=0.83)
337 {
338 offsetc=mypol3(costcos,m_thetapara[1][0],m_thetapara[1][1],m_thetapara[1][2],m_thetapara[1][3]);
339 sigcos=mypol3(costcos,m_thetapara[1][4],m_thetapara[1][5],m_thetapara[1][6],m_thetapara[1][7]);
340 }
341
342
343 offset=offsetc+sigcos*offsetp;
344 //offset=offsetc;
345 offset=offsetp+sigp*offsetc;
346 break;
347 }
348
349 case 2: {// Pion
350 double ptemp = ptrk;
351 double costm = cost;
352 if(rundedx2>0)
353 { if(ptrk < 0.3) ptemp = 0.3;
354 if(ptrk > 1.6) ptemp = 1.6;
355 }
356 else
357 { if(ptrk < 0.3) ptemp = 0.3;
358 if(ptrk > 1.6) ptemp = 1.6;
359 }
360
361 double plog = log(ptemp);
362 double costcos = cos(costm);
363 offsetp= mypol5(plog,m_momentpara[2][0],m_momentpara[2][1],m_momentpara[2][2],m_momentpara[2][3],m_momentpara[2][4],m_momentpara[2][5]);
364 sigp=mypol5(plog,m_momentpara[2][6],m_momentpara[2][7],m_momentpara[2][8],m_momentpara[2][9],m_momentpara[2][10],m_momentpara[2][11]);
365
366 if(costm<-0.83) {
367 offsetc=m_endcappara[2][0];
368 sigcos=m_endcappara[2][2];
369 }
370 if(costm>0.83) {
371 offsetc=m_endcappara[2][1];
372 sigcos=m_endcappara[2][3];
373 }
374 if(fabs(costm)<=0.83)
375 {
376 offsetc=mypol3(costcos,m_thetapara[2][0],m_thetapara[2][1],m_thetapara[2][2],m_thetapara[2][3]);
377 sigcos=mypol3(costcos,m_thetapara[2][4],m_thetapara[2][5],m_thetapara[2][6],m_thetapara[2][7]);
378 }
379
380
381 offset=offsetc+sigcos*offsetp;
382 //offset=offsetc;
383 offset=offsetp+sigp*offsetc;
384 break;
385 }
386
387 case 3: {// Kaon
388 double ptemp = ptrk;
389 double costm = cost;
390 if(rundedx2>0)
391 { if(ptrk < 0.4) ptemp = 0.4;
392 if(ptrk > 1.3) ptemp = 1.3;
393 }
394 else
395 { if(ptrk < 0.4) ptemp = 0.4;
396 if(ptrk > 1.3) ptemp = 1.3;
397 }
398 double plog = log(ptemp);
399 double costcos = cos(costm);
400 offsetp= mypol5(plog,m_momentpara[3][0],m_momentpara[3][1],m_momentpara[3][2],m_momentpara[3][3],m_momentpara[3][4],m_momentpara[3][5]);
401 sigp=mypol5(plog,m_momentpara[3][6],m_momentpara[3][7],m_momentpara[3][8],m_momentpara[3][9],m_momentpara[3][10],m_momentpara[3][11]);
402
403 if(costm<-0.83) {
404 offsetc=m_endcappara[3][0];
405 sigcos=m_endcappara[3][2];
406 }
407 if(costm>0.83) {
408 offsetc=m_endcappara[3][1];
409 sigcos=m_endcappara[3][3];
410 }
411 if(fabs(costm)<=0.83)
412 {
413 offsetc=mypol3(costcos,m_thetapara[3][0],m_thetapara[3][1],m_thetapara[3][2],m_thetapara[3][3]);
414 sigcos=mypol3(costcos,m_thetapara[3][4],m_thetapara[3][5],m_thetapara[3][6],m_thetapara[3][7]);
415 }
416
417
418 offset=offsetc+sigcos*offsetp;
419 //offset=offsetc;
420 offset=offsetp+sigp*offsetc;
421 break;
422 }
423
424 case 4 : { // Proton
425 double ptemp = ptrk;
426 double costm = cost;
427 if(rundedx2>0)
428 { if(ptrk < 0.5) ptemp = 0.5;
429 if(ptrk > 1.3) ptemp = 1.3;
430 }
431 else
432 { if(ptrk < 0.5) ptemp = 0.5;
433 if(ptrk > 1.3) ptemp = 1.3;
434 }
435 double plog = log(ptemp);
436 double costcos = cos(costm);
437 offsetp= mypol5(plog,m_momentpara[4][0],m_momentpara[4][1],m_momentpara[4][2],m_momentpara[4][3],m_momentpara[4][4],m_momentpara[4][5]);
438 sigp=mypol5(plog,m_momentpara[4][6],m_momentpara[4][7],m_momentpara[4][8],m_momentpara[4][9],m_momentpara[4][10],m_momentpara[4][11]);
439
440 if(costm<-0.83) {
441 offsetc=m_endcappara[4][0];
442 sigcos=m_endcappara[4][2];
443 }
444 if(costm>0.83) {
445 offsetc=m_endcappara[4][1];
446 sigcos=m_endcappara[4][3];
447 }
448 if(fabs(costm)<=0.83)
449 {
450 offsetc=mypol3(costcos,m_thetapara[4][0],m_thetapara[4][1],m_thetapara[4][2],m_thetapara[4][3]);
451 sigcos=mypol3(costcos,m_thetapara[4][4],m_thetapara[4][5],m_thetapara[4][6],m_thetapara[4][7]);
452 }
453
454 offset=offsetc+sigcos*offsetp;
455 //offset=offsetc;
456 offset=offsetp+sigp*offsetc;
457 break;
458 }
459
460 default:
461 offset = 0.0;
462 break;
463 }
464 // offset = 0.0;
465 return offset;
466}
467
468
469
470double TofCPID::sigmaTofC(int n, double ptrk, double cost) {
471 int rundedx3 = getRunNo();
472 double sigma = 1.0;
473 double sigmap = 1.0;
474 double sigmac = 1.0;
475 // double gb = ptrk/xmass(n);
476 switch(n) {
477
478 case 0: {// Electron
479 double ptemp = ptrk;
480 double costm = cost;
481 if(rundedx3>0)
482 { if(ptrk < 0.3) ptemp = 0.3;
483 if(ptrk > 1.3) ptemp = 1.3;
484 }
485 else
486 { if(ptrk < 0.3) ptemp = 0.3;
487 if(ptrk > 1.3) ptemp = 1.3;
488 }
489
490 double plog = log(ptemp);
491 double costcos = cos(costm);
492
493 sigmap=mypol5(plog,m_momentpara[0][6],m_momentpara[0][7],m_momentpara[0][8],m_momentpara[0][9],m_momentpara[0][10],m_momentpara[0][11]);
494
495 if(costm<-0.83) {
496 sigmac=m_endcappara[0][2];
497 }
498 if(costm>0.83) {
499 sigmac=m_endcappara[0][3];
500 }
501 if(fabs(costm)<0.83)
502 {
503 sigmac=mypol3(costcos,m_thetapara[0][4],m_thetapara[0][5],m_thetapara[0][6],m_thetapara[0][7]);
504 }
505
506 sigma=sigmap*sigmac;
507 //sigma=sigmac;
508 break;
509 }
510
511 case 1: {// Muon
512 double ptemp = ptrk;
513 double costm = cost;
514 if(rundedx3>0)
515 { if(ptrk < 0.3) ptemp = 0.3;
516 if(ptrk > 1.3) ptemp = 1.3;
517 }
518 else
519 { if(ptrk < 0.3) ptemp = 0.3;
520 if(ptrk > 1.3) ptemp = 1.3;
521 }
522
523 double plog = log(ptemp);
524 double costcos = cos(costm);
525
526 sigmap=mypol5(plog,m_momentpara[1][6],m_momentpara[1][7],m_momentpara[1][8],m_momentpara[1][9],m_momentpara[1][10],m_momentpara[1][11]);
527
528 if(costm<-0.83) {
529 sigmac=m_endcappara[1][2];
530 }
531 if(costm>0.83) {
532 sigmac=m_endcappara[1][3];
533 }
534 if(fabs(costm)<0.83)
535 {
536 sigmac=mypol3(costcos,m_thetapara[1][4],m_thetapara[1][5],m_thetapara[1][6],m_thetapara[1][7]);
537 }
538
539
540 sigma=sigmap*sigmac;
541 //sigma=sigmac;
542 break;
543 }
544
545 case 2: {// Pion
546 double ptemp = ptrk;
547 double costm = cost;
548 if(rundedx3>0)
549 { if(ptrk < 0.3) ptemp = 0.3;
550 if(ptrk > 1.6) ptemp = 1.6;
551 }
552 else
553 { if(ptrk < 0.3) ptemp = 0.3;
554 if(ptrk > 1.6) ptemp = 1.6;
555 }
556
557 double plog = log(ptemp);
558 double costcos = cos(costm);
559 sigmap=mypol5(plog,m_momentpara[2][6],m_momentpara[2][7],m_momentpara[2][8],m_momentpara[2][9],m_momentpara[2][10],m_momentpara[2][11]);
560
561 if(costm<-0.83) {
562 sigmac=m_endcappara[2][2];
563 }
564 if(costm>0.83) {
565 sigmac=m_endcappara[2][3];
566 }
567 if(fabs(costm)<0.83)
568 {
569 sigmac=mypol3(costcos,m_thetapara[2][4],m_thetapara[2][5],m_thetapara[2][6],m_thetapara[2][7]);
570 }
571
572 sigma=sigmap*sigmac;
573 //sigma=sigmac;
574
575 break;
576 }
577
578 case 3: { // Kaon
579 double ptemp = ptrk;
580 double costm = cost;
581
582 if(rundedx3>0)
583 { if(ptrk < 0.4) ptemp = 0.4;
584 if(ptrk > 1.3) ptemp = 1.3;
585 }
586 else
587 { if(ptrk < 0.4) ptemp = 0.4;
588 if(ptrk > 1.3) ptemp = 1.3;
589 }
590 double plog = log(ptemp);
591 double costcos = cos(costm);
592 sigmap=mypol5(plog,m_momentpara[3][6],m_momentpara[3][7],m_momentpara[3][8],m_momentpara[3][9],m_momentpara[3][10],m_momentpara[3][11]);
593
594 if(costm<-0.83) {
595 sigmac=m_endcappara[3][2];
596 }
597 if(costm>0.83) {
598 sigmac=m_endcappara[3][3];
599 }
600 if(fabs(costm)<0.83)
601 {
602 sigmac=mypol3(costcos,m_thetapara[3][4],m_thetapara[3][5],m_thetapara[3][6],m_thetapara[3][7]);
603 }
604
605 sigma=sigmap*sigmac;
606 //sigma=sigmac;
607 break;
608 }
609
610
611 case 4: {// Proton
612 double ptemp = ptrk;
613 double costm = cost;
614 if(rundedx3>0)
615 { if(ptrk < 0.5) ptemp = 0.5;
616 if(ptrk > 1.3) ptemp = 1.3;
617 }
618 else
619 { if(ptrk < 0.5) ptemp = 0.5;
620 if(ptrk > 1.3) ptemp = 1.3;
621 }
622 double plog = log(ptemp);
623 double costcos = cos(costm);
624 sigmap=mypol5(plog,m_momentpara[4][6],m_momentpara[4][7],m_momentpara[4][8],m_momentpara[4][9],m_momentpara[4][10],m_momentpara[4][11]);
625
626 if(costm<-0.83) {
627 sigmac=m_endcappara[4][2];
628 }
629 if(costm>0.83) {
630 sigmac=m_endcappara[4][3];
631 }
632 if(fabs(costm)<0.83)
633 {
634 sigmac=mypol3(costcos,m_thetapara[4][4],m_thetapara[4][5],m_thetapara[4][6],m_thetapara[4][7]);
635 }
636
637 sigma=sigmap*sigmac;
638 //sigma=sigmac;
639 break;
640 }
641
642 default:
643 sigma = 1.0;
644 break;
645 }
646
647 // sigma =1.0;
648 return sigma;
649}
650
651double TofCPID::mypol3(double x, double par0, double par1, double par2, double par3)
652{
653 double y = x;
654 return par0 + (par1 * y) +(par2 * y * y) + (par3 * y * y * y);
655
656}
657
658double TofCPID::mypol5(double x, double par0, double par1, double par2, double par3, double par4, double par5)
659{
660 double y = x;
661 return par0 + (par1 * y) +(par2 * y * y) + (par3 * y * y * y) + (par4 * y * y * y *y)+ (par5 * y * y * y * y * y);
662
663}
664
665
666
n
const Int_t n
Definition: DataBase/tau_mode.c:65
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
ParticleIDBase::probCalculate
double probCalculate(double chi2, int n)
Definition: ParticleIDBase.cxx:84
ParticleIDBase::pdfCalculate
double pdfCalculate(double offset, double sigma)
Definition: ParticleIDBase.cxx:92
TofCPID::calculate
void calculate()
Definition: TofCPID.cxx:55
TofCPID::particleIDCalculation
int particleIDCalculation()
Definition: TofCPID.cxx:161
TofCPID::init
void init()
Definition: TofCPID.cxx:43
TofCPID::mypol3
double mypol3(double x, double par0, double par1, double par2, double par3)
Definition: TofCPID.cxx:651
TofCPID::mypol5
double mypol5(double x, double par0, double par1, double par2, double par3, double par4, double par5)
Definition: TofCPID.cxx:658
TofCPID::instance
static TofCPID * instance()
Definition: TofCPID.cxx:19
TofCPID::offsetTofC
double offsetTofC(int n, double ptrk, double cost)
Definition: TofCPID.cxx:262
TofCPID::sigmaTofC
double sigmaTofC(int n, double ptrk, double cost)
Definition: TofCPID.cxx:470
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition: TofHitStatus.cxx:126