00001 #include "EventNavigator/NavigationTestAlg.h"
00002
00003 #include "EventNavigator/EventNavigator.h"
00004
00005 #include <cmath>
00006 #include "GaudiKernel/MsgStream.h"
00007 #include "GaudiKernel/SmartDataPtr.h"
00008
00009
00010
00011 #include "McTruth/McParticle.h"
00012 #include "McTruth/MdcMcHit.h"
00013
00014
00015 #include "MdcRecEvent/RecMdcTrack.h"
00016 #include "MdcRecEvent/RecMdcHit.h"
00017
00018
00019 #include "EmcRecEventModel/RecEmcShower.h"
00020
00021 #include "TH1F.h"
00022 #include "TH2F.h"
00023 #include "TFile.h"
00024
00025 using namespace std;
00026 using namespace Event;
00027
00028 NavigationTestAlg::NavigationTestAlg(const std::string& name, ISvcLocator* pSvcLocator) : Algorithm(name, pSvcLocator)
00029 {
00030 declareProperty("OutputFileName", m_fout="histo.root");
00031 }
00032
00033
00034
00035 StatusCode NavigationTestAlg::initialize(){
00036 MsgStream log(msgSvc(), name());
00037
00038 log << MSG::INFO << " Creating histograms" << endreq;
00039
00040
00041 m_histo.push_back(new TH1F("deltap","#delta p, true momentum - rec momentum, MeV/c",600,-100,500));
00042
00043 m_histo.push_back(new TH1F("n_parents","Number of parents for reconstructed track",20,0,20));
00044 m_histo.push_back(new TH1F("n_tracks","Number of tracks for one MC particle",20,0,20));
00045
00046 m_histo.push_back(new TH1F("hit_parent","Number of parent MC hits for reconstructed hit",20,-10,10));
00047 m_histo.push_back(new TH1F("hit_track","Number of reconstructed hits for MC hit",20,-10,10));
00048
00049 m_histo.push_back(new TH1F("true_p","true momentum, MeV/c",1500,0,1500));
00050 m_histo.push_back(new TH1F("rec_p","reconstructed momentum, MeV/c",1500,0,1500));
00051
00052 m_histo.push_back(new TH1F("n_emcparents","# of parents for reconstructed shower",20,0,20));
00053 m_histo.push_back(new TH1F("n_showers","# of showers for one true particle",20,0,20));
00054
00055 m_histo.push_back(new TH1F("pdg_miss","Pdg of particles with # trk==0 and #showers > 0",6000,-3000,3000));
00056
00057 m_histo.push_back(new TH1F("dE_true","True energy - shower energy",40,-2000,2000));
00058 m_histo.push_back(new TH1F("dE_rec","shower energy - sqrt(Prec**2+m_true**2)",40,-2000,2000));
00059
00060 m_histo.push_back(new TH1F("E_true","True energy, MeV",100,0,1000));
00061 m_histo.push_back(new TH1F("E_rec","Rec energy, MeV",100,0,1000));
00062
00063 m_histo2.push_back(new TH2F("deltap_p","#delta p/p vs. p, MeV",100,0,1000,20,0,2));
00064
00065 return StatusCode::SUCCESS;
00066 }
00067
00068
00069
00070 StatusCode NavigationTestAlg::execute() {
00071 MsgStream log(msgSvc(), name());
00072
00073
00074 SmartDataPtr<EventNavigator> navigator (eventSvc(),"/Event/Navigator");
00075 if( ! navigator )
00076 {
00077 log << MSG::ERROR << " Unable to retrieve EventNavigator" << endreq;
00078 return StatusCode::FAILURE;
00079 }
00080
00081 log << MSG::INFO << "EventNavigator object" << endreq;
00082 navigator->Print();
00083
00084
00085
00086 log << MSG::INFO << "=======================================================================" << endreq;
00087 log << MSG::INFO << "MC Particles ===============================================================" << endreq;
00088
00089
00090 SmartDataPtr<McParticleCol> mcParticles(eventSvc(),"/Event/MC/McParticleCol");
00091 if( ! mcParticles )
00092 {
00093 log << MSG::ERROR << " Unable to retrieve McParticleCol" << endreq;
00094 return StatusCode::FAILURE;
00095 }
00096
00097
00098 for( McParticleCol::iterator it = mcParticles->begin(); it != mcParticles->end(); it++ )
00099 {
00100
00101 int pdg_code = (*it)->particleProperty();
00102
00103
00104 double true_mom = (*it)->initialFourMomentum().vect().mag();
00105
00106
00107 m_histo[5]->Fill(true_mom);
00108
00109 log<<MSG::INFO<<"Retrieved McParticle # "<<(*it)->trackIndex()
00110 << " PDG " << pdg_code << " of true momentum "
00111 << true_mom <<" GeV/c " <<endreq;
00112
00113
00114 RecMdcTrackVector tracks = navigator->getMdcTracks(*it);
00115 RecMdcKalTrackVector ktracks = navigator->getMdcKalTracks(*it);
00116
00117
00118 m_histo[2]->Fill(tracks.size());
00119
00120 log << MSG::INFO << " Found " << tracks.size() << " tracks and " << ktracks.size() << " Kalman tracks" << endreq;
00121
00122
00123 for(unsigned int i=0; i<ktracks.size(); i++)
00124 {
00125
00126 double momentum = ktracks[i]->p();
00127
00128 log << MSG::INFO << "\t Track # " << i
00129 << " id = " << ktracks[i]->trackId()
00130 << " Pid " << ktracks[i]->getPidType()
00131 << " Ptot " << momentum << " GeV/c" << endreq;
00132
00133
00134 m_histo[0]->Fill(true_mom-momentum);
00135
00136
00137 m_histo2[0]->Fill(true_mom, fabs(true_mom-momentum)/true_mom);
00138 }
00139
00140
00141 RecEmcShowerVector showers = navigator->getEmcRecShowers(*it);
00142
00143 m_histo[8]->Fill(showers.size());
00144
00145 log << MSG::INFO << " Found " << showers.size() << " showers" << endreq;
00146
00147 for(unsigned int i=0; i<showers.size(); i++)
00148 {
00149 double true_energy = (*it)->initialFourMomentum().e();
00150 double rec_energy = showers[i]->energy()*1000;
00151
00152 log << MSG::INFO << "\t Shower # " << i
00153 << " Id = " << showers[i]->getShowerId().get_value()
00154 << " E = " << showers[i]->energy()*1000 << " MeV " << endreq;
00155
00156 m_histo[12]->Fill(true_energy);
00157 m_histo[13]->Fill(rec_energy);
00158 }
00159 }
00160
00161 log << MSG::INFO << "MDC Tracks ==============================================================" << endreq;
00162
00163
00164 SmartDataPtr<RecMdcKalTrackCol> mdcKalTracks(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
00165 if( ! mdcKalTracks )
00166 {
00167 log << MSG::ERROR << " Unable to retrieve MdcKalTrackCol" << endreq;
00168 return StatusCode::SUCCESS;
00169 }
00170
00171 for( RecMdcKalTrackCol::iterator it = mdcKalTracks->begin(); it != mdcKalTracks->end(); it++ )
00172 {
00173 McParticleVector particles = navigator->getMcParticles(*it);
00174
00175 log << MSG::INFO << "Retrieved " << particles.size() << " McParticles for for MdcKalTrack # "
00176 << (*it)->trackId() << " of reconstructed momentum " << (*it)->p() << " GeV/c (PID="
00177 << (*it)->getPidType() << ")" << endreq;
00178
00179
00180 m_histo[6]->Fill((*it)->p());
00181
00182
00183 m_histo[1]->Fill(particles.size());
00184
00185
00186 for(unsigned int i=0; i<particles.size(); i++)
00187 {
00188 int pdg_code = particles[i]->particleProperty();
00189
00190
00191 double true_mom = particles[i]->initialFourMomentum().vect().mag();
00192
00193
00194 int relevance = navigator->getMcParticleRelevance(*it, particles[i]);
00195
00196
00197 log << MSG::INFO << "\t" << pdg_code << " momentum "
00198 << true_mom << " relevance " << relevance << endreq;
00199 }
00200 }
00201
00202
00203 log << MSG::INFO << "EMC Showers ==============================================================" << endreq;
00204 SmartDataPtr<RecEmcShowerCol> emcShowers(eventSvc(),"/Event/Recon/RecEmcShowerCol");
00205 if( ! emcShowers )
00206 {
00207 log << MSG::ERROR << " Unable to retrieve EmcRecShowerCol" << endreq;
00208 return StatusCode::SUCCESS;
00209 }
00210
00211 int ij = 0;
00212 for( RecEmcShowerCol::iterator it = emcShowers->begin(); it != emcShowers->end(); it++ )
00213 {
00214 McParticleVector particles = navigator->getMcParticles(*it);
00215 log << MSG::INFO << "Retrieved McParticles for EmcShower # " << ij++
00216 << " Id " << (*it)->getShowerId().get_value()
00217 << " Energy = " << (*it)->energy() << endreq;
00218
00219 for(unsigned int i=0; i<particles.size(); i++)
00220 {
00221 int pdg_code = particles[i]->particleProperty();
00222
00223
00224 double true_e = particles[i]->initialFourMomentum().e();
00225
00226
00227 int relevance = navigator->getMcParticleRelevance(*it, particles[i]);
00228
00229 log << "\t Particle " << i << " PDG " << pdg_code << " E=" << true_e
00230 << " Relevance " << relevance << endreq;
00231 }
00232
00233
00234 m_histo[7]->Fill(particles.size());
00235 }
00236
00237 log << MSG::INFO << "=======================================================================" << endreq;
00238
00239 return StatusCode::SUCCESS;
00240 }
00241
00242
00243
00244 StatusCode NavigationTestAlg::finalize() {
00245
00246
00247 TFile of(m_fout.c_str(),"RECREATE");
00248 of.cd();
00249 for(vector<TH1F*>::iterator it=m_histo.begin(); it!=m_histo.end(); it++)
00250 (*it)->Write();
00251 for(vector<TH2F*>::iterator it=m_histo2.begin(); it!=m_histo2.end(); it++)
00252 (*it)->Write();
00253 of.Close();
00254
00255 return StatusCode::SUCCESS;
00256 }