1 #include "GaudiKernel/MsgStream.h"
2 #include "GaudiKernel/AlgFactory.h"
3 #include "GaudiKernel/SmartDataPtr.h"
4 #include "GaudiKernel/IDataProviderSvc.h"
5 #include "GaudiKernel/PropertyMgr.h"
17 #include "GaudiKernel/INTupleSvc.h"
18 #include "GaudiKernel/NTuple.h"
19 #include "GaudiKernel/Bootstrap.h"
20 #include "GaudiKernel/IHistogramSvc.h"
21 #include "CLHEP/Vector/ThreeVector.h"
22 #include "CLHEP/Vector/LorentzVector.h"
23 #include "CLHEP/Vector/TwoVector.h"
36 #include "GaudiKernel/Bootstrap.h"
37 #include "GaudiKernel/ISvcLocator.h"
38 using CLHEP::Hep3Vector;
39 using CLHEP::Hep2Vector;
40 using CLHEP::HepLorentzVector;
41 #include "CLHEP/Geometry/Point3D.h"
42 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
50 typedef std::vector<int>
Vint;
51 typedef std::vector<HepLorentzVector>
Vp4;
53 const double mpi = 0.13957;
59 double phi= trk->
fi0();
60 double kappa= trk->
kappa();
61 double tanl = trk->
tanl();
63 double px=-
sin(phi)/fabs(kappa);
68 double phi= trk->
fi0();
69 double kappa= trk->
kappa();
70 double tanl = trk->
tanl();
72 double py=
cos(phi)/fabs(kappa);
79 double phi= trk->
fi0();
80 double kappa= trk->
kappa();
81 double tanl = trk->
tanl();
83 double pz=tanl/fabs(kappa);
88 double phi= trk->
fi0();
89 double kappa= trk->
kappa();
90 double tanl = trk->
tanl();
92 double px=-
sin(phi)/fabs(kappa);
93 double py=
cos(phi)/fabs(kappa);
94 double pz=tanl/fabs(kappa);
96 return sqrt(px*px+py*py+pz*pz);
100 double phi0=trk->
fi0();
101 double kappa = trk->
kappa();
103 if(kappa!=0) pxy = 1.0/fabs(kappa);
106 double px = pxy * (-
sin(phi0));
107 double py = pxy *
cos(phi0);
114 int CalibEventSelect::idmax[43]={40,44,48,56,64,72,80,80,76,76,
115 88,88,100,100,112,112,128,128,140,140,
116 160,160,160,160,176,176,176,176,208,208,
117 208,208,240,240,240,240,256,256,256,256,
122 Algorithm(name, pSvcLocator) {
125 declareProperty(
"Output", m_output =
false);
126 declareProperty(
"Display", m_display =
false);
127 declareProperty(
"PrintMonitor", m_printmonitor=
false);
128 declareProperty(
"SelectRadBhabha", m_selectRadBhabha=
false);
129 declareProperty(
"SelectGGEE", m_selectGGEE=
false);
130 declareProperty(
"SelectGG4Pi", m_selectGG4Pi=
false);
131 declareProperty(
"SelectRadBhabhaT", m_selectRadBhabhaT=
false);
132 declareProperty(
"SelectRhoPi", m_selectRhoPi=
false);
133 declareProperty(
"SelectKstarK", m_selectKstarK=
false);
134 declareProperty(
"SelectPPPiPi", m_selectPPPiPi=
false);
135 declareProperty(
"SelectRecoJpsi", m_selectRecoJpsi=
false);
136 declareProperty(
"SelectBhabha", m_selectBhabha=
false);
137 declareProperty(
"SelectDimu", m_selectDimu=
false);
138 declareProperty(
"SelectCosmic", m_selectCosmic=
false);
139 declareProperty(
"SelectGenProton", m_selectGenProton=
false);
140 declareProperty(
"SelectPsipRhoPi", m_selectPsipRhoPi=
false);
141 declareProperty(
"SelectPsipKstarK", m_selectPsipKstarK=
false);
142 declareProperty(
"SelectPsipPPPiPi", m_selectPsipPPPiPi=
false);
143 declareProperty(
"SelectPsippCand", m_selectPsippCand=
false);
145 declareProperty(
"BhabhaScale", m_radbhabha_scale=1000);
146 declareProperty(
"RadBhabhaScale", m_bhabha_scale=1000);
147 declareProperty(
"DimuScale", m_dimu_scale=10);
148 declareProperty(
"CosmicScale", m_cosmic_scale=3);
149 declareProperty(
"ProtonScale", m_proton_scale=100);
151 declareProperty(
"CosmicLowp", m_cosmic_lowp=1.0);
153 declareProperty(
"WriteDst", m_writeDst=
false);
154 declareProperty(
"WriteRec", m_writeRec=
false);
155 declareProperty(
"Ecm", m_ecm=3.1);
156 declareProperty(
"Vr0cut", m_vr0cut=1.0);
157 declareProperty(
"Vz0cut", m_vz0cut=10.0);
158 declareProperty(
"Pt0HighCut", m_pt0HighCut=5.0);
159 declareProperty(
"Pt0LowCut", m_pt0LowCut=0.05);
160 declareProperty(
"EnergyThreshold", m_energyThreshold=0.05);
161 declareProperty(
"GammaPhiCut", m_gammaPhiCut=20.0);
162 declareProperty(
"GammaThetaCut", m_gammaThetaCut=20.0);
169 MsgStream log(
msgSvc(), name());
171 log << MSG::INFO <<
"in initialize()" << endmsg;
177 h_nGood=
histoSvc()->book(
"1D/nGoodtracks", 1,
"num of good chaged tracks", 20, 0, 20);
178 h_nCharge=
histoSvc()->book(
"1D/nCharge", 1,
"net charge", 20, -10, 10);
179 h_pmaxobeam=
histoSvc()->book(
"1D/pmax", 1,
"pmax over beam ratio", 100, 0, 3);
180 h_eopmax=
histoSvc()->book(
"1D/eopmax", 1,
"e over pmax ratio", 100, 0, 3);
181 h_eopsec=
histoSvc()->book(
"1D/eopsec", 1,
"e over psec ratio", 100, 0, 3);
182 h_deltap=
histoSvc()->book(
"1D/deltap", 1,
"pmax minus psec", 100, -3, 3);
183 h_esumoecm=
histoSvc()->book(
"1D/esumoverecm", 1,
"total energy over ecm ratio", 100, 0, 3);
184 h_egmax=
histoSvc()->book(
"1D/egmax", 1,
"most energetic photon energy", 100, 0, 3);
185 h_deltaphi=
histoSvc()->book(
"1D/deltaphi", 1,
"phi_pmax minus phi_sec", 150, -4, 4);
186 h_Pt=
histoSvc()->book(
"1D/Pt", 1,
"total Transverse Momentum", 200,-1, 4);
190 log << MSG::INFO <<
"creating sub-algorithms...." << endreq;
197 sc = createSubAlgorithm(
"EventWriter",
"WriteDst", m_subalg1);
198 if( sc.isFailure() ) {
199 log << MSG::ERROR <<
"Error creating Sub-Algorithm WriteDst" <<endreq;
202 log << MSG::INFO <<
"Success creating Sub-Algorithm WriteDst" <<endreq;
207 sc = createSubAlgorithm(
"EventWriter",
"WriteRec", m_subalg2);
208 if( sc.isFailure() ) {
209 log << MSG::ERROR <<
"Error creating Sub-Algorithm WriteRec" <<endreq;
212 log << MSG::INFO <<
"Success creating Sub-Algorithm WriteRec" <<endreq;
217 if(m_selectRadBhabha) {
218 sc = createSubAlgorithm(
"EventWriter",
"SelectRadBhabha", m_subalg3);
219 if( sc.isFailure() ) {
220 log << MSG::ERROR <<
"Error creating Sub-Algorithm SelectRadBhabha" <<endreq;
223 log << MSG::INFO <<
"Success creating Sub-Algorithm SelectRadBhabha" <<endreq;
228 sc = createSubAlgorithm(
"EventWriter",
"SelectGGEE", m_subalg4);
229 if( sc.isFailure() ) {
230 log << MSG::ERROR <<
"Error creating Sub-Algorithm SelectGGEE" <<endreq;
233 log << MSG::INFO <<
"Success creating Sub-Algorithm SelectGGEE" <<endreq;
238 sc = createSubAlgorithm(
"EventWriter",
"SelectGG4Pi", m_subalg5);
239 if( sc.isFailure() ) {
240 log << MSG::ERROR <<
"Error creating Sub-Algorithm SelectGG4Pi" <<endreq;
243 log << MSG::INFO <<
"Success creating Sub-Algorithm SelectGG4Pi" <<endreq;
248 if(m_selectRadBhabhaT) {
249 sc = createSubAlgorithm(
"EventWriter",
"SelectRadBhabhaT", m_subalg6);
250 if( sc.isFailure() ) {
251 log << MSG::ERROR <<
"Error creating Sub-Algorithm SelectRadBhabhaT" <<endreq;
254 log << MSG::INFO <<
"Success creating Sub-Algorithm SelectRadBhabhaT" <<endreq;
260 sc = createSubAlgorithm(
"EventWriter",
"SelectRhoPi", m_subalg7);
261 if( sc.isFailure() ) {
262 log << MSG::ERROR <<
"Error creating Sub-Algorithm SelectRhoPi" <<endreq;
265 log << MSG::INFO <<
"Success creating Sub-Algorithm SelectRhoPi" <<endreq;
271 sc = createSubAlgorithm(
"EventWriter",
"SelectKstarK", m_subalg8);
272 if( sc.isFailure() ) {
273 log << MSG::ERROR <<
"Error creating Sub-Algorithm SelectKstarK" <<endreq;
276 log << MSG::INFO <<
"Success creating Sub-Algorithm SelectKstarK" <<endreq;
283 sc = createSubAlgorithm(
"EventWriter",
"SelectPPPiPi", m_subalg9);
284 if( sc.isFailure() ) {
285 log << MSG::ERROR <<
"Error creating Sub-Algorithm SelectPPPiPi" <<endreq;
288 log << MSG::INFO <<
"Success creating Sub-Algorithm SelectPPPiPi" <<endreq;
293 if(m_selectRecoJpsi) {
294 sc = createSubAlgorithm(
"EventWriter",
"SelectRecoJpsi", m_subalg10);
295 if( sc.isFailure() ) {
296 log << MSG::ERROR <<
"Error creating Sub-Algorithm SelectRecoJpsi" <<endreq;
299 log << MSG::INFO <<
"Success creating Sub-Algorithm SelectRecoJpsi" <<endreq;
305 sc = createSubAlgorithm(
"EventWriter",
"SelectBhabha", m_subalg11);
306 if( sc.isFailure() ) {
307 log << MSG::ERROR <<
"Error creating Sub-Algorithm SelectBhabha" <<endreq;
310 log << MSG::INFO <<
"Success creating Sub-Algorithm SelectBhabha" <<endreq;
315 sc = createSubAlgorithm(
"EventWriter",
"SelectDimu", m_subalg12);
316 if( sc.isFailure() ) {
317 log << MSG::ERROR <<
"Error creating Sub-Algorithm SelectDimu" <<endreq;
320 log << MSG::INFO <<
"Success creating Sub-Algorithm SelectDimu" <<endreq;
325 sc = createSubAlgorithm(
"EventWriter",
"SelectCosmic", m_subalg13);
326 if( sc.isFailure() ) {
327 log << MSG::ERROR <<
"Error creating Sub-Algorithm SelectCosmic" <<endreq;
330 log << MSG::INFO <<
"Success creating Sub-Algorithm SelectCosmic" <<endreq;
334 if(m_selectGenProton) {
335 sc = createSubAlgorithm(
"EventWriter",
"SelectGenProton", m_subalg14);
336 if( sc.isFailure() ) {
337 log << MSG::ERROR <<
"Error creating Sub-Algorithm SelectGenProton" <<endreq;
340 log << MSG::INFO <<
"Success creating Sub-Algorithm SelectGenProton" <<endreq;
345 if(m_selectPsipRhoPi) {
346 sc = createSubAlgorithm(
"EventWriter",
"SelectPsipRhoPi", m_subalg15);
347 if( sc.isFailure() ) {
348 log << MSG::ERROR <<
"Error creating Sub-Algorithm SelectPsipRhoPi" <<endreq;
351 log << MSG::INFO <<
"Success creating Sub-Algorithm SelectPsipRhoPi" <<endreq;
356 if(m_selectPsipKstarK) {
357 sc = createSubAlgorithm(
"EventWriter",
"SelectPsipKstarK", m_subalg16);
358 if( sc.isFailure() ) {
359 log << MSG::ERROR <<
"Error creating Sub-Algorithm SelectPsipKstarK" <<endreq;
362 log << MSG::INFO <<
"Success creating Sub-Algorithm SelectPsipKstarK" <<endreq;
367 if(m_selectPsipPPPiPi) {
368 sc = createSubAlgorithm(
"EventWriter",
"SelectPsipPPPiPi", m_subalg17);
369 if( sc.isFailure() ) {
370 log << MSG::ERROR <<
"Error creating Sub-Algorithm SelectPsipPPPiPi" <<endreq;
373 log << MSG::INFO <<
"Success creating Sub-Algorithm SelectPsipPPPiPi" <<endreq;
378 if(m_selectPsippCand) {
379 sc = createSubAlgorithm(
"EventWriter",
"SelectPsippCand", m_subalg18);
380 if( sc.isFailure() ) {
381 log << MSG::ERROR <<
"Error creating Sub-Algorithm SelectPsippCand" <<endreq;
384 log << MSG::INFO <<
"Success creating Sub-Algorithm SelectPsippCand" <<endreq;
393 NTuplePtr nt0(
ntupleSvc(),
"FILE1/hadron");
394 if ( nt0 ) m_tuple0 = nt0;
396 m_tuple0 =
ntupleSvc()->book (
"FILE1/hadron", CLID_ColumnWiseTuple,
"N-Tuple example");
398 status = m_tuple0->addItem (
"esum", m_esum);
399 status = m_tuple0->addItem (
"eemc", m_eemc);
400 status = m_tuple0->addItem (
"etot", m_etot);
401 status = m_tuple0->addItem (
"nGood", m_nGood);
402 status = m_tuple0->addItem (
"nCharge", m_nCharge);
403 status = m_tuple0->addItem (
"nGam", m_nGam);
404 status = m_tuple0->addItem (
"ptot", m_ptot);
405 status = m_tuple0->addItem (
"pp", m_pp);
406 status = m_tuple0->addItem (
"pm", m_pm);
407 status = m_tuple0->addItem (
"run", m_runnb);
408 status = m_tuple0->addItem (
"event", m_evtnb);
409 status = m_tuple0->addItem (
"maxE", m_maxE);
410 status = m_tuple0->addItem (
"secE", m_secE);
411 status = m_tuple0->addItem (
"dthe", m_dthe);
412 status = m_tuple0->addItem (
"dphi", m_dphi);
413 status = m_tuple0->addItem (
"mdcHit1", m_mdcHit1);
414 status = m_tuple0->addItem (
"mdcHit2", m_mdcHit2);
417 log << MSG::ERROR <<
" Cannot book N-tuple:" << long(m_tuple0) << endmsg;
418 return StatusCode::FAILURE;
422 NTuplePtr nt1(
ntupleSvc(),
"FILE1/vxyz");
423 if ( nt1 ) m_tuple1 = nt1;
425 m_tuple1 =
ntupleSvc()->book (
"FILE1/vxyz", CLID_ColumnWiseTuple,
"ks N-Tuple example");
427 status = m_tuple1->addItem (
"vx0", m_vx0);
428 status = m_tuple1->addItem (
"vy0", m_vy0);
429 status = m_tuple1->addItem (
"vz0", m_vz0);
430 status = m_tuple1->addItem (
"vr0", m_vr0);
431 status = m_tuple1->addItem (
"theta0", m_theta0);
432 status = m_tuple1->addItem (
"p0", m_p0);
433 status = m_tuple1->addItem (
"pt0", m_pt0);
436 log << MSG::ERROR <<
" Cannot book N-tuple:" << long(m_tuple1) << endmsg;
437 return StatusCode::FAILURE;
441 NTuplePtr nt2(
ntupleSvc(),
"FILE1/photon");
442 if ( nt2 ) m_tuple2 = nt2;
444 m_tuple2 =
ntupleSvc()->book (
"FILE1/photon", CLID_ColumnWiseTuple,
"ks N-Tuple example");
446 status = m_tuple2->addItem (
"dthe", m_dthe);
447 status = m_tuple2->addItem (
"dphi", m_dphi);
448 status = m_tuple2->addItem (
"dang", m_dang);
449 status = m_tuple2->addItem (
"eraw", m_eraw);
452 log << MSG::ERROR <<
" Cannot book N-tuple:" << long(m_tuple2) << endmsg;
453 return StatusCode::FAILURE;
465 m_radBhabhaTNumber=0;
475 m_psipKstarKNumber=0;
476 m_psipPPPiPiNumber=0;
478 log << MSG::INFO <<
"successfully return from initialize()" <<endmsg;
479 return StatusCode::SUCCESS;
488 MsgStream log(
msgSvc(), name());
489 log << MSG::INFO <<
"in execute()" << endreq;
491 if( m_writeDst) m_subalg1->execute();
492 if( m_writeRec) m_subalg2->execute();
495 m_isRadBhabha =
false;
498 m_isRadBhabhaT =
false;
501 m_isRecoJpsi =
false;
506 m_isGenProton =
false;
507 m_isPsipRhoPi =
false;
508 m_isPsipKstarK =
false;
509 m_isPsipPPPiPi =
false;
510 m_isPsippCand =
false;
512 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
515 cout<<
" eventHeader "<<endl;
516 return StatusCode::FAILURE;
519 int run=eventHeader->runNumber();
520 int event=eventHeader->eventNumber();
527 cout<<
"new run is:"<<m_run<<endl;
528 cout<<
"beamE is:"<<beamE<<endl;
529 if(beamE>0 && beamE<3)
535 if(m_display && m_events%1000==0){
536 cout<< m_events <<
" -------- run,event: "<<run<<
","<<
event<<endl;
537 cout<<
"m_ecm="<<m_ecm<<endl;
543 cout<<
" evtRecEvent "<<endl;
544 return StatusCode::FAILURE;
547 log << MSG::DEBUG <<
"ncharg, nneu, tottks = "
548 << evtRecEvent->totalCharged() <<
" , "
549 << evtRecEvent->totalNeutral() <<
" , "
550 << evtRecEvent->totalTracks() <<endreq;
553 cout<<
" evtRecTrkCol "<<endl;
554 return StatusCode::FAILURE;
557 if(evtRecEvent->totalTracks()!=evtRecTrkCol->size())
return StatusCode::SUCCESS;
561 SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(),
"/Event/EvtRec/EvtRecPi0Col");
563 log << MSG::FATAL <<
"Could not find EvtRecPi0Col" << endreq;
564 return StatusCode::FAILURE;
568 int nPi0 = recPi0Col->size();
569 EvtRecPi0Col::iterator itpi0 = recPi0Col->begin();
571 double mpi0 = (*itpi0)->unconMass();
572 if ( fabs(
mpi0 - 0.135 )> 0.015 )
589 vector<int> iCP, iCN;
597 int nCosmicCharge =0;
599 for(
int i = 0;i < evtRecEvent->totalCharged(); i++)
606 if(!(*itTrk)->isMdcKalTrackValid())
continue;
611 double vx0 = mdcTrk->
x();
612 double vy0 = mdcTrk->
y();
613 double vz0 = mdcTrk->
z();
614 double vr0 = mdcTrk->
r();
615 double theta0 = mdcTrk->
theta();
616 double p0 =
P(mdcTrk);
617 double pt0 = sqrt( pow(
Px(mdcTrk),2)+pow(
Py(mdcTrk),2));
633 Hep3Vector xorigin(0,0,0);
635 Gaudi::svcLocator()->service(
"VertexDbSvc", vtxsvc);
639 xorigin.setX(dbv[0]);
640 xorigin.setY(dbv[1]);
641 xorigin.setZ(dbv[2]);
646 HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]);
649 HepVector vecipa = helixip.
a();
650 double db=fabs(vecipa[0]);
656 if(fabs(dz) >= m_vz0cut)
continue;
657 if(db >= m_vr0cut)
continue;
660 if(p0>m_cosmic_lowp && p0<20){
661 iCosmicGood.push_back((*itTrk)->trackId());
662 nCosmicCharge += mdcTrk->
charge();
667 if(pt0 >= m_pt0HighCut)
continue;
668 if(pt0 <= m_pt0LowCut)
continue;
670 iGood.push_back((*itTrk)->trackId());
671 nCharge += mdcTrk->
charge();
673 iCP.push_back((*itTrk)->trackId());
674 else if(mdcTrk->
charge()<0)
675 iCN.push_back((*itTrk)->trackId());
679 int nGood = iGood.size();
680 int nCosmicGood = iCosmicGood.size();
685 for(
int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
688 if(!(*itTrk)->isEmcShowerValid())
continue;
690 double eraw = emcTrk->
energy();
692 double costh =
cos(emcTrk->
theta());
693 if(fabs(costh)<0.83 && eraw<0.025)
continue;
694 if(fabs(costh)>=0.83 && eraw<0.05)
continue;
695 if(time<0 || time>14)
continue;
697 iGam.push_back((*itTrk)->trackId());
699 int nGam = iGam.size();
724 Hep3Vector Pt_charge(0,0,0);
726 for(
int i = 0; i < nGood; i++) {
733 if((*itTrk)->isMdcKalTrackValid()) {
740 double phi=
Phi(mdcTrk);
743 HepLorentzVector
ptrk;
755 Hep3Vector ptemp(
Px(mdcTrk),
Py(mdcTrk),0);
759 if((*itTrk)->isEmcShowerValid()) {
772 else if( p3 < pmax && p3> psec){
790 ppip.push_back(
ptrk);
793 ppim.push_back(
ptrk);
799 if((*itTrk)->isEmcShowerValid()) {
802 double eraw = emcTrk->
energy();
803 double phiemc = emcTrk->
phi();
804 double the = emcTrk->
theta();
806 HepLorentzVector pemc;
807 pemc.setPx(eraw*
sin(the)*
cos(phiemc));
808 pemc.setPy(eraw*
sin(the)*
sin(phiemc));
809 pemc.setPz(eraw*
cos(the));
811 pemc = pemc.boost(-0.011,0,0);
823 if(pmax!=0) eopmax=eccmax/pmax;
824 if(psec!=0) eopsec=eccsec/psec;
843 Hep3Vector Pt_neu(0,0,0);
845 for(
int i = 0; i < nGam; i++) {
848 double eraw = emcTrk->
energy();
849 double phi = emcTrk->
phi();
850 double the = emcTrk->
theta();
851 HepLorentzVector
ptrk;
857 pGam.push_back(
ptrk);
862 Hep3Vector ptemp(eraw*
sin(the)*
cos(phi), eraw*
sin(the)*
sin(phi),0);
870 double esum = echarge + eneu;
871 Hep3Vector Pt=Pt_charge+Pt_neu;
874 double phid=phimax-phisec;
875 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+
pi, CLHEP::twopi) - CLHEP::pi;
882 if( nGood == 2 && nCharge==0 && (m_selectBhabha || m_selectDimu) ){
885 if( m_events%m_bhabha_scale == 0 && m_selectBhabha && echarge/m_ecm>0.8 && eopmax>0.8 && eopsec>0.8){
892 if( iCP.size()==1 && iCN.size()==1 && m_events%m_dimu_scale == 0 && m_selectDimu && eemc/m_ecm<0.3){
950 double time1=-99,depth1=-99;
951 double time2=-99,depth2=-99;
952 double p1=-99,
p2=-99;
953 double emc1=0, emc2=0;
954 if( (*itp)->isTofTrackValid() ){
955 SmartRefVector<RecTofTrack> tofTrkCol= (*itp)->tofTrack();
956 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
958 for(;iter_tof!=tofTrkCol.end();iter_tof++){
960 status->
setStatus( (*iter_tof)->status() );
962 time1=(*iter_tof)->tof();
968 if( (*itp)->isMucTrackValid() ){
970 depth1=mucTrk->
depth();
974 if((*itp)->isMdcKalTrackValid()) {
980 if((*itp)->isEmcShowerValid()) {
985 if( (*itn)->isTofTrackValid() ){
986 SmartRefVector<RecTofTrack> tofTrkCol= (*itn)->tofTrack();
987 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
989 for(;iter_tof!=tofTrkCol.end();iter_tof++){
991 status->
setStatus( (*iter_tof)->status() );
993 time2=(*iter_tof)->tof();
999 if( (*itn)->isMucTrackValid() ){
1001 depth2=mucTrk->
depth();
1004 if((*itn)->isMdcKalTrackValid()) {
1011 if((*itn)->isEmcShowerValid()) {
1016 double ebeam=m_ecm/2.0;
1018 && (emc1/p1<0.4 || emc2/
p2<0.4) && (depth1>3 || depth2>3) &&
1019 emc1>0.15 && emc1<0.3 && emc2>0.15 && emc2<0.3 && emc1/p1<0.8 && emc2/
p2<0.8 &&
1020 ((depth1>80*p1-60 || depth1>40) || (depth2>80*
p2-60 || depth2>40)) ){
1033 if(m_selectCosmic && nCosmicGood == 2 && eemc/m_ecm<0.3){
1038 double time1=-99,depth1=-99;
1039 double time2=-99,depth2=-99;
1040 if( (*itp)->isTofTrackValid() ){
1041 SmartRefVector<RecTofTrack> tofTrkCol= (*itp)->tofTrack();
1042 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
1044 for(;iter_tof!=tofTrkCol.end();iter_tof++){
1046 status->
setStatus( (*iter_tof)->status() );
1048 time1=(*iter_tof)->tof();
1054 if( (*itp)->isMucTrackValid() ){
1056 depth1=mucTrk->
depth();
1059 if( (*itn)->isTofTrackValid() ){
1060 SmartRefVector<RecTofTrack> tofTrkCol= (*itn)->tofTrack();
1061 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
1063 for(;iter_tof!=tofTrkCol.end();iter_tof++){
1065 status->
setStatus( (*iter_tof)->status() );
1067 time2=(*iter_tof)->tof();
1073 if( (*itn)->isMucTrackValid() ){
1075 depth2=mucTrk->
depth();
1080 if( m_selectCosmic && time1!=-99 && time2!=-99 && fabs(time1-time2)>5 ){
1091 if(m_events%m_proton_scale ==0 ){
1093 bool protoncand=
false;
1096 for(
int i=0; i<nGood; i++){
1102 double pp =
P(mdcTrk);
1106 if((*iter)->isMdcDedxValid()){
1109 chiP=dedxTrk->
chiP();
1112 if( fabs(pp)<0.6 && dedx>1.2 && fabs(chiP)<6){
1113 ncharge+=mdcTrk->
charge();
1119 if( (nproton==1 && ncharge==-1) || (nproton>=2 && ncharge<=nproton-2))
1123 m_genProtonNumber++;
1133 h_nGood->fill(nGood);
1134 h_nCharge->fill(nCharge);
1135 h_pmaxobeam->fill(pmax/(m_ecm/2.0));
1136 h_eopmax->fill(eopmax);
1137 h_eopsec->fill(eopsec);
1138 h_deltap->fill(pmax-psec);
1139 h_esumoecm->fill(esum/m_ecm);
1140 h_egmax->fill(egmax);
1141 h_deltaphi->fill(phid);
1142 h_Pt->fill(Pt.mag());
1148 if(nGood>1 && pmax/(m_ecm/2.0)>0.4 && eopmax>0.5 && esum/m_ecm>0.75 &&
1149 egmax>0.08 && fabs( fabs(phid)-CLHEP::pi)*180.0/CLHEP::pi>2.86 )
1152 m_radBhabhaNumber++;
1158 if(nGood==2 && nCharge==0 && eemc/m_ecm>=0.7 && eopmax>=0.85 && eopmax<=1.15 && eopsec>=0.85 && eopsec<=1.15 ){
1160 if( fabs(fabs(pmax)-m_ecm/2.0)<0.1 && fabs(fabs(psec)-m_ecm/2.0)<0.1 ){
1161 if(m_events%1000==0){
1162 m_isRadBhabhaT=
true;
1163 m_radBhabhaTNumber++;
1167 m_isRadBhabhaT=
true;
1168 m_radBhabhaTNumber++;
1179 if(!m_isRadBhabhaT && nGood==2 && nCharge==0 && eopmax>=0.85 && eopmax<=1.15 && eopsec>=0.85 && eopsec<=1.15 && eemc/m_ecm<=0.8 && Pt.mag()<=0.2)
1186 if(m_selectGG4Pi && nGood==4 && nCharge==0 && pmax/(m_ecm/2.0)>0.04 && pmax/(m_ecm/2.0)<0.9 && esum/m_ecm>0.04 && esum/m_ecm<0.75 && Pt.mag()<=0.2)
1193 if( (m_selectRhoPi || m_selectKstarK) && nGood == 2 && nCharge==0 && nPi0 == 1){
1199 if( (*itone)->isMdcKalTrackValid() && (*ittwo)->isMdcKalTrackValid() ) {
1204 HepLorentzVector p4trk1;
1206 p4trk1.setPy(
Py(trk1));
1207 p4trk1.setPz(
Pz(trk1));
1210 HepLorentzVector p4trk2;
1211 p4trk2.setPx(
Px(trk2));
1212 p4trk2.setPy(
Py(trk2));
1213 p4trk2.setPz(
Pz(trk2));
1217 HepLorentzVector p4trk3;
1218 p4trk3.setPx(
Px(trk1));
1219 p4trk3.setPy(
Py(trk1));
1220 p4trk3.setPz(
Pz(trk1));
1221 p4trk3.setE(sqrt( pow(
P(trk1),2)+
mpi*
mpi) );
1223 HepLorentzVector p4trk4;
1224 p4trk4.setPx(
Px(trk2));
1225 p4trk4.setPy(
Py(trk2));
1226 p4trk4.setPz(
Pz(trk2));
1227 p4trk4.setE(sqrt( pow(
P(trk2),2)+
mpi*
mpi) );
1231 itpi0 = recPi0Col->begin();
1232 HepLorentzVector p4gam1 = (*itpi0)->hiPfit();
1233 HepLorentzVector p4gam2 = (*itpi0)->loPfit();
1234 HepLorentzVector p4pi0 = p4gam1+p4gam2;
1236 HepLorentzVector p4total = p4trk1 + p4trk2 + p4pi0;
1237 HepLorentzVector p4total2 = p4trk3 + p4trk4 + p4pi0;
1239 double rhopimass=p4total2.m();
1240 double kstarkmass=p4total.m();
1241 if( (kstarkmass > 3.0 && kstarkmass < 3.15) || (rhopimass > 3.0 && rhopimass < 3.15) ){
1246 double eop1=0.0,eop2=0.0;
1247 if( (*itone)->isEmcShowerValid() ){
1249 double etrk=emcTrk->
energy();
1252 eop1 = etrk/
P(trk1);
1257 if( (*ittwo)->isEmcShowerValid() ){
1259 double etrk=emcTrk->
energy();
1262 eop2 = etrk/
P(trk2);
1267 if(eop1<0.8 && eop2<0.8){
1269 if(rhopimass>3.0 && rhopimass<3.15){
1272 HepLorentzVector
ecms(0.034,0,0,m_ecm);
1292 bool oksq1 = kmfit->
Fit();
1293 double chi1 = kmfit->
chisq();
1296 if(oksq1 && chi1<=60){
1303 if(kstarkmass>3.0 && kstarkmass<3.15){
1306 double mkstarp=0, mkstarm=0;
1308 HepLorentzVector p4kstarp = p4trk1 + p4pi0;
1309 HepLorentzVector p4kstarm = p4trk2 + p4pi0;
1310 mkstarp = p4kstarp.m();
1311 mkstarm = p4kstarm.m();
1314 HepLorentzVector p4kstarm = p4trk1 + p4pi0;
1315 HepLorentzVector p4kstarp = p4trk2 + p4pi0;
1316 mkstarp = p4kstarp.m();
1317 mkstarm = p4kstarm.m();
1320 if ( (fabs(mkstarp-0.89166)<0.1 && fabs(mkstarm-0.89166)>0.1) || (fabs(mkstarm-0.89166)<0.1 && fabs(mkstarp-0.89166)>0.1 ) ){
1323 HepLorentzVector
ecms(0.034,0,0,m_ecm);
1342 bool oksq1 = kmfit->
Fit();
1343 double chi1 = kmfit->
chisq();
1346 if(oksq1 && chi1<=60){
1367 if(m_selectPPPiPi && nGood ==4 && nCharge == 0 && nPi0==0){
1378 HepLorentzVector p4trkpp;
1379 HepLorentzVector p4trkpm;
1380 HepLorentzVector p4trkpip;
1381 HepLorentzVector p4trkpim;
1386 p4trkpp.setPx(
Px(trk1));
1387 p4trkpp.setPy(
Py(trk1));
1388 p4trkpp.setPz(
Pz(trk1));
1392 p4trkpm.setPx(
Px(trk3));
1393 p4trkpm.setPy(
Py(trk3));
1394 p4trkpm.setPz(
Pz(trk3));
1398 p4trkpip.setPx(
Px(trk2));
1399 p4trkpip.setPy(
Py(trk2));
1400 p4trkpip.setPz(
Pz(trk2));
1401 p4trkpip.setE(sqrt( pow(
P(trk2),2)+
mpi*
mpi) );
1404 p4trkpim.setPx(
Px(trk4));
1405 p4trkpim.setPy(
Py(trk4));
1406 p4trkpim.setPz(
Pz(trk4));
1407 p4trkpim.setE(sqrt( pow(
P(trk4),2)+
mpi*
mpi) );
1409 double jpsimass1= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1415 p4trkpp.setPx(
Px(trk1));
1416 p4trkpp.setPy(
Py(trk1));
1417 p4trkpp.setPz(
Pz(trk1));
1421 p4trkpm.setPx(
Px(trk4));
1422 p4trkpm.setPy(
Py(trk4));
1423 p4trkpm.setPz(
Pz(trk4));
1427 p4trkpip.setPx(
Px(trk2));
1428 p4trkpip.setPy(
Py(trk2));
1429 p4trkpip.setPz(
Pz(trk2));
1430 p4trkpip.setE(sqrt( pow(
P(trk2),2)+
mpi*
mpi) );
1433 p4trkpim.setPx(
Px(trk3));
1434 p4trkpim.setPy(
Py(trk3));
1435 p4trkpim.setPz(
Pz(trk3));
1436 p4trkpim.setE(sqrt( pow(
P(trk3),2)+
mpi*
mpi) );
1438 double jpsimass2= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1445 p4trkpp.setPx(
Px(trk2));
1446 p4trkpp.setPy(
Py(trk2));
1447 p4trkpp.setPz(
Pz(trk2));
1451 p4trkpm.setPx(
Px(trk3));
1452 p4trkpm.setPy(
Py(trk3));
1453 p4trkpm.setPz(
Pz(trk3));
1457 p4trkpip.setPx(
Px(trk1));
1458 p4trkpip.setPy(
Py(trk1));
1459 p4trkpip.setPz(
Pz(trk1));
1460 p4trkpip.setE(sqrt( pow(
P(trk1),2)+
mpi*
mpi) );
1463 p4trkpim.setPx(
Px(trk4));
1464 p4trkpim.setPy(
Py(trk4));
1465 p4trkpim.setPz(
Pz(trk4));
1466 p4trkpim.setE(sqrt( pow(
P(trk4),2)+
mpi*
mpi) );
1468 double jpsimass3= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1475 p4trkpp.setPx(
Px(trk2));
1476 p4trkpp.setPy(
Py(trk2));
1477 p4trkpp.setPz(
Pz(trk2));
1481 p4trkpm.setPx(
Px(trk4));
1482 p4trkpm.setPy(
Py(trk4));
1483 p4trkpm.setPz(
Pz(trk4));
1487 p4trkpip.setPx(
Px(trk1));
1488 p4trkpip.setPy(
Py(trk1));
1489 p4trkpip.setPz(
Pz(trk1));
1490 p4trkpip.setE(sqrt( pow(
P(trk1),2)+
mpi*
mpi) );
1493 p4trkpim.setPx(
Px(trk3));
1494 p4trkpim.setPy(
Py(trk3));
1495 p4trkpim.setPz(
Pz(trk3));
1496 p4trkpim.setE(sqrt( pow(
P(trk3),2)+
mpi*
mpi) );
1498 double jpsimass4= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1503 if( fabs(jpsimass1-3.075)<=0.075 || fabs(jpsimass2-3.075)<=0.075 ||
1504 fabs(jpsimass3-3.075)<=0.075 || fabs(jpsimass4-3.075)<=0.075){
1509 double chi1, chi2, chi3, chi4;
1512 HepLorentzVector
ecms(0.034,0,0,m_ecm);
1531 bool oksq1 = kmfit->
Fit();
1532 chi1 = kmfit->
chisq();
1556 bool oksq1 = kmfit->
Fit();
1557 chi2 = kmfit->
chisq();
1563 else if(chi2<chimin){
1586 bool oksq1 = kmfit->
Fit();
1587 chi3 = kmfit->
chisq();
1593 else if(chi3<chimin){
1618 bool oksq1 = kmfit->
Fit();
1619 chi4 = kmfit->
chisq();
1626 else if(chi4<chimin){
1635 if(type!=0 && chimin<100){
1649 if( (m_selectPsipRhoPi || m_selectPsipKstarK || m_selectPsipPPPiPi) && (nGood==4 || nGood==6) && nCharge==0 ){
1654 double bestmass=1.0;
1656 vector<int> iJPsiP,iJPsiN;
1657 for(
int ip=0; ip<iCP.size(); ip++){
1658 for(
int in=0; in<iCN.size();in++){
1659 pione = evtRecTrkCol->begin() + iCP[ip];
1660 pitwo = evtRecTrkCol->begin() + iCN[in];
1661 pitrk1 = (*pione)->mdcKalTrack();
1662 pitrk2 = (*pitwo)->mdcKalTrack();
1663 Hep3Vector p1(
Px(pitrk1),
Py(pitrk1),
Pz(pitrk1));
1664 Hep3Vector
p2(
Px(pitrk2),
Py(pitrk2),
Pz(pitrk2));
1665 double E1=sqrt( pow(
P(pitrk1),2)+
mpi*
mpi);
1666 double E2=sqrt( pow(
P(pitrk2),2)+
mpi*
mpi);
1667 double recomass = sqrt(pow(3.686-E1-E2,2)- (-(p1+
p2)).mag2() );
1669 if(fabs(recomass-3.096)< fabs(bestmass-3.096)){
1679 pione = evtRecTrkCol->begin() + iCP[pindex];
1680 pitwo = evtRecTrkCol->begin() + iCN[nindex];
1681 pitrk1 = (*pione)->mdcKalTrack();
1682 pitrk2 = (*pitwo)->mdcKalTrack();
1686 double jpsicharge=0;
1687 for(
int ip=0; ip<iCP.size(); ip++){
1689 iJPsiP.push_back(iCP[ip]);
1692 jpsicharge+=trktmp->
charge();
1696 for(
int in=0; in<iCN.size(); in++){
1698 iJPsiN.push_back(iCN[in]);
1701 jpsicharge+=trktmp->
charge();
1706 HepLorentzVector
ecms(0.034,0,0,m_ecm);
1709 if( (m_selectPsipRhoPi || m_selectPsipKstarK) && (bestmass>3.075 && bestmass<3.120) && nGood==4 && jpsicharge==0 && nPi0==1){
1715 if( (*itone)->isMdcKalTrackValid() && (*ittwo)->isMdcKalTrackValid() ) {
1720 HepLorentzVector p4trk1;
1722 p4trk1.setPy(
Py(trk1));
1723 p4trk1.setPz(
Pz(trk1));
1726 HepLorentzVector p4trk2;
1727 p4trk2.setPx(
Px(trk2));
1728 p4trk2.setPy(
Py(trk2));
1729 p4trk2.setPz(
Pz(trk2));
1733 HepLorentzVector p4trk3;
1734 p4trk3.setPx(
Px(trk1));
1735 p4trk3.setPy(
Py(trk1));
1736 p4trk3.setPz(
Pz(trk1));
1737 p4trk3.setE(sqrt( pow(
P(trk1),2)+
mpi*
mpi) );
1739 HepLorentzVector p4trk4;
1740 p4trk4.setPx(
Px(trk2));
1741 p4trk4.setPy(
Py(trk2));
1742 p4trk4.setPz(
Pz(trk2));
1743 p4trk4.setE(sqrt( pow(
P(trk2),2)+
mpi*
mpi) );
1746 EvtRecPi0Col::iterator itpi0 = recPi0Col->begin();
1755 HepLorentzVector p4gam1 = (*itpi0)->hiPfit();
1756 HepLorentzVector p4gam2 = (*itpi0)->loPfit();
1757 HepLorentzVector p4pi0 = p4gam1 + p4gam2;
1758 HepLorentzVector p4total = p4trk1 + p4trk2 + p4pi0;
1759 HepLorentzVector p4total2 = p4trk3 + p4trk4 + p4pi0;
1762 HepLorentzVector p4kstarp = p4trk1 + p4pi0;
1763 HepLorentzVector p4kstarm = p4trk2 + p4pi0;
1764 double mkstarp = p4kstarp.m();
1765 double mkstarm = p4kstarm.m();
1767 if( (p4total.m() > 3.0 && p4total.m() < 3.15) || (p4total2.m() > 3.0 && p4total2.m() < 3.15) ){
1788 if(m_selectPsipRhoPi){
1799 bool oksq1 = kmfit->
Fit();
1800 double chi1 = kmfit->
chisq();
1803 if(oksq1 && chi1>0 && chi1<100){
1804 m_isPsipRhoPi =
true;
1805 m_psipRhoPiNumber++;
1808 if(m_selectPsipKstarK){
1820 bool oksq2 = kmfit2->
Fit();
1821 double chi2 = kmfit2->
chisq();
1823 if(oksq2 && chi2>0 && chi2<200 &&
1824 ( (fabs(mkstarp-0.89166)<0.1 && fabs(mkstarm-0.89166)>0.1)
1825 || (fabs(mkstarm-0.89166)<0.1 && fabs(mkstarp-0.89166)>0.1 ) )){
1826 m_isPsipKstarK =
true;
1827 m_psipKstarKNumber++;
1842 if(m_selectPsipPPPiPi && (bestmass>3.075 && bestmass<3.120) && nGood==6 && jpsicharge==0 && nPi0==0){
1854 HepLorentzVector p4trkpp;
1855 HepLorentzVector p4trkpm;
1856 HepLorentzVector p4trkpip;
1857 HepLorentzVector p4trkpim;
1862 p4trkpp.setPx(
Px(trk1));
1863 p4trkpp.setPy(
Py(trk1));
1864 p4trkpp.setPz(
Pz(trk1));
1868 p4trkpm.setPx(
Px(trk3));
1869 p4trkpm.setPy(
Py(trk3));
1870 p4trkpm.setPz(
Pz(trk3));
1875 p4trkpip.setPx(
Px(trk2));
1876 p4trkpip.setPy(
Py(trk2));
1877 p4trkpip.setPz(
Pz(trk2));
1878 p4trkpip.setE(sqrt( pow(
P(trk2),2)+
mpi*
mpi) );
1881 p4trkpim.setPx(
Px(trk4));
1882 p4trkpim.setPy(
Py(trk4));
1883 p4trkpim.setPz(
Pz(trk4));
1884 p4trkpim.setE(sqrt( pow(
P(trk4),2)+
mpi*
mpi) );
1886 double jpsimass1= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1892 p4trkpp.setPx(
Px(trk1));
1893 p4trkpp.setPy(
Py(trk1));
1894 p4trkpp.setPz(
Pz(trk1));
1898 p4trkpm.setPx(
Px(trk4));
1899 p4trkpm.setPy(
Py(trk4));
1900 p4trkpm.setPz(
Pz(trk4));
1904 p4trkpip.setPx(
Px(trk2));
1905 p4trkpip.setPy(
Py(trk2));
1906 p4trkpip.setPz(
Pz(trk2));
1907 p4trkpip.setE(sqrt( pow(
P(trk2),2)+
mpi*
mpi) );
1910 p4trkpim.setPx(
Px(trk3));
1911 p4trkpim.setPy(
Py(trk3));
1912 p4trkpim.setPz(
Pz(trk3));
1913 p4trkpim.setE(sqrt( pow(
P(trk3),2)+
mpi*
mpi) );
1915 double jpsimass2= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1922 p4trkpp.setPx(
Px(trk2));
1923 p4trkpp.setPy(
Py(trk2));
1924 p4trkpp.setPz(
Pz(trk2));
1928 p4trkpm.setPx(
Px(trk3));
1929 p4trkpm.setPy(
Py(trk3));
1930 p4trkpm.setPz(
Pz(trk3));
1934 p4trkpip.setPx(
Px(trk1));
1935 p4trkpip.setPy(
Py(trk1));
1936 p4trkpip.setPz(
Pz(trk1));
1937 p4trkpip.setE(sqrt( pow(
P(trk1),2)+
mpi*
mpi) );
1940 p4trkpim.setPx(
Px(trk4));
1941 p4trkpim.setPy(
Py(trk4));
1942 p4trkpim.setPz(
Pz(trk4));
1943 p4trkpim.setE(sqrt( pow(
P(trk4),2)+
mpi*
mpi) );
1945 double jpsimass3= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1952 p4trkpp.setPx(
Px(trk2));
1953 p4trkpp.setPy(
Py(trk2));
1954 p4trkpp.setPz(
Pz(trk2));
1958 p4trkpm.setPx(
Px(trk4));
1959 p4trkpm.setPy(
Py(trk4));
1960 p4trkpm.setPz(
Pz(trk4));
1964 p4trkpip.setPx(
Px(trk1));
1965 p4trkpip.setPy(
Py(trk1));
1966 p4trkpip.setPz(
Pz(trk1));
1967 p4trkpip.setE(sqrt( pow(
P(trk1),2)+
mpi*
mpi) );
1970 p4trkpim.setPx(
Px(trk3));
1971 p4trkpim.setPy(
Py(trk3));
1972 p4trkpim.setPz(
Pz(trk3));
1973 p4trkpim.setE(sqrt( pow(
P(trk3),2)+
mpi*
mpi) );
1975 double jpsimass4= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1977 if( fabs(jpsimass1-3.075)<=0.075 || fabs(jpsimass2-3.075)<=0.075 ||
1978 fabs(jpsimass3-3.075)<=0.075 || fabs(jpsimass4-3.075)<=0.075){
1982 double chi1, chi2, chi3, chi4;
1987 WTrackParameter wvpiTrk1, wvpiTrk2, wvpiTrk3, wvpiTrk4 , wvpiTrk5, wvpiTrk6 ;
2009 bool oksq1 = kmfit->
Fit();
2010 chi1 = kmfit->
chisq();
2039 bool oksq1 = kmfit->
Fit();
2040 chi2 = kmfit->
chisq();
2046 else if(chi2<chimin){
2075 bool oksq1 = kmfit->
Fit();
2076 chi3 = kmfit->
chisq();
2082 else if(chi3<chimin){
2109 bool oksq1 = kmfit->
Fit();
2110 chi4 = kmfit->
chisq();
2116 else if(chi4<chimin){
2125 if(chimin>0 && chimin<200){
2126 m_isPsipPPPiPi =
true;
2127 m_psipPPPiPiNumber++;
2143 if(m_selectPsippCand){
2146 if ( ! evtRecDTagCol ) {
2147 log << MSG::FATAL <<
"Could not find EvtRecDTagCol" << endreq;
2148 return StatusCode::FAILURE;
2150 if(evtRecDTagCol->size()>0){
2152 EvtRecDTagCol::iterator m_iterbegin=evtRecDTagCol->begin();
2153 EvtRecDTagCol::iterator m_iterend=evtRecDTagCol->end();
2154 for(EvtRecDTagCol::iterator
iter=m_iterbegin;
iter != m_iterend;
iter++){
2156 if( ((*iter)->decayMode()==
EvtRecDTag::kD0toKPiPi0 && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03) ||
2157 ((*iter)->decayMode()==
EvtRecDTag::kD0toKPiPi0Pi0 && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03) ||
2158 ((*iter)->decayMode()==
EvtRecDTag::kD0toKPiPiPi && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.03 && (*iter)->deltaE()<0.03) ||
2159 ((*iter)->decayMode()==
EvtRecDTag::kD0toKPiPiPiPi0 && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03) ||
2160 ((*iter)->decayMode()==
EvtRecDTag::kD0toKsPiPi && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.03 && (*iter)->deltaE()<0.03) ||
2161 ((*iter)->decayMode()==
EvtRecDTag::kD0toKsPiPiPi0 && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03)||
2162 ((*iter)->decayMode()==
EvtRecDTag::kDptoKPiPi && fabs((*iter)->mBC()-1.87)<0.006 && (*iter)->deltaE()>-0.03 && (*iter)->deltaE()<0.03) ||
2163 ((*iter)->decayMode()==
EvtRecDTag::kDptoKPiPiPi0 && fabs((*iter)->mBC()-1.87)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03)||
2164 ((*iter)->decayMode()==
EvtRecDTag::kDptoKsPiPi0 && fabs((*iter)->mBC()-1.87)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03) ||
2165 ((*iter)->decayMode()==
EvtRecDTag::kDptoKsPiPiPi&& fabs((*iter)->mBC()-1.87)<0.006 && (*iter)->deltaE()>-0.03 && (*iter)->deltaE()<0.03) ) {
2166 m_isPsippCand =
true;
2167 m_psippCandNumber++;
2181 if( m_selectRadBhabha && m_isRadBhabha ) m_subalg3->execute();
2182 if( m_selectGGEE && m_isGGEE ) m_subalg4->execute();
2183 if( m_selectGG4Pi && m_isGG4Pi ) m_subalg5->execute();
2184 if( m_selectRadBhabhaT && m_isRadBhabhaT ) m_subalg6->execute();
2185 if( m_selectRhoPi && m_isRhoPi ) m_subalg7->execute();
2186 if( m_selectKstarK && m_isKstarK ) m_subalg8->execute();
2187 if( m_selectPPPiPi && m_isPPPiPi ) m_subalg9->execute();
2188 if( m_selectRecoJpsi && m_isRecoJpsi ) m_subalg10->execute();
2189 if( m_selectBhabha && m_isBhabha ) m_subalg11->execute();
2190 if( m_selectDimu && m_isDimu ) m_subalg12->execute();
2191 if( m_selectCosmic && m_isCosmic ) m_subalg13->execute();
2192 if( m_selectGenProton && m_isGenProton ) m_subalg14->execute();
2193 if( m_selectPsipRhoPi && m_isPsipRhoPi ) m_subalg15->execute();
2194 if( m_selectPsipKstarK && m_isPsipKstarK ) m_subalg16->execute();
2195 if( m_selectPsipPPPiPi && m_isPsipPPPiPi ) m_subalg17->execute();
2196 if( m_selectPsippCand && m_isPsippCand ) m_subalg18->execute();
2220 return StatusCode::SUCCESS;
2227 MsgStream log(
msgSvc(), name());
2228 log << MSG::INFO <<
"in finalize()" << endmsg;
2230 cout<<
"Total events: "<<m_events<<endl;
2233 if(m_selectRadBhabha) {
2234 cout <<
"Selected rad bhabha: " << m_radBhabhaNumber << endl;
2239 cout <<
"Selected ggee events: " << m_GGEENumber << endl;
2243 cout <<
"Selected gg4pi events: " << m_GG4PiNumber << endl;
2246 if(m_selectRadBhabhaT) {
2247 cout <<
"Selected rad bhabha tight: " << m_radBhabhaTNumber << endl;
2251 cout <<
"Selected rhopi events: " << m_rhoPiNumber << endl;
2254 if(m_selectKstarK) {
2255 cout <<
"Selected kstark events: " << m_kstarKNumber << endl;
2258 if(m_selectPPPiPi) {
2259 cout <<
"Selected pppipi events: " << m_ppPiPiNumber << endl;
2262 if(m_selectRecoJpsi) {
2263 cout <<
"Selected recoil jsi events: " << m_recoJpsiNumber << endl;
2267 if(m_selectBhabha) {
2268 cout <<
"Selected bhabha events: " << m_bhabhaNumber << endl;
2273 cout <<
"Selected dimu events: " << m_dimuNumber << endl;
2277 if(m_selectCosmic) {
2278 cout <<
"Selected cosmic events: " << m_cosmicNumber << endl;
2281 if(m_selectGenProton) {
2282 cout <<
"Selected generic proton events: " << m_genProtonNumber << endl;
2285 if(m_selectPsipRhoPi) {
2286 cout <<
"Selected recoil rhopi events: " << m_psipRhoPiNumber << endl;
2289 if(m_selectPsipKstarK) {
2290 cout <<
"Selected recoil kstark events: " << m_psipKstarKNumber << endl;
2293 if(m_selectPsipPPPiPi) {
2294 cout <<
"Selected recoil pppipi events: " << m_psipPPPiPiNumber << endl;
2297 if(m_selectPsippCand) {
2298 cout <<
"Selected psi'' candi events: " << m_psippCandNumber << endl;
2301 return StatusCode::SUCCESS;
2307 double delta=0.0610865;
2312 if(
phi2>CLHEP::twopi)
phi2 -= CLHEP::twopi;
2323 if(ph<=phi2&&ph>=
phi1)
return true;
2327 if(ph>=
phi2||ph<=
phi1)
return true;
2336 const char host[] =
"202.122.37.69";
2337 const char user[] =
"guest";
2338 const char passwd[] =
"guestpass";
2339 const char db[] =
"run";
2340 unsigned int port_num = 3306;
2344 mysql = mysql_real_connect(mysql, host, user, passwd, db, port_num,
2348 if (mysql ==
NULL) {
2349 fprintf(stderr,
"can not open database: RunInfo for run %d lum\n",
runNo);
2354 snprintf(stmt, 1024,
2355 "select BER_PRB, BPR_PRB "
2356 "from RunParams where run_number = %d",
runNo);
2357 if (mysql_real_query(mysql, stmt, strlen(stmt))) {
2358 fprintf(stderr,
"query error\n");
2361 MYSQL_RES* result_set = mysql_store_result(mysql);
2362 MYSQL_ROW row = mysql_fetch_row(result_set);
2364 fprintf(stderr,
"cannot find data for RunNo %d\n",
runNo);
2368 double E_E=0, E_P=0;
2369 sscanf(row[0],
"%lf", &E_E);
2370 sscanf(row[1],
"%lf", &E_P);
2372 beamE=(E_E+E_P)/2.0;
double sin(const BesAngle a)
double cos(const BesAngle a)
double Px(RecMdcKalTrack *trk)
HepGeom::Point3D< double > HepPoint3D
double Pz(RecMdcKalTrack *trk)
std::vector< HepLorentzVector > Vp4
double Phi(RecMdcKalTrack *trk)
double P(RecMdcKalTrack *trk)
double Py(RecMdcKalTrack *trk)
EvtRecTrackCol::iterator EvtRecTrackIterator
std::vector< HepLorentzVector > Vp4
struct st_mysql_res MYSQL_RES
IHistogramSvc * histoSvc()
bool WhetherSector(double ph, double ph1, double ph2)
void readbeamEfromDb(int runNo, double &beamE)
const double theta() const
void setPx(const double px, const int pid)
const double tanl(void) const
static void setPidType(PidType pidType)
const double kappa(void) const
const double fi0(void) const
virtual bool isVertexValid()=0
virtual double * SigmaPrimaryVertex()=0
virtual double * PrimaryVertex()=0
static KinematicFit * instance()
void AddFourMomentum(int number, HepLorentzVector p4)
const HepVector & getZHelix() const
const HepSymMatrix & getZError() const
HepSymMatrix & getZErrorK()
void setStatus(unsigned int status)
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
const HepPoint3D & pivot(void) const
returns pivot position.
const HepVector & a(void) const
returns helix parameters.
_EXTERN_ std::string EvtRecEvent
_EXTERN_ std::string EvtRecDTagCol
_EXTERN_ std::string EvtRecTrackCol