Belle II Software development
TrepsInputModule.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/trepsinput/TrepsInputModule.h>
10#include <framework/geometry/B2Vector3.h>
11#include <string.h>
12#include <string>
13#include <TFile.h>
14#include <TRandom.h>
15
16#include <map>
17
18using namespace Belle2;
19
20//-----------------------------------------------------------------
21// Register the Module
22//-----------------------------------------------------------------
23REG_MODULE(TrepsInput);
24
25//-----------------------------------------------------------------
26// Implementation
27//-----------------------------------------------------------------
28
30{
31 // Set module properties
32 setDescription("Input from TREPS generator (No-tag), Input from TREPS generator for ee->ee hadrons");
33
34 // Parameter definitions
35 addParam("ParameterFile", m_parameterFile,
36 "parameter file for TREPS input", std::string("treps_par.dat"));
37 addParam("DifferentialCrossSectionFile", m_differentialCrossSectionFile,
38 "file name for differential cross section table input. If UseDiscreteAndSortedW is true, the file is used",
39 std::string("pipidcs.dat"));
40 addParam("WListTableFile", m_wListTableFile,
41 "file name for W-List table input. If UseDiscreteAndSortedW is false (default), the file is used", std::string("wlist_table.dat"));
42 addParam("UseDiscreteAndSortedW", m_useDiscreteAndSortedW,
43 "if true, use WListTable for discrete and sorted W. if false (default), use DifferentialCrossSection", false);
44
45 addParam("MaximalQ2", m_maximalQ2,
46 "Maximal :math:`Q^2 = -q^2`, where q is the difference between the initial "
47 "and final electron or positron momenta. Negative means no cut.",
48 -1.0);
49 addParam("MaximalAbsCosTheta", m_maximalAbsCosTheta,
50 "Maximal :math:`|\\cos(\\theta)|`, where :math:`\\theta` is the final-state particle "
51 "polar angle.", 1.01);
52 addParam("ApplyCosThetaCutCharged", m_applyCosThetaCutCharged,
53 "Whether to apply cut on :math:`|cos(theta)|` for charged particles only.",
54 true);
55 addParam("MinimalTransverseMomentum", m_minimalTransverseMomentum,
56 "Minimal transverse momentum of the final-state particles.", 0.0);
57 addParam("ApplyTransverseMomentumCutCharged",
59 "Whether to apply cut on the minimal transverse momentum for "
60 "charged particles only.", true);
61
62}
63
65{
66 m_initial.initialize();
67 m_mcparticles.registerInDataStore();
68}
69
71{
72 /*
73 * Check if the BeamParameters have changed (if they do, abort the job;
74 * otherwise cross section calculation will be a nightmare).
75 */
76 if (m_beamParams.hasChanged()) {
77 if (!m_initialized) {
79 B2INFO("Initialized the TREPS generator!");
80 } else {
81 B2FATAL("TrepsInputModule::event(): BeamParameters have changed within "
82 "a job, this is not supported for TREPS!");
83 }
84 }
85
86 /* Generation of the initial particle from beam parameters. */
87 const MCInitialParticles& initial = m_initial.generate();
88 ROOT::Math::XYZVector vertex = initial.getVertex();
89
91 if (m_generator.inmode != 0) return;
92
93 m_generator.wtcount++;
94 m_generator.wf = (float)(m_generator.wtable(1));
95
96 if (abs((double)m_generator.wf - m_generator.w) >= 0.001 && m_generator.wf > 0.01) {
97 B2INFO(" W value changed. " << m_generator.w << " to " << m_generator.wf);
98 m_generator.w = (double)m_generator.wf;
99 m_generator.updateW();
100 }
101 } else {
103 m_generator.updateW();
104 }
105
106 int idummy = 0;
107 int iret = m_generator.event_gen(idummy);
108 m_mpg.clear();
109
110 // fill data of the initial electron and positron
111 MCParticleGraph::GraphParticle& electron = m_mpg.addParticle();
113 electron.setPDG(11);
114 electron.setMomentum(ROOT::Math::XYZVector(m_generator.getElectronMomentum()));
115 electron.setMass(m_generator.me);
116 electron.setEnergy(sqrt(m_generator.me * m_generator.me + m_generator.getElectronMomentum().Mag2()));
117
118 MCParticleGraph::GraphParticle& positron = m_mpg.addParticle();
120 positron.setPDG(-11);
121 positron.setMomentum(ROOT::Math::XYZVector(m_generator.getPositronMomentum()));
122 positron.setMass(m_generator.me);
123 positron.setEnergy(sqrt(m_generator.me * m_generator.me + m_generator.getPositronMomentum().Mag2()));
124
125 if (iret >= 1) {
126
127 const Part_gen* part = m_generator.partgen;
128
129 // fill data of the final-state particles
130 for (int i = 0; i < m_generator.npart ; i++) {
131 auto& p = m_mpg.addParticle();
132 p.setPDG(part[i].part_prop.icode);
133 p.set4Vector(ROOT::Math::PxPyPzEVector(part[i].p.Px(), part[i].p.Py(), part[i].p.Pz(), part[i].p.E()));
134 p.setMass(part[i].part_prop.pmass);
136 p.setProductionVertex(vertex);
137 p.setValidVertex(true);
138 }
139 // fill data of the recoil electron and positron
140 auto& p1 = m_mpg.addParticle();
141 p1.setPDG(11);
142 p1.set4Vector(ROOT::Math::PxPyPzEVector(m_generator.pe.Px(), m_generator.pe.Py(), m_generator.pe.Pz(), m_generator.pe.E()));
143 p1.setMass(m_generator.me);
145 p1.setProductionVertex(vertex);
146 p1.setValidVertex(true);
147
148 auto& p2 = m_mpg.addParticle();
149 p2.setPDG(-11);
150 p2.set4Vector(ROOT::Math::PxPyPzEVector(m_generator.pp.Px(), m_generator.pp.Py(), m_generator.pp.Pz(), m_generator.pp.E()));
151 p2.setMass(m_generator.me);
153 p2.setProductionVertex(vertex);
154 p2.setValidVertex(true);
155
156 }
157 //Fill MCParticle List
158 m_mpg.generateList(m_mcparticles.getName(), MCParticleGraph::c_setDecayInfo);
159
160}
161
165
167{
168 if (m_generator.diffCrossSectionOfW.size() == 0) {
169 B2FATAL("Cross Section Table is empty !!!");
170 return 0.;
171 }
172
173 // lower_bound returns first iterator which meets >=W condition. --> upper side
174 auto it_upper = m_generator.diffCrossSectionOfW.lower_bound(W);
175 auto it_lower = it_upper;
176 it_lower--;
177
178 return (it_upper->second - it_lower->second) / (it_upper->first - it_lower->first) * (W - it_lower->first) + it_lower->second;
179}
180
182{
183
184 while (1) {
185 double crossSectionForMC = gRandom->Uniform(0.0, m_generator.totalCrossSectionForMC);
186
187 auto it_upper = m_generator.WOfCrossSectionForMC.lower_bound(crossSectionForMC);
188 auto it_lower = it_upper;
189 it_lower--;
190
191 double diffCrossSectionAtUpper = m_generator.diffCrossSectionOfW.at(it_upper->second);
192 double diffCrossSectionAtLower = m_generator.diffCrossSectionOfW.at(it_lower->second);
193 double limit = (diffCrossSectionAtUpper > diffCrossSectionAtLower) ? diffCrossSectionAtUpper * 1.01 : diffCrossSectionAtLower *
194 1.01;
195 // Higher diffCrossSection * 1.01 is set as limit. Same definition of limit is used in TrepsB::wtable().
196
197 double W = (crossSectionForMC - it_lower->first) / limit + it_lower->second ;
198 if (W < m_generator.diffCrossSectionOfW.begin()->first or W > m_generator.diffCrossSectionOfW.rbegin()->first)
199 B2FATAL("W has to be in [" << m_generator.diffCrossSectionOfW.begin()->first << ", "
200 << m_generator.diffCrossSectionOfW.rbegin()->first
201 << "] !!! W = " << W << ", crossSectionForMC = " << crossSectionForMC);
202
203 double trial = gRandom->Uniform(0.0, limit);
204 double crossSection = getCrossSection(W);
205
206 if (trial < crossSection)
207 return W;
208
209 }
210
211 return 0;
212}
213
214
216{
217 // Set parameter files
218 m_generator.setParameterFile(m_parameterFile);
219 m_generator.setDiffcrosssectionFile(m_differentialCrossSectionFile);
220 m_generator.setWlistFile(m_wListTableFile);
221
222 // Set cut parameters
223 m_generator.setMaximalQ2(m_maximalQ2);
224 m_generator.setMaximalAbsCosTheta(m_maximalAbsCosTheta);
225 m_generator.applyCosThetaCutCharged(m_applyCosThetaCutCharged);
226 m_generator.setMinimalTransverseMomentum(m_minimalTransverseMomentum);
227 m_generator.applyTransverseMomentumCutCharged(m_applyTransverseMomentumCutCharged);
228
229 // Initialize the initial particle information
230 const BeamParameters& nominalBeam = m_initial.getBeamParameters();
231 m_generator.setBeamEnergy(nominalBeam.getMass() / 2.);
232 B2Vector3D p3 = nominalBeam.getHER().Vect();
233 m_generator.setElectronMomentum(p3);
234 p3 = B2Vector3D(nominalBeam.getLER().Vect());
235 m_generator.setPositronMomentum(p3);
236
237 // Initialize generator;
238 m_generator.initp();
239
240 m_generator.create_hist();
241 m_generator.initg();
242
244 // Initialize wtable with WListTable
245 B2INFO("Discrete W-list is used !!!");
246
247 m_generator.wtcount = 0;
248 m_generator.wtable(0);
249
250 m_generator.w = (double)m_generator.wf;
251
252 m_generator.updateW();
253 } else {
254 // Initialize wtable with DifferentialCrossSection
255 m_generator.wtable();
256 }
257
258 m_initialized = true;
259}
This class contains the nominal beam parameters and the parameters used for smearing of the primary v...
This class contains the initial state for the given event.
const ROOT::Math::PxPyPzEVector & getHER() const
Get 4vector of the high energy beam.
const ROOT::Math::PxPyPzEVector & getLER() const
Get 4vector of the low energy beam.
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)
Class to represent Particle data in graph.
@ c_setDecayInfo
Set decay time and vertex.
@ c_Initial
bit 5: Particle is initial such as e+ or e- and not going to Geant4
Definition MCParticle.h:57
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition MCParticle.h:47
@ c_StableInGenerator
bit 1: Particle is stable, i.e., not decaying in the generator.
Definition MCParticle.h:49
void setMass(float mass)
Set particle mass.
Definition MCParticle.h:355
void addStatus(unsigned short int bitmask)
Add bitmask to current status.
Definition MCParticle.h:342
void setEnergy(float energy)
Set energy.
Definition MCParticle.h:361
void setPDG(int pdg)
Set PDG code of the particle.
Definition MCParticle.h:324
void setMomentum(const ROOT::Math::XYZVector &momentum)
Set particle momentum.
Definition MCParticle.h:406
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
std::string m_parameterFile
Parameter file which configures the setting of beam, production particle etc.
UtrepsB m_generator
Treps generator.
std::string m_differentialCrossSectionFile
Differential cross section table file.
double m_maximalQ2
Maximal Q^2 = -q^2, where q is the difference between the initial and final electron or positron mome...
DBObjPtr< BeamParameters > m_beamParams
BeamParameter.
bool m_initialized
True if generator has been initialized.
double getCrossSection(double W)
Returns the differential cross section for given W [GeV].
bool m_useDiscreteAndSortedW
if it is true, W-list table is used for discrete and sorted W
double m_minimalTransverseMomentum
Minimal transverse momentum of the final-state particles.
StoreArray< MCParticle > m_mcparticles
MCParticle collection.
virtual void terminate() override
initialization for trepsinput
bool m_applyCosThetaCutCharged
Whether to apply cut on |cos(theta)| for charged particles only.
virtual void generatorInitialize() override
initialization for trepsinput
double simulateW()
Simulate W distribution according to given input file of cross section table.
TrepsInputModule()
Constructor: Sets the description, the properties and the parameters of the module.
double m_maximalAbsCosTheta
Maximal |cos(theta)|, where theta is the final-state particle polar angle.
InitialParticleGeneration m_initial
initial particle used by BeamParameter class
MCParticleGraph m_mpg
An instance of the MCParticle graph.
std::string m_wListTableFile
W-List table file.
void initializeGenerator()
Initialize the TREPS generator.
bool m_applyTransverseMomentumCutCharged
Whether to apply cut on the minimal transverse momentum for charged particles only.
virtual void generatorEvent() override
input event from TREPS
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
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition B2Vector3.h:516
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Abstract base class for different kinds of events.