00001
00002
00003
00004
00005 #include "GaudiKernel/MsgStream.h"
00006 #include "GaudiKernel/AlgFactory.h"
00007 #include "GaudiKernel/ISvcLocator.h"
00008 #include "GaudiKernel/SmartDataPtr.h"
00009 #include "GaudiKernel/IDataProviderSvc.h"
00010 #include "GaudiKernel/PropertyMgr.h"
00011
00012 #include "EventModel/Event.h"
00013 #include "EvtRecEvent/EvtRecEvent.h"
00014 #include "EvtRecEvent/EvtRecTrack.h"
00015 #include "EvtRecEvent/EvtRecVeeVertex.h"
00016 #include "EventModel/EventHeader.h"
00017 #include "EventModel/Event.h"
00018 #include "EventModel/EventModel.h"
00019 #include "McTruth/McParticle.h"
00020 #include "ParticleID/ParticleID.h"
00021 #include "VertexFit/VertexFit.h"
00022 #include "VertexFit/SecondVertexFit.h"
00023 #include "VeeVertexAlg/KShortReconstruction.h"
00024 #include <vector>
00025
00026 typedef std::vector<int> Vint;
00027 const double mpi = 0.13957;
00028
00029 KShortReconstruction::KShortReconstruction(const std::string& name, ISvcLocator* pSvcLocator):
00030 Algorithm(name, pSvcLocator) {
00031
00032
00033 declareProperty("Vz0Cut", m_vzCut = 20.0);
00034 declareProperty("CosThetaCut", m_cosThetaCut=0.93);
00035 declareProperty("ChisqCut", m_chisqCut = 100.0);
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047 }
00048
00049
00050 StatusCode KShortReconstruction::initialize() {
00051 MsgStream log(msgSvc(), name());
00052 log << MSG::INFO << "in initialize()" <<endreq;
00053
00054 return StatusCode::SUCCESS;
00055 }
00056
00057
00058 StatusCode KShortReconstruction::finalize() {
00059 MsgStream log(msgSvc(), name());
00060 log << MSG::INFO << "in finalize()" << endreq;
00061
00062 return StatusCode::SUCCESS;
00063 }
00064
00065 StatusCode KShortReconstruction::registerEvtRecVeeVertexCol(
00066 EvtRecVeeVertexCol* aNewEvtRecVeeVertexCol, MsgStream& log) {
00067 StatusCode sc = eventSvc()->registerObject("/Event/EvtRec/EvtRecVeeVertexCol",
00068 aNewEvtRecVeeVertexCol);
00069 if (sc != StatusCode::SUCCESS) {
00070 log << MSG::FATAL << "Could not register EvtRecVeeVertexCol in TDS!" << endreq;
00071 }
00072 return sc;
00073 }
00074
00075
00076 StatusCode KShortReconstruction::execute() {
00077 MsgStream log(msgSvc(), name());
00078 log << MSG::INFO << "in execute()" << endreq;
00079
00080 StatusCode sc;
00081
00083
00085 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00086 SmartDataPtr<EvtRecEvent> recEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00087 SmartDataPtr<EvtRecTrackCol> recTrackCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00088 log << MSG::DEBUG << "run and event = " << eventHeader->runNumber()
00089 << " " << eventHeader->eventNumber() << endreq;
00090 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00091 << recEvent->totalCharged() << " , "
00092 << recEvent->totalNeutral() << " , "
00093 << recEvent->totalTracks() <<endreq;
00094 int evtNo = eventHeader->eventNumber();
00095
00096
00097
00098
00099 SmartDataPtr<EvtRecVeeVertexCol> veeVertexCol(eventSvc(),
00100 EventModel::EvtRec::EvtRecVeeVertexCol);
00101 if (!veeVertexCol) {
00102 veeVertexCol = new EvtRecVeeVertexCol;
00103 sc = registerEvtRecVeeVertexCol(veeVertexCol, log);
00104 if (sc != StatusCode::SUCCESS) {
00105 return sc;
00106 }
00107 }
00108
00110
00112 Vint ipip, ipim, iGood;
00113 for (unsigned int i = 0; i < recEvent->totalCharged(); i++) {
00114 EvtRecTrackIterator itTrk = recTrackCol->begin() + i;
00115
00116 if(!(*itTrk)->isMdcTrackValid()) continue;
00117 if(!(*itTrk)->isMdcKalTrackValid()) continue;
00118 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
00119 if (fabs(cos(mdcTrk->theta())) >= m_cosThetaCut) continue;
00120 if (fabs(mdcTrk->z()) >= m_vzCut) continue;
00121 iGood.push_back(i);
00122 if (mdcTrk->charge() > 0) ipip.push_back(i);
00123 if (mdcTrk->charge() < 0) ipim.push_back(i);
00124 }
00125
00126
00127 if (ipip.size() < 1 || ipim.size() < 1) return StatusCode::SUCCESS;
00128
00130
00132 HepPoint3D vx(0., 0., 0.);
00133 HepSymMatrix Evx(3, 0);
00134 double bx = 1E+6;
00135 double by = 1E+6;
00136 double bz = 1E+6;
00137 Evx[0][0] = bx*bx;
00138 Evx[1][1] = by*by;
00139 Evx[2][2] = bz*bz;
00140 VertexParameter vxpar;
00141 vxpar.setVx(vx);
00142 vxpar.setEvx(Evx);
00143
00144 VertexFit *vtxfit0 = VertexFit::instance();
00145
00146 for(unsigned int i1 = 0; i1 < ipip.size(); i1++) {
00147 int ip1 = ipip[i1];
00148 RecMdcKalTrack *pipKalTrk = (*(recTrackCol->begin()+ip1))->mdcKalTrack();
00149 RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
00150 WTrackParameter wpiptrk(mpi, pipKalTrk->helix(), pipKalTrk->err());
00151
00152 for(unsigned int i2 = 0; i2 < ipim.size(); i2++) {
00153 int ip2 = ipim[i2];
00154 RecMdcKalTrack *pimKalTrk = (*(recTrackCol->begin()+ip2))->mdcKalTrack();
00155 RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
00156 WTrackParameter wpimtrk(mpi, pimKalTrk->helix(), pimKalTrk->err());
00157
00158 vtxfit0->init();
00159 vtxfit0->AddTrack(0, wpiptrk);
00160 vtxfit0->AddTrack(1, wpimtrk);
00161 vtxfit0->AddVertex(0, vxpar, 0, 1);
00162 if(!(vtxfit0->Fit(0))) continue;
00163 vtxfit0->BuildVirtualParticle(0);
00164 if(vtxfit0->chisq(0) > m_chisqCut) continue;
00165 WTrackParameter wKshort = vtxfit0->wVirtualTrack(0);
00166 std::pair<int, int> pair;
00167 pair.first = 3;
00168 pair.second = 3;
00169
00170 EvtRecTrack* track0 = *(recTrackCol->begin() + ip1);
00171 EvtRecTrack* track1 = *(recTrackCol->begin() + ip2);
00172
00173 EvtRecVeeVertex* KsVertex = new EvtRecVeeVertex;
00174 KsVertex->setVertexId(310);
00175 KsVertex->setVertexType(0);
00176 KsVertex->setChi2(vtxfit0->chisq(0));
00177 KsVertex->setNdof(1);
00178 KsVertex->setMass(wKshort.p().m());
00179 KsVertex->setW(wKshort.w());
00180 KsVertex->setEw(wKshort.Ew());
00181 KsVertex->setPair(pair);
00182 KsVertex->setNCharge(0);
00183 KsVertex->setNTracks(2);
00184 KsVertex->addDaughter(track0, 0);
00185 KsVertex->addDaughter(track1, 1);
00186 veeVertexCol->push_back(KsVertex);
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196 }
00197 }
00198
00199 return StatusCode::SUCCESS;
00200 }