#include <NavigationEx.h>
Public Member Functions | |
| StatusCode | execute () |
| StatusCode | execute () |
| StatusCode | finalize () |
| StatusCode | finalize () |
| StatusCode | initialize () |
| StatusCode | initialize () |
| NavigationEx (const std::string &name, ISvcLocator *pSvcLocator) | |
| NavigationEx (const std::string &name, ISvcLocator *pSvcLocator) | |
Private Member Functions | |
| Hep3Vector | momentum (const RecMdcTrack *trk) |
| Hep3Vector | momentum (const RecMdcTrack *trk) |
Private Attributes | |
| std::string | m_fout |
| std::vector< TH1F * > | m_histo |
| std::vector< TH1F * > | m_histo |
| std::vector< TH2F * > | m_histo2 |
| std::vector< TH2F * > | m_histo2 |
| HepPDT::ParticleDataTable * | m_particleTable |
| HepPDT::ParticleDataTable * | m_particleTable |
|
||||||||||||
|
00030 : Algorithm(name, pSvcLocator)
00031 {
00032 declareProperty("OutputFileName", m_fout="histo.root");
00033 }
|
|
||||||||||||
|
|
|
|
|
|
|
00083 {
00084 MsgStream log(msgSvc(), name());
00085
00086 // Get EventNavigator from the TDS
00087 SmartDataPtr<EventNavigator> navigator (eventSvc(),"/Event/Navigator");
00088 if( ! navigator )
00089 {
00090 log << MSG::ERROR << " Unable to retrieve EventNavigator" << endreq;
00091 return StatusCode::FAILURE;
00092 }
00093
00094 // ======== Analyse Monte-Carlo information using navigator
00095 log << MSG::INFO << "=======================================================================" << endreq;
00096 log << MSG::INFO << "=======================================================================" << endreq;
00097
00098 // Get McParticle collection
00099 SmartDataPtr<McParticleCol> mcParticles(eventSvc(),"/Event/MC/McParticleCol");
00100 if( ! mcParticles )
00101 {
00102 log << MSG::ERROR << " Unable to retrieve McParticleCol" << endreq;
00103 return StatusCode::FAILURE;
00104 }
00105
00106 // For each McParticle ...
00107 for( McParticleCol::iterator it = mcParticles->begin(); it != mcParticles->end(); it++ )
00108 {
00109 // get charge of McParticles
00110 int pdg_code = (*it)->particleProperty();
00111
00112 const HepPDT::ParticleData* particle = m_particleTable->particle(abs(pdg_code));
00113 if( ! particle ){
00114 log<<MSG::WARNING<<"Exotic particle found with PDG code "<<pdg_code<<endreq;
00115 continue;
00116 }
00117 // skip neutrals
00118 // if( particle->charge() == 0 )
00119 // continue;
00120 if( particle->charge() != 0 )
00121 {
00122 // get true momentum value
00123 double true_mom = (*it)->initialFourMomentum().vect().mag();
00124
00125 // fill true momentum
00126 m_histo[5]->Fill(true_mom);
00127
00128 log<<MSG::INFO<<"Retrieved McParticle "<<particle->name()<<" TrackIndex "<<(*it)->trackIndex()
00129 << " of true momentum " << true_mom <<" MeV/c " <<endreq;
00130
00131 // get list of reconstructed tracks, which correspond to this particle, using EventNavigator
00132 RecMdcTrackVector tracks = navigator->getMdcTracks(*it);
00133
00134 // fill the distribution of number of rec. tracks corresponding to single particle
00135 m_histo[2]->Fill(tracks.size());
00136
00137 log << MSG::INFO << "Found " << tracks.size() << " tracks" << endreq;
00138
00139 // for each track
00140 for(unsigned int i=0; i<tracks.size(); i++)
00141 {
00142 Hep3Vector p = momentum(tracks[i]);
00143
00144 log << MSG::INFO << "\t Track # " << i << " id=" << tracks[i]->trackId()
00145 << " momentum " << p.mag() << " MeV/c" << endreq;
00146
00147 // fill difference of true momentum and rec momentum for the _same_ track
00148 m_histo[0]->Fill(true_mom-p.mag());
00149
00150 // fill delta_p/p wrt p for the _same_ track
00151 m_histo2[0]->Fill(true_mom, fabs(true_mom-p.mag())/true_mom);
00152 }
00153 }
00154
00155 // get list of reconstructed tracks, which correspond to this particle, using EventNavigator
00156 RecEmcShowerVector showers = navigator->getEmcRecShowers(*it);
00157
00158 m_histo[8]->Fill(showers.size());
00159 log << MSG::INFO << "Found " << showers.size() << " showers" << endreq;
00160
00161 if ( showers.size()==1 )
00162 {
00163 double true_energy = (*it)->initialFourMomentum().e();
00164 double rec_energy = showers[0]->energy()*1000; // MeV
00165
00166 m_histo[12]->Fill(true_energy);
00167 m_histo[13]->Fill(rec_energy);
00168
00169 RecMdcTrackVector tracks = navigator->getMdcTracks(*it);
00170 if (tracks.size() == 0)
00171 {
00172 m_histo[9]->Fill(pdg_code);
00173 }
00174 }
00175 }
00176
00177 log << MSG::INFO << "=======================================================================" << endreq;
00178 // ======== Analyse reconstructed track information using navigator
00179
00180 SmartDataPtr<RecMdcTrackCol> mdcTracks(eventSvc(),"/Event/Recon/RecMdcTrackCol");
00181 if( ! mdcTracks )
00182 {
00183 log << MSG::ERROR << " Unable to retrieve MdcTrackCol" << endreq;
00184 return StatusCode::SUCCESS;
00185 }
00186
00187 for( RecMdcTrackCol::iterator it = mdcTracks->begin(); it != mdcTracks->end(); it++ )
00188 {
00189 double rec_mom = momentum(*it).mag();
00190
00191 // fill rec momentum
00192 m_histo[6]->Fill(rec_mom);
00193
00194 log << MSG::INFO << "Retrieved McParticles for for MdcTrack # " << (*it)->trackId()
00195 << " of reconstructed momentum " << rec_mom << " MeV/c" << endreq;
00196
00197 McParticleVector particles = navigator->getMcParticles(*it);
00198
00199 // fill the distribution of number of parent particles corresponding to this track
00200 m_histo[1]->Fill(particles.size());
00201
00202 // for each parent particle
00203 for(unsigned int i=0; i<particles.size(); i++)
00204 {
00205 int pdg_code = particles[i]->particleProperty();
00206
00207 const HepPDT::ParticleData* particle = m_particleTable->particle(abs(pdg_code));
00208 if( ! particle )
00209 {
00210 log<<MSG::WARNING<<"Exotic particle found with PDG code "<<pdg_code<<endreq;
00211 continue;
00212 }
00213
00214 // get true momentum value
00215 double true_mom = particles[i]->initialFourMomentum().vect().mag();
00216
00217 // dump particle names and momenta
00218 log << MSG::INFO << "\t" << particle->name() << " momentum " << true_mom << endreq;
00219 }
00220 }
00221
00222 // ======== Analyse reconstructed shower information using navigator
00223 SmartDataPtr<RecEmcShowerCol> emcShowers(eventSvc(),"/Event/Recon/RecEmcShowerCol");
00224 if( ! emcShowers )
00225 {
00226 log << MSG::ERROR << " Unable to retrieve EmcRecShowerCol" << endreq;
00227 return StatusCode::SUCCESS;
00228 }
00229
00230 for( RecEmcShowerCol::iterator it = emcShowers->begin(); it != emcShowers->end(); it++ )
00231 {
00232 McParticleVector particles = navigator->getMcParticles(*it);
00233
00234 // fill the distribution of number of parent particles corresponding to this track
00235 m_histo[7]->Fill(particles.size());
00236 }
00237
00238 SmartDataPtr<EmcMcHitCol> emcMcHits(eventSvc(),"/Event/MC/EmcMcHitCol");
00239 if( ! emcMcHits )
00240 {
00241 log << MSG::ERROR << " Unable to retrieve EmcMcHits" << endreq;
00242 return StatusCode::SUCCESS;
00243 }
00244
00245 for( EmcMcHitCol::iterator it = emcMcHits->begin(); it != emcMcHits->end(); it++ )
00246 {
00247 log << MSG::DEBUG << "HitMap size = " << (*it)->getHitMap().size() << endreq;
00248 }
00249
00250 log << MSG::INFO << "=======================================================================" << endreq;
00251
00252 return StatusCode::SUCCESS;
00253 }
|
|
|
|
|
|
00257 {
00258
00259 // Save histograms to ROOT file
00260 TFile of(m_fout.c_str(),"RECREATE");
00261 of.cd();
00262 for(vector<TH1F*>::iterator it=m_histo.begin(); it!=m_histo.end(); it++)
00263 (*it)->Write();
00264 for(vector<TH2F*>::iterator it=m_histo2.begin(); it!=m_histo2.end(); it++)
00265 (*it)->Write();
00266 of.Close();
00267
00268 return StatusCode::SUCCESS;
00269 }
|
|
|
|
|
|
00037 {
00038 MsgStream log(msgSvc(), name());
00039
00040 // Get the Particle Properties Service
00041 IPartPropSvc* p_PartPropSvc;
00042 static const bool CREATEIFNOTTHERE(true);
00043 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
00044 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
00045 log << MSG::ERROR << " Could not initialize Particle Properties Service" << endreq;
00046 return PartPropStatus;
00047 }
00048
00049 m_particleTable = p_PartPropSvc->PDT();
00050
00051 log << MSG::INFO << " Creating histograms" << endreq;
00052
00053 // create histograms
00054 m_histo.push_back(new TH1F("deltap","#delta p, true momentum - rec momentum, MeV/c",600,-100,500)); //0
00055
00056 m_histo.push_back(new TH1F("n_parents","Number of parents for reconstructed track",20,0,20)); // 1
00057 m_histo.push_back(new TH1F("n_tracks","Number of tracks for one MC particle",20,0,20)); // 2
00058
00059 m_histo.push_back(new TH1F("hit_parent","Number of parent MC hits for reconstructed hit",20,-10,10)); // 3
00060 m_histo.push_back(new TH1F("hit_track","Number of reconstructed hits for MC hit",20,-10,10)); // 4
00061
00062 m_histo.push_back(new TH1F("true_p","true momentum, MeV/c",1500,0,1500)); // 5
00063 m_histo.push_back(new TH1F("rec_p","reconstructed momentum, MeV/c",1500,0,1500)); // 6
00064
00065 m_histo.push_back(new TH1F("n_emcparents","# of parents for reconstructed shower",20,0,20)); // 7
00066 m_histo.push_back(new TH1F("n_showers","# of showers for one true particle",20,0,20)); // 8
00067
00068 m_histo.push_back(new TH1F("pdg_miss","Pdg of particles with # trk==0 and #showers > 0",6000,-3000,3000)); //9
00069
00070 m_histo.push_back(new TH1F("dE_true","True energy - shower energy",40,-2000,2000)); //10
00071 m_histo.push_back(new TH1F("dE_rec","shower energy - sqrt(Prec**2+m_true**2)",40,-2000,2000)); //11
00072
00073 m_histo.push_back(new TH1F("E_true","True energy, MeV",100,0,1000)); //12
00074 m_histo.push_back(new TH1F("E_rec","Rec energy, MeV",100,0,1000)); //13
00075
00076 m_histo2.push_back(new TH2F("deltap_p","#delta p/p vs. p, MeV",100,0,1000,20,0,2));
00077
00078 return StatusCode::SUCCESS;
00079 }
|
|
|
|
|
|
00273 {
00274 // double dr = trk->getDr();
00275 double fi0 = trk->helix(1);
00276 double cpa = trk->helix(2);
00277 double tanl = trk->helix(4);
00278
00279 double pxy = 0;
00280
00281 if(cpa != 0) pxy = 1/fabs(cpa);
00282 else pxy = 0;
00283
00284 double px = pxy * (-sin(fi0));
00285 double py = pxy * cos(fi0);
00286 double pz = pxy * tanl;
00287
00288 Hep3Vector p;
00289 p.set(px*1000, py*1000, pz*1000); // MeV/c
00290 return p;
00291 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1.3.9.1