00126 {
00127 StatusCode sc = StatusCode::SUCCESS;
00128
00129 MsgStream log(msgSvc(), name());
00130 log << MSG::INFO << "in execute()" << endreq;
00131
00132
00133
00134
00135 setFilterPassed(false);
00136
00137 m_pass[0] += 1;
00138
00139 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00140 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00141 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00142
00143 Vint ilam, ip, ipi, iGood;
00144 iGood.clear();
00145 ilam.clear();
00146 ip.clear();
00147 ipi.clear();
00148
00149 Vp4 ppi, pp;
00150 ppi.clear();
00151 pp.clear();
00152
00153 int TotCharge = 0;
00154 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00155 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00156 if(!(*itTrk)->isMdcTrackValid()) continue;
00157 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00158 if(fabs(mdcTrk->z()) >= m_vz0cut) continue;
00159 if(mdcTrk->r() >= m_vr0cut) continue;
00160 iGood.push_back(i);
00161 TotCharge += mdcTrk->charge();
00162 }
00163
00164
00165
00166 int nGood = iGood.size();
00167
00168
00169
00170
00171
00172 if((nGood < 2) || (TotCharge!=0)) return sc;
00173
00174 m_pass[1] += 1;
00175
00176
00177
00178 ParticleID *pid = ParticleID::instance();
00179 for(int i = 0; i < nGood; i++) {
00180 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00181
00182 pid->init();
00183 pid->setMethod(pid->methodProbability());
00184 pid->setChiMinCut(4);
00185 pid->setRecTrack(*itTrk);
00186 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());
00187 pid->identify(pid->onlyPion() | pid->onlyProton());
00188 pid->calculate();
00189 if(!(pid->IsPidInfoValid())) continue;
00190
00191 if(!(*itTrk)->isMdcKalTrackValid()) continue ;
00192 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
00193
00194 if(mdcKalTrk->charge() == 0) continue;
00195
00196 if(pid->probPion() > 0.001 && (pid->probPion() > pid->probProton())){
00197 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
00198 ipi.push_back(iGood[i]);
00199 HepLorentzVector ptrk;
00200 ptrk.setPx(mdcKalTrk->px());
00201 ptrk.setPy(mdcKalTrk->py());
00202 ptrk.setPz(mdcKalTrk->pz());
00203 double p3 = ptrk.mag();
00204 ptrk.setE(sqrt(p3*p3+mpi*mpi));
00205 ppi.push_back(ptrk);
00206 }
00207 if(pid->probProton() > 0.001 && (pid->probProton() > pid->probPion())){
00208 mdcKalTrk->setPidType(RecMdcKalTrack::proton);
00209 ip.push_back(iGood[i]);
00210 HepLorentzVector ptrk;
00211 ptrk.setPx(mdcKalTrk->px());
00212 ptrk.setPy(mdcKalTrk->py());
00213 ptrk.setPz(mdcKalTrk->pz());
00214 double p3 = ptrk.mag();
00215 ptrk.setE(sqrt(p3*p3+mp*mp));
00216 pp.push_back(ptrk);
00217 }
00218 }
00219
00220 m_pass[2] += 1;
00221
00222 int npi = ipi.size();
00223 int np = ip.size();
00224 m_npi = npi;
00225 m_np = np;
00226 if(npi < 1 || np <1) return sc;
00227
00228 m_pass[3] += 1;
00229
00230
00231
00232 double chi=999.;
00233 double delm;
00234
00235 HepPoint3D vx(0., 0., 0.);
00236 HepSymMatrix Evx(3, 0);
00237 double bx = 1E+6;
00238 double by = 1E+6;
00239 double bz = 1E+6;
00240 Evx[0][0] = bx*bx;
00241 Evx[1][1] = by*by;
00242 Evx[2][2] = bz*bz;
00243 VertexParameter vxpar;
00244 vxpar.setVx(vx);
00245 vxpar.setEvx(Evx);
00246
00247 int ipion=-1, iproton=-1;
00248
00249 VertexFit *vtxfit0 = VertexFit::instance();
00250 SecondVertexFit *vtxfit = SecondVertexFit::instance();
00251
00252 for(unsigned int i1 = 0; i1 < ipi.size(); i1++) {
00253 RecMdcKalTrack *piTrk = (*(evtRecTrkCol->begin()+ipi[i1]))->mdcKalTrack();
00254 piTrk->setPidType(RecMdcKalTrack::pion);
00255 WTrackParameter wpitrk(mpi, piTrk->helix(), piTrk->err());
00256
00257 for(unsigned int i2 = 0; i2 < ip.size(); i2++) {
00258 RecMdcKalTrack *protonTrk = (*(evtRecTrkCol->begin()+ip[i2]))->mdcKalTrack();
00259 protonTrk->setPidType(RecMdcKalTrack::proton);
00260 WTrackParameter wptrk(mp, protonTrk->helix(), protonTrk->err());
00261
00262 int NCharge = 0;
00263 NCharge =piTrk->charge() + protonTrk->charge();
00264
00265
00266
00267 if(NCharge != 0) continue;
00268
00269 vtxfit0->init();
00270 vtxfit0->AddTrack(0, wpitrk);
00271 vtxfit0->AddTrack(1, wptrk);
00272 vtxfit0->AddVertex(0, vxpar, 0, 1);
00273 if(!(vtxfit0->Fit(0))) continue;
00274 vtxfit0->BuildVirtualParticle(0);
00275
00276 vtxfit->init();
00277
00278 vtxfit->AddTrack(0, vtxfit0->wVirtualTrack(0));
00279 vtxfit->setVpar(vtxfit0->vpar(0));
00280 if(!vtxfit->Fit()) continue;
00281
00282
00283
00284
00285 if(vtxfit->chisq() > chi) continue;
00286
00287 delm = fabs((vtxfit->p4par()).m()-mlam);
00288 chi = vtxfit->chisq();
00289 ipion=ipi[i1];
00290 iproton=ip[i2];
00291 }
00292 }
00293
00294
00295
00296 if(ipion==-1 || iproton==-1) return sc;
00297 m_pass[4] += 1;
00298
00299 RecMdcKalTrack *pionTrk = (*(evtRecTrkCol->begin()+ipion))->mdcKalTrack();
00300 pionTrk->setPidType(RecMdcKalTrack::pion);
00301 WTrackParameter wpiontrk(mpi, pionTrk->helix(), pionTrk->err());
00302
00303 RecMdcKalTrack *protonTrk = (*(evtRecTrkCol->begin()+iproton))->mdcKalTrack();
00304 protonTrk->setPidType(RecMdcKalTrack::proton);
00305 WTrackParameter wprotontrk(mp, protonTrk->helix(), protonTrk->err());
00306 vtxfit0->init();
00307 vtxfit0->AddTrack(0, wpiontrk);
00308 vtxfit0->AddTrack(1, wprotontrk);
00309 vtxfit0->AddVertex(0, vxpar, 0, 1);
00310 if(!(vtxfit0->Fit(0))) return sc;
00311 vtxfit0->BuildVirtualParticle(0);
00312
00313 vtxfit->init();
00314
00315 vtxfit->AddTrack(0, vtxfit0->wVirtualTrack(0));
00316 vtxfit->setVpar(vtxfit0->vpar(0));
00317 if(!vtxfit->Fit()) return sc;
00318
00319 m_lamRawMass = vtxfit->p4par().m();
00320 m_ctau = vtxfit->ctau();
00321 m_lxyz = vtxfit->decayLength();
00322 m_elxyz = vtxfit->decayLengthError();
00323 m_chis = vtxfit->chisq();
00324 m_plam = vtxfit->p4par().rho();
00325
00326 if(m_lxyz>3.5){
00327
00328 TH1 *h(0);
00329 if (m_thsvc->getHist("/DQAHist/InclLam/InclLam_mass", h).isSuccess()) {
00330 h->Fill(m_lamRawMass);
00331 } else {
00332 log << MSG::ERROR << "Couldn't retrieve incllam_mass" << endreq;
00333 }
00334
00335 }
00336 m_tuple1->write();
00338
00339
00340
00341 (*(evtRecTrkCol->begin()+ipion))->setPartId(3);
00342 (*(evtRecTrkCol->begin()+iproton))->setPartId(5);
00343
00344
00345
00346
00347
00348 (*(evtRecTrkCol->begin()+ipion))->setQuality(1);
00349 (*(evtRecTrkCol->begin()+iproton))->setQuality(1);
00350
00351
00352
00353
00354 setFilterPassed(true);
00355
00356 return StatusCode::SUCCESS;
00357 }