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