114 {
115 StatusCode sc = StatusCode::SUCCESS;
116
117 MsgStream log(
msgSvc(), name() );
118 log << MSG::INFO << "in execute()" << endmsg;
119
120
121
122
123 setFilterPassed( false );
124
125 m_pass[0] += 1;
126
127 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
130
131 Vint iks, ipip, ipim, iGood;
132 iGood.clear();
133 iks.clear();
134 ipip.clear();
135 ipim.clear();
136
138 ppip.clear();
139 ppim.clear();
140
141 int TotCharge = 0;
142 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
143 {
145 if ( !( *itTrk )->isMdcTrackValid() ) continue;
146 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
147 if ( fabs( mdcTrk->
z() ) >= m_vz0cut )
continue;
148 if ( mdcTrk->
r() >= m_vr0cut )
continue;
149 iGood.push_back( i );
150 TotCharge += mdcTrk->
charge();
151 }
152
153
154
155 int nGood = iGood.size();
156
157
158
159
160
161 if ( ( nGood < 2 ) || ( TotCharge != 0 ) ) return sc;
162
163 m_pass[1] += 1;
164
165
166
168 for ( int i = 0; i < nGood; i++ )
169 {
171
178
179
182 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
183 RecMdcKalTrack* mdcKalTrk = ( *itTrk )->mdcKalTrack();
184
185
187
189 HepLorentzVector ptrk;
190 ptrk.setPx( mdcKalTrk->
px() );
191 ptrk.setPy( mdcKalTrk->
py() );
192 ptrk.setPz( mdcKalTrk->
pz() );
193 double p3 = ptrk.mag();
194 ptrk.setE( sqrt( p3 * p3 +
xmass[2] *
xmass[2] ) );
195
196 if ( mdcKalTrk->
charge() > 0 )
197 {
198 ipip.push_back( iGood[i] );
199 ppip.push_back( ptrk );
200 }
201 else
202 {
203 ipim.push_back( iGood[i] );
204 ppim.push_back( ptrk );
205 }
206 }
207
208 m_pass[2] += 1;
209 int npip = ipip.size();
210 int npim = ipim.size();
211 m_npip = npip;
212 m_npim = npim;
213
214 if ( npip < 1 || npim < 1 ) return sc;
215
216 m_pass[3] += 1;
217
218
219
220
221 double chi, delm;
222 double chisq = 999.;
223
224
226 HepSymMatrix Evx( 3, 0 );
227 double bx = 1E+6;
228 double by = 1E+6;
229 double bz = 1E+6;
230 Evx[0][0] = bx * bx;
231 Evx[1][1] = by * by;
232 Evx[2][2] = bz * bz;
233 VertexParameter vxpar;
236
237
238
239
240
241
242
243
244
245
246 int ip1 = -1, ip2 = -1;
247
250
251 for ( int i1 = 0; i1 < ipip.size(); i1++ )
252 {
253 RecMdcKalTrack* pipTrk = ( *( evtRecTrkCol->begin() + ipip[i1] ) )->mdcKalTrack();
255 WTrackParameter wpiptrk(
xmass[2], pipTrk->
helix(), pipTrk->
err() );
256
257 for ( int i2 = 0; i2 < ipim.size(); i2++ )
258 {
259 RecMdcKalTrack* pimTrk = ( *( evtRecTrkCol->begin() + ipim[i2] ) )->mdcKalTrack();
261 WTrackParameter wpimtrk(
xmass[2], pimTrk->
helix(), pimTrk->
err() );
262
267 if ( !( vtxfit0->
Fit( 0 ) ) )
continue;
269
273 if ( !vtxfit->
Fit() )
continue;
274 chi = vtxfit->
chisq();
275
276
277 if ( chi > chisq ) continue;
278 delm = fabs( ( vtxfit->
p4par().m() ) -
mk0 );
279 chisq = chi;
280 ip1 = ipip[i1];
281 ip2 = ipim[i2];
282 }
283 }
284
285
286
287 if ( ip1 == -1 || ip2 == -1 ) return sc;
288 m_pass[4] += 1;
289
290 RecMdcKalTrack* pi1Trk = ( *( evtRecTrkCol->begin() + ip1 ) )->mdcKalTrack();
292 WTrackParameter wpi1trk(
xmass[2], pi1Trk->
helix(), pi1Trk->
err() );
293
294 RecMdcKalTrack* pi2Trk = ( *( evtRecTrkCol->begin() + ip2 ) )->mdcKalTrack();
296 WTrackParameter wpi2trk(
xmass[2], pi2Trk->
helix(), pi2Trk->
err() );
301 if ( !( vtxfit0->
Fit( 0 ) ) )
return sc;
303
307 if ( !vtxfit->
Fit() )
return sc;
308
309 m_ksRawMass = vtxfit->
p4par().m();
310 m_ctau = vtxfit->
ctau();
313 m_chis = vtxfit->
chisq();
314 m_pk0 = vtxfit->
p4par().rho();
315
316
317 if ( m_lxyz > 0.4 && m_chis < 20. )
318 {
319
320 TH1* h( 0 );
321 if ( m_thsvc->getHist( "/DQAHist/InclKs/InclKs_mass", h ).isSuccess() )
322 { h->Fill( m_ksRawMass ); }
323 else { log << MSG::ERROR << "Couldn't retrieve inclks_mass" << endmsg; }
324 }
325
326 m_tuple1->write().ignore();
327
328
329
330
331
332
333 ( *( evtRecTrkCol->begin() + ip1 ) )->tagPion();
334 ( *( evtRecTrkCol->begin() + ip2 ) )->tagPion();
335
336
337
338
339
340 ( *( evtRecTrkCol->begin() + ip1 ) )->setQuality( 2 );
341 ( *( evtRecTrkCol->begin() + ip2 ) )->setQuality( 2 );
342
343
344
345 setFilterPassed( true );
346
347 return StatusCode::SUCCESS;
348}
HepGeom::Point3D< double > HepPoint3D
EvtRecTrackCol::iterator EvtRecTrackIterator
std::vector< HepLorentzVector > Vp4
const HepVector & helix() const
static void setPidType(PidType pidType)
const HepSymMatrix & err() const
int methodProbability() const
void setChiMinCut(const double chi=4)
void setRecTrack(EvtRecTrack *trk)
void setMethod(const int method)
void identify(const int pidcase)
void usePidSys(const int pidsys)
static ParticleID * instance()
bool IsPidInfoValid() const
HepLorentzVector p4par() const
double decayLength() const
double decayLengthError() const
static SecondVertexFit * instance()
void setVpar(const VertexParameter vpar)
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
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
_EXTERN_ std::string EvtRecTrackCol