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