Belle II Software  release-05-02-19
PhokharaInputModule.cc
1 /************************************************************************
2 * BASF2 (Belle Analysis Framework 2) *
3 * Copyright(C) 2014 Belle II Collaboration *
4 * *
5 * Author: The Belle II Collaboration *
6 * Contributors: Torben Ferber *
7 * *
8 * This software is provided "as is" without any warranty. *
9 **************************************************************************/
10 
11 #include <generators/modules/phokharainput/PhokharaInputModule.h>
12 
13 #include <framework/datastore/StoreArray.h>
14 #include <framework/utilities/FileSystem.h>
15 #include <framework/utilities/IOIntercept.h>
16 
17 using namespace std;
18 using namespace Belle2;
19 
20 //-----------------------------------------------------------------
21 // Register the Module
22 //-----------------------------------------------------------------
23 REG_MODULE(PhokharaInput)
24 
25 //-----------------------------------------------------------------
26 // Implementation
27 //-----------------------------------------------------------------
29 {
30  //Set module properties
31  setDescription("Generates radiative return events with PHOKHARA.");
32 
33  //Parameter definition
34  addParam("FinalState", m_finalState,
35  "Final state: mu+mu-(0, default), 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)",
36  1);
37  addParam("ReplaceMuonsByVirtualPhoton", m_replaceMuonsByVirtualPhoton,
38  "Replace muons by a virtual photon (for FinalState == 0 only).",
39  false);
40  addParam("SearchMax", m_nSearchMax, "Number of events used to search for maximum of differential cross section", 100000);
41  addParam("Epsilon", m_epsilon, "Soft/hard photon separator", 0.0001);
42  addParam("nMaxTrials", m_nMaxTrials, "Maximum trials per event", 10000);
43 
44  addParam("LO", m_LO, "LO: 1ph(0, default), Born: 0ph(1), only Born: 0ph(-1)", 0);
45  addParam("NLO", m_NLO, "NLO, for 1ph only: off(0, default), on(1)", 0);
46  addParam("QED", m_QED, "ISR only(0, default), ISR+FSR(1), ISR+INT+FSR(2)", 0);
47  addParam("NLOIFI", m_NLOIFI, "IFI, NLO only: off(0, default), on(1)", 0);
48  addParam("Alpha", m_alpha, "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 
58  addParam("ScatteringAngleRangePhoton", m_ScatteringAngleRangePhoton,
59  "Min [0] and Max [1] value for the scattering angle of photons [deg], default (0, 180)", make_vector(0.0, 180.0));
60  addParam("ScatteringAngleRangeFinalStates", m_ScatteringAngleRangeFinalStates,
61  "Min [0] and Max [1] value for the scattering angle of pions(muons,nucleons,kaons) [deg], default (0, 180)", make_vector(0.0,
62  180.0));
63 
64  addParam("MinInvMassHadronsGamma", m_MinInvMassHadronsGamma, "Minimal hadrons/muons-gamma invariant mass squared [GeV^2]", 0.0);
65  addParam("MinInvMassHadrons", m_MinInvMassHadrons, "Minimal hadrons/muons invariant mass squared [GeV^2]", 0.2);
66  addParam("ForceMinInvMassHadronsCut", m_ForceMinInvMassHadronsCut,
67  "Force application of the MinInvMassHadrons cut "
68  "It is ignored by PHOKHARA with LO = 1, NLO = 1.",
69  false);
70  addParam("MaxInvMassHadrons", m_MaxInvMassHadrons, "Maximal hadrons/muons invariant mass squared [GeV^2]", 112.0);
71  addParam("MinEnergyGamma", m_MinEnergyGamma, "Minimal photon energy/missing energy, must be greater than 0.0098 * CMS energy [GeV]",
72  0.15);
73 
74  addParam("ParameterFile", m_ParameterFile, "File that contain particle properties",
75  FileSystem::findFile("/data/generators/phokhara/const_and_model_paramall9.1.dat"));
76  addParam("BeamEnergySpread", m_BeamEnergySpread,
77  "Simulate beam-energy spread (initializes PHOKHARA for every "
78  "event - very slow).", false);
79 // addParam("InputFile", m_InputFile, "File that contain input configuration", FileSystem::findFile("/data/generators/phokhara/input_9.1.dat"));
80 
81 }
82 
83 //-----------------------------------------------------------------
84 // Destructor
85 //-----------------------------------------------------------------
86 PhokharaInputModule::~PhokharaInputModule()
87 {
88 
89 }
90 
91 //-----------------------------------------------------------------
92 // Initialize
93 //-----------------------------------------------------------------
94 void PhokharaInputModule::initialize()
95 {
96  StoreArray<MCParticle> mcparticle;
97  mcparticle.registerInDataStore();
98 
99  if (m_replaceMuonsByVirtualPhoton) {
100  if (m_finalState != 0) {
101  B2FATAL("You requested to replace muons by a virtual photon, but the "
102  "final state is not mu+ mu-.");
103  }
104  if (m_QED != 0) {
105  B2FATAL("You requested to replace muons by a virtual photon. In this "
106  "mode, PHOKHARA works as an ISR generator. The parameter QED "
107  "should be set to 0 (ISR only). If FSR is taken into account "
108  "(QED = 1 or 2), the results will be incorrect.");
109  }
110  if (m_NLOIFI != 0) {
111  B2FATAL("You requested to replace muons by a virtual photon. In this "
112  "mode, PHOKHARA works as an ISR generator. The parameter NLOIFI "
113  "should be set to 0 (off). If simultaneous emission of initial "
114  "and final-state photons is taken into account (NLOIFI = 1), "
115  "the results will be incorrect.");
116  }
117  }
118  //Beam Parameters, initial particle - PHOKHARA cannot handle beam energy spread
119  if (m_BeamEnergySpread) {
120  m_initial.setAllowedFlags(BeamParameters::c_smearVertex |
121  BeamParameters::c_smearBeam);
122  }
123  m_initial.initialize();
124 
125 }
126 
127 //-----------------------------------------------------------------
128 // Event
129 //-----------------------------------------------------------------
130 void PhokharaInputModule::event()
131 {
132 
133  // Check if the BeamParameters have changed (if they do, abort the job! otherwise cross section calculation will be a nightmare.)
134  if (m_beamParams.hasChanged()) {
135  if (!m_initialized) {
136  initializeGenerator();
137  } else {
138  B2FATAL("PhokharaInputModule::event(): BeamParameters have changed within a job, this is not supported for PHOKHARA!");
139  }
140  }
141 
142  // initial particle from beam parameters
143  const MCInitialParticles& initial = m_initial.generate();
144 
145  // true boost
146  TLorentzRotation boost = initial.getCMSToLab();
147 
148  // vertex
149  TVector3 vertex = initial.getVertex();
150 
151  m_mcGraph.clear();
152  if (m_BeamEnergySpread) {
153  // PHOKHARA does not support beam-energy spread. To simulate it,
154  // the generator is initialized with the new energy for every event.
155  IOIntercept::OutputToLogMessages initLogCapture(
156  "PHOKHARA", LogConfig::c_Debug, LogConfig::c_Debug, 100, 100);
157  initLogCapture.start();
158  m_generator.setCMSEnergy(initial.getMass());
159  m_generator.init(m_ParameterFile);
160  initLogCapture.finish();
161  }
162  m_generator.generateEvent(m_mcGraph, vertex, boost);
163  m_mcGraph.generateList("", MCParticleGraph::c_setDecayInfo | MCParticleGraph::c_checkCyclic);
164 
165 }
166 
167 
168 //-----------------------------------------------------------------
169 // Terminate
170 //-----------------------------------------------------------------
171 void PhokharaInputModule::terminate()
172 {
173 
174  m_generator.term();
175 }
176 
177 void PhokharaInputModule::initializeGenerator()
178 {
179 
180  const BeamParameters& nominal = m_initial.getBeamParameters();
181  m_cmsEnergy = nominal.getMass();
182 
183  m_generator.setCMSEnergy(m_cmsEnergy);
184 
185  m_generator.setNSearchMax(m_nSearchMax);
186  m_generator.setFinalState(m_finalState);
187  m_generator.setReplaceMuonsByVirtualPhoton(m_replaceMuonsByVirtualPhoton);
188  m_generator.setNMaxTrials(m_nMaxTrials);
189  m_generator.setEpsilon(m_epsilon);
190  m_generator.setLO(m_LO);
191  m_generator.setNLO(m_NLO);
192  m_generator.setQED(m_QED);
193  m_generator.setNLOIFI(m_NLOIFI);
194  m_generator.setAlpha(m_alpha);
195  m_generator.setPionFF(m_pionff);
196  m_generator.setKaonFF(m_kaonff);
197  m_generator.setPionStructure(m_pionstructure);
198  m_generator.setNarrowRes(m_narres);
199  m_generator.setProtonFF(m_protonff);
200 
201  m_generator.setMinInvMassHadronsGamma(m_MinInvMassHadronsGamma);
202  m_generator.setm_MinInvMassHadrons(m_MinInvMassHadrons);
203  m_generator.setForceMinInvMassHadronsCut(m_ForceMinInvMassHadronsCut);
204  m_generator.setm_MaxInvMassHadrons(m_MaxInvMassHadrons);
205  m_generator.setMinEnergyGamma(m_MinEnergyGamma);
206  m_generator.setScatteringAngleRangePhoton(vectorToPair(m_ScatteringAngleRangePhoton, "ScatteringAngleRangePhoton"));
207  m_generator.setScatteringAngleRangeFinalStates(vectorToPair(m_ScatteringAngleRangeFinalStates, "ScatteringAngleRangeFinalStates"));
208 
209  m_generator.init(m_ParameterFile);
210 
211  m_initialized = true;
212 
213 }
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::IOIntercept::OutputToLogMessages::finish
bool finish()
Finish the capture and emit the message if output has appeard on stdout or stderr.
Definition: IOIntercept.cc:245
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::MCInitialParticles
This class contains the initial state for the given event.
Definition: MCInitialParticles.h:35
Belle2::MCInitialParticles::getMass
double getMass() const
Get the invariant mass of the collision (= energy in CMS)
Definition: MCInitialParticles.h:152
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::IOIntercept::InterceptOutput::start
bool start()
Start intercepting the output.
Definition: IOIntercept.h:165
Belle2::MCInitialParticles::getVertex
const TVector3 & getVertex() const
Get the position of the collision.
Definition: MCInitialParticles.h:143
Belle2::PhokharaInputModule
The Phokhara Generator module.
Definition: PhokharaInputModule.h:46
Belle2::IOIntercept::OutputToLogMessages
Capture stdout and stderr and convert into log messages.
Definition: IOIntercept.h:236
Belle2::MCInitialParticles::getCMSToLab
const TLorentzRotation & getCMSToLab() const
Return the LorentzRotation to convert from CMS to lab frame.
Definition: MCInitialParticles.h:161
Belle2::StoreArray< MCParticle >
Belle2::BeamParameters
This class contains the nominal beam parameters and the parameters used for smearing of the primary v...
Definition: BeamParameters.h:33