161 {
162
163
164
165 MsgStream log(
msgSvc(), name() );
166 log << MSG::INFO << "in execute()" << endmsg;
167
168 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
169 int runNo = eventHeader->runNumber();
170 int event = eventHeader->eventNumber();
171 log << MSG::DEBUG <<
"run, evtnum = " <<
runNo <<
" , " <<
event << endmsg;
173 m_nrec = event;
174
175 setFilterPassed( false );
176
177 if ( ( event % 100000 ) == 0 ) { cout << "event: " << event << endl; }
178 Ncut0++;
179
181
182 log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
183 << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
184
186
187 SmartDataPtr<Event::McParticleCol> mcParticleCol( eventSvc(), "/Event/MC/McParticleCol" );
188 double temp_p_pp, temp_p_pm;
189
190 int t_j = 0;
191
192
193
194
195
196 Vint iGood, ipip, ipim, ikaonp, ikaonm, iprotonp, iprotonm;
197 iGood.clear();
198 ipip.clear();
199 ipim.clear();
200 ikaonp.clear();
201 ikaonm.clear();
202 iprotonp.clear();
203 iprotonm.clear();
205 ppip.clear();
206 ppim.clear();
207
208 int nCharge = 0;
209
210 Hep3Vector xorigin( 0, 0, 0 );
211
212
213 IVertexDbSvc* vtxsvc;
214 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc ).ignore();
216 {
219
220
221 xorigin.setX( dbv[0] );
222 xorigin.setY( dbv[1] );
223 xorigin.setZ( dbv[2] );
224 }
225
226 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
227 {
229 if ( !( *itTrk )->isMdcTrackValid() ) continue;
230 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
231 double pch = mdcTrk->
p();
232 double x0 = mdcTrk->
x();
233 double y0 = mdcTrk->
y();
234 double z0 = mdcTrk->
z();
235 double phi0 = mdcTrk->
helix( 1 );
236 double xv = xorigin.x();
237 double yv = xorigin.y();
238 double Rxy = ( x0 - xv ) *
cos( phi0 ) + ( y0 - yv ) *
sin( phi0 );
239
240
241
242
243
244 HepVector a = mdcTrk->
helix();
245 HepSymMatrix Ea = mdcTrk->
err();
247 HepPoint3D IP( xorigin[0], xorigin[1], xorigin[2] );
248 VFHelix helixip( point0, a, Ea );
249 helixip.pivot( IP );
250 HepVector vecipa = helixip.a();
251 double Rvxy0 = fabs( vecipa[0] );
252 double Rvz0 = vecipa[3];
253 double Rvphi0 = vecipa[1];
254
255 if ( fabs( Rvz0 ) >= m_vz0cut ) continue;
256 if ( fabs( Rvxy0 ) >= m_vr0cut ) continue;
257 if (
cos( mdcTrk->
theta() ) > 0.93 )
continue;
258 if ( pch > 5 ) continue;
259
260 iGood.push_back( i );
261 nCharge += mdcTrk->
charge();
262 if ( t_j < 4 )
263 {
264 if ( ( *itTrk )->isMdcDedxValid() )
265 {
266 RecMdcDedx* dedxTrk = ( *itTrk )->mdcDedx();
267
268 m_dedxchi_e[t_j] = dedxTrk->
chiE();
269 m_dedxchi_mu[t_j] = dedxTrk->
chiMu();
270 m_dedxchi_pi[t_j] = dedxTrk->
chiPi();
271 m_dedxchi_kaon[t_j] = dedxTrk->
chiK();
272 m_dedxchi_proton[t_j] = dedxTrk->
chiP();
273
275 }
276 if ( ( *itTrk )->isTofTrackValid() )
277 {
278 SmartRefVector<RecTofTrack> tofTrkCol = ( *itTrk )->tofTrack();
279 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
280 for ( ; iter_tof != tofTrkCol.end(); iter_tof++ )
281 {
282 TofHitStatus* status = new TofHitStatus;
283 status->
setStatus( ( *iter_tof )->status() );
285 {
286 double time = ( *iter_tof )->tof();
287 double sigma = 1.1 * ( *iter_tof )->sigma( 0 ) / 1000;
288 double expE_tof = ( *iter_tof )->texpElectron();
289 double expMu_tof = ( *iter_tof )->texpMuon();
290 double expPi_tof = ( *iter_tof )->texpPion();
291 double expK_tof = ( *iter_tof )->texpKaon();
292 double expP_tof = ( *iter_tof )->texpProton();
293
294 if ( sigma != 0 )
295 {
296
297 m_tofchi_e[t_j] = (
time - expE_tof );
298 m_tofchi_mu[t_j] = (
time - expMu_tof );
299 m_tofchi_pi[t_j] = (
time - expPi_tof );
300 m_tofchi_kaon[t_j] = (
time - expK_tof );
301 m_tofchi_proton[t_j] = (
time - expP_tof );
302 }
303 }
304 delete status;
305 }
306 }
307
308 t_j++;
309 }
310 }
311
312
313
314
315 int nGood = iGood.size();
316 m_nGood = nGood;
317 m_nCharge = nCharge;
318 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endmsg;
319 if ( ( nGood != 4 ) || ( nCharge != 0 ) ) { return StatusCode::SUCCESS; }
320 Ncut1++;
321
322 for ( int i = 0; i < nGood; i++ )
323 {
324 m_eop[i] = 5.;
326 if ( !( *itTrk )->isMdcTrackValid() ) continue;
327 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
328 double p = mdcTrk->
p();
329 m_eop[i] = 6.;
330 if ( !( *itTrk )->isEmcShowerValid() ) continue;
331 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
332 double eraw = emcTrk->
energy();
333 m_eop[i] = eraw / p;
334 }
335
337 iGam.clear();
338 for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
339 {
340 if ( i >= evtRecTrkCol->size() ) break;
342 if ( !( *itTrk )->isEmcShowerValid() ) continue;
343 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
344 Hep3Vector emcpos( emcTrk->
x(), emcTrk->
y(), emcTrk->
z() );
345
346 double dthe = 200.;
347 double dphi = 200.;
348 double dang = 200.;
349 for ( int j = 0; j < evtRecEvent->totalCharged(); j++ )
350 {
352 if ( !( *jtTrk )->isExtTrackValid() ) continue;
353 RecExtTrack* extTrk = ( *jtTrk )->extTrack();
356
357 double angd = extpos.angle( emcpos );
358 double thed = extpos.theta() - emcpos.theta();
359 double phid = extpos.deltaPhi( emcpos );
360 thed = fmod( thed + CLHEP::twopi + CLHEP::twopi +
pi, CLHEP::twopi ) - CLHEP::pi;
361 phid = fmod( phid + CLHEP::twopi + CLHEP::twopi +
pi, CLHEP::twopi ) - CLHEP::pi;
362 if ( angd < dang )
363 {
364 dang = angd;
365 dthe = thed;
366 dphi = phid;
367 }
368 }
369 if ( dang >= 200 ) continue;
370 double eraw = emcTrk->
energy();
371 dthe = dthe * 180 / ( CLHEP::pi );
372 dphi = dphi * 180 / ( CLHEP::pi );
373 dang = dang * 180 / ( CLHEP::pi );
374 if ( eraw < m_energyThreshold ) continue;
375
376 if ( fabs( dang ) < m_gammaAngleCut ) continue;
377 double dtime = emcTrk->
time();
378
379 if ( dtime < 0 ) continue;
380 if ( dtime > 14 ) continue;
381
382
383
384
385 iGam.push_back( ( *itTrk )->trackId() );
386 }
387
388
389
390
391 int nGam = iGam.size();
392 m_nGam = nGam;
393
394 log << MSG::DEBUG << "num Good Photon " << nGam << " , " << evtRecEvent->totalNeutral()
395 << endmsg;
396
397
398
399 Ncut2++;
400
401
402
403
404
405
406
407
409 pCh.clear();
411 iPid.clear();
412 iCh.clear();
413 int npi = 0, npip = 0, npim = 0;
414 int nkaon = 0, nkaonp = 0, nkaonm = 0;
415 int nproton = 0, nprotonp = 0, nprotonm = 0;
417
418 for ( int i = 0; i < nGood; i++ )
419 {
421
424
425
429
431
432
433
435 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
436
437
438 iCh.push_back( mdcTrk->
charge() );
440 {
441 iPid.push_back( 0 );
442 HepLorentzVector ptrk;
443 ptrk.setPx( mdcTrk->
px() );
444 ptrk.setPy( mdcTrk->
py() );
445 ptrk.setPz( mdcTrk->
pz() );
446 double p3 = ptrk.vect().mag();
447 ptrk.setE( sqrt( p3 * p3 +
mpi *
mpi ) );
448 pCh.push_back( ptrk );
449 }
451
452
453
454
455
456
457 RecMdcKalTrack* mdcKalTrk =
458 ( *itTrk )->mdcKalTrack();
459
460
461
462
463
464 HepLorentzVector ptrk;
465 ptrk.setPx( mdcTrk->
px() );
466 ptrk.setPy( mdcTrk->
py() );
467 ptrk.setPz( mdcTrk->
pz() );
468
470 {
471 npi = npi + 1;
472 iPid.push_back( 2 );
473 if ( ( *itTrk )->isMdcKalTrackValid() )
474 {
476 ptrk.setPx( mdcKalTrk->
px() );
477 ptrk.setPy( mdcKalTrk->
py() );
478 ptrk.setPz( mdcKalTrk->
pz() );
479 }
480 double p3 = ptrk.vect().mag();
481 ptrk.setE( sqrt( p3 * p3 +
mpi *
mpi ) );
482 pCh.push_back( ptrk );
483 if ( mdcTrk->
charge() > 0 ) npip++;
484 if ( mdcTrk->
charge() < 0 ) npim++;
485 }
488 {
489 nkaon = nkaon + 1;
490 iPid.push_back( 3 );
491 if ( ( *itTrk )->isMdcKalTrackValid() )
492 {
494 ptrk.setPx( mdcKalTrk->
px() );
495 ptrk.setPy( mdcKalTrk->
py() );
496 ptrk.setPz( mdcKalTrk->
pz() );
497 }
498 double p3 = ptrk.vect().mag();
500 pCh.push_back( ptrk );
501 }
504 {
505 nproton = nproton + 1;
506 iPid.push_back( 4 );
507 if ( ( *itTrk )->isMdcKalTrackValid() )
508 {
510 ptrk.setPx( mdcKalTrk->
px() );
511 ptrk.setPy( mdcKalTrk->
py() );
512 ptrk.setPz( mdcKalTrk->
pz() );
513 }
514 double p3 = ptrk.vect().mag();
516 pCh.push_back( ptrk );
517 if ( mdcTrk->
charge() > 0 ) nprotonp++;
518 if ( mdcTrk->
charge() < 0 ) nprotonm++;
519 }
520 else
521 {
522 iPid.push_back( 0 );
523 HepLorentzVector ptrk;
524 ptrk.setPx( mdcTrk->
px() );
525 ptrk.setPy( mdcTrk->
py() );
526 ptrk.setPz( mdcTrk->
pz() );
527 double p3 = ptrk.vect().mag();
528 ptrk.setE( sqrt( p3 * p3 +
mpi *
mpi ) );
529 pCh.push_back( ptrk );
530
531 }
532 }
533 m_npi = npi;
534 m_nkaon = nkaon;
535 m_nproton = nproton;
536
537
538
539 m_istat_pmiss = 0;
540 m_istat_pbmiss = 0;
541 m_istat_pipmiss = 0;
542 m_istat_pimmiss = 0;
543 for ( int i = 0; i < 4; i++ )
544 {
545 m_ptrk_pmiss[i] = 5;
546 m_ptrk_pbmiss[i] = 5;
547 m_ptrk_pipmiss[i] = 5;
548 m_ptrk_pimmiss[i] = 5;
549 }
550 for ( int i = 0; i < 3; i++ )
551 {
552 m_id_pmiss[i] = 0;
553 m_id_pbmiss[i] = 0;
554 m_id_pipmiss[i] = 0;
555 m_id_pimmiss[i] = 0;
556 }
557
558
559
560 HepLorentzVector
ecms( 0.03406837, 0, 0, 3.0971873 );
561 HepLorentzVector pmiss;
562 m_mpmiss = 5.;
563 m_mpbmiss = 5.;
564 m_mpipmiss = 5.;
565 m_mpimmiss = 5.;
566 m_ppmiss = 5.;
567 m_ppbmiss = 5.;
568 m_ppipmiss = 5.;
569 m_ppimmiss = 5.;
570
571 int istat_nhit;
572
573 if ( ( npip == 1 ) && ( npim == 1 ) && ( nprotonm == 1 ) )
574 {
575 pmiss.setPx( 0 );
576 pmiss.setPy( 0 );
577 pmiss.setPz( 0 );
578 pmiss.setE( 0 );
579 istat_nhit = 0;
580
581 int j = 0;
582 for ( int i = 0; i < nGood; i++ )
583 {
585
586 RecMdcKalTrack* mdcKalTrk =
587 ( *itTrk )->mdcKalTrack();
588
589
590
591 double p;
592 double eraw;
593 if ( ( *itTrk )->isEmcShowerValid() )
594 {
595 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
597 }
598
599
600 if ( ( iPid[i] * iCh[i] ) == 2 )
601 {
602 m_index_pmiss[j] = i;
603 pmiss = pmiss + pCh[i];
606 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
607 j++;
608 }
609 else if ( ( iPid[i] * iCh[i] ) == -2 )
610 {
611 m_index_pmiss[j] = i;
612 pmiss = pmiss + pCh[i];
615 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
616
617 j++;
618 }
619 else if ( ( iPid[i] * iCh[i] ) == -4 )
620 {
621 m_index_pmiss[j] = i;
622 pmiss = pmiss + pCh[i];
625 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
626
627 j++;
628 }
629 else
630 {
631 if ( nGood == 4 )
632 {
633 m_index_pmiss[3] = i;
636 HepLorentzVector ptrk;
637 ptrk.setPx( mdcKalTrk->
px() );
638 ptrk.setPy( mdcKalTrk->
py() );
639 ptrk.setPz( mdcKalTrk->
pz() );
640 double p3 = ptrk.vect().mag();
642 m_ptrk_pmiss[0] = ptrk.px();
643 m_ptrk_pmiss[1] = ptrk.py();
644 m_ptrk_pmiss[2] = ptrk.pz();
645 m_ptrk_pmiss[3] = ptrk.e();
646 }
647 }
648 }
649
650 for ( int i = 0; i < nGood; i++ )
651 {
652 if ( ( iPid[i] * iCh[i] ) == 2 ) continue;
653 if ( ( iPid[i] * iCh[i] ) == -2 ) continue;
654 if ( ( iPid[i] * iCh[i] ) == -4 ) continue;
655 if ( iCh[i] < 0 ) continue;
656
658 for ( int ii = 0; ii < 3; ii++ )
659 {
665 else if ( ii == 1 )
667 else if ( ii == 2 )
668 {
670 }
671
673
678 m_id_pmiss[ii]++;
679 }
680 }
681
682 pmiss =
ecms - pmiss;
683 m_mpmiss = pmiss.m();
684 m_ppmiss = pmiss.rho();
685
686 if ( fabs( m_mpmiss -
mproton ) < 0.02 && istat_nhit == 0 )
687 {
688 m_istat_pmiss = 1;
689 ( *( evtRecTrkCol->begin() + iGood[m_index_pmiss[3]] ) )->setPartId( 3 );
690 ( *( evtRecTrkCol->begin() + iGood[m_index_pmiss[3]] ) )->setQuality( 0 );
691 }
692
693 Ncut3++;
694 }
695 if ( ( npip == 1 ) && ( npim == 1 ) && ( nprotonp == 1 ) )
696 {
697 pmiss.setPx( 0 );
698 pmiss.setPy( 0 );
699 pmiss.setPz( 0 );
700 pmiss.setE( 0 );
701 istat_nhit = 0;
702
703 int j = 0;
704 for ( int i = 0; i < nGood; i++ )
705 {
707
708 RecMdcKalTrack* mdcKalTrk =
709 ( *itTrk )->mdcKalTrack();
710
711
712
713 double p;
714 double eraw;
715 if ( ( *itTrk )->isEmcShowerValid() )
716 {
717 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
719 }
720
721
722
723 if ( ( iPid[i] * iCh[i] ) == 2 )
724 {
725 m_index_pbmiss[j] = i;
726 pmiss = pmiss + pCh[i];
729 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
730
731 j++;
732 }
733 else if ( ( iPid[i] * iCh[i] ) == -2 )
734 {
735 m_index_pbmiss[j] = i;
736
737 pmiss = pmiss + pCh[i];
740 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
741
742 j++;
743 }
744 else if ( ( iPid[i] * iCh[i] ) == 4 )
745 {
746 m_index_pbmiss[j] = i;
747 pmiss = pmiss + pCh[i];
750 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
751
752 j++;
753 }
754 else
755 {
756 if ( nGood == 4 )
757 {
758 m_index_pbmiss[3] = i;
761 HepLorentzVector ptrk;
762 ptrk.setPx( mdcKalTrk->
px() );
763 ptrk.setPy( mdcKalTrk->
py() );
764 ptrk.setPz( mdcKalTrk->
pz() );
765 double p3 = ptrk.vect().mag();
767 m_ptrk_pbmiss[0] = ptrk.px();
768 m_ptrk_pbmiss[1] = ptrk.py();
769 m_ptrk_pbmiss[2] = ptrk.pz();
770 m_ptrk_pbmiss[3] = ptrk.e();
771 }
772 }
773 }
774
775 for ( int i = 0; i < nGood; i++ )
776 {
777 if ( ( iPid[i] * iCh[i] ) == 2 ) continue;
778 if ( ( iPid[i] * iCh[i] ) == -2 ) continue;
779 if ( ( iPid[i] * iCh[i] ) == 4 ) continue;
780 if ( iCh[i] > 0 ) continue;
781
783 for ( int ii = 0; ii < 3; ii++ )
784 {
790 else if ( ii == 1 )
792 else if ( ii == 2 )
793 {
795 }
796
798
803 m_id_pbmiss[ii]++;
804 }
805 }
806 pmiss =
ecms - pmiss;
807 m_mpbmiss = pmiss.m();
808 m_ppbmiss = pmiss.rho();
809 if ( fabs( m_mpbmiss -
mproton ) < 0.02 && istat_nhit == 0 )
810 {
811 m_istat_pbmiss = 1;
812 ( *( evtRecTrkCol->begin() + iGood[m_index_pbmiss[3]] ) )->setPartId( 3 );
813 ( *( evtRecTrkCol->begin() + iGood[m_index_pbmiss[3]] ) )->setQuality( 0 );
814 }
815 Ncut4++;
816 }
817 int initial_pip, initial_pim;
818
819
820
821
822 if ( ( npim == 1 ) && ( nprotonp == 1 ) && ( nprotonm == 1 ) )
823 {
824 pmiss.setPx( 0 );
825 pmiss.setPy( 0 );
826 pmiss.setPz( 0 );
827 pmiss.setE( 0 );
828 istat_nhit = 0;
829
830 int j = 0;
831 for ( int i = 0; i < nGood; i++ )
832 {
834
835 RecMdcKalTrack* mdcKalTrk =
836 ( *itTrk )->mdcKalTrack();
837
838
839
840 double p;
841 double eraw;
842 if ( ( *itTrk )->isEmcShowerValid() )
843 {
844 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
846 }
847
848
849
850 if ( ( iPid[i] * iCh[i] ) == 4 )
851 {
852 m_index_pipmiss[j] = i;
853 pmiss = pmiss + pCh[i];
856 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
857
858 j++;
859 }
860 else if ( ( iPid[i] * iCh[i] ) == -4 )
861 {
862 m_index_pipmiss[j] = i;
863
864 pmiss = pmiss + pCh[i];
867 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
868
869 j++;
870 }
871 else if ( ( iPid[i] * iCh[i] ) == -2 )
872 {
873 m_index_pipmiss[j] = i;
874
875 pmiss = pmiss + pCh[i];
878 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
879
880 j++;
881 }
882 else
883 {
884 if ( nGood == 4 )
885 {
886 m_index_pipmiss[3] = i;
887
890 HepLorentzVector ptrk;
891 ptrk.setPx( mdcKalTrk->
px() );
892 ptrk.setPy( mdcKalTrk->
py() );
893 ptrk.setPz( mdcKalTrk->
pz() );
894 double p3 = ptrk.vect().mag();
895 ptrk.setE( sqrt( p3 * p3 +
mpi *
mpi ) );
896 m_ptrk_pipmiss[0] = ptrk.px();
897 m_ptrk_pipmiss[1] = ptrk.py();
898 m_ptrk_pipmiss[2] = ptrk.pz();
899 m_ptrk_pipmiss[3] = ptrk.e();
900 }
901 }
902 }
903
904 for ( int i = 0; i < nGood; i++ )
905 {
906 if ( ( iPid[i] * iCh[i] ) == 4 ) continue;
907 if ( ( iPid[i] * iCh[i] ) == -4 ) continue;
908 if ( ( iPid[i] * iCh[i] ) == -2 ) continue;
909 if ( iCh[i] < 0 ) continue;
910
912 for ( int ii = 0; ii < 3; ii++ )
913 {
919 else if ( ii == 1 )
921 else if ( ii == 2 )
922 {
924 }
925
927
932 m_id_pipmiss[ii]++;
933 }
934 }
935 pmiss =
ecms - pmiss;
936 m_mpipmiss = pmiss.m();
937 m_ppipmiss = pmiss.rho();
938 if ( fabs( m_mpipmiss -
mpi ) < 0.05 && istat_nhit == 0 )
939 {
940 m_istat_pipmiss = 1;
941 ( *( evtRecTrkCol->begin() + iGood[m_index_pipmiss[3]] ) )->setPartId( 2 );
942 ( *( evtRecTrkCol->begin() + iGood[m_index_pipmiss[3]] ) )->setQuality( 0 );
943 }
944 Ncut5++;
945 }
946 if ( ( npip == 1 ) && ( nprotonp == 1 ) && ( nprotonm == 1 ) )
947 {
948 pmiss.setPx( 0 );
949 pmiss.setPy( 0 );
950 pmiss.setPz( 0 );
951 pmiss.setE( 0 );
952 istat_nhit = 0;
953
954 int j = 0;
955 for ( int i = 0; i < nGood; i++ )
956 {
958
959 RecMdcKalTrack* mdcKalTrk =
960 ( *itTrk )->mdcKalTrack();
961
962
963
964 double p;
965 double eraw;
966 if ( ( *itTrk )->isEmcShowerValid() )
967 {
968 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
970 }
971
972
973
974
975 if ( ( iPid[i] * iCh[i] ) == 4 )
976 {
977 m_index_pimmiss[j] = i;
978 pmiss = pmiss + pCh[i];
981 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
982
983 j++;
984 }
985 else if ( ( iPid[i] * iCh[i] ) == -4 )
986 {
987 m_index_pimmiss[j] = i;
988
989 pmiss = pmiss + pCh[i];
992 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
993
994 j++;
995 }
996 else if ( ( iPid[i] * iCh[i] ) == 2 )
997 {
998 m_index_pimmiss[j] = i;
999 pmiss = pmiss + pCh[i];
1002 if ( m_dedxngoodhit[i] < 20 ) istat_nhit = 1;
1003
1004 j++;
1005 }
1006 else
1007 {
1008 if ( nGood == 4 )
1009 {
1010 m_index_pimmiss[3] = i;
1013 HepLorentzVector ptrk;
1014 ptrk.setPx( mdcKalTrk->
px() );
1015 ptrk.setPy( mdcKalTrk->
py() );
1016 ptrk.setPz( mdcKalTrk->
pz() );
1017 double p3 = ptrk.vect().mag();
1018 ptrk.setE( sqrt( p3 * p3 +
mpi *
mpi ) );
1019 m_ptrk_pimmiss[0] = ptrk.px();
1020 m_ptrk_pimmiss[1] = ptrk.py();
1021 m_ptrk_pimmiss[2] = ptrk.pz();
1022 m_ptrk_pimmiss[3] = ptrk.e();
1023 }
1024 }
1025 }
1026
1027 for ( int i = 0; i < nGood; i++ )
1028 {
1029 if ( ( iPid[i] * iCh[i] ) == 4 ) continue;
1030 if ( ( iPid[i] * iCh[i] ) == -4 ) continue;
1031 if ( ( iPid[i] * iCh[i] ) == 2 ) continue;
1032 if ( iCh[i] > 0 ) continue;
1033
1035 for ( int ii = 0; ii < 3; ii++ )
1036 {
1042 else if ( ii == 1 )
1044 else if ( ii == 2 )
1045 {
1047 }
1048
1050
1055 m_id_pimmiss[ii]++;
1056 }
1057 }
1058 pmiss =
ecms - pmiss;
1059 m_mpimmiss = pmiss.m();
1060 m_ppimmiss = pmiss.rho();
1061 if ( fabs( m_mpimmiss -
mpi ) < 0.05 && istat_nhit == 0 )
1062 {
1063 m_istat_pimmiss = 1;
1064 ( *( evtRecTrkCol->begin() + iGood[m_index_pimmiss[3]] ) )->setPartId( 2 );
1065 ( *( evtRecTrkCol->begin() + iGood[m_index_pimmiss[3]] ) )->setQuality( 0 );
1066 }
1067 Ncut6++;
1068 }
1069
1070 if ( m_istat_pmiss != 1 && m_istat_pbmiss != 1 && m_istat_pipmiss != 1 &&
1071 m_istat_pimmiss != 1 )
1072 { return StatusCode::SUCCESS; }
1073
1074
1075 if ( m_saventuple ) m_tuple0->write().ignore();
1076
1077 setFilterPassed( true );
1078
1079 return StatusCode::SUCCESS;
1080}
HepGeom::Point3D< double > HepPoint3D
EvtRecTrackCol::iterator EvtRecTrackIterator
std::vector< HepLorentzVector > Vp4
double sin(const BesAngle a)
double cos(const BesAngle a)
const Hep3Vector emcPosition() const
const int emcVolumeNumber() const
static void setPidType(PidType pidType)
const double theta() const
const HepSymMatrix err() const
const HepVector helix() const
......
virtual bool isVertexValid()=0
virtual double * SigmaPrimaryVertex()=0
virtual double * PrimaryVertex()=0
int methodProbability() const
void setChiMinCut(const double chi=4)
void setRecTrack(EvtRecTrack *trk)
void setMethod(const int method)
void identify(const int pidcase)
void usePidSys(const int pidsys)
static ParticleID * instance()
bool IsPidInfoValid() const
double probProton() const
void setStatus(unsigned int status)
_EXTERN_ std::string EvtRecEvent
_EXTERN_ std::string EvtRecTrackCol