diff --git a/Generator/src/GtBeamBackgroundTool.cpp b/Generator/src/GtBeamBackgroundTool.cpp index fba5b0b4da7c995d8c09ceb7581cfdbf9de79ba0..c21a12919febfb722d38f6e228e38af714f625fa 100644 --- a/Generator/src/GtBeamBackgroundTool.cpp +++ b/Generator/src/GtBeamBackgroundTool.cpp @@ -4,6 +4,8 @@ #include "BeamBackgroundFileParserV0.h" #include "GuineaPigPairsFileParser.h" +#include "CLHEP/Random/RandPoisson.h" + #include "TVector3.h" // for rotation DECLARE_COMPONENT(GtBeamBackgroundTool) @@ -51,37 +53,47 @@ bool GtBeamBackgroundTool::mutate(MyHepMC::GenEvent& event) { // TODO: should sample according to the rates // dummy: get one and stop to generate for (auto& [label, parser]: m_beaminputs) { - IBeamBackgroundFileParser::BeamBackgroundData beamdata; - auto isok = parser->load(beamdata); - if (not isok) { - error() << "Failed to load beam background data from the parser " << label << endmsg; - return false; - } - // fill the value - float charge = beamdata.pdgid == 11 ? -1: 1; - - TVector3 pos(beamdata.x,beamdata.y,beamdata.z); - TVector3 mom(beamdata.px,beamdata.py,beamdata.pz); + int nparticles = 1; // by default, only one particle - auto itrot = m_rotYmaps.find(label); - if (itrot != m_rotYmaps.end() ) { - info() << "Apply rotation along Y " << itrot->second << endmsg; + if (auto itRate = m_ratemaps.find(label); itRate != m_ratemaps.end()) { + nparticles = itRate->second * m_timewindow; - pos.RotateY(itrot->second); - mom.RotateY(itrot->second); + // now sampling the number of particles + nparticles = CLHEP::RandPoisson::shoot(nparticles); } - // create the MC particle - auto mcp = event.m_mc_vec.create(); - mcp.setPDG(beamdata.pdgid); - mcp.setGeneratorStatus(1); - mcp.setSimulatorStatus(1); - mcp.setCharge(static_cast<float>(charge)); - mcp.setTime(beamdata.t); - mcp.setMass(beamdata.mass); - mcp.setVertex(edm4hep::Vector3d(pos.X(), pos.Y(), pos.Z())); - mcp.setMomentum(edm4hep::Vector3f(mom.X(), mom.Y(), mom.Z())); - + for (auto iparticle = 0; iparticle < nparticles; ++iparticle) { + IBeamBackgroundFileParser::BeamBackgroundData beamdata; + auto isok = parser->load(beamdata); + if (not isok) { + error() << "Failed to load beam background data from the parser " << label << endmsg; + return false; + } + // fill the value + float charge = beamdata.pdgid == 11 ? -1: 1; + + TVector3 pos(beamdata.x,beamdata.y,beamdata.z); + TVector3 mom(beamdata.px,beamdata.py,beamdata.pz); + + auto itrot = m_rotYmaps.find(label); + if (itrot != m_rotYmaps.end() ) { + info() << "Apply rotation along Y " << itrot->second << endmsg; + + pos.RotateY(itrot->second); + mom.RotateY(itrot->second); + } + + // create the MC particle + auto mcp = event.m_mc_vec.create(); + mcp.setPDG(beamdata.pdgid); + mcp.setGeneratorStatus(1); + mcp.setSimulatorStatus(1); + mcp.setCharge(static_cast<float>(charge)); + mcp.setTime(beamdata.t); + mcp.setMass(beamdata.mass); + mcp.setVertex(edm4hep::Vector3d(pos.X(), pos.Y(), pos.Z())); + mcp.setMomentum(edm4hep::Vector3f(mom.X(), mom.Y(), mom.Z())); + } } return true; diff --git a/Generator/src/GtBeamBackgroundTool.h b/Generator/src/GtBeamBackgroundTool.h index 632f3851fea0e8468d53d9e63ccc42ad34b65dff..b3f2254888ea3ad871288ce7ce39a0a17c61ce42 100644 --- a/Generator/src/GtBeamBackgroundTool.h +++ b/Generator/src/GtBeamBackgroundTool.h @@ -54,7 +54,8 @@ private: private: Gaudi::Property<std::map<std::string, std::string>> m_inputmaps{this, "InputFileMap"}; Gaudi::Property<std::map<std::string, std::string>> m_formatmaps{this, "InputFormatMap"}; - Gaudi::Property<std::map<std::string, double>> m_ratemaps {this, "InputRateMap"}; + Gaudi::Property<std::map<std::string, double>> m_ratemaps {this, "InputRateMap"}; // unit: Hz + Gaudi::Property<double> m_timewindow{this, "TimeWindow", 1e-6}; // unit: s // unit of beam energy: GeV Gaudi::Property<std::map<std::string, double>> m_Ebeammaps{this, "InputBeamEnergyMap"};