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..f8551a9e5e8ba 100644 --- a/Generators/CMakeLists.txt +++ b/Generators/CMakeLists.txt @@ -29,6 +29,8 @@ o2_add_library(Generators src/GeneratorExternalParam.cxx src/GeneratorFromFile.cxx src/GeneratorFromO2KineParam.cxx + src/GeneratorFileOrCmd.cxx + src/GeneratorFileOrCmdParam.cxx src/PrimaryGenerator.cxx src/PrimaryGeneratorParam.cxx src/TriggerExternalParam.cxx @@ -38,6 +40,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> @@ -71,6 +75,8 @@ set(headers include/Generators/GeneratorExternalParam.h include/Generators/GeneratorFromFile.h include/Generators/GeneratorFromO2KineParam.h + include/Generators/GeneratorFileOrCmd.h + include/Generators/GeneratorFileOrCmdParam.h include/Generators/PrimaryGenerator.h include/Generators/PrimaryGeneratorParam.h include/Generators/TriggerExternalParam.h @@ -79,6 +85,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/GeneratorFileOrCmd.h b/Generators/include/Generators/GeneratorFileOrCmd.h new file mode 100644 index 0000000000000..0c9c618928549 --- /dev/null +++ b/Generators/include/Generators/GeneratorFileOrCmd.h @@ -0,0 +1,243 @@ +// Copyright 2019-2020 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_GENERATORFILEORCMD_H_ +#define ALICEO2_EVENTGEN_GENERATORFILEORCMD_H_ +#include +#include +#include + +namespace o2 +{ +namespace conf +{ +class SimConfig; +} +namespace eventgen +{ + +/** Service class for either reading from a file or executing a + program writing to a specific file */ +struct GeneratorFileOrCmd { + /** + * Configure the generator from parameters and the general + * simulation configuration. This is implemented as a member + * function so as to better facilitate changes. */ + void setup(const GeneratorFileOrCmdParam& param, + const conf::SimConfig& config); + /** + * Set command to execute in bacground rather than reading from + * existing file(s) + * + * @param cmd Command line. Can include options for the program to + * execute, but should not include pipes. + */ + void setCmd(const std::string& cmd) { mCmd = cmd; } + /** + * Set the number of events that a background command should + * generate. This should come from @c SimConfig::getNEents. + * + * @param nev Number of events to generate. This is passed via @c + * mNEventsSwitch to the command line. + */ + void setNEvents(unsigned int nev) { mNEvents = nev; } + /** + * Set the random number seed that a background command should use. + * This should come from @c SimConfig::getStartSeed + * + * @param seed Random number seed. Will be passed to the + * commandline using the @c mSeedSwitch. + */ + void setSeed(unsigned long seed) { mSeed = seed; } + /** + * Set the maximum impact parameter to sample by a background + * command. This should come from @c SimConfig::getBMax() + * + * @param bmax Maximum impact parameter, in fm, to sample. This is + * passed to the command line via the @c mBmaxSwitch. + */ + void setBmax(float bmax) { mBmax = bmax; } + /** + * Set the file names to read from. + * + * @param filenames A comma seperated list of files to read from. + */ + void setFileNames(const std::string& filenames); + /** + * Set the output switch. + * + * @param opt Command line switch (e.g., @c -o) to specify output + * file name. */ + void setOutputSwitch(const std::string& opt) { mOutputSwitch = opt; } + /** + * Set the seed switch. + * + * @param opt Command line switch (e.g., @c -s) to specify the + * random number seed to use when generating events. + */ + void setSeedSwitch(const std::string& opt) { mSeedSwitch = opt; } + /** + * Set the nevents switch. + * + * @param opt Command line switch (e.g., @c -n) to specify the + * number of events to generate. + */ + void setNEventsSwitch(const std::string& opt) { mNEventsSwitch = opt; } + /** + * Set the maximum impact parameter switch. + * + * @param opt Command line switch (e.g., @c -b) to specify the + * maximum impact parameter (in fm) to sample when generating + * events. + */ + void setBmaxSwitch(const std::string& opt) { mBmaxSwitch = opt; } + /** + * Set the background switch. + * + * @param opt Command line switch (e.g., @c &) to detach and send + * the event generator program into the background. + */ + void setBackgroundSwitch(const std::string& opt) { mBackgroundSwitch = opt; } + /** Set the wait time, in miliseconds, when waiting for data */ + void setWait(int miliseconds = 500) { mWait = miliseconds; } + + protected: + /** + * Format a command line using the set command line, option flags, + * and option values. + * + * @return formatted command line. + */ + virtual std::string makeCmdLine() const; + /** + * Execute a command line (presumably formatted by @c makeCmdLine). + * If the command failed to execute, then make it a fatal error. + * + * @param cmd Command line to execute, presumabley formatted by @c + * makeCmdLine. + + * @return true if the background command line was executed, false + * otherwise. + */ + virtual bool executeCmdLine(const std::string& cmd) const; + /** + * Create a temporary file (and close it immediately). On success, + * the list of file names is cleared and the name of the temporary + * file set as the sole element in that list. + * + * @return true if the temporary file name was generated + * successfully. + */ + virtual bool makeTemp(); + /** + * Remove the temporary file if it was set and it exists. + * + * @return true if the temporary file was removed. + */ + virtual bool removeTemp() const; + /** + * Make a fifo at the location of the first element of the list of + * file names (presumably a temporary file as created by + * makeTemp). + * + * @return true if the FIFo was made correctly + */ + virtual bool makeFifo() const; + /** + * Ensure that all files in the list of file names exists. If @e + * any of te file names point to a net resource (e.g., @c alien://or + * @c https://) then this member function should @e not be called + * + * The file names registered are replaced with their canonical path + * equivalents. + * + * @return true if all currently registered file names can be found. + */ + virtual bool ensureFiles(); + /** + * Wait for data to be available in case we're executing a + * background command + */ + virtual void waitForData(const std::string& filename) const; + /** + * Possible command line to execute. The command executed must + * accept the switches defined below. Note if @c mOutputSwitch is + * set to @c ">", then the program @e must write data to standard + * output + */ + std::string mCmd = ""; + /** + * List of file names to read. In case we're executing a command, + * then there will only be one file name in the list + */ + std::list mFileNames; + /** + * Name of temporary file, if it was created. + */ + std::string mTemporary; + /** + * Number of events to generate in case we're executing a command. + * This is passed to the program via the switch @c + * mNEventsSwitch + */ + unsigned int mNEvents = 0; + /** + * Random number seed to use in case we're executing a command. + * This is passed to the program via the switch @c + * mSeedSwitch + */ + unsigned long mSeed = 0; + /** + * Maximum impact parameter to sample in case we're executing a + * command. IF negative, then it is not passed to the command. + * This is passed to the command via the switch @c mBmaxSwitch + */ + float mBmax = -1.f; + /** + * Switch to direct output to specified file. In case of a fifo, + * this should often be @c ">" - i.e., the standard output of the + * program is redirected to the fifo. + */ + std::string mOutputSwitch = ">"; + /** + * The switch specify the random number seed to the program + * executed + */ + std::string mSeedSwitch = "-s"; + /** + * The switch to specify the number of events to generate to the + * program being executed + */ + std::string mNEventsSwitch = "-n"; + /** + * The switch to specify maximum impact parameter to sample by the + * program being executed. + */ + std::string mBmaxSwitch = "-b"; + /** + * The "switch" to put the program being executed in the + * background. + */ + std::string mBackgroundSwitch = "&"; + /** + * Time in miliseconds between each wait for data + */ + int mWait = 500; +}; + +} // namespace eventgen +} // namespace o2 +#endif +// Local Variables: +// mode: C++ +// End: diff --git a/Generators/include/Generators/GeneratorFileOrCmdParam.h b/Generators/include/Generators/GeneratorFileOrCmdParam.h new file mode 100644 index 0000000000000..fe6dfa3a80722 --- /dev/null +++ b/Generators/include/Generators/GeneratorFileOrCmdParam.h @@ -0,0 +1,45 @@ +// 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_GENERATORFILEORCMDPARAM_H_ +#define ALICEO2_EVENTGEN_GENERATORFILEORCMDPARAM_H_ + +#include "CommonUtils/ConfigurableParam.h" +#include "CommonUtils/ConfigurableParamHelper.h" +#include + +namespace o2 +{ +namespace eventgen +{ + +/** + ** a parameter class/struct to keep the settings of + ** the Fileorcmd event generator and + ** allow the user to modify them + **/ +struct GeneratorFileOrCmdParam : public o2::conf::ConfigurableParamHelper { + std::string fileNames = ""; + std::string cmd = ""; // Program command line to spawn + std::string outputSwitch = ">"; + std::string seedSwitch = "-s"; + std::string bMaxSwitch = "-b"; + std::string nEventsSwitch = "-n"; + std::string backgroundSwitch = "&"; + O2ParamDef(GeneratorFileOrCmdParam, "GeneratorFileOrCmd"); +}; + +} // end namespace eventgen +} // end namespace o2 + +#endif // ALICEO2_EVENTGEN_GENERATORFILEORCMDPARAM_H_ diff --git a/Generators/include/Generators/GeneratorHepMC.h b/Generators/include/Generators/GeneratorHepMC.h index 341f4b8a0e38a..9a5c607f209da 100644 --- a/Generators/include/Generators/GeneratorHepMC.h +++ b/Generators/include/Generators/GeneratorHepMC.h @@ -15,7 +15,8 @@ #define ALICEO2_EVENTGEN_GENERATORHEPMC_H_ #include "Generators/Generator.h" -#include +#include "Generators/GeneratorFileOrCmd.h" +#include "Generators/GeneratorHepMCParam.h" #ifdef GENERATORS_WITH_HEPMC3_DEPRECATED namespace HepMC @@ -35,33 +36,51 @@ class FourVector; namespace o2 { +namespace conf +{ +class SimConfig; +} namespace eventgen { /*****************************************************************/ /*****************************************************************/ -class GeneratorHepMC : public Generator +class GeneratorHepMC : public Generator, public GeneratorFileOrCmd { public: /** default constructor **/ GeneratorHepMC(); /** constructor **/ - GeneratorHepMC(const Char_t* name, const Char_t* title = "ALICEo2 HepMC Generator"); + GeneratorHepMC(const Char_t* name, + const Char_t* title = "ALICEo2 HepMC Generator"); /** destructor **/ ~GeneratorHepMC() override; - /** Initialize the generator if needed **/ + /** Initialize the generator. **/ Bool_t Init() override; - /** methods to override **/ + /** + * Configure the generator from parameters and the general + * simulation configuration. This is implemented as a member + * function so as to better facilitate changes. */ + void setup(const GeneratorFileOrCmdParam& param0, + const GeneratorHepMCParam& param, + const conf::SimConfig& config); + /** + * Generate a single event. The event is read in from the current + * input file. Returns false if a new event could not be read. + **/ Bool_t generateEvent() override; + /** + * Import particles from the last read event into a vector + * TParticle. Returns false if no particles could be exported to + * the vector. + */ Bool_t importParticles() override; /** setters **/ - void setVersion(Int_t val) { mVersion = val; }; - void setFileName(std::string val) { mFileName = val; }; void setEventsToSkip(uint64_t val) { mEventsToSkip = val; }; protected: @@ -77,18 +96,17 @@ 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; + /** Make our reader */ + bool makeReader(); + /** HepMC interface **/ - std::ifstream mStream; //! - std::string mFileName; - Int_t mVersion; - uint64_t mEventsToSkip; -#ifdef GENERATORS_WITH_HEPMC3_DEPRECATED - HepMC::Reader* mReader; //! - HepMC::GenEvent* mEvent; //! -#else - HepMC3::Reader* mReader; //! - HepMC3::GenEvent* mEvent; //! -#endif + uint64_t mEventsToSkip = 0; + int mVersion = 0; + std::shared_ptr mReader; + /** Event structure */ + HepMC3::GenEvent* mEvent = nullptr; ClassDefOverride(GeneratorHepMC, 1); diff --git a/Generators/include/Generators/GeneratorHepMCParam.h b/Generators/include/Generators/GeneratorHepMCParam.h index 4c53918c70634..5df4013489f95 100644 --- a/Generators/include/Generators/GeneratorHepMCParam.h +++ b/Generators/include/Generators/GeneratorHepMCParam.h @@ -30,9 +30,9 @@ namespace eventgen **/ struct GeneratorHepMCParam : public o2::conf::ConfigurableParamHelper { - std::string fileName = ""; - int version = 2; + int version = 0; uint64_t eventsToSkip = 0; + std::string fileName = ""; O2ParamDef(GeneratorHepMCParam, "HepMC"); }; diff --git a/Generators/include/Generators/GeneratorTParticle.h b/Generators/include/Generators/GeneratorTParticle.h new file mode 100644 index 0000000000000..e4ddb5fa1f340 --- /dev/null +++ b/Generators/include/Generators/GeneratorTParticle.h @@ -0,0 +1,119 @@ +// 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 +#include + +// Forward decls +class TChain; +class TParticle; +class TClonesArray; + +namespace o2 +{ +namespace conf +{ +class SimConfig; +} +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 GeneratorFileOrCmd +{ + public: + /** CTOR */ + GeneratorTParticle(); + /** 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; + + /** + * Configure the generator from parameters and the general + * simulation configuration. This is implemented as a member + * function so as to better facilitate changes. */ + void setup(const GeneratorFileOrCmdParam& param0, + const GeneratorTParticleParam& param, + const conf::SimConfig& config); + /** 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 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; } + + protected: + /** Name of the tree to read */ + std::string mTreeName = "T"; + /** Name of branch containing a TClonesArray of TParticle */ + std::string mBranchName = "Particles"; + /** Current entry */ + unsigned int mEntry = 0; + /** Chain of files */ + TChain* mChain = 0; + /** Array to read particles into */ + TClonesArray* mTParticles; + + 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..98059f3b0e1e6 --- /dev/null +++ b/Generators/include/Generators/GeneratorTParticleParam.h @@ -0,0 +1,39 @@ +// 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"; + 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..6780315974c2c 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,28 @@ 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& param0 = GeneratorFileOrCmdParam::Instance(); + auto& param = GeneratorTParticleParam::Instance(); + LOG(info) << "Init 'GeneratorTParticle' with the following parameters"; + LOG(info) << param0; + LOG(info) << param; + auto tgen = new o2::eventgen::GeneratorTParticle(); + tgen->setup(param0, param, conf); + 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& param0 = GeneratorFileOrCmdParam::Instance(); auto& param = GeneratorHepMCParam::Instance(); LOG(info) << "Init \'GeneratorHepMC\' with following parameters"; + LOG(info) << param0; LOG(info) << param; auto hepmcGen = new o2::eventgen::GeneratorHepMC(); - hepmcGen->setFileName(param.fileName); - hepmcGen->setVersion(param.version); - hepmcGen->setEventsToSkip(param.eventsToSkip); + hepmcGen->setup(param0, param, conf); primGen->AddGenerator(hepmcGen); #endif #ifdef GENERATORS_WITH_PYTHIA6 diff --git a/Generators/src/GeneratorFileOrCmd.cxx b/Generators/src/GeneratorFileOrCmd.cxx new file mode 100644 index 0000000000000..9ff45bcd9c867 --- /dev/null +++ b/Generators/src/GeneratorFileOrCmd.cxx @@ -0,0 +1,214 @@ +// Copyright 2019-2020 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 "SimulationDataFormat/MCUtils.h" +#include "Generators/GeneratorFileOrCmd.h" +// For fifo's and system call +#include +#include // POSIX only +#include // POISX only +#include +// For filesystem operations +#include +// Waits +#include +// Log messages +#include +#include + +namespace +{ +std::string ltrim(const std::string& s, const std::string& what = "\" ' ") +{ + auto start = s.find_first_not_of(what); + if (start == std::string::npos) { + return ""; + } + + return s.substr(start); +} + +std::string rtrim(const std::string& s, const std::string& what = "\"' ") +{ + auto end = s.find_last_not_of(what); + return s.substr(0, end + 1); +} + +std::string trim(const std::string& s, const std::string& what = "\"' ") +{ + return rtrim(ltrim(s, what), what); +} +} // namespace + +namespace o2 +{ +namespace eventgen +{ +// ----------------------------------------------------------------- +void GeneratorFileOrCmd::setup(const GeneratorFileOrCmdParam& param, + const conf::SimConfig& config) +{ + setFileNames(param.fileNames); + setCmd(param.cmd); + setOutputSwitch(trim(param.outputSwitch)); + setSeedSwitch(trim(param.seedSwitch)); + setBmaxSwitch(trim(param.bMaxSwitch)); + setNEventsSwitch(trim(param.nEventsSwitch)); + setBackgroundSwitch(trim(param.backgroundSwitch)); + setSeed(config.getStartSeed()); + setNEvents(config.getNEvents()); + setBmax(config.getBMax()); +} +// ----------------------------------------------------------------- +void GeneratorFileOrCmd::setFileNames(const std::string& filenames) +{ + std::stringstream s(filenames); + std::string f; + while (std::getline(s, f, ',')) { + mFileNames.push_back(f); + } +} +// ----------------------------------------------------------------- +std::string GeneratorFileOrCmd::makeCmdLine() const +{ + std::string fileName = mFileNames.front(); + std::stringstream s; + s << mCmd << " "; + if (not mSeedSwitch.empty() and mSeedSwitch != "none") { + s << mSeedSwitch << " " << mSeed << " "; + } + if (not mNEventsSwitch.empty() and mNEventsSwitch != "none") { + s << mNEventsSwitch << " " << mNEvents << " "; + } + if (not mBmaxSwitch.empty() and mBmax >= 0 and mBmaxSwitch != "none") { + s << mBmaxSwitch.c_str() << " " << mBmax << " "; + } + + s << mOutputSwitch << " " << fileName << " " + << mBackgroundSwitch; + return s.str(); +} +// ----------------------------------------------------------------- +bool GeneratorFileOrCmd::executeCmdLine(const std::string& cmd) const +{ + LOG(info) << "Command line to execute: \"" << cmd << "\""; + int ret = std::system(cmd.c_str()); + if (ret != 0) { + LOG(fatal) << "Failed to spawn \"" << cmd << "\""; + return false; + } + return true; +} +// ----------------------------------------------------------------- +bool GeneratorFileOrCmd::makeTemp() +{ + mFileNames.clear(); + char buf[] = "generatorFifoXXXXXX"; + auto fp = mkstemp(buf); + if (fp < 0) { + LOG(fatal) << "Failed to make temporary file: " + << std::strerror(errno); + return false; + } + mTemporary = std::string(buf); + mFileNames.push_back(mTemporary); + close(fp); + return true; +} +// ----------------------------------------------------------------- +bool GeneratorFileOrCmd::removeTemp() const +{ + if (mTemporary.empty()) { + LOG(info) << "Temporary file name empty, nothing to remove"; + return false; + } + + // Get the file we're reading from + std::filesystem::path p(mTemporary); + + // Check if the file exists + if (not std::filesystem::exists(p)) { + LOG(info) << "Temporary file " << p << " does not exist"; + return true; + } + + // Remove temporary file + std::error_code ec; + std::filesystem::remove(p, ec); + if (ec) { + LOG(warn) << "When removing " << p << ": " << ec.message(); + } + + // Ignore errors when removing the temporary file + return true; +} +// ----------------------------------------------------------------- +bool GeneratorFileOrCmd::makeFifo() const +{ + // First remove the temporary file if it exists, + // otherwise we may not be able to make the FIFO + removeTemp(); + + std::string fileName = mFileNames.front(); + + int ret = mkfifo(fileName.c_str(), 0600); + if (ret != 0) { + LOG(fatal) << "Failed to make fifo \"" << fileName << "\": " + << std::strerror(errno); + return false; + } + + return true; +} +// ----------------------------------------------------------------- +bool GeneratorFileOrCmd::ensureFiles() +{ + try { + for (auto& f : mFileNames) { + auto c = std::filesystem::canonical(std::filesystem::path(f)); + f = c.c_str(); + } + } catch (std::exception& e) { + LOG(error) << e.what(); + return false; + } + return true; +} +// ----------------------------------------------------------------- +void GeneratorFileOrCmd::waitForData(const std::string& filename) const +{ + using namespace std::chrono_literals; + + // Get the file we're reading from + std::filesystem::path p(filename); + + LOG(debug) << "Waiting for data on " << p; + + // Wait until child process creates the file + while (not std::filesystem::exists(p)) { + std::this_thread::sleep_for(mWait * 1ms); + } + + // 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(mWait * 1ms); + } + + // Give the child process 1 second to post the data to the file + LOG(debug) << "Got data in " << p << ", sleeping for a while"; + std::this_thread::sleep_for(mWait * 2ms); +} + +} /* namespace eventgen */ +} /* namespace o2 */ diff --git a/Generators/src/GeneratorFileOrCmdParam.cxx b/Generators/src/GeneratorFileOrCmdParam.cxx new file mode 100644 index 0000000000000..e452f272c3735 --- /dev/null +++ b/Generators/src/GeneratorFileOrCmdParam.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/GeneratorFileOrCmdParam.h" +O2ParamImpl(o2::eventgen::GeneratorFileOrCmdParam); +// +// EOF +// diff --git a/Generators/src/GeneratorHepMC.cxx b/Generators/src/GeneratorHepMC.cxx index 6c062268bfacc..3fec13ae490d2 100644 --- a/Generators/src/GeneratorHepMC.cxx +++ b/Generators/src/GeneratorHepMC.cxx @@ -14,14 +14,15 @@ #include "SimulationDataFormat/MCUtils.h" #include "Generators/GeneratorHepMC.h" #include "Generators/GeneratorHepMCParam.h" -#include "HepMC3/ReaderAscii.h" -#include "HepMC3/ReaderAsciiHepMC2.h" +#include "SimulationDataFormat/MCEventHeader.h" +#include "SimConfig/SimConfig.h" +#include "HepMC3/ReaderFactory.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 "FairPrimaryGenerator.h" @@ -36,19 +37,14 @@ namespace eventgen /*****************************************************************/ GeneratorHepMC::GeneratorHepMC() - : Generator("ALICEo2", "ALICEo2 HepMC Generator"), mStream(), mFileName(), mVersion(3), mReader(nullptr), mEvent(nullptr) + : GeneratorHepMC("ALICEo2", "ALICEo2 HepMC Generator") { - /** default constructor **/ - - mEvent = new HepMC3::GenEvent(); - mInterface = reinterpret_cast(mEvent); - mInterfaceName = "hepmc"; } /*****************************************************************/ GeneratorHepMC::GeneratorHepMC(const Char_t* name, const Char_t* title) - : Generator(name, title), mStream(), mFileName(), mVersion(3), mReader(nullptr), mEvent(nullptr) + : Generator(name, title) { /** constructor **/ @@ -62,36 +58,66 @@ GeneratorHepMC::GeneratorHepMC(const Char_t* name, const Char_t* title) GeneratorHepMC::~GeneratorHepMC() { /** default destructor **/ - - if (mStream.is_open()) { - mStream.close(); - } + LOG(info) << "Destructing GeneratorHepMC"; if (mReader) { mReader->close(); - delete mReader; } if (mEvent) { delete mEvent; } + removeTemp(); } /*****************************************************************/ +void GeneratorHepMC::setup(const GeneratorFileOrCmdParam& param0, + const GeneratorHepMCParam& param, + const conf::SimConfig& config) +{ + if (not param.fileName.empty()) { + LOG(warn) << "The use of the key \"HepMC.fileName\" is " + << "deprecated, use \"GeneratorFileOrCmd.fileNames\" instead"; + } + + GeneratorFileOrCmd::setup(param0, config); + if (not param.fileName.empty()) { + setFileNames(param.fileName); + } + mVersion = param.version; + setEventsToSkip(param.eventsToSkip); + + if (param.version != 0 and mCmd.empty()) { + LOG(warn) << "The key \"HepMC.version\" is no longer used when " + << "reading from files. The format version of the input files " + << "are automatically deduced."; + } +} + +/*****************************************************************/ Bool_t GeneratorHepMC::generateEvent() { + LOG(debug) << "Generating an event"; /** generate event **/ + int tries = 0; + do { + LOG(debug) << " try # " << ++tries; + if (not mReader and not makeReader()) { + return false; + } - /** clear and read event **/ - mEvent->clear(); - mReader->read_event(*mEvent); - if (mReader->failed()) { - return kFALSE; - } - /** set units to desired output **/ - mEvent->set_units(HepMC3::Units::GEV, HepMC3::Units::MM); + /** clear and read event **/ + mEvent->clear(); + mReader->read_event(*mEvent); + if (not mReader->failed()) { + /** set units to desired output **/ + mEvent->set_units(HepMC3::Units::GEV, HepMC3::Units::MM); + LOG(debug) << "Read one event " << mEvent->event_number(); + return true; + } + } while (true); - /** success **/ - return kTRUE; + /** failure **/ + return false; } /*****************************************************************/ @@ -143,6 +169,213 @@ Bool_t GeneratorHepMC::importParticles() return kTRUE; } +namespace +{ +template +bool putAttributeInfoImpl(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()); + return true; + } + return false; +} + +void putAttributeInfo(o2::dataformats::MCEventHeader* eventHeader, + const std::string& name, + const std::shared_ptr& a) +{ + using IntAttribute = HepMC3::IntAttribute; + using LongAttribute = HepMC3::LongAttribute; + using FloatAttribute = HepMC3::FloatAttribute; + using DoubleAttribute = HepMC3::DoubleAttribute; + using StringAttribute = HepMC3::StringAttribute; + using CharAttribute = HepMC3::CharAttribute; + using LongLongAttribute = HepMC3::LongLongAttribute; + using LongDoubleAttribute = HepMC3::LongDoubleAttribute; + using UIntAttribute = HepMC3::UIntAttribute; + using ULongAttribute = HepMC3::ULongAttribute; + using ULongLongAttribute = HepMC3::ULongLongAttribute; + using BoolAttribute = HepMC3::BoolAttribute; + + if (putAttributeInfoImpl(eventHeader, name, a)) { + return; + } + if (putAttributeInfoImpl(eventHeader, name, a)) { + return; + } + if (putAttributeInfoImpl(eventHeader, name, a)) { + return; + } + if (putAttributeInfoImpl(eventHeader, name, a)) { + return; + } + if (putAttributeInfoImpl(eventHeader, name, a)) { + return; + } + if (putAttributeInfoImpl(eventHeader, name, a)) { + return; + } + if (putAttributeInfoImpl(eventHeader, name, a)) { + return; + } + if (putAttributeInfoImpl(eventHeader, name, a)) { + return; + } + if (putAttributeInfoImpl(eventHeader, name, a)) { + return; + } + if (putAttributeInfoImpl(eventHeader, name, a)) { + return; + } + if (putAttributeInfoImpl(eventHeader, name, a)) { + return; + } + if (putAttributeInfoImpl(eventHeader, name, a)) { + return; + } +} +} // 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 GeneratorHepMC::makeReader() +{ + // Reset the reader smart pointer + LOG(debug) << "Reseting the reader"; + mReader.reset(); + + // Check that we have any file names left + if (mFileNames.size() < 1) { + LOG(debug) << "No more files to read, return false"; + return false; + } + + // If we have file names left, pop the top of the list (LIFO) + auto filename = mFileNames.front(); + mFileNames.pop_front(); + + LOG(debug) << "Next file to read: \"" << filename << "\" " + << mFileNames.size() << " left"; + + if (not mCmd.empty()) { + // For FIFO reading, we assume straight ASCII output always. + // Unfortunately, the HepMC3::deduce_reader `stat`s the filename + // which isn't supported on a FIFO, so we have to use the reader + // directly. Here, we allow for version 2 formats if the user + // specifies that + LOG(info) << "Creating ASCII reader of " << filename; + if (mVersion == 2) { + mReader.reset(new HepMC3::ReaderAsciiHepMC2(filename)); + } else { + mReader.reset(new HepMC3::ReaderAscii(filename)); + } + } else { + LOG(info) << "Deduce a reader of " << filename; + mReader = HepMC3::deduce_reader(filename); + } + + bool ret = bool(mReader) and not mReader->failed(); + LOG(info) << "Reader is " << mReader.get() << " " << ret; + return ret; +} + /*****************************************************************/ Bool_t GeneratorHepMC::Init() @@ -152,41 +385,81 @@ Bool_t GeneratorHepMC::Init() /** init base class **/ Generator::Init(); - /** open file **/ - std::string filename = gSystem->ExpandPathName(mFileName.c_str()); - mStream.open(filename); - if (!mStream.is_open()) { - LOG(fatal) << "Cannot open input file: " << filename << std::endl; - return kFALSE; - } - - /** create reader according to HepMC version **/ - switch (mVersion) { - case 2: - mStream.close(); - mReader = new HepMC3::ReaderAsciiHepMC2(filename); - break; - case 3: - mReader = new HepMC3::ReaderAscii(mStream); - break; - default: - LOG(fatal) << "Unsupported HepMC version: " << mVersion << std::endl; - return kFALSE; - } - - // skip events at the beginning - if (!mReader->failed()) { - LOGF(info, "%i events to skip.", mEventsToSkip); - for (auto ind = 0; ind < mEventsToSkip; ind++) { - if (!generateEvent()) { - LOGF(error, "The file %s only contains %i events!", mFileName, ind); - break; - } + // 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 + // following command line argument + // + // `-n NEVENTS` to set the number of events to produce. + // + // Optionally, the command line should also accept + // + // `-s SEED` to set the random number seed + // `-b FM` to set the maximum impact parameter to sample + // `-o OUTPUT` to set the output file name + // + // All of this can conviniently be achieved via a wrapper script + // around the actual EG program. + if (not mCmd.empty()) { + // Set filename to be a temporary name + if (not makeTemp()) { + return false; + } + + // Make a fifo + if (not makeFifo()) { + return false; + } + + // Build command line, rediret stdout to our fifo and put + std::string cmd = makeCmdLine(); + LOG(debug) << "EG command line is \"" << cmd << "\""; + + // Execute the command line + if (not executeCmdLine(cmd)) { + LOG(fatal) << "Failed to spawn \"" << cmd << "\""; + return false; + } + } else { + // If no command line was given, ensure that all files are present + // on the system. Note, in principle, HepMC3 can read from remote + // files + // + // root:// XRootD served + // http[s]:// Web served + // gsidcap:// DCap served + // + // These will all be handled in HepMC3 via ROOT's TFile protocol + // and the files are assumed to contain a TTree named + // `hepmc3_tree` and that tree has the branches + // + // `hepmc3_event` with object of type `HepMC3::GenEventData` + // `GenRunInfo` with object of type `HepMC3::GenRunInfoData` + // + // where the last branch is optional. + // + // However, here we will assume system local files. If _any_ of + // the listed files do not exist, then we fail. + if (not ensureFiles()) { + return false; } } - /** success **/ - return !mReader->failed(); + // Create reader for current (first) file + return true; } /*****************************************************************/ diff --git a/Generators/src/GeneratorPythia8.cxx b/Generators/src/GeneratorPythia8.cxx index 13c17003ee612..7038647ee4ba0 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,17 @@ 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..ab68f7f39b1bf --- /dev/null +++ b/Generators/src/GeneratorTParticle.cxx @@ -0,0 +1,154 @@ +// 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 + +namespace o2 +{ +namespace eventgen +{ +/*****************************************************************/ +GeneratorTParticle::GeneratorTParticle() +{ + setOutputSwitch("-o"); +} + +/*****************************************************************/ +GeneratorTParticle::~GeneratorTParticle() +{ + if (mChain) { + TFile* file = mChain->GetCurrentFile(); + if (file) { + mChain->RecursiveRemove(file); + } + delete mChain; + } + if (mCmd.empty()) { + return; + } + + removeTemp(); +} +/*****************************************************************/ +Bool_t GeneratorTParticle::Init() +{ + mChain = new TChain(mTreeName.c_str()); + mTParticles = new TClonesArray("TParticle"); + mChain->SetBranchAddress(mBranchName.c_str(), &mTParticles); + + if (not mCmd.empty()) { + // Set filename to be a temporary name + if (not makeTemp()) { + return false; + } + + // Build command line, Assumes command line parameter + std::string cmd = makeCmdLine(); + LOG(info) << "EG command line is \"" << cmd << "\""; + + // Execute the background command + if (not executeCmdLine(cmd)) { + LOG(fatal) << "Failed to spawn \"" << cmd << "\""; + return false; + } + } + for (auto filename : mFileNames) { + mChain->AddFile(filename.c_str()); + } + + // Clear the array of file names + mFileNames.clear(); + + return true; +} + +/*****************************************************************/ +void GeneratorTParticle::setup(const GeneratorFileOrCmdParam& param0, + const GeneratorTParticleParam& param, + const conf::SimConfig& config) +{ + GeneratorFileOrCmd::setup(param0, config); + setTreeName(param.treeName); + setBranchName(param.branchName); +} + +/*****************************************************************/ +Bool_t GeneratorTParticle::generateEvent() +{ + // If this is the first entry, and we're executing a command, then + // wait until the input file exists and actually contain some data. + if (mEntry == 0 and not mCmd.empty()) { + waitForData(mTemporary); + } + + // Read in the next entry in the chain + int read = mChain->GetEntry(mEntry); + mEntry++; + + // If we got an error while reading, then give error message + if (read < 0) { + LOG(error) << "Failed to read entry " << mEntry << " of chain"; + } + + // If we had an error or nothing was read back, then return false + 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..35e77febee172 100644 --- a/Generators/src/GeneratorsLinkDef.h +++ b/Generators/src/GeneratorsLinkDef.h @@ -67,5 +67,10 @@ #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> + ; +#pragma link C++ class o2::eventgen::GeneratorFileOrCmdParam + ; +#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::eventgen::GeneratorFileOrCmdParam> + ; #endif diff --git a/run/CMakeLists.txt b/run/CMakeLists.txt index 07f23328230dd..f9b0f5e220ea6 100644 --- a/run/CMakeLists.txt +++ b/run/CMakeLists.txt @@ -271,7 +271,7 @@ o2_add_test_command(NAME o2sim_hepmc -g hepmc --configKeyValues -"HepMC.fileName=${CMAKE_SOURCE_DIR}/Generators/share/data/pythia.hepmc;HepMC.version=2;align-geom.mDetectors=none" +"GeneratorFileOrCmd.fileNames=${CMAKE_SOURCE_DIR}/Generators/share/data/pythia.hepmc;HepMC.version=2;align-geom.mDetectors=none" -o o2simhepmc LABELS long sim hepmc3 diff --git a/run/SimExamples/HepMC/README.md b/run/SimExamples/HepMC/README.md new file mode 100644 index 0000000000000..0369ee9d90bdd --- /dev/null +++ b/run/SimExamples/HepMC/README.md @@ -0,0 +1,169 @@ + + +Here are pointers on how to use `GeneratorHepMC` selected by the +option `-g hepmc` for `o2-sim`. + +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). + +## 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 "GeneratorFileOrCmd.fileNames=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 "GeneratorFileOrCmd.cmd=eg" + +See also [`child.sh`](child.sh). + +There are some requirements on the program `eg`: + +- The EG program _must_ be able to write the HepMC event structures to + a specified file. The option passed to the program is specified via + the key `GeneratorFileOrCmd.outputSwitch`. This defaults to `>` + which means the EG program is assumed to write the HepMC event + structures to standard output, _and_ that nothing else is printed on + standard output. +- It _must_ accept an option to set the number of events to generate. + This is controlled by the configuration key + `GeneratorFileOrCmd.nEventsSwitch` and defaults to `-n`. Thus, the + EG application should accept `-n 10` to mean that it should generate + `10` events, for example. +- The EG application should accept a command line switch to set the + random number generator seed. This option is specified via the + configuration key `GeneratorFileOrCmd.seedSwitch` and defaults to + `-s`. Thus, the EG application must accept `-s 123456` to mean to + set the random number seed to `123456` for example. +- The EG application should accept a command line switch to set the + maximum impact parameter (in Fermi-metre) sampled. This is set via + the configuration key `GeneratorFileOrCmd.bMaxSwithc` and defaults + to `-b`. Thus, the EG application should take the command line + argument `-b 10` to mean that it should only generate events with an + impact parameter between 0fm and 10fm. + +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 "GeneratorFileOrCmd.cmd=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 configuration keys + +The `GeneratorHepMC` (and sister generator `GeneratorTParticle`) +allows customisation of the execution via configuration keys passed +via `--configKeyValues` + +- `HepMC.eventsToSkip=number` a number events to skip at the beginning + of each file read. + +- `HepMC.version` - when reading the events from files, this option is + no longer needed. The code itself figures out which format version + the input file is in. If executing a child process through + `GeneratorFileOrCmd.cmd` and the EG writes out HepMC2 format, then + this _must_ be set to `2`. Otherwise, HepMC3 is assumed. + +- `GeneratorFileOrCmd.fileNames=list` a comma separated list of HepMC + files to read. + +- `GeneratorFileOrCmd.cmd=command line` a command line to execute as a + background child process. If this is set (not the empty string), + then `GeneratorFileOrCmd.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. + + - `GeneratorFileOrCmd.outputSwitch=switch` (default `>`) to specify + output file. The default of `>` assumes that the program write + HepMC events, and _only_ those, to standard output. + + - `GeneratorFileOrCmd.seedSwitch=switch` (default `-s`) to specify + the random number generator seed. The value passed is selected by + the `o2-sim` option `--seed` + + - `GeneratorFileOrCmd.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` + + - `GeneratorFileOrCmd.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`) + + - `GeneratorFileOrCmd.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 deprecated + + - `HepMC.fileName` - use `GeneratorFileOrCmd.fileNames` + +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..89b563ec8d550 --- /dev/null +++ b/run/SimExamples/HepMC/child.sh @@ -0,0 +1,53 @@ +#!/usr/bin/env bash + +cmd="./crmc.sh" +more="GeneratorFileOrCmd.bMaxSwitch=none" +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 "GeneratorFileOrCmd.cmd=$cmd;${more}" \ + --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..1c2d2f7827dbc --- /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 hepmc -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..e64f5561a775c --- /dev/null +++ b/run/SimExamples/HepMC/read.sh @@ -0,0 +1,48 @@ +#!/usr/bin/env 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 "GeneratorFileOrCmd.fileNames=$inp" \ + --outPrefix "$out" --seed $seed --nEvents $nev $@ diff --git a/run/SimExamples/HepMC_STARlight/run.sh b/run/SimExamples/HepMC_STARlight/run.sh index b75cee52b7117..3eb5b91f24a6d 100755 --- a/run/SimExamples/HepMC_STARlight/run.sh +++ b/run/SimExamples/HepMC_STARlight/run.sh @@ -10,9 +10,9 @@ set -x # PART a) env -i HOME="$HOME" USER="$USER" PATH="/bin:/usr/bin:/usr/local/bin" \ - ALIBUILD_WORK_DIR="$ALIBUILD_WORK_DIR" ./run-starlight.sh + ALIBUILD_WORK_DIR="$ALIBUILD_WORK_DIR" ./run-starlight.sh # PART b) NEV=1000 o2-sim -j 20 -n ${NEV} -g hepmc -m PIPE ITS -o sim \ - --configKeyValues "HepMC.fileName=starlight.hepmc;HepMC.version=2;Diamond.position[2]=0.1;Diamond.width[2]=0.05" + --configKeyValues "GeneratorFileOrCmd.fileNames=starlight.hepmc;Diamond.position[2]=0.1;Diamond.width[2]=0.05" 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..c9df40665d829 --- /dev/null +++ b/run/SimExamples/TParticle/README.md @@ -0,0 +1,265 @@ + + +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 "GeneratorFileOrCmd.fileNames=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 "GeneratorFileOrCmd.fileNames=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 "GeneratorFileOrCmd.cmd=eg" + +See also [`child.sh`](child.sh). Do + + ./child.sh --help + +for a list of options. + +There are some requirements on the program `eg`: + +- The EG program _must_ be able to write the HepMC event structures to + a specified file. The option passed to the program is specified via + the key `GeneratorFileOrCmd.outputSwitch`. This defaults to `-o`. +- It _must_ accept an option to set the number of events to generate. + This is controlled by the configuration key + `GeneratorFileOrCmd.nEventsSwitch` and defaults to `-n`. Thus, the + EG application should accept `-n 10` to mean that it should generate + `10` events, for example. +- The EG application should accept a command line switch to set the + random number generator seed. This option is specified via the + configuration key `GeneratorFileOrCmd.seedSwitch` and defaults to + `-s`. Thus, the EG application must accept `-s 123456` to mean to + set the random number seed to `123456` for example. +- The EG application should accept a command line switch to set the + maximum impact parameter (in Fermi-metre) sampled. This is set via + the configuration key `GeneratorFileOrCmd.bMaxSwithc` and defaults + to `-b`. Thus, the EG application should take the command line + argument `-b 10` to mean that it should only generate events with an + impact parameter between 0fm and 10fm. + +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 "GeneratorFileOrCmd.cmd=./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`) is +configured through configuration keys set via `--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. + +- `GeneratorFileOrCmd.fileNames=list` a comma separated list of HepMC + files to read + +- `GeneratorFileOrCmd.cmd=command line` a command line to execute as a + background child process. If this is set (not the empty string), + then `GeneratorFileOrCmd.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 or special value `none`, then that switch + and corresponding option value is not passed to the child program. + + - `GeneratorFileOrCmd.outputSwitch=switch` (default `>`) to specify + output file. The default of `>` assumes that the program write + events, and _only_ those, to standard output. + + - `GeneratorFileOrCmd.seedSwitch=switch` (default `-s`) to specify + the random number generator seed. The value passed is selected by + the `o2-sim` option `--seed`. + + - `GeneratorFileOrCmd.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`. + + - `GeneratorFileOrCmd.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`). + + - `GeneratorFileOrCmd.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. + +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_ + + +## TODO + +### 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..ae188024808d0 --- /dev/null +++ b/run/SimExamples/TParticle/child.sh @@ -0,0 +1,57 @@ +#!/usr/bin/env 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 ' ' '_'` + +set -x +export VMCWORKDIR=${O2_ROOT}/share +o2-sim -g tparticle \ + --configKeyValues "GeneratorFileOrCmd.cmd=$cmd $opt;GeneratorFileOrCmd.outputSwitch=-o" \ + --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 +export VMCWORKDIR=${O2_ROOT}/share +o2-sim -g tparticle --configKeyValues "GeneratorFileOrCmd.fileNames=${inp}" \ + --outPrefix "$out" --seed $seed --nEvents $nev $@ +