Belle II Software development
Public Member Functions | Private Member Functions | Private Attributes | List of all members
KKGenInterface Class Reference
Packages » generators

Interface class for using the KKMC generator. More...

#include <KKGenInterface.h>

Public Member Functions

 KKGenInterface ()
 Constructor.
 
 ~KKGenInterface ()
 Destructor.
 
 KKGenInterface (const KKGenInterface &m)=delete
 Copy constructor, explicitly deleted.
 
KKGenInterface & operator= (const KKGenInterface &m)=delete
 Assignment operator, explicitly deleted.
 
int setup (const std::string &KKdefaultFileName, const std::string &tauinputFileName, const std::string &taudecaytableFileName, const std::string &KKMCOutputFileName)
 Setup for KKMC and TAUOLA.
 
void set_beam_info (double Ecms0, double Ecms0Spread)
 Setup for beams information.
 
int simulateEvent (MCParticleGraph &graph, const ConditionalGaussGenerator &lorentzGenerator, ROOT::Math::XYZVector vertex)
 Simulate the events.
 
void term ()
 Terminate the generator.
 

Private Member Functions

int addParticles2Graph (MCParticleGraph &graph, ROOT::Math::XYZVector vertex)
 Add particles to the MCParticleGraph.
 
void updateGraphParticle (int, MCParticleGraph::GraphParticle *gParticle, ROOT::Math::XYZVector vertex)
 Update the MCParticleGraph.
 

Private Attributes

std::unordered_map< int, int > m_mapPythiaIDtoPDG
 Map between PYTHIA id and PDG codes.
 

Detailed Description

Interface class for using the KKMC generator.

Definition at line 68 of file KKGenInterface.h.

Constructor & Destructor Documentation

◆ KKGenInterface()

KKGenInterface ( )
inline

Constructor.

Definition at line 75 of file KKGenInterface.h.

75{}

◆ ~KKGenInterface()

~KKGenInterface ( )
inline

Destructor.

Definition at line 80 of file KKGenInterface.h.

80{}

Member Function Documentation

◆ addParticles2Graph()

int addParticles2Graph ( MCParticleGraph &  graph,
ROOT::Math::XYZVector  vertex 
)
private

Add particles to the MCParticleGraph.

Definition at line 179 of file KKGenInterface.cc.

180{
181 // KKMC generates at least five particles:
182 // beam (e+ e-), intermediate gamma/Z, f+ f- (f=mu, tau, ...)
183 if (hepevt_.nhep < 5) {
184 B2ERROR("KKMC-generated event has not been produced correctty!");
185 return hepevt_.nhep;
186 }
187
188 std::vector <MCParticleGraph::GraphParticle*> MCPList;
189
190 //Fill top particle in the tree & starting the queue:
191 for (int i = 1; i <= hepevt_.nhep; ++i) {
192 int position = graph.size();
193 graph.addParticle();
194 MCParticleGraph::GraphParticle* p = &graph[position];
195 updateGraphParticle(i, p, vertex);
196 MCPList.push_back(p);
197 }
198
199 // From produced particles, its mother is assigned
200 for (int i = 4; i <= hepevt_.nhep; ++i) {
201 MCParticleGraph::GraphParticle* p = MCPList[i - 1];
202 if (hepevt_.jmohep[i - 1][0] > 0 && hepevt_.jmohep[i - 1][0] <= hepevt_.nhep) {
203 int j = hepevt_.jmohep[i - 1][0];
204 MCParticleGraph::GraphParticle* q = MCPList[j - 1];
205 p->comesFrom((*q));
206 }
207 }
208 return hepevt_.nhep;
209}
Belle2::KKGenInterface::updateGraphParticle
void updateGraphParticle(int, MCParticleGraph::GraphParticle *gParticle, ROOT::Math::XYZVector vertex)
Update the MCParticleGraph.
Definition: KKGenInterface.cc:212
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:75
Belle2::MCParticleGraph::size
size_t size() const
Return the number of particles in the graph.
Definition: MCParticleGraph.h:242
Belle2::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:294
hepevt_type::jmohep
int jmohep[nmxhep][2]
parent particles.
Definition: KKGenInterface.h:29
hepevt_type::nhep
int nhep
number of particles.
Definition: KKGenInterface.h:26

◆ set_beam_info()

void set_beam_info ( double  Ecms0,
double  Ecms0Spread 
)

Setup for beams information.

Definition at line 84 of file KKGenInterface.cc.

