Belle II Software  release-05-01-25
KKGenInterface.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Kiyoshi Hayasaka
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 /* Own header. */
12 #include <generators/kkmc/KKGenInterface.h>
13 
14 /* Belle 2 headers. */
15 #include <framework/gearbox/Const.h>
16 #include <framework/gearbox/Unit.h>
17 #include <framework/logging/Logger.h>
18 #include <framework/particledb/EvtGenDatabasePDG.h>
19 #include <mdst/dataobjects/MCParticleGraph.h>
20 
21 /* ROOT headers. */
22 #include <THashList.h>
23 #include <TLorentzVector.h>
24 
25 /* C++ headers. */
26 #include <cmath>
27 #include <string>
28 #include <utility>
29 
30 using namespace Belle2;
31 
32 int KKGenInterface::setup(const std::string& KKdefaultFileName, const std::string& tauinputFileName,
33  const std::string& taudecaytableFileName, const std::string& KKMCOutputFileName)
34 {
35  B2DEBUG(20, "Begin initialisation of KKGen Interface.");
36 
37  // This is to avoid character corruption of TAUOLA output
38  char DatX_d[132], DatX_u[132], DatX_p[132], DatX_o[132];
39  for (int i = 0; i < 132; ++i) {
40  DatX_d[i] = ' ';
41  DatX_u[i] = ' ';
42  DatX_p[i] = ' ';
43  DatX_o[i] = ' ';
44  }
45  strcpy(DatX_d, KKdefaultFileName.c_str());
46  int length = strlen(DatX_d);
47  DatX_d[length] = ' ';
48  strcpy(DatX_u, tauinputFileName.c_str());
49  length = strlen(DatX_u);
50  DatX_u[length] = ' ';
51  strcpy(DatX_p, taudecaytableFileName.c_str());
52  length = strlen(DatX_p);
53  DatX_p[length] = ' ';
54  strcpy(DatX_o, KKMCOutputFileName.c_str());
55  length = strlen(DatX_o);
56  DatX_o[length] = ' ';
57  int irand = 0;
58  kk_init_(DatX_d, DatX_u, DatX_p, &irand, DatX_o);
59 
60  // seed of random generator should be set here
61  kk_init_seed_();
62 
63  // create mapping from pythia id to pdg code
64  auto particlelist = EvtGenDatabasePDG::Instance()->ParticleList();
65  m_mapPythiaIDtoPDG.reserve(particlelist->GetEntries());
66  for (TObject* object : *particlelist) {
67  EvtGenParticlePDG* particle = dynamic_cast<EvtGenParticlePDG*>(object);
68  if (!particle) {
69  B2FATAL("Something is wrong with the EvtGenDatabasePDG, got object not inheriting from EvtGenParticlePDG");
70  }
71  m_mapPythiaIDtoPDG[particle->PythiaID()] = particle->PdgCode();
72  }
73 
74  B2DEBUG(20, "End initialisation of KKGen Interface.");
75 
76  return 0;
77 }
78 
79 void KKGenInterface::set_beam_info(TLorentzVector P4_LER, double Espread_LER, TLorentzVector P4_HER, double Espread_HER)
80 {
81 
82  // Beam 4 momenta settings
83 // double crossing_angle = 0.;
84  double ph = P4_HER.Vect().Mag();
85  double pl = P4_LER.Vect().Mag();
86  double eh = P4_HER.E();
87  double el = P4_LER.E();
88  if (ph > 0. && pl > 0. && eh > 0. && el > 0.) {
89 
90  double pxh, pyh, pzh, pxl, pyl, pzl;
91  pxh = P4_HER.Px();
92  pyh = P4_HER.Py();
93  pzh = P4_HER.Pz();
94 
95  pxl = P4_LER.Px();
96  pyl = P4_LER.Py();
97  pzl = P4_LER.Pz();
98 
99  char buf[200];
100  sprintf(buf,
101  "Set Beam info: (%9.4f, %9.4f, %9.4f, %9.4f), (%9.4f, %9.4f, %9.4f, %9.4f)", pxh, pyh, pzh, eh, pxl, pyl, pzl, el);
102  B2DEBUG(100, buf);
103 
104  kk_putbeam_(&pxh, &pyh, &pzh, &eh, &pxl, &pyl, &pzl, &el);
105 
106  B2DEBUG(20, "Espread_LER=" << Espread_LER);
107  B2DEBUG(20, "Espread_HER=" << Espread_HER);
108  double Espread_CM = 0.0;
109 
110  sprintf(buf,
111  "Set Beam Energy spread: %9.4f", Espread_CM);
112  B2DEBUG(100, buf);
113  kk_begin_run_(&Espread_CM);
114  } else {
115  char buf[200];
116  sprintf(buf,
117  "Wrongly Set Beam info: Eh=%9.4f, Ph=%9.4f, El=%9.4f, Pl=%9.4f",
118  eh, ph, el, pl);
119  B2DEBUG(100, buf);
120  }
121 
122 }
123 
124 int KKGenInterface::simulateEvent(MCParticleGraph& graph, TVector3 vertex)
125 {
126  B2DEBUG(20, "Start simulation of KKGen Interface.");
127  int status = 0;
128  kk_event_(&status);
129 
130  // before storing event to MCParticle, check /hepevt/ common block
131  B2DEBUG(100, "HepEVT table:");
132 
133  for (int i = 0; i < hepevt_.nhep; ++i) {
134  char buf[200];
135  sprintf(buf,
136  "IntA: %3d %4d %8d %4d %4d %4d %9.4f %9.4f %9.4f %9.4f %9.4f",
137  i + 1, hepevt_.isthep[i], hepevt_.idhep[i],
138  hepevt_.jmohep[i][0], hepevt_.jdahep[i][0],
139  hepevt_.jdahep[i][1], hepevt_.phep[i][0],
140  hepevt_.phep[i][1], hepevt_.phep[i][2],
141  hepevt_.phep[i][3], hepevt_.phep[i][4]);
142  B2DEBUG(100, buf);
143  }
144 
145  int npar = addParticles2Graph(graph, vertex);
147  B2DEBUG(100, "GraphParticles:");
148 
149  // check MCParticleGraph
150  for (int i = 0; i < npar; ++i) {
151  MCParticleGraph::GraphParticle* p = &graph[i];
152  int moID = 0;
153  char buf[200];
154  sprintf(buf, "IntB: %3d %4u %8d %4d %4d %4d %9.4f %9.4f %9.4f %9.4f",
155  p->getIndex() , p->getStatus() , p->getPDG() , moID ,
156  p->getFirstDaughter() , p->getLastDaughter() ,
157  p->get4Vector().Px() , p->get4Vector().Py() ,
158  p->get4Vector().Pz() , p->get4Vector().E());
159  B2DEBUG(100, buf);
160 
161  }
162 
163  B2DEBUG(20, "End simulation of KKGen Interface.");
164  return hepevt_.nhep; //returns the number of generated particles from KKMC
165 }
166 
167 
168 
170 {
171  // KKMC generates at least five particles:
172  // beam (e+ e-), intermediate gamma/Z, f+ f- (f=mu, tau, ...)
173  if (hepevt_.nhep < 5) {
174  B2ERROR("KKMC-generated event has not been produced correctty!");
175  return hepevt_.nhep;
176  }
177 
178  std::vector <MCParticleGraph::GraphParticle*> MCPList;
179 
180  //Fill top particle in the tree & starting the queue:
181  for (int i = 1; i <= hepevt_.nhep; ++i) {
182  int position = graph.size();
183  graph.addParticle();
184  MCParticleGraph::GraphParticle* p = &graph[position];
185  updateGraphParticle(i, p, vertex);
186  MCPList.push_back(p);
187  }
188 
189  // From produced particles, its mother is assigned
190  for (int i = 4; i <= hepevt_.nhep; ++i) {
191  MCParticleGraph::GraphParticle* p = MCPList[i - 1];
192  if (hepevt_.jmohep[i - 1][0] > 0 && hepevt_.jmohep[i - 1][0] <= hepevt_.nhep) {
193  int j = hepevt_.jmohep[i - 1][0];
194  MCParticleGraph::GraphParticle* q = MCPList[j - 1];
195  p->comesFrom((*q));
196  }
197  }
198  return hepevt_.nhep;
199 }
200 
201 
202 void KKGenInterface::updateGraphParticle(int index, MCParticleGraph::GraphParticle* gParticle, TVector3 vertex)
203 {
204  if (index < 1 || index > hepevt_.nhep)
205  return;
206  //updating the GraphParticle information from /hepevt/ common block information
207 
208  // convert PYTHIA ID to PDG ID
209  auto iter = m_mapPythiaIDtoPDG.find(hepevt_.idhep[index - 1]);
210  if (iter != end(m_mapPythiaIDtoPDG)) {
211  gParticle->setPDG(iter->second);
212  } else {
213  //not in the map, set pythia id directly
214  gParticle->setPDG(hepevt_.idhep[index - 1]);
215  }
216 
217  //all(!) particles from the generator have to be primary
219 
220  // initial beam electron and positron should be Initial.
221  if (abs(hepevt_.idhep[index - 1]) == 11 &&
222  hepevt_.jmohep[index - 1][0] == 0 &&
223  hepevt_.jmohep[index - 1][1] == 0 &&
224  hepevt_.isthep[index - 1] == 3 &&
225  index > 0 && index < 3) {
226  gParticle->addStatus(MCParticle::c_Initial);
227  }
228 
229  // Z0 or W+/W- must be flagged as virtual
230  if (hepevt_.idhep[index - 1] == 23 || std::abs(hepevt_.idhep[index - 1]) == 24) {
232  } else if (hepevt_.isthep[index - 1] == 1) {
234  }
235 
236  // set photon flags (for now: set ISR and FSR (CEEX is undefined by definition, unsure about IFI))
237  // PHOTOS could be called by TAUOLA, not sure what to do about that
238  if (hepevt_.idhep[index - 1] == 22) {
239  if (hepevt_.jmohep[index - 1][0] == 1) {
242  }
243  }
244 
245  TLorentzVector p4(hepevt_.phep[index - 1][0],
246  hepevt_.phep[index - 1][1],
247  hepevt_.phep[index - 1][2],
248  hepevt_.phep[index - 1][3]);
249  gParticle->setMass(hepevt_.phep[index - 1][4]);
250  gParticle->set4Vector(p4);
251 
252  //set vertex including smearing (if user requested)
253  TVector3 pProductionVertex(hepevt_.vhep[index - 1][0]*Unit::mm,
254  hepevt_.vhep[index - 1][1]*Unit::mm,
255  hepevt_.vhep[index - 1][2]*Unit::mm);
256  if (!gParticle->hasStatus(MCParticle::c_Initial)) {
257  pProductionVertex = pProductionVertex + vertex;
258  }
259  gParticle->setProductionVertex(pProductionVertex);
260  gParticle->setProductionTime(hepevt_.vhep[index - 1][3]*Unit::mm / Const::speedOfLight);
261  gParticle->setValidVertex(true);
262 }
263 
265 {
266  double xsec = 0.;
267  double xsecerr = 0.;
268  kk_term_(&xsec, &xsecerr);
269 }
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:211
Belle2::MCParticleGraph::size
size_t size() const
Return the number of particles in the graph.
Definition: MCParticleGraph.h:253
Belle2::MCParticle::c_IsFSRPhoton
@ c_IsFSRPhoton
bit 7: Particle is from finial state radiation
Definition: MCParticle.h:72
Belle2::MCParticleGraph
Class to build, validate and sort a particle decay chain.
Definition: MCParticleGraph.h:48
Belle2::MCParticle::c_IsVirtual
@ c_IsVirtual
bit 4: Particle is virtual and not going to Geant4.
Definition: MCParticle.h:66
Belle2::MCParticleGraph::c_checkCyclic
@ c_checkCyclic
Check for cyclic dependencies.
Definition: MCParticleGraph.h:75
Belle2::MCParticle::c_IsISRPhoton
@ c_IsISRPhoton
bit 6: Particle is from initial state radiation
Definition: MCParticle.h:70
Belle2::KKGenInterface::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: KKGenInterface.cc:32
Belle2::KKGenInterface::term
void term()
Terminate the generator.
Definition: KKGenInterface.cc:264
Belle2::MCParticleGraph::c_setDecayInfo
@ c_setDecayInfo
Set decay time and vertex.
Definition: MCParticleGraph.h:74
Belle2::KKGenInterface::addParticles2Graph
int addParticles2Graph(MCParticleGraph &graph, TVector3 vertex)
Add particles to the MCParticleGraph.
Definition: KKGenInterface.cc:169
Belle2::MCParticle::hasStatus
bool hasStatus(unsigned short int bitmask) const
Return if specific status bit is set.
Definition: MCParticle.h:140
hepevt_type::isthep
int isthep[nmxhep]
status code.
Definition: KKGenInterface.h:28
Belle2::Const::speedOfLight
static const double speedOfLight
[cm/ns]
Definition: Const.h:568
Belle2::MCParticle::setProductionTime
void setProductionTime(float time)
Set production time.
Definition: MCParticle.h:392
Belle2::KKGenInterface::m_mapPythiaIDtoPDG
std::unordered_map< int, int > m_mapPythiaIDtoPDG
Map between PYTHIA id and PDG codes.
Definition: KKGenInterface.h:134
hepevt_type::jdahep
int jdahep[nmxhep][2]
childreen particles.
Definition: KKGenInterface.h:31
hepevt_type::phep
double phep[nmxhep][5]
four-momentum, mass [GeV].
Definition: KKGenInterface.h:32
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::MCParticle::setPDG
void setPDG(int pdg)
Set PDG code of the particle.
Definition: MCParticle.h:343
hepevt_type::jmohep
int jmohep[nmxhep][2]
parent particles.
Definition: KKGenInterface.h:30
Belle2::MCParticle::c_Initial
@ c_Initial
bit 5: Particle is initial such as e+ or e- and not going to Geant4
Definition: MCParticle.h:68
Belle2::KKGenInterface::set_beam_info
void set_beam_info(TLorentzVector P4_LER, double Espread_LER, TLorentzVector P4_HER, double Espread_HER)
Setup for beams information.
Definition: KKGenInterface.cc:79
Belle2::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:305
Belle2::MCParticle::setProductionVertex
void setProductionVertex(const TVector3 &vertex)
Set production vertex position.
Definition: MCParticle.h:404
Belle2::MCParticle::addStatus
void addStatus(unsigned short int bitmask)
Add bitmask to current status.
Definition: MCParticle.h:361
Belle2::MCParticle::c_StableInGenerator
@ c_StableInGenerator
bit 1: Particle is stable, i.e., not decaying in the generator.
Definition: MCParticle.h:60
Belle2::Unit::mm
static const double mm
[millimeters]
Definition: Unit.h:80
hepevt_type::nhep
int nhep
number of particles.
Definition: KKGenInterface.h:27
Belle2::MCParticle::setMass
void setMass(float mass)
Set particle mass.
Definition: MCParticle.h:374
Belle2::MCParticle::set4Vector
void set4Vector(const TLorentzVector &p4)
Sets the 4Vector of particle.
Definition: MCParticle.h:446
Belle2::KKGenInterface::simulateEvent
int simulateEvent(MCParticleGraph &graph, TVector3 vertex)
Simulate the events.
Definition: KKGenInterface.cc:124
Belle2::EvtGenParticlePDG
Helper class for setting additional TParticlePDG members and storing the ones it doesn't have yet.
Definition: EvtGenParticlePDG.h:32
hepevt_type::vhep
double vhep[nmxhep][4]
vertex [mm].
Definition: KKGenInterface.h:33
Belle2::MCParticle::setValidVertex
void setValidVertex(bool valid)
Set indication wether vertex and time information is valid or just default.
Definition: MCParticle.h:386
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:86
Belle2::MCParticle::c_PrimaryParticle
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition: MCParticle.h:58
Belle2::EvtGenDatabasePDG::Instance
static EvtGenDatabasePDG * Instance()
Instance method that loads the EvtGen table.
Definition: EvtGenDatabasePDG.cc:27
Belle2::KKGenInterface::updateGraphParticle
void updateGraphParticle(int, MCParticleGraph::GraphParticle *gParticle, TVector3 vertex)
Update the MCParticleGraph.
Definition: KKGenInterface.cc:202
hepevt_type::idhep
int idhep[nmxhep]
particle ident KF.
Definition: KKGenInterface.h:29