BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
SimplePIDSvc.cxx
Go to the documentation of this file.
1
2#include "GaudiKernel/MsgStream.h"
3#include "GaudiKernel/SmartDataPtr.h"
4#include "TFile.h"
5#include "TMath.h"
6#include <cmath>
7#include <cstdlib>
8#include <fstream>
9#include <iostream>
10
11#include "DstEvent/TofHitStatus.h"
12#include "EventModel/EventHeader.h"
13
14#include "SimplePIDSvc.h"
15
16using namespace std;
17DECLARE_COMPONENT( SimplePIDSvc )
18
19SimplePIDSvc::SimplePIDSvc( const std::string& name, ISvcLocator* svcLoc )
20 : base_class( name, svcLoc ) {
21 declareProperty( "DedxChiCut", m_dedx_chi_cut = 4 );
22 declareProperty( "TofChiCut", m_tof_chi_cut = 4 );
23 declareProperty( "IsTofCorr", m_tof_corr = true );
24 declareProperty( "IsDedxCorr", m_dedx_corr = true );
25 declareProperty( "EidRatio", m_eid_ratio = 0.80 );
26}
27
28SimplePIDSvc::~SimplePIDSvc() { ; }
29
30StatusCode SimplePIDSvc::initialize() {
31 MsgStream log( msgSvc(), name() );
32 log << MSG::INFO << "in SimplePIDSvc initialize()" << endmsg;
33
34 StatusCode sc = Service::initialize();
35
36 sc = serviceLocator()->service( "EventDataSvc", eventSvc_, true );
37
38 loadHistogram();
39 loadSecondPar(); // the second correct factor
40
41 return sc;
42}
43
44StatusCode SimplePIDSvc::finalize() {
45 MsgStream log( msgSvc(), name() );
46 log << MSG::INFO << "in SimplePIDSvc finalize()" << endmsg;
47
48 StatusCode sc = Service::finalize();
49
50 for ( unsigned int i = 0; i < 2; i++ )
51 {
52 for ( unsigned int j = 0; j < 4; j++ )
53 {
54 f_dedx[i][j]->Close();
55 f_tof_q[i][j]->Close();
56 f_tof_bgcost[i][j]->Close();
57 f_tof_wgt[i][j]->Close();
58 f_tof_final[i][j]->Close();
59 f_tofec_q[i][j]->Close();
60 f_tofec_bg[i][j]->Close();
61 f_tofec_cost[i][j]->Close();
62 }
63 }
64 for ( unsigned int i = 0; i < 3; i++ )
65 {
66 for ( unsigned int j = 0; j < 4; j++ ) { f_emc[i][j]->Close(); }
67 }
68
69 return sc;
70}
71
72/*StatusCode SimplePIDSvc::queryInterface(const InterfaceID& riid, void** ppvIF)
73 {
74 if (ISimplePIDSvc::interfaceID().versionMatch(riid))
75 {
76 *ppvIF = dynamic_cast<ISimplePIDSvc*>(this);
77 }
78 else
79 {
80 return Service::queryInterface(riid, ppvIF);
81 }
82 addRef();
83 return StatusCode::SUCCESS;
84 }
85 */
86
87void SimplePIDSvc::preparePID( EvtRecTrack* track ) {
88
89 SmartDataPtr<Event::EventHeader> eventHeaderpid( eventSvc_, "/Event/EventHeader" );
90 m_run = eventHeaderpid->runNumber();
91
92 if ( track->isMdcKalTrackValid() )
93 {
94 RecMdcKalTrack* mdckalTrk = track->mdcKalTrack();
95 RecMdcKalTrack::PidType trk_type[5] = {
96 RecMdcKalTrack::electron, RecMdcKalTrack::muon, RecMdcKalTrack::pion,
97 RecMdcKalTrack::kaon, RecMdcKalTrack::proton,
98 };
99 double mass[5] = {
100 0.000511, 0.105658, 0.13957, 0.493677, 0.938272,
101 };
102 for ( unsigned int pid = 0; pid < 5; pid++ )
103 {
104 mdckalTrk->setPidType( trk_type[pid] );
105 m_p[pid] = mdckalTrk->p();
106 m_betagamma[pid] = m_p[pid] / mass[pid];
107 m_charge[pid] = mdckalTrk->charge();
108 m_cost[pid] = cos( mdckalTrk->theta() );
109 }
110 }
111 else
112 {
113 for ( unsigned int i = 0; i < 5; i++ )
114 {
115 m_p[i] = -99;
116 m_betagamma[i] = -99;
117 m_cost[i] = -99;
118 m_charge[i] = 0;
119 }
120 }
121
122 // dE/dx PID
123 loadDedxInfo( track );
124 if ( m_dedx_corr )
125 {
126 dedxCorrection();
127 dedxSecondCorrection();
128 }
129 // TOF PID
130 loadTOFInfo( track );
131 if ( m_tof_corr )
132 {
133 if ( m_tof_barrel == 1 )
134 {
135 tofBarrelCorrection();
136 tofBarrelSecondCorrection();
137 }
138 else if ( m_tof_barrel == 0 )
139 {
140 tofEndcapCorrection();
141 tofEndcapSecondCorrection();
142 }
143 }
144 // EMC
145 loadEMCInfo( track );
146
147 calprob();
148}
149
150void SimplePIDSvc::calprob() {
151 bool usededx = false;
152 bool usetof = false;
153
154 for ( unsigned int i = 0; i < 5; i++ )
155 {
156 if ( !usededx && fabs( m_dedx_chi[i] ) < m_dedx_chi_cut ) usededx = true;
157 if ( !usetof && fabs( m_tof_chi[i] ) < m_tof_chi_cut ) usetof = true;
158
159 m_dedx_only[i] = false;
160 }
161 if ( !usededx )
162 {
163 for ( unsigned int i = 0; i < 5; i++ ) m_dedx_chi[i] = -99;
164 }
165 if ( !usetof )
166 {
167 for ( unsigned int i = 0; i < 5; i++ ) m_tof_chi[i] = -99;
168 }
169
170 for ( unsigned int i = 0; i < 5; i++ )
171 {
172 m_prob[i] = -99;
173 double chi2 = 0;
174 int ndf = 0;
175
176 if ( usededx && usetof )
177 {
178 chi2 = pow( m_dedx_chi[i], 2 ) + pow( m_tof_chi[i], 2 );
179 ndf = 2;
180 }
181 else if ( usededx && !usetof )
182 {
183 chi2 = pow( m_dedx_chi[i], 2 );
184 ndf = 1;
185 m_dedx_only[i] = true;
186 }
187 else if ( !usededx && usetof )
188 {
189 chi2 = pow( m_tof_chi[i], 2 );
190 ndf = 1;
191 }
192 if ( ndf > 0 && chi2 > 0 ) m_prob[i] = TMath::Prob( chi2, ndf );
193 }
194}
195
196int SimplePIDSvc::getRunIdx( int run_no ) {
197 // -1: beyond correction region
198 // 0: 2010 psi(3770) data
199 // 1: 2011 psi(3770) data
200 // 2: 2010 psi(3770) mc
201 // 3: 2011 psi(3770) mc
202 const int RUN_BEGIN_DATA_10 = 11414;
203 const int RUN_END_DATA_10 = 14604;
204 const int RUN_BEGIN_MC_10 = -14604;
205 const int RUN_END_MC_10 = -11414;
206 const int RUN_BEGIN_DATA_11 = 20448;
207 const int RUN_END_DATA_11 = 23454;
208 const int RUN_BEGIN_MC_11 = -23454;
209 const int RUN_END_MC_11 = -20448;
210
211 const int RUN_BEGIN_DATA_22 = 70521;
212 const int RUN_END_DATA_22 = 73930;
213 const int RUN_BEGIN_MC_22 = -73930;
214 const int RUN_END_MC_22 = -70521;
215
216 const int RUN_BEGIN_DATA_23 = 74031;
217 const int RUN_END_DATA_23 = 78536;
218 const int RUN_BEGIN_MC_23 = -78536;
219 const int RUN_END_MC_23 = -74031;
220
221 const int RUN_BEGIN_DATA_24 = 78615;
222 const int RUN_END_DATA_24 = 81094;
223 const int RUN_BEGIN_MC_24 = -81094;
224 const int RUN_END_MC_24 = -78615;
225
226 const int RUN_BEGIN_DATA_OffRes_24 = 81095;
227 const int RUN_END_DATA_OffRes_24 = 81631;
228 const int RUN_BEGIN_MC_OffRes_24 = -81631;
229 const int RUN_END_MC_OffRes_24 = -81095;
230
231 if ( run_no >= RUN_BEGIN_DATA_10 && run_no <= RUN_END_DATA_10 ) return 0;
232 else if ( run_no >= RUN_BEGIN_DATA_11 && run_no <= RUN_END_DATA_11 ) return 1;
233 else if ( run_no >= RUN_BEGIN_MC_10 && run_no <= RUN_END_MC_10 ) return 2;
234 else if ( run_no >= RUN_BEGIN_MC_11 && run_no <= RUN_END_MC_11 ) return 3;
235
236 else if ( run_no >= RUN_BEGIN_DATA_22 && run_no <= RUN_END_DATA_22 ) return 1;
237 else if ( run_no >= RUN_BEGIN_MC_22 && run_no <= RUN_END_MC_22 ) return 3;
238 else if ( run_no >= RUN_BEGIN_DATA_23 && run_no <= RUN_END_DATA_23 ) return 1;
239 else if ( run_no >= RUN_BEGIN_MC_23 && run_no <= RUN_END_MC_23 ) return 3;
240 else if ( run_no >= RUN_BEGIN_DATA_24 && run_no <= RUN_END_DATA_24 ) return 1;
241 else if ( run_no >= RUN_BEGIN_MC_24 && run_no <= RUN_END_MC_24 ) return 3;
242 else if ( run_no >= RUN_BEGIN_DATA_OffRes_24 && run_no <= RUN_END_DATA_OffRes_24 ) return 1;
243 else if ( run_no >= RUN_BEGIN_MC_OffRes_24 && run_no <= RUN_END_MC_OffRes_24 ) return 3;
244 else return -1;
245}
246
247void SimplePIDSvc::loadTOFInfo( EvtRecTrack* track ) {
248 // Initialization
249 for ( unsigned int i = 0; i < 8; i++ )
250 {
251 for ( unsigned int j = 0; j < 5; j++ ) m_tof_dt[i][j] = -99.;
252 m_tof_ph[i] = -99.;
253 }
254 for ( unsigned int i = 0; i < 2; i++ )
255 {
256 m_tof_zr[i] = -9999.;
257 m_tof_counter[i] = -1;
258 }
259 for ( unsigned int i = 0; i < 5; i++ ) { m_tof_chi[i] = -99.; }
260 m_tof_barrel = -1;
261
262 if ( !track->isExtTrackValid() || !track->isTofTrackValid() ) return;
263
264 SmartRefVector<RecTofTrack> tofTrk = track->tofTrack();
265 SmartRefVector<RecTofTrack>::iterator it;
266 RecExtTrack* extTrk = track->extTrack();
267 double zrhit[2];
268 zrhit[0] = extTrk->tof1Position().z();
269 zrhit[1] = extTrk->tof2Position().z();
270
271 TofHitStatus* hitst = new TofHitStatus;
272
273 for ( it = tofTrk.begin(); it != tofTrk.end(); it++ )
274 {
275 unsigned int st = ( *it )->status();
276 hitst->setStatus( st );
277 if ( hitst->is_raw() ) continue; // empty TOF hit
278 bool barrel = hitst->is_barrel();
279 bool readout = hitst->is_readout();
280 bool counter = hitst->is_counter();
281 bool cluster = hitst->is_cluster();
282 int layer = hitst->layer();
283 double tof = ( *it )->tof();
284 double ph = ( *it )->ph();
285 m_tof_counter[layer - 1] = ( *it )->tofID();
286
287 if ( barrel ) { m_tof_barrel = 1; }
288 else
289 {
290 m_tof_barrel = 0;
291 zrhit[0] = extTrk->tof1Position().rho();
292 }
293 m_tof_zr[0] = zrhit[0];
294 m_tof_zr[1] = zrhit[1];
295
296 int idx = -1;
297 if ( readout )
298 {
299 if ( barrel ) idx = ( ( st & 0xC0 ) >> 5 ) + ( ( ( st ^ 0x20 ) & 0x20 ) >> 5 ) - 2;
300 else idx = 7;
301 }
302 else if ( counter ) { idx = layer + 3; }
303 else if ( cluster ) { idx = 6; }
304 if ( idx == -1 ) continue;
305 m_tof_ph[idx] = ph;
306 for ( unsigned int i = 0; i < 5; i++ )
307 {
308 double offset = ( *it )->toffset( i );
309 double texp = ( *it )->texp( i );
310 if ( texp < 0.0 ) continue;
311 double dt = tof - offset - texp;
312 m_tof_dt[idx][i] = dt;
313 }
314 }
315 delete hitst;
316}
317
318void SimplePIDSvc::tofBarrelCorrection() {
319 const double EPS = 1e-4;
320 const double BG_LOW = 0.20;
321 const double BG_UP = 7.40;
322 const double COST_LOW = -0.81;
323 const double COST_UP = 0.81;
324 const double Q_LOW = 0.;
325 const double Q_UP = 9000.;
326 const double P_LOW = 0.2;
327 const double P_UP = 1.3;
328 const int BIN_BG = 15;
329 const int BIN_COST = 15;
330 const int BIN_P = 15;
331 double BG[BIN_BG + 1] = { 0.20, 0.87, 1.11, 1.35, 1.55, 1.72, 1.91, 2.17,
332 2.63, 3.05, 3.47, 3.93, 4.50, 5.27, 6.00, 7.40 };
333 double COST[BIN_COST + 1] = { -0.81, -0.64, -0.53, -0.43, -0.33, -0.23, -0.13, -0.04,
334 0.05, 0.14, 0.24, 0.34, 0.44, 0.54, 0.65, 0.81 };
335 double P[BIN_P + 1] = { 0.20, 0.47, 0.56, 0.65, 0.72, 0.79, 0.86, 0.92,
336 0.98, 1.03, 1.08, 1.13, 1.17, 1.22, 1.26, 1.30 };
337 int idx = getRunIdx( m_run );
338
339 if ( idx != -1 )
340 {
341 for ( unsigned int i = 0; i < 4; i++ )
342 { // only correct e, pi, K
343 double bg;
344 int bin_bg, bin_cost, bin_wgt;
345 int pid;
346 if ( i == 0 )
347 {
348 bg = max( P_LOW + EPS, min( P_UP - EPS, m_p[i] ) );
349 bin_bg = findBin( P, BIN_P, bg );
350 pid = 0;
351 }
352 else if ( i == 2 || i == 3 )
353 {
354 bg = max( BG_LOW + EPS, min( BG_UP - EPS, m_betagamma[i] ) );
355 bin_bg = findBin( BG, BIN_BG, bg );
356 pid = 1;
357 }
358 else { continue; }
359 double cost = m_cost[i];
360 int charge = m_charge[i];
361 double t[5], q;
362 double offset, sigma;
363 double offset_q, offset_bgcost;
364 int flag[4] = {
365 0,
366 0,
367 0,
368 0,
369 };
370 cost = max( COST_LOW + EPS, min( COST_UP - EPS, cost ) );
371 bin_cost = findBin( COST, BIN_COST, cost );
372 if ( bin_bg == -1 || bin_cost == -1 ) continue;
373
374 // corrections
375 for ( unsigned int j = 0; j < 4; j++ )
376 {
377 t[j] = m_tof_dt[j][i];
378 if ( fabs( t[j] + 99. ) < EPS ) // no readout
379 flag[j] = 0;
380 else flag[j] = 1;
381 q = m_tof_ph[j];
382 q = max( Q_LOW + EPS, min( Q_UP - EPS, q ) );
383 if ( charge == 1 )
384 {
385 offset_q = h_tof_p_q_offset[pid][idx][j]->Interpolate( q );
386 offset_bgcost = h_tof_p_bgcost_offset[pid][idx][j]->Interpolate( bg, cost );
387 t[j] = t[j] - offset_q - offset_bgcost;
388 }
389 else
390 {
391 offset_q = h_tof_m_q_offset[pid][idx][j]->Interpolate( q );
392 offset_bgcost = h_tof_m_bgcost_offset[pid][idx][j]->Interpolate( bg, cost );
393 t[j] = t[j] - offset_q - offset_bgcost;
394 }
395 }
396 bin_wgt = flag[0] * 8 + flag[1] * 4 + flag[2] * 2 + flag[3] - 1;
397 if ( bin_wgt == -1 ) continue;
398 t[4] = 0;
399 for ( unsigned int j = 0; j < 4; j++ )
400 {
401 if ( charge == 1 )
402 t[4] += t[j] *
403 h_tof_p_wgt[pid][idx][bin_wgt][j]->GetBinContent( bin_bg + 1, bin_cost + 1 );
404 else
405 t[4] += t[j] *
406 h_tof_m_wgt[pid][idx][bin_wgt][j]->GetBinContent( bin_bg + 1, bin_cost + 1 );
407 }
408 if ( charge == 1 )
409 {
410 t[4] /= h_tof_p_wgt[pid][idx][bin_wgt][4]->GetBinContent( bin_bg + 1, bin_cost + 1 );
411 offset = h_tof_p_final_offset[pid][idx][bin_wgt]->Interpolate( bg, cost );
412 sigma = h_tof_p_final_sigma[pid][idx][bin_wgt]->Interpolate( bg, cost );
413 m_tof_chi[i] = ( t[4] - offset ) / sigma;
414 }
415 else
416 {
417 t[4] /= h_tof_m_wgt[pid][idx][bin_wgt][4]->GetBinContent( bin_bg + 1, bin_cost + 1 );
418 offset = h_tof_m_final_offset[pid][idx][bin_wgt]->Interpolate( bg, cost );
419 sigma = h_tof_m_final_sigma[pid][idx][bin_wgt]->Interpolate( bg, cost );
420 m_tof_chi[i] = ( t[4] - offset ) / sigma;
421 }
422 }
423 }
424}
425
426void SimplePIDSvc::tofEndcapCorrection() {
427 const double EPS = 1e-4;
428 const double BG_LOW = 0.30;
429 const double BG_UP = 7.40;
430 const double Q_LOW = 0.;
431 const double Q_UP = 6000.;
432 const double COST_EAST_LOW = 0.720;
433 const double COST_EAST_UP = 0.930;
434 const double COST_WEST_LOW = -0.930;
435 const double COST_WEST_UP = -0.720;
436 const double P_LOW = 0.2;
437 const double P_UP = 1.3;
438
439 int idx = getRunIdx( m_run );
440
441 if ( idx != -1 )
442 {
443 for ( unsigned int i = 0; i < 4; i++ )
444 { // only correct e, pi, K
445 int pid;
446 double bg;
447 if ( i == 0 )
448 {
449 bg = max( P_LOW + EPS, min( P_UP - EPS, m_p[i] ) );
450 pid = 0;
451 }
452 else if ( i == 2 || i == 3 )
453 {
454 bg = max( BG_LOW + EPS, min( BG_UP - EPS, m_betagamma[i] ) );
455 pid = 1;
456 }
457 else { continue; }
458
459 int flag; // 0:east, 1:west
460 double cost = m_cost[i];
461 int charge = m_charge[i];
462 double t = m_tof_dt[7][i];
463 double q = m_tof_ph[7];
464 double off_q, off_bg, off_cost;
465 double sg_q, sg_bg, sg_cost;
466 if ( cost > 0 )
467 {
468 flag = 0;
469 cost = max( COST_EAST_LOW + EPS, min( COST_EAST_UP - EPS, cost ) );
470 }
471 else
472 {
473 flag = 1;
474 cost = max( COST_WEST_LOW + EPS, min( COST_WEST_UP - EPS, cost ) );
475 }
476 q = max( Q_LOW + EPS, min( Q_UP - EPS, q ) );
477
478 // corrections
479 if ( charge == 1 )
480 {
481 off_q = h_tofec_p_q_offset[pid][idx][flag]->Interpolate( q );
482 sg_q = h_tofec_p_q_sigma[pid][idx][flag]->Interpolate( q );
483 off_bg = h_tofec_p_bg_offset[pid][idx][flag]->Interpolate( bg );
484 sg_bg = h_tofec_p_bg_sigma[pid][idx][flag]->Interpolate( bg );
485 off_cost = h_tofec_p_cost_offset[pid][idx][flag]->Interpolate( cost );
486 sg_cost = h_tofec_p_cost_sigma[pid][idx][flag]->Interpolate( cost );
487 m_tof_chi[i] = ( ( ( t - off_q ) / sg_q - off_bg ) / sg_bg - off_cost ) / sg_cost;
488 }
489 else
490 {
491 off_q = h_tofec_m_q_offset[pid][idx][flag]->Interpolate( q );
492 sg_q = h_tofec_m_q_sigma[pid][idx][flag]->Interpolate( q );
493 off_bg = h_tofec_m_bg_offset[pid][idx][flag]->Interpolate( bg );
494 sg_bg = h_tofec_m_bg_sigma[pid][idx][flag]->Interpolate( bg );
495 off_cost = h_tofec_m_cost_offset[pid][idx][flag]->Interpolate( cost );
496 sg_cost = h_tofec_m_cost_sigma[pid][idx][flag]->Interpolate( cost );
497 m_tof_chi[i] = ( ( ( t - off_q ) / sg_q - off_bg ) / sg_bg - off_cost ) / sg_cost;
498 }
499 }
500 }
501}
502
503void SimplePIDSvc::loadDedxInfo( EvtRecTrack* track ) {
504 if ( track->isMdcDedxValid() )
505 {
506 RecMdcDedx* dedx_trk = track->mdcDedx();
507 for ( unsigned int i = 0; i < 5; i++ ) m_dedx_chi[i] = dedx_trk->chi( i );
508 }
509 else
510 {
511 for ( unsigned int i = 0; i < 5; i++ ) m_dedx_chi[i] = -99;
512 }
513}
514
515void SimplePIDSvc::dedxCorrection() {
516 int idx = getRunIdx( m_run );
517 const double EPS = 1e-4;
518 const double BG_LOW = 0.20;
519 const double BG_UP = 7.40;
520 const double COST_LOW = -0.93;
521 const double COST_UP = 0.93;
522 const double P_LOW = 0.2;
523 const double P_UP = 1.3;
524 if ( idx != -1 )
525 {
526 double offset, sigma;
527 for ( unsigned int i = 0; i < 4; i++ )
528 { // only correct e, pi, K
529 double bg;
530 int pid;
531 if ( i == 0 )
532 {
533 bg = max( P_LOW + EPS, min( P_UP - EPS, m_p[i] ) );
534 pid = 0;
535 }
536 else if ( i == 2 || i == 3 )
537 {
538 bg = max( BG_LOW + EPS, min( BG_UP - EPS, m_betagamma[i] ) );
539 pid = 1;
540 }
541 else { continue; }
542 double cost = m_cost[i];
543 double charge = m_charge[i];
544 cost = max( COST_LOW + EPS, min( COST_UP - EPS, cost ) );
545 if ( charge == 1 )
546 {
547 offset = h_dedx_p_offset[pid][idx]->Interpolate( bg, cost );
548 sigma = h_dedx_p_sigma[pid][idx]->Interpolate( bg, cost );
549 m_dedx_chi[i] = ( m_dedx_chi[i] - offset ) / sigma;
550 }
551 else
552 {
553 offset = h_dedx_m_offset[pid][idx]->Interpolate( bg, cost );
554 sigma = h_dedx_m_sigma[pid][idx]->Interpolate( bg, cost );
555 m_dedx_chi[i] = ( m_dedx_chi[i] - offset ) / sigma;
556 }
557 }
558 }
559}
560
561void SimplePIDSvc::tofBarrelSecondCorrection() {
562
563 int idx = getRunIdx( m_run );
564 const double EPS = 1e-4;
565 const double P_LOW = 0.0;
566 const double P_UP = 1.3;
567 const int BIN_P = 10;
568
569 const int RUN_BEGIN_DATA_10 = 11414;
570 const int RUN_END_DATA_10 = 14604;
571 const int RUN_BEGIN_MC_10 = -14604;
572 const int RUN_END_MC_10 = -11414;
573 const int RUN_BEGIN_DATA_11 = 20448;
574 const int RUN_END_DATA_11 = 23454;
575 const int RUN_BEGIN_MC_11 = -23454;
576 const int RUN_END_MC_11 = -20448;
577
578 const int RUN_BEGIN_DATA_22 = 70521;
579 const int RUN_END_DATA_22 = 73930;
580 const int RUN_BEGIN_MC_22 = -73930;
581 const int RUN_END_MC_22 = -70521;
582
583 const int RUN_BEGIN_DATA_23 = 74031;
584 const int RUN_END_DATA_23 = 78536;
585 const int RUN_BEGIN_MC_23 = -78536;
586 const int RUN_END_MC_23 = -74031;
587
588 const int RUN_BEGIN_DATA_24 = 78615;
589 const int RUN_END_DATA_24 = 81094;
590 const int RUN_BEGIN_MC_24 = -81094;
591 const int RUN_END_MC_24 = -78615;
592
593 const int RUN_BEGIN_DATA_OffRes_24 = 81095;
594 const int RUN_END_DATA_OffRes_24 = 81631;
595 const int RUN_BEGIN_MC_OffRes_24 = -81631;
596 const int RUN_END_MC_OffRes_24 = -81095;
597
598 double P[BIN_P + 1] = { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.3 };
599 if ( idx != -1 )
600 {
601
602 for ( unsigned int i = 2; i < 4; i++ )
603 { // second correct for tof of pi k
604
605 int aa = 99, bb = 99;
606
607 int bin_p;
608 double ptk = m_p[i];
609 ptk = max( P_LOW + EPS, min( P_UP - EPS, m_p[i] ) );
610 bin_p = findBin( P, BIN_P, ptk );
611
612 if ( i == 2 )
613 { // pi
614
615 if ( ( m_run >= RUN_BEGIN_DATA_10 && m_run <= RUN_END_DATA_10 && idx == 0 ) ||
616 ( m_run >= RUN_BEGIN_DATA_11 && m_run <= RUN_END_DATA_11 && idx == 1 ) )
617 {
618 aa = 2;
619 bb = 2;
620 } // round0304 data
621
622 if ( ( m_run >= RUN_BEGIN_MC_10 && m_run <= RUN_END_MC_10 && idx == 2 ) ||
623 ( m_run >= RUN_BEGIN_MC_11 && m_run <= RUN_END_MC_11 && idx == 3 ) )
624 {
625 aa = 2;
626 bb = 3;
627 } // round0304 inc
628
629 if ( m_run >= RUN_BEGIN_DATA_22 && m_run <= RUN_END_DATA_22 && idx == 1 )
630 {
631 aa = 3;
632 bb = 2;
633 } // round15 data
634
635 if ( m_run >= RUN_BEGIN_MC_22 && m_run <= RUN_END_MC_22 && idx == 3 )
636 {
637 aa = 3;
638 bb = 3;
639 } // round15 inc
640
641 if ( m_run >= RUN_BEGIN_DATA_23 && m_run <= RUN_END_DATA_23 && idx == 1 )
642 {
643 aa = 7;
644 bb = 2;
645 } // round16 data added by Yijia Zeng
646
647 if ( m_run >= RUN_BEGIN_MC_23 && m_run <= RUN_END_MC_23 && idx == 3 )
648 {
649 aa = 7;
650 bb = 3;
651 } // round16 inc added by Yijia Zeng
652
653 if ( m_run >= RUN_BEGIN_DATA_24 && m_run <= RUN_END_DATA_24 && idx == 1 )
654 {
655 aa = 10;
656 bb = 2;
657 } // round17 data added by Yijia Zeng
658
659 if ( m_run >= RUN_BEGIN_MC_24 && m_run <= RUN_END_MC_24 && idx == 3 )
660 {
661 aa = 10;
662 bb = 3;
663 } // round17 inc added by Yijia Zeng
664
665 if ( m_run >= RUN_BEGIN_DATA_OffRes_24 && m_run <= RUN_END_DATA_OffRes_24 && idx == 1 )
666 {
667 aa = 10;
668 bb = 2;
669 } // psipp off-resonance data added by Yijia Zeng
670
671 if ( m_run >= RUN_BEGIN_MC_OffRes_24 && m_run <= RUN_END_MC_OffRes_24 && idx == 3 )
672 {
673 aa = 10;
674 bb = 3;
675 } // psipp off-resonance inc added by Yijia Zeng
676 }
677
678 if ( i == 3 )
679 { // k
680
681 if ( ( m_run >= RUN_BEGIN_DATA_10 && m_run <= RUN_END_DATA_10 && idx == 0 ) ||
682 ( m_run >= RUN_BEGIN_DATA_11 && m_run <= RUN_END_DATA_11 && idx == 1 ) )
683 {
684 aa = 2;
685 bb = 0;
686 }
687
688 if ( ( m_run >= RUN_BEGIN_MC_10 && m_run <= RUN_END_MC_10 && idx == 2 ) ||
689 ( m_run >= RUN_BEGIN_MC_11 && m_run <= RUN_END_MC_11 && idx == 3 ) )
690 {
691 aa = 2;
692 bb = 1;
693 }
694
695 if ( m_run >= RUN_BEGIN_DATA_22 && m_run <= RUN_END_DATA_22 && idx == 1 )
696 {
697 aa = 3;
698 bb = 0;
699 }
700
701 if ( m_run >= RUN_BEGIN_MC_22 && m_run <= RUN_END_MC_22 && idx == 3 )
702 {
703 aa = 3;
704 bb = 1;
705 }
706
707 if ( m_run >= RUN_BEGIN_DATA_23 && m_run <= RUN_END_DATA_23 && idx == 1 )
708 {
709 aa = 7;
710 bb = 0;
711 } // round16 data added by Yijia Zeng
712
713 if ( m_run >= RUN_BEGIN_MC_23 && m_run <= RUN_END_MC_23 && idx == 3 )
714 {
715 aa = 7;
716 bb = 1;
717 } // round16 inc added by Yijia Zeng
718
719 if ( m_run >= RUN_BEGIN_DATA_24 && m_run <= RUN_END_DATA_24 && idx == 1 )
720 {
721 aa = 10;
722 bb = 0;
723 } // round17 data added by Yijia Zeng
724
725 if ( m_run >= RUN_BEGIN_MC_24 && m_run <= RUN_END_MC_24 && idx == 3 )
726 {
727 aa = 10;
728 bb = 1;
729 } // round17 inc added by Yijia Zeng
730
731 if ( m_run >= RUN_BEGIN_DATA_OffRes_24 && m_run <= RUN_END_DATA_OffRes_24 && idx == 1 )
732 {
733 aa = 10;
734 bb = 0;
735 } // psipp off-resonance data added by Yijia Zeng
736
737 if ( m_run >= RUN_BEGIN_MC_OffRes_24 && m_run <= RUN_END_MC_OffRes_24 && idx == 3 )
738 {
739 aa = 10;
740 bb = 1;
741 } // psipp off-resonance inc added by Yijia Zeng
742 }
743
744 if ( m_tof_chi[i] != -99 && m_run > 0 )
745 {
746
747 // NO Correction on the higher FIVE momentum region from Yijia for round 16 & 17 data
748 // sample NO Correction on the higher SIX momentum region from Kaikai He for round03-04
749 // and round15 The Corrected Distribution is somehow worse than the Original one
750 if ( !( ( m_run >= RUN_BEGIN_DATA_23 && m_run <= RUN_END_DATA_23 ) ||
751 ( m_run >= RUN_BEGIN_DATA_24 && m_run <= RUN_END_DATA_24 ) ||
752 ( m_run >= RUN_BEGIN_DATA_OffRes_24 && m_run <= RUN_END_DATA_OffRes_24 ) ) )
753 {
754 if ( bin_p < 4 )
755 {
756 m_tof_chi[i] =
757 ( m_tof_chi[i] - m_gaussion_mean[aa][bb][bin_p] ) /
758 ( m_gaussion_sigmab[aa][bb][bin_p] / m_gaussion_sigmab[aa][bb + 1][bin_p] );
759 }
760 if ( bin_p >= 4 )
761 { m_tof_chi[i] = ( m_tof_chi[i] - m_gaussion_mean[aa][bb][bin_p] ); }
762 }
763 else
764 {
765 if ( bin_p < 5 )
766 {
767 m_tof_chi[i] =
768 ( m_tof_chi[i] - m_gaussion_mean[aa][bb][bin_p] ) /
769 ( m_gaussion_sigmab[aa][bb][bin_p] / m_gaussion_sigmab[aa][bb + 1][bin_p] );
770 }
771 if ( bin_p >= 5 )
772 { m_tof_chi[i] = ( m_tof_chi[i] - m_gaussion_mean[aa][bb][bin_p] ); }
773 }
774 }
775 if ( m_tof_chi[i] != -99 && m_run < 0 )
776 {
777 m_tof_chi[i] = ( m_tof_chi[i] - m_gaussion_mean[aa][bb][bin_p] ) /
778 ( m_gaussion_sigmab[aa][bb][bin_p] / m_gaussion_sigmab[aa][bb][bin_p] );
779 }
780 }
781
782 } // end idx!= -1
783}
784
785void SimplePIDSvc::tofEndcapSecondCorrection() {
786
787 int idx = getRunIdx( m_run );
788 const double EPS = 1e-4;
789 const double P_LOW = 0.0;
790 const double P_UP = 1.3;
791 const int BIN_P = 10;
792
793 const int RUN_BEGIN_DATA_10 = 11414;
794 const int RUN_END_DATA_10 = 14604;
795 const int RUN_BEGIN_MC_10 = -14604;
796 const int RUN_END_MC_10 = -11414;
797 const int RUN_BEGIN_DATA_11 = 20448;
798 const int RUN_END_DATA_11 = 23454;
799 const int RUN_BEGIN_MC_11 = -23454;
800 const int RUN_END_MC_11 = -20448;
801
802 const int RUN_BEGIN_DATA_22 = 70521;
803 const int RUN_END_DATA_22 = 73930;
804 const int RUN_BEGIN_MC_22 = -73930;
805 const int RUN_END_MC_22 = -70521;
806
807 const int RUN_BEGIN_DATA_23 = 74031;
808 const int RUN_END_DATA_23 = 78536;
809 const int RUN_BEGIN_MC_23 = -78536;
810 const int RUN_END_MC_23 = -74031;
811
812 const int RUN_BEGIN_DATA_24 = 78615;
813 const int RUN_END_DATA_24 = 81094;
814 const int RUN_BEGIN_MC_24 = -81094;
815 const int RUN_END_MC_24 = -78615;
816
817 const int RUN_BEGIN_DATA_OffRes_24 = 81095;
818 const int RUN_END_DATA_OffRes_24 = 81631;
819 const int RUN_BEGIN_MC_OffRes_24 = -81631;
820 const int RUN_END_MC_OffRes_24 = -81095;
821
822 double P[BIN_P + 1] = { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.3 };
823 if ( idx != -1 )
824 {
825
826 for ( unsigned int i = 2; i < 4; i++ )
827 { // second correct for tof of pi k
828
829 int aa = 99, bb = 99;
830
831 int bin_p;
832 double ptk = m_p[i];
833 ptk = max( P_LOW + EPS, min( P_UP - EPS, m_p[i] ) );
834 bin_p = findBin( P, BIN_P, ptk );
835
836 if ( i == 2 )
837 { // pi
838
839 if ( ( m_run >= RUN_BEGIN_DATA_10 && m_run <= RUN_END_DATA_10 && idx == 0 ) ||
840 ( m_run >= RUN_BEGIN_DATA_11 && m_run <= RUN_END_DATA_11 && idx == 1 ) )
841 {
842 aa = 4;
843 bb = 2;
844 } // round0304 data endcap
845
846 if ( ( m_run >= RUN_BEGIN_MC_10 && m_run <= RUN_END_MC_10 && idx == 2 ) ||
847 ( m_run >= RUN_BEGIN_MC_11 && m_run <= RUN_END_MC_11 && idx == 3 ) )
848 {
849 aa = 4;
850 bb = 3;
851 } // round0304 inc endcap
852
853 if ( m_run >= RUN_BEGIN_DATA_22 && m_run <= RUN_END_DATA_22 && idx == 1 )
854 {
855 aa = 5;
856 bb = 2;
857 } // round15 data endcap
858
859 if ( m_run >= RUN_BEGIN_MC_22 && m_run <= RUN_END_MC_22 && idx == 3 )
860 {
861 aa = 5;
862 bb = 3;
863 } // round15 inc endcap
864
865 if ( m_run >= RUN_BEGIN_DATA_23 && m_run <= RUN_END_DATA_23 && idx == 1 )
866 {
867 aa = 8;
868 bb = 2;
869 } // round16 data endcap added by Yijia Zeng
870
871 if ( m_run >= RUN_BEGIN_MC_23 && m_run <= RUN_END_MC_23 && idx == 3 )
872 {
873 aa = 8;
874 bb = 3;
875 } // round16 inc endcap added by Yijia Zeng
876
877 if ( m_run >= RUN_BEGIN_DATA_24 && m_run <= RUN_END_DATA_24 && idx == 1 )
878 {
879 aa = 11;
880 bb = 2;
881 } // round17 data endcap added by Yijia Zeng
882
883 if ( m_run >= RUN_BEGIN_MC_24 && m_run <= RUN_END_MC_24 && idx == 3 )
884 {
885 aa = 11;
886 bb = 3;
887 } // round17 inc endcap added by Yijia Zeng
888
889 if ( m_run >= RUN_BEGIN_DATA_OffRes_24 && m_run <= RUN_END_DATA_OffRes_24 && idx == 1 )
890 {
891 aa = 11;
892 bb = 2;
893 } // psipp off-resonance data endcap added by Yijia Zeng
894
895 if ( m_run >= RUN_BEGIN_MC_OffRes_24 && m_run <= RUN_END_MC_OffRes_24 && idx == 3 )
896 {
897 aa = 11;
898 bb = 3;
899 } // psipp off-resonance inc endcap added by Yijia Zeng
900 }
901
902 if ( i == 3 )
903 { // k
904
905 if ( ( m_run >= RUN_BEGIN_DATA_10 && m_run <= RUN_END_DATA_10 && idx == 0 ) ||
906 ( m_run >= RUN_BEGIN_DATA_11 && m_run <= RUN_END_DATA_11 && idx == 1 ) )
907 {
908 aa = 4;
909 bb = 0;
910 }
911
912 if ( ( m_run >= RUN_BEGIN_MC_10 && m_run <= RUN_END_MC_10 && idx == 2 ) ||
913 ( m_run >= RUN_BEGIN_MC_11 && m_run <= RUN_END_MC_11 && idx == 3 ) )
914 {
915 aa = 4;
916 bb = 1;
917 }
918
919 if ( m_run >= RUN_BEGIN_DATA_22 && m_run <= RUN_END_DATA_22 && idx == 1 )
920 {
921 aa = 5;
922 bb = 0;
923 }
924
925 if ( m_run >= RUN_BEGIN_MC_22 && m_run <= RUN_END_MC_22 && idx == 3 )
926 {
927 aa = 5;
928 bb = 1;
929 }
930
931 if ( m_run >= RUN_BEGIN_DATA_23 && m_run <= RUN_END_DATA_23 && idx == 1 )
932 {
933 aa = 8;
934 bb = 0;
935 } // round16 data endcap added by Yijia Zeng
936
937 if ( m_run >= RUN_BEGIN_MC_23 && m_run <= RUN_END_MC_23 && idx == 3 )
938 {
939 aa = 8;
940 bb = 1;
941 } // round16 inc endcap added by Yijia Zeng
942
943 if ( m_run >= RUN_BEGIN_DATA_24 && m_run <= RUN_END_DATA_24 && idx == 1 )
944 {
945 aa = 11;
946 bb = 0;
947 } // round17 data endcap added by Yijia Zeng
948
949 if ( m_run >= RUN_BEGIN_MC_24 && m_run <= RUN_END_MC_24 && idx == 3 )
950 {
951 aa = 11;
952 bb = 1;
953 } // round17 inc endcap added by Yijia Zeng
954
955 if ( m_run >= RUN_BEGIN_DATA_OffRes_24 && m_run <= RUN_END_DATA_OffRes_24 && idx == 1 )
956 {
957 aa = 11;
958 bb = 0;
959 } // psipp off-resonance data endcap added by Yijia Zeng
960
961 if ( m_run >= RUN_BEGIN_MC_OffRes_24 && m_run <= RUN_END_MC_OffRes_24 && idx == 3 )
962 {
963 aa = 11;
964 bb = 1;
965 } // psipp off-resonance inc endcap added by Yijia Zeng
966 }
967
968 if ( m_tof_chi[i] != -99 && ( aa == 5 || aa == 8 || aa == 11 ) )
969 { // for round15 and round16 and round17 data and inc
970
971 m_tof_chi[i] = ( m_tof_chi[i] - m_gaussion_mean[aa][bb][bin_p] ) / ( 60. / 110. );
972 }
973
974 if ( m_tof_chi[i] != -99 && aa == 4 )
975 { // for round0304 data and inc
976
977 m_tof_chi[i] = ( m_tof_chi[i] - m_gaussion_mean[aa][bb][bin_p] );
978 }
979 }
980
981 } // end idx!= -1
982}
983
984void SimplePIDSvc::dedxSecondCorrection()
985
986{
987
988 int idx = getRunIdx( m_run );
989 const double EPS = 1e-4;
990 const double P_LOW = 0.0;
991 const double P_UP = 1.3;
992 const int BIN_P = 10;
993
994 const int RUN_BEGIN_DATA_10 = 11414;
995 const int RUN_END_DATA_10 = 14604;
996 const int RUN_BEGIN_MC_10 = -14604;
997 const int RUN_END_MC_10 = -11414;
998 const int RUN_BEGIN_DATA_11 = 20448;
999 const int RUN_END_DATA_11 = 23454;
1000 const int RUN_BEGIN_MC_11 = -23454;
1001 const int RUN_END_MC_11 = -20448;
1002
1003 const int RUN_BEGIN_DATA_22 = 70521;
1004 const int RUN_END_DATA_22 = 73930;
1005 const int RUN_BEGIN_MC_22 = -73930;
1006 const int RUN_END_MC_22 = -70521;
1007
1008 const int RUN_BEGIN_DATA_23 = 74031;
1009 const int RUN_END_DATA_23 = 78536;
1010 const int RUN_BEGIN_MC_23 = -78536;
1011 const int RUN_END_MC_23 = -74031;
1012
1013 const int RUN_BEGIN_DATA_24 = 78615;
1014 const int RUN_END_DATA_24 = 81094;
1015 const int RUN_BEGIN_MC_24 = -81094;
1016 const int RUN_END_MC_24 = -78615;
1017
1018 const int RUN_BEGIN_DATA_OffRes_24 = 81095;
1019 const int RUN_END_DATA_OffRes_24 = 81631;
1020 const int RUN_BEGIN_MC_OffRes_24 = -81631;
1021 const int RUN_END_MC_OffRes_24 = -81095;
1022
1023 double P[BIN_P + 1] = { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.3 };
1024 if ( idx != -1 )
1025 {
1026
1027 for ( unsigned int i = 2; i < 4; i++ )
1028 { // second correct for dedx of pi k
1029
1030 int aa = 99, bb = 99;
1031
1032 int bin_p;
1033 double ptk = m_p[i];
1034 ptk = max( P_LOW + EPS, min( P_UP - EPS, m_p[i] ) );
1035 bin_p = findBin( P, BIN_P, ptk );
1036
1037 if ( i == 2 )
1038 { // pi
1039
1040 if ( ( m_run >= RUN_BEGIN_DATA_10 && m_run <= RUN_END_DATA_10 && idx == 0 ) ||
1041 ( m_run >= RUN_BEGIN_DATA_11 && m_run <= RUN_END_DATA_11 && idx == 1 ) )
1042 {
1043 aa = 0;
1044 bb = 2;
1045 } // round0304 data
1046
1047 if ( ( m_run >= RUN_BEGIN_MC_10 && m_run <= RUN_END_MC_10 && idx == 2 ) ||
1048 ( m_run >= RUN_BEGIN_MC_11 && m_run <= RUN_END_MC_11 && idx == 3 ) )
1049 {
1050 aa = 0;
1051 bb = 3;
1052 } // round0304 inc
1053
1054 if ( m_run >= RUN_BEGIN_DATA_22 && m_run <= RUN_END_DATA_22 && idx == 1 )
1055 {
1056 aa = 1;
1057 bb = 2;
1058 } // round15 data
1059
1060 if ( m_run >= RUN_BEGIN_MC_22 && m_run <= RUN_END_MC_22 && idx == 3 )
1061 {
1062 aa = 1;
1063 bb = 3;
1064 } // round15 inc
1065
1066 if ( m_run >= RUN_BEGIN_DATA_23 && m_run <= RUN_END_DATA_23 && idx == 1 )
1067 {
1068 aa = 6;
1069 bb = 2;
1070 } // round16 data added by Yijia Zeng
1071
1072 if ( m_run >= RUN_BEGIN_MC_23 && m_run <= RUN_END_MC_23 && idx == 3 )
1073 {
1074 aa = 6;
1075 bb = 3;
1076 } // round16 inc added by Yijia Zeng
1077
1078 if ( m_run >= RUN_BEGIN_DATA_24 && m_run <= RUN_END_DATA_24 && idx == 1 )
1079 {
1080 aa = 9;
1081 bb = 2;
1082 } // round17 data added by Yijia Zeng
1083
1084 if ( m_run >= RUN_BEGIN_MC_24 && m_run <= RUN_END_MC_24 && idx == 3 )
1085 {
1086 aa = 9;
1087 bb = 3;
1088 } // round17 inc added by Yijia Zeng
1089
1090 if ( m_run >= RUN_BEGIN_DATA_OffRes_24 && m_run <= RUN_END_DATA_OffRes_24 && idx == 1 )
1091 {
1092 aa = 9;
1093 bb = 2;
1094 } // psipp off-resonance data added by Yijia Zeng
1095
1096 if ( m_run >= RUN_BEGIN_MC_OffRes_24 && m_run <= RUN_END_MC_OffRes_24 && idx == 3 )
1097 {
1098 aa = 9;
1099 bb = 3;
1100 } // psipp off-resonance inc added by Yijia Zeng
1101 }
1102
1103 if ( i == 3 )
1104 { // k
1105
1106 if ( ( m_run >= RUN_BEGIN_DATA_10 && m_run <= RUN_END_DATA_10 && idx == 0 ) ||
1107 ( m_run >= RUN_BEGIN_DATA_11 && m_run <= RUN_END_DATA_11 && idx == 1 ) )
1108 {
1109 aa = 0;
1110 bb = 0;
1111 }
1112
1113 if ( ( m_run >= RUN_BEGIN_MC_10 && m_run <= RUN_END_MC_10 && idx == 2 ) ||
1114 ( m_run >= RUN_BEGIN_MC_11 && m_run <= RUN_END_MC_11 && idx == 3 ) )
1115 {
1116 aa = 0;
1117 bb = 1;
1118 }
1119
1120 if ( m_run >= RUN_BEGIN_DATA_22 && m_run <= RUN_END_DATA_22 && idx == 1 )
1121 {
1122 aa = 1;
1123 bb = 0;
1124 }
1125
1126 if ( m_run >= RUN_BEGIN_MC_22 && m_run <= RUN_END_MC_22 && idx == 3 )
1127 {
1128 aa = 1;
1129 bb = 1;
1130 }
1131
1132 if ( m_run >= RUN_BEGIN_DATA_23 && m_run <= RUN_END_DATA_23 && idx == 1 )
1133 {
1134 aa = 6;
1135 bb = 0;
1136 } // round16 data added by Yijia Zeng
1137
1138 if ( m_run >= RUN_BEGIN_MC_23 && m_run <= RUN_END_MC_23 && idx == 3 )
1139 {
1140 aa = 6;
1141 bb = 1;
1142 } // round16 inc added by Yijia Zeng
1143
1144 if ( m_run >= RUN_BEGIN_DATA_24 && m_run <= RUN_END_DATA_24 && idx == 1 )
1145 {
1146 aa = 9;
1147 bb = 0;
1148 } // round17 data added by Yijia Zeng
1149
1150 if ( m_run >= RUN_BEGIN_MC_24 && m_run <= RUN_END_MC_24 && idx == 3 )
1151 {
1152 aa = 9;
1153 bb = 1;
1154 } // round17 inc added by Yijia Zeng
1155
1156 if ( m_run >= RUN_BEGIN_DATA_OffRes_24 && m_run <= RUN_END_DATA_OffRes_24 && idx == 1 )
1157 {
1158 aa = 9;
1159 bb = 0;
1160 } // psipp off-resonance data added by Yijia Zeng
1161
1162 if ( m_run >= RUN_BEGIN_MC_OffRes_24 && m_run <= RUN_END_MC_OffRes_24 && idx == 3 )
1163 {
1164 aa = 9;
1165 bb = 1;
1166 } // psipp off-resonance inc added by Yijia Zeng
1167 }
1168
1169 if ( m_dedx_chi[i] != -99 && m_run > 0 )
1170 {
1171 m_dedx_chi[i] =
1172 ( m_dedx_chi[i] - m_gaussion_mean[aa][bb][bin_p] ) /
1173 ( m_gaussion_sigmab[aa][bb][bin_p] / m_gaussion_sigmab[aa][bb + 1][bin_p] );
1174 }
1175 if ( m_dedx_chi[i] != -99 && m_run < 0 )
1176 {
1177 m_dedx_chi[i] =
1178 ( m_dedx_chi[i] - m_gaussion_mean[aa][bb][bin_p] ) /
1179 ( m_gaussion_sigmab[aa][bb][bin_p] / m_gaussion_sigmab[aa][bb][bin_p] );
1180 }
1181 }
1182
1183 } // end idx!= -1
1184}
1185
1186void SimplePIDSvc::loadSecondPar() {
1187 string path = getenv( "SIMPLEPIDSVCROOT" );
1188
1189 for ( int i = 0; i < 12; i++ )
1190 {
1191
1192 const char* dir;
1193
1194 if ( i == 0 ) dir = "round0304_dedx";
1195 else if ( i == 1 ) dir = "round15_dedx";
1196 else if ( i == 2 ) dir = "round0304_tof";
1197 else if ( i == 3 ) dir = "round15_tof";
1198 else if ( i == 4 ) dir = "round0304_tof_endcap";
1199 else if ( i == 5 ) dir = "round15_tof_endcap";
1200 // Added by Yijia Zeng for round 16 data sample
1201 else if ( i == 6 ) dir = "round16_dedx";
1202 else if ( i == 7 ) dir = "round16_tof";
1203 else if ( i == 8 ) dir = "round16_tof_endcap";
1204 // Added by Yijia Zeng for round 17 data sample
1205 else if ( i == 9 ) dir = "round17_dedx";
1206 else if ( i == 10 ) dir = "round17_tof";
1207 else if ( i == 11 ) dir = "round17_tof_endcap";
1208
1209 else
1210 {
1211 cout << "Boundary Error! " << endl;
1212 exit( 1 );
1213 }
1214
1215 for ( int j = 0; j < 4; j++ )
1216 {
1217
1218 const char* name;
1219
1220 if ( j == 0 ) name = "data_K";
1221 else if ( j == 1 ) name = "inc_K";
1222 else if ( j == 2 ) name = "data_pi";
1223 else if ( j == 3 ) name = "inc_pi";
1224 else
1225 {
1226 cout << "Boundary Error! " << endl;
1227 exit( 1 );
1228 }
1229
1230 ifstream second_cor( path +
1231 Form( "/share/second_correct/%s/%s_%s.dat", dir, dir, name ) );
1232
1233 for ( int m = 0; m < 10; m++ )
1234 {
1235
1236 second_cor >> m_gaussion_mean[i][j][m] >> m_gaussion_sigma[i][j][m] >>
1237 m_gaussion_sigmab[i][j][m];
1238 }
1239
1240 second_cor.close();
1241 }
1242 }
1243}
1244
1245void SimplePIDSvc::loadHistogram() {
1246 string path = getenv( "SIMPLEPIDSVCROOT" );
1247 vector<string> filename;
1248 for ( unsigned int idx = 0; idx < 2; idx++ )
1249 {
1250 const char* dir;
1251 if ( idx == 0 ) dir = "electron";
1252 else if ( idx == 1 ) dir = "kpi";
1253 else
1254 {
1255 cout << "Boundary Error! " << endl;
1256 exit( 1 );
1257 }
1258
1259 // dedx
1260 filename.clear();
1261 filename.push_back( path + Form( "/share/%s/dedx/dedx_d10.root", dir ) );
1262 filename.push_back( path + Form( "/share/%s/dedx/dedx_d11.root", dir ) );
1263 filename.push_back( path + Form( "/share/%s/dedx/dedx_m10.root", dir ) );
1264 filename.push_back( path + Form( "/share/%s/dedx/dedx_m11.root", dir ) );
1265 for ( unsigned int i = 0; i < filename.size(); i++ )
1266 {
1267 f_dedx[idx][i] = new TFile( filename[i].c_str(), "READ" );
1268 const char* name;
1269 if ( i == 0 ) name = "d10";
1270 else if ( i == 1 ) name = "d11";
1271 else if ( i == 2 ) name = "m10";
1272 else if ( i == 3 ) name = "m11";
1273 else
1274 {
1275 cout << "Boundary Error! " << endl;
1276 exit( 1 );
1277 }
1278 h_dedx_p_offset[idx][i] =
1279 (TH2D*)f_dedx[idx][i]->Get( Form( "h_dedx_p_offset_%s", name ) );
1280 h_dedx_p_sigma[idx][i] = (TH2D*)f_dedx[idx][i]->Get( Form( "h_dedx_p_sigma_%s", name ) );
1281 h_dedx_m_offset[idx][i] =
1282 (TH2D*)f_dedx[idx][i]->Get( Form( "h_dedx_m_offset_%s", name ) );
1283 h_dedx_m_sigma[idx][i] = (TH2D*)f_dedx[idx][i]->Get( Form( "h_dedx_m_sigma_%s", name ) );
1284 }
1285 // tof_barrel q
1286 filename.clear();
1287 filename.push_back( path + Form( "/share/%s/tof_barrel/tof_q_d10.root", dir ) );
1288 filename.push_back( path + Form( "/share/%s/tof_barrel/tof_q_d11.root", dir ) );
1289 filename.push_back( path + Form( "/share/%s/tof_barrel/tof_q_m10.root", dir ) );
1290 filename.push_back( path + Form( "/share/%s/tof_barrel/tof_q_m11.root", dir ) );
1291 for ( unsigned int i = 0; i < filename.size(); i++ )
1292 {
1293 f_tof_q[idx][i] = new TFile( filename[i].c_str(), "READ" );
1294 const char* name;
1295 if ( i == 0 ) name = "d10";
1296 else if ( i == 1 ) name = "d11";
1297 else if ( i == 2 ) name = "m10";
1298 else if ( i == 3 ) name = "m11";
1299 else
1300 {
1301 cout << "Boundary Error! " << endl;
1302 exit( 1 );
1303 }
1304 for ( unsigned int j = 0; j < 4; j++ )
1305 {
1306 h_tof_p_q_offset[idx][i][j] =
1307 (TH1D*)f_tof_q[idx][i]->Get( Form( "h_tof_p_q_offset_%s_%d", name, j ) );
1308 h_tof_m_q_offset[idx][i][j] =
1309 (TH1D*)f_tof_q[idx][i]->Get( Form( "h_tof_m_q_offset_%s_%d", name, j ) );
1310 h_tof_p_q_sigma[idx][i][j] =
1311 (TH1D*)f_tof_q[idx][i]->Get( Form( "h_tof_p_q_sigma_%s_%d", name, j ) );
1312 h_tof_m_q_sigma[idx][i][j] =
1313 (TH1D*)f_tof_q[idx][i]->Get( Form( "h_tof_m_q_sigma_%s_%d", name, j ) );
1314 }
1315 }
1316 // tof_barrel bg&cost
1317 filename.clear();
1318 filename.push_back( path + Form( "/share/%s/tof_barrel/tof_bg_cost_d10.root", dir ) );
1319 filename.push_back( path + Form( "/share/%s/tof_barrel/tof_bg_cost_d11.root", dir ) );
1320 filename.push_back( path + Form( "/share/%s/tof_barrel/tof_bg_cost_m10.root", dir ) );
1321 filename.push_back( path + Form( "/share/%s/tof_barrel/tof_bg_cost_m11.root", dir ) );
1322 for ( unsigned int i = 0; i < filename.size(); i++ )
1323 {
1324 f_tof_bgcost[idx][i] = new TFile( filename[i].c_str(), "READ" );
1325 const char* name;
1326 if ( i == 0 ) name = "d10";
1327 else if ( i == 1 ) name = "d11";
1328 else if ( i == 2 ) name = "m10";
1329 else if ( i == 3 ) name = "m11";
1330 else
1331 {
1332 cout << "Boundary Error! " << endl;
1333 exit( 1 );
1334 }
1335 for ( unsigned int j = 0; j < 4; j++ )
1336 {
1337 h_tof_p_bgcost_offset[idx][i][j] =
1338 (TH2D*)f_tof_bgcost[idx][i]->Get( Form( "h_tof_p_bgcost_offset_%s_%d", name, j ) );
1339 h_tof_m_bgcost_offset[idx][i][j] =
1340 (TH2D*)f_tof_bgcost[idx][i]->Get( Form( "h_tof_m_bgcost_offset_%s_%d", name, j ) );
1341 h_tof_p_bgcost_sigma[idx][i][j] =
1342 (TH2D*)f_tof_bgcost[idx][i]->Get( Form( "h_tof_p_bgcost_sigma_%s_%d", name, j ) );
1343 h_tof_m_bgcost_sigma[idx][i][j] =
1344 (TH2D*)f_tof_bgcost[idx][i]->Get( Form( "h_tof_m_bgcost_sigma_%s_%d", name, j ) );
1345 }
1346 }
1347 // tof_barrel wgt
1348 filename.clear();
1349 filename.push_back( path + Form( "/share/%s/tof_barrel/tof_wgt_d10.root", dir ) );
1350 filename.push_back( path + Form( "/share/%s/tof_barrel/tof_wgt_d11.root", dir ) );
1351 filename.push_back( path + Form( "/share/%s/tof_barrel/tof_wgt_m10.root", dir ) );
1352 filename.push_back( path + Form( "/share/%s/tof_barrel/tof_wgt_m11.root", dir ) );
1353 for ( unsigned int i = 0; i < filename.size(); i++ )
1354 {
1355 f_tof_wgt[idx][i] = new TFile( filename[i].c_str(), "READ" );
1356 const char* name;
1357 if ( i == 0 ) name = "d10";
1358 else if ( i == 1 ) name = "d11";
1359 else if ( i == 2 ) name = "m10";
1360 else if ( i == 3 ) name = "m11";
1361 else
1362 {
1363 cout << "Boundary Error! " << endl;
1364 exit( 1 );
1365 }
1366 for ( unsigned int j = 0; j < 15; j++ )
1367 {
1368 for ( unsigned int k = 0; k < 5; k++ )
1369 {
1370 h_tof_p_wgt[idx][i][j][k] =
1371 (TH2D*)f_tof_wgt[idx][i]->Get( Form( "h_tof_p_wgt_%s_%d_%d", name, j, k ) );
1372 h_tof_m_wgt[idx][i][j][k] =
1373 (TH2D*)f_tof_wgt[idx][i]->Get( Form( "h_tof_m_wgt_%s_%d_%d", name, j, k ) );
1374 }
1375 }
1376 }
1377 // tof_barrel corr
1378 filename.clear();
1379 filename.push_back( path + Form( "/share/%s/tof_barrel/tof_final_d10.root", dir ) );
1380 filename.push_back( path + Form( "/share/%s/tof_barrel/tof_final_d11.root", dir ) );
1381 filename.push_back( path + Form( "/share/%s/tof_barrel/tof_final_m10.root", dir ) );
1382 filename.push_back( path + Form( "/share/%s/tof_barrel/tof_final_m11.root", dir ) );
1383 for ( unsigned int i = 0; i < filename.size(); i++ )
1384 {
1385 f_tof_final[idx][i] = new TFile( filename[i].c_str(), "READ" );
1386 const char* name;
1387 if ( i == 0 ) name = "d10";
1388 else if ( i == 1 ) name = "d11";
1389 else if ( i == 2 ) name = "m10";
1390 else if ( i == 3 ) name = "m11";
1391 else
1392 {
1393 cout << "Boundary Error! " << endl;
1394 exit( 1 );
1395 }
1396 for ( unsigned int j = 0; j < 15; j++ )
1397 {
1398 h_tof_p_final_offset[idx][i][j] =
1399 (TH2D*)f_tof_final[idx][i]->Get( Form( "h_tof_p_final_offset_%s_%d", name, j ) );
1400 h_tof_m_final_offset[idx][i][j] =
1401 (TH2D*)f_tof_final[idx][i]->Get( Form( "h_tof_m_final_offset_%s_%d", name, j ) );
1402 h_tof_p_final_sigma[idx][i][j] =
1403 (TH2D*)f_tof_final[idx][i]->Get( Form( "h_tof_p_final_sigma_%s_%d", name, j ) );
1404 h_tof_m_final_sigma[idx][i][j] =
1405 (TH2D*)f_tof_final[idx][i]->Get( Form( "h_tof_m_final_sigma_%s_%d", name, j ) );
1406 }
1407 }
1408 // tof_endcap q
1409 filename.clear();
1410 filename.push_back( path + Form( "/share/%s/tof_endcap/tofec_q_d10.root", dir ) );
1411 filename.push_back( path + Form( "/share/%s/tof_endcap/tofec_q_d11.root", dir ) );
1412 filename.push_back( path + Form( "/share/%s/tof_endcap/tofec_q_m10.root", dir ) );
1413 filename.push_back( path + Form( "/share/%s/tof_endcap/tofec_q_m11.root", dir ) );
1414 for ( unsigned int i = 0; i < filename.size(); i++ )
1415 {
1416 f_tofec_q[idx][i] = new TFile( filename[i].c_str(), "READ" );
1417 const char* name;
1418 if ( i == 0 ) name = "d10";
1419 else if ( i == 1 ) name = "d11";
1420 else if ( i == 2 ) name = "m10";
1421 else if ( i == 3 ) name = "m11";
1422 else
1423 {
1424 cout << "Boundary Error! " << endl;
1425 exit( 1 );
1426 }
1427 for ( unsigned int j = 0; j < 2; j++ )
1428 {
1429 h_tofec_p_q_offset[idx][i][j] =
1430 (TH1D*)f_tofec_q[idx][i]->Get( Form( "h_tofec_p_q_offset_%s_%d", name, j ) );
1431 h_tofec_m_q_offset[idx][i][j] =
1432 (TH1D*)f_tofec_q[idx][i]->Get( Form( "h_tofec_m_q_offset_%s_%d", name, j ) );
1433 h_tofec_p_q_sigma[idx][i][j] =
1434 (TH1D*)f_tofec_q[idx][i]->Get( Form( "h_tofec_p_q_sigma_%s_%d", name, j ) );
1435 h_tofec_m_q_sigma[idx][i][j] =
1436 (TH1D*)f_tofec_q[idx][i]->Get( Form( "h_tofec_m_q_sigma_%s_%d", name, j ) );
1437 }
1438 }
1439 // tof_endcap bg
1440 filename.clear();
1441 filename.push_back( path + Form( "/share/%s/tof_endcap/tofec_bg_d10.root", dir ) );
1442 filename.push_back( path + Form( "/share/%s/tof_endcap/tofec_bg_d11.root", dir ) );
1443 filename.push_back( path + Form( "/share/%s/tof_endcap/tofec_bg_m10.root", dir ) );
1444 filename.push_back( path + Form( "/share/%s/tof_endcap/tofec_bg_m11.root", dir ) );
1445 for ( unsigned int i = 0; i < filename.size(); i++ )
1446 {
1447 f_tofec_bg[idx][i] = new TFile( filename[i].c_str(), "READ" );
1448 const char* name;
1449 if ( i == 0 ) name = "d10";
1450 else if ( i == 1 ) name = "d11";
1451 else if ( i == 2 ) name = "m10";
1452 else if ( i == 3 ) name = "m11";
1453 else
1454 {
1455 cout << "Boundary Error! " << endl;
1456 exit( 1 );
1457 }
1458 for ( unsigned int j = 0; j < 2; j++ )
1459 {
1460 h_tofec_p_bg_offset[idx][i][j] =
1461 (TH1D*)f_tofec_bg[idx][i]->Get( Form( "h_tofec_p_bg_offset_%s_%d", name, j ) );
1462 h_tofec_m_bg_offset[idx][i][j] =
1463 (TH1D*)f_tofec_bg[idx][i]->Get( Form( "h_tofec_m_bg_offset_%s_%d", name, j ) );
1464 h_tofec_p_bg_sigma[idx][i][j] =
1465 (TH1D*)f_tofec_bg[idx][i]->Get( Form( "h_tofec_p_bg_sigma_%s_%d", name, j ) );
1466 h_tofec_m_bg_sigma[idx][i][j] =
1467 (TH1D*)f_tofec_bg[idx][i]->Get( Form( "h_tofec_m_bg_sigma_%s_%d", name, j ) );
1468 }
1469 }
1470 // tof_endcap cost
1471 filename.clear();
1472 filename.push_back( path + Form( "/share/%s/tof_endcap/tofec_cost_d10.root", dir ) );
1473 filename.push_back( path + Form( "/share/%s/tof_endcap/tofec_cost_d11.root", dir ) );
1474 filename.push_back( path + Form( "/share/%s/tof_endcap/tofec_cost_m10.root", dir ) );
1475 filename.push_back( path + Form( "/share/%s/tof_endcap/tofec_cost_m11.root", dir ) );
1476 for ( unsigned int i = 0; i < filename.size(); i++ )
1477 {
1478 f_tofec_cost[idx][i] = new TFile( filename[i].c_str(), "READ" );
1479 const char* name;
1480 if ( i == 0 ) name = "d10";
1481 else if ( i == 1 ) name = "d11";
1482 else if ( i == 2 ) name = "m10";
1483 else if ( i == 3 ) name = "m11";
1484 else
1485 {
1486 cout << "Boundary Error! " << endl;
1487 exit( 1 );
1488 }
1489 for ( unsigned int j = 0; j < 2; j++ )
1490 {
1491 h_tofec_p_cost_offset[idx][i][j] =
1492 (TH1D*)f_tofec_cost[idx][i]->Get( Form( "h_tofec_p_cost_offset_%s_%d", name, j ) );
1493 h_tofec_m_cost_offset[idx][i][j] =
1494 (TH1D*)f_tofec_cost[idx][i]->Get( Form( "h_tofec_m_cost_offset_%s_%d", name, j ) );
1495 h_tofec_p_cost_sigma[idx][i][j] =
1496 (TH1D*)f_tofec_cost[idx][i]->Get( Form( "h_tofec_p_cost_sigma_%s_%d", name, j ) );
1497 h_tofec_m_cost_sigma[idx][i][j] =
1498 (TH1D*)f_tofec_cost[idx][i]->Get( Form( "h_tofec_m_cost_sigma_%s_%d", name, j ) );
1499 }
1500 }
1501 }
1502 for ( unsigned int idx = 0; idx < 3; idx++ )
1503 {
1504 const char* dir;
1505 if ( idx == 0 ) dir = "electron/emc";
1506 else if ( idx == 1 ) dir = "kpi/emc_pion";
1507 else if ( idx == 2 ) dir = "kpi/emc_kaon";
1508 else
1509 {
1510 cout << "Boundary Error! " << endl;
1511 exit( 1 );
1512 }
1513 // emc
1514 filename.clear();
1515 filename.push_back( path + Form( "/share/%s/emc_d10.root", dir ) );
1516 filename.push_back( path + Form( "/share/%s/emc_d11.root", dir ) );
1517 filename.push_back( path + Form( "/share/%s/emc_m10.root", dir ) );
1518 filename.push_back( path + Form( "/share/%s/emc_m11.root", dir ) );
1519 for ( unsigned int i = 0; i < filename.size(); i++ )
1520 {
1521 f_emc[idx][i] = new TFile( filename[i].c_str(), "READ" );
1522 const char* name;
1523 if ( i == 0 ) name = "d10";
1524 else if ( i == 1 ) name = "d11";
1525 else if ( i == 2 ) name = "m10";
1526 else if ( i == 3 ) name = "m11";
1527 else
1528 {
1529 cout << "Boundary Error! " << endl;
1530 exit( 1 );
1531 }
1532 for ( unsigned int j = 0; j < 15; j++ )
1533 {
1534 for ( unsigned int k = 0; k < 25; k++ )
1535 {
1536 h_emc_ep[idx][i][j][k] =
1537 (TH1D*)f_emc[idx][i]->Get( Form( "h_ep_%s_%d_%d", name, j, k ) );
1538 h_emc_e35[idx][i][j][k] =
1539 (TH1D*)f_emc[idx][i]->Get( Form( "h_e35_%s_%d_%d", name, j, k ) );
1540 }
1541 }
1542 }
1543 }
1544 cout << "Successfully Return from Loading Initializations by package SimplePIDSvc ... "
1545 << endl;
1546}
1547
1548int SimplePIDSvc::findBin( double* a, int length, double value ) {
1549 for ( int i = 0; i < length; i++ )
1550 {
1551 if ( value > a[i] && value <= a[i + 1] ) { return i; }
1552 }
1553 if ( value < a[0] ) { return 0; }
1554 else { return length; }
1555}
1556
1557double SimplePIDSvc::getChi2( int i ) {
1558 return pow( m_dedx_chi[i], 2 ) + pow( m_tof_chi[i], 2 );
1559}
1560
1561void SimplePIDSvc::loadEMCInfo( EvtRecTrack* track ) {
1562 // Initialization
1563 for ( unsigned int i = 0; i < 5; i++ )
1564 {
1565 m_emc_eop[i] = -99.;
1566 m_emc_likelihood[i] = -99.;
1567 }
1568 m_emc_e = -99.;
1569 m_emc_e13 = -99.;
1570 m_emc_e35 = -99.;
1571 m_emc_sec = -99.;
1572 m_emc_lat = -99.;
1573 m_lh_electron = -99.;
1574
1575 if ( !track->isEmcShowerValid() ) return;
1576
1577 RecEmcShower* emc_trk = track->emcShower();
1578 m_emc_e = emc_trk->energy();
1579 for ( unsigned int i = 0; i < 5; i++ ) { m_emc_eop[i] = m_emc_e / fabs( m_p[i] ); }
1580 double eseed = emc_trk->eSeed();
1581 double e3 = emc_trk->e3x3();
1582 double e5 = emc_trk->e5x5();
1583 m_emc_sec = emc_trk->secondMoment() / 1000.;
1584 m_emc_lat = emc_trk->latMoment();
1585 if ( e3 != 0 ) { m_emc_e13 = eseed / e3; }
1586 if ( e5 != 0 ) { m_emc_e35 = e3 / e5; }
1587}
1588
1589bool SimplePIDSvc::calEMCLikelihood() {
1590 if ( m_emc_eop[0] < 0 ) return false;
1591
1592 int idx = getRunIdx( m_run );
1593 const Int_t BIN_P = 15;
1594 const Int_t BIN_COST = 25;
1595 const Int_t BIN_PID = 3;
1596 const double EPS = 1e-4;
1597 // electron
1598 double P[BIN_PID][BIN_P + 1] = {
1599 { 0.20, 0.47, 0.56, 0.65, 0.72, 0.79, 0.86, 0.92, 0.98, 1.03, 1.08, 1.13, 1.17, 1.22,
1600 1.26, 1.30 },
1601 { 0.20, 0.26, 0.31, 0.35, 0.39, 0.42, 0.46, 0.49, 0.53, 0.57, 0.62, 0.67, 0.73, 0.80,
1602 0.88, 1.05 },
1603 { 0.20, 0.33, 0.39, 0.43, 0.48, 0.52, 0.56, 0.61, 0.67, 0.73, 0.76, 0.81, 0.85, 0.90,
1604 0.96, 1.05 },
1605 };
1606 double COST[BIN_PID][BIN_COST + 1] = {
1607 { -0.930, -0.910, -0.905, -0.897, -0.890, -0.881, -0.871, -0.858, -0.775,
1608 -0.732, -0.669, -0.561, -0.330, 0.199, 0.515, 0.645, 0.718, 0.766,
1609 0.804, 0.870, 0.882, 0.891, 0.898, 0.906, 0.913, 0.930 },
1610 { -0.930, -0.810, -0.728, -0.648, -0.574, -0.501, -0.431, -0.364, -0.295,
1611 -0.228, -0.161, -0.096, -0.031, 0.035, 0.100, 0.167, 0.234, 0.301,
1612 0.370, 0.439, 0.510, 0.580, 0.655, 0.733, 0.813, 0.930 },
1613 { -0.930, -0.804, -0.721, -0.643, -0.568, -0.497, -0.429, -0.362, -0.293,
1614 -0.228, -0.161, -0.096, -0.029, 0.035, 0.100, 0.166, 0.233, 0.298,
1615 0.365, 0.432, 0.500, 0.571, 0.644, 0.722, 0.805, 0.930 },
1616 };
1617
1618 double vp, vcost;
1619 int pid;
1620 int bin_p, bin_cost;
1621 for ( unsigned int i = 0; i < 4; i++ )
1622 {
1623 // only e, pi ,K
1624 if ( i == 0 ) pid = 0;
1625 else if ( i == 2 ) pid = 1;
1626 else if ( i == 3 ) pid = 2;
1627 else continue;
1628
1629 vp = max( P[pid][0] + EPS, min( P[pid][BIN_P] - EPS, m_p[i] ) );
1630 vcost = max( COST[pid][0] + EPS, min( COST[pid][BIN_COST] - EPS, m_cost[i] ) );
1631 bin_p = findBin( P[pid], BIN_P, vp );
1632 bin_cost = findBin( COST[pid], BIN_COST, vcost );
1633
1634 m_emc_likelihood[i] = h_emc_ep[pid][idx][bin_p][bin_cost]->Interpolate( m_emc_eop[i] ) *
1635 h_emc_e35[pid][idx][bin_p][bin_cost]->Interpolate( m_emc_e35 );
1636 }
1637 double a = m_prob[0] > 0 ? m_prob[0] : 0;
1638 double b = m_prob[2] > 0 ? m_prob[2] : 0;
1639 double c = m_prob[3] > 0 ? m_prob[3] : 0;
1640 double sum = a * m_emc_likelihood[0] + b * m_emc_likelihood[2] + c * m_emc_likelihood[3];
1641
1642 if ( sum > 0 && m_prob[0] > 0 )
1643 {
1644 m_lh_electron = m_prob[0] * m_emc_likelihood[0] / sum;
1645 return true;
1646 }
1647 else { return false; }
1648}
1649
1650bool SimplePIDSvc::ispion() {
1651 if ( m_prob[2] > 0.00 && m_prob[2] > m_prob[3] ) return true;
1652 else return false;
1653}
1654
1655bool SimplePIDSvc::iskaon() {
1656 if ( m_prob[3] > 0.00 && m_prob[3] > m_prob[2] ) return true;
1657 else return false;
1658}
1659
1660// bool SimplePIDSvc::iselectron_(bool emc)
1661//{
1662// if (m_prob[0] > 0 && m_prob[0] > m_prob[2] && m_prob[0] > m_prob[3])
1663// {
1664// if (!emc)
1665// return true;
1666// else if (fabs(m_cost[0]) < 0.7 && m_emc_eop[0] > 0 && m_emc_eop[0] <
1667//0.8) return false; else if (fabs(m_cost[0]) >= 0.7 && fabs(m_cost[0])<0.8 && m_emc_eop[0]
1668//> 0 && m_emc_eop[0] < -7.5*fabs(m_cost[0])+6.05) return false; else if (fabs(m_cost[0]) >
1669//0.85 && m_emc_eop[0] > 0 && m_emc_eop[0] < 0.6) return false; else return true;
1670// }
1671// else
1672// return false;
1673// }
1674
1675bool SimplePIDSvc::iselectron( bool emc ) {
1676 if ( !emc )
1677 {
1678 if ( m_prob[0] > 0 && m_prob[0] > m_prob[2] && m_prob[0] > m_prob[3] ) return true;
1679 else return false;
1680 }
1681 else
1682 {
1683 if ( calEMCLikelihood() )
1684 {
1685 if ( m_lh_electron > m_eid_ratio ) return true;
1686 else return false;
1687 }
1688 else
1689 {
1690 if ( m_prob[0] > 0 && m_prob[0] > m_prob[2] && m_prob[0] > m_prob[3] ) return true;
1691 else return false;
1692 }
1693 }
1694}
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
P
double P(RecMdcKalTrack *trk)
Definition CalibEventSelect.cxx:84
mass
double mass
Definition CosmicGenerator.cxx:128
min
#define min(a, b)
Definition Eepipi/src/ee2eepp/basesv5.1/f2c.h:157
max
#define max(a, b)
Definition Eepipi/src/ee2eepp/basesv5.1/f2c.h:158
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition KKsem.h:33
dt
TGraph2DErrors * dt
Definition McCor.cxx:41
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
EPS
const double EPS
Definition TRunge.cxx:42
EvtRecTrack::tofTrack
SmartRefVector< RecTofTrack > tofTrack()
Definition Event/EvtRecEvent/include/EvtRecEvent/EvtRecTrack.h:57
SimplePIDSvc::SimplePIDSvc
SimplePIDSvc(const std::string &name, ISvcLocator *svcLoc)
Definition SimplePIDSvc.cxx:19
SimplePIDSvc::~SimplePIDSvc
virtual ~SimplePIDSvc()
Definition SimplePIDSvc.cxx:28
SimplePIDSvc::getChi2
double getChi2(int i)
Definition SimplePIDSvc.cxx:1557
SimplePIDSvc::preparePID
void preparePID(EvtRecTrack *track)
Definition SimplePIDSvc.cxx:87
SimplePIDSvc::initialize
virtual StatusCode initialize()
Definition SimplePIDSvc.cxx:30
SimplePIDSvc::iselectron
bool iselectron(bool emc=true)
Definition SimplePIDSvc.cxx:1675
SimplePIDSvc::finalize
virtual StatusCode finalize()
Definition SimplePIDSvc.cxx:44
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition TofHitStatus.cxx:115
EvtCyclic3::c_str
char * c_str(Index i)
Definition EvtCyclic3.cc:240
t
int t()
Definition t.c:1