125 {
126
127 MsgStream log(
msgSvc(), name());
128 log << MSG::INFO << "in execute()" << endreq;
129
130
131
132
133 setFilterPassed(false);
134
135 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
136 m_runNo=eventHeader->runNumber();
137 m_event=eventHeader->eventNumber();
138 log << MSG::DEBUG <<"run, evtnum = "
139 << m_runNo << " , "
140 << m_event <<endreq;
141
142
144 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
145 << evtRecEvent->totalCharged() << " , "
146 << evtRecEvent->totalNeutral() << " , "
147 << evtRecEvent->totalTracks() <<endreq;
148
150
151 if(evtRecEvent->totalNeutral()>100) {
152 return StatusCode::SUCCESS;
153 }
154
155 Vint iGood, iplus, iminus;
156 iGood.clear();
157 iplus.clear();
158 iminus.clear();
160 ppip.clear();
161 ppim.clear();
162
163 Hep3Vector xorigin(0,0,0);
164
166 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
170 xorigin.setX(dbv[0]);
171 xorigin.setY(dbv[1]);
172 xorigin.setZ(dbv[2]);
173 }
174
175 int nCharge = 0;
176 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
178 if(!(*itTrk)->isMdcTrackValid()) continue;
179 if (!(*itTrk)->isMdcKalTrackValid()) continue;
181
182 double pch =mdcTrk->
p();
183 double x0 =mdcTrk->
x();
184 double y0 =mdcTrk->
y();
185 double z0 =mdcTrk->
z();
186 double phi0=mdcTrk->
helix(1);
187 double xv=xorigin.x();
188 double yv=xorigin.y();
189 double Rxy=fabs((x0-xv)*
cos(phi0)+(y0-yv)*
sin(phi0));
190
191 if(fabs(z0) >= m_vz0cut) continue;
192 if(Rxy >= m_vr0cut) continue;
193
194
195
196 TH1 *h(0);
197 if (m_thsvc->getHist("/DQAHist/JsiLL/hrxy", h).isSuccess()) {
198 h->Fill(Rxy);
199 } else {
200 log << MSG::ERROR << "Couldn't retrieve hrxy" << endreq;
201 }
202 if (m_thsvc->getHist("/DQAHist/JsiLL/hz", h).isSuccess()) {
203 h->Fill(z0);
204 } else {
205 log << MSG::ERROR << "Couldn't retrieve hz" << endreq;
206 }
207
208
209 iGood.push_back(i);
210 nCharge += mdcTrk->
charge();
211 if (mdcTrk->
charge() > 0) {
212 iplus.push_back(i);
213 } else {
214 iminus.push_back(i);
215 }
216 }
217
218
219
220
221 int nGood = iGood.size();
222 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
223 if((nGood != 4)||(nCharge!=0)){
224 return StatusCode::SUCCESS;
225 }
226
227
228
229
230
231
232
233
234 int pidp0 = 5, pidp1 = 3, pidm0 = 5, pidm1 = 3;
235
236 RecMdcKalTrack *itTrkp = (*(evtRecTrkCol->begin() + iplus[0]))->mdcKalTrack();
237 RecMdcKalTrack *itTrkpip = (*(evtRecTrkCol->begin() + iplus[1]))->mdcKalTrack();
238 RecMdcKalTrack *itTrkpb = (*(evtRecTrkCol->begin() + iminus[0]))->mdcKalTrack();
239 RecMdcKalTrack *itTrkpim = (*(evtRecTrkCol->begin() + iminus[1]))->mdcKalTrack();
240
241
242
243 if (itTrkp->
p() < itTrkpip->
p()){
244 itTrkp = (*(evtRecTrkCol->begin() + iplus[1]))->mdcKalTrack();
245 itTrkpip = (*(evtRecTrkCol->begin() + iplus[0]))->mdcKalTrack();
246 pidp0 = 3;
247 pidp1 = 5;
248 }
249 if (itTrkpb->
p() < itTrkpim->
p()){
250 itTrkpb = (*(evtRecTrkCol->begin() + iminus[1]))->mdcKalTrack();
251 itTrkpim = (*(evtRecTrkCol->begin() + iminus[0]))->mdcKalTrack();
252 pidm0 = 3;
253 pidm1 = 5;
254 }
255 if (itTrkp->
p() < 0.7 || itTrkp->
p() >1.1)
return StatusCode::SUCCESS;
256 if (itTrkpb->
p() < 0.7 || itTrkpb->
p() >1.1)
return StatusCode::SUCCESS;
257 if (itTrkpip->
p() > 0.35)
return StatusCode::SUCCESS;
258 if (itTrkpim->
p() > 0.35)
return StatusCode::SUCCESS;
259
260
261
266
267
268 m_chisq = 9999.0;
269
271 HepSymMatrix Evx(3, 0);
272 double bx = 1E+6;
273 double by = 1E+6;
274 double bz = 1E+6;
275 Evx[0][0] = bx*bx;
276 Evx[1][1] = by*by;
277 Evx[2][2] = bz*bz;
278
282
288
293
298 if(!vtxfita0->
Fit(0))
return StatusCode::SUCCESS;
304 if(!vtxfita->
Fit())
return StatusCode::SUCCESS;
305
307
312 if(!vtxfitb0->
Fit(0))
return StatusCode::SUCCESS;
315
319 if(!vtxfitb->
Fit())
return StatusCode::SUCCESS;
320
322
327 if(!vtxfit->
Fit(0))
return StatusCode::SUCCESS;
331
332
333
334
336
337
338 HepLorentzVector
ecms(0.034065,0.0,0.0,3.0969);
339 const Hep3Vector
u_cms = -
ecms.boostVector();
340
345
346 if(!kmfit->
Fit())
return StatusCode::SUCCESS;
347 m_chisq = kmfit->
chisq();
348 if(m_chisq > 50) return StatusCode::SUCCESS;
349 HepLorentzVector kmf_pLambda = kmfit->
pfit(0);
350 HepLorentzVector kmf_pLambdabar = kmfit->
pfit(1);
351
352 kmf_pLambda.boost(
u_cms);
353 kmf_pLambdabar.boost(
u_cms);
354 m_mLambda = kmf_pLambda.m();
355 m_mLambdabar = kmf_pLambdabar.m();
356 m_pLambda = kmf_pLambda.rho();
357 m_pLambdabar = kmf_pLambdabar.rho();
358
359 if(fabs(m_mLambda-1.1157)>0.003||fabs(m_mLambdabar-1.1157)>0.003) return StatusCode::SUCCESS;
360
361
362 m_tuple->write();
363
364
365
366
367
368
369
370 (*(evtRecTrkCol->begin()+iplus[0]))->setPartId(pidp0);
371 (*(evtRecTrkCol->begin()+iplus[1]))->setPartId(pidp1);
372 (*(evtRecTrkCol->begin()+iminus[0]))->setPartId(pidm0);
373 (*(evtRecTrkCol->begin()+iminus[1]))->setPartId(pidm1);
374
375
376
377
378
379 (*(evtRecTrkCol->begin()+iplus[0]))->setQuality(0);
380 (*(evtRecTrkCol->begin()+iplus[1]))->setQuality(0);
381 (*(evtRecTrkCol->begin()+iminus[0]))->setQuality(0);
382 (*(evtRecTrkCol->begin()+iminus[1]))->setQuality(0);
383
384
385
386
387 setFilterPassed(true);
388
389
390 return StatusCode::SUCCESS;
391
392}
EvtRecTrackCol::iterator EvtRecTrackIterator
std::vector< HepLorentzVector > Vp4
double sin(const BesAngle a)
double cos(const BesAngle a)
static void setPidType(PidType pidType)
const HepVector helix() const
......
virtual bool isVertexValid()=0
virtual double * SigmaPrimaryVertex()=0
virtual double * PrimaryVertex()=0
static KinematicFit * instance()
void AddFourMomentum(int number, HepLorentzVector p4)
HepLorentzVector pfit(int n) const
const HepVector & getZHelix() const
const HepSymMatrix & getZError() const
HepSymMatrix & getZErrorP()
static SecondVertexFit * instance()
void setVpar(const VertexParameter vpar)
WTrackParameter wpar() const
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
_EXTERN_ std::string EvtRecTrackCol