Warning, /include/ThePEG/Vectors/HepMCConverter.tcc is written in an unsupported language. File is not indexed.
0001 // -*- C++ -*-
0002 //
0003 // HepMCConverter.tcc is a part of ThePEG - Toolkit for HEP Event Generation
0004 // Copyright (C) 1999-2019 Leif Lonnblad
0005 //
0006 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 //
0010 // This is the implementation of the non-inlined, non-templated member
0011 // functions of the HepMCConverter class.
0012 //
0013
0014 #include "ThePEG/StandardModel/StandardModelBase.h"
0015 #include "ThePEG/Repository/EventGenerator.h"
0016 #include "ThePEG/EventRecord/Particle.h"
0017 #include "ThePEG/EventRecord/StandardSelectors.h"
0018 #include "ThePEG/EventRecord/Collision.h"
0019 #include "ThePEG/EventRecord/Step.h"
0020 #include "ThePEG/EventRecord/SubProcess.h"
0021 #include "ThePEG/Handlers/XComb.h"
0022 #include "ThePEG/Handlers/EventHandler.h"
0023 #include "ThePEG/PDF/PartonExtractor.h"
0024 #include "ThePEG/PDF/PDF.h"
0025 #include "ThePEG/PDT/StandardMatchers.h"
0026 #include "ThePEG/Utilities/Throw.h"
0027 namespace ThePEG {
0028
0029
0030 template <typename HepMCEventT, typename Traits>
0031 typename HepMCConverter<HepMCEventT,Traits>::GenEvent *
0032 HepMCConverter<HepMCEventT,Traits>::
0033 convert(const Event & ev, bool nocopies, Energy eunit, Length lunit) {
0034 HepMCConverter<HepMCEventT,Traits> converter(ev, nocopies, eunit, lunit);
0035 return converter.geneve;
0036 }
0037
0038 template <typename HepMCEventT, typename Traits>
0039 void HepMCConverter<HepMCEventT,Traits>::
0040 convert(const Event & ev, GenEvent & gev, bool nocopies) {
0041 HepMCConverter<HepMCEventT,Traits>
0042 converter(ev, gev, nocopies,
0043 Traits::momentumUnit(gev), Traits::lengthUnit(gev));
0044 }
0045
0046 template <typename HepMCEventT, typename Traits>
0047 void HepMCConverter<HepMCEventT,Traits>::
0048 convert(const Event & ev, GenEvent & gev, bool nocopies,
0049 Energy eunit, Length lunit) {
0050 HepMCConverter<HepMCEventT,Traits> converter(ev, gev, nocopies, eunit, lunit);
0051 }
0052
0053 template <typename HepMCEventT, typename Traits>
0054 HepMCConverter<HepMCEventT,Traits>::
0055 HepMCConverter(const Event & ev, bool nocopies, Energy eunit, Length lunit)
0056 : energyUnit(eunit), lengthUnit(lunit) {
0057
0058 geneve = Traits::newEvent(ev.number(), ev.weight(), ev.optionalWeights());
0059
0060 init(ev, nocopies);
0061
0062 }
0063
0064 template <typename HepMCEventT, typename Traits>
0065 HepMCConverter<HepMCEventT,Traits>::
0066 HepMCConverter(const Event & ev, GenEvent & gev, bool nocopies,
0067 Energy eunit, Length lunit)
0068 : energyUnit(eunit), lengthUnit(lunit) {
0069
0070 geneve = &gev;
0071 Traits::resetEvent(geneve, ev.number(), ev.weight(), ev.optionalWeights());
0072
0073 init(ev, nocopies);
0074
0075 }
0076
0077 struct ParticleOrderNumberCmp {
0078 bool operator()(tcPPtr a, tcPPtr b) const {
0079 return a->number() < b->number();
0080 }
0081 };
0082
0083 template <typename HepMCEventT, typename Traits>
0084 void HepMCConverter<HepMCEventT,Traits>::init(const Event & ev, bool nocopies) {
0085
0086 if ( lengthUnit != millimeter && lengthUnit != centimeter )
0087 Throw<HepMCConverterException>()
0088 << "Length unit used for HepMC::GenEvent was not MM nor CM."
0089 << Exception::runerror;
0090 if ( energyUnit != GeV && energyUnit != MeV )
0091 Throw<HepMCConverterException>()
0092 << "Momentum unit used for HepMC::GenEvent was not GEV nor MEV."
0093 << Exception::runerror;
0094 Traits::setUnits(*geneve, energyUnit, lengthUnit);
0095
0096 tcEHPtr eh;
0097 if ( ev.primaryCollision() && ( eh =
0098 dynamic_ptr_cast<tcEHPtr>(ev.primaryCollision()->handler()) ) ) {
0099
0100 // Get general event info if present.
0101 Traits::setScaleAndAlphas(*geneve, eh->lastScale(),
0102 eh->lastAlphaS(),eh->lastAlphaEM(),
0103 energyUnit);
0104 }
0105
0106 // Extract all particles and order them.
0107 tcPVector all;
0108 ev.select(back_inserter(all), SelectAll());
0109 stable_sort(all.begin(), all.end(), ParticleOrderNumberCmp());
0110 vertices.reserve(all.size()*2);
0111
0112 // Create GenParticle's and map them to the ThePEG particles.
0113 for ( int i = 0, N = all.size(); i < N; ++i ) {
0114 tcPPtr p = all[i];
0115 if ( nocopies && p->next() ) continue;
0116 if ( pmap.find(p) != pmap.end() ) continue;
0117 pmap[p] = createParticle(p);
0118 if ( !p->children().empty() || p->next() ) {
0119 // If the particle has children it should have a decay vertex:
0120 vertices.push_back(Vertex());
0121 decv[p] = &vertices.back();
0122 vertices.back().in.insert(p);
0123 }
0124
0125 if ( !p->parents().empty() || p->previous() ||
0126 (p->children().empty() && !p->next()) ) {
0127 // If the particle has parents it should have a production
0128 // vertex. If neither parents or children it should still have a
0129 // dummy production vertex.
0130 vertices.push_back(Vertex());
0131 prov[p] = &vertices.back();
0132 vertices.back().out.insert(p);
0133 }
0134 }
0135
0136 // Now go through the the particles again, and join the vertices.
0137 for ( int i = 0, N = all.size(); i < N; ++i ) {
0138 tcPPtr p = all[i];
0139 if ( nocopies ) {
0140 if ( p->next() ) continue;
0141 for ( int i = 0, N = p->children().size(); i < N; ++i )
0142 join(p, p->children()[i]->final());
0143 tcPPtr pp = p;
0144 while ( pp->parents().empty() && pp->previous() ) pp = pp->previous();
0145 for ( int i = 0, N = pp->parents().size(); i < N; ++i )
0146 join(pp->parents()[i]->final(), p);
0147 } else {
0148 for ( int i = 0, N = p->children().size(); i < N; ++i )
0149 join(p, p->children()[i]);
0150 if ( p->next() ) join(p, p->next());
0151 for ( int i = 0, N = p->parents().size(); i < N; ++i )
0152 join(p->parents()[i], p);
0153 if ( p->previous() ) join(p->previous(), p);
0154 }
0155 }
0156
0157 // Time to create the GenVertex's
0158 for ( typename VertexMap::iterator it = prov.begin(); it != prov.end(); ++it )
0159 if ( !member(vmap, it->second) )
0160 vmap[it->second] = createVertex(it->second);
0161 for ( typename VertexMap::iterator it = decv.begin(); it != decv.end(); ++it )
0162 if ( !member(vmap, it->second) )
0163 vmap[it->second] = createVertex(it->second);
0164
0165 // First add the decay vertices for the first beam particle to avoid issue in HepMC3
0166 const Vertex * prim = decv[ev.incoming().first];
0167 if(prim) {
0168 Traits::addVertex(*geneve, vmap[prim]);
0169 vmap.erase(prim);
0170 }
0171
0172 // Now find the primary signal process vertex defined to be the
0173 // decay vertex of the first parton coming into the primary hard
0174 // sub-collision.
0175 tSubProPtr sub = ev.primarySubProcess();
0176 if ( sub && sub->incoming().first && sub->incoming().first!=ev.incoming().first) {
0177 const Vertex * prim = decv[sub->incoming().first];
0178 Traits::setSignalProcessVertex(*geneve, vmap[prim]);
0179 vmap.erase(prim);
0180 }
0181
0182 // Then add the rest of the vertices.
0183 for ( typename GenVertexMap::iterator it = vmap.begin();
0184 it != vmap.end(); ++it )
0185 Traits::addVertex(*geneve, it->second);
0186
0187 // and the incoming beam particles
0188 Traits::setBeamParticles(*geneve,pmap[ev.incoming().first],
0189 pmap[ev.incoming().second]);
0190
0191 // and the PDF info
0192 setPdfInfo(ev);
0193
0194 // and the cross section info
0195 Traits::setCrossSection(*geneve,
0196 eh->integratedXSec()/picobarn,
0197 eh->integratedXSecErr()/picobarn);
0198 for ( int i = 0, N = all.size(); i < N; ++i ) {
0199 tcPPtr p = all[i];
0200 if ( pmap.find(p) == pmap.end() ) continue;
0201 GenParticlePtrT gp = pmap[p];
0202 if ( p->hasColourInfo() ) {
0203 // Check if the particle is connected to colour lines, in which
0204 // case the lines are mapped to an integer and set in the
0205 // GenParticle's Flow info.
0206 tcColinePtr l;
0207 if ( (l = p->colourLine()) ) {
0208 if ( !member(flowmap, l) ) flowmap[l] = flowmap.size() + 500;
0209 Traits::setColourLine(*gp, 1, flowmap[l]);
0210 }
0211 if ( (l = p->antiColourLine()) ) {
0212 if ( !member(flowmap, l) ) flowmap[l] = flowmap.size() + 500;
0213 Traits::setColourLine(*gp, 2, flowmap[l]);
0214 }
0215 }
0216 if ( p->spinInfo() && p->spinInfo()->hasPolarization() ) {
0217 DPair pol = p->spinInfo()->polarization();
0218 Traits::setPolarization(*gp, pol.first, pol.second);
0219 }
0220 }
0221 }
0222
0223 template <typename HepMCEventT, typename Traits>
0224 typename HepMCConverter<HepMCEventT,Traits>::GenParticlePtrT
0225 HepMCConverter<HepMCEventT,Traits>::createParticle(tcPPtr p) const {
0226 int status = 1;
0227 size_t nChildren = p->children().size();
0228 if ( nChildren > 0 || p->next() ) status = 11;
0229 if ( nChildren > 1 ) {
0230 long id = p->data().id();
0231 if ( BaryonMatcher::Check(id) || MesonMatcher::Check(id) ||
0232 id == ParticleID::muminus || id == ParticleID::muplus ||
0233 id == ParticleID::tauminus || id == ParticleID::tauplus ) {
0234 bool child = false;
0235 for(unsigned int ix=0;ix<nChildren;++ix) {
0236 if(p->children()[ix]->id()==id) {
0237 child = true;
0238 break;
0239 }
0240 }
0241 if ( !child ) {
0242 if(p->data().widthCut()!=ZERO) {
0243 if(p->mass() <= p->data().massMax() &&
0244 p->mass() >= p->data().massMin() )
0245 status = 2;
0246 }
0247 else {
0248 status = 2;
0249 }
0250 }
0251 }
0252 }
0253 GenParticlePtrT gp =
0254 Traits::newParticle(p->momentum(), p->id(), p->status() ? p->status() : status, energyUnit);
0255
0256 if ( p->spinInfo() && p->spinInfo()->hasPolarization() ) {
0257 DPair pol = p->spinInfo()->polarization();
0258 Traits::setPolarization(*gp, pol.first, pol.second);
0259 }
0260
0261 return gp;
0262
0263 }
0264
0265 template <typename HepMCEventT, typename Traits>
0266 void HepMCConverter<HepMCEventT,Traits>::join(tcPPtr parent, tcPPtr child) {
0267 Vertex * dec = decv[parent];
0268 Vertex * pro = prov[child];
0269 if ( !pro || !dec ) Throw<HepMCConverterException>()
0270 << "Found a reference to a ThePEG::Particle which was not in the Event."
0271 << Exception::eventerror;
0272 if ( pro == dec ) return;
0273 while ( !pro->in.empty() ) {
0274 dec->in.insert(*(pro->in.begin()));
0275 decv[*(pro->in.begin())] = dec;
0276 pro->in.erase(pro->in.begin());
0277 }
0278 while ( !pro->out.empty() ) {
0279 dec->out.insert(*(pro->out.begin()));
0280 prov[*(pro->out.begin())] = dec;
0281 pro->out.erase(pro->out.begin());
0282 }
0283 }
0284
0285 template <typename HepMCEventT, typename Traits>
0286 typename HepMCConverter<HepMCEventT,Traits>::GenVertexPtrT
0287 HepMCConverter<HepMCEventT,Traits>::createVertex(Vertex * v) {
0288 if ( !v ) Throw<HepMCConverterException>()
0289 << "Found internal null Vertex." << Exception::abortnow;
0290
0291 GenVertexPtrT gv = Traits::newVertex();
0292
0293 // We assume that the vertex position is the average of the decay
0294 // vertices of all incoming and the creation vertices of all
0295 // outgoing particles in the lab. Note that this will probably not
0296 // be useful information for very small distances.
0297 LorentzPoint p;
0298 for ( tcParticleSet::iterator it = v->in.begin();
0299 it != v->in.end(); ++it ) {
0300 p += (**it).labDecayVertex();
0301 Traits::addIncoming(*gv, pmap[*it]);
0302 }
0303 for ( tcParticleSet::iterator it = v->out.begin();
0304 it != v->out.end(); ++it ) {
0305 p += (**it).labVertex();
0306 Traits::addOutgoing(*gv, pmap[*it]);
0307 }
0308
0309 p /= double(v->in.size() + v->out.size());
0310 Traits::setPosition(*gv, p, lengthUnit);
0311
0312 return gv;
0313 }
0314
0315 template <typename HepMCEventT, typename Traits>
0316 void HepMCConverter<HepMCEventT,Traits>::setPdfInfo(const Event & e) {
0317 // ids of the partons going into the primary sub process
0318 tSubProPtr sub = e.primarySubProcess();
0319 int id1 = sub->incoming().first ->id();
0320 int id2 = sub->incoming().second->id();
0321 // get the event handler
0322 tcEHPtr eh = dynamic_ptr_cast<tcEHPtr>(e.handler());
0323 // get the values of x
0324 double x1 = eh->lastX1();
0325 double x2 = eh->lastX2();
0326 // get the pdfs
0327 pair<PDF,PDF> pdfs;
0328 pdfs.first = eh->pdf<PDF>(sub->incoming().first );
0329 pdfs.second = eh->pdf<PDF>(sub->incoming().second);
0330 // get the scale
0331 Energy2 scale = eh->lastScale();
0332 // get the values of the pdfs
0333 double xf1 = pdfs.first.xfx(sub->incoming().first->dataPtr(), scale, x1);
0334 double xf2 = pdfs.second.xfx(sub->incoming().second->dataPtr(), scale, x2);
0335
0336 Traits::setPdfInfo(*geneve, id1, id2, x1, x2, sqrt(scale/GeV2), xf1, xf2);
0337
0338 }
0339
0340
0341
0342 }
0343
0344
