00001 #include <fstream> 00002 00003 #include <math.h> 00004 //#include <string.h> 00005 00006 #include "TruthExamples/HistSample.h" 00007 // #include "GeneratorObject/McEventCollection.h" 00008 00009 #include "GaudiKernel/MsgStream.h" 00010 #include "GaudiKernel/AlgFactory.h" 00011 00012 #include "GaudiKernel/SmartDataPtr.h" 00013 #include "GaudiKernel/DataSvc.h" 00014 00015 #include "GaudiKernel/IHistogramSvc.h" 00016 00017 #include "AIDA/IHistogram1D.h" 00018 #include "AIDA/IHistogram2D.h" 00019 00020 #include "GaudiKernel/PropertyMgr.h" 00021 #include "GaudiKernel/INTupleSvc.h" 00022 #include "GaudiKernel/ISvcLocator.h" 00023 #include "GaudiKernel/IDataProviderSvc.h" 00024 #include "GaudiKernel/MsgStream.h" 00025 #include "GaudiKernel/ObjectList.h" 00026 00027 // #include "StoreGate/StoreGateSvc.h" 00028 00029 #include "HepMC/GenEvent.h" 00030 #include "HepMC/GenParticle.h" 00031 #include "HepMC/GenVertex.h" 00032 // #include "CLHEP/HepMC/GenEvent.h" 00033 // #include "CLHEP/HepMC/GenParticle.h" 00034 // #include "CLHEP/HepMC/GenVertex.h" 00035 // #include "CLHEP/HepMC/ParticleDataTableConfig.h" 00036 00037 // #include "PartPropSvc/PartPropSvc.h" 00038 #include "HepPDT/ParticleData.hh" 00039 00040 #include "GeneratorModule/GeneratorName.h" 00041 00042 # 00043 //static const AlgFactory<HistSample> Factory; 00044 //const IAlgFactory& HistSampleFactory = Factory; 00045 00046 HistSample::HistSample(const std::string& name, ISvcLocator* pSvcLocator) : 00047 Algorithm(name, pSvcLocator), 00048 m_hgenerated(0), m_hfinal(0), m_ncharged(0), 00049 m_hChargedPt(0), m_hChargedEta(0), 00050 m_hZPtall(0), m_hZPt(0), m_hZPte(0), m_hZPtm(0), m_hZPtt(0), 00051 m_massZall(0), m_massZ(0), m_massZe(0), m_massZm(0), m_massZt(0), 00052 m_hPtPaire(0), m_hPtPairm(0), m_hPtPairt(0), 00053 m_massPaire(0), m_massPairm(0), m_massPairt(0), 00054 m_rapidity(0), m_pseudorapidity(0), m_hpte() 00055 { 00056 //Declare the algorithm's properties 00057 declareProperty("HistogramFlag", m_produceHistogram = true ); 00058 // declareProperty("HistogramFlag", m_produceHistogram = false ); 00059 00060 } 00061 00062 StatusCode HistSample::initialize(){ 00063 00064 StatusCode result = StatusCode::SUCCESS; 00065 MsgStream msglog(messageService(), name()); 00066 msglog << MSG::INFO << ">>> HistSample from Initialize" << endreq; 00067 // cout << "----- HistSample World From initialize" << endl; 00068 00069 /* 00070 StatusCode sc = service("StoreGateSvc", m_sgSvc); 00071 if (sc.isFailure()) { 00072 msglog << MSG::ERROR << "Could not find StoreGateSvc" << endreq; 00073 return sc; 00074 } 00075 */ 00076 00077 /* 00078 // Get the Particle Properties Service 00079 IPartPropSvc* p_PartPropSvc; 00080 static const bool CREATEIFNOTTHERE(true); 00081 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE); 00082 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) { 00083 msglog << MSG::ERROR << " Could not initialize Particle Properties Service" << endreq; 00084 return PartPropStatus; 00085 } 00086 00087 m_particleTable = p_PartPropSvc->PDT(); 00088 */ 00089 00090 00091 // Register the histograms 00092 // m_hgenerated = histoSvc()->book("/stat/1Dhist/1","Generated",100,0,1000); 00093 // if( 0 == m_hgenerated ) { 00094 // msglog << MSG::ERROR << "Cannot register histo Generated" << endreq; 00095 // } 00096 00097 m_hgenerated = histoSvc()->book("/stat/1Dhist/1","Generated",100,0,1200); 00098 if (0 == m_hgenerated) { 00099 msglog << MSG::ERROR << " ERROR booking histogram" << endreq; 00100 result = StatusCode::FAILURE; 00101 } 00102 m_hfinal = histoSvc()->book("/stat/1Dhist/2","Final state",100,0,800); 00103 m_ncharged = histoSvc()->book("/stat/1Dhist/3","number of charged tracks",100,0,500); 00104 m_hChargedPt = histoSvc()->book("/stat/1Dhist/4","Pt of charged tracks",100,0,100); 00105 m_hChargedEta = histoSvc()->book("/stat/1Dhist/5","Pseudorapidity of charged tracks",100,-12,12); 00106 m_hZPtall = histoSvc()->book("/stat/1Dhist/6","ZPt, all",100,0,200); 00107 m_hZPt = histoSvc()->book("/stat/1Dhist/7","ZPt, charged leptons",100,0,200); 00108 m_hZPte = histoSvc()->book("/stat/1Dhist/8","ZPt electrons",100,0,200); 00109 m_hZPtm = histoSvc()->book("/stat/1Dhist/9","ZPt muons",100,0,200); 00110 m_hZPtt = histoSvc()->book("/stat/1Dhist/10","ZPt taus",100,0,200); 00111 m_massZall = histoSvc()->book("/stat/1Dhist/11","Z mass, all",40,70,110); 00112 m_massZ = histoSvc()->book("/stat/1Dhist/12","Z mass, charged leptons",40,70,110); 00113 m_massZe = histoSvc()->book("/stat/1Dhist/13","Z mass, electrons",40,70,110); 00114 m_massZm = histoSvc()->book("/stat/1Dhist/14","Z mass, muons",40,70,110); 00115 m_massZt = histoSvc()->book("/stat/1Dhist/15","Z mass, taus",40,70,110); 00116 m_hPtPaire = histoSvc()->book("/stat/1Dhist/16","Pt electron pairs",100,0,200); 00117 m_hPtPairm = histoSvc()->book("/stat/1Dhist/17","Pt muon pairs",100,0,200); 00118 m_hPtPairt = histoSvc()->book("/stat/1Dhist/18","Pt tau pairs",100,0,200); 00119 m_massPaire = histoSvc()->book("/stat/1Dhist/19","mass, electron pairs",40,70,110); 00120 m_massPairm = histoSvc()->book("/stat/1Dhist/20","mass, muon pairs",40,70,110); 00121 m_massPairt = histoSvc()->book("/stat/1Dhist/21","mass, tau pairs",40,70,110); 00122 m_rapidity = histoSvc()->book("/stat/1Dhist/22","Rapidity of charged tracks",100,-12,12); 00123 m_pseudorapidity = histoSvc()->book("/stat/1Dhist/23","Pseudorapidity of charged tracks",100,-12,12); 00124 m_hpte = histoSvc()->book("/stat/1Dhist/24","electon pt",100,0,100); 00125 00126 // Initialization terminated 00127 // return StatusCode::SUCCESS; 00128 return result; 00129 } 00130 StatusCode HistSample::execute() { 00131 00132 00133 MsgStream msglog(messageService(), name()); 00134 msglog << MSG::INFO << ">>> HistSample from execute" << endreq; 00135 // cout << "----- HistSample World From execute" << endl; 00136 00137 // Read Data from Transient Store 00138 /* 00139 const McEventCollection* mcCollptr; 00140 if ( m_sgSvc->retrieve(mcCollptr).isFailure() ) { 00141 msglog << MSG::ERROR << "Could not retrieve McEventCollection" 00142 << endreq; 00143 return StatusCode::FAILURE; 00144 } 00145 00146 // cout << " ---- Begin Iterating Over McEventCollection --- " << endl; 00147 McEventCollection::const_iterator it; 00148 for(it=mcCollptr->begin(); it!=mcCollptr->end(); it++) { 00149 // cout << " Generator: " << (*it)->generatorName() << endl; 00151 // (*it)->pGenEvt()->print(); 00152 00153 // fix the STDHEP mess for the Z status 00154 int properStatus = 2; 00155 // switch ( (*it)->generatorName()[0] ) { 00156 // case 'P': 00157 // properStatus = 2; 00158 // break; 00159 // case 'I': 00160 // properStatus = 12; 00161 // break; 00162 // case 'H': 00163 // properStatus = 195; 00164 // break; 00165 // default: 00166 // properStatus = 2; 00167 // break; 00168 // } 00169 HepMC::GenEvent* genEvt = (*it); 00170 int g_id = genEvt->signal_process_id(); 00171 switch ( first_generator(g_id) ) { 00172 case PYTHIA: 00173 properStatus = 2; 00174 break; 00175 case ISAJET: 00176 properStatus = 12; 00177 break; 00178 case HERWIG: 00179 properStatus = 195; 00180 break; 00181 default: 00182 properStatus = 2; 00183 break; 00184 } 00185 00186 00187 int number_particles = 0; 00188 int final_state = 0; 00189 int number_charged = 0; 00190 // Iterate over MC particles 00191 for(HepMC::GenEvent::particle_iterator pitr=genEvt->particles_begin(); 00192 pitr!=genEvt->particles_end(); ++pitr ){ 00193 ++number_particles; 00194 // cout << "particle " << ((*pitr)->pdg_id()) << " status " << ((*pitr)->status()) << endl; 00195 00196 // Z decays 00197 if( ((*pitr)->pdg_id()) == 23 && ((*pitr)->status()) == properStatus ){ 00198 HepMC::GenVertex::particle_iterator firstChild = 00199 (*pitr)->end_vertex()->particles_begin(HepMC::children); 00200 HepMC::GenVertex::particle_iterator endChild = 00201 (*pitr)->end_vertex()->particles_end(HepMC::children); 00202 // plot Pt and mass of the Z (generator values) 00203 if( ((*firstChild)->pdg_id()) != 23 ){ 00204 double ZPt = (*pitr)->momentum().perp(); 00205 m_hZPtall->fill( ZPt, 1.); 00206 double Zmass = (*pitr)->momentum().m(); 00207 m_massZall->fill(Zmass, 1.); 00208 } 00209 // Z decays to l+l- 00210 double sumx = 0.0; 00211 double sumy = 0.0; 00212 double sumz = 0.0; 00213 double sume = 0.0; 00214 int nelds = 0; 00215 int nmuds = 0; 00216 int ntauds = 0; 00217 int Zchild = 0; 00218 HepMC::GenVertex::particle_iterator thisChild = firstChild; 00219 for(; thisChild != endChild++; ++thisChild){ 00220 if( ((*thisChild)->pdg_id()) != 23 ){ 00221 // cout << "Zchild id " << (*thisChild)->pdg_id() << endl; 00222 ++Zchild; 00223 if( abs((*thisChild)->pdg_id()) == 11 || abs((*thisChild)->pdg_id()) == 13 || 00224 abs((*thisChild)->pdg_id()) == 15 ){ 00225 if( abs((*thisChild)->pdg_id()) == 11 )++nelds; 00226 if( abs((*thisChild)->pdg_id()) == 13 )++nmuds; 00227 if( abs((*thisChild)->pdg_id()) == 15 )++ntauds; 00228 sumx += (*thisChild)->momentum().x(); 00229 sumy += (*thisChild)->momentum().y(); 00230 sumz += (*thisChild)->momentum().z(); 00231 sume += (*thisChild)->momentum().e(); 00232 } 00233 } 00234 } 00235 // cout << "Zchildren " << Zchild << " nelds " << nelds << " nmuds " << nmuds << " ntauds " << ntauds << endl; 00236 if( Zchild != 0 ){ 00237 double PtZ2 = sumx*sumx + sumy*sumy; 00238 double PtZ = 0.; 00239 if(PtZ2 > 0.) PtZ = sqrt(PtZ2); 00240 if(PtZ != 0.)m_hZPt->fill( PtZ, 1.); 00241 double massZ2 = sume*sume - sumx*sumx - sumy*sumy - sumz*sumz; 00242 double massZ = 0.; 00243 if(massZ2 > 0.) massZ = sqrt(massZ2); 00244 if(massZ != 0.)m_massZ->fill( massZ, 1.); 00245 if(nelds == 2){ 00246 m_hZPte->fill( PtZ, 1.); 00247 m_massZe->fill( massZ, 1.); 00248 } 00249 if(nmuds == 2){ 00250 m_hZPtm->fill( PtZ, 1.); 00251 m_massZm->fill( massZ, 1.); 00252 } 00253 if(ntauds == 2){ 00254 m_hZPtt->fill( PtZ, 1.); 00255 m_massZt->fill( massZ, 1.); 00256 } 00257 } 00258 } 00259 // all decays to l+l- pairs in the event 00260 if( (*pitr)->end_vertex() ){ 00261 HepMC::GenVertex::particle_iterator fstChild = 00262 (*pitr)->end_vertex()->particles_begin(HepMC::children); 00263 HepMC::GenVertex::particle_iterator lstChild = 00264 (*pitr)->end_vertex()->particles_end(HepMC::children); 00265 double sumpx = 0.0; 00266 double sumpy = 0.0; 00267 double sumpz = 0.0; 00268 double sumpe = 0.0; 00269 int nel = 0; 00270 int nmu = 0; 00271 int ntau = 0; 00272 HepMC::GenVertex::particle_iterator aChild = fstChild; 00273 for(; aChild != lstChild++; ++aChild){ 00274 if( abs((*aChild)->pdg_id()) == 11 || abs((*aChild)->pdg_id()) == 13 || 00275 abs((*aChild)->pdg_id()) == 15 ){ 00276 if( abs((*aChild)->pdg_id()) == 11 )++nel; 00277 if( abs((*aChild)->pdg_id()) == 13 )++nmu; 00278 if( abs((*aChild)->pdg_id()) == 15 )++ntau; 00279 sumpx += (*aChild)->momentum().x(); 00280 sumpy += (*aChild)->momentum().y(); 00281 sumpz += (*aChild)->momentum().z(); 00282 sumpe += (*aChild)->momentum().e(); 00283 } 00284 } 00285 if(nel == 2 || nmu == 2 || ntau == 2 ){ 00286 double PtPair2 = sumpx*sumpx + sumpy*sumpy; 00287 double PtPair = 0.; 00288 if(PtPair2 > 0.) PtPair = sqrt(PtPair2); 00289 double massPair2 = sumpe*sumpe - sumpx*sumpx - sumpy*sumpy - sumpz*sumpz; 00290 double massPair = 0.; 00291 if(massPair2 > 0.) massPair = sqrt(massPair2); 00292 if(nel == 2){ 00293 m_hPtPaire->fill( PtPair, 1.); 00294 m_massPaire->fill( massPair, 1.); 00295 } 00296 if(nmu == 2){ 00297 m_hPtPairm->fill( PtPair, 1.); 00298 m_massPairm->fill( massPair, 1.); 00299 } 00300 if(ntau == 2){ 00301 m_hPtPairt->fill( PtPair, 1.); 00302 m_massPairt->fill( massPair, 1.); 00303 } 00304 } 00305 } 00306 // charged tracks 00307 double c = 0.; 00308 HepPDT::ParticleData* ap = m_particleTable->particle( abs( (*pitr)->pdg_id() ) ); 00309 00310 if(!ap){ 00311 // std::cout << "ChargeService WARNING: abs(id) " 00312 // << abs((*pitr)->pdg_id()) 00313 // << " is not in particle data table" << std::endl; 00314 } 00315 else 00316 { 00317 c = ap->charge(); 00318 } 00319 00320 if( ((*pitr)->status() == 1) && ( ! (*pitr)->end_vertex()) ) { 00321 ++final_state; 00322 if( ((*pitr)->pdg_id()) == 11 ){ 00323 double ePt = (*pitr)->momentum().perp(); 00324 m_hpte->fill( ePt, 1.); 00325 } 00326 if(c!=0.){ 00327 ++number_charged; 00328 double chargedPt = (*pitr)->momentum().perp(); 00329 double chargedEta = (*pitr)->momentum().pseudoRapidity(); 00330 m_hChargedPt->fill( chargedPt, 1.); 00331 m_hChargedEta->fill( chargedEta, 1.); 00332 00333 double px = (*pitr)->momentum().x(); 00334 double py = (*pitr)->momentum().y(); 00335 double pz = (*pitr)->momentum().z(); 00336 double pe = (*pitr)->momentum().e(); 00337 double pp2 = px*px + py*py + pz*pz; 00338 double pp3 = 0.; 00339 if(pp2 > 0.) pp3 = sqrt(pp2); 00340 double y = -999.; 00341 if( (pe-pz) != 0. && (pe+pz)/(pe-pz) > 0. ) y = (1./2.)*log((pe+pz)/(pe-pz)); 00342 double eta = -999.; 00343 if( (pp3-pz) != 0. && (pp3+pz)/(pp3-pz) > 0. ) eta = (1./2.)*log((pp3+pz)/(pp3-pz)); 00344 m_rapidity->fill( y, 1.); 00345 m_pseudorapidity->fill( eta, 1.); 00346 } 00347 } 00348 } 00349 m_hgenerated->fill( number_particles, 1.); 00350 m_hfinal->fill( final_state, 1.); 00351 m_ncharged->fill( number_charged, 1.); 00352 00353 } 00354 */ 00355 // End of execution for each event 00356 return StatusCode::SUCCESS; 00357 } 00358 StatusCode HistSample::finalize() { 00359 00360 MsgStream msglog(messageService(), name()); 00361 msglog << MSG::INFO << ">>> HistSample from finalize" << endreq; 00362 // cout << "----- HistSample World From finalize" << endl; 00363 00364 // End of finalization step 00365 return StatusCode::SUCCESS; 00366 } 00367 00368
1.4.7