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"
6#include "GaudiKernel/Service.h"
7#include "GaudiKernel/ISvcLocator.h"
8#include "GaudiKernel/SvcFactory.h"
20#include "GaudiKernel/INTupleSvc.h"
21#include "GaudiKernel/NTuple.h"
22#include "GaudiKernel/Bootstrap.h"
23#include "GaudiKernel/IHistogramSvc.h"
35 m_iterstag(0), m_iterdtag1(0), m_iterdtag2(0),
36 m_chargebegin(0),m_chargeend(0),m_showerbegin(0),m_showerend(0),
37 m_tag1desigma(1.0),m_tag2desigma(1.0){
39 SmartDataPtr<EvtRecEvent> evtRecEvent(
eventSvc(),
"/Event/EvtRec/EvtRecEvent");
40 if ( ! evtRecEvent ) {
41 cout << MSG::FATAL <<
"Could not find EvtRecEvent" << endl;
45 SmartDataPtr<EvtRecTrackCol> evtRecTrackCol(
eventSvc(),
"/Event/EvtRec/EvtRecTrackCol");
46 if ( ! evtRecTrackCol ) {
47 cout << MSG::FATAL <<
"Could not find EvtRecTrackCol" << endl;
52 StatusCode sc = Gaudi::svcLocator()->service(
"SimplePIDSvc", m_simplePIDSvc);
61 m_chargebegin=evtRecTrackCol->begin();
62 m_chargeend=evtRecTrackCol->begin()+evtRecEvent->totalCharged();
63 m_showerbegin=evtRecTrackCol->begin()+evtRecEvent->totalCharged();
64 m_showerend=evtRecTrackCol->begin()+evtRecEvent->totalTracks();
68 if ( ! evtRecDTagCol ) {
69 cout <<
"Could not find EvtRecDTagCol" << endl;
74 SmartDataPtr<EvtRecVeeVertexCol> evtRecVeeVertexCol(
eventSvc(),
"/Event/EvtRec/EvtRecVeeVertexCol");
75 if ( ! evtRecVeeVertexCol ) {
76 cout<<
"Could not find EvtRecVeeVertexCol" << endl;
81 SmartDataPtr<EvtRecPi0Col> recPi0Col(
eventSvc(),
"/Event/EvtRec/EvtRecPi0Col");
83 cout<<
"Could not find EvtRecPi0Col" << endl;
89 SmartDataPtr<EvtRecEtaToGGCol> recEtaToGGCol(
eventSvc(),
"/Event/EvtRec/EvtRecEtaToGGCol");
90 if ( ! recEtaToGGCol ) {
91 cout<<
"Could not find EvtRecEtaToGGCol" << endl;
97 m_iterbegin=evtRecDTagCol->begin();
98 m_iterend=evtRecDTagCol->end();
99 m_pi0iterbegin=recPi0Col->begin();
100 m_pi0iterend=recPi0Col->end();
101 m_etaiterbegin=recEtaToGGCol->begin();
102 m_etaiterend=recEtaToGGCol->end();
103 m_ksiterbegin=evtRecVeeVertexCol->begin();
104 m_ksiterend=evtRecVeeVertexCol->end();
106 if(evtRecDTagCol->size() > 0)
107 m_isdtaglistempty=
false;
109 m_isdtaglistempty=
true;
117 for(EvtRecDTagCol::iterator
iter=m_iterbegin;
iter != m_iterend;
iter++){
119 if( (
int)( (*iter)->decayMode())< 200 )
120 m_d0modes.push_back( index );
121 else if( (
int)( (*iter)->decayMode())< 400 )
122 m_dpmodes.push_back( index );
123 else if( (
int)( (*iter)->decayMode())< 1000 )
124 m_dsmodes.push_back( index );
147bool DTagTool::compare(EvtRecDTagCol::iterator pair1_iter1,EvtRecDTagCol::iterator pair1_iter2,EvtRecDTagCol::iterator pair2_iter1,EvtRecDTagCol::iterator pair2_iter2,
double mD,
string smass){
152 value1 = fabs(0.5*((*pair1_iter1)->mBC()+(*pair1_iter2)->mBC())-mD);
153 value2 = fabs(0.5*((*pair2_iter1)->mBC()+(*pair2_iter2)->mBC())-mD);
155 else if(smass==
"de"){
156 value1 = pow((*pair1_iter1)->deltaE()/m_tag1desigma,2)+pow((*pair1_iter2)->deltaE()/m_tag2desigma,2);
157 value2 = pow((*pair2_iter1)->deltaE()/m_tag1desigma,2)+pow((*pair2_iter2)->deltaE()/m_tag2desigma,2);
159 else if(smass==
"inv"){
160 value1 = fabs(0.5*((*pair1_iter1)->mass()+(*pair1_iter2)->mass())-mD);
161 value2 = fabs(0.5*((*pair2_iter1)->mass()+(*pair2_iter2)->mass())-mD);
165 if( value1 <= value2)
175 for(EvtRecDTagCol::iterator
iter=m_iterbegin;
iter != m_iterend;
iter++){
178 if( (*iter)->decayMode() == decaymode && (*iter)->type() == 1 )
179 mode.push_back( index );
182 if( (*iter)->decayMode() == decaymode )
183 mode.push_back( index );
197 for(EvtRecDTagCol::iterator
iter=m_iterbegin;
iter != m_iterend;
iter++){
200 if( (*iter)->decayMode() == decaymode && (*iter)->type() == 1 )
201 mode.push_back( index );
204 if( (*iter)->decayMode() == decaymode )
205 mode.push_back( index );
240 EvtRecDTagCol::iterator iter_dtag=m_iterbegin;
241 for ( ; iter_dtag != m_iterend; iter_dtag++){
244 if( (*iter_dtag)->type()!=1 || (*iter_dtag)->decayMode() !=
mode || (*iter_dtag)->charm() != tagcharm )
248 if( (*iter_dtag)->decayMode() !=
mode || (*iter_dtag)->charm() != tagcharm )
251 if(fabs((*iter_dtag)->deltaE())<fabs(de_min)){
253 m_iterstag=iter_dtag;
254 de_min=(*iter_dtag)->deltaE();
270 EvtRecDTagCol::iterator iter_dtag=m_iterbegin;
271 for ( ; iter_dtag != m_iterend; iter_dtag++){
274 if( (*iter_dtag)->type()!=1 || (*iter_dtag)->decayMode() !=
mode )
278 if( (*iter_dtag)->decayMode() !=
mode )
282 if(fabs((*iter_dtag)->deltaE())<fabs(de_min)){
284 m_iterstag=iter_dtag;
285 de_min=(*iter_dtag)->deltaE();
303 return findDTag(
static_cast<int>(mode1),
static_cast<int>(mode2), smass);
309 return findDTag(
static_cast<int>(mode1), tagcharm1,
static_cast<int>(mode2), tagcharm2, smass);
317 return findADTag(
static_cast<int>(mode1),
static_cast<int>(mode2));
323 return findADTag(
static_cast<int>(mode1), tagcharm1,
static_cast<int>(mode2), tagcharm2);
333 int tagcharm1= (mode1<10 || mode1>=200)?+1:0;
334 int tagcharm2= (mode2<10 || mode2>=200)?-1:0;
336 if(tagcharm1*tagcharm2>0){
337 cout<<
"double tagging warning! two modes can't have same nonzero charmness"<<endl;
343 if( mode1 < 200 && mode2 < 200)
345 else if ( mode1>=200 && mode1 < 400 && mode2>=200 && mode2 < 400)
347 else if ( mode1>=400 && mode1 < 1000 && mode2>=400 && mode2 < 1000)
350 cout<<
"double tag modes are not from same particle ! "<<endl;
355 vector<int> igood1, igood2;
356 igood1.clear(),igood2.clear();
358 EvtRecDTagCol::iterator iter_dtag=m_iterbegin;
361 for ( ; iter_dtag != m_iterend; iter_dtag++){
362 int iter_mode=(*iter_dtag)->decayMode();
363 int iter_charm=(*iter_dtag)->charm();
364 int iter_type=(*iter_dtag)->type();
367 if( iter_mode == mode1 && iter_charm == tagcharm1 && iter_type==1 )
368 igood1.push_back(index);
369 if( tagcharm1!=0 && iter_mode == mode1 && iter_charm == -tagcharm1 && iter_type==1 )
370 igood1.push_back(index);
372 if( iter_mode == mode2 && iter_charm == tagcharm2 && iter_type==1 )
373 igood2.push_back(index);
374 if( tagcharm2!=0 && iter_mode == mode2 && iter_charm == -tagcharm2 && iter_type==1 )
375 igood2.push_back(index);
378 if( iter_mode == mode1 && iter_charm == tagcharm1 )
379 igood1.push_back(index);
380 if( tagcharm1!=0 && iter_mode == mode1 && iter_charm == -tagcharm1 )
381 igood1.push_back(index);
383 if( iter_mode == mode2 && iter_charm == tagcharm2 )
384 igood2.push_back(index);
385 if( tagcharm2!=0 && iter_mode == mode2 && iter_charm == -tagcharm2 )
386 igood2.push_back(index);
393 if(igood1.size()<1 || igood2.size()<1)
397 int index_i=0, index_j=0;
399 EvtRecDTagCol::iterator iter_i, iter_j;
401 for(
int i=0; i<igood1.size(); i++){
403 iter_i=m_iterbegin+igood1[i];
405 int charm_i=(*iter_i)->charm();
406 for(
int j=0;j<igood2.size();j++){
407 iter_j=m_iterbegin+igood2[j];
409 int charm_j=(*iter_j)->charm();
410 if( charm_i*charm_j>0 || igood2[j] == igood1[i] )
continue;
419 if(
compare(iter_i,iter_j,m_iterdtag1,m_iterdtag2,mDcand,smass) ){
435 if(tagcharm1*tagcharm2>0){
436 cout<<
"double tagging warning! two modes can't have same nonzero charmness"<<endl;
442 if( mode1 < 200 && mode2 < 200)
444 else if (mode1>=200 && mode1 < 400 && mode2>=200 && mode2 < 400)
446 else if (mode1>=400 && mode1 < 1000 && mode2>=400 && mode2 < 1000)
449 cout<<
"double tag modes are not from same particle ! "<<endl;
454 vector<int> igood1, igood2;
455 igood1.clear(),igood2.clear();
457 EvtRecDTagCol::iterator iter_dtag=m_iterbegin;
459 for ( ; iter_dtag != m_iterend; iter_dtag++){
460 int iter_mode=(*iter_dtag)->decayMode();
461 int iter_charm=(*iter_dtag)->charm();
462 int iter_type=(*iter_dtag)->type();
465 if( iter_mode == mode1 && iter_charm == tagcharm1 && iter_type==1 )
466 igood1.push_back(index);
468 if( iter_mode == mode2 && iter_charm == tagcharm2 && iter_type==1)
469 igood2.push_back(index);
472 if( iter_mode == mode1 && iter_charm == tagcharm1)
473 igood1.push_back(index);
475 if( iter_mode == mode2 && iter_charm == tagcharm2)
476 igood2.push_back(index);
485 if(igood1.size()<1 || igood2.size()<1)
489 int index_i=0, index_j=0;
494 EvtRecDTagCol::iterator iter_i, iter_j;
496 for(
int i=0; i<igood1.size(); i++){
498 iter_i=m_iterbegin+igood1[i];
499 int charm_i=(*iter_i)->charm();
500 for(
int j=0;j<igood2.size();j++){
501 iter_j=m_iterbegin+igood2[j];
502 int charm_j=(*iter_j)->charm();
503 if( charm_i*charm_j>0 || igood2[j] == igood1[i] )
continue;
512 if(
compare(iter_i,iter_j,m_iterdtag1,m_iterdtag2,mDcand,smass) ){
530 int tagcharm1= (mode1<10 || mode1>=200)?+1:0;
531 int tagcharm2= (mode2<10 || mode2>=200)?-1:0;
533 if(tagcharm1*tagcharm2>0){
534 cout<<
"double tagging warning! two modes can't have same nonzero charmness"<<endl;
540 if( mode1 < 200 && mode2 < 200)
542 else if ( mode1>=200 && mode1 < 400 && mode2>=200 && mode2 < 400)
544 else if ( mode1>=400 && mode1 < 1000 && mode2>=400 && mode2 < 1000)
547 cout<<
"double tag modes are not from same particle ! "<<endl;
552 vector<int> igood1, igood2;
553 igood1.clear(),igood2.clear();
555 EvtRecDTagCol::iterator iter_dtag=m_iterbegin;
558 for ( ; iter_dtag != m_iterend; iter_dtag++){
559 int iter_mode=(*iter_dtag)->decayMode();
560 int iter_charm=(*iter_dtag)->charm();
561 int iter_type=(*iter_dtag)->type();
564 if( iter_mode == mode1 && iter_charm == tagcharm1 && iter_type==1 )
565 igood1.push_back(index);
566 if( tagcharm1!=0 && iter_mode == mode1 && iter_charm == -tagcharm1 && iter_type==1 )
567 igood1.push_back(index);
569 if( iter_mode == mode2 && iter_charm == tagcharm2 && iter_type==1 )
570 igood2.push_back(index);
571 if( tagcharm2!=0 && iter_mode == mode2 && iter_charm == -tagcharm2 && iter_type==1 )
572 igood2.push_back(index);
575 if( iter_mode == mode1 && iter_charm == tagcharm1 )
576 igood1.push_back(index);
577 if( tagcharm1!=0 && iter_mode == mode1 && iter_charm == -tagcharm1 )
578 igood1.push_back(index);
580 if( iter_mode == mode2 && iter_charm == tagcharm2 )
581 igood2.push_back(index);
582 if( tagcharm2!=0 && iter_mode == mode2 && iter_charm == -tagcharm2 )
583 igood2.push_back(index);
592 EvtRecDTagCol::iterator iter_i, iter_j;
594 for(
int i=0; i<igood1.size(); i++){
596 iter_i=m_iterbegin+igood1[i];
597 int charm_i=(*iter_i)->charm();
599 for(
int j=0;j<igood2.size();j++){
600 iter_j=m_iterbegin+igood2[j];
601 int charm_j=(*iter_j)->charm();
602 if( charm_i*charm_j>0 || igood2[j] == igood1[i] )
continue;
605 m_viterdtag1.push_back(m_iterbegin+igood1[i]);
606 m_viterdtag2.push_back(m_iterbegin+igood2[j]);
612 if(m_viterdtag1.size()>0){
621 if(tagcharm1*tagcharm2>0){
622 cout<<
"double tagging warning! two modes can't have same nonzero charmness"<<endl;
628 if( mode1 < 200 && mode2 < 200)
630 else if (mode1>=200 && mode1 < 400 && mode2>=200 && mode2 < 400)
632 else if (mode1>=400 && mode1 < 1000 && mode2>=400 && mode2 < 1000)
635 cout<<
"double tag modes are not from same particle ! "<<endl;
640 vector<int> igood1, igood2;
641 igood1.clear(),igood2.clear();
643 EvtRecDTagCol::iterator iter_dtag=m_iterbegin;
645 for ( ; iter_dtag != m_iterend; iter_dtag++){
646 int iter_mode=(*iter_dtag)->decayMode();
647 int iter_charm=(*iter_dtag)->charm();
648 int iter_type=(*iter_dtag)->type();
651 if( iter_mode == mode1 && iter_charm == tagcharm1 && iter_type==1 )
652 igood1.push_back(index);
654 if( iter_mode == mode2 && iter_charm == tagcharm2 && iter_type==1)
655 igood2.push_back(index);
658 if( iter_mode == mode1 && iter_charm == tagcharm1)
659 igood1.push_back(index);
661 if( iter_mode == mode2 && iter_charm == tagcharm2)
662 igood2.push_back(index);
673 int index_i=0, index_j=0;
674 EvtRecDTagCol::iterator iter_i, iter_j;
676 for(
int i=0; i<igood1.size(); i++){
678 iter_i=m_iterbegin+igood1[i];
679 int charm_i=(*iter_i)->charm();
680 for(
int j=0;j<igood2.size();j++){
681 iter_j=m_iterbegin+igood2[j];
682 int charm_j=(*iter_j)->charm();
683 if( charm_i*charm_j>0 || igood2[j] == igood1[i] )
continue;
687 m_viterdtag1.push_back(m_iterbegin+igood1[i]);
688 m_viterdtag2.push_back(m_iterbegin+igood2[j]);
693 if(m_viterdtag1.size()>0)
710 cout<<
" print mode:"<< (*iter)->decayMode()<<endl;
711 cout<<
"beam energy is:"<< (*iter)->beamE()<<endl;
712 cout<<
"mBC is:"<< (*iter)->mBC()<<endl;
713 cout<<
"deltaE is:"<< (*iter)->deltaE()<<endl;
714 cout<<
"inv mass is:"<< (*iter)->mass()<<endl;
715 cout<<
"charge is:"<< (*iter)->charge()<<endl;
716 cout<<
"charm is:"<< (*iter)->charm()<<endl;
717 cout<<
"num of children is:"<< (*iter)->numOfChildren()<<endl;
719 cout<<
"found "<< (*iter)->tracks().size()<<
" same side tracks."<<endl;
720 cout<<
"found "<< (*iter)->otherTracks().size()<<
" other side tracks."<<endl;
721 cout<<
"found "<< (*iter)->showers().size()<<
" same side showers."<<endl;
722 cout<<
"found "<< (*iter)->otherShowers().size()<<
" other side showers."<<endl;
730 HepLorentzVector p41= (*pi0Itr)->hiPfit();
731 HepLorentzVector p42= (*pi0Itr)->loPfit();
741 HepLorentzVector p41=
p4shower(emctrk1);
742 HepLorentzVector p42=
p4shower(emctrk2);
752 HepLorentzVector p41= (*etaItr)->hiPfit();
753 HepLorentzVector p42= (*etaItr)->loPfit();
763 HepLorentzVector p41=
p4shower(emctrk1);
764 HepLorentzVector p42=
p4shower(emctrk2);
774 SmartRefVector<EvtRecTrack> showers=(*iter_dtag)->showers();
775 if(showers.size()<2*numpi0){
776 cout<<
"too many pi0 required"<<endl;
783 for(EvtRecPi0Col::iterator pi0Itr = m_pi0iterbegin;
784 pi0Itr < m_pi0iterend; pi0Itr++){
790 int heGammaTrkId = heGammaTrk->
trackId();
791 int leGammaTrkId = leGammaTrk->
trackId();
793 for(
int index=0; index<numpi0; ++index){
797 for(
int i=index*2; i<index*2+2; ++i){
799 if(!isheid && heGammaTrkId == showers[i]->trackId())
801 if(!isleid && leGammaTrkId == showers[i]->trackId())
806 canid.push_back( pi0Itr - m_pi0iterbegin);
810 if(canid.size()==numpi0){
825 SmartRefVector<EvtRecTrack> showers=(*iter_dtag)->showers();
826 if(showers.size()<2*numeta){
827 cout<<
"too many eta required"<<endl;
834 for(EvtRecEtaToGGCol::iterator etaItr = m_etaiterbegin;
835 etaItr < m_etaiterend; etaItr++){
841 int heGammaTrkId = heGammaTrk->
trackId();
842 int leGammaTrkId = leGammaTrk->
trackId();
844 for(
int index=0; index<numeta; ++index){
848 for(
int i=index*2; i<index*2+2; ++i){
850 if(!isheid && heGammaTrkId == showers[i]->trackId())
852 if(!isleid && leGammaTrkId == showers[i]->trackId())
857 canid.push_back( etaItr - m_etaiterbegin);
861 if(canid.size()==numeta){
874 SmartRefVector<EvtRecTrack> tracks=(*iter_dtag)->tracks();
875 if(tracks.size()<2*numks){
876 cout<<
"too many kshort required"<<endl;
886 for(EvtRecVeeVertexCol::iterator ksItr = m_ksiterbegin;
887 ksItr < m_ksiterend; ksItr++){
890 if((*ksItr)->vertexId() != 310)
continue;
895 int ksChild1TrkId = aKsChild1Trk->
trackId();
896 int ksChild2TrkId = aKsChild2Trk->
trackId();
898 for(
int index=0; index<numks; ++index){
902 for(
int i=index*2; i<index*2+2; ++i){
903 if(!isc1id && ksChild1TrkId == tracks[i]->trackId())
905 if(!isc2id && ksChild2TrkId == tracks[i]->trackId())
910 canid.push_back( ksItr - m_ksiterbegin);
913 if(canid.size()==numks){
924 double Eemc=shower->
energy();
925 double phi=shower->
phi();
926 double theta=shower->
theta();
927 HepLorentzVector
p4(Eemc*
sin(theta)*
cos(phi),Eemc*
sin(theta)*
sin(phi),Eemc*
cos(theta),Eemc);
957 double dr(0),phi0(0),kappa(0),dz(0),tanl(0);
966 if (kappa > 0.0000000001)
968 else if (kappa < -0.0000000001)
972 if(kappa!=0) pxy = 1.0/fabs(kappa);
974 double px = pxy * (-
sin(phi0));
975 double py = pxy *
cos(phi0);
976 double pz = pxy * tanl;
978 double e = sqrt( pxy*pxy + pz*pz +
mass*
mass );
980 return HepLorentzVector(px, py, pz, e);
987 SmartRefVector<EvtRecTrack> pionid=(*m_iterbegin)->pionId();
989 for(
int i=0; i < pionid.size() ;i++){
990 if( trk->
trackId() == pionid[i]->trackId()){
1001 SmartRefVector<EvtRecTrack> kaonid=(*m_iterbegin)->kaonId();
1003 for(
int i=0; i < kaonid.size() ;i++){
1004 if( trk->
trackId() == kaonid[i]->trackId()){
1059 dedxchi=dedxTrk->
chiMu();
1064 ptrk= mdcKalTrk->
p();
1071 double eop = Eemc/
ptrk;
1075 depth=mucTrk->
depth();
1078 if( fabs(dedxchi)<5 && fabs(Eemc)>0.15 && fabs(Eemc)<0.3 && (depth>=80*
ptrk-60 || depth>40))
1089 Hep3Vector xorigin(0,0,0);
1091 Gaudi::svcLocator()->service(
"VertexDbSvc", vtxsvc);
1095 xorigin.setX(dbv[0]);
1096 xorigin.setY(dbv[1]);
1097 xorigin.setZ(dbv[2]);
1104 HepSymMatrix Ea = mdcKalTrk->
getZError();
1106 HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]);
1109 HepVector
vec = helixp.
a();
1110 double vrl =
vec[0];
1111 double vzl =
vec[3];
1114 if(fabs(vrl)<1 && fabs(vzl)<10 && fabs(
costheta)<0.93)
1125 Hep3Vector xorigin(0,0,0);
1127 Gaudi::svcLocator()->service(
"VertexDbSvc", vtxsvc);
1131 xorigin.setX(dbv[0]);
1132 xorigin.setY(dbv[1]);
1133 xorigin.setZ(dbv[2]);
1140 HepSymMatrix Ea = mdcKalTrk->
getZError();
1142 HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]);
1145 HepVector
vec = helixp.
a();
1146 double vrl =
vec[0];
1147 double vzl =
vec[3];
1150 if(fabs(vzl)<20 && fabs(
costheta)<0.93){
1160 double eraw = emcTrk->
energy();
1161 double time = emcTrk->
time();
1162 double costh =
cos(emcTrk->
theta());
1164 (fabs(costh)<0.80 && eraw>0.025)
1165 || (fabs(costh)>0.86 &&fabs(costh)<0.92 && eraw>0.05)
1166 ) && time>0 && time<14 )
1175 Hep3Vector xorigin(0,0,0);
1177 Gaudi::svcLocator()->service(
"VertexDbSvc", vtxsvc);
1181 xorigin.setX(dbv[0]);
1182 xorigin.setY(dbv[1]);
1183 xorigin.setZ(dbv[2]);
1189 Hep3Vector emcpos(emcTrk->
x(), emcTrk->
y(), emcTrk->
z());
1190 Hep3Vector Gm_Vec(emcTrk->
x(), emcTrk->
y(), emcTrk->
z());
1191 Hep3Vector Gm_Mom = Gm_Vec - xorigin;
1192 Gm_Mom.setMag(emcTrk->
energy());
1193 HepLorentzVector shP4(Gm_Mom, emcTrk->
energy());
1195 double costh = shP4.vect().cosTheta();
1196 double eraw = emcTrk->
energy();
1197 double time = emcTrk->
time();
1198 if(time<0 || time>14)
return false;
1199 if(!((fabs(costh)<0.80&&eraw>0.025) || (fabs(costh)>0.86&&fabs(costh)<0.92&&eraw>0.05)))
return false;
1201 for(
int j = 0; j < evtRecEvent->totalCharged(); j++){
1203 if(!(*jtTrk)->isExtTrackValid())
continue;
1207 double angd = extpos.angle(emcpos);
1208 if(angd < dang) dang = angd;
1210 if(dang>=200)
return false;
1211 dang = dang * 180 / (CLHEP::pi);
1212 if(fabs(dang) < shwDang)
return false;
1220 vector<EvtRecTrackIterator> iGood;
1224 iGood.push_back(
iter);
1227 if(iGood.size() != 2)
1231 double time1=-99,time2=-99;
1232 for(vector<EvtRecTrackIterator>::size_type i=0;i<2;i++){
1233 if( (*iGood[i])->isTofTrackValid() ){
1234 SmartRefVector<RecTofTrack> tofTrkCol= (*iGood[i])->tofTrack();
1235 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
1237 for(;iter_tof!=tofTrkCol.end();iter_tof++){
1239 status->
setStatus( (*iter_tof)->status() );
1242 time1=(*iter_tof)->tof();
1244 time2=(*iter_tof)->tof();
1250 if( time1!=-99 && time2!=-99 && fabs(time1-time2)>5)
1265 bool gotgoodshower=
false;
1269 if(!(*iter)->isEmcShowerValid())
continue;
1272 double eraw = emcTrk->
energy();
1273 double time = emcTrk->
time();
1274 if( !(eraw>0.05 && time>0 && time<14 ))
1281 if(angle1>20 && angle2>20){
1289 return gotgoodshower;
1301 Hep3Vector emcpos(emcTrk->
x(), emcTrk->
y(), emcTrk->
z());
1308 if(extpos==(-99,-99,-99))
1311 return extpos.angle(emcpos)*180/CLHEP::pi;
1317 SmartRefVector<EvtRecTrack> tracks1=(*iter1)->tracks();
1318 SmartRefVector<EvtRecTrack> showers1=(*iter1)->showers();
1319 SmartRefVector<EvtRecTrack> tracks2=(*iter2)->tracks();
1320 SmartRefVector<EvtRecTrack> showers2=(*iter2)->showers();
1323 for(
int i=0; i<tracks1.size(); i++){
1324 for(
int j=0; j<tracks2.size(); j++){
1325 if(tracks1[i]==tracks2[j])
1331 for(
int i=0; i<showers1.size(); i++){
1332 for(
int j=0; j<showers2.size(); j++){
1333 if(showers1[i]==showers2[j])
1345 StatusCode sc = Gaudi::svcLocator()->service (
"EventDataSvc", m_evtSvc,
true);
1346 if( sc.isFailure() ) {
double sin(const BesAngle a)
double cos(const BesAngle a)
DOUBLE_PRECISION count[3]
EvtRecTrackCol::iterator EvtRecTrackIterator
const Hep3Vector emcPosition() const
const int emcVolumeNumber() const
const double theta() const
static void setPidType(PidType pidType)
bool isMdcKalTrackValid()
RecEmcShower * emcShower()
RecMdcKalTrack * mdcKalTrack()
virtual void preparePID(EvtRecTrack *track)=0
virtual bool iselectron(bool emc=false)=0
virtual bool isVertexValid()=0
virtual double * SigmaPrimaryVertex()=0
virtual double * PrimaryVertex()=0
const HepVector & getZHelix() const
HepVector & getZHelixMu()
const HepSymMatrix & getZError() const
void setStatus(unsigned int status)
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