00001 #include "CLHEP/Vector/LorentzVector.h"
00002 #include "CLHEP/Geometry/Point3D.h"
00003 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00004 typedef HepGeom::Point3D<double> HepPoint3D;
00005 #endif
00006 #include "GaudiKernel/MsgStream.h"
00007 #include "GaudiKernel/SmartDataPtr.h"
00008 #include "GaudiKernel/IDataManagerSvc.h"
00009 #include "EventModel/EventModel.h"
00010 #include "EventModel/Event.h"
00011 #include "EventModel/EventHeader.h"
00012 #include "EvtRecEvent/EvtRecEvent.h"
00013 #include "EvtRecEvent/EvtRecTrack.h"
00014 #include "EvtRecEvent/EvtRecPrimaryVertex.h"
00015 #include "VertexFit/KinematicFit.h"
00016 #include "VertexFit/VertexFit.h"
00017 #include "VertexFit/HTrackParameter.h"
00018 #include "VertexFit/KalmanVertexFit.h"
00019 #include "VertexFit/BField.h"
00020 #include "VertexFit/FastVertexFit.h"
00021 #include "ParticleID/ParticleID.h"
00022 #include "PrimaryVertexAlg/PrimaryVertex.h"
00023
00024 #include <assert.h>
00025 #include "TMath.h"
00026 #include "TH2D.h"
00027 #include "TH1D.h"
00028 #include "TF1.h"
00029 #include <map>
00030 #include <iostream>
00031 #include <fstream>
00032
00033 using CLHEP::HepLorentzVector;
00034 using namespace std;
00035
00036 const double xmass[5] = {0.000511, 0.105658, 0.139570, 0.493677, 0.938272};
00037 const double ecms = 3.097;
00038 const double mpi0 = 0.134976;
00039 const double momega = 0.78265;
00040
00041 typedef std::vector<int> Vint;
00042 typedef std::vector<HepLorentzVector> Vp4;
00043
00044
00045 PrimaryVertex::PrimaryVertex(const std::string& name, ISvcLocator* pSvcLocator) :
00046 Algorithm (name, pSvcLocator)
00047 {
00048
00049 declareProperty("Output", m_output = 0);
00050 declareProperty("TrackNumberCut", m_trackNumberCut = 1);
00051 declareProperty("CosThetaCut", m_cosThetaCut=0.93);
00052 declareProperty("Vz0Cut", m_vz0Cut = 40.0);
00053
00054 declareProperty("FitMethod", m_fitMethod = 0);
00055
00056 declareProperty("Chi2CutOfGlobal", m_globalChisqCut = 200.0);
00057
00058
00059 declareProperty("ChisqCut", m_chisqCut = 200.0);
00060 declareProperty("TrackIteration", m_trackIteration = 5);
00061 declareProperty("VertexIteration", m_vertexIteration = 5);
00062 declareProperty("Chi2CutforTrkIter", m_chi2CutforTrkIter = 0.1);
00063 declareProperty("FreedomCut", m_freedomCut = 1);
00064 declareProperty("Chi2CutforSmooth", m_chi2CutforSmooth = 200.0);
00065
00066
00067
00068 }
00069
00070
00071 StatusCode PrimaryVertex::initialize()
00072 {
00073 MsgStream log(msgSvc(), name());
00074 log << MSG::INFO << "in initialize()" << endmsg;
00075
00076 for(int i = 0; i < 15; i++){
00077 m_sel_number[i] = 0;
00078 }
00079
00080 StatusCode status;
00081
00082 if (m_output == 1) {
00083 NTuplePtr nt1(ntupleSvc(), "FILE999/glbvtx");
00084 if(nt1) m_tuple1 = nt1;
00085 else {
00086 m_tuple1 = ntupleSvc()->book ("FILE999/glbvtx", CLID_ColumnWiseTuple, "global vertex");
00087 if(m_tuple1) {
00088 status = m_tuple1->addItem("gvx", m_gvx);
00089 status = m_tuple1->addItem("gvy", m_gvy);
00090 status = m_tuple1->addItem("gvz", m_gvz);
00091 status = m_tuple1->addItem("chig", m_chig);
00092 status = m_tuple1->addItem("ndofg", m_ndofg);
00093 status = m_tuple1->addItem("probg", m_probg);
00094 } else {
00095 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_tuple1) << endmsg;
00096 return StatusCode::FAILURE;
00097 }
00098 }
00099
00100 NTuplePtr nt2(ntupleSvc(), "FILE999/chisq");
00101 if(nt2) m_tuple2 = nt2;
00102 else {
00103 m_tuple2 = ntupleSvc()->book ("FILE999/chisq", CLID_ColumnWiseTuple, "chi-square of ");
00104 if(m_tuple2) {
00105 status = m_tuple2->addItem("chis", m_chis);
00106 status = m_tuple2->addItem("probs", m_probs);
00107 status = m_tuple2->addItem("chif", m_chif);
00108 status = m_tuple2->addItem("probf", m_probf);
00109 } else {
00110 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_tuple2) << endmsg;
00111 return StatusCode::FAILURE;
00112 }
00113 }
00114
00115 NTuplePtr nt3(ntupleSvc(), "FILE999/Pull");
00116 if(nt3) m_tuple3 = nt3;
00117 else {
00118 m_tuple3 = ntupleSvc()->book ("FILE999/Pull", CLID_ColumnWiseTuple, "Pull");
00119 if(m_tuple3) {
00120 status = m_tuple3->addItem("pull_drho", m_pull_drho);
00121 status = m_tuple3->addItem("pull_phi", m_pull_phi);
00122 status = m_tuple3->addItem("pull_kapha", m_pull_kapha);
00123 status = m_tuple3->addItem("pull_dz", m_pull_dz);
00124 status = m_tuple3->addItem("pull_lamb", m_pull_lamb);
00125 status = m_tuple3->addItem("pull_momentum", m_pull_momentum);
00126 } else {
00127 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_tuple3) << endmsg;
00128 return StatusCode::FAILURE;
00129 }
00130 }
00131
00132 NTuplePtr nt4(ntupleSvc(), "FILE999/kalvtx");
00133 if(nt4) m_tuple4 = nt4;
00134 else {
00135 m_tuple4 = ntupleSvc()->book ("FILE999/kalvtx", CLID_ColumnWiseTuple, "kalman vertex");
00136 if(m_tuple4) {
00137 status = m_tuple4->addItem("kvx", m_kvx);
00138 status = m_tuple4->addItem("kvy", m_kvy);
00139 status = m_tuple4->addItem("kvz", m_kvz);
00140 status = m_tuple4->addItem("chik", m_chik);
00141 status = m_tuple4->addItem("ndofk", m_ndofk);
00142 status = m_tuple4->addItem("probk", m_probk);
00143 } else {
00144 log << MSG::ERROR << "Cannot book N-tuple:" << long(m_tuple4) << endmsg;
00145 return StatusCode::FAILURE;
00146 }
00147 }
00148 }
00149
00150 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00151 return StatusCode::SUCCESS;
00152 }
00153
00154
00155 StatusCode PrimaryVertex::RegisterEvtRecPrimaryVertex(
00156 EvtRecPrimaryVertex* aNewEvtRecPrimaryVertex, MsgStream& log) {
00157 DataObject *aEvtRecEvent;
00158 eventSvc()->findObject("/Event/EvtRec",aEvtRecEvent);
00159 if (aEvtRecEvent==NULL) {
00160 aEvtRecEvent = new EvtRecEvent();
00161 StatusCode sc = eventSvc()->registerObject("/Event/EvtRec",aEvtRecEvent);
00162 if(sc!=StatusCode::SUCCESS) {
00163 log << MSG::FATAL << "Could not register EvtRecEvent" <<endreq;
00164 return StatusCode::FAILURE;
00165 }
00166 }
00167 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
00168 DataObject *aEvtRecPrimaryVertex;
00169 eventSvc()->findObject("/Event/EvtRec/EvtRecPrimaryVertex",aEvtRecPrimaryVertex);
00170 if(aEvtRecPrimaryVertex != NULL) {
00171
00172 dataManSvc->clearSubTree("/Event/EvtRec/EvtRecPrimaryVertex");
00173 eventSvc()->unregisterObject("/Event/EvtRec/EvtRecPrimaryVertex");
00174 }
00175 StatusCode sc = eventSvc()->registerObject("/Event/EvtRec/EvtRecPrimaryVertex",
00176 aNewEvtRecPrimaryVertex);
00177 if (sc != StatusCode::SUCCESS) {
00178 log << MSG::FATAL << "Could not register EvtRecPrimaryVertex in TDS!" << endreq;
00179 return StatusCode::FAILURE;
00180 }
00181
00182 return StatusCode::SUCCESS;
00183 }
00184
00186
00188 void PrimaryVertex::SelectGoodChargedTracks(
00189 SmartDataPtr<EvtRecEvent>& recEvent,
00190 SmartDataPtr<EvtRecTrackCol>& recTrackCol,
00191 Vint& icp, Vint& icm, Vint& iGood) {
00192 assert(recEvent->totalCharged() <= recTrackCol->size());
00193 for (unsigned int i = 0; i < recEvent->totalCharged(); i++) {
00194 EvtRecTrackIterator itTrk = recTrackCol->begin() + i;
00195 if(!(*itTrk)->isMdcTrackValid()) continue;
00196 if(!(*itTrk)->isMdcKalTrackValid()) continue;
00197 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
00198 if (fabs(cos(mdcTrk->theta())) >= m_cosThetaCut) continue;
00199 if (fabs(mdcTrk->z()) >= m_vz0Cut) continue;
00200 iGood.push_back((*itTrk)->trackId());
00201 if (mdcTrk->charge() > 0) {
00202 icp.push_back((*itTrk)->trackId());
00203 } else if (mdcTrk->charge() < 0) {
00204 icm.push_back((*itTrk)->trackId());
00205 }
00206 }
00207 }
00208
00209 void InitVertexParameter(VertexParameter& vxpar) {
00210 HepPoint3D vx(0., 0., 0.);
00211 HepSymMatrix Evx(3, 0);
00212 double bx = 1E+6;
00213 double by = 1E+6;
00214 double bz = 1E+6;
00215 Evx[0][0] = bx*bx;
00216 Evx[1][1] = by*by;
00217 Evx[2][2] = bz*bz;
00218 vxpar.setVx(vx);
00219 vxpar.setEvx(Evx);
00220 }
00221
00222
00223 StatusCode PrimaryVertex::execute()
00224 {
00225 MsgStream log(msgSvc(), name());
00226 log << MSG::INFO << "in execute()" << endmsg;
00227
00228 cout.precision(20);
00230
00232 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00233 SmartDataPtr<EvtRecEvent> recEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00234 SmartDataPtr<EvtRecTrackCol> recTrackCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00235 log << MSG::DEBUG << "run and event = " << eventHeader->runNumber()
00236 << " " << eventHeader->eventNumber() << endreq;
00237 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00238 << recEvent->totalCharged() << " , " << recEvent->totalNeutral() << " , "
00239 << recEvent->totalTracks() <<endreq;
00240 m_sel_number[0]++;
00241
00242
00243
00244
00245
00246 EvtRecPrimaryVertex* aNewEvtRecPrimaryVertex = new EvtRecPrimaryVertex;
00247
00248 StatusCode sc = RegisterEvtRecPrimaryVertex(aNewEvtRecPrimaryVertex, log);
00249 if (sc != StatusCode::SUCCESS) {
00250 return sc;
00251 }
00252
00253
00254 if (recEvent->totalCharged() < m_trackNumberCut) {
00255 return StatusCode::SUCCESS;
00256 }
00257 m_sel_number[1]++;
00258
00259 Vint icp, icm, iGood;
00260 SelectGoodChargedTracks(recEvent, recTrackCol, icp, icm, iGood);
00261
00262
00263 if (iGood.size() < m_trackNumberCut) {
00264 return StatusCode::SUCCESS;
00265 }
00266 m_sel_number[2]++;
00267
00268
00269 VertexParameter vxpar;
00270 InitVertexParameter(vxpar);
00271
00272
00273 if (m_fitMethod == 0) {
00274 VertexFit* vtxfit = VertexFit::instance();
00275 vtxfit->init();
00276 Vint tlis, trkidlis;
00277 for (int i = 0; i < iGood.size(); i++) {
00278 int trk_id = iGood[i];
00279 EvtRecTrackIterator itTrk = recTrackCol->begin() + trk_id;
00280 RecMdcKalTrack *kalTrk = (*itTrk)->mdcKalTrack();
00281 RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
00282 WTrackParameter wtrk(xmass[2], kalTrk->helix(), kalTrk->err());
00283 vtxfit->AddTrack(i ,wtrk);
00284 tlis.push_back(i);
00285 trkidlis.push_back(trk_id);
00286 }
00287 vtxfit->AddVertex(0, vxpar, tlis);
00288 if(!vtxfit->Fit(0)) return StatusCode::SUCCESS;
00289 if(vtxfit->chisq(0) > m_globalChisqCut) return StatusCode::SUCCESS;
00290 if (m_output == 1) {
00291 m_chig = vtxfit->chisq(0);
00292 m_ndofg = 2 * iGood.size() - 3;
00293 m_probg = TMath::Prob(m_chig, m_ndofg);
00294 HepVector glvtx = vtxfit->Vx(0);
00295 m_gvx = glvtx[0];
00296 m_gvy = glvtx[1];
00297 m_gvz = glvtx[2];
00298 m_tuple1->write();
00299 }
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315 aNewEvtRecPrimaryVertex->setNTracks(iGood.size());
00316 aNewEvtRecPrimaryVertex->setTrackIdList(trkidlis);
00317 aNewEvtRecPrimaryVertex->setChi2(vtxfit->chisq(0));
00318 aNewEvtRecPrimaryVertex->setNdof(2 * iGood.size() - 3);
00319 aNewEvtRecPrimaryVertex->setFitMethod(0);
00320 aNewEvtRecPrimaryVertex->setVertex(vtxfit->Vx(0));
00321 aNewEvtRecPrimaryVertex->setErrorVertex(vtxfit->Evx(0));
00322 aNewEvtRecPrimaryVertex->setIsValid(true);
00323
00324 } else if (m_fitMethod == 1) {
00325
00326 KalmanVertexFit *kalvtx = KalmanVertexFit::instance();
00327 kalvtx->init();
00328 kalvtx->initVertex(vxpar);
00329 kalvtx->setChisqCut(m_chisqCut);
00330 kalvtx->setTrackIteration(m_trackIteration);
00331 kalvtx->setVertexIteration(m_vertexIteration);
00332 kalvtx->setTrackIterationCut(m_chi2CutforTrkIter);
00333 for(int i = 0; i < iGood.size(); i++) {
00334 int trk_id = iGood[i];
00335 EvtRecTrackIterator itTrk = recTrackCol->begin()+trk_id;
00336 RecMdcKalTrack *kalTrk = (*itTrk)->mdcKalTrack();
00337 RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
00338 HTrackParameter htrk(kalTrk->helix(), kalTrk->err(), trk_id, 2);
00339 kalvtx->addTrack(htrk);
00340 }
00341 kalvtx->filter();
00342
00343
00344 if(kalvtx->ndof() >= m_freedomCut) {
00345
00346 for(int i = 0; i < iGood.size(); i++) {
00347 if (m_output == 1) {
00348 m_chif = kalvtx->chiF(i);
00349 m_chis = kalvtx->chiS(i);
00350 m_probs = TMath::Prob(m_chis, 2);
00351 m_probf = TMath::Prob(m_chif, 2);
00352 m_tuple2->write();
00353 }
00354 if(kalvtx->chiS(i) > m_chi2CutforSmooth) kalvtx->remove(i);
00355 }
00356 }
00357 if(kalvtx->ndof() >= m_freedomCut) {
00358 for(int i = 0; i < iGood.size(); i++) {
00359 kalvtx->smooth(i);
00360
00361 HepVector Pull(5, 0);
00362 Pull = kalvtx->pull(i);
00363 if (m_output == 1) {
00364 m_pull_drho = Pull[0];
00365 m_pull_phi = Pull[1];
00366 m_pull_kapha = Pull[2];
00367 m_pull_dz = Pull[3];
00368 m_pull_lamb = Pull[4];
00369 m_pull_momentum = kalvtx->pullmomentum(i);
00370 m_tuple3->write();
00371 }
00372 }
00373 if (m_output == 1) {
00374 m_chik = kalvtx->chisq();
00375 m_ndofk = kalvtx->ndof();
00376 m_probk = TMath::Prob(m_chik, m_ndofk);
00377 HepVector xp(3, 0);
00378 xp = kalvtx->x();
00379 m_kvx = xp[0];
00380 m_kvy = xp[1];
00381 m_kvz = xp[2];
00382 m_tuple4->write();
00383 }
00384
00385 m_sel_number[3]++;
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399 aNewEvtRecPrimaryVertex->setNTracks(kalvtx->numTrack());
00400 aNewEvtRecPrimaryVertex->setTrackIdList(kalvtx->trackIDList());
00401 aNewEvtRecPrimaryVertex->setChi2(kalvtx->chisq());
00402 aNewEvtRecPrimaryVertex->setNdof(kalvtx->ndof());
00403 aNewEvtRecPrimaryVertex->setFitMethod(1);
00404 aNewEvtRecPrimaryVertex->setVertex(kalvtx->x());
00405 aNewEvtRecPrimaryVertex->setErrorVertex(kalvtx->Ex());
00406 aNewEvtRecPrimaryVertex->setIsValid(true);
00407
00408 }
00409 }
00410 return StatusCode::SUCCESS;
00411 }
00412
00413
00414
00415 StatusCode PrimaryVertex::finalize()
00416 {
00417 MsgStream log(msgSvc(), name());
00418 log << MSG::INFO << "in finalize()" << endmsg;
00419
00420 log << MSG::ALWAYS << "==================================" << endreq;
00421 log << MSG::ALWAYS << "survived event :";
00422 for(int i = 0; i < 10; i++){
00423 log << MSG::ALWAYS << m_sel_number[i] << " ";
00424 }
00425 log << MSG::ALWAYS << endreq;
00426 log << MSG::ALWAYS << "==================================" << endreq;
00427 return StatusCode::SUCCESS;
00428 }
00429