diff --git a/DataFormats/simulation/include/SimulationDataFormat/MCEventHeader.h b/DataFormats/simulation/include/SimulationDataFormat/MCEventHeader.h index 1532973778eed..b0dcd6891bb9c 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/MCEventHeader.h +++ b/DataFormats/simulation/include/SimulationDataFormat/MCEventHeader.h @@ -27,6 +27,67 @@ namespace dataformats class GeneratorHeader; +/** Common keys for information in MC event header */ +struct MCInfoKeys { + /** @{ +@name HepMC3 heavy-ion fields */ + static constexpr const char* impactParameter = "Bimpact"; + static constexpr const char* nPart = "Npart"; + static constexpr const char* nPartProjectile = "Npart_proj"; + static constexpr const char* nPartTarget = "Npart_targ"; + static constexpr const char* nColl = "Ncoll"; + static constexpr const char* nCollHard = "Ncoll_hard"; + static constexpr const char* nCollNNWounded = "NColl_NNw"; + static constexpr const char* nCollNWoundedN = "NColl_NwN"; + static constexpr const char* nCollNWoundedNwounded = "NColl_NwNW"; + static constexpr const char* planeAngle = "eventPsi"; + static constexpr const char* sigmaInelNN = "sigmaInelNN"; + static constexpr const char* centrality = "centrality"; + static constexpr const char* nSpecProjectileProton = "Nspec_proj_p"; + static constexpr const char* nSpecProjectileNeutron = "Nspec_proj_n"; + static constexpr const char* nSpecTargetProton = "Nspec_targ_p"; + static constexpr const char* nSpecTargetNeutron = "Nspec_targ_n"; + /** @} */ + /** @{ +@name HepMC3 PDF information + +In principle a header can have many of these. In that case, +each set should be prefixed with "_" where "" is a +serial number. + */ + static constexpr const char* pdfParton1Id = "pdf_parton_1_id"; + static constexpr const char* pdfParton2Id = "pdf_parton_2_id"; + static constexpr const char* pdfX1 = "pdf_x1"; + static constexpr const char* pdfX2 = "pdf_x2"; + static constexpr const char* pdfScale = "pdf_scale"; + static constexpr const char* pdfXF1 = "pdf_par_x1"; + static constexpr const char* pdfXF2 = "pdf_par_x2"; + static constexpr const char* pdfCode1 = "pdf_lhc_1_id"; + static constexpr const char* pdfCode2 = "pdf_lhc_2_id"; + /** @} */ + /** @{ +@name HepMC3 cross-section information + +In principle we can have one cross section per weight. In that +case, each should be post-fixed by "_" where "" is a +serial number. These should then matcht possible names of +weights. + */ + static constexpr const char* acceptedEvents = "accepted_events"; + static constexpr const char* attemptedEvents = "attempted_events"; + static constexpr const char* xSection = "cross_section"; + static constexpr const char* xSectionError = "cross_section_error"; + /** @} */ + /** @{ +@name Common fields */ + static constexpr const char* generator = "generator"; + static constexpr const char* generatorVersion = "version"; + static constexpr const char* processName = "processName"; + static constexpr const char* processCode = "processCode"; + static constexpr const char* weight = "weight"; + /** @} */ +}; + /*****************************************************************/ /*****************************************************************/ diff --git a/Generators/CMakeLists.txt b/Generators/CMakeLists.txt index 0b0f92955a651..6e4ebc1c2da37 100644 --- a/Generators/CMakeLists.txt +++ b/Generators/CMakeLists.txt @@ -38,6 +38,8 @@ o2_add_library(Generators src/GenCosmicsParam.cxx src/GeneratorFactory.cxx src/GeneratorGeantinos.cxx + src/GeneratorTParticle.cxx + src/GeneratorTParticleParam.cxx $<$:src/GeneratorPythia6.cxx> $<$:src/GeneratorPythia6Param.cxx> $<$:src/GeneratorPythia8.cxx> @@ -79,6 +81,8 @@ set(headers include/Generators/QEDGenParam.h include/Generators/GenCosmicsParam.h include/Generators/GeneratorGeantinos.h + include/Generators/GeneratorTParticle.h + include/Generators/GeneratorTParticleParam.h ) if (pythia6_FOUND) diff --git a/Generators/include/Generators/GeneratorHepMC.h b/Generators/include/Generators/GeneratorHepMC.h index 341f4b8a0e38a..8bcc567dd768b 100644 --- a/Generators/include/Generators/GeneratorHepMC.h +++ b/Generators/include/Generators/GeneratorHepMC.h @@ -62,7 +62,9 @@ class GeneratorHepMC : public Generator /** setters **/ void setVersion(Int_t val) { mVersion = val; }; void setFileName(std::string val) { mFileName = val; }; + void setProgCmd(std::string val) { mProgCmd = val; }; void setEventsToSkip(uint64_t val) { mEventsToSkip = val; }; + void setNEvents(unsigned int val) { mNEvents = val; } protected: /** copy constructor **/ @@ -77,11 +79,16 @@ class GeneratorHepMC : public Generator const HepMC3::FourVector getBoostedVector(const HepMC3::FourVector& vector, Double_t boost); #endif + /** methods that can be overridded **/ + void updateHeader(o2::dataformats::MCEventHeader* eventHeader) override; + /** HepMC interface **/ std::ifstream mStream; //! std::string mFileName; + std::string mProgCmd; Int_t mVersion; uint64_t mEventsToSkip; + unsigned int mNEvents; #ifdef GENERATORS_WITH_HEPMC3_DEPRECATED HepMC::Reader* mReader; //! HepMC::GenEvent* mEvent; //! diff --git a/Generators/include/Generators/GeneratorHepMCParam.h b/Generators/include/Generators/GeneratorHepMCParam.h index 4c53918c70634..b44841df3cb61 100644 --- a/Generators/include/Generators/GeneratorHepMCParam.h +++ b/Generators/include/Generators/GeneratorHepMCParam.h @@ -31,6 +31,7 @@ namespace eventgen struct GeneratorHepMCParam : public o2::conf::ConfigurableParamHelper { std::string fileName = ""; + std::string progCmd = ""; // Program command line to spawn, must write HepMC on stdout int version = 2; uint64_t eventsToSkip = 0; O2ParamDef(GeneratorHepMCParam, "HepMC"); diff --git a/Generators/include/Generators/GeneratorTParticle.h b/Generators/include/Generators/GeneratorTParticle.h new file mode 100644 index 0000000000000..8887ab344434f --- /dev/null +++ b/Generators/include/Generators/GeneratorTParticle.h @@ -0,0 +1,114 @@ +// Copyright 2023-2099 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// @author: Christian Holm Christensen +#ifndef ALICEO2_GENERATORTPARTICLE_H_ +#define ALICEO2_GENERATORTPARTICLE_H_ +#include +#include +#include + +// Forward decls +class TChain; +class TParticle; +class TClonesArray; + +namespace o2 +{ +namespace eventgen +{ +/// A class that reads in particles of class @c TParticle from a +/// branch in a @c TChain. +/// +/// Optionally, a program that generates such a @c TTree can be +/// spawn, and the @c TParticles written to a file from which this +/// object reads them in. This is done with +/// +/// --configKeyValues "TParticle.progCmd=" +/// +/// which will execute the specified EG program with the given +/// options. The EG program _must_ support the command line options +/// +/// -n NUMBER Number of events to generate +/// -o FILENAME Name of file to write to +/// +/// The tree name and particle branch names can be configured. +/// +/// --configKeyValues "TParticle.treeName=T,TParticle.branchName=P" +/// +/// File(s) to read are also configurable +/// +/// --configKeyValues "TParticle.fileNames=foo.root,bar.root" +/// +class GeneratorTParticle : public Generator +{ + public: + /** CTOR */ + GeneratorTParticle() = default; + /** CTOR */ + GeneratorTParticle(const std::string& name) + : Generator(name.c_str(), "ALICEo2 TParticle Generator") + { + } + /** DTOR */ + virtual ~GeneratorTParticle(); + + /** Initialize this generator. This will set up the chain. +Optionally, if a command line was specified by @c +TParticle.progCmd then that command line is executed in the +background and events are read from the output file of that +program */ + Bool_t Init() override; + + /** Read in the next entry from the chain. Returns false in +case of errors or no more entries to read. */ + Bool_t generateEvent() override; + + /** Import the read-in particles into the steer particle +stack */ + Bool_t importParticles() override; + + /** Set the names of files to read, separated by commas */ + void setFileNames(const std::string& val); + /** Set the name of the tree in the files. The tree _must_ +reside in the top-level directory of the files. */ + void setTreeName(const std::string& val) { mTreeName = val; } + /** Set the branch name of the branch that holds a @c +TClonesArray of @c TParticle objects */ + void setBranchName(const std::string& val) { mBranchName = val; } + /** Set child program command line to (optionally) execute */ + void setProgCmd(const std::string& val) { mProgCmd = val; } + /** Set the number of events to generate. */ + void setNEvents(unsigned int nev) { mNEvents = nev; } + + protected: + std::string mTreeName = "T"; + std::string mBranchName = "Particles"; + std::string mProgCmd = ""; + std::list mFileNames; + unsigned int mNEvents = 0; + unsigned int mEntry = 0; + TChain* mChain; + TClonesArray* mTParticles; + + void waitForData(); + + ClassDefOverride(GeneratorTParticle, 1); +}; +} // namespace eventgen +} // namespace o2 +#endif +// Local Variables: +// mode: C++ +// End: +// +// EOF +// diff --git a/Generators/include/Generators/GeneratorTParticleParam.h b/Generators/include/Generators/GeneratorTParticleParam.h new file mode 100644 index 0000000000000..2acfb2c13c623 --- /dev/null +++ b/Generators/include/Generators/GeneratorTParticleParam.h @@ -0,0 +1,41 @@ +// Copyright 2023-2099 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +// @author: Christian Holm Christensen +#ifndef ALICEO2_EVENTGEN_GENERATORTPARTICLEPARAM_H_ +#define ALICEO2_EVENTGEN_GENERATORTPARTICLEPARAM_H_ + +#include "CommonUtils/ConfigurableParam.h" +#include "CommonUtils/ConfigurableParamHelper.h" +#include + +namespace o2 +{ +namespace eventgen +{ + +/** + a parameter class/struct to keep the settings of the TGenerator + event generator and allow the user to modify them */ +struct GeneratorTParticleParam : public o2::conf::ConfigurableParamHelper { + std::string treeName = "T"; + std::string branchName = "Particles"; + std::string fileNames = "tparticle.root"; + std::string progCmd = ""; + O2ParamDef(GeneratorTParticleParam, "GeneratorTParticle"); +}; +} // end namespace eventgen +} // end namespace o2 + +#endif // ALICEO2_EVENTGEN_GENERATORHEPMCPARAM_H_ +// Local Variables: +// mode: C++ +// End: diff --git a/Generators/src/GeneratorFactory.cxx b/Generators/src/GeneratorFactory.cxx index ae47ede14f4a3..e572b5594e71e 100644 --- a/Generators/src/GeneratorFactory.cxx +++ b/Generators/src/GeneratorFactory.cxx @@ -18,6 +18,8 @@ #include #include #include +#include +#include #ifdef GENERATORS_WITH_PYTHIA6 #include #include @@ -86,6 +88,7 @@ void GeneratorFactory::setPrimaryGenerator(o2::conf::SimConfig const& conf, Fair o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance()); auto genconfig = conf.getGenerator(); + LOG(info) << "** Generator to use: '" << genconfig << "'"; if (genconfig.compare("boxgen") == 0) { // a simple "box" generator configurable via BoxGunparam auto& boxparam = BoxGunParam::Instance(); @@ -156,16 +159,32 @@ void GeneratorFactory::setPrimaryGenerator(o2::conf::SimConfig const& conf, Fair } } LOG(info) << "using external O2 kinematics"; + } else if (genconfig.compare("tparticle") == 0) { + // External ROOT file(s) with tree of TParticle in clones array, + // or external program generating such a file + auto& param = GeneratorTParticleParam::Instance(); + LOG(info) << "Init 'GeneratorTParticle' with the following parameters"; + LOG(info) << param; + auto tgen = new o2::eventgen::GeneratorTParticle(); + tgen->setFileNames(param.fileNames); + tgen->setProgCmd(param.progCmd); + tgen->setTreeName(param.treeName); + tgen->setBranchName(param.branchName); + tgen->setNEvents(conf.getNEvents()); + primGen->AddGenerator(tgen); #ifdef GENERATORS_WITH_HEPMC3 } else if (genconfig.compare("hepmc") == 0) { - // external HepMC file + // external HepMC file, or external program writing HepMC event + // records to standard output. auto& param = GeneratorHepMCParam::Instance(); LOG(info) << "Init \'GeneratorHepMC\' with following parameters"; LOG(info) << param; auto hepmcGen = new o2::eventgen::GeneratorHepMC(); hepmcGen->setFileName(param.fileName); + hepmcGen->setProgCmd(param.progCmd); hepmcGen->setVersion(param.version); hepmcGen->setEventsToSkip(param.eventsToSkip); + hepmcGen->setNEvents(conf.getNEvents()); primGen->AddGenerator(hepmcGen); #endif #ifdef GENERATORS_WITH_PYTHIA6 diff --git a/Generators/src/GeneratorHepMC.cxx b/Generators/src/GeneratorHepMC.cxx index 6c062268bfacc..a5fb510ba479f 100644 --- a/Generators/src/GeneratorHepMC.cxx +++ b/Generators/src/GeneratorHepMC.cxx @@ -14,15 +14,22 @@ #include "SimulationDataFormat/MCUtils.h" #include "Generators/GeneratorHepMC.h" #include "Generators/GeneratorHepMCParam.h" +#include "SimulationDataFormat/MCEventHeader.h" #include "HepMC3/ReaderAscii.h" #include "HepMC3/ReaderAsciiHepMC2.h" #include "HepMC3/GenEvent.h" #include "HepMC3/GenParticle.h" #include "HepMC3/GenVertex.h" #include "HepMC3/FourVector.h" +#include "HepMC3/Version.h" #include "TParticle.h" #include "TSystem.h" +#include +#include // POSIX only +#include // POISX only +#include + #include #include "FairPrimaryGenerator.h" #include @@ -85,6 +92,7 @@ Bool_t GeneratorHepMC::generateEvent() mEvent->clear(); mReader->read_event(*mEvent); if (mReader->failed()) { + LOG(error) << "Failed to read one event from input"; return kFALSE; } /** set units to desired output **/ @@ -143,6 +151,132 @@ Bool_t GeneratorHepMC::importParticles() return kTRUE; } +namespace +{ +void putAttributeInfo(o2::dataformats::MCEventHeader* eventHeader, + const std::string& name, + const std::shared_ptr& a) +{ + if (auto* p = dynamic_cast(a.get())) + eventHeader->putInfo(name, p->value()); + if (auto* p = dynamic_cast(a.get())) + eventHeader->putInfo(name, p->value()); + if (auto* p = dynamic_cast(a.get())) + eventHeader->putInfo(name, p->value()); + if (auto* p = dynamic_cast(a.get())) + eventHeader->putInfo(name, p->value()); + if (auto* p = dynamic_cast(a.get())) + eventHeader->putInfo(name, p->value()); + if (auto* p = dynamic_cast(a.get())) + eventHeader->putInfo(name, p->value()); + if (auto* p = dynamic_cast(a.get())) + eventHeader->putInfo(name, p->value()); + if (auto* p = dynamic_cast(a.get())) + eventHeader->putInfo(name, p->value()); + if (auto* p = dynamic_cast(a.get())) + eventHeader->putInfo(name, p->value()); + if (auto* p = dynamic_cast(a.get())) + eventHeader->putInfo(name, p->value()); + if (auto* p = dynamic_cast(a.get())) + eventHeader->putInfo(name, p->value()); + if (auto* p = dynamic_cast(a.get())) + eventHeader->putInfo(name, p->value()); +} +} // namespace + +/*****************************************************************/ + +void GeneratorHepMC::updateHeader(o2::dataformats::MCEventHeader* eventHeader) +{ + /** update header **/ + using Key = o2::dataformats::MCInfoKeys; + + eventHeader->putInfo(Key::generator, "hepmc"); + eventHeader->putInfo(Key::generatorVersion, HEPMC3_VERSION_CODE); + + auto xSection = mEvent->cross_section(); + auto pdfInfo = mEvent->pdf_info(); + auto hiInfo = mEvent->heavy_ion(); + + // Set default cross-section + if (xSection) { + eventHeader->putInfo(Key::xSection, xSection->xsec()); + eventHeader->putInfo(Key::xSectionError, xSection->xsec_err()); + eventHeader->putInfo(Key::acceptedEvents, + xSection->get_accepted_events()); + eventHeader->putInfo(Key::attemptedEvents, + xSection->get_attempted_events()); + } + + // Set weights and cross sections + size_t iw = 0; + for (auto w : mEvent->weights()) { + std::string post = (iw > 0 ? "_" + std::to_string(iw) : ""); + eventHeader->putInfo(Key::weight + post, w); + if (xSection) { + eventHeader->putInfo(Key::xSection, xSection->xsec(iw)); + eventHeader->putInfo(Key::xSectionError, xSection->xsec_err(iw)); + } + iw++; + } + + // Set the PDF information + if (pdfInfo) { + eventHeader->putInfo(Key::pdfParton1Id, pdfInfo->parton_id[0]); + eventHeader->putInfo(Key::pdfParton2Id, pdfInfo->parton_id[1]); + eventHeader->putInfo(Key::pdfX1, pdfInfo->x[0]); + eventHeader->putInfo(Key::pdfX2, pdfInfo->x[1]); + eventHeader->putInfo(Key::pdfScale, pdfInfo->scale); + eventHeader->putInfo(Key::pdfXF1, pdfInfo->xf[0]); + eventHeader->putInfo(Key::pdfXF2, pdfInfo->xf[1]); + eventHeader->putInfo(Key::pdfCode1, pdfInfo->pdf_id[0]); + eventHeader->putInfo(Key::pdfCode2, pdfInfo->pdf_id[1]); + } + + // Set heavy-ion information + if (hiInfo) { + eventHeader->putInfo(Key::impactParameter, + hiInfo->impact_parameter); + eventHeader->putInfo(Key::nPart, + hiInfo->Npart_proj + hiInfo->Npart_targ); + eventHeader->putInfo(Key::nPartProjectile, hiInfo->Npart_proj); + eventHeader->putInfo(Key::nPartTarget, hiInfo->Npart_targ); + eventHeader->putInfo(Key::nColl, hiInfo->Ncoll); + eventHeader->putInfo(Key::nCollHard, hiInfo->Ncoll_hard); + eventHeader->putInfo(Key::nCollNNWounded, + hiInfo->N_Nwounded_collisions); + eventHeader->putInfo(Key::nCollNWoundedN, + hiInfo->Nwounded_N_collisions); + eventHeader->putInfo(Key::nCollNWoundedNwounded, + hiInfo->Nwounded_Nwounded_collisions); + eventHeader->putInfo(Key::planeAngle, + hiInfo->event_plane_angle); + eventHeader->putInfo(Key::sigmaInelNN, + hiInfo->sigma_inel_NN); + eventHeader->putInfo(Key::centrality, hiInfo->centrality); + eventHeader->putInfo(Key::nSpecProjectileProton, hiInfo->Nspec_proj_p); + eventHeader->putInfo(Key::nSpecProjectileNeutron, hiInfo->Nspec_proj_n); + eventHeader->putInfo(Key::nSpecTargetProton, hiInfo->Nspec_targ_p); + eventHeader->putInfo(Key::nSpecTargetNeutron, hiInfo->Nspec_targ_n); + } + + for (auto na : mEvent->attributes()) { + std::string name = na.first; + if (name == "GenPdfInfo" || + name == "GenCrossSection" || + name == "GenHeavyIon") + continue; + + for (auto ia : na.second) { + int no = ia.first; + auto at = ia.second; + std::string post = (no == 0 ? "" : std::to_string(no)); + + putAttributeInfo(eventHeader, name + post, at); + } + } +} + /*****************************************************************/ Bool_t GeneratorHepMC::Init() @@ -152,14 +286,77 @@ Bool_t GeneratorHepMC::Init() /** init base class **/ Generator::Init(); - /** open file **/ std::string filename = gSystem->ExpandPathName(mFileName.c_str()); + + // If a EG command line is given, then we make a fifo on a temporary + // file, and directs the EG to write to that fifo. We will then set + // up the HepMC3 reader to read from that fifo. + // + // o2-sim -g hepmc --configKeyValues "HepMC.progCmd=" ... + // + // where is the command line to run an event generator. The + // event generator should output HepMC event records to standard + // output. Nothing else, but the HepMC event record may be output + // to standard output. If the EG has other output to standard + // output, then a filter can be set-up. For example + // + // crmc -n 3 -o hepmc3 -c /optsw/inst/etc/crmc.param -f /dev/stdout \ + // | sed -n 's/^\(HepMC::\|[EAUWVP] \)/\1/p' + // + // What's more, the event generator program must accept the command + // line argument `-n NEVENTS` to set the number of events to + // produce. + // + // Perhaps we should consider a way to set a seed on the EG. It + // could be another configuration parameter. Of course, if the EG + // program accepts a seed option, say `-s SEED`, then one could + // simply pass + // + // -s \$RANDOM + // + // to as part of the command line in `progCmd`. + // + // All of this can conviniently be achieved via a wrapper script + // around the actual EG program. + if (not mProgCmd.empty()) { + // Set filename to be a temporary name + // Should perhaps use + // + // TString base("xxxxxx"); + // auto fp = gSystem->TempFileName(base); + // fclose(fp); + // + filename = std::tmpnam(nullptr); + + // Make a fifo + int ret = mkfifo(filename.c_str(), 0600); + if (ret != 0) { + LOG(fatal) << "Failed to make fifo \"" << filename << "\""; + return false; + } + + // Build command line, rediret stdout to our fifo and put + // in the background. + std::string cmd = + mProgCmd + + " -n " + std::to_string(mNEvents) + + " > " + filename + " &"; + LOG(info) << "EG command line is \"" << cmd << "\""; + + ret = std::system(cmd.c_str()); + if (ret != 0) { + LOG(fatal) << "Failed to spawn \"" << cmd << "\""; + return false; + } + } + /** open file **/ mStream.open(filename); if (!mStream.is_open()) { LOG(fatal) << "Cannot open input file: " << filename << std::endl; return kFALSE; } + LOG(info) << "Set up reader to read from \"" << filename << "\"" << std::endl; /** create reader according to HepMC version **/ switch (mVersion) { case 2: diff --git a/Generators/src/GeneratorPythia8.cxx b/Generators/src/GeneratorPythia8.cxx index 13c17003ee612..574bed41773d8 100644 --- a/Generators/src/GeneratorPythia8.cxx +++ b/Generators/src/GeneratorPythia8.cxx @@ -202,12 +202,39 @@ Bool_t void GeneratorPythia8::updateHeader(o2::dataformats::MCEventHeader* eventHeader) { /** update header **/ - - eventHeader->putInfo("generator", "pythia8"); - eventHeader->putInfo("version", PYTHIA_VERSION_INTEGER); - eventHeader->putInfo("processName", mPythia.info.name()); - eventHeader->putInfo("processCode", mPythia.info.code()); - eventHeader->putInfo("weight", mPythia.info.weight()); + using Key = o2::dataformats::MCInfoKeys; + + eventHeader->putInfo(Key::generator, "pythia8"); + eventHeader->putInfo(Key::generatorVersion, PYTHIA_VERSION_INTEGER); + eventHeader->putInfo(Key::processName, mPythia.info.name()); + eventHeader->putInfo(Key::processCode, mPythia.info.code()); + eventHeader->putInfo(Key::weight, mPythia.info.weight()); + + auto& info = mPythia.info; + + // Set PDF information + eventHeader->putInfo(Key::pdfParton1Id, info.id1pdf()); + eventHeader->putInfo(Key::pdfParton2Id, info.id2pdf()); + eventHeader->putInfo(Key::pdfX1, info.x1pdf()); + eventHeader->putInfo(Key::pdfX2, info.x2pdf()); + eventHeader->putInfo(Key::pdfScale, info.QFac()); + eventHeader->putInfo(Key::pdfXF1, info.pdf1()); + eventHeader->putInfo(Key::pdfXF2, info.pdf2()); + + // Set cross section + eventHeader->putInfo(Key::xSection, info.sigmaGen() * 1e9); + eventHeader->putInfo(Key::xSectionError, info.sigmaErr() * 1e9); + + // Set weights (overrides cross-section for each weight) + size_t iw = 0; + auto xsecErr = info.weightContainerPtr->getTotalXsecErr(); + for (auto w : info.weightContainerPtr->getTotalXsec()) { + std::string post = (iw == 0 ? "" : "_" + std::to_string(iw)); + eventHeader->putInfo(Key::weight + post, info.weightValueByIndex(iw)); + eventHeader->putInfo(Key::xSection + post, w * 1e9); + eventHeader->putInfo(Key::xSectionError + post, xsecErr[iw] * 1e9); + iw++; + } #if PYTHIA_VERSION_INTEGER < 8300 auto hiinfo = mPythia.info.hiinfo; @@ -218,7 +245,7 @@ void GeneratorPythia8::updateHeader(o2::dataformats::MCEventHeader* eventHeader) if (hiinfo) { /** set impact parameter **/ eventHeader->SetB(hiinfo->b()); - eventHeader->putInfo("Bimpact", hiinfo->b()); + eventHeader->putInfo(Key::impactParameter, hiinfo->b()); auto bImp = hiinfo->b(); /** set Ncoll, Npart and Nremn **/ int nColl, nPart; @@ -230,7 +257,8 @@ void GeneratorPythia8::updateHeader(o2::dataformats::MCEventHeader* eventHeader) getNpart(nPartProtonProj, nPartNeutronProj, nPartProtonTarg, nPartNeutronTarg); getNremn(nRemnProtonProj, nRemnNeutronProj, nRemnProtonTarg, nRemnNeutronTarg); getNfreeSpec(nFreeNeutronProj, nFreeProtonProj, nFreeNeutronTarg, nFreeProtonTarg); - eventHeader->putInfo("Ncoll", nColl); + eventHeader->putInfo(Key::nColl, nColl); + // These are all non-HepMC3 fields - of limited use eventHeader->putInfo("Npart", nPart); eventHeader->putInfo("Npart_proj_p", nPartProtonProj); eventHeader->putInfo("Npart_proj_n", nPartNeutronProj); @@ -244,6 +272,18 @@ void GeneratorPythia8::updateHeader(o2::dataformats::MCEventHeader* eventHeader) eventHeader->putInfo("Nfree_proj_p", nFreeProtonProj); eventHeader->putInfo("Nfree_targ_n", nFreeNeutronTarg); eventHeader->putInfo("Nfree_targ_p", nFreeProtonTarg); + + // --- HepMC3 conforming information --- + // This is how the Pythia authors define Ncoll + // eventHeader->putInfo(Key::nColl, + // hiinfo->nAbsProj() + hiinfo->nDiffProj() + + // hiinfo->nAbsTarg() + hiinfo->nDiffTarg() - + // hiiinfo->nCollND() - hiinfo->nCollDD()); + eventHeader->putInfo(Key::nPartProjectile, + hiinfo->nAbsProj() + hiinfo->nDiffProj()); + eventHeader->putInfo(Key::nPartTarget, + hiinfo->nAbsTarg() + hiinfo->nDiffTarg()); + eventHeader->putInfo(Key::nCollHard, hiinfo->nCollNDTot()); } } @@ -302,6 +342,10 @@ void GeneratorPythia8::getNcoll(const Pythia8::Info& info, int& nColl) auto hiinfo = info.hiInfo; #endif + // This is how the Pythia authors define Ncoll + nColl = (hiinfo->nAbsProj() + hiinfo->nDiffProj() + + hiinfo->nAbsTarg() + hiinfo->nDiffTarg() - + hiinfo->nCollND() - hiinfo->nCollDD()); nColl = 0; if (!hiinfo) { @@ -330,6 +374,16 @@ void GeneratorPythia8::getNpart(const Pythia8::Info& info, int& nPart) /** compute number of participants as the sum of all participants nucleons **/ + // This is how the Pythia authors calculate Npart +#if PYTHIA_VERSION_INTEGER < 8300 + auto hiinfo = info.hiinfo; +#else + auto hiinfo = info.hiInfo; +#endif + if (hiinfo) + nPart = (hiinfo->nAbsProj() + hiinfo->nDiffProj() + + hiinfo->nAbsTarg() + hiinfo->nDiffTarg()); + int nProtonProj, nNeutronProj, nProtonTarg, nNeutronTarg; getNpart(info, nProtonProj, nNeutronProj, nProtonTarg, nNeutronTarg); nPart = nProtonProj + nNeutronProj + nProtonTarg + nNeutronTarg; diff --git a/Generators/src/GeneratorTParticle.cxx b/Generators/src/GeneratorTParticle.cxx new file mode 100644 index 0000000000000..d04f0116bfe98 --- /dev/null +++ b/Generators/src/GeneratorTParticle.cxx @@ -0,0 +1,186 @@ +// Copyright 2023-2099 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// @author Christian Holm Christensen +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // POSIX only +#include // POISX only +#include + +namespace o2 +{ +namespace eventgen +{ +GeneratorTParticle::~GeneratorTParticle() +{ + if (mChain) { + TFile* file = mChain->GetCurrentFile(); + if (file) + mChain->RecursiveRemove(file); + delete mChain; + } + if (mProgCmd.empty()) + return; + + // Get the file we're reading from + std::filesystem::path p(mFileNames.front()); + + // Wait until child process creates the file + if (not std::filesystem::exists(p)) + return; + + // Remove temporary file + std::error_code ec; + std::filesystem::remove(p, ec); +} + +void GeneratorTParticle::setFileNames(const std::string& val) +{ + std::stringstream s; + std::string f; + while (std::getline(s, f, ',')) + mFileNames.push_back(f); +} + +Bool_t GeneratorTParticle::Init() +{ + mChain = new TChain(mTreeName.c_str()); + mTParticles = new TClonesArray("TParticle"); + mChain->SetBranchAddress(mBranchName.c_str(), &mTParticles); + + if (not mProgCmd.empty()) { + // Set filename to be a temporary name + // Should perhaps use + // + // TString base("xxxxxx"); + // auto fp = gSystem->TempFileName(base); + // fclose(fp); + // + std::string filename = std::tmpnam(nullptr); + + // Build command line, Assumes command line parameter + // + // -n NUMBER of events to produce + // -o FILENAME of output file + // + // A script can be wrapped around existing EGs to ensure these + // options are observed. + std::string cmd = + mProgCmd + + " -n " + std::to_string(mNEvents) + + " -o " + filename + " &"; + LOG(info) << "EG command line is \"" << cmd << "\""; + + int ret = std::system(cmd.c_str()); + if (ret != 0) { + LOG(fatal) << "Failed to spawn \"" << cmd << "\""; + return false; + } + + mFileNames.clear(); + mFileNames.push_back(filename); + } + for (auto filename : mFileNames) + mChain->AddFile(filename.c_str()); + + return true; +} + +void GeneratorTParticle::waitForData() +{ + if (mProgCmd.empty()) + return; // Not from child process + + using namespace std::chrono_literals; + + // Get the file we're reading from + std::filesystem::path p(mFileNames.front()); + + LOG(info) << "Waiting for data on " << p; + + // Wait until child process creates the file + while (not std::filesystem::exists(p)) + std::this_thread::sleep_for(500ms); + + // Wait until we have more data in the file than just the file + // header + while (std::filesystem::file_size(p) <= 256) + std::this_thread::sleep_for(500ms); + + // Give the child process 1 second to post the data to the file + LOG(info) << "Got data in " << p << ", sleeping for a while"; + std::this_thread::sleep_for(1s); +} + +Bool_t GeneratorTParticle::generateEvent() +{ + if (mEntry == 0) + waitForData(); + + int read = mChain->GetEntry(mEntry); + mEntry++; + + if (read < 0) + LOG(error) << "Failed to read entry " << mEntry << " of chain"; + + if (read <= 0) + return false; + + return true; +} + +Bool_t GeneratorTParticle::importParticles() +{ + for (auto* object : *mTParticles) { + TParticle* particle = static_cast(object); + auto statusCode = particle->GetStatusCode(); + if (!mcgenstatus::isEncoded(statusCode)) + statusCode = mcgenstatus::MCGenStatusEncoding(statusCode, 0) + .fullEncoding; + + mParticles.emplace_back(particle->GetPdgCode(), + statusCode, + particle->GetFirstMother(), + particle->GetSecondMother(), + particle->GetFirstDaughter(), + particle->GetLastDaughter(), + particle->Px(), + particle->Py(), + particle->Pz(), + particle->Energy(), + particle->Vx(), + particle->Vy(), + particle->Vz(), + particle->T()); + auto& tgt = mParticles[mParticles.size() - 1]; + tgt.SetPolarTheta(particle->GetPolarTheta()); + tgt.SetPolarPhi(particle->GetPolarPhi()); + tgt.SetCalcMass(particle->GetCalcMass()); + tgt.SetWeight(particle->GetWeight()); + } + return true; +} +} // namespace eventgen +} // namespace o2 +// +// EOF +// diff --git a/Generators/src/GeneratorTParticleParam.cxx b/Generators/src/GeneratorTParticleParam.cxx new file mode 100644 index 0000000000000..abb2d39681a3e --- /dev/null +++ b/Generators/src/GeneratorTParticleParam.cxx @@ -0,0 +1,18 @@ +// Copyright 2023-2099 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// @author Christian Holm Christensen + +#include "Generators/GeneratorTParticleParam.h" +O2ParamImpl(o2::eventgen::GeneratorTParticleParam); +// +// EOF +// diff --git a/Generators/src/GeneratorsLinkDef.h b/Generators/src/GeneratorsLinkDef.h index 0d42b2ce505b6..669b1f37019c1 100644 --- a/Generators/src/GeneratorsLinkDef.h +++ b/Generators/src/GeneratorsLinkDef.h @@ -67,5 +67,8 @@ #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::eventgen::QEDGenParam> + ; #pragma link C++ class o2::eventgen::GenCosmicsParam + ; #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::eventgen::GenCosmicsParam> + ; +#pragma link C++ class o2::eventgen::GeneratorTParticle + ; +#pragma link C++ class o2::eventgen::GeneratorTParticleParam + ; +#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::eventgen::GeneratorTParticleParam> + ; #endif diff --git a/run/SimExamples/HepMC/README.md b/run/SimExamples/HepMC/README.md new file mode 100644 index 0000000000000..fec93402e9ea2 --- /dev/null +++ b/run/SimExamples/HepMC/README.md @@ -0,0 +1,148 @@ + + +Here are pointers on how to use `GeneratorHepMC` selected by the +option `-g hepmc` for `o2-sim`. + +## Reading HepMC files + +The generator `GeneratorHepMC` can read events from a +[HepMC(3)](http://hepmc.web.cern.ch/hepmc/) formatted file. These +files can be produced by a standalone event generator program (EG). +Examples of such programs are + +- [Pythia8](https://pythia.org) +- The [CRMC](https://gitlab.iap.kit.edu/AirShowerPhysics/crmc) suite +- [Herwig](https://herwig.hepforge.org/) +- [SMASH](https://smash-transport.github.io/) +- ... and many others + +Please refer to the documentation of these for more on how to make +event files in the HepMC format. + +To make a simulation reading from the file `events.hepmc`, do + + o2-sim -g hepmc --configKeyValues "HepMC.fileName=events.hepmc" ... + +See also [`read.sh`](read.sh). + +## Reading HepMC events from child process +`GeneratorHepMC` can not only read HepMC events from a file, but can +also spawn an child EG to produce events. Suppose we have a program +named `eg` which is some EG that writes HepMC event records to the +standard output. Then we can execute a simulation using this external +EG by + + o2-sim -g hepmc --configKeyValues "HepMC.progCmd=eg" + +See also [`child.sh`](child.sh). + +There are some requirements on the program `eg`: + +- It _must_ write the HepMC event structures to standard output + (`/dev/stdout`). +- It may _not_ write other information to standard output. +- It _must_ accept the option `-n n-events` to set the number of + events to produce to `n-events`. + +If a program does not adhere to these requirements, it will often be +simple enough to make a small wrapper script that enforce this. For +example, `crmc` will write a lot of information to standard output. +We can filter that out via a shell script ([`crmc.sh`](crmc.sh)) like + + #!/bin/sh + crmc $@ -o hepmc3 -f /dev/stdout | sed -n 's/^\(HepMC::\|[EAUWVP] \)/\1/p' + +The `sed` command selects lines that begin with `HepMC::`, or one +of single characters `E` (event), `A` (attribute), `U` (units), `W` +(weight), `V` (vertex), or `P` (particle) followed by a space. This +should in most cases be enough to filter out extra stuff written on +standard output. + +The script above also passes any additional command line options on to +`crmc` via `$@`. We can utilise this with `o2-sim` to set options to +the CRMC suite. For example, if we want to simulate p-Pb collisions +using DpmJET, we can do + + o2-sim -g hepmc --configKeyValues "HepMC.progCmd=crmc.sh -m 12 -i2212 -I 1002080820" + + +### Implementation details + +Internally `GeneratorHepMC` + +1. creates a unique temporary file name in the working directory, +2. then creates a FIFO (or named pipe, see + [Wikipedia](https://en.wikipedia.org/wiki/Named_pipe)), +3. builds a command line, e.g., + + eg options > fifo-name & + +4. and executes that command line + +## The future + +The `GeneratorHepMC` (and sister generator `GeneratorTParticle`) will +in the not so distant future be upgraded with new functionality to +more easily customise reading files and executing a child process. In +particular + +- HepMC event structures can be read from any file format supported by + HepMC it self (see + [here](http://hepmc.web.cern.ch/hepmc/group__IO.html) and + [here](http://hepmc.web.cern.ch/hepmc/group__factory.html). + +- New options that can be specified in `--configKeyValues` + + - `HepMC.eventsToSkip=number` a number events to skip at the + beginning of each file read. + + - `FileOrCmd.fileNames=list` a comma separated list of HepMC files + to read. + + - `FileOrCmd.cmd=command line` a command line to execute as a + background child process. If this is set (not the empty string), + then `FileOrCmd.fileNames` is ignored. + + - A number of keys that specifies the command line option switch + that the child program accepts for certain things. If any of + these are set to the empty string, then that switch and + corresponding option value is not passed to the child program. + + - `FileOrCmd.outputSwitch=switch` (default `>`) to specify output + file. The default of `>` assumes that the program write HepMC + events, and _only_ those, to standard output. + + - `FileOrCmd.seedSwitch=switch` (default `-s`) to specify the + random number generator seed. The value passed is selected by + the `o2-sim` option `--seed` + + - `FileOrCmd.bMaxSwitch=switch` (default `-b`) to specify the + upper limit on the impact parameters sampled. The value passed + is selected by the `o2-sim` option `--bMax` + + - `FileOrCmd.nEventsSwitch=switch` (default `-n`) to specify the + number of events to generate. The value passed is selected by + the `o2-sim` option `--nEvents` or (`-n`) + + - `FileOrCmd.backgroundSwitch=switch` (default `&`) to specify how + the program is put in the background. Typically this should be + `&`, but a program may itself fork to the background. + +- Some options are no longer available + - `HepMC.fileName` - use `FileOrCmd.fileNames` + - `HepMC.progCmd` - use `FileOrCmd.cmd` + +The command line build will now be + +> _commandLine_ _nEventsSwitch_ _nEvents_ _seedSwitch_ _seed_ +> _bMaxSwitch_ _bMax_ _outputSwitch_ _output_ _backgroundSwitch_ + +If any of the `Switch` keys are empty, then the corresponding option +is not propagated to the command line. For example, if _bMaxSwitch_ +is empty, then the build command line will be + +> _commandLine_ _nEventsSwitch_ _nEvents_ _seedSwitch_ _seed_ +> _outputSwitch_ _output_ _backgroundSwitch_ + diff --git a/run/SimExamples/HepMC/child.sh b/run/SimExamples/HepMC/child.sh new file mode 100755 index 0000000000000..3a4f9ba3e40c7 --- /dev/null +++ b/run/SimExamples/HepMC/child.sh @@ -0,0 +1,51 @@ +#!/usr/bin/bash + +cmd="./crmc.sh" +seed=$RANDOM +nev=1 +out= + +usage() +{ + cat </dev/stderr + exit 1 + ;; + esac + shift +done + +if test "x$out" = "x" ; then + out=`echo $cmd | sed 's,^\./,,' | tr '[$/. ]' '_'` +fi +out=`echo "$out" | tr ' ' '_'` + +export VMCWORKDIR=${O2_ROOT}/share +o2-sim -g hepmc --configKeyValues "HepMC.progCmd=$cmd" \ + --outPrefix "$out" --seed $seed --nEvents $nev $@ diff --git a/run/SimExamples/HepMC/crmc.sh b/run/SimExamples/HepMC/crmc.sh new file mode 100755 index 0000000000000..347de0446e5fd --- /dev/null +++ b/run/SimExamples/HepMC/crmc.sh @@ -0,0 +1,7 @@ +#!/bin/sh +# This script _may not_ write to standard out, except for the HepMC +# event record. + +crmcParam=$(dirname $(dirname `which crmc`))/etc/crmc.param +exec crmc -c $crmcParam $@ -o hepmc3 -f /dev/stdout | \ + sed -n 's/^\(HepMC::\|[EAUWVP] \)/\1/p' diff --git a/run/SimExamples/HepMC/read.sh b/run/SimExamples/HepMC/read.sh new file mode 100755 index 0000000000000..96f0415fa13ac --- /dev/null +++ b/run/SimExamples/HepMC/read.sh @@ -0,0 +1,48 @@ +#!/usr/bin/bash + +inp=events.hepmc +seed=$RANDOM +nev=1 +out= + +usage() +{ + cat </dev/stderr + exit 1 + ;; + esac + shift +done + +if test "x$out" = "x" ; then + out=`basename $inp .hepmc` +fi +out=`echo "$out" | tr ' ' '_'` + +export VMCWORKDIR=${O2_ROOT}/share +o2-sim -g hepmc --configKeyValues "HepMC.fileName=$inp" \ + --outPrefix "$out" --seed $seed --nEvents $nev $@ diff --git a/run/SimExamples/README.md b/run/SimExamples/README.md index 2613f6bb2773c..ae2d15d47316c 100644 --- a/run/SimExamples/README.md +++ b/run/SimExamples/README.md @@ -11,6 +11,7 @@ * \subpage refrunSimExamplesAdaptive_Pythia8 * \subpage refrunSimExamplesAliRoot_Hijing * \subpage refrunSimExamplesAliRoot_AMPT +* \subpage refrunSimExamplesHepMC * \subpage refrunSimExamplesHepMC_STARlight * \subpage refrunSimExamplesJet_Embedding_Pythia8 * \subpage refrunSimExamplesStepMonitoringSimple1 @@ -20,4 +21,5 @@ * \subpage refrunSimExamplesSelective_Transport_pi0 * \subpage refrunSimExamplesCustom_EventInfo * \subpage refrunSimExamplesMCTrackToDPL +* \subpage refrunSimExamplesTParticle /doxy --> diff --git a/run/SimExamples/TParticle/MyEG.macro b/run/SimExamples/TParticle/MyEG.macro new file mode 100644 index 0000000000000..a067b4c1051b7 --- /dev/null +++ b/run/SimExamples/TParticle/MyEG.macro @@ -0,0 +1,128 @@ +#include +#include +#include +#include +#include +#include + +//-------------------------------------------------------------------- +// Our generator class. Really simple. +class MyGenerator : public TGenerator +{ + public: + Long_t projectilePDG; + Long_t targetPDG; + Double_t sqrts; + MyGenerator() {} + void Initialize(Long_t projectile, + Long_t target, + Double_t sqrts) + { + this->projectilePDG = projectile; + this->targetPDG = target; + this->sqrts = sqrts; + } + void GenerateEvent() + { /* Do something */ + } + TObjArray* ImportParticles(Option_t* option = "") { return 0; } + Int_t ImportParticles(TClonesArray* particles, Option_t* option = "") + { + Int_t nParticles = 10; + Int_t iParticle = 0; + // Make beam particles + new ((*particles)[iParticle++]) TParticle(projectilePDG, 4, -1, -1, + 2, nParticles - 1, + 0, 0, sqrts / 2, + TMath::Sqrt(1 + sqrts * sqrts), + 0, 0, 0, 0); + new ((*particles)[iParticle++]) TParticle(projectilePDG, 4, -1, -1, + 2, nParticles - 1, + 0, 0, -sqrts / 2, + TMath::Sqrt(1 + sqrts * sqrts), + 0, 0, 0, 0); + for (; iParticle < nParticles; iParticle++) + new ((*particles)[iParticle]) + TParticle(211, 1, 0, 1, + -1, -1, + 0.1 * iParticle, + 0.1 * iParticle, + 0.1 * iParticle, + TMath::Sqrt(0.03 * iParticle * iParticle + 0.14 * 0.14), + 0, 0, 0, 0); + + return nParticles; + } +}; +//-------------------------------------------------------------------- +// Our steering class +struct MySteer { + TGenerator* generator; + TFile* file; + TTree* tree; + TClonesArray* particles; + Int_t every; + MySteer(TGenerator* generator, const TString& output, Int_t every) + : generator(generator), + file(TFile::Open(output, "RECREATE")), + tree(new TTree("T", "T")), + particles(new TClonesArray("TParticle")), + every(every) + { + tree->SetDirectory(file); + tree->Branch("Particles", &particles); + } + ~MySteer() + { + close(); + } + void event() + { + particles->Clear(); + generator->GenerateEvent(); + generator->ImportParticles(particles); + tree->Fill(); + } + void sync() + { + // Important so that GeneratorTParticle picks up the events as + // they come. + tree->AutoSave("SaveSelf FlushBaskets Overwrite"); + } + void run(Int_t nev) + { + for (Int_t iev = 0; iev < nev; iev++) { + event(); + + if (every > 0 and (iev % every == 0) and iev != 0) + sync(); + } + } + void close() + { + if (not file) + return; + file->Write(); + file->Close(); + file = nullptr; + } +}; + +//-------------------------------------------------------------------- +// Our steering function +void MyEG(Int_t nev, const TString& out, Int_t seed, Int_t every = 1) +{ + gRandom->SetSeed(seed); + + MyGenerator* eg = new MyGenerator(); + eg->Initialize(2212, 2212, 5200); + + MySteer steer(eg, out, every); + steer.run(nev); +} +// Local Variables: +// mode: C++ +// End: +// +// EOF +// diff --git a/run/SimExamples/TParticle/README.md b/run/SimExamples/TParticle/README.md new file mode 100644 index 0000000000000..b07080d51f044 --- /dev/null +++ b/run/SimExamples/TParticle/README.md @@ -0,0 +1,258 @@ + + +Here are pointers on how to use `GeneratorTParticle` selected by the +option `-g tparticle` for `o2-sim`. + +## Reading TParticle files + +The generator `GeneratorTParticle` can read events from a ROOT file +containing a `TTree` with a branch holding a `TClonesArray` of +`TParticle` objects. These files can be produced by a standalone event +generator program (EG). + +To make a simulation reading from the file `particles.root`, do + + o2-sim -g tparticle --configKeyValues "GeneratorTParticle.fileName=particles.root" ... + +See also [`read.sh`](read.sh). Do + + ./read.sh --help + +for a list of options. This expects an input file with a `TTree` +with a single `TBranch` holding a `TClonesArray` of `TParticle` +objects. One such example file can be made with [`myeg.sh`]. Do + + ./myeg.sh --help + +for a list of options. + +### Configurations + +- The name of the `TTree` to read can be set by the configuration key + `GeneratorTParticle.treeName`. The default is `T` +- The name of the `TBranch` holding the `TClonesArray` of `TParticle` + objects can be set by the configuration key + `GeneratorTParticle.branchName`. The default is `Particles` + +For example + + o2-sim -g tparticle --configKeyValues "GeneratorTParticle.fileName=particles.root;GeneratorTParticle.treeName=Events;GeneratorTParticle.branchName=Tracks" ... + +## Reading TParticle events from child process + +`GeneratorTParticle` can not only read events from a file, but can +also spawn an child EG to produce events. Suppose we have a program +named `eg` which is some EG that writes `TParticle` event records to a +file . Then we can execute a simulation using this external +EG by + + o2-sim -g tgenerator --configKeyValues "GeneratorTParticle.progCmd=eg" + +See also [`child.sh`](child.sh). Do + + ./child.sh --help + +for a list of options. + +There are some requirements on the program `eg`: + +- It _must_ accept the option `-n n-events` to set the number of + events to produce to `n-events`. + +- It _must_ accept the option `-o output` to set the output file name. + +If a program does not adhere to these requirements, it will often be +simple enough to make a small wrapper script that enforce this. + +### Configurations + +Same as above. + +### Example EG + +The child-process feature allows us to use almost any EG for which we +have a `TGenerator` interface without compiling it in to O2. Suppose +we have defined the class `MyGenerator` to produce events. + + class MyGenerator : public TGenerator { + public: + MyGenerator(); + void Initialize(Long_t projectile, + Long_t target, + Double_t sqrts); + void GenerateEvent(); + Int_t ImportParticles(TClonesArray* particles,Option_t*option=""); + }; + +and a steering class + + struct MySteer { + TGenerator* generator; + TFile* file; + TTree* tree; + TClonesArray* particle; + Int_t flushEvery; + MySteer(TGenerator* generator, + const TString& output, + Int_t flushEvery) + : generator(generator) + file(TFile::Open(output,"RECREATE")), + tree("T","Particle tree"), + particles(new TClonesArray("TParticle")), + flushEvery(flushEvery) + { + tree->SetDirectory(file); + tree->Branch("Particles",&particles); + } + ~MySteer() { close(); } + void event() { + particles->Clear(); + generator->GenerateEvent(); + generator->ImportParticles(particles); + tree->Fill(); + } + void sync() { + tree->AutoSave("SaveSelf FlushBaskets Overwrite"); + } + void run(Int_t nev) { + for (Int_t iev = 0; iev < nev; iev++) { + event(); + if (flushEvery > 0 and (iev % flushEvery == 0) and iev != 0) + sync(); + } + } + void close() { + if (not file) return; + file->Write(); + file->Close(); + file = nullptr; + } + }; + +Then we could make the script [`MyEG.macro` (complete code)](MyEG.macro) like + + void MyEG(Int_t nev,const TString& out,Int_t every=1) + { + MyGenerator* eg = new MyGenerator(); + eg->Initialize(2212, 2212, 5200); + + MySteer steer(eg, out, every); + steer.run(nev); + } + +and a simple shell-script [`myeg.sh`](myeg.sh) to pass arguments to +the `MyEG.macro` script + + #!/bin/sh + + nev=1 + out=particles.root + + while test $# -gt 0 ; do + case $1 in + -n) nev=$2 ; shift ;; + -o) out=$2 ; shift ;; + *) ;; + esac + shift + done + + root -l MyEG.macro -- $nev \"$out\" + +We can then do + + o2-sim -g tgenerator --configKeyValues "GeneratorTParticle.progCmd=./myeg.sh" + +to produce events with our generator `MyGenerator`. + + +### Implementation details + +Internally `GeneratorTParticle` + +1. creates a unique temporary file name in the working directory, +2. builds a command line, e.g., + + eg options -o temporary-name & + +3. and executes that command line + +## The future + +The `GeneratorTParticle` (and sister generator `GeneratorHepMC`) will +in the not so distant future be upgraded with new functionality to +more easily customise reading files and executing a child process. In +particular + +- New options that can be specified in `--configKeyValues` + + - `GeneratorTParticle.treeName=name` the name of the `TTree` in the + input files. + + - `GeneratorTParticle.branchName=name` the name of the `TBranch` in + the `TTree` that holds the `TClonesArray` of `TParticle` objects. + + - `FileOrCmd.fileNames=list` a comma separated list of HepMC files + to read + + - `FileOrCmd.cmd=command line` a command line to execute as a + background child process. If this is set (not the empty string), + then `FileOrCmd.fileNames` is ignored. + + - A number of keys that specifies the command line option switch + that the child program accepts for certain things. If any of + these are set to the empty string, then that switch and + corresponding option value is not passed to the child program. + + - `FileOrCmd.outputSwitch=switch` (default `>`) to specify output + file. The default of `>` assumes that the program write HepMC + events, and _only_ those, to standard output. + + - `FileOrCmd.seedSwitch=switch` (default `-s`) to specify the + random number generator seed. The value passed is selected by + the `o2-sim` option `--seed` + + - `FileOrCmd.bMaxSwitch=switch` (default `-b`) to specify the + upper limit on the impact parameters sampled. The value passed + is selected by the `o2-sim` option `--bMax` + + - `FileOrCmd.nEventsSwitch=switch` (default `-n`) to specify the + number of events to generate. The value passed is selected by + the `o2-sim` option `--nEvents` or (`-n`) + + - `FileOrCmd.backgroundSwitch=switch` (default `&`) to specify how + the program is put in the background. Typically this should be + `&`, but a program may itself fork to the background. + +- Some options are no longer available + + - `GeneratorTParticle.fileName` - use `FileOrCmd.fileNames` + - `GeneratorTParticle.progCmd` - use `FileOrCmd.cmd` + +The command line build will now be + +> _commandLine_ _nEventsSwitch_ _nEvents_ _seedSwitch_ _seed_ +> _bMaxSwitch_ _bMax_ _outputSwitch_ _output_ _backgroundSwitch_ + +If any of the `Switch` keys are empty or set to `none`, then the +corresponding option is not propagated to the command line. For +example, if _bMaxSwitch_ is empty, then the build command line will be + +> _commandLine_ _nEventsSwitch_ _nEvents_ _seedSwitch_ _seed_ +> _outputSwitch_ _output_ _backgroundSwitch_ + + +### Header information + +The class `GeneratorTParticle` will take a key parameter, say +`headerName` which will indicate a branch that contains header +information. Under that branch, the class will then search for leaves +(`TLeaf`) that correspond to standard header information keys (see +o2::dataformats::MCInfoKeys). If any of those leaves are present, +then the corresponding keys will be set on the generated event header. + +Thus, as long as the generator observes the convention used, we can +also import auxiliary information (impact parameter, Npart, ...) from +the input files in addition to the particle information. diff --git a/run/SimExamples/TParticle/child.sh b/run/SimExamples/TParticle/child.sh new file mode 100755 index 0000000000000..802ebe4a3689a --- /dev/null +++ b/run/SimExamples/TParticle/child.sh @@ -0,0 +1,58 @@ +#!/usr/bin/bash + +cmd="./myeg.sh" +seed=$RANDOM +nev=1 +out= +opt= + +usage() +{ + cat </dev/stderr + exit 1 + ;; + esac + shift +done + +if test "x$out" = "x" ; then + out=`echo $cmd | sed 's,^\./,,' | tr '[$/. ]' '_'` +fi +out=`echo "$out" | tr ' ' '_'` + +keys="Generator.progCmd=$cmd $opt" +# keys="FileOrCmd.cmd=$cmd $opt;FileOrCmd.outputSwitch=-o" +set -x +export VMCWORKDIR=${O2_ROOT}/share +o2-sim -g tparticle --configKeyValues "$keys" \ + --outPrefix "$out" --seed $seed --nEvents $nev $@ + diff --git a/run/SimExamples/TParticle/myeg.sh b/run/SimExamples/TParticle/myeg.sh new file mode 100755 index 0000000000000..4e76c8a8a02d8 --- /dev/null +++ b/run/SimExamples/TParticle/myeg.sh @@ -0,0 +1,37 @@ +#!/bin/sh + + +nev=1 +out=particles.root +seed=0 +opt= + +usage() +{ + cat </dev/stderr + exit 1 + ;; + esac + shift +done + +if test "x$out" = "x" ; then + out=`basename $inp .root` +fi +out=`echo "$out" | tr ' ' '_'` + +set +e + +# Future FileOrCmd.fileNames=${inp} +export VMCWORKDIR=${O2_ROOT}/share +o2-sim -g tparticle --configKeyValues "GeneratorTParticle.fileName=${inp}" \ + --outPrefix "$out" --seed $seed --nEvents $nev $@ +