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