85{
86 if (Ecms0 > 0. && Ecms0Spread >= 0.) {
87 // KKMC requires spread of energy of a single beam in CMS system
88 double E0spread = Ecms0Spread / std::sqrt(2);
89 kk_begin_run_(&Ecms0, &E0spread);
90 } else {
91 B2DEBUG(100, "Wrong beam info");
92 }
93}

◆ setup()

int setup ( const std::string &  KKdefaultFileName,
const std::string &  tauinputFileName,
const std::string &  taudecaytableFileName,
const std::string &  KKMCOutputFileName 
)

Setup for KKMC and TAUOLA.

Definition at line 36 of file KKGenInterface.cc.

38{
39 B2DEBUG(20, "Begin initialisation of KKGen Interface.");
40
41 // This is to avoid character corruption of TAUOLA output
42 char DatX_d[132], DatX_u[132], DatX_p[132], DatX_o[132];
43 for (int i = 0; i < 132; ++i) {
44 DatX_d[i] = ' ';
45 DatX_u[i] = ' ';
46 DatX_p[i] = ' ';
47 DatX_o[i] = ' ';
48 }
49 strcpy(DatX_d, KKdefaultFileName.c_str());
50 int length = strlen(DatX_d);
51 DatX_d[length] = ' ';
52 strcpy(DatX_u, tauinputFileName.c_str());
53 length = strlen(DatX_u);
54 DatX_u[length] = ' ';
55 strcpy(DatX_p, taudecaytableFileName.c_str());
56 length = strlen(DatX_p);
57 DatX_p[length] = ' ';
58 strcpy(DatX_o, KKMCOutputFileName.c_str());
59 length = strlen(DatX_o);
60 DatX_o[length] = ' ';
61 int irand = 0;
62 kk_init_(DatX_d, DatX_u, DatX_p, &irand, DatX_o);
63
64 // seed of random generator should be set here
65 kk_init_seed_();
66
67 // create mapping from pythia id to pdg code
68 auto particlelist = EvtGenDatabasePDG::Instance()->ParticleList();
69 m_mapPythiaIDtoPDG.reserve(particlelist->GetEntries());
70 for (TObject* object : *particlelist) {
71 EvtGenParticlePDG* particle = dynamic_cast<EvtGenParticlePDG*>(object);
72 if (!particle) {
73 B2FATAL("Something is wrong with the EvtGenDatabasePDG, got object not inheriting from EvtGenParticlePDG");
74 }
75 m_mapPythiaIDtoPDG[particle->PythiaID()] = particle->PdgCode();
76 }
77
78 B2DEBUG(20, "End initialisation of KKGen Interface.");
79
80 return 0;
81}
Belle2::EvtGenDatabasePDG::Instance
static EvtGenDatabasePDG * Instance()
Instance method that loads the EvtGen table.
Definition: EvtGenDatabasePDG.cc:25
Belle2::EvtGenParticlePDG
Helper class for setting additional TParticlePDG members and storing the ones it doesn't have yet.
Definition: EvtGenParticlePDG.h:22
Belle2::KKGenInterface::m_mapPythiaIDtoPDG
std::unordered_map< int, int > m_mapPythiaIDtoPDG
Map between PYTHIA id and PDG codes.
Definition: KKGenInterface.h:129

◆ simulateEvent()

int simulateEvent ( MCParticleGraph &  graph,
const ConditionalGaussGenerator &  lorentzGenerator,
ROOT::Math::XYZVector  vertex 
)

Simulate the events.

Definition at line 96 of file KKGenInterface.cc.

