33 {
34
35 int nGood = iGood.size();
36 int nGam = iGam.size();
37
38 iGoodtag.clear();
39 iGamtag.clear();
40
41 double mass_bcgg, delE_tag_temp;
42 int m_chargetag, m_chargepi1, m_chargepi2;
43 int ipi1_temp, ipi2_temp, ipi3_temp;
44 HepLorentzVector kmfit1, kmfit2, kmfit3, kmfit4, pddd;
45
46 HepLorentzVector pddd_temp;
47 IDataProviderSvc* eventSvc = NULL;
48 Gaudi::svcLocator()->service( "EventDataSvc", eventSvc );
50 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc, "/Event/EventHeader" );
51
52 int runNo = eventHeader->runNumber();
53 int rec = eventHeader->eventNumber();
54
55 double xecm = 2 * Ebeam;
56
57 pipimd = false;
58 double tagmode = 0;
59
60 if ( ( evtRecEvent->totalCharged() < 2 ) ) { return; }
61
63
64 ISimplePIDSvc* simple_pid;
65 Gaudi::svcLocator()->service( "SimplePIDSvc", simple_pid );
66
67 double deltaE_tem = 0.20;
68 int ncount1 = 0;
69
70 Hep3Vector xorigin( 0, 0, 0 );
71 IVertexDbSvc* vtxsvc;
72 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc );
74 {
77 xorigin.setX( dbv[0] );
78 xorigin.setY( dbv[1] );
79 xorigin.setZ( dbv[2] );
80 }
81
82 double xv = xorigin.x();
83 double yv = xorigin.y();
84 double zv = xorigin.z();
85
87 HepPoint3D IP( xorigin[0], xorigin[1], xorigin[2] );
88
89 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
90 {
92
93 int ipi1 = ( *itTrk )->trackId();
94
95 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
96 RecMdcKalTrack* mdcKalTrk1 = ( *itTrk )->mdcKalTrack();
98
99 m_chargepi1 = mdcKalTrk1->
charge();
100 if (
abs( m_chargepi1 ) != 1 )
continue;
101
102
104 HepSymMatrix Ea1 = mdcKalTrk1->
getZError();
105
106 VFHelix helixip3_1( point0, a1, Ea1 );
107 helixip3_1.pivot( IP );
108 HepVector vecipa1 = helixip3_1.a();
109
110 double dr1 = fabs( vecipa1[0] );
111 double dz1 = fabs( vecipa1[3] );
112 double costheta1 =
cos( mdcKalTrk1->
theta() );
113
114 if ( dr1 >= 1.0 ) continue;
115 if ( dz1 >= 10.0 ) continue;
116 if ( fabs( costheta1 ) >= 0.93 ) continue;
117
118 if ( PID_flag == 5 )
119 {
122 continue;
123 }
124
126
127
128
129
130 for ( int j = 0; j < evtRecEvent->totalCharged(); j++ )
131 {
133
134 int ipi2 = ( *itTrk )->trackId();
135 if ( ipi1 == ipi2 ) continue;
136
137 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
138 RecMdcKalTrack* mdcKalTrk2 = ( *itTrk )->mdcKalTrack();
140
141 m_chargepi2 = mdcKalTrk2->
charge();
142 if ( ( m_chargepi1 + m_chargepi2 ) != 0 ) continue;
143
144
146 HepSymMatrix Ea2 = mdcKalTrk2->
getZError();
147 VFHelix helixip3_2( point0, a2, Ea2 );
148 helixip3_2.pivot( IP );
149 HepVector vecipa2 = helixip3_2.a();
150
151 double dr2 = fabs( vecipa2[0] );
152 double dz2 = fabs( vecipa2[3] );
153 double costheta2 =
cos( mdcKalTrk2->
theta() );
154 if ( dr2 >= 1.0 ) continue;
155 if ( dz2 >= 10.0 ) continue;
156 if ( fabs( costheta2 ) >= 0.93 ) continue;
157
158 if ( PID_flag == 5 )
159 {
162 continue;
163 }
164
166
167
168 HepPoint3D vx( xorigin.x(), xorigin.y(), xorigin.z() );
169 HepSymMatrix Evx( 3, 0 );
170 double bx = 1E+6;
171 Evx[0][0] = bx * bx;
172 double by = 1E+6;
173 Evx[1][1] = by * by;
174 double bz = 1E+6;
175 Evx[2][2] = bz * bz;
176 VertexParameter vxpar;
179
180
186 if ( !vtxfit->
Fit( 0 ) )
continue;
188
189 WTrackParameter wpip = vtxfit->
wtrk( 0 );
190 WTrackParameter wpim = vtxfit->
wtrk( 1 );
191
192 HepVector pip_val = HepVector( 7, 0 );
194 HepVector pim_val = HepVector( 7, 0 );
196
197 HepLorentzVector P_PIP( pip_val[0], pip_val[1], pip_val[2], pip_val[3] );
198 HepLorentzVector P_PIM( pim_val[0], pim_val[1], pim_val[2], pim_val[3] );
199
200 P_PIP.boost( -0.011, 0, 0 );
201 P_PIM.boost( -0.011, 0, 0 );
202 pddd = P_PIP + P_PIM;
203
204 double ppipi = pddd.rho();
205
206 double temp1 = (
ecms / 2 ) * (
ecms / 2 ) - ppipi * ppipi;
207 if ( temp1 < 0 ) temp1 = 0;
208 double mass_bc_tem = sqrt( temp1 );
209 if ( mass_bc_tem < 1.82 || mass_bc_tem > 1.89 ) continue;
210
211 double delE_tag_tag =
ecms / 2 - pddd.e();
212
213 if ( fabs( delE_tag_tag ) < deltaE_tem )
214 {
215 deltaE_tem = fabs( delE_tag_tag );
216 delE_tag_temp = delE_tag_tag;
217 mass_bcgg = mass_bc_tem;
218
219 pddd_temp = pddd;
220
221 ipi1_temp = ipi1;
222 ipi2_temp = ipi2;
223
224 ncount1 = 1;
225 }
226 }
227 }
228
229 if ( ncount1 == 1 )
230 {
231 tagmode = 19;
232 if ( m_chargetag < 0 ) tagmode = -19;
233 tagmd = tagmode;
234 mass_bc = mass_bcgg;
235 delE_tag = delE_tag_temp;
236 cqtm = -0.0;
237
238 iGoodtag.push_back( ipi1_temp );
239 iGoodtag.push_back( ipi2_temp );
240
241 iGamtag.push_back( 9999 );
242 iGamtag.push_back( 9999 );
243 iGamtag.push_back( 9999 );
244 iGamtag.push_back( 9999 );
245
246 ptag = pddd_temp;
247
248 pipimd = true;
249 }
250}
HepGeom::Point3D< double > HepPoint3D
EvtRecTrackCol::iterator EvtRecTrackIterator
double cos(const BesAngle a)
const double theta() const
static void setPidType(PidType pidType)
virtual double probKaon()=0
virtual void preparePID(EvtRecTrack *track)=0
virtual double probPion()=0
virtual bool isVertexValid()=0
virtual double * SigmaPrimaryVertex()=0
virtual double * PrimaryVertex()=0
const HepVector & getZHelix() const
const HepSymMatrix & getZError() const
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
WTrackParameter wtrk(int n) const
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
static VertexFit * instance()
void setEvx(const HepSymMatrix &eVx)
void setVx(const HepPoint3D &vx)
_EXTERN_ std::string EvtRecEvent
1.16.1