12 #include <generators/koralw/KoralW.h>
15 #include <framework/gearbox/Unit.h>
16 #include <framework/logging/Logger.h>
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);
62 void kw_getmomdec_(
double* p1,
double* p2,
double* p3,
double* p4);
63 void kw_getxsecmc_(
double* xSecMC,
double* xErrMC);
65 void koralw_warning_ecm_(
double* ecmconfig,
double* ecm)
67 B2WARNING(
"KORALW: Different center of mass energy in config file (obsolete), E=" << *ecmconfig <<
", and from beam parameters, E="
73 void KoralW::init(
const std::string& dataPath,
const std::string& userDataFile)
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 !");
79 string dataPathNew = dataPath;
80 if (dataPath[dataPath.length() - 1] !=
'/') dataPathNew +=
"/";
83 size_t pathLength = dataPathNew.size();
84 kw_setdatapath_(dataPathNew.c_str(), &pathLength);
88 const string defaultDataFile = dataPathNew +
"KoralW_Default.data";
89 pathLength = defaultDataFile.size();
90 kw_readatax_(defaultDataFile.c_str(), &pathLength, &reset, &m_numXPar, m_XPar);
94 pathLength = userDataFile.size();
95 kw_readatax_(userDataFile.c_str(), &pathLength, &reset, &m_numXPar, m_XPar);
98 kw_initialize_(&m_cmsEnergy, m_XPar);
101 unsigned int mar1 = gRandom->Integer(m_seed1);
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);
109 void KoralW::generateEvent(
MCParticleGraph& mcGraph, TVector3 vertex, TLorentzRotation boost)
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);
118 storeParticle(mcGraph, hepevt_.
phep[iPart], hepevt_.
vhep[iPart], hepevt_.
idhep[iPart], vertex, boost);
129 kw_getxsecmc_(&m_crossSection, &m_crossSectionError);
136 void KoralW::storeParticle(
MCParticleGraph& mcGraph,
const float* mom,
const float* vtx,
int pdg, TVector3 vertex,
137 TLorentzRotation boost,
bool isVirtual,
bool isInitial)
150 part.setStatus(MCParticle::c_IsVirtual);
151 }
else if (isInitial) {
152 part.setStatus(MCParticle::c_Initial);
156 part.addStatus(MCParticle::c_PrimaryParticle);
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);
166 TLorentzVector p4 = part.get4Vector();
172 TVector3 v3 = part.getProductionVertex();
174 part.setProductionVertex(v3);
175 part.setValidVertex(
true);