00001 #include "GaudiKernel/MsgStream.h"
00002 #include "GaudiKernel/AlgFactory.h"
00003 #include "GaudiKernel/ISvcLocator.h"
00004 #include "GaudiKernel/SmartDataPtr.h"
00005 #include "GaudiKernel/IDataProviderSvc.h"
00006 #include "GaudiKernel/PropertyMgr.h"
00007 #include "EventModel/EventModel.h"
00008 #include "EventModel/Event.h"
00009 #include "EventModel/EventHeader.h"
00010 #include "EvtRecEvent/EvtRecEvent.h"
00011 #include "EvtRecEvent/EvtRecTrack.h"
00012 #include "EvtRecEvent/EvtRecPi0.h"
00013 #include "GaudiKernel/INTupleSvc.h"
00014 #include "GaudiKernel/NTuple.h"
00015 #include "GaudiKernel/Bootstrap.h"
00016 #include "GaudiKernel/IHistogramSvc.h"
00017 #include "CLHEP/Vector/ThreeVector.h"
00018 #include "CLHEP/Vector/LorentzVector.h"
00019 #include "CLHEP/Vector/TwoVector.h"
00020 #include "CLHEP/Geometry/Point3D.h"
00021 #include "MdcRecEvent/RecMdcKalTrack.h"
00022
00023 #include "BFieldCorr/BFieldCorr.h"
00024
00025 #include "CLHEP/Matrix/Vector.h"
00026 #include "CLHEP/Vector/LorentzVector.h"
00027 #include "CLHEP/Vector/ThreeVector.h"
00028 #include "CLHEP/Matrix/SymMatrix.h"
00029 #include "CLHEP/Matrix/Matrix.h"
00030 using CLHEP::HepVector;
00031 using CLHEP::HepLorentzVector;
00032 using CLHEP::Hep3Vector;
00033 using CLHEP::HepMatrix;
00034 using CLHEP::HepSymMatrix;
00035
00036 #include "CLHEP/Geometry/Point3D.h"
00037 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00038 typedef HepGeom::Point3D<double> HepPoint3D;
00039 #endif
00040
00041 #include <cmath>
00042 #include <vector>
00043 using namespace std;
00044
00045 BFieldCorr::BFieldCorr(const std::string& name, ISvcLocator* pSvcLocator):Algorithm(name,pSvcLocator)
00046 {
00047 declareProperty("CorrectionFactor", m_factor = 1.0000);
00048 }
00049
00050 StatusCode BFieldCorr::initialize()
00051 {
00052 m_Ew = HepSymMatrix(5, 0);
00053 m_Ew[0][0] = 1.0;
00054 m_Ew[1][1] = 1.0;
00055 m_Ew[2][2] = fabs(m_factor) < 1e-6 ? 1.0 : 1.0 / m_factor;
00056 m_Ew[3][3] = 1.0;
00057 m_Ew[4][4] = 1.0;
00058 RUN_BEGIN_10 = 11414;
00059 RUN_END_10 = 14604;
00060 RUN_BEGIN_11 = 20448;
00061 RUN_END_11 = 23454;
00062
00063 MsgStream log(msgSvc(), name());
00064 log << MSG::INFO << "in MagCorr::initialize()" <<endmsg;
00065 return StatusCode::SUCCESS;
00066 }
00067
00068 StatusCode BFieldCorr::execute()
00069 {
00070 MsgStream log(msgSvc(), name());
00071 log << MSG::INFO << "in BFieldCorr::execute()" << endreq;
00072 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00073 int run_no = eventHeader->runNumber();
00074 if (run_no < 0) return StatusCode::SUCCESS;
00075
00076 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00077 if (!evtRecEvent)
00078 {
00079 log << MSG::FATAL << "Could not find EvtRecEvent" << endreq;
00080 return StatusCode::SUCCESS;
00081 }
00082
00083 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00084 if (!evtRecTrkCol)
00085 {
00086 log << MSG::FATAL << "Could not find EvtRecTrackCol" << endreq;
00087 return StatusCode::SUCCESS;
00088 }
00089
00090 if ((unsigned int)evtRecEvent->totalTracks() > evtRecTrkCol->size()) return StatusCode::SUCCESS;
00091
00092
00093
00094 double factor;
00095 if (fabs(m_factor - 1.000) < 1e-6)
00096 {
00097 if (run_no >= RUN_BEGIN_10 && run_no <= RUN_END_10)
00098 factor = 1.0004;
00099 else if (run_no >= RUN_BEGIN_11 && run_no <= RUN_END_11)
00100 factor = 1.0002;
00101 else
00102 factor = 1.0000;
00103 }
00104 else
00105 {
00106 factor = m_factor;
00107 }
00108 m_Ew[2][2] = fabs(factor) < 1e-6 ? 1.0 : 1.0 / factor;
00109
00110 for (int i = 0; i < evtRecEvent->totalCharged(); i++)
00111 {
00112 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
00113 if (!(*itTrk)->isMdcTrackValid()) continue;
00114 if (!(*itTrk)->isMdcKalTrackValid()) continue;
00115
00116
00117 RecMdcTrack *mdc_trk = (*itTrk)->mdcTrack();
00118 mdc_trk->setHelix( m_Ew * mdc_trk->helix() );
00119 mdc_trk->setError( mdc_trk->err().similarity(m_Ew) );
00120 mdc_trk->setP( factor * mdc_trk->p() );
00121 mdc_trk->setPxy( factor * mdc_trk->pxy() );
00122 mdc_trk->setPx( factor * mdc_trk->px() );
00123 mdc_trk->setPy( factor * mdc_trk->py() );
00124 mdc_trk->setPz( factor * mdc_trk->pz() );
00125
00126
00127 RecMdcKalTrack *kal_trk = (*itTrk)->mdcKalTrack();
00128 RecMdcKalTrack::PidType trk_type[] = {
00129 RecMdcKalTrack::electron,
00130 RecMdcKalTrack::muon,
00131 RecMdcKalTrack::pion,
00132 RecMdcKalTrack::kaon,
00133 RecMdcKalTrack::proton,
00134 };
00135
00136 for (int j = 0; j < 5; j++)
00137 {
00138 kal_trk->setPidType( trk_type[j] );
00139 kal_trk->setZHelix( m_Ew * kal_trk->helix(), j );
00140 kal_trk->setZError( kal_trk->err().similarity(m_Ew), j );
00141 kal_trk->setP( factor * kal_trk->p(), j );
00142 kal_trk->setPxy( factor * kal_trk->pxy(), j );
00143 kal_trk->setPx( factor * kal_trk->px(), j );
00144 kal_trk->setPy( factor * kal_trk->py(), j );
00145 kal_trk->setPz( factor * kal_trk->pz(), j );
00146 }
00147 }
00148 return StatusCode::SUCCESS;
00149 }
00150
00151 StatusCode BFieldCorr::finalize()
00152 {
00153 MsgStream log(msgSvc(), name());
00154 log << MSG::INFO << "in BFieldCorr::finalize()" << endreq;
00155 return StatusCode::SUCCESS;
00156 }