12 #include <generators/kkmc/KKGenInterface.h>
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>
22 #include <THashList.h>
23 #include <TLorentzVector.h>
33 const std::string& taudecaytableFileName,
const std::string& KKMCOutputFileName)
35 B2DEBUG(20,
"Begin initialisation of KKGen Interface.");
38 char DatX_d[132], DatX_u[132], DatX_p[132], DatX_o[132];
39 for (
int i = 0; i < 132; ++i) {
45 strcpy(DatX_d, KKdefaultFileName.c_str());
46 int length = strlen(DatX_d);
48 strcpy(DatX_u, tauinputFileName.c_str());
49 length = strlen(DatX_u);
51 strcpy(DatX_p, taudecaytableFileName.c_str());
52 length = strlen(DatX_p);
54 strcpy(DatX_o, KKMCOutputFileName.c_str());
55 length = strlen(DatX_o);
58 kk_init_(DatX_d, DatX_u, DatX_p, &irand, DatX_o);
66 for (TObject*
object : *particlelist) {
69 B2FATAL(
"Something is wrong with the EvtGenDatabasePDG, got object not inheriting from EvtGenParticlePDG");
74 B2DEBUG(20,
"End initialisation of KKGen Interface.");
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.) {
90 double pxh, pyh, pzh, pxl, pyl, pzl;
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);
104 kk_putbeam_(&pxh, &pyh, &pzh, &eh, &pxl, &pyl, &pzl, &el);
106 B2DEBUG(20,
"Espread_LER=" << Espread_LER);
107 B2DEBUG(20,
"Espread_HER=" << Espread_HER);
108 double Espread_CM = 0.0;
111 "Set Beam Energy spread: %9.4f", Espread_CM);
113 kk_begin_run_(&Espread_CM);
117 "Wrongly Set Beam info: Eh=%9.4f, Ph=%9.4f, El=%9.4f, Pl=%9.4f",
126 B2DEBUG(20,
"Start simulation of KKGen Interface.");
131 B2DEBUG(100,
"HepEVT table:");
133 for (
int i = 0; i < hepevt_.
nhep; ++i) {
136 "IntA: %3d %4d %8d %4d %4d %4d %9.4f %9.4f %9.4f %9.4f %9.4f",
140 hepevt_.
phep[i][1], hepevt_.
phep[i][2],
141 hepevt_.
phep[i][3], hepevt_.
phep[i][4]);
147 B2DEBUG(100,
"GraphParticles:");
150 for (
int i = 0; i < npar; ++i) {
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());
163 B2DEBUG(20,
"End simulation of KKGen Interface.");
173 if (hepevt_.
nhep < 5) {
174 B2ERROR(
"KKMC-generated event has not been produced correctty!");
178 std::vector <MCParticleGraph::GraphParticle*> MCPList;
181 for (
int i = 1; i <= hepevt_.
nhep; ++i) {
182 int position = graph.
size();
186 MCPList.push_back(p);
190 for (
int i = 4; i <= hepevt_.
nhep; ++i) {
192 if (hepevt_.
jmohep[i - 1][0] > 0 && hepevt_.
jmohep[i - 1][0] <= hepevt_.
nhep) {
193 int j = hepevt_.
jmohep[i - 1][0];
204 if (index < 1 || index > hepevt_.
nhep)
211 gParticle->
setPDG(iter->second);
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) {
230 if (hepevt_.
idhep[index - 1] == 23 || std::abs(hepevt_.
idhep[index - 1]) == 24) {
232 }
else if (hepevt_.
isthep[index - 1] == 1) {
238 if (hepevt_.
idhep[index - 1] == 22) {
239 if (hepevt_.
jmohep[index - 1][0] == 1) {
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]);
253 TVector3 pProductionVertex(hepevt_.
vhep[index - 1][0]*
Unit::mm,
257 pProductionVertex = pProductionVertex + vertex;
268 kk_term_(&xsec, &xsecerr);