128 {
129 StatusCode sc = StatusCode::SUCCESS;
130
131 MsgStream log(
msgSvc(), name() );
132 log << MSG::INFO << "in execute()" << endmsg;
133
134
135
136
137 setFilterPassed( false );
138
139 m_pass[0] += 1;
140
141 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
144
145 Vint ilam, ip, ipi, iGood;
146 iGood.clear();
147 ilam.clear();
148 ip.clear();
149 ipi.clear();
150
152 ppi.clear();
153 pp.clear();
154
155 int TotCharge = 0;
156 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
157 {
159 if ( !( *itTrk )->isMdcTrackValid() ) continue;
160 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
161 if ( fabs( mdcTrk->
z() ) >= m_vz0cut )
continue;
162 if ( mdcTrk->
r() >= m_vr0cut )
continue;
163 iGood.push_back( i );
164 TotCharge += mdcTrk->
charge();
165 }
166
167
168
169 int nGood = iGood.size();
170
171
172
173
174
175 if ( ( nGood < 2 ) || ( TotCharge != 0 ) ) return sc;
176
177 m_pass[1] += 1;
178
179
180
182 for ( int i = 0; i < nGood; i++ )
183 {
185
194
195 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
196 RecMdcKalTrack* mdcKalTrk = ( *itTrk )->mdcKalTrack();
197
198 if ( mdcKalTrk->
charge() == 0 )
continue;
199
201 {
203 ipi.push_back( iGood[i] );
204 HepLorentzVector ptrk;
205 ptrk.setPx( mdcKalTrk->
px() );
206 ptrk.setPy( mdcKalTrk->
py() );
207 ptrk.setPz( mdcKalTrk->
pz() );
208 double p3 = ptrk.mag();
209 ptrk.setE( sqrt( p3 * p3 +
mpi *
mpi ) );
210 ppi.push_back( ptrk );
211 }
213 {
215 ip.push_back( iGood[i] );
216 HepLorentzVector ptrk;
217 ptrk.setPx( mdcKalTrk->
px() );
218 ptrk.setPy( mdcKalTrk->
py() );
219 ptrk.setPz( mdcKalTrk->
pz() );
220 double p3 = ptrk.mag();
221 ptrk.setE( sqrt( p3 * p3 +
mp *
mp ) );
222 pp.push_back( ptrk );
223 }
224 }
225
226 m_pass[2] += 1;
227
228 int npi = ipi.size();
229 int np = ip.size();
230 m_npi = npi;
231 m_np = np;
232 if ( npi < 1 || np < 1 ) return sc;
233
234 m_pass[3] += 1;
235
236
237
238 double chi = 999.;
239 double delm;
240
242 HepSymMatrix Evx( 3, 0 );
243 double bx = 1E+6;
244 double by = 1E+6;
245 double bz = 1E+6;
246 Evx[0][0] = bx * bx;
247 Evx[1][1] = by * by;
248 Evx[2][2] = bz * bz;
249 VertexParameter vxpar;
252
253 int ipion = -1, iproton = -1;
254
257
258 for ( unsigned int i1 = 0; i1 < ipi.size(); i1++ )
259 {
260 RecMdcKalTrack* piTrk = ( *( evtRecTrkCol->begin() + ipi[i1] ) )->mdcKalTrack();
262 WTrackParameter wpitrk(
mpi, piTrk->
helix(), piTrk->
err() );
263
264 for ( unsigned int i2 = 0; i2 < ip.size(); i2++ )
265 {
266 RecMdcKalTrack* protonTrk = ( *( evtRecTrkCol->begin() + ip[i2] ) )->mdcKalTrack();
268 WTrackParameter wptrk(
mp, protonTrk->
helix(), protonTrk->
err() );
269
270 int NCharge = 0;
272
273
274
275 if ( NCharge != 0 ) continue;
276
281 if ( !( vtxfit0->
Fit( 0 ) ) )
continue;
283
285
288 if ( !vtxfit->
Fit() )
continue;
289
290
291
292
293 if ( vtxfit->
chisq() > chi )
continue;
294
295 delm = fabs( ( vtxfit->
p4par() ).m() -
mlam );
296 chi = vtxfit->
chisq();
297 ipion = ipi[i1];
298 iproton = ip[i2];
299 }
300 }
301
302
303
304 if ( ipion == -1 || iproton == -1 ) return sc;
305 m_pass[4] += 1;
306
307 RecMdcKalTrack* pionTrk = ( *( evtRecTrkCol->begin() + ipion ) )->mdcKalTrack();
309 WTrackParameter wpiontrk(
mpi, pionTrk->
helix(), pionTrk->
err() );
310
311 RecMdcKalTrack* protonTrk = ( *( evtRecTrkCol->begin() + iproton ) )->mdcKalTrack();
313 WTrackParameter wprotontrk(
mp, protonTrk->
helix(), protonTrk->
err() );
318 if ( !( vtxfit0->
Fit( 0 ) ) )
return sc;
320
322
325 if ( !vtxfit->
Fit() )
return sc;
326
327 m_lamRawMass = vtxfit->
p4par().m();
328 m_ctau = vtxfit->
ctau();
331 m_chis = vtxfit->
chisq();
332 m_plam = vtxfit->
p4par().rho();
333
334 if ( m_lxyz > 3.5 )
335 {
336
337 TH1* h( 0 );
338 if ( m_thsvc->getHist( "/DQAHist/InclLam/InclLam_mass", h ).isSuccess() )
339 { h->Fill( m_lamRawMass ); }
340 else { log << MSG::ERROR << "Couldn't retrieve incllam_mass" << endmsg; }
341 }
342 m_tuple1->write();
343
344
345
346
347
348
349 ( *( evtRecTrkCol->begin() + ipion ) )->tagPion();
350 ( *( evtRecTrkCol->begin() + iproton ) )->tagProton();
351
352
353
354
355
356 ( *( evtRecTrkCol->begin() + ipion ) )->setQuality( 1 );
357 ( *( evtRecTrkCol->begin() + iproton ) )->setQuality( 1 );
358
359
360
361
362 setFilterPassed( true );
363
364 return StatusCode::SUCCESS;
365}
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
double probProton() 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