Belle II Software development
BabayagaNLOInputModule.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/babayaganloinput/BabayagaNLOInputModule.h>
10
11#include <framework/datastore/StoreArray.h>
12#include <framework/utilities/FileSystem.h>
13#include <framework/utilities/IOIntercept.h>
14
15#include <Math/Vector3D.h>
16
17using namespace std;
18using namespace Belle2;
19
20//-----------------------------------------------------------------
21// Register the Module
22//-----------------------------------------------------------------
23REG_MODULE(BabayagaNLOInput);
24
25//-----------------------------------------------------------------
26// Implementation
27//-----------------------------------------------------------------
28
29extern "C" {
30 double babayaganlo_getrandomcmsenergy_()
31 {
33 }
34}
35
36// BabayagaNLOInputModule::BabayagaNLOInputModule() : Module(), s_initial(BeamParameters::c_smearALL)
38
40{
41 //Set module properties
42 setDescription("Generates radiative Bhabha scattering and exclusive two-photon events with the high precision QED generator called BabaYaga@NLO.");
43
44 //Parameter definition
45 addParam("VacuumPolarization", m_vacPol, "Vacuum polarization: off, hadr5 (Jegerlehner) or hlmnt (Teubner, default)",
46 std::string("hlmnt"));
47 addParam("Model", m_model, "Model: ps or matched (default)", std::string("matched"));
48 addParam("Order", m_order, "Order: born, alpha, exp (default)", std::string("exp"));
49 addParam("Mode", m_mode, "Mode: weighted or unweighted (default)", std::string("unweighted"));
50 addParam("FinalState", m_finalState, "Final state: ee (default), mm or gg", std::string("ee"));
51 addParam("MinEnergyFrac", m_eMinFrac,
52 "Fractional minimum energy for leptons (ee or mm mode) or photons (gg mode) in the final state [fraction of ECMS]", -1.0);
53 addParam("MinEnergy", m_eMin, "Minimum energy for leptons (ee or mm mode) or photons (gg mode) in the final state [GeV]", 0.10);
54 addParam("Epsilon", m_epsilon, "Soft/hard photon separator [fraction of ECMS/2], must be <=1e-7 for ee and mumu, <=1e-5 for gg",
55 1.e-7);
56 addParam("MaxAcollinearity", m_maxAcollinearity, "Maximum acollinearity angle between finale state leptons/photons [degree]",
57 180.0);
58 addParam("FMax", m_fMax, "Maximum of differential cross section weight (fmax)", -1.);
59 addParam("NPhotons", m_nPhot, "Fixed number of (hard) photons are generated, -1 for any number", -1);
60 addParam("SearchMax", m_nSearchMax, "Number of events used to search for maximum of differential cross section", 500000);
61 addParam("ScatteringAngleRange", m_ScatteringAngleRange, "Min [0] and Max [1] value for the scattering angle [deg].",
62 make_vector(15.0, 165.0));
63 addParam("ExtraFile", m_fileNameExtraInfo, "ROOT file that contains additional information.", std::string(""));
64 addParam("DebugEnergySpread", m_Spread, "TEMPORARY SOLUTION! Approximate energy spread per beam (CMS)", 5.e-3);
65 addParam("VPUncertainty", m_Uncertainty, "Calculate VP uncertainty by internal reweighting", false);
66 addParam("NSKDataFile", m_NSKDataFile, "File that contain VP data from Novosibirsk (nsk)",
67 FileSystem::findFile("/data/generators/babayaganlo/vpol_novosibirsk.dat"));
68
69 // user cuts
70 addParam("UserMode", m_userMode, "User mode similar to TEEGG: ETRON, EGAMMA, GAMMA or PRESCALE or NONE (default)",
71 std::string("NONE"));
72 addParam("EEMIN", m_eemin, "Minimum CMS energy of the tagged e-/e+ (GeV)", -1.0);
73 addParam("TEMIN", m_temin, "Minimum CMS angle between the tagged e-/e+ and -z axis (deg)", -1.0);
74 addParam("EGMIN", m_egmin, "Minimum CMS energy of the gamma (GeV)", -1.0);
75 addParam("TGMIN", m_tgmin, "Minimum CMS angle between the gamma and -z axis (deg)", -1.0);
76 addParam("EEVETO", m_eeveto, "Minimum CMS energy to veto e-/e+ (GeV)", -1.0);
77 addParam("TEVETO", m_teveto, "Maximum CMS theta of e-/e+ in final state (deg)", -1.0);
78 addParam("EGVETO", m_egveto, "Minimum CMS energy to veto gamma (GeV)", -1.0);
79 addParam("TGVETO", m_tgveto, "Maximum CMS angle between the gamma and -z axis (deg)", -1.0);
80 addParam("MaxPrescale", m_maxprescale, "Maximum prescale factor (used for maximum differential cross section)", 1.0);
81
82 //initialize member variables
84 m_th1dSDif = 0;
85
86}
87
91
93{
94 //Initialize MCParticle collection
95 StoreArray<MCParticle> mcparticle;
96 mcparticle.registerInDataStore();
97
98 //open extrafile
99 if (m_fileNameExtraInfo != "") {
100 m_fileExtraInfo = new TFile(m_fileNameExtraInfo.c_str(), "RECREATE") ;
101 m_fileExtraInfo->cd();
102 m_th1dSDif = new TH1D("sdif", "sdif", 1000, 0., 100000.);
103 }
104
105 //Initialize initial particle for beam parameters.
106 s_initial.initialize();
107
108 // Initialize ExtraInfo (hold prescale values)
109 m_generator.initExtraInfo();
110}
111
113{
114 // Check if the BeamParameters have changed (if they do, abort the job! otherwise cross section calculation will be a nightmare.)
115 if (m_beamParams.hasChanged()) {
116 if (!m_initialized) {
118 } else {
119 B2FATAL("BabayagaNLOInputModule::event(): BeamParameters have changed within a job, this is not supported for BABAYAGA!");
120 }
121 }
122
123 // initial particle from beam parameters
124 const MCInitialParticles& initial = s_initial.generate();
125
126 // CM energy
127 double ecm = initial.getMass();
128
129 // true boost (per event!)
130 ROOT::Math::LorentzRotation boost = initial.getCMSToLab();
131
132 // vertex
133 ROOT::Math::XYZVector vertex = initial.getVertex();
134
135 m_mcGraph.clear();
136 m_generator.generateEvent(m_mcGraph, ecm, vertex, boost); // actual generator call
137
139}
140
145
147{
148 // generator parameters
149 m_generator.setScatAngle(vectorToPair(m_ScatteringAngleRange, "ScatteringAngleRange"));
150 m_generator.setMaxAcollinearity(m_maxAcollinearity);
151 m_generator.setFinalState(m_finalState);
152 m_generator.setEpsilon(m_epsilon);
153 m_generator.setVacPol(m_vacPol);
154 m_generator.setOrder(m_order);
155 m_generator.setModel(m_model);
156 m_generator.setMode(m_mode);
157 m_generator.setEnergySpread(m_Spread);
158 m_generator.setVPUncertainty(m_Uncertainty);
159 m_generator.setNPhotons(m_nPhot);
160 m_generator.setNSKDataFile(m_NSKDataFile);
161 m_generator.setEEMIN(m_eemin);
162 m_generator.setTEMIN(m_temin);
163 m_generator.setEGMIN(m_egmin);
164 m_generator.setTGMIN(m_tgmin);
165 m_generator.setEEVETO(m_eeveto);
166 m_generator.setTEVETO(m_teveto);
167 m_generator.setEGVETO(m_egveto);
168 m_generator.setTGVETO(m_tgveto);
169 m_generator.setUserMode(m_userMode);
170 m_generator.setMaxPrescale(m_maxprescale);
171
172 // set a nominal value for some initialization inside the generator
173 const BeamParameters& nominal = s_initial.getBeamParameters();
174 double ecmnominal = nominal.getMass();
175 m_generator.setCmsEnergyNominal(ecmnominal);
176
177 // check which of the emin values are used:
178 // its either a fixed value or a fractional one
179 // if both are provided, the fractional has priority
180 if (m_eMinFrac >= 0.) {
181 B2INFO("Setting EMIN using fraction " << m_eMinFrac << " of CMS energy: cut=" << m_eMinFrac * ecmnominal << " GeV");
182 m_generator.setMinEnergy(m_eMinFrac * ecmnominal);
183 } else m_generator.setMinEnergy(m_eMin);
184
185 // check if a maximum weight is provided
186 if (m_fMax > 0. && (m_mode == "unweighted" || m_mode == "uw")) {
187 B2INFO("Setting FMax manually (no maximum search is performed)." <<
188 LogVar("FMax", m_fMax));
189 m_generator.setFMax(m_fMax);
190 m_generator.setNSearchMax(-1);
191 } else if (m_mode == "weighted" || m_mode == "w") {
192 m_generator.setNSearchMax(-1);
193 } else {
194 m_generator.setNSearchMax(m_nSearchMax);
195 }
196
197 IOIntercept::OutputToLogMessages initLogCapture("BabaYaga@NLO", LogConfig::c_Info, LogConfig::c_Info, 100, 100);
198 initLogCapture.start();
199 m_generator.init();
200 initLogCapture.finish();
201
202 m_initialized = true;
203}
double m_eeveto
Minimum CMS energy to veto e-/e+ (GeV).
std::string m_NSKDataFile
data file for the NSK VP data.
double m_tgveto
Maximum CMS angle between the gamma and -z axis (deg).
static InitialParticleGeneration s_initial
Initial particle for beam parameters, static because we need to call it from FORTRAN here.
double m_maxAcollinearity
maximum acollinearity angle between finale state leptons in degrees.
double m_maxprescale
Maximum prescale value.
DBObjPtr< BeamParameters > m_beamParams
BeamParameter.
bool m_initialized
True if generator has been initialized.
double m_Spread
TEMPORARY SOLUTION!
std::string m_order
Module parameters.
double m_tgmin
Minimum CMS angle between the gamma and -z axis (deg).
double m_epsilon
Soft/hard photon separator in units of CMS/2.
double m_teveto
Maximum CMS theta of e-/e+ in final state (deg).
virtual void terminate() override
Method is called at the end of the event processing.
double m_temin
Minimum CMS angle between the tagged e-/e+ and -z axis (deg).
std::string m_fileNameExtraInfo
Extra ROOT file that contains the weight distribution to check overweight bias.
std::string m_userMode
User mode similar to TEEGG: ETRON, EGAMMA, GAMMA or PRESCALE or NONE.
int m_nSearchMax
Events used to search maximum of differential cross section.
int m_nPhot
fixed number of nphot (hard) photons are generated.
virtual void generatorInitialize() override
Initializes the module.
bool m_Uncertainty
vary all VP related parameters and extracted total uncertainty.
std::string m_vacPol
Vacuum polarization: off, hadr5 or hmnt.
double m_egmin
Minimum CMS energy of the gamma (GeV).
double m_eMinFrac
Fractional energy (of cms energy) for leptons in the final state.
static double generateCMSEnergy()
Static method to get a random CMS energy (used via extern C from fortran).
double m_eemin
Minimum CMS energy of the tagged e-/e+ (GeV).
std::string m_finalState
Final state: ee, mm or gg.
MCParticleGraph m_mcGraph
The MCParticle graph object.
double m_egveto
Minimum CMS energy to veto gamma (GeV).
double m_eMin
Minimum energy for leptons in the final state, in GeV.
std::vector< double > m_ScatteringAngleRange
Min [0] and Max [1] value for the scat.
void initializeGenerator()
Method is called to initialize the generator.
std::string m_model
model: matched or ps.
std::string m_mode
mode: weighted or unweighted.
virtual void generatorEvent() override
Method is called for each event.
TH1D * m_th1dSDif
Histograms with the event weights.
double m_fMax
Maximum differential cross section weight.
This class contains the nominal beam parameters and the parameters used for smearing of the primary v...
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_Info
Info: for informational messages, e.g.
Definition LogConfig.h:27
This class contains the initial state for the given event.
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
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
Class to store variables with their name which were sent to the logging service.
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.