81 MsgStream log(
msgSvc(), name());
83 log << MSG::INFO <<
"in initialize()" << endmsg;
84 Ncut0=Ncut1=Ncut2=Ncut3=Ncut4=Ncut5=Ncut6=0;
89 NTuplePtr nt0(
ntupleSvc(),
"DQAFILE/pi2p2");
90 if ( nt0 ) m_tuple0 = nt0;
92 m_tuple0 =
ntupleSvc()->book (
"DQAFILE/pi2p2", CLID_ColumnWiseTuple,
"ks N-Tuple example");
94 status = m_tuple0->addItem (
"nrun", m_nrun);
95 status = m_tuple0->addItem (
"nrec", m_nrec);
97 status = m_tuple0->addItem (
"nGam", m_nGam);
98 status = m_tuple0->addItem (
"nGood", m_nGood);
99 status = m_tuple0->addItem (
"nCharge", m_nCharge);
100 status = m_tuple0->addItem (
"npi", m_npi);
101 status = m_tuple0->addItem (
"nkaon", m_nkaon);
102 status = m_tuple0->addItem (
"nproton", m_nproton);
104 status = m_tuple0->addItem (
"dedxchi_e", 4, m_dedxchi_e);
105 status = m_tuple0->addItem (
"dedxchi_mu", 4, m_dedxchi_mu);
106 status = m_tuple0->addItem (
"dedxchi_pi", 4, m_dedxchi_pi);
107 status = m_tuple0->addItem (
"dedxchi_kaon", 4, m_dedxchi_kaon);
108 status = m_tuple0->addItem (
"dedxchi_proton", 4, m_dedxchi_proton);
110 status = m_tuple0->addItem (
"tofchi_e", 4, m_tofchi_e);
111 status = m_tuple0->addItem (
"tofchi_mu", 4, m_tofchi_mu);
112 status = m_tuple0->addItem (
"tofchi_pi", 4, m_tofchi_pi);
113 status = m_tuple0->addItem (
"tofchi_kaon", 4, m_tofchi_kaon);
114 status = m_tuple0->addItem (
"tofchi_proton", 4, m_tofchi_proton);
116 status = m_tuple0->addItem (
"trackfitchi", 4, m_trackfitchi);
117 status = m_tuple0->addItem (
"dedxngoodhit", 4, m_dedxngoodhit);
118 status = m_tuple0->addItem (
"trackfitndof", 4, m_trackfitndof);
120 status = m_tuple0->addItem (
"index_pmiss", 4, m_index_pmiss);
121 status = m_tuple0->addItem (
"index_pbmiss", 4, m_index_pbmiss);
122 status = m_tuple0->addItem (
"index_pipmiss", 4, m_index_pipmiss);
123 status = m_tuple0->addItem (
"index_pimmiss", 4, m_index_pimmiss);
126 status = m_tuple0->addItem (
"istat_pmiss", m_istat_pmiss);
127 status = m_tuple0->addItem (
"istat_pbmiss", m_istat_pbmiss);
128 status = m_tuple0->addItem (
"istat_pipmiss", m_istat_pipmiss);
129 status = m_tuple0->addItem (
"istat_pimmiss", m_istat_pimmiss);
131 status = m_tuple0->addItem (
"mpmiss", m_mpmiss);
132 status = m_tuple0->addItem (
"mpbmiss", m_mpbmiss);
133 status = m_tuple0->addItem (
"mpipmiss", m_mpipmiss);
134 status = m_tuple0->addItem (
"mpimmiss", m_mpimmiss);
136 status = m_tuple0->addItem (
"ppmiss", m_ppmiss);
137 status = m_tuple0->addItem (
"ppbmiss", m_ppbmiss);
138 status = m_tuple0->addItem (
"ppipmiss", m_ppipmiss);
139 status = m_tuple0->addItem (
"ppimmiss", m_ppimmiss);
142 status = m_tuple0->addItem (
"ptrk_pmiss", 4, m_ptrk_pmiss);
143 status = m_tuple0->addItem (
"ptrk_pbmiss", 4, m_ptrk_pbmiss);
144 status = m_tuple0->addItem (
"ptrk_pipmiss", 4, m_ptrk_pipmiss);
145 status = m_tuple0->addItem (
"ptrk_pimmiss", 4, m_ptrk_pimmiss);
147 status = m_tuple0->addItem (
"id_pmiss", 3, m_id_pmiss);
148 status = m_tuple0->addItem (
"id_pbmiss", 3, m_id_pbmiss);
149 status = m_tuple0->addItem (
"id_pipmiss", 3, m_id_pipmiss);
150 status = m_tuple0->addItem (
"id_pimmiss", 3, m_id_pimmiss);
153 status = m_tuple0->addItem (
"eop", 4, m_eop);
159 log << MSG::ERROR <<
" Cannot book N-tuple:" << long(m_tuple0) << endmsg;
160 return StatusCode::FAILURE;
170 log << MSG::INFO <<
"successfully return from initialize()" <<endmsg;
171 return StatusCode::SUCCESS;
180 MsgStream log(
msgSvc(), name());
181 log << MSG::INFO <<
"in execute()" << endreq;
183 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
184 int runNo=eventHeader->runNumber();
185 int event=eventHeader->eventNumber();
186 log << MSG::DEBUG <<
"run, evtnum = "
192 setFilterPassed(
false);
194 if ((event%100000)==0) {cout <<
"event: "<<
event << endl ; }
199 log << MSG::DEBUG <<
"ncharg, nneu, tottks = "
200 << evtRecEvent->totalCharged() <<
" , "
201 << evtRecEvent->totalNeutral() <<
" , "
202 << evtRecEvent->totalTracks() <<endreq;
206 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),
"/Event/MC/McParticleCol");
207 double temp_p_pp,temp_p_pm;
215 Vint iGood, ipip, ipim , ikaonp, ikaonm, iprotonp,iprotonm;
229 Hep3Vector xorigin(0,0,0);
233 Gaudi::svcLocator()->service(
"VertexDbSvc", vtxsvc);
239 xorigin.setX(dbv[0]);
240 xorigin.setY(dbv[1]);
241 xorigin.setZ(dbv[2]);
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);
262 HepVector a = mdcTrk->
helix();
263 HepSymMatrix Ea = mdcTrk->
err();
265 HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]);
268 HepVector vecipa = helixip.
a();
269 double Rvxy0=fabs(vecipa[0]);
270 double Rvz0=vecipa[3];
271 double Rvphi0=vecipa[1];
273 if(fabs(Rvz0) >= m_vz0cut)
continue;
274 if(fabs(Rvxy0) >= m_vr0cut)
continue;
275 if(
cos(mdcTrk->
theta())>0.93)
continue;
279 nCharge += mdcTrk->
charge();
282 if((*itTrk)->isMdcDedxValid())
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();
294 if((*itTrk)->isTofTrackValid())
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();
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 );
332 int nGood = iGood.size();
335 log << MSG::DEBUG <<
"ngood, totcharge = " << nGood <<
" , " << nCharge << endreq;
336 if((nGood != 4)||(nCharge!=0)){
337 return StatusCode::SUCCESS;
341 for(
int i = 0; i < nGood; i++) {
344 if(!(*itTrk)->isMdcTrackValid())
continue;
346 double p = mdcTrk->
p();
348 if(!(*itTrk)->isEmcShowerValid())
continue;
350 double eraw = emcTrk->
energy();
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());
366 for(
int j = 0; j < evtRecEvent->totalCharged(); j++) {
368 if(!(*jtTrk)->isExtTrackValid())
continue;
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;
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;
391 if(fabs(dang) < m_gammaAngleCut)
continue;
392 double dtime = emcTrk->
time();
394 if(dtime<0)
continue;
395 if(dtime>14)
continue;
400 iGam.push_back((*itTrk)->trackId());
406 int nGam = iGam.size();
409 log << MSG::DEBUG <<
"num Good Photon " << nGam <<
" , " <<evtRecEvent->totalNeutral()<<endreq;
427 int npi=0,npip=0,npim=0;
428 int nkaon=0,nkaonp=0,nkaonm=0;
429 int nproton=0,nprotonp=0,nprotonm=0;
432 for(
int i = 0; i < nGood; i++) {
449 iCh.push_back(mdcTrk->
charge());
453 HepLorentzVector
ptrk;
457 double p3 =
ptrk.vect().mag();
472 HepLorentzVector
ptrk;
480 if((*itTrk)->isMdcKalTrackValid())
483 ptrk.setPx(mdcKalTrk->
px());
484 ptrk.setPy(mdcKalTrk->
py());
485 ptrk.setPz(mdcKalTrk->
pz());
487 double p3 =
ptrk.vect().mag();
490 if(mdcTrk->
charge()>0) npip++;
491 if(mdcTrk->
charge()<0) npim++;
497 if((*itTrk)->isMdcKalTrackValid())
500 ptrk.setPx(mdcKalTrk->
px());
501 ptrk.setPy(mdcKalTrk->
py());
502 ptrk.setPz(mdcKalTrk->
pz());
504 double p3 =
ptrk.vect().mag();
511 if((*itTrk)->isMdcKalTrackValid())
514 ptrk.setPx(mdcKalTrk->
px());
515 ptrk.setPy(mdcKalTrk->
py());
516 ptrk.setPz(mdcKalTrk->
pz());
518 double p3 =
ptrk.vect().mag();
521 if(mdcTrk->
charge()>0) nprotonp++;
522 if(mdcTrk->
charge()<0) nprotonm++;
528 HepLorentzVector
ptrk;
532 double p3 =
ptrk.vect().mag();
556 for(
int i=0; i<3; i++)
566 HepLorentzVector
ecms(0.03406837,0,0,3.0971873);
567 HepLorentzVector pmiss;
579 if((npip==1)&&(npim==1)&&(nprotonm==1))
588 for(
int i = 0; i < nGood; i++)
596 if((*itTrk)->isEmcShowerValid())
603 if((iPid[i]*iCh[i])== 2 )
609 if(m_dedxngoodhit[i]<20) istat_nhit = 1;
612 else if((iPid[i]*iCh[i])== -2 )
618 if(m_dedxngoodhit[i]<20) istat_nhit = 1;
622 else if((iPid[i]*iCh[i])== -4 )
628 if(m_dedxngoodhit[i]<20) istat_nhit = 1;
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();
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;
661 for(
int ii=0; ii<3; ii++)
681 pmiss =
ecms - pmiss;
683 m_ppmiss = pmiss.rho();
685 if(fabs(m_mpmiss-
mproton)<0.02&&istat_nhit==0)
688 (*(evtRecTrkCol->begin()+ iGood[m_index_pmiss[3]]))->setPartId(3);
689 (*(evtRecTrkCol->begin()+ iGood[m_index_pmiss[3]]))->setQuality(0);
695 if((npip==1)&&(npim==1)&&(nprotonp==1))
704 for(
int i = 0; i < nGood; i++)
712 if((*itTrk)->isEmcShowerValid())
720 if((iPid[i]*iCh[i])== 2 )
722 m_index_pbmiss[j] = i;
723 pmiss = pmiss + pCh[i];
726 if(m_dedxngoodhit[i]<20) istat_nhit=1;
730 else if((iPid[i]*iCh[i])== -2 )
732 m_index_pbmiss[j] = i;
734 pmiss = pmiss + pCh[i];
737 if(m_dedxngoodhit[i]<20) istat_nhit=1;
741 else if((iPid[i]*iCh[i])== 4 )
743 m_index_pbmiss[j] = i;
744 pmiss = pmiss + pCh[i];
747 if(m_dedxngoodhit[i]<20) istat_nhit=1;
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();
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;
780 for(
int ii=0; ii<3; ii++)
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)
805 (*(evtRecTrkCol->begin()+ iGood[m_index_pbmiss[3]]))->setPartId(3);
806 (*(evtRecTrkCol->begin()+ iGood[m_index_pbmiss[3]]))->setQuality(0);
810 int initial_pip, initial_pim;
815 if((npim==1)&&(nprotonp==1)&&(nprotonm==1))
824 for(
int i = 0; i < nGood; i++)
832 if((*itTrk)->isEmcShowerValid())
840 if((iPid[i]*iCh[i])== 4 )
842 m_index_pipmiss[j] = i;
843 pmiss = pmiss +pCh[i];
846 if(m_dedxngoodhit[i]<20) istat_nhit =1;
850 else if((iPid[i]*iCh[i])== -4 )
852 m_index_pipmiss[j] = i;
854 pmiss = pmiss +pCh[i];
857 if(m_dedxngoodhit[i]<20) istat_nhit=1;
861 else if((iPid[i]*iCh[i])== -2 )
863 m_index_pipmiss[j] = i;
865 pmiss = pmiss +pCh[i];
868 if(m_dedxngoodhit[i]<20) istat_nhit=1;
876 m_index_pipmiss[3] = i;
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();
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;
902 for(
int ii=0; ii<3; ii++)
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)
927 (*(evtRecTrkCol->begin()+ iGood[m_index_pipmiss[3]]))->setPartId(2);
928 (*(evtRecTrkCol->begin()+ iGood[m_index_pipmiss[3]]))->setQuality(0);
933 if((npip==1)&&(nprotonp==1)&&(nprotonm==1))
942 for(
int i = 0; i < nGood; i++)
950 if((*itTrk)->isEmcShowerValid())
960 if((iPid[i]*iCh[i])== 4 )
962 m_index_pimmiss[j]=i;
963 pmiss = pmiss + pCh[i];
966 if(m_dedxngoodhit[i]<20) istat_nhit=1;
970 else if((iPid[i]*iCh[i])== -4 )
972 m_index_pimmiss[j]=i;
974 pmiss = pmiss + pCh[i];
977 if(m_dedxngoodhit[i]<20) istat_nhit =1;
981 else if((iPid[i]*iCh[i])== 2 )
983 m_index_pimmiss[j]=i;
984 pmiss = pmiss + pCh[i];
987 if(m_dedxngoodhit[i]<20) istat_nhit = 1;
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();
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;
1020 for(
int ii=0; ii<3; ii++)
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)
1045 (*(evtRecTrkCol->begin()+ iGood[m_index_pimmiss[3]]))->setPartId(2);
1046 (*(evtRecTrkCol->begin()+ iGood[m_index_pimmiss[3]]))->setQuality(0);
1054 if(m_istat_pmiss!=1 &&m_istat_pbmiss!=1&&m_istat_pipmiss!=1 && m_istat_pimmiss!=1)
1055 {
return StatusCode::SUCCESS;}
1059 if(m_saventuple) m_tuple0->write();
1061 setFilterPassed(
true);
1063 return StatusCode::SUCCESS;