176 {
177
178
179
180 MsgStream log(
msgSvc(), name());
181 log << MSG::INFO << "in execute()" << endreq;
182
183 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
184 int runNo=eventHeader->runNumber();
185 int event=eventHeader->eventNumber();
186 log << MSG::DEBUG <<"run, evtnum = "
188 << event <<endreq;
190 m_nrec= event;
191
192 setFilterPassed(false);
193
194 if ((event%100000)==0) {cout <<"event: "<< event << endl ; }
195 Ncut0++;
196
198
199 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
200 << evtRecEvent->totalCharged() << " , "
201 << evtRecEvent->totalNeutral() << " , "
202 << evtRecEvent->totalTracks() <<endreq;
203
205
206 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
207 double temp_p_pp,temp_p_pm;
208
209 int t_j=0;
210
211
212
213
214
215 Vint iGood, ipip, ipim , ikaonp, ikaonm, iprotonp,iprotonm;
216 iGood.clear();
217 ipip.clear();
218 ipim.clear();
219 ikaonp.clear();
220 ikaonm.clear();
221 iprotonp.clear();
222 iprotonm.clear();
224 ppip.clear();
225 ppim.clear();
226
227 int nCharge = 0;
228
229 Hep3Vector xorigin(0,0,0);
230
231
233 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
237
238
239 xorigin.setX(dbv[0]);
240 xorigin.setY(dbv[1]);
241 xorigin.setZ(dbv[2]);
242 }
243
244 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
246 if(!(*itTrk)->isMdcTrackValid()) continue;
248 double pch=mdcTrk->
p();
249 double x0=mdcTrk->
x();
250 double y0=mdcTrk->
y();
251 double z0=mdcTrk->
z();
252 double phi0=mdcTrk->
helix(1);
253 double xv=xorigin.x();
254 double yv=xorigin.y();
255 double Rxy=(x0-xv)*
cos(phi0)+(y0-yv)*
sin(phi0);
256
257
258
259
260
261
262 HepVector a = mdcTrk->
helix();
263 HepSymMatrix Ea = mdcTrk->
err();
265 HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]);
267 helixip.pivot(IP);
268 HepVector vecipa = helixip.a();
269 double Rvxy0=fabs(vecipa[0]);
270 double Rvz0=vecipa[3];
271 double Rvphi0=vecipa[1];
272
273 if(fabs(Rvz0) >= m_vz0cut) continue;
274 if(fabs(Rvxy0) >= m_vr0cut) continue;
275 if(
cos(mdcTrk->
theta())>0.93)
continue;
276 if(pch>5) continue;
277
278 iGood.push_back(i);
279 nCharge += mdcTrk->
charge();
280 if(t_j<4)
281 {
282 if((*itTrk)->isMdcDedxValid())
283 {
285
286 m_dedxchi_e[t_j] = dedxTrk->
chiE();
287 m_dedxchi_mu[t_j] = dedxTrk->
chiMu();
288 m_dedxchi_pi[t_j] = dedxTrk->
chiPi();
289 m_dedxchi_kaon[t_j] = dedxTrk->
chiK();
290 m_dedxchi_proton[t_j] = dedxTrk->
chiP();
291
293 }
294 if((*itTrk)->isTofTrackValid())
295 {
296 SmartRefVector<RecTofTrack> tofTrkCol=(*itTrk)->tofTrack();
297 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
298 for(;iter_tof!=tofTrkCol.end();iter_tof++){
300 status->
setStatus( (*iter_tof)->status() );
302 double time=(*iter_tof)->tof();
303 double sigma=1.1*(*iter_tof)->sigma(0)/1000;
304 double expE_tof=(*iter_tof)->texpElectron();
305 double expMu_tof=(*iter_tof)->texpMuon();
306 double expPi_tof=(*iter_tof)->texpPion();
307 double expK_tof=(*iter_tof)->texpKaon();
308 double expP_tof=(*iter_tof)->texpProton();
309
311
312 m_tofchi_e[t_j] = ( time- expE_tof );
313 m_tofchi_mu[t_j] = ( time- expMu_tof);
314 m_tofchi_pi[t_j] = ( time- expPi_tof );
315 m_tofchi_kaon[t_j] = ( time- expK_tof );
316 m_tofchi_proton[t_j] = ( time- expP_tof );
317 }
318
319 }
320 delete status;
321 }
322 }
323
324 t_j++;
325 }
326
327 }
328
329
330
331
332 int nGood = iGood.size();
333 m_nGood=nGood;
334 m_nCharge=nCharge;
335 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
336 if((nGood != 4)||(nCharge!=0)){
337 return StatusCode::SUCCESS;
338 }
339 Ncut1++;
340
341 for(int i = 0; i < nGood; i++) {
342 m_eop[i] = 5.;
344 if(!(*itTrk)->isMdcTrackValid()) continue;
346 double p = mdcTrk->
p();
347 m_eop[i] = 6.;
348 if(!(*itTrk)->isEmcShowerValid()) continue;
350 double eraw = emcTrk->
energy();
351 m_eop[i]=eraw/p;
352 }
353
355 iGam.clear();
356 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
357 if(i>=evtRecTrkCol->size())break;
359 if(!(*itTrk)->isEmcShowerValid()) continue;
361 Hep3Vector emcpos(emcTrk->
x(), emcTrk->
y(), emcTrk->
z());
362
363 double dthe = 200.;
364 double dphi = 200.;
365 double dang = 200.;
366 for(int j = 0; j < evtRecEvent->totalCharged(); j++) {
368 if(!(*jtTrk)->isExtTrackValid()) continue;
372
373 double angd = extpos.angle(emcpos);
374 double thed = extpos.theta() - emcpos.theta();
375 double phid = extpos.deltaPhi(emcpos);
376 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+
pi, CLHEP::twopi) - CLHEP::pi;
377 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+
pi, CLHEP::twopi) - CLHEP::pi;
378 if(angd < dang){
379 dang = angd;
380 dthe = thed;
381 dphi = phid;
382 }
383 }
384 if(dang>=200) continue;
385 double eraw = emcTrk->
energy();
386 dthe = dthe * 180 / (CLHEP::pi);
387 dphi = dphi * 180 / (CLHEP::pi);
388 dang = dang * 180 / (CLHEP::pi);
389 if(eraw < m_energyThreshold) continue;
390
391 if(fabs(dang) < m_gammaAngleCut) continue;
392 double dtime = emcTrk->
time();
393
394 if(dtime<0) continue;
395 if(dtime>14) continue;
396
397
398
399
400 iGam.push_back((*itTrk)->trackId());
401 }
402
403
404
405
406 int nGam = iGam.size();
407 m_nGam=nGam;
408
409 log << MSG::DEBUG << "num Good Photon " << nGam << " , " <<evtRecEvent->totalNeutral()<<endreq;
410
411
412
413 Ncut2++;
414
415
416
417
418
419
420
421
423 pCh.clear();
425 iPid.clear();
426 iCh.clear();
427 int npi=0,npip=0,npim=0;
428 int nkaon=0,nkaonp=0,nkaonm=0;
429 int nproton=0,nprotonp=0,nprotonm=0;
431
432 for(int i = 0; i < nGood; i++) {
434
437
438
443
444
447
448
449 iCh.push_back(mdcTrk->
charge());
451 {
452 iPid.push_back(0);
453 HepLorentzVector
ptrk;
457 double p3 =
ptrk.vect().mag();
460
461 }
463
464
465
466
467
469
470
471
472 HepLorentzVector
ptrk;
476
478 {npi=npi+1;
479 iPid.push_back(2);
480 if((*itTrk)->isMdcKalTrackValid())
481 {
483 ptrk.setPx(mdcKalTrk->
px());
484 ptrk.setPy(mdcKalTrk->
py());
485 ptrk.setPz(mdcKalTrk->
pz());
486 }
487 double p3 =
ptrk.vect().mag();
490 if(mdcTrk->
charge()>0) npip++;
491 if(mdcTrk->
charge()<0) npim++;
492
493 }
495 {nkaon=nkaon+1;
496 iPid.push_back(3);
497 if((*itTrk)->isMdcKalTrackValid())
498 {
500 ptrk.setPx(mdcKalTrk->
px());
501 ptrk.setPy(mdcKalTrk->
py());
502 ptrk.setPz(mdcKalTrk->
pz());
503 }
504 double p3 =
ptrk.vect().mag();
507 }
509 {nproton=nproton+1;
510 iPid.push_back(4);
511 if((*itTrk)->isMdcKalTrackValid())
512 {
514 ptrk.setPx(mdcKalTrk->
px());
515 ptrk.setPy(mdcKalTrk->
py());
516 ptrk.setPz(mdcKalTrk->
pz());
517 }
518 double p3 =
ptrk.vect().mag();
521 if(mdcTrk->
charge()>0) nprotonp++;
522 if(mdcTrk->
charge()<0) nprotonm++;
523
524 }
525 else
526 {
527 iPid.push_back(0);
528 HepLorentzVector
ptrk;
532 double p3 =
ptrk.vect().mag();
535
536 }
537
538 }
539 m_npi = npi;
540 m_nkaon = nkaon;
541 m_nproton = nproton;
542
543
544
545 m_istat_pmiss=0;
546 m_istat_pbmiss=0;
547 m_istat_pipmiss=0;
548 m_istat_pimmiss=0;
549 for(int i=0;i<4;i++)
550 {
551 m_ptrk_pmiss[i]=5;
552 m_ptrk_pbmiss[i]=5;
553 m_ptrk_pipmiss[i]=5;
554 m_ptrk_pimmiss[i]=5;
555 }
556 for(int i=0; i<3; i++)
557 {
558 m_id_pmiss[i]=0;
559 m_id_pbmiss[i]=0;
560 m_id_pipmiss[i]=0;
561 m_id_pimmiss[i]=0;
562 }
563
564
565
566 HepLorentzVector
ecms(0.03406837,0,0,3.0971873);
567 HepLorentzVector pmiss;
568 m_mpmiss=5.;
569 m_mpbmiss=5.;
570 m_mpipmiss=5.;
571 m_mpimmiss=5.;
572 m_ppmiss=5.;
573 m_ppbmiss=5.;
574 m_ppipmiss=5.;
575 m_ppimmiss=5.;
576
577 int istat_nhit;
578
579 if((npip==1)&&(npim==1)&&(nprotonm==1))
580 {
581 pmiss.setPx(0);
582 pmiss.setPy(0);
583 pmiss.setPz(0);
584 pmiss.setE(0);
585 istat_nhit=0;
586
587 int j=0;
588 for(int i = 0; i < nGood; i++)
589 {
591
593
594 double p;
595 double eraw;
596 if((*itTrk)->isEmcShowerValid())
597 {
600 }
601
602
603 if((iPid[i]*iCh[i])== 2 )
604 {
605 m_index_pmiss[j]=i;
606 pmiss=pmiss+pCh[i];
609 if(m_dedxngoodhit[i]<20) istat_nhit = 1;
610 j++;
611 }
612 else if((iPid[i]*iCh[i])== -2 )
613 {
614 m_index_pmiss[j]=i;
615 pmiss=pmiss+pCh[i];
618 if(m_dedxngoodhit[i]<20) istat_nhit = 1;
619
620 j++;
621 }
622 else if((iPid[i]*iCh[i])== -4 )
623 {
624 m_index_pmiss[j]=i;
625 pmiss=pmiss+pCh[i];
628 if(m_dedxngoodhit[i]<20) istat_nhit = 1;
629
630 j++;
631 }
632 else
633 {
634 if(nGood==4)
635 {
636 m_index_pmiss[3]=i;
639 HepLorentzVector
ptrk;
640 ptrk.setPx(mdcKalTrk->
px());
641 ptrk.setPy(mdcKalTrk->
py());
642 ptrk.setPz(mdcKalTrk->
pz());
643 double p3 =
ptrk.vect().mag();
645 m_ptrk_pmiss[0]=
ptrk.px();
646 m_ptrk_pmiss[1]=
ptrk.py();
647 m_ptrk_pmiss[2]=
ptrk.pz();
648 m_ptrk_pmiss[3]=
ptrk.e();
649
650 }
651 }
652 }
653
654 for(int i = 0; i < nGood; i++) {
655 if((iPid[i]*iCh[i])== 2 ) continue;
656 if((iPid[i]*iCh[i])== -2 ) continue;
657 if((iPid[i]*iCh[i])== -4 ) continue;
658 if(iCh[i]<0) continue;
659
661 for(int ii=0; ii<3; ii++)
662 {
667 if(ii==0)
669 else if(ii==1)
671 else if(ii==2)
673
678 }
679 }
680
681 pmiss =
ecms - pmiss;
682 m_mpmiss= pmiss.m();
683 m_ppmiss = pmiss.rho();
684
685 if(fabs(m_mpmiss-
mproton)<0.02&&istat_nhit==0)
686 {
687 m_istat_pmiss=1;
688 (*(evtRecTrkCol->begin()+ iGood[m_index_pmiss[3]]))->setPartId(3);
689 (*(evtRecTrkCol->begin()+ iGood[m_index_pmiss[3]]))->setQuality(0);
690
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
709
710 double p;
711 double eraw;
712 if((*itTrk)->isEmcShowerValid())
713 {
716 }
717
718
719
720 if((iPid[i]*iCh[i])== 2 )
721 {
722 m_index_pbmiss[j] = i;
723 pmiss = pmiss + pCh[i];
726 if(m_dedxngoodhit[i]<20) istat_nhit=1;
727
728 j++;
729 }
730 else if((iPid[i]*iCh[i])== -2 )
731 {
732 m_index_pbmiss[j] = i;
733
734 pmiss = pmiss + pCh[i];
737 if(m_dedxngoodhit[i]<20) istat_nhit=1;
738
739 j++;
740 }
741 else if((iPid[i]*iCh[i])== 4 )
742 {
743 m_index_pbmiss[j] = i;
744 pmiss = pmiss + pCh[i];
747 if(m_dedxngoodhit[i]<20) istat_nhit=1;
748
749 j++;
750 }
751 else
752 {
753 if(nGood==4)
754 {
755 m_index_pbmiss[3] = i;
758 HepLorentzVector
ptrk;
759 ptrk.setPx(mdcKalTrk->
px());
760 ptrk.setPy(mdcKalTrk->
py());
761 ptrk.setPz(mdcKalTrk->
pz());
762 double p3 =
ptrk.vect().mag();
764 m_ptrk_pbmiss[0]=
ptrk.px();
765 m_ptrk_pbmiss[1]=
ptrk.py();
766 m_ptrk_pbmiss[2]=
ptrk.pz();
767 m_ptrk_pbmiss[3]=
ptrk.e();
768
769 }
770 }
771 }
772
773 for(int i = 0; i < nGood; i++) {
774 if((iPid[i]*iCh[i])== 2 ) continue;
775 if((iPid[i]*iCh[i])== -2 ) continue;
776 if((iPid[i]*iCh[i])== 4 ) continue;
777 if(iCh[i]>0) continue;
778
780 for(int ii=0; ii<3; ii++)
781 {
786 if(ii==0)
788 else if(ii==1)
790 else if(ii==2)
792
797 }
798 }
799 pmiss =
ecms - pmiss;
800 m_mpbmiss= pmiss.m();
801 m_ppbmiss = pmiss.rho();
802 if(fabs(m_mpbmiss-
mproton)<0.02&&istat_nhit==0)
803 {
804 m_istat_pbmiss=1;
805 (*(evtRecTrkCol->begin()+ iGood[m_index_pbmiss[3]]))->setPartId(3);
806 (*(evtRecTrkCol->begin()+ iGood[m_index_pbmiss[3]]))->setQuality(0);
807 }
808 Ncut4++;
809 }
810 int initial_pip, initial_pim;
811
812
813
814
815 if((npim==1)&&(nprotonp==1)&&(nprotonm==1))
816 {
817 pmiss.setPx(0);
818 pmiss.setPy(0);
819 pmiss.setPz(0);
820 pmiss.setE(0);
821 istat_nhit=0;
822
823 int j=0;
824 for(int i = 0; i < nGood; i++)
825 {
827
829
830 double p;
831 double eraw;
832 if((*itTrk)->isEmcShowerValid())
833 {
836 }
837
838
839
840 if((iPid[i]*iCh[i])== 4 )
841 {
842 m_index_pipmiss[j] = i;
843 pmiss = pmiss +pCh[i];
846 if(m_dedxngoodhit[i]<20) istat_nhit =1;
847
848 j++;
849 }
850 else if((iPid[i]*iCh[i])== -4 )
851 {
852 m_index_pipmiss[j] = i;
853
854 pmiss = pmiss +pCh[i];
857 if(m_dedxngoodhit[i]<20) istat_nhit=1;
858
859 j++;
860 }
861 else if((iPid[i]*iCh[i])== -2 )
862 {
863 m_index_pipmiss[j] = i;
864
865 pmiss = pmiss +pCh[i];
868 if(m_dedxngoodhit[i]<20) istat_nhit=1;
869
870 j++;
871 }
872 else
873 {
874 if(nGood==4)
875 {
876 m_index_pipmiss[3] = i;
877
880 HepLorentzVector
ptrk;
881 ptrk.setPx(mdcKalTrk->
px());
882 ptrk.setPy(mdcKalTrk->
py());
883 ptrk.setPz(mdcKalTrk->
pz());
884 double p3 =
ptrk.vect().mag();
886 m_ptrk_pipmiss[0]=
ptrk.px();
887 m_ptrk_pipmiss[1]=
ptrk.py();
888 m_ptrk_pipmiss[2]=
ptrk.pz();
889 m_ptrk_pipmiss[3]=
ptrk.e();
890
891 }
892 }
893 }
894
895 for(int i = 0; i < nGood; i++) {
896 if((iPid[i]*iCh[i])== 4 ) continue;
897 if((iPid[i]*iCh[i])== -4 ) continue;
898 if((iPid[i]*iCh[i])== -2 ) continue;
899 if(iCh[i]<0) continue;
900
902 for(int ii=0; ii<3; ii++)
903 {
908 if(ii==0)
910 else if(ii==1)
912 else if(ii==2)
914
919 }
920 }
921 pmiss =
ecms - pmiss;
922 m_mpipmiss = pmiss.m();
923 m_ppipmiss = pmiss.rho();
924 if(fabs(m_mpipmiss-
mpi)<0.05&&istat_nhit==0)
925 {
926 m_istat_pipmiss=1;
927 (*(evtRecTrkCol->begin()+ iGood[m_index_pipmiss[3]]))->setPartId(2);
928 (*(evtRecTrkCol->begin()+ iGood[m_index_pipmiss[3]]))->setQuality(0);
929
930 }
931 Ncut5++;
932 }
933 if((npip==1)&&(nprotonp==1)&&(nprotonm==1))
934 {
935 pmiss.setPx(0);
936 pmiss.setPy(0);
937 pmiss.setPz(0);
938 pmiss.setE(0);
939 istat_nhit=0;
940
941 int j=0;
942 for(int i = 0; i < nGood; i++)
943 {
945
947
948 double p;
949 double eraw;
950 if((*itTrk)->isEmcShowerValid())
951 {
954 }
955
956
957
958
959
960 if((iPid[i]*iCh[i])== 4 )
961 {
962 m_index_pimmiss[j]=i;
963 pmiss = pmiss + pCh[i];
966 if(m_dedxngoodhit[i]<20) istat_nhit=1;
967
968 j++;
969 }
970 else if((iPid[i]*iCh[i])== -4 )
971 {
972 m_index_pimmiss[j]=i;
973
974 pmiss = pmiss + pCh[i];
977 if(m_dedxngoodhit[i]<20) istat_nhit =1;
978
979 j++;
980 }
981 else if((iPid[i]*iCh[i])== 2 )
982 {
983 m_index_pimmiss[j]=i;
984 pmiss = pmiss + pCh[i];
987 if(m_dedxngoodhit[i]<20) istat_nhit = 1;
988
989 j++;
990 }
991 else
992 {
993 if(nGood==4)
994 {
995 m_index_pimmiss[3]=i;
998 HepLorentzVector
ptrk;
999 ptrk.setPx(mdcKalTrk->
px());
1000 ptrk.setPy(mdcKalTrk->
py());
1001 ptrk.setPz(mdcKalTrk->
pz());
1002 double p3 =
ptrk.vect().mag();
1004 m_ptrk_pimmiss[0]=
ptrk.px();
1005 m_ptrk_pimmiss[1]=
ptrk.py();
1006 m_ptrk_pimmiss[2]=
ptrk.pz();
1007 m_ptrk_pimmiss[3]=
ptrk.e();
1008
1009 }
1010 }
1011 }
1012
1013 for(int i = 0; i < nGood; i++) {
1014 if((iPid[i]*iCh[i])== 4 ) continue;
1015 if((iPid[i]*iCh[i])== -4 ) continue;
1016 if((iPid[i]*iCh[i])== 2 ) continue;
1017 if(iCh[i]>0) continue;
1018
1020 for(int ii=0; ii<3; ii++)
1021 {
1026 if(ii==0)
1028 else if(ii==1)
1030 else if(ii==2)
1032
1037 }
1038 }
1039 pmiss =
ecms - pmiss;
1040 m_mpimmiss = pmiss.m();
1041 m_ppimmiss = pmiss.rho();
1042 if(fabs(m_mpimmiss-
mpi)<0.05&&istat_nhit==0)
1043 {
1044 m_istat_pimmiss=1;
1045 (*(evtRecTrkCol->begin()+ iGood[m_index_pimmiss[3]]))->setPartId(2);
1046 (*(evtRecTrkCol->begin()+ iGood[m_index_pimmiss[3]]))->setQuality(0);
1047
1048 }
1049 Ncut6++;
1050 }
1051
1052
1053
1054 if(m_istat_pmiss!=1 &&m_istat_pbmiss!=1&&m_istat_pipmiss!=1 && m_istat_pimmiss!=1)
1055 {return StatusCode::SUCCESS;}
1056
1057
1058
1059 if(m_saventuple) m_tuple0->write();
1060
1061 setFilterPassed(true);
1062
1063 return StatusCode::SUCCESS;
1064}
double sin(const BesAngle a)
double cos(const BesAngle a)
EvtRecTrackCol::iterator EvtRecTrackIterator
std::vector< HepLorentzVector > Vp4
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