ATLAS Offline Software
BeamHaloGeneratorAlg.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
9 
10 #include "GaudiKernel/MsgStream.h"
11 #include "AtlasHepMC/GenEvent.h"
13 #include <cmath>
14 
15 //--------------------------------------------------------------------------
16 
17 BeamHaloGeneratorAlg::BeamHaloGeneratorAlg(const std::string& name, ISvcLocator* svcLocator)
18  : GenModule(name,svcLocator)
19 {
20  for (int i=0;i<NPLOTS;i++) {
21  m_validationPlots[i] = nullptr;
22  }
23 }
24 
25 //--------------------------------------------------------------------------
26 
28  ATH_MSG_INFO("Initialising this BeamHaloGeneratorAlg: " << name());
29 
30  ATH_MSG_INFO("================ Generator Settings =================");
31  ATH_MSG_INFO(" input type ------------------ " << m_inputTypeStr);
32  ATH_MSG_INFO(" input file ------------------ " << m_inputFile);
33  ATH_MSG_INFO(" interface plane ------------- " << m_interfacePlane << " (mm)");
34  ATH_MSG_INFO(" enable flip about z=0 ------- " << (m_enableFlip ? "True" : "False"));
35  ATH_MSG_INFO(" flip probability ----------- " << m_flipProbability);
36  ATH_MSG_INFO(" enable sampling ------------- " << (m_enableSampling ? "True" : "False"));
37  ATH_MSG_INFO(" binary buffer file ---------- " << m_bufferFileName);
40  for(;itr!=itr_end;++itr) ATH_MSG_INFO((*itr));
41  ATH_MSG_INFO(" Produce monitoring plots ---- " << (m_doMonitoringPlots ? "True" : "False"));
42  ATH_MSG_INFO(" Random stream name ---------- " << m_randomStream);
43  ATH_MSG_INFO("=====================================================");
44 
45  // Check the flip probability
46  if(m_flipProbability <= 0) {
47  ATH_MSG_INFO("Event flipping will be disabled.");
48  m_enableFlip = false;
49  }
50  else if(m_flipProbability > 1) {
51  ATH_MSG_WARNING("Flip probability " << m_flipProbability << " is out of range. Event flipping will be disabled.");
52  m_enableFlip = false;
53  }
54  else {
55  m_enableFlip = true;
56  }
57 
58  // Retrieve pointer to THistSvc if monitoring plots are requested
61  StatusCode sc = service("THistSvc", m_tHistSvc);
62  if (sc.isFailure() || !m_tHistSvc) {
63  ATH_MSG_FATAL("Unable to locate Service THistSvc");
64  return sc;
65  }
66 
67  // Create log10 bins for x-axis of E plots.
68  const Int_t nbins_E = 60;
69  double xmin = 1e-2;
70  double xmax = 3.5e3;
71  double logxmin = std::log10(xmin);
72  double logxmax = std::log10(xmax);
73  double binwidth = (logxmax-logxmin)/nbins_E;
74  Double_t xbins[nbins_E+1];
75  xbins[0] = xmin;
76  for (Int_t i=1;i<=nbins_E;i++) {
77  xbins[i] = xmin + std::pow(10,logxmin+i*binwidth);
78  }
79 
80  // Create the monitoring plots
81  m_validationPlots[PRI_R] = new TH1F("primaryR",";Radius [m];Events/[m]",100,0., 0.2);
82  m_validationPlots[PRI_Z] = new TH1F("primaryZ",";z [m];Events/[m]",1100,0., 550.); // Beam gas and TCT
83  m_validationPlots[PRI_Z_TCT] = new TH1F("primaryZ_TCT",";z [m];Events/[m]",100,144., 149.); // TCT region
84  m_validationPlots[SP_R_ALL] = new TH1F("scoringPlaneR_all",";Radius [m];Particles/[m]",240,0., 24.);
85  m_validationPlots[SP_E_ALL] = new TH1F("scoringPlaneE_all",";Energy [GeV];Particles/[GeV]",nbins_E,xbins);
86  m_validationPlots[SP_PZ_ALL] = new TH1F("scoringPlanePz_all",";p_{z} [GeV];Particles/[GeV]",350,0., 3500.);
87  m_validationPlots[SP_PT_ALL] = new TH1F("scoringPlanePt_all",";p_{T} [GeV];Particles/[GeV]",500,0., 50.);
88  m_validationPlots[SP_R_PROTONS] = new TH1F("scoringPlaneR_protons",";Radius [m];Particles/[m]",240,0., 24.);
89  m_validationPlots[SP_E_PROTONS] = new TH1F("scoringPlaneE_protons",";Energy [GeV];Particles/[GeV]",nbins_E,xbins);
90  m_validationPlots[SP_PZ_PROTONS] = new TH1F("scoringPlanePz_protons",";p_{z} [GeV];Particles/[GeV]",350,0., 3500.);
91  m_validationPlots[SP_PT_PROTONS] = new TH1F("scoringPlanePt_protons",";p_{T} [GeV];Particles/[GeV]",500,0., 50.);
92  m_validationPlots[SP_R_MUONS] = new TH1F("scoringPlaneR_muons",";Radius [m];Particles/[m]",240,0., 24.);
93  m_validationPlots[SP_E_MUONS] = new TH1F("scoringPlaneE_muons",";Energy [GeV];Particles/[GeV]",nbins_E,xbins);
94  m_validationPlots[SP_PZ_MUONS] = new TH1F("scoringPlanePz_muons",";p_{z} [GeV];Particles/[GeV]",350,0., 3500.);
95  m_validationPlots[SP_PT_MUONS] = new TH1F("scoringPlanePt_muons",";p_{T} [GeV];Particles/[GeV]",500,0., 50.);
96 
97  if((sc = m_tHistSvc->regHist("/BeamHalo/primaryR", m_validationPlots[PRI_R])) == StatusCode::FAILURE) return sc;
98  if((sc = m_tHistSvc->regHist("/BeamHalo/primaryZ", m_validationPlots[PRI_Z])) == StatusCode::FAILURE) return sc;
99  if((sc = m_tHistSvc->regHist("/BeamHalo/primaryZ_TCT", m_validationPlots[PRI_Z_TCT])) == StatusCode::FAILURE) return sc;
100  if((sc = m_tHistSvc->regHist("/BeamHalo/scoringPlaneR_all", m_validationPlots[SP_R_ALL])) == StatusCode::FAILURE) return sc;
101  if((sc = m_tHistSvc->regHist("/BeamHalo/scoringPlaneE_all", m_validationPlots[SP_E_ALL])) == StatusCode::FAILURE) return sc;
102  if((sc = m_tHistSvc->regHist("/BeamHalo/scoringPlanePz_all", m_validationPlots[SP_PZ_ALL])) == StatusCode::FAILURE) return sc;
103  if((sc = m_tHistSvc->regHist("/BeamHalo/scoringPlanePt_all", m_validationPlots[SP_PT_ALL])) == StatusCode::FAILURE) return sc;
104  if((sc = m_tHistSvc->regHist("/BeamHalo/scoringPlaneR_protons", m_validationPlots[SP_R_PROTONS])) == StatusCode::FAILURE) return sc;
105  if((sc = m_tHistSvc->regHist("/BeamHalo/scoringPlaneE_protons", m_validationPlots[SP_E_PROTONS])) == StatusCode::FAILURE) return sc;
106  if((sc = m_tHistSvc->regHist("/BeamHalo/scoringPlanePz_protons", m_validationPlots[SP_PZ_PROTONS])) == StatusCode::FAILURE) return sc;
107  if((sc = m_tHistSvc->regHist("/BeamHalo/scoringPlanePt_protons", m_validationPlots[SP_PT_PROTONS])) == StatusCode::FAILURE) return sc;
108  if((sc = m_tHistSvc->regHist("/BeamHalo/scoringPlaneR_muons", m_validationPlots[SP_R_MUONS])) == StatusCode::FAILURE) return sc;
109  if((sc = m_tHistSvc->regHist("/BeamHalo/scoringPlaneE_muons", m_validationPlots[SP_E_MUONS])) == StatusCode::FAILURE) return sc;
110  if((sc = m_tHistSvc->regHist("/BeamHalo/scoringPlanePz_muons", m_validationPlots[SP_PZ_MUONS])) == StatusCode::FAILURE) return sc;
111  if((sc = m_tHistSvc->regHist("/BeamHalo/scoringPlanePt_muons", m_validationPlots[SP_PT_MUONS])) == StatusCode::FAILURE) return sc;
112  }
113 
114  // Check the input type string
115  if (m_inputTypeStr == "MARS-NM") {
117  m_inputFile,
119  }
120  else if (m_inputTypeStr == "FLUKA-VT") {
122  m_inputFile,
124  }
125  else if (m_inputTypeStr == "FLUKA-RB") {
127  m_inputFile,
129  }
130  else {
132  ATH_MSG_FATAL("Input type " << m_inputTypeStr << " is not known. Available types are: MARS-NM or FLUKA-VT");
133  return StatusCode::FAILURE;
134  }
135 
136  // Set the options using those specified to this algorithm
143 
144  // Initialise the generator.
145  if(m_beamHaloGenerator->genInitialize() != 0) return StatusCode::FAILURE;
146 
147  return StatusCode::SUCCESS;
148 }
149 
150 //---------------------------------------------------------------------------
151 
153 
154  const EventContext& ctx = Gaudi::Hive::currentContext();
155  CLHEP::HepRandomEngine* rndmEngine = this->getRandomEngine(m_randomStream, ctx);
156 
157 
158  // Clear the event ready for it to be filled with the next beam halo event.
159  m_evt.clear();
160 
161  // Fill an event with particles from the converted ASCII input file.
162  // (If the end of file has already been reached this will return a
163  // non-zero value.)
164  if(m_beamHaloGenerator->fillEvt(&m_evt, rndmEngine) != 0) return StatusCode::FAILURE;
165 
166  // Fill monitoring plots if requested
167  if(m_doMonitoringPlots) {
168  auto weightContainer = m_evt.weights();
169  if(weightContainer.size() != 5) {
170  ATH_MSG_WARNING("The number of weights for this event is not equal to 5.");
171  return StatusCode::SUCCESS;
172  }
173  double weight = weightContainer[0];
174  HepMC::FourVector primaryPosition(weightContainer[1],weightContainer[2],weightContainer[3],weightContainer[4]);
175  m_validationPlots[PRI_R]->Fill(primaryPosition.perp()/1000.0,weight);
176  m_validationPlots[PRI_Z]->Fill(weightContainer[3]/1000.0,weight);
177  m_validationPlots[PRI_Z_TCT]->Fill(std::fabs(weightContainer[3]/1000.0),weight);
178 
179  double values[4];
180  int pdgId;
181  for (auto hepmc_part: m_evt) {
182  auto prodVertex = hepmc_part->production_vertex();
183  if(!prodVertex) continue;
184 
185  // Store the values for use in the if conditions that follow
186  values[0] = prodVertex->position().perp()/1000.;
187  values[1] = hepmc_part->momentum().e()/1000.;
188  values[2] = hepmc_part->momentum().pz()/1000.;
189  values[3] = hepmc_part->momentum().perp()/1000.;
190 
191  pdgId = hepmc_part->pdg_id(); if(pdgId<0) pdgId = -pdgId;
196  if(pdgId == 2212) {
201  }
202  else if(pdgId == 13) {
207  }
208  }
209  }
210 
211  return StatusCode::SUCCESS;
212 }
213 
214 //---------------------------------------------------------------------------
215 
217  *event = m_evt;
218  return StatusCode::SUCCESS;
219 }
220 
221 //---------------------------------------------------------------------------
222 
224  if(m_beamHaloGenerator->genFinalize() != 0) return StatusCode::FAILURE;
225  return StatusCode::SUCCESS;
226 }
227 
228 //---------------------------------------------------------------------------
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
BeamHaloGeneratorAlg::m_inputFile
StringProperty m_inputFile
Input file name.
Definition: BeamHaloGeneratorAlg.h:50
FlukaHaloGenerator.h
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
BeamHaloGenerator::setEnableSampling
void setEnableSampling(bool enableSampling)
Turn on or off sampling from the input ASCII file.
Definition: BeamHaloGenerator.h:72
GenEvent.h
BeamHaloGeneratorAlg::m_flipProbability
FloatProperty m_flipProbability
Flip probability.
Definition: BeamHaloGeneratorAlg.h:59
BeamHaloGeneratorAlg::m_enableFlip
BooleanProperty m_enableFlip
Flag for flipping event.
Definition: BeamHaloGeneratorAlg.h:56
BeamHaloGenerator.h
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
BeamHaloGeneratorAlg::m_enableSampling
BooleanProperty m_enableSampling
Flag to enable or disable sampling.
Definition: BeamHaloGeneratorAlg.h:62
BeamHaloGenerator::fillEvt
virtual int fillEvt(HepMC::GenEvent *evt, CLHEP::HepRandomEngine *engine)=0
A function to create one event in HepMC format.
BeamHaloGeneratorAlg::genInitialize
virtual StatusCode genInitialize()
For initializing the generator, if required.
Definition: BeamHaloGeneratorAlg.cxx:27
binwidth
bool binwidth
Definition: listroot.cxx:58
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
BeamHaloGeneratorAlg::SP_PZ_MUONS
@ SP_PZ_MUONS
Definition: BeamHaloGeneratorAlg.h:101
BeamHaloGeneratorAlg::fillEvt
virtual StatusCode fillEvt(HepMC::GenEvent *evt)
Fill the GenEvent pointer with the contents of the GenEvent cache.
Definition: BeamHaloGeneratorAlg.cxx:216
BeamHaloGeneratorAlg::m_bufferFileName
StringProperty m_bufferFileName
The name of the binary buffer file, needed for sampling from a converted file.
Definition: BeamHaloGeneratorAlg.h:66
AthCommonMsg< Algorithm >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition: AthCommonMsg.h:30
BeamHaloGeneratorAlg::m_generatorSettings
StringArrayProperty m_generatorSettings
A vector of strings defining generator settings.
Definition: BeamHaloGeneratorAlg.h:69
BeamHaloGeneratorAlg::m_tHistSvc
ITHistSvc * m_tHistSvc
A pointer to the THist service for validation plots.
Definition: BeamHaloGeneratorAlg.h:75
BeamHaloGeneratorAlg::SP_PZ_ALL
@ SP_PZ_ALL
Definition: BeamHaloGeneratorAlg.h:93
MarsHaloGenerator
A class to provide conversion from a MARS format ASCII input record into HepMC format,...
Definition: MarsHaloGenerator.h:18
MarsHaloGenerator.h
BeamHaloGeneratorAlg::SP_R_MUONS
@ SP_R_MUONS
Definition: BeamHaloGeneratorAlg.h:99
BeamHaloGeneratorAlg::PRI_Z
@ PRI_Z
Definition: BeamHaloGeneratorAlg.h:89
PowhegPy8EG_H2a.pdgId
dictionary pdgId
Definition: PowhegPy8EG_H2a.py:128
BeamHaloGeneratorAlg::m_interfacePlane
DoubleProperty m_interfacePlane
The position of the interface plane in mm.
Definition: BeamHaloGeneratorAlg.h:53
GenModule::getRandomEngine
CLHEP::HepRandomEngine * getRandomEngine(const std::string &streamName, const EventContext &ctx) const
Definition: GenModule.cxx:34
BeamHaloGenerator::setDebugEnable
void setDebugEnable(bool debug)
A function to turn on or off debug print out.
Definition: BeamHaloGenerator.h:79
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
python.Bindings.values
values
Definition: Control/AthenaPython/python/Bindings.py:797
BeamHaloGeneratorAlg::SP_R_ALL
@ SP_R_ALL
Definition: BeamHaloGeneratorAlg.h:91
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:200
BeamHaloGeneratorAlg.h
GenModule
Base class for common behaviour of generator interfaces.
Definition: GenModule.h:39
BeamHaloGeneratorAlg::m_beamHaloGenerator
BeamHaloGenerator * m_beamHaloGenerator
A pointer to the beam halo generator.
Definition: BeamHaloGeneratorAlg.h:81
BeamHaloGeneratorAlg::SP_PT_ALL
@ SP_PT_ALL
Definition: BeamHaloGeneratorAlg.h:94
FlukaHaloGenerator
A class to provide conversion from a FLUKA format ASCII input record into HepMC format,...
Definition: FlukaHaloGenerator.h:18
BeamHaloGenerator::genFinalize
virtual int genFinalize()
A function to finalise the generator.
Definition: BeamHaloGenerator.cxx:57
BeamHaloGeneratorAlg::m_inputTypeStr
StringProperty m_inputTypeStr
Input file type and therefore associated beam halo generator that should be used.
Definition: BeamHaloGeneratorAlg.h:47
event
POOL::TEvent event(POOL::TEvent::kClassAccess)
lumiFormat.i
int i
Definition: lumiFormat.py:92
xmin
double xmin
Definition: listroot.cxx:60
BeamHaloGeneratorAlg::genFinalize
virtual StatusCode genFinalize()
For finalising the generator, if required.
Definition: BeamHaloGeneratorAlg.cxx:223
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
BeamHaloGeneratorAlg::NPLOTS
@ NPLOTS
Definition: BeamHaloGeneratorAlg.h:103
BeamHaloGeneratorAlg::m_validationPlots
TH1F * m_validationPlots[NPLOTS]
An array of TH1F pointers for validation plots.
Definition: BeamHaloGeneratorAlg.h:106
BeamHaloGeneratorAlg::m_randomStream
StringProperty m_randomStream
Name of the random number stream.
Definition: BeamHaloGeneratorAlg.h:78
BeamHaloGeneratorAlg::callGenerator
virtual StatusCode callGenerator()
Read one event from the selected input and convert it into GenEvent format.
Definition: BeamHaloGeneratorAlg.cxx:152
BeamHaloGeneratorAlg::PRI_Z_TCT
@ PRI_Z_TCT
Definition: BeamHaloGeneratorAlg.h:90
BeamHaloGeneratorAlg::SP_R_PROTONS
@ SP_R_PROTONS
Definition: BeamHaloGeneratorAlg.h:95
BeamHaloGeneratorAlg::PRI_R
@ PRI_R
Definition: BeamHaloGeneratorAlg.h:88
BeamHaloGeneratorAlg::SP_E_MUONS
@ SP_E_MUONS
Definition: BeamHaloGeneratorAlg.h:100
GenBase::particleTable
const HepPDT::ParticleDataTable & particleTable() const
Get a particle data table.
Definition: GenBase.h:118
BeamHaloGeneratorAlg::SP_E_ALL
@ SP_E_ALL
Definition: BeamHaloGeneratorAlg.h:92
BeamHaloGeneratorAlg::SP_E_PROTONS
@ SP_E_PROTONS
Definition: BeamHaloGeneratorAlg.h:96
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:195
BeamHaloGenerator::setFlipProbability
void setFlipProbability(float flipProbability)
Set probability of flipping an event about Z=0.
Definition: BeamHaloGenerator.h:69
BeamHaloGenerator::setBufferFileName
void setBufferFileName(const std::string &bufferFileName)
Set the name of the binary buffer file, needed for sampling from a converted file.
Definition: BeamHaloGenerator.h:76
BeamHaloGeneratorAlg::SP_PZ_PROTONS
@ SP_PZ_PROTONS
Definition: BeamHaloGeneratorAlg.h:97
LArCellBinning.xbins
int xbins
Definition: LArCellBinning.py:163
BeamHaloGenerator::setEnableFlip
void setEnableFlip(bool enableFlip)
Turn on or off the event flipping code.
Definition: BeamHaloGenerator.h:66
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
WeightContainer.h
DEBUG
#define DEBUG
Definition: page_access.h:11
BeamHaloGenerator::genInitialize
virtual int genInitialize()
A function to initialise the generator.
Definition: BeamHaloGenerator.cxx:48
BeamHaloGeneratorAlg::SP_PT_PROTONS
@ SP_PT_PROTONS
Definition: BeamHaloGeneratorAlg.h:98
BeamHaloGeneratorAlg::m_evt
HepMC::GenEvent m_evt
An empty GenEvent to cache the generate output between callGenerator and fillEvt.
Definition: BeamHaloGeneratorAlg.h:85
xmax
double xmax
Definition: listroot.cxx:61
BeamHaloGeneratorAlg::SP_PT_MUONS
@ SP_PT_MUONS
Definition: BeamHaloGeneratorAlg.h:102
BeamHaloGenerator::setInterfacePlane
void setInterfacePlane(double interfacePlane)
Set the position of the interface plane in mm.
Definition: BeamHaloGenerator.h:63
BeamHaloGeneratorAlg::m_doMonitoringPlots
BooleanProperty m_doMonitoringPlots
A flag to allow monitoring plots to be turned on or off.
Definition: BeamHaloGeneratorAlg.h:72
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
BeamHaloGeneratorAlg::BeamHaloGeneratorAlg
BeamHaloGeneratorAlg(const std::string &name, ISvcLocator *svcLocator)
Definition: BeamHaloGeneratorAlg.cxx:17