Belle II Software  release-05-01-25
KoralW.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2011 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Andreas Moll *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 /* Own header. */
12 #include <generators/koralw/KoralW.h>
13 
14 /* Belle 2 headers. */
15 #include <framework/gearbox/Unit.h>
16 #include <framework/logging/Logger.h>
17 
18 /* ROOT headers. */
19 #include <TRandom.h>
20 
21 using namespace std;
22 using namespace Belle2;
23 
24 //Declaration of the KoralW FORTRAN methods and common blocks.
25 extern"C" {
26 
27  extern struct {
28  double q1[4];
29  double q2[4];
30  double p1[4];
31  double p2[4];
32  double p3[4];
33  double p4[4];
34  } momdec_;
35 
36  extern struct {
37  double qeff1[4];
38  double qeff2[4];
39  double sphum[4];
40  double sphot[4][100];
41  int nphot;
42  } momset_;
43 
44  extern struct {
45  int nevhep;
46  int nhep;
47  int isthep[2000];
48  int idhep[2000];
49  int jmohep[2000][2];
50  int jdahep[2000][2];
51  float phep[2000][5];
52  float vhep[2000][4];
53  } hepevt_;
54 
55  void kw_setdatapath_(const char* filename, size_t* length);
56  void kw_readatax_(const char* filename, size_t* length, int* reset, const int* max, double* xpar);
57  void kw_initialize_(double* ecm, double* xpar);
58  void marini_(unsigned int* mar1, unsigned int* mar2, unsigned int* mar3);
59  void rmarin_(unsigned int* mar1, unsigned int* mar2, unsigned int* mar3);
60  void kw_make_();
61  void kw_finalize_();
62  void kw_getmomdec_(double* p1, double* p2, double* p3, double* p4);
63  void kw_getxsecmc_(double* xSecMC, double* xErrMC);
64 
65  void koralw_warning_ecm_(double* ecmconfig, double* ecm)
66  {
67  B2WARNING("KORALW: Different center of mass energy in config file (obsolete), E=" << *ecmconfig << ", and from beam parameters, E="
68  << *ecm);
69  }
70 
71 }
72 
73 void KoralW::init(const std::string& dataPath, const std::string& userDataFile)
74 {
75  if (dataPath.empty()) B2FATAL("KoralW: The specified data path is empty !");
76  if (userDataFile.empty()) B2FATAL("KoralW: The specified user data file is empty !");
77 
78  //Make sure the dataPath ends with an "/"
79  string dataPathNew = dataPath;
80  if (dataPath[dataPath.length() - 1] != '/') dataPathNew += "/";
81 
82  //Set the path to the data files
83  size_t pathLength = dataPathNew.size();
84  kw_setdatapath_(dataPathNew.c_str(), &pathLength);
85 
86  //Load the default KoralW input data
87  int reset = 1;
88  const string defaultDataFile = dataPathNew + "KoralW_Default.data";
89  pathLength = defaultDataFile.size();
90  kw_readatax_(defaultDataFile.c_str(), &pathLength, &reset, &m_numXPar, m_XPar);
91 
92  //Load the user KoralW input data
93  reset = 0;
94  pathLength = userDataFile.size();
95  kw_readatax_(userDataFile.c_str(), &pathLength, &reset, &m_numXPar, m_XPar);
96 
97  //Initialize KoralW
98  kw_initialize_(&m_cmsEnergy, m_XPar);
99 
100  //Initialize random number generator
101  unsigned int mar1 = gRandom->Integer(m_seed1); //The range for this seed seems to be [0, 900000000]
102  unsigned int mar2 = gRandom->Integer(m_seed2);
103  unsigned int mar3 = gRandom->Integer(m_seed3);
104  marini_(&mar1, &mar2, &mar3);
105  rmarin_(&mar1, &mar2, &mar3);
106 }
107 
108 
109 void KoralW::generateEvent(MCParticleGraph& mcGraph, TVector3 vertex, TLorentzRotation boost)
110 {
111  kw_make_();
112 
113  //Store the particles with status id 3 and 2 as virtual particles, the particles with status id 1 as real particles
114  for (int iPart = 0; iPart < hepevt_.nhep; ++iPart) {
115  if (hepevt_.isthep[iPart] > 1) {
116  storeParticle(mcGraph, hepevt_.phep[iPart], hepevt_.vhep[iPart], hepevt_.idhep[iPart], vertex, boost, true);
117  } else {
118  storeParticle(mcGraph, hepevt_.phep[iPart], hepevt_.vhep[iPart], hepevt_.idhep[iPart], vertex, boost);
119  }
120  }
121 }
122 
123 
124 void KoralW::term()
125 {
126  kw_finalize_();
127 
128  //Get the cross section
129  kw_getxsecmc_(&m_crossSection, &m_crossSectionError);
130 }
131 
132 //=========================================================================
133 // Protected methods
134 //=========================================================================
135 
136 void KoralW::storeParticle(MCParticleGraph& mcGraph, const float* mom, const float* vtx, int pdg, TVector3 vertex,
137  TLorentzRotation boost, bool isVirtual, bool isInitial)
138 {
139  // //Create particle
140  //MCParticleGraph::GraphParticle& part = mcGraph.addParticle();
141  //if (!isVirtual) {
142  // part.setStatus(MCParticle::c_PrimaryParticle);
143  //} else {
144  // part.setStatus(MCParticle::c_IsVirtual);
145  //}
146 
147  // RG 6/25/14 Add new flag for ISR "c_Initial"
149  if (isVirtual) {
150  part.setStatus(MCParticle::c_IsVirtual);
151  } else if (isInitial) {
152  part.setStatus(MCParticle::c_Initial);
153  }
154 
155  //every particle from a generator is primary, TF
156  part.addStatus(MCParticle::c_PrimaryParticle);
157  part.setPDG(pdg);
158  part.setFirstDaughter(0);
159  part.setLastDaughter(0);
160  part.setMomentum(TVector3(mom[0], mom[1], mom[2]));
161  part.setEnergy(mom[3]);
162  part.setMass(mom[4]);
163  part.setProductionVertex(vtx[0]*Unit::mm, vtx[1]*Unit::mm, vtx[2]*Unit::mm);
164 
165  //boost, TF
166  TLorentzVector p4 = part.get4Vector();
167  p4 = boost * p4;
168  part.set4Vector(p4);
169 
170  //set vertex, TF
171  if (!isInitial) {
172  TVector3 v3 = part.getProductionVertex();
173  v3 = v3 + vertex;
174  part.setProductionVertex(v3);
175  part.setValidVertex(true);
176  }
177 }
Belle2::MCParticleGraph
Class to build, validate and sort a particle decay chain.
Definition: MCParticleGraph.h:48
hepevt_type::isthep
int isthep[nmxhep]
status code.
Definition: KKGenInterface.h:28
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::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:305
hepevt_type::nhep
int nhep
number of particles.
Definition: KKGenInterface.h:27
hepevt_type::vhep
double vhep[nmxhep][4]
vertex [mm].
Definition: KKGenInterface.h:33
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:86
hepevt_type::idhep
int idhep[nmxhep]
particle ident KF.
Definition: KKGenInterface.h:29