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, m_chargek1, m_chargek2;
43 int ik1_temp, ik2_temp, ipi1_temp, ipi2_temp;
44 HepLorentzVector pddd;
45 HepLorentzVector pddd_temp;
46
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 k0kkmd = false;
58 double tagmode = 0;
59
60 if ( ( evtRecEvent->totalCharged() < 4 ) ) { 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 HepSymMatrix xoriginEx( 3, 0 );
72 IVertexDbSvc* vtxsvc;
73 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc );
75 {
78 xorigin.setX( dbv[0] );
79 xorigin.setY( dbv[1] );
80 xorigin.setZ( dbv[2] );
81
82 xoriginEx[0][0] = vv[0] * vv[0];
83 xoriginEx[1][1] = vv[1] * vv[1];
84 xoriginEx[2][2] = vv[2] * vv[2];
85 }
86
87 double xv = xorigin.x();
88 double yv = xorigin.y();
89 double zv = xorigin.z();
90
92 HepPoint3D IP( xorigin[0], xorigin[1], xorigin[2] );
93
94
95 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
96 {
98
99 int ipi1 = ( *itTrk )->trackId();
100
101 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
102 RecMdcKalTrack* mdcKalTrk1 = ( *itTrk )->mdcKalTrack();
104
105 m_chargepi1 = mdcKalTrk1->
charge();
106 if (
abs( m_chargepi1 ) != 1 )
continue;
107
108
110 HepSymMatrix Ea1 = mdcKalTrk1->
getZError();
111
112 VFHelix helixip3_1( point0, a1, Ea1 );
113 helixip3_1.pivot( IP );
114 HepVector vecipa1 = helixip3_1.a();
115
116 double dr1 = fabs( vecipa1[0] );
117 double dz1 = fabs( vecipa1[3] );
118 double costheta1 =
cos( mdcKalTrk1->
theta() );
119
120 if ( dr1 >= 15.0 ) continue;
121 if ( dz1 >= 25.0 ) continue;
122 if ( fabs( costheta1 ) >= 0.93 ) continue;
123
125
126
127
128
129 for ( int j = 0; j < evtRecEvent->totalCharged(); j++ )
130 {
132
133 int ipi2 = ( *itTrk )->trackId();
134 if ( ipi1 == ipi2 ) continue;
135
136 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
137 RecMdcKalTrack* mdcKalTrk2 = ( *itTrk )->mdcKalTrack();
139
140 m_chargepi2 = mdcKalTrk2->
charge();
141 if ( ( m_chargepi2 + m_chargepi1 ) != 0 ) continue;
142
143
145 HepSymMatrix Ea2 = mdcKalTrk2->
getZError();
146 VFHelix helixip3_2( point0, a2, Ea2 );
147 helixip3_2.pivot( IP );
148 HepVector vecipa2 = helixip3_2.a();
149
150 double dr2 = fabs( vecipa2[0] );
151 double dz2 = fabs( vecipa2[3] );
152 double costheta2 =
cos( mdcKalTrk2->
theta() );
153 if ( dr2 >= 15.0 ) continue;
154 if ( dz2 >= 25.0 ) continue;
155 if ( fabs( costheta2 ) >= 0.93 ) continue;
156
158
159 HepPoint3D vx( xorigin.x(), xorigin.y(), xorigin.z() );
160 HepSymMatrix Evx( 3, 0 );
161 double bx = 1E+6;
162 Evx[0][0] = bx * bx;
163 double by = 1E+6;
164 Evx[1][1] = by * by;
165 double bz = 1E+6;
166 Evx[2][2] = bz * bz;
167 VertexParameter vxpar;
170
171 HepVector pip_val = HepVector( 7, 0 );
172 HepVector pim_val = HepVector( 7, 0 );
173 pip_val = pip.w();
174 pim_val = pim.w();
175 HepLorentzVector ptrktagk0( pip_val[0] + pim_val[0], pip_val[1] + pim_val[1],
176 pip_val[2] + pim_val[2], pip_val[3] + pim_val[3] );
177 double m_xmtagk0_tem = ptrktagk0.mag();
178 if ( fabs( ptrktagk0.m() - 0.498 ) > 0.1 ) continue;
179
185 if ( !( vtxfit0->
Fit( 0 ) ) )
continue;
188 WTrackParameter wksp = vtxfit0->
wtrk( 0 );
189 WTrackParameter wksm = vtxfit0->
wtrk( 1 );
191 VertexParameter wks_var = vtxfit0->
vpar( 0 );
192
193
194
195
196 for ( int k = 0; k < evtRecEvent->totalCharged(); k++ )
197 {
199
200 int ik1 = ( *itTrk )->trackId();
201 if ( ipi2 == ik1 || ipi1 == ik1 ) continue;
202
203 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
204 RecMdcKalTrack* mdcKalTrk3 = ( *itTrk )->mdcKalTrack();
206
207 m_chargek1 = mdcKalTrk3->
charge();
208 if (
abs( m_chargek1 ) != 1 )
continue;
209
210
213 VFHelix helixip3_3( point0, a3, Ea3 );
214 helixip3_3.pivot( IP );
215 HepVector vecipa3 = helixip3_3.a();
216
217 double dr3 = fabs( vecipa3[0] );
218 double dz3 = fabs( vecipa3[3] );
219 double costheta3 =
cos( mdcKalTrk3->
theta() );
220 if ( dr3 >= 1.0 ) continue;
221 if ( dz3 >= 10.0 ) continue;
222 if ( fabs( costheta3 ) >= 0.93 ) continue;
223
224 if ( PID_flag == 5 )
225 {
227 if ( simple_pid->
probKaon() < 0.0 ||
229 continue;
230 }
231
233
234
235
236
237 for ( int l = 0; l < evtRecEvent->totalCharged(); l++ )
238 {
240
241 int ik2 = ( *itTrk )->trackId();
242 if ( ik2 == ik1 || ik2 == ipi2 || ik2 == ipi1 ) continue;
243
244 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
245 RecMdcKalTrack* mdcKalTrk4 = ( *itTrk )->mdcKalTrack();
247
248 m_chargek2 = mdcKalTrk4->
charge();
249 if ( ( m_chargek1 + m_chargek2 ) != 0 ) continue;
250
251
254 VFHelix helixip3_4( point0, a4, Ea4 );
255 helixip3_4.pivot( IP );
256 HepVector vecipa4 = helixip3_4.a();
257
258 double dr4 = fabs( vecipa4[0] );
259 double dz4 = fabs( vecipa4[3] );
260 double costheta4 =
cos( mdcKalTrk4->
theta() );
261 if ( dr4 >= 1.0 ) continue;
262 if ( dz4 >= 10.0 ) continue;
263 if ( fabs( costheta4 ) >= 0.93 ) continue;
264
265 if ( PID_flag == 5 )
266 {
268 if ( simple_pid->
probKaon() < 0.0 ||
270 continue;
271 }
272
273
275
281 if ( !vtxfit_2->
Fit( 0 ) )
continue;
283
284 WTrackParameter wkap = vtxfit_2->
wtrk( 0 );
285 WTrackParameter wkam = vtxfit_2->
wtrk( 1 );
286
289 vxpar.
setEvx( xoriginEx );
293 if ( !vtxfit->
Fit() )
continue;
294
295 if ( vtxfit->
chisq() > 999. )
continue;
297
298 double m_massks1_tem = vtxfit->
p4par().m();
299 if ( m_massks1_tem < 0.485 || m_massks1_tem > 0.515 ) continue;
300 HepLorentzVector p4kstag = vtxfit->
p4par();
302
303 HepVector kap_val = HepVector( 7, 0 );
304 HepVector kam_val = HepVector( 7, 0 );
305 HepVector ksp_val = HepVector( 7, 0 );
306 HepVector ksm_val = HepVector( 7, 0 );
307
312
313 HepLorentzVector P_KAP( kap_val[0], kap_val[1], kap_val[2], kap_val[3] );
314 HepLorentzVector P_KAM( kam_val[0], kam_val[1], kam_val[2], kam_val[3] );
315 HepLorentzVector P_KSP( ksp_val[0], ksp_val[1], ksp_val[2], ksp_val[3] );
316 HepLorentzVector P_KSM( ksm_val[0], ksm_val[1], ksm_val[2], ksm_val[3] );
317
318 P_KAM.boost( -0.011, 0, 0 );
319 P_KAP.boost( -0.011, 0, 0 );
320 P_KSP.boost( -0.011, 0, 0 );
321 P_KSM.boost( -0.011, 0, 0 );
322 p4kstag.boost( -0.011, 0, 0 );
323
324 pddd = P_KAM + P_KAP + p4kstag;
325
326 double pk0kk = pddd.rho();
327
328 double temp1 = (
ecms / 2 ) * (
ecms / 2 ) - pk0kk * pk0kk;
329 if ( temp1 < 0 ) temp1 = 0;
330 double mass_bc_tem = sqrt( temp1 );
331 if ( mass_bc_tem < 1.82 || mass_bc_tem > 1.89 ) continue;
332
333 double delE_tag_tag =
ecms / 2 - pddd.e();
334
335 if ( fabs( delE_tag_tag ) < deltaE_tem )
336 {
337 deltaE_tem = fabs( delE_tag_tag );
338 delE_tag_temp = delE_tag_tag;
339 mass_bcgg = mass_bc_tem;
340
341 pddd_temp = pddd;
342
343 ipi1_temp = ipi1;
344 ipi2_temp = ipi2;
345 ik1_temp = ik1;
346 ik2_temp = ik2;
347
348 ncount1 = 1;
349 }
350 }
351 }
352 }
353 }
354
355 if ( ncount1 == 1 )
356 {
357 tagmode = 22;
358 if ( m_chargetag < 0 ) tagmode = -22;
359 tagmd = tagmode;
360 mass_bc = mass_bcgg;
361 delE_tag = delE_tag_temp;
362 cqtm = 0.0;
363
364 iGoodtag.push_back( ipi1_temp );
365 iGoodtag.push_back( ipi2_temp );
366 iGoodtag.push_back( ik1_temp );
367 iGoodtag.push_back( ik2_temp );
368
369 iGamtag.push_back( 9999 );
370 iGamtag.push_back( 9999 );
371 iGamtag.push_back( 9999 );
372 iGamtag.push_back( 9999 );
373
374 ptag = pddd_temp;
375
376 k0kkmd = true;
377 }
378}
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
HepSymMatrix & getZErrorK()
HepLorentzVector p4par() const
void setPrimaryVertex(const VertexParameter vpar)
double decayLength() const
static SecondVertexFit * instance()
void setVpar(const VertexParameter vpar)
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
WTrackParameter wtrk(int n) const
WTrackParameter wVirtualTrack(int n) const
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
static VertexFit * instance()
VertexParameter vpar(int n) const
void BuildVirtualParticle(int number)
void setEvx(const HepSymMatrix &eVx)
void setVx(const HepPoint3D &vx)
_EXTERN_ std::string EvtRecEvent