Belle II Software prerelease-10-00-00a
PhokharaInputModule.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <generators/modules/phokharainput/PhokharaInputModule.h>
10
11#include <framework/datastore/StoreArray.h>
12#include <framework/utilities/FileSystem.h>
13#include <framework/utilities/IOIntercept.h>
14
15using namespace std;
16using namespace Belle2;
17
18//-----------------------------------------------------------------
19// Register the Module
20//-----------------------------------------------------------------
21REG_MODULE(PhokharaInput);
22
23//-----------------------------------------------------------------
24// Implementation
25//-----------------------------------------------------------------
27{
28 //Set module properties
29 setDescription("Generates radiative return events with PHOKHARA 10.");
30
31 //Parameter definition
32 addParam("FinalState", m_finalState,
33 "Final state: mu+mu-(0), pi+pi-(1), 2pi0pi+pi-(2), 2pi+2pi-(3), ppbar(4), nnbar(5), K+K-(6), K0K0bar(7), pi+pi-pi0(8), lamb(->pi-p)lambbar(->pi+pbar)(9), eta pi+ pi- (10), chi_c1 production (11), chi_c2 production (12), pi0 g (13), eta g (14), etaP g (15)",
34 1);
35 addParam("ReplaceMuonsByVirtualPhoton", m_replaceMuonsByVirtualPhoton,
36 "Replace muons by a virtual photon (for FinalState == 0 only).",
37 false);
38 addParam("SearchMax", m_nSearchMax, "Number of events used to search for maximum of differential cross section", 100000);
39 addParam("Epsilon", m_epsilon, "Soft/hard photon separator", 0.0001);
40 addParam("nMaxTrials", m_nMaxTrials, "Maximum trials per event", 10000);
41 addParam("Weighted", m_weighted, "generate weighted events", 0);
42 addParam("LO", m_LO, "ph0 Born: 1ph(0), Born: 0ph(1), only Born: 0ph(-1)", 0);
43 addParam("NLO", m_NLO, "1 photon : Born(0), NLO(1)", 0);
44 addParam("FullNLO", m_fullNLO, "full NLO : No(0), Yes(1). Allowed only for ph0=1, nlo=1, fsr=2, fsrnlo=1", 0);
45 addParam("QED", m_QED, "ISR only(0, default), ISR+FSR(1), ISR+INT+FSR(2)", 0);
46 addParam("IFSNLO", m_IFSNLO, "IFSNLO: no(0), yes(1)", 0);
47 addParam("Alpha", m_alpha,
48 "Vacuum polarization switch: off (0), on (1, [by Fred Jegerlehner], default), on (2, [by Thomas Teubner])",
49 1);
50 addParam("PionFF", m_pionff, "Pion FF: KS PionFormFactor(0, default), GS old 1), GS new(2)", 0);
51 addParam("PionStructure", m_pionstructure,
52 "For pi+pi- only: f0+f0(600): K+K- model(0, default), no structure model(1), no f0+f0(600)(2), f0 KLOE(3)", 0);
53 addParam("KaonFF", m_kaonff,
54 "Kaon FF: KaonFormFactor constrained(0, default), KaonFormFactor unconstrained(1) KaonFormFactor old(2)", 0);
55 addParam("NarrowRes", m_narres, "Only for m_finalState = 0,1,6,7: No narrow resonances (0, default), J/Psi (1) and Psi(2S) (2)", 0);
56 addParam("ProtonFF", m_protonff, "ProtonFormFactor old(0), ProtonFormFactor new(1)", 1);
57 addParam("ChiSW", m_chi_sw,
58 "chi_sw: Radiative return(0), Chi production(1), Radiative return + Chi production (2). Works only for pion=11 and pion=12.", 0);
59 addParam("SwitchBeamResolution", m_be_r,
60 "be_r: without beam resolution(0), with beam resolution(1). Works only for pion=11 and pion=12; (0) assumes exact CMS-Energy, (1) each beam ennergy (=CMS-Energy/2) is smeared with Gaussian distribution of the given variance=BeamResolution**2.",
61 0);
62 addParam("BeamResolution", m_beamres, "beamres - beam resolution for pion==11 and pion==12 only", 0.);
63
64 addParam("ScatteringAngleRangePhoton", m_ScatteringAngleRangePhoton,
65 "Min [0] and Max [1] value for the scattering angle of photons [deg], default (0, 180)", make_vector(0.0, 180.0));
66 addParam("ScatteringAngleRangeFinalStates", m_ScatteringAngleRangeFinalStates,
67 "Min [0] and Max [1] value for the scattering angle of pions(muons,nucleons,kaons) [deg], default (0, 180)", make_vector(0.0,
68 180.0));
69
70 addParam("MinInvMassHadronsGamma", m_MinInvMassHadronsGamma, "Minimal hadrons/muons-gamma invariant mass squared [GeV^2]", 0.0);
71 addParam("MinInvMassHadrons", m_MinInvMassHadrons, "Minimal hadrons/muons invariant mass squared [GeV^2]", 0.01);
72 addParam("ForceMinInvMassHadronsCut", m_ForceMinInvMassHadronsCut,
73 "Force application of the MinInvMassHadrons cut. It is ignored by PHOKHARA with LO = 1, NLO = 1.",
74 false);
75 addParam("MaxInvMassHadrons", m_MaxInvMassHadrons, "Maximal hadrons/muons invariant mass squared [GeV^2]", 112.0);
76 addParam("MinEnergyGamma", m_MinEnergyGamma, "Minimal photon energy/missing energy, must be greater than 0.0098 * CMS energy [GeV]",
77 0.15);
78 const std::string defaultParameterFile =
79 FileSystem::findFile("/data/generators/phokhara/const_and_model_paramall10.0.dat");
80 addParam("ParameterFile", m_ParameterFile,
81 "File that contains particle properties.",
82 defaultParameterFile);
83 addParam("BeamEnergySpread", m_BeamEnergySpread,
84 "Simulate beam-energy spread (initializes PHOKHARA for every event - very slow).", false);
85
86}
87
88//-----------------------------------------------------------------
89// Destructor
90//-----------------------------------------------------------------
95
96//-----------------------------------------------------------------
97// Initialize
98//-----------------------------------------------------------------
100{
101 StoreArray<MCParticle> mcparticle;
102 mcparticle.registerInDataStore();
103
105 if (m_finalState != 0) {
106 B2FATAL("You requested to replace muons by a virtual photon, but the "
107 "final state is not mu+ mu-.");
108 }
109 if (m_QED != 0) {
110 B2FATAL("You requested to replace muons by a virtual photon. In this "
111 "mode, PHOKHARA works as an ISR generator. The parameter QED "
112 "should be set to 0 (ISR only). If FSR is taken into account "
113 "(QED = 1 or 2), the results will be incorrect.");
114 }
115 if (m_IFSNLO != 0) {
116 B2FATAL("You requested to replace muons by a virtual photon. In this "
117 "mode, PHOKHARA works as an ISR generator. The parameter IFSNLO "
118 "should be set to 0 (off). If simultaneous emission of initial "
119 "and final-state photons is taken into account (IFSNLO = 1), "
120 "the results will be incorrect.");
121 }
122 }
123 if (m_finalState == 11 || m_finalState == 12) {
124 if (m_narres != 1 || m_LO != 0 || m_NLO != 0) {
125 B2FATAL("Generation of the chi_c1 or chi_c2 (final states 11 and 12, "
126 "respectively) requires to turn on the J/psi (NarrowRes = 1) "
127 "and use LO radiative-return mode (LO = 0, NLO = 0)");
128 }
129 }
130 //Beam Parameters, initial particle - PHOKHARA cannot handle beam energy spread
131 if (m_BeamEnergySpread) {
134 }
135 m_initial.initialize();
136
137}
138
139//-----------------------------------------------------------------
140// Event
141//-----------------------------------------------------------------
143{
144
145 // Check if the BeamParameters have changed (if they do, abort the job! otherwise cross section calculation will be a nightmare.)
146 if (m_beamParams.hasChanged()) {
147 if (!m_initialized) {
149 } else {
150 B2FATAL("PhokharaInputModule::event(): BeamParameters have changed within a job, this is not supported for PHOKHARA!");
151 }
152 }
153
154 StoreObjPtr<EventMetaData> evtMetaData("EventMetaData", DataStore::c_Event);
155
156 // initial particle from beam parameters
157 const MCInitialParticles& initial = m_initial.generate();
158
159 // true boost
160 ROOT::Math::LorentzRotation boost = initial.getCMSToLab();
161
162 // vertex
163 ROOT::Math::XYZVector vertex = initial.getVertex();
164
165 m_mcGraph.clear();
166 if (m_BeamEnergySpread) {
167 // PHOKHARA does not support beam-energy spread. To simulate it,
168 // the generator is initialized with the new energy for every event.
170 "PHOKHARA", LogConfig::c_Debug, LogConfig::c_Debug, 100, 100);
171 initLogCapture.start();
172 m_generator.setCMSEnergy(initial.getMass());
174 initLogCapture.finish();
175 }
176 double weight = m_generator.generateEvent(m_mcGraph, vertex, boost);
178
179 //store the weight (1.0 for unweighted events)
180 evtMetaData->setGeneratedWeight(weight);
181}
182
183
184//-----------------------------------------------------------------
185// Terminate
186//-----------------------------------------------------------------
188{
189
190 m_generator.term();
191}
192
194{
195
196 const BeamParameters& nominal = m_initial.getBeamParameters();
197 m_cmsEnergy = nominal.getMass();
198
199 m_generator.setCMSEnergy(m_cmsEnergy);
200
201 m_generator.setNSearchMax(m_nSearchMax);
202 m_generator.setWeighted(m_weighted);
203 m_generator.setFinalState(m_finalState);
204 m_generator.setReplaceMuonsByVirtualPhoton(m_replaceMuonsByVirtualPhoton);
205 m_generator.setNMaxTrials(m_nMaxTrials);
206 m_generator.setEpsilon(m_epsilon);
207 m_generator.setLO(m_LO);
208 m_generator.setNLO(m_NLO);
209 m_generator.setFullNLO(m_fullNLO);
210 m_generator.setQED(m_QED);
211 m_generator.setIFSNLO(m_IFSNLO);
212 m_generator.setAlpha(m_alpha);
213 m_generator.setPionFF(m_pionff);
214 m_generator.setKaonFF(m_kaonff);
215 m_generator.setPionStructure(m_pionstructure);
216 m_generator.setNarrowRes(m_narres);
217 m_generator.setProtonFF(m_protonff);
218 m_generator.setChiSW(m_chi_sw);
219 m_generator.setSwitchBeamResolution(m_be_r);
220 m_generator.setBeamResolution(m_beamres);
221
222
223 m_generator.setMinInvMassHadronsGamma(m_MinInvMassHadronsGamma);
224 m_generator.setm_MinInvMassHadrons(m_MinInvMassHadrons);
225 m_generator.setForceMinInvMassHadronsCut(m_ForceMinInvMassHadronsCut);
226 m_generator.setm_MaxInvMassHadrons(m_MaxInvMassHadrons);
227 m_generator.setMinEnergyGamma(m_MinEnergyGamma);
228 m_generator.setScatteringAngleRangePhoton(vectorToPair(m_ScatteringAngleRangePhoton, "ScatteringAngleRangePhoton"));
229 m_generator.setScatteringAngleRangeFinalStates(vectorToPair(m_ScatteringAngleRangeFinalStates, "ScatteringAngleRangeFinalStates"));
230
232
233 m_initialized = true;
234
235}
This class contains the nominal beam parameters and the parameters used for smearing of the primary v...
@ c_Event
Different object in each event, all objects/arrays are invalidated after event() function has been ca...
Definition DataStore.h:59
static std::string findFile(const std::string &path, bool silent=false)
Search for given file or directory in local or central release directory, and return absolute path if...
bool start()
Start intercepting the output.
Capture stdout and stderr and convert into log messages.
bool finish()
Finish the capture and emit the message if output has appeard on stdout or stderr.
@ c_Debug
Debug: for code development.
Definition LogConfig.h:26
This class contains the initial state for the given event.
@ c_smearBeam
smear the full beam momentum (energy and direction)
const ROOT::Math::LorentzRotation & getCMSToLab() const
Return the LorentzRotation to convert from CMS to lab frame.
const ROOT::Math::XYZVector & getVertex() const
Get the position of the collision.
double getMass() const
Get the invariant mass of the collision (= energy in CMS)
@ c_checkCyclic
Check for cyclic dependencies.
@ c_setDecayInfo
Set decay time and vertex.
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
double m_beamres
beam resolution for chi2 studies
int m_pionstructure
for pi+pi- only: f0+f0(600): K+K- model(0), "no structure" model(1), no f0+f0(600)(2),...
int m_LO
LO: ph0 Born: 1ph(0), Born: 0ph(1), only Born: 0ph(-1)
std::vector< double > m_ScatteringAngleRangePhoton
Minimal/Maximal photon angle/missing momentum angle.
bool m_ForceMinInvMassHadronsCut
Force application of the above cut.
int m_finalState
Module parameters.
int m_QED
ISR only(0), ISR+FSR(1), ISR+INT+FSR(2)
DBObjPtr< BeamParameters > m_beamParams
BeamParameter.
bool m_initialized
True if generator has been initialized.
int m_narres
narr_res: no narrow resonances (0), J/psi (1) and psi(2S) (2) only for m_finalState = 0,...
int m_alpha
vacuum polarization switch: off (0), on (1,[by Fred Jegerlehner]), on (2,[by Thomas Teubner])
int m_chi_sw
chi_sw: Radiative return(0), Chi production(1), Radiative return + Chi production (2)
double m_epsilon
Soft/hard photon separator.
int m_pionff
FF_pion: KS PionFormFactor(0),GS old (1),GS new (2)
virtual ~PhokharaInputModule()
Destructor.
virtual void terminate() override
Method is called at the end of the event processing.
int m_protonff
ProtonFormFactor old(0), ProtonFormFactor new(1)
int m_IFSNLO
IFSNLO: no(0), yes(1)
int m_nSearchMax
Events used to search maximum of differential cross section.
int m_nMaxTrials
Events before loop is aborted.
int m_be_r
be_r: without beam resolution(0), with beam resolution(1)
virtual void generatorInitialize() override
Initializes the module.
double m_MinInvMassHadronsGamma
m_MinInvMassHadronsGamma [GeV^2].
std::string m_ParameterFile
file that holds all resonance parameters
MCParticleGraph m_mcGraph
The MCParticle graph object.
std::vector< double > m_ScatteringAngleRangeFinalStates
Minimal/Maximal pions(muons,nucleons,kaons) momentum angle.
int m_kaonff
FF_kaon: KaonFormFactor constrained(0), KaonFormFactor unconstrained(1) KaonFormFactor old(2)
InitialParticleGeneration m_initial
initial particle used by BeamParameter class
int m_weighted
generate weighted events
int m_NLO
1 photon : Born(0), NLO(1)
double m_MaxInvMassHadrons
m_MaxInvMassHadrons [GeV^2]
int m_fullNLO
full NLO : No(0), Yes(1)
double m_MinEnergyGamma
m_MinEnergyGamma [GeV].
double m_MinInvMassHadrons
m_MinInvMassHadrons [GeV^2].
void initializeGenerator()
Method is called to initialize the generator.
bool m_BeamEnergySpread
Simulate beam-energy spread.
virtual void generatorEvent() override
Method is called for each event.
bool m_replaceMuonsByVirtualPhoton
Replace muons by a virtual photon.
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
Accessor to arrays stored in the data store.
Definition StoreArray.h:113
Type-safe access to single objects in the data store.
Definition StoreObjPtr.h:96
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
std::pair< T, T > vectorToPair(std::vector< T > &vec, const std::string &name="")
std::vector< T > make_vector(T const &t1, T const &t2)
make_vector.
Abstract base class for different kinds of events.
STL namespace.