diff --git a/Generators/TruthIO/TruthIO/DumpMC.h b/Generators/TruthIO/TruthIO/DumpMC.h new file mode 100644 index 0000000000000000000000000000000000000000..8587f83806880008b5f8a93b385f2f6335fea8ee --- /dev/null +++ b/Generators/TruthIO/TruthIO/DumpMC.h @@ -0,0 +1,25 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TRUTHIO_DUMPMC_H +#define TRUTHIO_DUMPMC_H + +#include "GeneratorModules/GenBase.h" + + +/// Dump MC event record +class DumpMC : public GenBase { +public: + + DumpMC(const std::string& name, ISvcLocator* pSvcLocator); + StatusCode initialize(); + StatusCode execute(); + + std::string m_keyout; + bool m_VerboseOutput; + bool m_DeepCopy; + bool m_EtaPhi; +}; + +#endif diff --git a/Generators/TruthIO/TruthIO/HepMCReadFromFile.h b/Generators/TruthIO/TruthIO/HepMCReadFromFile.h new file mode 100644 index 0000000000000000000000000000000000000000..f938340cb5a54e152f6780f65bd0533c64208c03 --- /dev/null +++ b/Generators/TruthIO/TruthIO/HepMCReadFromFile.h @@ -0,0 +1,27 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TRUTHIO_HEPMCREADFROMFILE_H +#define TRUTHIO_HEPMCREADFROMFILE_H + +#include "GeneratorModules/GenBase.h" +#include "HepMC/IO_GenEvent.h" +#include <memory> + +class HepMCReadFromFile : public GenBase { +public: + + HepMCReadFromFile(const std::string& name, ISvcLocator* pSvcLocator); + StatusCode initialize(); + StatusCode execute(); + +private: + + std::string m_input_file; + int m_event_number; + + std::auto_ptr<HepMC::IO_GenEvent> m_hepmcio; +}; + +#endif diff --git a/Generators/TruthIO/TruthIO/PrintHijingPars.h b/Generators/TruthIO/TruthIO/PrintHijingPars.h new file mode 100644 index 0000000000000000000000000000000000000000..8acbf61943d55ea06f04b74e436201cf11cc6db7 --- /dev/null +++ b/Generators/TruthIO/TruthIO/PrintHijingPars.h @@ -0,0 +1,27 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TRUTHIO_PRINTHIJINGPARS_H +#define TRUTHIO_PRINTHIJINGPARS_H + +#include "GeneratorModules/GenBase.h" + +class PrintHijingPars:public GenBase { +public: + PrintHijingPars(const std::string& name, ISvcLocator* pSvcLocator); + StatusCode initialize(); + StatusCode execute(); + StatusCode finalize(); + +private: + + // Setable Properties:- + std::string m_key; + bool m_VerboseOutput; + int m_firstEvt; + int m_lastEvt; +}; + +#endif + diff --git a/Generators/TruthIO/TruthIO/PrintMC.h b/Generators/TruthIO/TruthIO/PrintMC.h new file mode 100644 index 0000000000000000000000000000000000000000..567ce752f98640d2b4c0d83a8fa8be3edbfddfe5 --- /dev/null +++ b/Generators/TruthIO/TruthIO/PrintMC.h @@ -0,0 +1,31 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TRUTHIO_PRINTMC_H +#define TRUTHIO_PRINTMC_H + +#include "GeneratorModules/GenBase.h" + + +/// Print MC event details for a range of event numbers +class PrintMC : public GenBase { +public: + + PrintMC(const std::string& name, ISvcLocator* pSvcLocator); + StatusCode initialize(); + StatusCode execute(); + +private: + + std::string m_keyout; + bool m_VerboseOutput; + std::string m_printsty; + bool m_vertexinfo; + int m_firstEvt; + int m_lastEvt; + bool m_trustHepMC; + +}; + +#endif diff --git a/Generators/TruthIO/TruthIO/WriteHepMC.h b/Generators/TruthIO/TruthIO/WriteHepMC.h new file mode 100644 index 0000000000000000000000000000000000000000..d1dcdaaffb38a4caaa6dd846073be0dbac95d2d4 --- /dev/null +++ b/Generators/TruthIO/TruthIO/WriteHepMC.h @@ -0,0 +1,28 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TRUTHIO_WRITEHEPMC_H +#define TRUTHIO_WRITEHEPMC_H + +#include "GeneratorModules/GenBase.h" +#include "HepMC/IO_GenEvent.h" +#include <memory> + + +/// Write the MC event record to file in IO_GenEvent text format +class WriteHepMC : public GenBase { +public: + + WriteHepMC(const std::string& name, ISvcLocator* pSvcLocator); + StatusCode initialize(); + StatusCode execute(); + + std::string m_outfile; + int m_precision; + + std::auto_ptr<HepMC::IO_GenEvent> m_hepmcio; + +}; + +#endif diff --git a/Generators/TruthIO/cmt/requirements b/Generators/TruthIO/cmt/requirements new file mode 100644 index 0000000000000000000000000000000000000000..0bdf33a9402ce7358c72998166940bbb640d48df --- /dev/null +++ b/Generators/TruthIO/cmt/requirements @@ -0,0 +1,24 @@ +package TruthIO + +author Andy Buckley <andy.buckley@cern.ch> + +use AtlasPolicy AtlasPolicy-* +use GaudiInterface GaudiInterface-* External +use GeneratorModules GeneratorModules-* Generators +use AtlasHepMC AtlasHepMC-* External + +private +use GeneratorObjects GeneratorObjects-* Generators +use EventInfo EventInfo-* Event +use HepPDT v* LCG_Interfaces +end_private + +private +macro_append HepMC_linkopts " $(HepMC_IO_linkopts) " +end_private + +public +library TruthIO files=" *.cxx " -s=components *.cxx +apply_pattern component_library + +apply_pattern declare_joboptions files="*.py" diff --git a/Generators/TruthIO/share/jobOptions.read_hepmc.py b/Generators/TruthIO/share/jobOptions.read_hepmc.py new file mode 100644 index 0000000000000000000000000000000000000000..f3a6237913851f2849caece3079de6173117fe8b --- /dev/null +++ b/Generators/TruthIO/share/jobOptions.read_hepmc.py @@ -0,0 +1,50 @@ +############################################################### +# +# Job options file +# +#============================================================== +#-------------------------------------------------------------- +# General Application Configuration options +#-------------------------------------------------------------- +import AthenaCommon.AtlasUnixGeneratorJob + +from AthenaCommon.AppMgr import theApp +from AthenaCommon.AppMgr import ServiceMgr + +#-------------------------------------------------------------- +# Private Application Configuration options +#-------------------------------------------------------------- +from AthenaCommon.AlgSequence import AlgSequence +job=AlgSequence() +from ReadHepMCFromFile.ReadHepMCFromFileConf import HepMCReadFromFile +job += HepMCReadFromFile() + +# Set output level threshold (2=DEBUG, 3=INFO, 4=WARNING, 5=ERROR, 6=FATAL ) +ServiceMgr.MessageSvc.OutputLevel = INFO +#-------------------------------------------------------------- +# Event related parameters +#-------------------------------------------------------------- +# Number of events to be processed (default is 10) +theApp.EvtMax = 5 +#-------------------------------------------------------------- +# Algorithms Private Options +#-------------------------------------------------------------- + +re = job.HepMCReadFromFile +re.AsciiFile= "hepmc.out" + +from AthenaPoolCnvSvc.WriteAthenaPool import AthenaPoolOutputStream +Stream1 = AthenaPoolOutputStream( "Stream1" ) +Stream1.ItemList += [ "2101#*", "133273#*" ] +# This line is not really needed. Just to show that PoolSvc is already covered by the ServiceMgr +PoolSvc = ServiceMgr.PoolSvc +Stream1.OutputFile = "McEvent.root" + +#--------------------------------------------------------------- +# Ntuple service output +#--------------------------------------------------------------- +#============================================================== +# +# End of job options file +# +############################################################### diff --git a/Generators/TruthIO/src/DumpMC.cxx b/Generators/TruthIO/src/DumpMC.cxx new file mode 100644 index 0000000000000000000000000000000000000000..e2e4805f836b5f8413bc5758ec5624dc5cbbc86d --- /dev/null +++ b/Generators/TruthIO/src/DumpMC.cxx @@ -0,0 +1,145 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "TruthIO/DumpMC.h" +#include "GaudiKernel/PhysicalConstants.h" +using namespace Gaudi::Units; +using namespace std; + + +DumpMC::DumpMC(const string& name, ISvcLocator* pSvcLocator) + : GenBase(name, pSvcLocator) +{ + declareProperty("McEventOutputKey", m_keyout="GEN_EVENT"); + declareProperty("VerboseOutput", m_VerboseOutput=true); + declareProperty("DeepCopy", m_DeepCopy=false); + declareProperty("EtaPhi", m_EtaPhi=false); +} + + +StatusCode DumpMC::initialize() { + CHECK(GenBase::initialize()); + if (m_mcEventKey == m_keyout && m_DeepCopy) { + ATH_MSG_FATAL("Input and output MC event keys cannot be the same (" << m_mcEventKey << ")"); + return StatusCode::FAILURE; + } + return StatusCode::SUCCESS; +} + + +StatusCode DumpMC::execute() { + if (m_DeepCopy) { + McEventCollection* mcCollptra = new McEventCollection(); + // Fill the new McEventCollection with a copy of the initial HepMC::GenEvent + for (McEventCollection::const_iterator evt = events()->begin(); evt != events()->end(); ++evt) { + mcCollptra->push_back(new HepMC::GenEvent(*(*evt))); + } + // Loop over all events in McEventCollection + for (McEventCollection::iterator evt = mcCollptra->begin(); evt != mcCollptra->end(); ++evt) { + // Loop over the vertices of the event + set<int> Lambdas; + set<int> vtx_to_delete; + for (HepMC::GenEvent::vertex_const_iterator vtx = (*evt)->vertices_begin(); vtx != (*evt)->vertices_end(); ++vtx) { + // Loop over the particles that produced the vertex + if ( vtx_to_delete.find((*vtx)->barcode()) == vtx_to_delete.end() && + (*vtx)->particles_in_const_begin() != (*vtx)->particles_in_const_end() ) { + HepMC::GenVertex::particles_in_const_iterator part = (*vtx)->particles_in_const_begin(); + bool lambda_not_found = true; + const int abspid = abs((*part)->pdg_id()); + do { + if (abspid == 310 || abspid == 3122 || abspid == 3222 || abspid == 3112 || + abspid == 3322 || abspid == 3312 || abspid == 3334 ) { + lambda_not_found = false; + Lambdas.insert((*part)->barcode()); + vtx_to_delete.insert((*vtx)->barcode()); + // In case a lambda was found, store in vertices to be deleted all the vertices + // that the products of Lambda created + for (HepMC::GenVertex::vertex_iterator desc = (*vtx)->vertices_begin(HepMC::descendants); + desc != (*vtx)->vertices_end(HepMC::descendants); ++desc) { + vtx_to_delete.insert((*desc)->barcode()); + } + } + ++part; + } while (part != (*vtx)->particles_in_const_end() && lambda_not_found); + } + } + + // Set Lambda's status to stable + for (set<int>::iterator l = Lambdas.begin(); l != Lambdas.end(); ++l) { + HepMC::GenParticle* lam = (*evt)->barcode_to_particle(*l); + lam->set_status(1); + } + // Delete all Lambda vertices from the event + for (set<int>::iterator v = vtx_to_delete.begin(); v != vtx_to_delete.end(); ++v) { + HepMC::GenVertex* vdel = (*evt)->barcode_to_vertex(*v); + (*evt)->remove_vertex(vdel); + delete vdel; + } + } + + if (evtStore()->record(mcCollptra, m_keyout).isFailure()){ + ATH_MSG_ERROR("Could not record McEventCollection"); + return StatusCode::FAILURE; + } + } + + // Loop over all events in McEventCollection + for (McEventCollection::const_iterator itr = events()->begin(); itr != events()->end(); ++itr) { + //int g_id = (*itr)->signal_process_id(); + //GeneratorName_print(g_id); + HepMC::PdfInfo *pdfinfo = (*itr)->pdf_info(); + if (pdfinfo) { + cout << "PdfInfo: " + << pdfinfo->id1() << ", " + << pdfinfo->id2() << ", " + << pdfinfo->x1() << ", " + << pdfinfo->x2() << ", " + << pdfinfo->scalePDF() << ", " + << pdfinfo->pdf1() << ", " + << pdfinfo->pdf2() << ", " + << pdfinfo->pdf_id1() << ", " + << pdfinfo->pdf_id2() + << endl; + } + + HepMC::HeavyIon *ion = (*itr)->heavy_ion(); + if (ion) { + std::cout << std::endl; + std::cout << "Heavy Ion: " + << ion->Ncoll_hard() <<", " + << ion->Npart_proj() <<" , " + << ion->Npart_targ()<< ", " + << ion->Ncoll()<< ", " + << ion->spectator_neutrons() << ", " + << ion->spectator_protons() << ", " + << ion->N_Nwounded_collisions() << ", " + << ion->Nwounded_N_collisions() << ", " + << ion->Nwounded_Nwounded_collisions() << ", " + << ion->impact_parameter() << ", " + << ion->event_plane_angle() << ", " + << ion->eccentricity() << ", " + << ion->sigma_inel_NN() + << std::endl; + } + + if (m_VerboseOutput) { + if (!m_EtaPhi) { + (*itr)->print(); // standard HepMc dump + } else { // sort particles by rapidity and then dump + // Loop over all particles in the event and build up the grid + const HepMC::GenEvent* genEvt = (*itr); + for (HepMC::GenEvent::particle_const_iterator pitr=genEvt->particles_begin(); pitr!=genEvt->particles_end(); ++pitr) { + // if( (*pitr)->status() == 1 ){// stables only + double rapid =(*pitr)->momentum().pseudoRapidity(); + double phi = (*pitr)->momentum().phi(); //phi is in range -pi to pi + double et=(*pitr)->momentum().perp(); + int p_stat=(*pitr)->status(); + int p_id = (*pitr)->pdg_id(); + cout << " eta = " << rapid<< " Phi = " << phi << " Et = " <<et/GeV << " Status= " << p_stat << " PDG ID= "<< p_id << endl; + } + } + } + } + return StatusCode::SUCCESS; +} diff --git a/Generators/TruthIO/src/HepMCReadFromFile.cxx b/Generators/TruthIO/src/HepMCReadFromFile.cxx new file mode 100644 index 0000000000000000000000000000000000000000..406bf87bdf0ade6bfcb89ca23bd507095887f93b --- /dev/null +++ b/Generators/TruthIO/src/HepMCReadFromFile.cxx @@ -0,0 +1,36 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "TruthIO/HepMCReadFromFile.h" + + +HepMCReadFromFile::HepMCReadFromFile(const std::string& name, ISvcLocator* pSvcLocator) : + GenBase(name, pSvcLocator) +{ + declareProperty("InputFile", m_input_file="events.hepmc"); + m_event_number = 0; +} + + +StatusCode HepMCReadFromFile::initialize() { + CHECK(GenBase::initialize()); + // Initialize input file + m_hepmcio.reset( new HepMC::IO_GenEvent(m_input_file.c_str(), std::ios::in) ); + // Initialize event number + m_event_number = 0; + return StatusCode::SUCCESS; +} + + +StatusCode HepMCReadFromFile::execute() { + /// @todo Should be a do-while until the read is successful or end of file? + HepMC::GenEvent* evt = m_hepmcio->read_next_event(); + if (evt) { + ++m_event_number; + evt->set_event_number(m_event_number); + GeVToMeV(evt); + events()->push_back(evt); + } + return StatusCode::SUCCESS; +} diff --git a/Generators/TruthIO/src/PrintHijingPars.cxx b/Generators/TruthIO/src/PrintHijingPars.cxx new file mode 100644 index 0000000000000000000000000000000000000000..f4f900a15f2dd5926911e047848b9b2251057757 --- /dev/null +++ b/Generators/TruthIO/src/PrintHijingPars.cxx @@ -0,0 +1,66 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "TruthIO/PrintHijingPars.h" + +#include "GeneratorObjects/HijingEventParams.h" + +PrintHijingPars::PrintHijingPars(const std::string& name, ISvcLocator* pSvcLocator) : + GenBase(name, pSvcLocator) +{ + // Set users' request + declareProperty("HijingParamsKey", m_key="Hijing_event_params"); + declareProperty("VerboseOutput", m_VerboseOutput=true); + declareProperty("FirstEvent", m_firstEvt=1); + declareProperty("LastEvent", m_lastEvt=100000000); +} + +StatusCode PrintHijingPars::initialize(){ + + ATH_MSG_INFO(">>> PrintHijingPars from Initialize"); + std::cout << "----- PrintHijingPars From initialize" << std::endl; + + if (m_lastEvt<m_firstEvt) { m_lastEvt = m_firstEvt; } + + // Initialization terminated + return StatusCode::SUCCESS; +} + +StatusCode PrintHijingPars::execute() +{ + + const HijingEventParams *hijing_pars; + const StatusCode sc = evtStore()->retrieve(hijing_pars, m_key); + if (!sc.isSuccess()) { + ATH_MSG_ERROR("Could not retrieve Hijing_event_params"); + return StatusCode::FAILURE; + } + + // printout hijing event params + // ---------------------------- + std::cout << "----- PrintHijingPars -----" << std::endl; + std::cout << "np: " << hijing_pars->get_np() << std::endl; + std::cout << "nt: " << hijing_pars->get_nt() << std::endl; + std::cout << "n0: " << hijing_pars->get_n0() << std::endl; + std::cout << "n01: " << hijing_pars->get_n01() << std::endl; + std::cout << "n10: " << hijing_pars->get_n10() << std::endl; + std::cout << "n11: " << hijing_pars->get_n11() << std::endl; + std::cout << "natt: " << hijing_pars->get_natt() << std::endl; + std::cout << "jatt: " << hijing_pars->get_jatt() << std::endl; + std::cout << "b: " << hijing_pars->get_b() << std::endl; + std::cout << "bphi: " << hijing_pars->get_bphi() << std::endl; + std::cout << "---------------------------" << std::endl; + + // End of execution for each event + return StatusCode::SUCCESS; +} + +StatusCode PrintHijingPars::finalize() { + + ATH_MSG_INFO(">>> PrintHijingPars from finalize"); + std::cout << "----- PrintHijingPars From finalize" << std::endl; + + // End of finalization step + return StatusCode::SUCCESS; +} diff --git a/Generators/TruthIO/src/PrintMC.cxx b/Generators/TruthIO/src/PrintMC.cxx new file mode 100644 index 0000000000000000000000000000000000000000..3237aa5bdd73a0155da0c58bbe58ca04594b8626 --- /dev/null +++ b/Generators/TruthIO/src/PrintMC.cxx @@ -0,0 +1,254 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// PrintMC - dump details of MC events +// Updated 29/10/2009 by Andy Buckley <andy.buckley@cern.ch> + +#include "TruthIO/PrintMC.h" +#include "GeneratorObjects/McEventCollection.h" + +#include "EventInfo/EventInfo.h" +#include "EventInfo/EventID.h" + +#include "HepPDT/ParticleData.hh" +#include "HepPDT/ParticleDataTable.hh" +#include "HepMC/GenEvent.h" + + +inline void drawLine(std::ostream& os) { + os << std::string(80, '_') << "\n" << std::endl; +} + + +PrintMC::PrintMC(const std::string& name, ISvcLocator* pSvcLocator) + : GenBase(name, pSvcLocator) +{ + // Declare JO interface + declareProperty("VerboseOutput", m_VerboseOutput=true, "Print out detailed information or not"); + declareProperty("PrintStyle", m_printsty="Barcode", "Type of info to print out: 'Barcode' or 'Vertex'"); + declareProperty("VertexInfo", m_vertexinfo=false, "Print out vertex information"); + declareProperty("TrustHepMC", m_trustHepMC=false, "Use the event numbers in the HepMC record to choose the events to print"); + declareProperty("FirstEvent", m_firstEvt=1, "Lowest event number to print"); + declareProperty("LastEvent", m_lastEvt=100000000, "Highest event number to print"); +} + + +StatusCode PrintMC::initialize() { + CHECK(GenBase::initialize()); + // Check settings + if (m_printsty != "Barcode" && m_printsty != "Vertex") { + ATH_MSG_ERROR("Unknown PrintStyle: " << m_printsty << "! Choose Barcode/Vertex"); + return StatusCode::FAILURE; + } + + // Check that last > first + if (m_lastEvt < m_firstEvt) m_lastEvt = m_firstEvt; + + return StatusCode::SUCCESS; +} + + + +/// @todo Avoid use of unprotected std::cout: stringstream + MsgStream would be better +StatusCode PrintMC::execute() { + // If output already turned off by passing last dumped event, just return + /// @todo I get the feeling VerboseOutput is being abused here... + if (!m_VerboseOutput) return StatusCode::SUCCESS; + + // Loop over all events in McEventCollection + + for (McEventCollection::const_iterator itr = events()->begin(); itr != events()->end(); ++itr) { + + // Get event number from HepMC + HepMC::GenEvent* evt = *itr; + int evtnum = evt->event_number(); + // Override with evtnum from Athena if enabled and functional + if (!m_trustHepMC) { + const EventInfo* evtinfo; + const StatusCode sc = evtStore()->retrieve(evtinfo); + if (sc.isSuccess()) { + evtnum = evtinfo->event_ID()->event_number(); + } + } + + // Return immediately if event number is outside the print range + if (evtnum < m_firstEvt) { + return StatusCode::SUCCESS; + } + if (evtnum > m_lastEvt) { + /// @todo Looks like misuse of the verbose output flag + m_VerboseOutput = false; + return StatusCode::SUCCESS; + } + + // Printout header + ATH_MSG_INFO(">>> PrintMC from execute"); + //std::cout << "Next event in the bag" << std::endl; + // int g_id = evt->signal_process_id(); + // GeneratorName_print(g_id); + // std::cout << std::endl; + + // VERTEX format + /// @todo Isn't this if (m_VerboseOutput... redundant? + if (m_VerboseOutput && m_printsty == "Vertex") { + evt->print(); // standard HepMC dump of vertex list + } + // BARCODE format + /// @todo Isn't this if (m_VerboseOutput... redundant? + else if (m_VerboseOutput && m_printsty == "Barcode") { + drawLine(std::cout); + std::cout << "GenEvent: #" << evt->event_number() + << " ID=" << evt->signal_process_id() + << " SignalProcessGenVertex Barcode: " + << ( evt->signal_process_vertex() ? evt->signal_process_vertex()->barcode() : 0 ) << "\n"; + std::cout << " Entries this event: " << evt->vertices_size() << " vertices, " + << evt->particles_size() << " particles.\n"; + + // Random State + // std::cout << " RndmState(" << evt->random_states().size() << ")="; + //for ( std::vector<long int>::const_iterator rs + // = evt->random_states().begin(); + // rs != evt->random_states().end(); rs++ ) { std::cout << *rs << " "; } + //std::cout << "\n"; + + // Weights + std::cout << " Weights(" << evt->weights().size() << ")="; + for ( HepMC::WeightContainer::const_iterator wgt = evt->weights().begin(); + wgt != evt->weights().end(); wgt++ ) { std::cout << *wgt << " "; } + std::cout << "\n"; + std::cout << " EventScale " << evt->event_scale() + << " [energy] \t alphaQCD=" << evt->alphaQCD() + << "\t alphaQED=" << evt->alphaQED() << std::endl; + + if (evt->pdf_info()) { + std::cout << "PdfInfo: id1=" << evt->pdf_info()->id1() + << " id2=" << evt->pdf_info()->id2() + << " x1=" << evt->pdf_info()->x1() + << " x2=" << evt->pdf_info()->x2() + << " q=" << evt->pdf_info()->scalePDF() + << " xpdf1=" << evt->pdf_info()->pdf1() + << " xpdf2=" << evt->pdf_info()->pdf2() + << std::endl; + } + else { + std::cout << "PdfInfo: EMPTY"; + } + + // Print a legend to describe the particle info + char particle_legend[120]; + sprintf( particle_legend," %9s %8s %-15s %4s %8s %8s (%9s,%9s,%9s,%9s,%9s)", + "Barcode","PDG ID","Name","Stat","ProdVtx","DecayVtx","Px","Py","Pz","E ","m"); + std::cout << std::endl; + std::cout << " GenParticle Legend\n" << particle_legend << "\n"; + if (m_vertexinfo) { + sprintf( particle_legend," %60s (%9s,%9s,%9s,%9s)"," ","Vx","Vy","Vz","Vct "); + std::cout << particle_legend << std::endl; + } + drawLine(std::cout); + + // Print all particles + for (HepMC::GenEvent::particle_const_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p) { + int p_bcode = (*p)->barcode(); + int p_pdg_id = (*p)->pdg_id(); + double p_px = (*p)->momentum().px(); + double p_py = (*p)->momentum().py(); + double p_pz = (*p)->momentum().pz(); + double p_pe = (*p)->momentum().e(); + int p_stat = (*p)->status(); + int p_prodvtx = 0; + if ((*p)->production_vertex() && (*p)->production_vertex()->barcode() != 0) { + p_prodvtx = (*p)->production_vertex()->barcode(); + } + int p_endvtx = 0; + if ((*p)->end_vertex() && (*p)->end_vertex()->barcode() != 0) { + p_endvtx=(*p)->end_vertex()->barcode(); + } + double v_x = 0; + double v_y = 0; + double v_z = 0; + double v_ct = 0; + if ((*p)->production_vertex()) { + v_x = (*p)->production_vertex()->position().x(); + v_y = (*p)->production_vertex()->position().y(); + v_z = (*p)->production_vertex()->position().z(); + v_ct = (*p)->production_vertex()->position().t(); + } + + // Access the PDG table to get the particle name (and mass?) + std::string sname; + double p_mass = (*p)->generated_mass(); + const HepPDT::ParticleData* ap = particleData(abs(p_pdg_id)); + if (!ap) { + ATH_MSG_DEBUG("PID " << abs(p_pdg_id) << " is not in particle data table"); + } else { + const double p_charge = ap->charge() * (p_pdg_id < 0 ? -1 : 1); // assuming that charged leptons are in the PDT... + /// @todo Only do this for hadrons? Or only if there is no generated_mass entry? Commenting out for now... + // p_mass = ap->mass().value(); + // Build particle name string + sname = ap->name(); + if (p_charge < 0) { + const size_t plusidx = sname.rfind("+"); + if (plusidx != std::string::npos) { + sname.replace(plusidx, 1, "-"); + } + } + } + // PDG table missing or exceptions + // quarks & gluons + if (p_pdg_id == 21) sname="g"; + else if (p_pdg_id == 1) sname="d"; + else if (p_pdg_id == -1) sname="d~"; + else if (p_pdg_id == 2) sname="u"; + else if (p_pdg_id == -2) sname="u~"; + else if (p_pdg_id == 3) sname="s"; + else if (p_pdg_id == -3) sname="s~"; + else if (p_pdg_id == 4) sname="c"; + else if (p_pdg_id == -4) sname="c~"; + else if (p_pdg_id == 5) sname="b"; + else if (p_pdg_id == -5) sname="b~"; + else if (p_pdg_id == 6) sname="t"; + else if (p_pdg_id == -6) sname="t~"; + // shower-specific + else if (p_pdg_id == 91) sname="cluster"; + else if (p_pdg_id == 92) sname="string"; + else if (p_pdg_id == 9922212) sname="remn"; + else if (p_pdg_id == 2101) sname="ud"; + else if (p_pdg_id == 2203) sname="uu"; + + // Calc mass if missing + /// @todo Is there no better way to detect unspecified masses than == 0? Final state photons *should* be 0 + if (p_mass == 0 && (p_stat == 2 || (p_stat != 1 && p_pdg_id != 22))) { + p_mass = (*p)->momentum().m(); + //p_mass = sqrt(p_pe*p_pe - p_px*p_px - p_py*p_py - p_pz*p_pz); + } + + const char* p_name = sname.c_str() ; + char particle_entries[120]; + sprintf(particle_entries, " %9i %8i %-15s %4i %8i %8i (%+9.3g,%+9.3g,%+9.3g,%+9.3g,%9.3g)", + p_bcode, p_pdg_id, p_name, p_stat, p_prodvtx, p_endvtx, p_px, p_py, p_pz, p_pe, p_mass); + std::cout << particle_entries << "\n"; + if (m_vertexinfo) { + sprintf(particle_entries," %60s (%+9.3g,%+9.3g,%+9.3g,%+9.3g)"," ",v_x, v_y, v_z, v_ct); + std::cout << particle_entries << "\n"; + } + } + + + // while ( !evt->vertices_empty() ) { + // HepMC::GenVertex* vtx = ( evt->m_vertex_barcodes.begin() )->second; + // int bcode = ( evt->m_vertex_barcodes.begin() )->first; + // std::cout << "vtx= "<< vtx << " bc1=" << vtx->barcode() << " bc2=" << bcode << std::endl; + // std::cout << "evt= "<< evt << " pev=" << vtx->parent_event() << std::endl; + // evt->m_vertex_barcodes.erase( evt->m_vertex_barcodes.begin() ); + // delete vtx; + // } + // evt->clear(); + // drawLine(std::cout); + + } + + } + // End of execution for each event + return StatusCode::SUCCESS; +} diff --git a/Generators/TruthIO/src/WriteHepMC.cxx b/Generators/TruthIO/src/WriteHepMC.cxx new file mode 100644 index 0000000000000000000000000000000000000000..5faf1be996ea2b2136797ea1cbdeeba28e950821 --- /dev/null +++ b/Generators/TruthIO/src/WriteHepMC.cxx @@ -0,0 +1,28 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "TruthIO/WriteHepMC.h" + + +WriteHepMC::WriteHepMC(const std::string& name, ISvcLocator* pSvcLocator) + : GenBase(name, pSvcLocator) +{ + declareProperty("OutputFile", m_outfile="events.hepmc"); + declareProperty("Precision", m_precision=8); +} + + +StatusCode WriteHepMC::initialize() { + CHECK(GenBase::initialize()); + m_hepmcio.reset( new HepMC::IO_GenEvent(m_outfile) ); + m_hepmcio->precision(m_precision); + return StatusCode::SUCCESS; +} + + +StatusCode WriteHepMC::execute() { + // Just write out the first (i.e. signal) event in the collection + m_hepmcio->write_event(event()); + return StatusCode::SUCCESS; +} diff --git a/Generators/TruthIO/src/components/TruthIO_entries.cxx b/Generators/TruthIO/src/components/TruthIO_entries.cxx new file mode 100644 index 0000000000000000000000000000000000000000..9a344e720fd26cddb274bd533a7197eb66c482c6 --- /dev/null +++ b/Generators/TruthIO/src/components/TruthIO_entries.cxx @@ -0,0 +1,20 @@ +#include "GaudiKernel/DeclareFactoryEntries.h" +#include "TruthIO/DumpMC.h" +#include "TruthIO/HepMCReadFromFile.h" +#include "TruthIO/WriteHepMC.h" +#include "TruthIO/PrintMC.h" +#include "TruthIO/PrintHijingPars.h" + +DECLARE_ALGORITHM_FACTORY(DumpMC) +DECLARE_ALGORITHM_FACTORY(HepMCReadFromFile) +DECLARE_ALGORITHM_FACTORY(PrintMC) +DECLARE_ALGORITHM_FACTORY(WriteHepMC) +DECLARE_ALGORITHM_FACTORY(PrintHijingPars) + +DECLARE_FACTORY_ENTRIES(GeneratorModules) { + DECLARE_ALGORITHM(DumpMC) + DECLARE_ALGORITHM(HepMCReadFromFile) + DECLARE_ALGORITHM(PrintMC) + DECLARE_ALGORITHM(WriteHepMC) + DECLARE_ALGORITHM(PrintHijingPars) +} diff --git a/Generators/TruthIO/src/components/TruthIO_load.cxx b/Generators/TruthIO/src/components/TruthIO_load.cxx new file mode 100644 index 0000000000000000000000000000000000000000..b4ac6f565535e248811443a7c9d95600c51abe3e --- /dev/null +++ b/Generators/TruthIO/src/components/TruthIO_load.cxx @@ -0,0 +1,3 @@ +#include "GaudiKernel/LoadFactoryEntries.h" + +LOAD_FACTORY_ENTRIES(TruthIO)