98{
99 B2DEBUG(20, "Start simulation of KKGen Interface.");
100 int status = 0;
101 kk_event_(&status);
102
103 ROOT::Math::PxPyPzEVector pHERorg(hepevt_.phep[0][0], hepevt_.phep[0][1], hepevt_.phep[0][2], hepevt_.phep[0][3]);
104 ROOT::Math::PxPyPzEVector pLERorg(hepevt_.phep[1][0], hepevt_.phep[1][1], hepevt_.phep[1][2], hepevt_.phep[1][3]);
105 ROOT::Math::PxPyPzEVector pTotOrg = pHERorg + pLERorg;
106
107 // KKMC allows generation of events with E-spread of beams, where
108 // without spread both energies are equal and momenta aligned along z-asis
109 // When spread it on, the system is no more CMS, so we transform to it
110 ROOT::Math::LorentzRotation rotKKMC(ROOT::Math::Boost(pTotOrg.BoostToCM()));
111
112 // CMS energy from KKMC used for conditional generator
113 double EcmsNow = pTotOrg.M();
114
115 // Calculate Lorentz transformation to LAB for this event
116 Eigen::VectorXd transVec = lorentzGenerator.generate(EcmsNow);
117 ROOT::Math::LorentzRotation rot = MCInitialParticles::cmsToLab(transVec[1], transVec[2], transVec[3], transVec[4], transVec[5]);
118
119 // Total Lorentz transformation
120 ROOT::Math::LorentzRotation rotTot = rot * rotKKMC;
121
122 for (int i = 0; i < hepevt_.nhep; ++i) {
123 ROOT::Math::PxPyPzEVector p4cms(hepevt_.phep[i][0], hepevt_.phep[i][1], hepevt_.phep[i][2], hepevt_.phep[i][3]);
124
125 // transform to LAB
126 ROOT::Math::PxPyPzEVector p4lab = rotTot * p4cms;
127
128 hepevt_.phep[i][0] = p4lab.Px();
129 hepevt_.phep[i][1] = p4lab.Py();
130 hepevt_.phep[i][2] = p4lab.Pz();
131 hepevt_.phep[i][3] = p4lab.E();
132 //hepevt_.phep[i][4] = p4lab.M();
133 }
134
135
136 // Shift vertex point due to tau life time
137 kk_shifttaudecayvtx_();
138
139
140 // before storing event to MCParticle, check /hepevt/ common block
141 B2DEBUG(100, "HepEVT table:");
142
143 for (int i = 0; i < hepevt_.nhep; ++i) {
144 char buf[200];
145 sprintf(buf,
146 "IntA: %3d %4d %8d %4d %4d %4d %9.4f %9.4f %9.4f %9.4f %9.4f",
147 i + 1, hepevt_.isthep[i], hepevt_.idhep[i],
148 hepevt_.jmohep[i][0], hepevt_.jdahep[i][0],
149 hepevt_.jdahep[i][1], hepevt_.phep[i][0],
150 hepevt_.phep[i][1], hepevt_.phep[i][2],
151 hepevt_.phep[i][3], hepevt_.phep[i][4]);
152 B2DEBUG(100, buf);
153 }
154
155 int npar = addParticles2Graph(graph, vertex);
156 graph.generateList("", MCParticleGraph::c_setDecayInfo | MCParticleGraph::c_checkCyclic);
157 B2DEBUG(100, "GraphParticles:");
158
159 // check MCParticleGraph
160 for (int i = 0; i < npar; ++i) {
161 MCParticleGraph::GraphParticle* p = &graph[i];
162 int moID = 0;
163 char buf[200];
164 sprintf(buf, "IntB: %3d %4u %8d %4d %4d %4d %9.4f %9.4f %9.4f %9.4f",
165 p->getIndex(), p->getStatus(), p->getPDG(), moID,
166 p->getFirstDaughter(), p->getLastDaughter(),
167 p->get4Vector().Px(), p->get4Vector().Py(),
168 p->get4Vector().Pz(), p->get4Vector().E());
169 B2DEBUG(100, buf);
170
171 }
172
173 B2DEBUG(20, "End simulation of KKGen Interface.");
174 return hepevt_.nhep; //returns the number of generated particles from KKMC
175}
Belle2::ConditionalGaussGenerator::generate
Eigen::VectorXd generate(double x0) const
generate random vector based on the provided first component x0
Definition: ConditionalGaussGenerator.cc:99
Belle2::KKGenInterface::addParticles2Graph
int addParticles2Graph(MCParticleGraph &graph, ROOT::Math::XYZVector vertex)
Add particles to the MCParticleGraph.
Definition: KKGenInterface.cc:179
Belle2::MCInitialParticles::cmsToLab
static ROOT::Math::LorentzRotation cmsToLab(double bX, double bY, double bZ, double angleXZ, double angleYZ)
Return the LorentzRotation from CMS to LAB based on the following parameters.
Definition: MCInitialParticles.cc:31
Belle2::MCParticleGraph::c_checkCyclic
@ c_checkCyclic
Check for cyclic dependencies.
Definition: MCParticleGraph.h:64
Belle2::MCParticleGraph::c_setDecayInfo
@ c_setDecayInfo
Set decay time and vertex.
Definition: MCParticleGraph.h:63
Belle2::MCParticleGraph::generateList
void generateList(const std::string &name="", int options=c_setNothing)
Generates the MCParticle list and stores it in the StoreArray with the given name.
Definition: MCParticleGraph.cc:209
hepevt_type::isthep
int isthep[nmxhep]
status code.
Definition: KKGenInterface.h:27
hepevt_type::jdahep
int jdahep[nmxhep][2]
childreen particles.
Definition: KKGenInterface.h:30
hepevt_type::phep
double phep[nmxhep][5]
four-momentum, mass [GeV].
Definition: KKGenInterface.h:31
hepevt_type::idhep
int idhep[nmxhep]
particle ident KF.
Definition: KKGenInterface.h:28

