BOSS_7.0.2: /home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Reconstruction/VeeVertexAlg/VeeVertexAlg-00-00-05/src/LambdaReconstruction.cxx Source File

00001 //
00002 //      Lambda -> p+ pi- Reconstruction
00003 // Anti-Lambda -> p- pi+
00004 //
00005 
00006 #include "GaudiKernel/MsgStream.h"
00007 #include "GaudiKernel/AlgFactory.h"
00008 #include "GaudiKernel/ISvcLocator.h"
00009 #include "GaudiKernel/SmartDataPtr.h"
00010 #include "GaudiKernel/IDataProviderSvc.h"
00011 #include "GaudiKernel/PropertyMgr.h"
00012 
00013 #include "EventModel/Event.h"
00014 #include "EvtRecEvent/EvtRecEvent.h"
00015 #include "EvtRecEvent/EvtRecTrack.h"
00016 #include "EvtRecEvent/EvtRecVeeVertex.h"
00017 #include "EventModel/EventHeader.h"
00018 #include "EventModel/Event.h"
00019 #include "EventModel/EventModel.h"
00020 #include "McTruth/McParticle.h"
00021 #include "ParticleID/ParticleID.h"
00022 #include "VertexFit/VertexFit.h"
00023 #include "VertexFit/SecondVertexFit.h"
00024 #include "VeeVertexAlg/LambdaReconstruction.h"
00025 #include <vector>
00026 
00027 typedef std::vector<int> Vint;
00028 const double mpi =  0.13957;
00029 const double mp = 0.938272;
00030 //*******************************************************************************************
00031 LambdaReconstruction::LambdaReconstruction(const std::string& name, ISvcLocator* pSvcLocator):
00032 Algorithm(name, pSvcLocator)  {
00033   //Declare the properties
00034   
00035   declareProperty("Vz0Cut", m_vzCut = 20.0);
00036   declareProperty("CosThetaCut", m_cosThetaCut=0.93);
00037   declareProperty("ChisqCut", m_chisqCut = 100.0);
00038   /*
00039   declareProperty("PidUseDedx", m_useDedx = true);
00040   declareProperty("PidUseTof1", m_useTof1 = true);
00041   declareProperty("PidUseTof2", m_useTof2 = true);
00042   declareProperty("PidUseTofE", m_useTofE = false);
00043   declareProperty("PidUseTofQ", m_useTofQ = false);
00044   declareProperty("PidUseEmc", m_useEmc = false);
00045   declareProperty("PidProbCut", m_PidProbCut= 0.001);
00046   declareProperty("MassRangeLower", m_massRangeLower = 0.45);
00047   declareProperty("MassRangeHigher", m_massRangeHigher = 0.55);
00048   */
00049 }
00050 
00051 //******************************************************************************************
00052 StatusCode LambdaReconstruction::initialize() {
00053   MsgStream log(msgSvc(), name());
00054   log << MSG::INFO << "in initialize()" <<endreq;
00055   
00056   return StatusCode::SUCCESS;
00057 }
00058 
00059 //********************************************************************************************
00060 StatusCode LambdaReconstruction::finalize() {
00061   MsgStream log(msgSvc(), name());
00062   log << MSG::INFO << "in finalize()" << endreq;
00063 
00064   return StatusCode::SUCCESS;
00065 }
00066 
00067 StatusCode LambdaReconstruction::registerEvtRecVeeVertexCol(
00068     EvtRecVeeVertexCol* aNewEvtRecVeeVertexCol, MsgStream& log) {
00069   StatusCode sc = eventSvc()->registerObject("/Event/EvtRec/EvtRecVeeVertexCol",
00070                                              aNewEvtRecVeeVertexCol);
00071   if (sc != StatusCode::SUCCESS) {
00072     log << MSG::FATAL << "Could not register EvtRecVeeVertexCol in TDS!" << endreq;
00073   }
00074   return sc;
00075 }
00076 
00077 //*********************************************************************************************
00078 StatusCode LambdaReconstruction::execute() {
00079   MsgStream log(msgSvc(), name());
00080   log << MSG::INFO << "in execute()" << endreq;
00081 
00082   StatusCode sc;
00083 
00085   // Read REC data
00087   SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00088   SmartDataPtr<EvtRecEvent> recEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00089   SmartDataPtr<EvtRecTrackCol> recTrackCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00090   log << MSG::DEBUG << "run and event = " << eventHeader->runNumber()
00091                                         << " " << eventHeader->eventNumber() << endreq;
00092   log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00093           << recEvent->totalCharged() << " , "
00094           << recEvent->totalNeutral() << " , "
00095       << recEvent->totalTracks() <<endreq;
00096   int evtNo = eventHeader->eventNumber();
00097   //if (eventHeader->eventNumber() % 1000 == 0) {
00098   //  cout << "Event number = " << evtNo << endl;
00099   //}
00100 
00101   SmartDataPtr<EvtRecVeeVertexCol> veeVertexCol(eventSvc(),
00102                                                 EventModel::EvtRec::EvtRecVeeVertexCol);
00103   if (!veeVertexCol) {
00104     veeVertexCol = new EvtRecVeeVertexCol;
00105     sc = registerEvtRecVeeVertexCol(veeVertexCol, log);
00106     if (sc != StatusCode::SUCCESS) {
00107       return sc;
00108     }
00109   }
00110 
00112   //  No PID : identify positive and negative
00114   Vint icp, icm, iGood;
00115   for (unsigned int i = 0; i < recEvent->totalCharged(); i++) {
00116     EvtRecTrackIterator itTrk = recTrackCol->begin() + i;
00117     if(!(*itTrk)->isMdcTrackValid()) continue;
00118     if(!(*itTrk)->isMdcKalTrackValid()) continue;
00119     RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
00120     if (fabs(cos(mdcTrk->theta())) >= m_cosThetaCut) continue;
00121     if (fabs(mdcTrk->z()) >= m_vzCut) continue;
00122     iGood.push_back(i);
00123     if (mdcTrk->charge() > 0) icp.push_back(i);
00124     if (mdcTrk->charge() < 0) icm.push_back(i);
00125   }
00126   // Cut on good charged tracks 
00127   //if (iGood.size() < 2) return StatusCode::SUCCESS;
00128   if (icp.size() < 1 || icm.size() < 1) return StatusCode::SUCCESS;
00129 
00131   //  Vertex Fit
00133   HepPoint3D vx(0., 0., 0.);
00134   HepSymMatrix Evx(3, 0);
00135   double bx = 1E+6;
00136   double by = 1E+6;
00137   double bz = 1E+6;
00138   Evx[0][0] = bx*bx;
00139   Evx[1][1] = by*by;
00140   Evx[2][2] = bz*bz;
00141   VertexParameter vxpar;
00142   vxpar.setVx(vx);
00143   vxpar.setEvx(Evx);
00144  
00145   VertexFit *vtxfit0 = VertexFit::instance();
00146 
00147   // For Lambda reconstruction
00148   for(unsigned int i1 = 0; i1 < icp.size(); i1++) {
00149     int ip1 = icp[i1];
00150     RecMdcKalTrack *ppKalTrk = (*(recTrackCol->begin()+ip1))->mdcKalTrack();
00151     RecMdcKalTrack::setPidType(RecMdcKalTrack::proton);
00152     WTrackParameter wpptrk(mp, ppKalTrk->helix(), ppKalTrk->err());
00153 
00154     for(unsigned int i2 = 0; i2 < icm.size(); i2++) {
00155       int ip2 = icm[i2];
00156       RecMdcKalTrack *pimKalTrk = (*(recTrackCol->begin()+ip2))->mdcKalTrack();
00157       RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
00158       WTrackParameter wpimtrk(mpi, pimKalTrk->helix(), pimKalTrk->err());
00159    
00160       vtxfit0->init();
00161       vtxfit0->AddTrack(0, wpptrk);
00162       vtxfit0->AddTrack(1, wpimtrk);
00163       vtxfit0->AddVertex(0, vxpar, 0, 1);
00164       if(!(vtxfit0->Fit(0))) continue;
00165       vtxfit0->BuildVirtualParticle(0);
00166       if(vtxfit0->chisq(0) > m_chisqCut) continue;
00167       WTrackParameter wLamb = vtxfit0->wVirtualTrack(0); //FIXME
00168       std::pair<int, int> pair;
00169       pair.first = 5;
00170       pair.second = 3;
00171 
00172       EvtRecTrack* track0 = *(recTrackCol->begin() + ip1);
00173       EvtRecTrack* track1 = *(recTrackCol->begin() + ip2);
00174 
00175       EvtRecVeeVertex* LambdaVertex = new EvtRecVeeVertex;
00176       LambdaVertex->setVertexId(3122);
00177       LambdaVertex->setVertexType(1);
00178       LambdaVertex->setChi2(vtxfit0->chisq(0));
00179       LambdaVertex->setNdof(1);
00180       LambdaVertex->setMass(wLamb.p().m());
00181       LambdaVertex->setW(wLamb.w());
00182       LambdaVertex->setEw(wLamb.Ew());
00183       LambdaVertex->setPair(pair);
00184       LambdaVertex->setNCharge(0);
00185       LambdaVertex->setNTracks(2);
00186       LambdaVertex->addDaughter(track0, 0);
00187       LambdaVertex->addDaughter(track1, 1);
00188       veeVertexCol->push_back(LambdaVertex);
00189   /* 
00190       cout << "============= After Vertex fitting ===========" << endl; 
00191       cout << "Event number = " << evtNo << endl;
00192       cout << "chi2 = " << vtxfit0->chisq(0) << endl;
00193       cout << "mass = " << wLamb.p().m() << endl;
00194       cout << "w = " << wLamb.w() << endl;
00195       cout << "Ew = " << wLamb.Ew() << endl;
00196       cout << "track 0 : Z : " << track0->mdcTrack()->z() << endl;
00197       cout << "track 1 : Z : " << track1->mdcTrack()->z() << endl; */
00198     }
00199   }
00200 
00201   // For Anti-Lambda reconstruction
00202   for(unsigned int i1 = 0; i1 < icp.size(); i1++) {
00203     int ip1 = icp[i1];
00204     RecMdcKalTrack *pipKalTrk = (*(recTrackCol->begin()+ip1))->mdcKalTrack();
00205     RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);
00206     WTrackParameter wpiptrk(mpi, pipKalTrk->helix(), pipKalTrk->err());
00207 
00208     for(unsigned int i2 = 0; i2 < icm.size(); i2++) {
00209       int ip2 = icm[i2];
00210       RecMdcKalTrack *pmKalTrk = (*(recTrackCol->begin()+ip2))->mdcKalTrack();
00211       RecMdcKalTrack::setPidType(RecMdcKalTrack::proton);
00212       WTrackParameter wpmtrk(mp, pmKalTrk->helix(), pmKalTrk->err());
00213 
00214       vtxfit0->init();
00215       vtxfit0->AddTrack(0, wpiptrk);
00216       vtxfit0->AddTrack(1, wpmtrk);
00217       vtxfit0->AddVertex(0, vxpar, 0, 1);
00218       if(!(vtxfit0->Fit(0))) continue;
00219       vtxfit0->BuildVirtualParticle(0);
00220       if(vtxfit0->chisq(0) > m_chisqCut) continue;
00221       WTrackParameter wALamb = vtxfit0->wVirtualTrack(0); //FIXME
00222       std::pair<int, int> pair;
00223       pair.first = 3;
00224       pair.second = 5;
00225 
00226       EvtRecTrack* track0 = *(recTrackCol->begin() + ip1);
00227       EvtRecTrack* track1 = *(recTrackCol->begin() + ip2);
00228 
00229       EvtRecVeeVertex* ALambdaVertex = new EvtRecVeeVertex;
00230       ALambdaVertex->setVertexId(-3122);
00231       ALambdaVertex->setVertexType(1);
00232       ALambdaVertex->setChi2(vtxfit0->chisq(0));
00233       ALambdaVertex->setNdof(1);
00234       ALambdaVertex->setMass(wALamb.p().m());
00235       ALambdaVertex->setW(wALamb.w());
00236       ALambdaVertex->setEw(wALamb.Ew());
00237       ALambdaVertex->setPair(pair);
00238       ALambdaVertex->setNCharge(0);
00239       ALambdaVertex->setNTracks(2);
00240       ALambdaVertex->addDaughter(track0, 0);
00241       ALambdaVertex->addDaughter(track1, 1);
00242       veeVertexCol->push_back(ALambdaVertex);
00243   /* 
00244       cout << "============= After Vertex fitting ===========" << endl; 
00245       cout << "Event number = " << evtNo << endl;
00246       cout << "chi2 = " << vtxfit0->chisq(0) << endl;
00247       cout << "mass = " << wALamb.p().m() << endl;
00248       cout << "w = " << wALamb.w() << endl;
00249       cout << "Ew = " << wALamb.Ew() << endl;
00250       cout << "track 0 : Z : " << track0->mdcTrack()->z() << endl;
00251       cout << "track 1 : Z : " << track1->mdcTrack()->z() << endl; */
00252     }
00253   }
00254 
00255 
00256 
00257   return StatusCode::SUCCESS;
00258 }