11 #include <framework/gearbox/Unit.h>
12 #include <framework/gearbox/Const.h>
13 #include <framework/logging/Logger.h>
14 #include <framework/particledb/EvtGenDatabasePDG.h>
15 #include <framework/utilities/FileSystem.h>
16 #include <generators/evtgen/EvtGenInterface.h>
17 #include <mdst/dataobjects/MCParticleGraph.h>
18 #include <generators/evtgen/EvtGenModelRegister.h>
24 #include <EvtGenExternal/EvtExternalGenList.hh>
25 #include <EvtGenBase/EvtAbsRadCorr.hh>
26 #include <EvtGenBase/EvtCPUtil.hh>
27 #include <EvtGenBase/EvtDecayTable.hh>
28 #include <EvtGenBase/EvtDecayBase.hh>
29 #include <EvtGenBase/EvtParticleFactory.hh>
30 #include <EvtGenBase/EvtPDL.hh>
31 #include <EvtGenBase/EvtRandom.hh>
33 #include <TLorentzVector.h>
40 EvtGenInterface::~EvtGenInterface()
42 EvtDecayTable* evtDecayTable = EvtDecayTable::getInstance();
43 for (
unsigned int i = 0; i < EvtPDL::entries(); ++i) {
44 for (
int j = 0; j < evtDecayTable->getNMode(i); ++j) {
45 delete evtDecayTable->getDecay(i, j);
48 if (m_Generator)
delete m_Generator;
51 EvtGen* EvtGenInterface::createEvtGen(
const std::string& DECFileName,
const bool coherentMixing)
54 initLogCapture.
start();
55 EvtRandom::setRandomEngine(&EvtGenInterface::m_eng);
58 std::list<EvtDecayBase*> extraModels = EvtGenModelRegister::getModels();
61 EvtExternalGenList genList;
62 EvtAbsRadCorr* radCorrEngine = genList.getPhotosModel();
63 list<EvtDecayBase*> modelList = genList.getListOfModels();
64 extraModels.insert(extraModels.end(), modelList.begin(), modelList.end());
67 EvtGenDatabasePDG::Instance()->WriteEvtGenTable(tmp);
69 auto mixingMode = EvtCPUtil::Incoherent;
71 mixingMode = EvtCPUtil::Coherent;
73 B2WARNING(
"Evtgen has been set to decay the B mesons incoherently. This is useful as a workaround only to generate Y(5S, 6S) -> BBar for QCD studies.");
74 B2WARNING(
"If you are generating Y(4S) events, you _really_ must use the coherent decay mode.");
77 return new EvtGen(DECFileName.c_str(), tmp.getName().c_str(), &EvtGenInterface::m_eng,
78 radCorrEngine, &extraModels, mixingMode);
83 int EvtGenInterface::setup(
const std::string& DECFileName,
const std::string& parentParticle,
84 const std::string& userFileName,
const bool coherentMixing)
86 B2DEBUG(150,
"Begin initialisation of EvtGen Interface.");
90 initLogCapture.
start();
92 m_Generator = createEvtGen(DECFileName, coherentMixing);
94 if (!userFileName.empty()) {
95 m_Generator->readUDecay(userFileName.c_str());
99 if (parentParticle !=
"") {
100 m_ParentInitialized =
true;
101 m_ParentParticle = EvtPDL::getId(parentParticle);
105 B2DEBUG(150,
"End initialisation of EvtGen Interface.");
111 int EvtGenInterface::simulateEvent(
MCParticleGraph& graph, TLorentzVector pParentParticle, TVector3 pPrimaryVertex,
112 int inclusiveType,
const std::string& inclusiveParticle)
114 EvtId inclusiveParticleID, inclusiveAntiParticleID;
116 if (!m_ParentInitialized)
117 B2FATAL(
"Parent particle is not initialized.");
119 m_pinit.set(pParentParticle.E(), pParentParticle.X(), pParentParticle.Y(), pParentParticle.Z());
121 if (inclusiveType != 0) {
122 inclusiveParticleID = EvtPDL::getId(inclusiveParticle);
123 if (inclusiveParticleID.getId() < 0)
124 B2FATAL(
"Incorrect inclusive particle " << inclusiveParticle);
125 inclusiveAntiParticleID = EvtPDL::chargeConj(inclusiveParticleID);
126 if (inclusiveAntiParticleID.getId() < 0) {
127 B2FATAL(
"Cannot find the charge-conjugate particle for "
128 << inclusiveParticle);
132 bool we_got_inclusive_particle =
false;
134 m_logCapture.start();
135 m_parent = EvtParticleFactory::particleFactory(m_ParentParticle, m_pinit);
136 m_parent->setVectorSpinDensity();
137 m_Generator->generateDecay(m_parent);
138 m_logCapture.finish();
140 if (inclusiveType != 0) {
141 EvtParticle* p = m_parent;
147 if (p->getId() == inclusiveParticleID ||
148 (inclusiveType == 2 && p->getId() == inclusiveAntiParticleID)) {
149 we_got_inclusive_particle =
true;
152 p = p->nextIter(m_parent);
155 if (!we_got_inclusive_particle) {
156 m_parent->deleteTree();
160 we_got_inclusive_particle =
true;
162 }
while (!we_got_inclusive_particle);
166 int iPart = addParticles2Graph(m_parent, graph, pPrimaryVertex, NULL);
167 graph.
generateList(
"", MCParticleGraph::c_setDecayInfo | MCParticleGraph::c_checkCyclic);
169 m_parent->deleteTree();
179 TLorentzVector momentum = parent.
get4Vector();
181 m_pinit.set(momentum.E(), momentum.X(), momentum.Y(), momentum.Z());
182 m_logCapture.start();
185 parent.
setDecayTime(-std::numeric_limits<float>::infinity());
187 id = EvtPDL::evtIdFromStdHep(pdg);
188 m_parent = EvtParticleFactory::particleFactory(
id, m_pinit);
190 m_parent->setVectorSpinDensity();
192 m_parent->setDiagonalSpinDensity();
193 m_Generator->generateDecay(m_parent);
194 m_logCapture.finish();
195 int iPart = addParticles2Graph(m_parent, graph, vertex, &parent, parent.
getProductionTime());
196 m_parent->deleteTree();
200 int EvtGenInterface::addParticles2Graph(EvtParticle* top,
MCParticleGraph& graph, TVector3 pPrimaryVertex,
204 const int existingParticles = graph.
size();
210 updateGraphParticle(top, p, pPrimaryVertex);
212 typedef pair<MCParticleGraph::GraphParticle*, EvtParticle*> halfFamily;
213 halfFamily currFamily;
214 queue < halfFamily > heritancesQueue;
216 for (uint idaughter = 0; idaughter < top->getNDaug(); idaughter++) {
217 currFamily.first = p;
218 currFamily.second = top->getDaug(idaughter);
219 heritancesQueue.push(currFamily);
223 while (!heritancesQueue.empty()) {
224 currFamily = heritancesQueue.front();
225 heritancesQueue.pop();
228 EvtParticle* currDaughter = currFamily.second;
232 updateGraphParticle(currDaughter, graphDaughter, pPrimaryVertex, timeOffset);
237 int nGrandChildren = currDaughter->getNDaug();
239 if (nGrandChildren == 0)
240 graphDaughter->
addStatus(MCParticle::c_StableInGenerator);
242 for (
int igrandchild = 0; igrandchild < nGrandChildren; igrandchild++) {
243 currFamily.first = graphDaughter;
244 currFamily.second = currDaughter->getDaug(igrandchild);
245 heritancesQueue.push(currFamily);
250 return graph.
size() - existingParticles;
255 TVector3 pPrimaryVertex,
double timeOffset)
259 gParticle->
setStatus(MCParticle::c_PrimaryParticle);
260 gParticle->
setMass(eParticle->mass());
261 gParticle->
setPDG(EvtPDL::getStdHep(eParticle->getId()));
263 EvtVector4R EvtP4 = eParticle->getP4Lab();
264 TLorentzVector p4(EvtP4.get(1), EvtP4.get(2), EvtP4.get(3), EvtP4.get(0));
267 EvtVector4R Evtpos = eParticle->get4Pos();
269 TVector3 pVertex(Evtpos.get(1)*Unit::mm, Evtpos.get(2)*Unit::mm, Evtpos.get(3)*Unit::mm);
270 pVertex = pVertex + pPrimaryVertex;
273 gParticle->
setProductionTime((Evtpos.get(0)*Unit::mm / Const::speedOfLight) + timeOffset);
277 if (eParticle->getAttribute(
"FSR")) {
278 gParticle->
addStatus(MCParticle::c_IsPHOTOSPhoton);