◆ term()

void term ( )

Terminate the generator.

Definition at line 310 of file KKGenInterface.cc.

311{
312 double xsec = 0.;
313 double xsecerr = 0.;
314 kk_term_(&xsec, &xsecerr);
315}

◆ updateGraphParticle()

void updateGraphParticle ( int  index,
MCParticleGraph::GraphParticle *  gParticle,
ROOT::Math::XYZVector  vertex 
)
private

Update the MCParticleGraph.

Definition at line 212 of file KKGenInterface.cc.

213{
214 if (index < 1 || index > hepevt_.nhep)
215 return;
216 //updating the GraphParticle information from /hepevt/ common block information
217
218 // convert PYTHIA ID to PDG ID
219 auto iter = m_mapPythiaIDtoPDG.find(hepevt_.idhep[index - 1]);
220 if (iter != end(m_mapPythiaIDtoPDG)) {
221 gParticle->setPDG(iter->second);
222 } else {
223 //not in the map, set pythia id directly
224 gParticle->setPDG(hepevt_.idhep[index - 1]);
225 }
226
227 //all(!) particles from the generator have to be primary
228 gParticle->addStatus(MCParticleGraph::GraphParticle::c_PrimaryParticle);
229
230 // initial beam electron and positron should be Initial.
231 if (abs(hepevt_.idhep[index - 1]) == 11 &&
232 hepevt_.jmohep[index - 1][0] == 0 &&
233 hepevt_.jmohep[index - 1][1] == 0 &&
234 hepevt_.isthep[index - 1] == 3 &&
235 index < 3) {
236 gParticle->addStatus(MCParticle::c_Initial);
237 }
238
239 // Z0 or W+/W- must be flagged as virtual
240 if (hepevt_.idhep[index - 1] == 23 || std::abs(hepevt_.idhep[index - 1]) == 24) {
241 gParticle->addStatus(MCParticleGraph::GraphParticle::c_IsVirtual);
242 } else if (hepevt_.isthep[index - 1] == 1) {
243 gParticle->addStatus(MCParticleGraph::GraphParticle::c_StableInGenerator);
244 }
245
246 // set photon flags
247 if (hepevt_.idhep[index - 1] == 22) {
248 // ISR and FSR from incoming beams AND at the production vertex of outgoing fermions
249 // added by KKMC using CEEX model which cannot be split into ISR and FSR due to interference
250 if (hepevt_.jmohep[index - 1][0] == 1) {
251 gParticle->addStatus(MCParticleGraph::GraphParticle::c_IsISRPhoton);
252 gParticle->addStatus(MCParticleGraph::GraphParticle::c_IsFSRPhoton);
253 } else {
254 // FSR from leptonic tau decay added by Tauola
255 // FSR from hadronic tau decay added by PHOTOS, but electromagnetic decay of hadrons by Tauola
256 // No interference with radiation from CEEX because tau decay vertices are separted due to lifetime
257 int moth = hepevt_.jmohep[index - 1][0];
258 int mothid = hepevt_.idhep[moth - 1];
259 int gmoth = hepevt_.jmohep[moth - 1][0];
260 int gmothid = hepevt_.idhep[gmoth - 1];
261 if (abs(mothid) == 15 || abs(mothid) == 24) {
262 B2DEBUG(200, "moth, mothid, gmoth, gmothid = " << moth << " " << mothid << " " << gmoth << " " << gmothid);
263 }
264 int leptonicdecay = 0;
265 int mothdau1 = hepevt_.jdahep[moth - 1][0];
266 int mothdau2 = hepevt_.jdahep[moth - 1][1];
267 for (int idau = mothdau1; idau <= mothdau2; ++idau) {
268 int dauid = abs(hepevt_.idhep[idau - 1]);
269 if (dauid == 11 || dauid == 12 || dauid == 13 || dauid == 14 || dauid == 16) {
270 leptonicdecay++;
271 }
272 }
273 if (abs(mothid) == 15 && leptonicdecay > 1) {
274 B2DEBUG(200, "mothdau1, mothdau2, leptonicdecay = " << mothdau1 << " " << mothdau2 << " " << leptonicdecay);
275 }
276 if ((abs(mothid) == 15) || // from tau decays
277 (abs(mothid) == 24 && abs(gmothid) == 15)) { // from W in tau -> W nu decays
278 if (leptonicdecay == 3) { // leptonic tau decay
279 gParticle->addStatus(MCParticleGraph::GraphParticle::c_IsFSRPhoton);
280 } else {
281 // PHOTOSPhoton flag to mark FSR photons added by PHOTOS, which cannot be separated otherwise
282 if (hepevt_.jdahep[index - 1][0] == 0 && hepevt_.jdahep[index - 1][1] == -1) {
283 gParticle->addStatus(MCParticleGraph::GraphParticle::c_IsFSRPhoton);
284 gParticle->addStatus(MCParticleGraph::GraphParticle::c_IsPHOTOSPhoton);
285 }
286 }
287 }
288 }
289 }
290
291 ROOT::Math::PxPyPzEVector p4(hepevt_.phep[index - 1][0],
292 hepevt_.phep[index - 1][1],
293 hepevt_.phep[index - 1][2],
294 hepevt_.phep[index - 1][3]);
295 gParticle->setMass(hepevt_.phep[index - 1][4]);
296 gParticle->set4Vector(p4);
297
298 //set vertex including smearing (if user requested)
299 ROOT::Math::XYZVector pProductionVertex(hepevt_.vhep[index - 1][0]*Unit::mm,
300 hepevt_.vhep[index - 1][1]*Unit::mm,
301 hepevt_.vhep[index - 1][2]*Unit::mm);
302 if (!gParticle->hasStatus(MCParticle::c_Initial)) {
303 pProductionVertex = pProductionVertex + vertex;
304 }
305 gParticle->setProductionVertex(pProductionVertex);
306 gParticle->setProductionTime(hepevt_.vhep[index - 1][3]*Unit::mm / Const::speedOfLight);
307 gParticle->setValidVertex(true);
308}
Belle2::Const::speedOfLight
static const double speedOfLight
[cm/ns]
Definition: Const.h:695
Belle2::MCParticle::c_IsFSRPhoton
@ c_IsFSRPhoton
bit 7: Particle is from finial state radiation
Definition: MCParticle.h:61
Belle2::MCParticle::c_Initial
@ c_Initial
bit 5: Particle is initial such as e+ or e- and not going to Geant4
Definition: MCParticle.h:57
Belle2::MCParticle::c_IsPHOTOSPhoton
@ c_IsPHOTOSPhoton
bit 8: Particle is an radiative photon from PHOTOS
Definition: MCParticle.h:63
Belle2::MCParticle::c_PrimaryParticle
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition: MCParticle.h:47
Belle2::MCParticle::c_IsVirtual
@ c_IsVirtual
bit 4: Particle is virtual and not going to Geant4.
Definition: MCParticle.h:55
Belle2::MCParticle::c_StableInGenerator
@ c_StableInGenerator
bit 1: Particle is stable, i.e., not decaying in the generator.
Definition: MCParticle.h:49
Belle2::MCParticle::c_IsISRPhoton
@ c_IsISRPhoton
bit 6: Particle is from initial state radiation
Definition: MCParticle.h:59
Belle2::MCParticle::setMass
void setMass(float mass)
Set particle mass.
Definition: MCParticle.h:366
Belle2::MCParticle::addStatus
void addStatus(unsigned short int bitmask)
Add bitmask to current status.
Definition: MCParticle.h:353
Belle2::MCParticle::hasStatus
bool hasStatus(unsigned short int bitmask) const
Return if specific status bit is set.
Definition: MCParticle.h:129
Belle2::MCParticle::setValidVertex
void setValidVertex(bool valid)
Set indication wether vertex and time information is valid or just default.
Definition: MCParticle.h:378
Belle2::MCParticle::setProductionVertex
void setProductionVertex(const ROOT::Math::XYZVector &vertex)
Set production vertex position.
Definition: MCParticle.h:396
Belle2::MCParticle::setPDG
void setPDG(int pdg)
Set PDG code of the particle.
Definition: MCParticle.h:335
Belle2::MCParticle::set4Vector
void set4Vector(const ROOT::Math::PxPyPzEVector &p4)
Sets the 4Vector of particle.
Definition: MCParticle.h:438
Belle2::MCParticle::setProductionTime
void setProductionTime(float time)
Set production time.
Definition: MCParticle.h:384
Belle2::Unit::mm
static const double mm
[millimeters]
Definition: Unit.h:70
hepevt_type::vhep
double vhep[nmxhep][4]
vertex [mm].
Definition: KKGenInterface.h:32

Member Data Documentation

◆ m_mapPythiaIDtoPDG

std::unordered_map<int, int> m_mapPythiaIDtoPDG
private

Map between PYTHIA id and PDG codes.

Definition at line 129 of file KKGenInterface.h.

The documentation for this class was generated from the following files: