diff --git a/EventRecord/Particle.cc b/EventRecord/Particle.cc --- a/EventRecord/Particle.cc +++ b/EventRecord/Particle.cc @@ -1,524 +1,525 @@ // -*- C++ -*- // // Particle.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined member functions of // the Particle class. // #include "Particle.h" #include "ThePEG/EventRecord/Step.h" #include "ThePEG/EventRecord/Event.h" #include "ThePEG/EventRecord/ColourLine.h" #include "ThePEG/Utilities/Rebinder.h" #include "ThePEG/Config/algorithm.h" #include "ThePEG/EventRecord/ParticleTraits.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/DecayMode.h" #include #include #include #ifdef ThePEG_TEMPLATES_IN_CC_FILE #include "Particle.tcc" #endif using namespace ThePEG; Particle::ParticleRep::ParticleRep(const ParticleRep & p) : theParents(p.theParents), theChildren(p.theChildren), thePrevious(p.thePrevious), theNext(p.theNext), theBirthStep(p.theBirthStep), theVertex(p.theVertex), theLifeLength(p.theLifeLength), theScale(p.theScale), theVetoScale(p.theVetoScale), theNumber(p.theNumber), theExtraInfo(p.theExtraInfo.size()) { if ( p.theColourInfo ) theColourInfo = dynamic_ptr_cast(p.theColourInfo->clone()); if ( p.theSpinInfo ) theSpinInfo = dynamic_ptr_cast(p.theSpinInfo->clone()); for ( int i = 0, N = p.theExtraInfo.size(); i < N; ++i ) theExtraInfo[i] = p.theExtraInfo[i]->clone(); } Particle::Particle(const Particle & p) - : Base(p), theData(p.theData), theMomentum(p.theMomentum), theRep(p.theRep) { + : Base(p), theData(p.theData), theMomentum(p.theMomentum), + theRep(p.theRep), theStatus(p.theStatus) { if ( p.theRep ) { theRep = new ParticleRep(*p.theRep); theRep->theParents.clear(); } } Particle::~Particle() { if ( theRep ) { if ( colourLine() ) colourLine()->removeColoured(this); if ( antiColourLine() ) antiColourLine()->removeAntiColoured(this); delete theRep; } theRep = 0; theData = cEventPDPtr(); } void Particle::initFull() { if ( theRep ) return; theRep = new ParticleRep; Energy width = data().generateWidth(mass()); if ( width > ZERO ) { Time lifetime = data().generateLifeTime(mass(), width); theRep->theLifeLength.setTau(lifetime); theRep->theLifeLength. setVect((momentum().vect()*(lifetime / max(mass(), Constants::epsilon*GeV)))); theRep->theLifeLength.rescaleEnergy(); } } PPtr Particle::clone() const { return ptr_new(*this); } void Particle::rebind(const EventTranslationMap & trans) { for ( ParticleVector::iterator pit = rep().theChildren.begin(); pit != rep().theChildren.end(); ++pit ) *pit = trans.translate(*pit); for ( tParticleVector::iterator pit = rep().theParents.begin(); pit != rep().theParents.end(); ++pit ) *pit = trans.translate(*pit); rep().thePrevious = trans.translate(rep().thePrevious); rep().theNext = trans.translate(rep().theNext); if ( hasColourInfo() ) colourInfo()->rebind(trans); if ( spinInfo() ) spinInfo()->rebind(trans); rep().theBirthStep = trans.translate(rep().theBirthStep); for ( EIVector::const_iterator ie = rep().theExtraInfo.begin(); ie != rep().theExtraInfo.end(); ++ie ) (**ie).rebind(trans); } tParticleSet Particle::siblings() const { tParticleSet theSiblings; for ( tParticleVector::const_iterator pit = parents().begin(); pit != parents().end(); ++pit ) theSiblings.insert((*pit)->children().begin(), (*pit)->children().end()); theSiblings.erase(const_cast(this)); return theSiblings; } void Particle::colourNeighbour(tPPtr p, bool anti) { tColinePtr line = colourLine(!anti); if ( !line ) line = ColourLine::create(this, !anti); line->addColoured(p, anti); } void Particle::outgoingColour(tPPtr p, bool anti) { tColinePtr line = colourLine(anti); if ( !line ) line = ColourLine::create(this, anti); line->addColoured(p, anti); } tPPtr Particle::incomingColour(bool anti) const { if ( !hasColourInfo() ) return tPPtr(); tColinePtr line = colourLine(anti); if ( !line ) return tPPtr(); for ( int i = 0, N = parents().size(); i < N; ++i ) if ( parents()[i]->hasColourLine(line, anti) ) return parents()[i]; return tPPtr(); } tPPtr Particle::outgoingColour(bool anti) const { if ( !hasColourInfo() ) return tPPtr(); tColinePtr line = colourLine(anti); if ( !line ) return tPPtr(); for ( int i = 0, N = children().size(); i < N; ++i ) if ( children()[i]->hasColourLine(line, anti) ) return children()[i]; return tPPtr(); } LorentzPoint Particle::labVertex() const { LorentzPoint r(rep().theBirthStep && rep().theBirthStep->collision()? vertex() + rep().theBirthStep->collision()->vertex(): vertex()); return r; } void Particle::setLabVertex(const LorentzPoint & p) { rep().theVertex = ( rep().theBirthStep && rep().theBirthStep->collision()? p - rep().theBirthStep->collision()->vertex() : p ); } void Particle::transform(const LorentzRotation & r) { if ( hasRep() && spinInfo() ) spinInfo()->transform(momentum(), r); theMomentum.transform(r); if ( !hasRep() ) return; rep().theVertex.transform(r); rep().theLifeLength.transform(r); } void Particle::deepTransform(const LorentzRotation & r) { transform(r); if ( !theRep ) return; for ( int i = 0, N = children().size(); i < N; ++i ) rep().theChildren[i]->deepTransform(r); if ( rep().theNext ) rep().theNext->deepTransform(r); } void Particle::rotateX(double a) { LorentzRotation r; r.rotateX(a); transform(r); } void Particle::deepRotateX(double a) { LorentzRotation r; r.rotateX(a); deepTransform(r); } void Particle::rotateY(double a) { LorentzRotation r; r.rotateY(a); transform(r); } void Particle::deepRotateY(double a) { LorentzRotation r; r.rotateY(a); deepTransform(r); } void Particle::rotateZ(double a) { LorentzRotation r; r.rotateZ(a); transform(r); } void Particle::deepRotateZ(double a) { LorentzRotation r; r.rotateZ(a); deepTransform(r); } void Particle::rotate(double a, const Axis & axis) { LorentzRotation r; r.rotate(a, axis); transform(r); } void Particle::deepRotate(double a, const Axis & axis) { LorentzRotation r; r.rotate(a, axis); deepTransform(r); } string Particle::outputFormat = "%n3%s10 %i7 %p[,]0 %c(,) %^^0%vv0 %>>0%<>0 %l{,}0\n" " %x10.3%y10.3%z10.3%e10.3%m10.3\n"; int getNumber(string::const_iterator & pos, int def) { if ( !isdigit(*pos) ) return def; def = *pos++ - '0'; while ( isdigit(*pos) ) def = 10*def + *pos++ - '0'; return def; } void writePrecision(ostream & os, string::const_iterator & pos, int defw, int defp, double x) { defw = getNumber(pos, defw); if ( *pos == '.' ) defp = getNumber(++pos, defp); int oldp = os.precision(); os << setprecision(defp) << setw(defw) << x << setprecision(oldp); } void writeStringAdjusted(ostream & os, bool left, int w, string str) { while ( !left && w-- > int(str.size()) ) os << ' '; os << str; while ( left && w-- > int(str.size()) ) os << ' '; } template void writeParticleRanges(ostream & os, const Container & co, char sep, int w) { set cnum; for ( typename Container::const_iterator it = co.begin(); it != co.end(); ++it) cnum.insert((**it).number()); bool elipsis = false; int last = -10; for ( set::iterator it = cnum.begin(); it != cnum.end(); ++it) { int n = *it; int next = 0; set::iterator itn = it; if ( ++itn != cnum.end() ) next = *itn; bool writeit = true; bool writesep = false; if ( elipsis && ( n != last + 1 || n != next - 1 ) ) elipsis = false; else if ( !elipsis && n == last + 1 && n == next -1 ) { os << ".."; elipsis = true; writeit = false; } else if ( elipsis && n == last + 1 && n == next -1 ) writeit = false; else if ( it != cnum.begin() ) writesep = true; if ( writeit ) { if ( writesep ) os << sep; os << setw(w) << n; } last = n; } } ostream & ThePEG::operator<<(ostream & os, const Particle & p) { return p.print(os, p.birthStep()); } ostream & Particle::print(ostream & os, tcStepPtr step) const { if ( !step ) step = birthStep(); tCollPtr coll = step? step->collision(): tCollPtr(); tEventPtr event = coll? coll->event(): tEventPtr(); string::const_iterator pos = Particle::outputFormat.begin(); ios::fmtflags saveflags = os.setf(ios::fixed, ios::floatfield); while ( pos != Particle::outputFormat.end() ) { if ( *pos == '%' && ++pos != Particle::outputFormat.end() ) { bool left = false; if ( *pos == '-' ) { left = true; os.setf(ios::left, ios::adjustfield); ++pos; } else { os.setf(ios::right, ios::adjustfield); } char mark; char open; char close; char sep; int w; string str; string fill; if ( pos == Particle::outputFormat.end() ) break; bool fullColour = false; switch ( *pos ) { case 'n': os << setw(getNumber(++pos, 3)) << number(); break; case 'i': os << setw(getNumber(++pos, 8)) << id(); break; case 's': writeStringAdjusted(os, left, getNumber(++pos, 8), PDGName()); break; case 'x': writePrecision(os, ++pos, 10, 3, momentum().x()/GeV); break; case 'y': writePrecision(os, ++pos, 10, 3, momentum().y()/GeV); break; case 'z': writePrecision(os, ++pos, 10, 3, momentum().z()/GeV); break; case 'e': writePrecision(os, ++pos, 10, 3, momentum().e()/GeV); break; case 'm': writePrecision(os, ++pos, 10, 3, momentum().mass()/GeV); break; case 'P': fullColour = true; [[fallthrough]]; case 'p': open = *++pos; sep = *++pos; close = *++pos; w = getNumber(++pos, 0); if ( parents().empty() ) break; if ( open ) os << open; writeParticleRanges(os, parents(), sep, w); if ( fullColour && hasColourInfo() && ( incomingColour() || incomingAntiColour() ) ) { if ( close ) os << open; if ( incomingColour() ) os << "+" << incomingColour()->number(); if ( incomingAntiColour() ) os << "-" << incomingAntiColour()->number(); if ( close ) os << close; } if ( close ) os << close; break; case 'l': open = *++pos; sep = *++pos; close = *++pos; w = getNumber(++pos, 0); if ( hasColourInfo() && ( colourLine() || antiColourLine() ) && event) { if ( open ) os << open; vector clines = colourInfo()->colourLines(); for ( int i = 0, N = clines.size(); i < N; ++i ) { if ( i > 0 && sep ) os << sep; clines[i]->write(os, event, false); } vector aclines = colourInfo()->antiColourLines(); for ( int i = 0, N = aclines.size(); i < N; ++i ) { if ( ( i > 0 || clines.size() ) && sep ) os << sep; aclines[i]->write(os, event, true); } if ( close ) os << close; } break; case 'C': fullColour = true; [[fallthrough]]; case 'c': open = *++pos; sep = *++pos; close = *++pos; w = getNumber(++pos, 0); if ( children().empty() ) break; if ( open ) os << open; writeParticleRanges(os, children(), sep, w); if ( fullColour && hasColourInfo() && ( outgoingColour() || outgoingAntiColour() ) ) { if ( close ) os << open; if ( outgoingColour() ) os << "+" << outgoingColour()->number(); if ( outgoingAntiColour() ) os << "-" << outgoingAntiColour()->number(); if ( close ) os << close; } if ( close ) os << close; break; case '>': mark = *++pos; w = getNumber(++pos, 0); if ( hasColourInfo() && step && step->colourNeighbour(this) ) { os << setw(w-1) << step->colourNeighbour(this)->number() << mark; } break; case '<': mark = *++pos; w = getNumber(++pos, 0); if ( hasColourInfo() && step && step->antiColourNeighbour(this) ) { int n = step->antiColourNeighbour(this)->number(); ostringstream oss; oss << mark << n; writeStringAdjusted(os, left, w, oss.str()); } break; case 'v': mark = *++pos; w = getNumber(++pos, 0); if ( next() ) { if ( left && mark ) os << mark; os << setw(w) << next()->number(); if ( !left && mark ) os << mark; } break; case '^': mark = *++pos; w = getNumber(++pos, 0); if ( previous() ) { if ( left && mark ) os << mark; os << setw(w) << previous()->number(); if ( !left && mark ) os << mark; } break; case 'd': switch ( *++pos ) { case 'x': writePrecision(os, ++pos, 10, 3, lifeLength().x()/mm); break; case 'y': writePrecision(os, ++pos, 10, 3, lifeLength().y()/mm); break; case 'z': writePrecision(os, ++pos, 10, 3, lifeLength().z()/mm); break; case 't': writePrecision(os, ++pos, 10, 3, lifeLength().e()/mm); break; case 'T': writePrecision(os, ++pos, 10, 3, lifeLength().tau()/mm); break; } break; case 'V': switch ( *++pos ) { case 'x': writePrecision(os, ++pos, 10, 3, vertex().x()/mm); break; case 'y': writePrecision(os, ++pos, 10, 3, vertex().y()/mm); break; case 'z': writePrecision(os, ++pos, 10, 3, vertex().z()/mm); break; case 't': writePrecision(os, ++pos, 10, 3, vertex().e()/mm); break; } break; case 'L': switch ( *++pos ) { case 'x': writePrecision(os, ++pos, 10, 3, labVertex().x()/mm); break; case 'y': writePrecision(os, ++pos, 10, 3, labVertex().y()/mm); break; case 'z': writePrecision(os, ++pos, 10, 3, labVertex().z()/mm); break; case 't': writePrecision(os, ++pos, 10, 3, labVertex().e()/mm); break; } break; default: os << *pos++; } } else { if ( pos != Particle::outputFormat.end() ) os << *pos++; } } os.flags(saveflags); return os; } void Particle::debugme() const { cerr << *this; EventRecordBase::debugme(); } void Particle::persistentOutput(PersistentOStream & os) const { EventConfig::putParticleData(os, theData); - os << ounit(theMomentum, GeV) << bool( theRep != 0 ); + os << ounit(theMomentum, GeV) << theStatus << bool( theRep != 0 ); if ( !theRep ) return; os << rep().theParents << rep().theChildren << rep().thePrevious << rep().theNext << rep().theBirthStep << ounit(rep().theVertex, mm) << ounit(rep().theLifeLength, mm) << ounit(rep().theScale, GeV2) << ounit(rep().theVetoScale, GeV2) << rep().theNumber << rep().theDecayMode << rep().theColourInfo << rep().theSpinInfo << rep().theExtraInfo; } void Particle::persistentInput(PersistentIStream & is, int) { bool readRep; EventConfig::getParticleData(is, theData); - is >> iunit(theMomentum, GeV) >> readRep; + is >> iunit(theMomentum, GeV) >> theStatus >> readRep; if ( !readRep ) return; if ( !hasRep() ) theRep = new ParticleRep; is >> rep().theParents >> rep().theChildren >> rep().thePrevious >> rep().theNext >> rep().theBirthStep >> iunit(rep().theVertex, mm) >> iunit(rep().theLifeLength, mm) >> iunit(rep().theScale, GeV2) >> iunit(rep().theVetoScale, GeV2) >> rep().theNumber >> rep().theDecayMode >> rep().theColourInfo >> rep().theSpinInfo >> rep().theExtraInfo; } ClassDescription Particle::initParticle; void Particle::Init() {} diff --git a/EventRecord/Particle.h b/EventRecord/Particle.h --- a/EventRecord/Particle.h +++ b/EventRecord/Particle.h @@ -1,1180 +1,1195 @@ // -*- C++ -*- // // Particle.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef ThePEG_Particle_H #define ThePEG_Particle_H // This is the decalaration of the Particle class. #include "EventConfig.h" #include "ThePEG/Vectors/Lorentz5Vector.h" #include "ThePEG/Vectors/LorentzRotation.h" #include "ThePEG/Utilities/ClassDescription.h" #include "ThePEG/EventRecord/MultiColour.h" #include "ThePEG/EventRecord/SpinInfo.h" #include "ThePEG/PDT/ParticleData.h" namespace ThePEG { /** * The Particle class is used to describe an instance of a * particle. Properties of the corresponding particle type can be * accessed through a pointer to a ParticleData object. * * A Particle object contains pointers to other particles, such as a * list of parents and a list of children. It may also contain a * pointer to the previous or next instance of the same physical * particle if the properties of a given particle has been changed * during the generation. Coloured particles contains pointers to * ColourLine defining the colour connections to other particles. * * The Particle also has a pointer to the Step object where it was * first introduced in the Event. * * When printing a particle the format of the output is governed by * the static outputFormat string. When a particle is sent * to an ostream, the format string is written but with * keys prefixed by the \% character replaced with * infromation about the particle as follows:
\%\% is * replaced by a singel \%
\%C sets a flag so * that subsequent output of children and parents etc. will contain * colour information.
\%n is replaced by the particles * number in a fied ow width
\%s is replaced by the name * of the particle
\%i is replaced by the id number of * the particle type
\%x, \%y, \%z, \%e, \%m is replaced by * the x-, y-, z-, energy- and mass-component of the particles * momentum respectively
\%dx, \%dy, \%dz, \%dt, \%dT is * replaced by the x-, y-, z-, time- and invariant time-component of * the particles lifeLength respectively
\%Vx, \%Vy, \%Vz, * \%Vt is replaced by the x-, y-, z- and time-component of the * creation point relative to the vertex of the * collision.
\%Lx, \%Ly, \%Lz, \%Lt is replaced by the x-, * y-, z- and time-component of the creation point in the current * lab-system.
\%p[,] is replaced by a list (of numbers) * of the particles parents enclosed by [ and ] * and separated by ,, The parent from which the particle * inherits its (anti-) colour is marked by a (-)+
* \%c(,) is replaced by a list (of numbers) of the * particles children enclosed by ( and ) and * separated by ,, The child which inherits the particles * (anti-) colour is marked by a (-)+
\%> is replaced * by the number of the colour neighbor
\%< is * replaced by the number of the anti-colour neighbor
* \%^ is replaced by the number of the previous instance of * the same physical particle
\%v is replaced by the * number of the next instance of the same physical particle.
* \%l{,} is replaced by the indices of the colour lines to * which this particle is connected enclosed by { and * } and separated by ,, The line corresponding * to the (anti-) colour of the particle is prefixed by a (-)+ * * @see Event * @see Collision * @see Step * @see SubProcess * @see Lorentz5Vector * @see ColourLine * @see ColourBase */ class Particle: public EventRecordBase { public: /** Most of the Event classes are friends with each other. */ friend class Event; /** Most of the Event classes are friends with each other. */ friend class Collision; /** Most of the Event classes are friends with each other. */ friend class Step; /** Most of the Event classes are friends with each other. */ friend class SubProcess; /** ParticleData needs to be a friend. */ friend class ParticleData; struct ParticleRep; public: /** * Standard Constructor. Note that the default constructor is * private - there is no particle without a pointer to a * ParticleData object. */ - Particle(tcEventPDPtr newData) : theData(newData), theRep(0) {} + Particle(tcEventPDPtr newData) : theData(newData), theRep(0), theStatus(0) {} /** * Copy constructor. */ Particle(const Particle &); /** * Destructor. */ virtual ~Particle(); //@} /** @name Functions relating to ancestry of particles. */ //@{ /** * Returns true if and only if this particle has decayed. */ bool decayed() const { return hasRep() && !rep().theChildren.empty(); } /** * The list of decay products. */ const ParticleVector & children() const { static const ParticleVector null; return hasRep() ? rep().theChildren : null; } /** * Add a child (the childs parent pointer will be set accordingly). */ void addChild(tPPtr c) { rep().theChildren.push_back(c); (c->rep()).theParents.push_back(this); } /** * Remove the given child from the list of children of this particle * (the corresponding parent pointer of the child will also be * removed). */ void abandonChild(tPPtr child) { removeChild(child); child->removeParent(this); } /** * The list of parent particles. */ const tParticleVector & parents() const { static const tParticleVector null; return hasRep() ? rep().theParents : null; } /** * Return a set of neighboring particles coming from the same decay * as this one. The return value is a newly recalculated set * every time. It must be stored to be used further, do not directly call * e.g. siblings().begin() or siblings().end()! */ tParticleSet siblings() const; /** * Undo the decay of this particle, removing all children (and * grand children ...) from the event record */ void undecay() { if ( hasRep() ) { rep().theChildren.clear(); rep().theNext = tPPtr(); } } /** * If this particle has decayed set the corresponding decay mode. */ void decayMode(tDMPtr dm) { rep().theDecayMode = dm; } /** * If this particle has decayed get the corresponding decay mode. */ tDMPtr decayMode() const { return hasRep() ? rep().theDecayMode : tDMPtr(); } /** * Next instance. Pointer to another instance of the same * physical particle in later steps. */ tPPtr next() const { return hasRep() ? rep().theNext : PPtr(); } /** * Previous instance. Pointer to another instance of the same * physical particle in earlier steps. */ tPPtr previous() const { return hasRep() ? rep().thePrevious : tPPtr(); } /** * Original instance. If there exists another previous instance of * this particle return this instance (recursively). */ tcPPtr original() const { return previous() ? tcPPtr(previous()->original()) : tcPPtr(this); } /** * Original instance. If there exists another previous instance of * this particle return this instance (recursively). */ tPPtr original() { return previous() ? previous()->original() : tPPtr(this); } /** * Final instance. If there exists another subsequent instance of * this particle return this instance (recursively). */ tcPPtr final() const { return next() ? tcPPtr(next()->final()) : tcPPtr(this); } /** * Final instance. If there exists another subsequent instance of * this particle return this instance (recursively). */ tPPtr final() { return next() ? next()->final() : tPPtr(this); } //@} /** @name Relations to the Event and Step. */ //@{ /** * Get the first Step object where this particle occurred. */ tStepPtr birthStep() const { return hasRep() ? rep().theBirthStep : tStepPtr(); } /** * Get the order-number for this particle in the current event. */ int number() const { return hasRep() ? rep().theNumber : 0; } + + /** + * Get the status code of the particle + */ + int status() const { return theStatus; } + + /** + * Set the status code of the particle + */ + void status(int n) { theStatus = n; } //@} /** @name Access the underlying ParticleData object. */ //@{ /** * Access the ParticleData object of this particle type */ const ParticleDataClass & data() const { return *theData; } /** * Access the ParticleData object of this particle type */ tcEventPDPtr dataPtr() const { return theData; } /** * Return the PDG name of this particle. */ const string & PDGName() const { return data().PDGName(); } /** * Return the PDG id number of this particle. */ long id() const { return data().id(); } //@} /** @name Functions to access the momentum. */ //@{ /** * Return the momentum of this particle. */ const Lorentz5Momentum & momentum() const { return theMomentum; } /** * Set the 3-momentum of this particle. The energy is set to be * consistent with the mass. */ void set3Momentum(const Momentum3 & p) { theMomentum.setVect(p); theMomentum.rescaleEnergy(); } /** * Set the momentum of this particle. Afterwards, the underlying * Lorentz5Momentum may have inconsistent mass. */ void setMomentum(const LorentzMomentum & p) { theMomentum = p; } /** * Set the momentum and mass. */ void set5Momentum(const Lorentz5Momentum & p) { theMomentum = p; } /** * Acces the mass of this particle. */ Energy mass() const { return momentum().mass(); } /** * Acces the mass of this particle type. */ Energy nominalMass() const { return data().mass(); } /** * Get the scale at which this particle is considered resolved. */ Energy2 scale() const { return hasRep() ? rep().theScale : -1.0*GeV2; } /** * Set the scale at which this particle is considered resolved. */ void scale(Energy2 q2) { rep().theScale = q2; } /** * Get the scale above which this particle should * not radiate. */ Energy2 vetoScale() const { return hasRep() ? rep().theVetoScale : -1.0*GeV2; } /** * Set the scale above which this particle should * not radiate. */ void vetoScale(Energy2 q2) { rep().theVetoScale = q2; } /** * Return the transverse mass (squared), calculated from the energy * and the longitudinal momentum. */ Energy2 mt2() const { return sqr(momentum().t()) - sqr(momentum().z()); } /** * Return the transverse mass (squared), calculated from the energy * and the longitudinal momentum. */ Energy mt() const { return sqrt(mt2()); } /** * Return the transverse mass (squared), calculated from the mass * and the transverse momentum. */ Energy2 perpmass2() const { return momentum().perp2() + momentum().mass2(); } /** * Return the transverse mass (squared), calculated from the mass * and the transverse momentum. */ Energy perpmass() const { return sqrt(perpmass2()); } /** * Return the (pseudo) rapidity. */ double rapidity() const { return ( Pplus() > ZERO && Pminus() > ZERO )? 0.5*log(Pplus()/Pminus()) : Constants::MaxFloat; } /** * Return the (pseudo) rapidity. */ double eta() const { Energy rho = momentum().rho(); return rho > abs(momentum().z())? 0.5*log((rho+momentum().z())/(rho-momentum().z())) : Constants::MaxFloat; } /** * Return the positive and negative light-cone momenta. */ Energy Pplus() const { return momentum().plus(); } /** * Return the positive and negative light-cone momenta. */ Energy Pminus() const { return momentum().minus(); } //@} /** @name Functions to access the position. */ //@{ /** * The creation vertex of this particle. The point is given * relative to the collision vertex. */ const LorentzPoint & vertex() const { static const LorentzPoint null; return hasRep() ? rep().theVertex : null; } /** * The creation vertex of this particle. The absolute * position in the lab is given. */ LorentzPoint labVertex() const; /** * The decay vertex of this particle. The point is given * relative to the collision vertex. */ LorentzPoint decayVertex() const { return vertex() + lifeLength(); } /** * The decay vertex of this particle. The absolute * position in the lab is given. */ LorentzPoint labDecayVertex() const { return labVertex() + lifeLength(); } /** * The life time/length. Return the Lorentz vector connecting the * creation to the decay vertes. */ const Lorentz5Distance & lifeLength() const { static const Lorentz5Distance null; return hasRep() ? rep().theLifeLength : null; } /** * Set the creation vertex relative to the collision vertex. */ void setVertex(const LorentzPoint & p) { rep().theVertex = p; } /** * Set the creation vertex in the lab frame of this particle. */ void setLabVertex(const LorentzPoint &); /** * Set the life length of this particle. The life time will be * automatically rescaled to be consistent with the invariant * distance. */ void setLifeLength(const Distance & d) { rep().theLifeLength.setVect(d); rep().theLifeLength.rescaleEnergy(); } /** * Set the life time/length of a particle. The invariant distance * may become inconsistent. */ void setLifeLength(const LorentzDistance & d) { rep().theLifeLength = d; } /** * Set the life time/length of a particle. */ void setLifeLength(const Lorentz5Distance & d) { rep().theLifeLength = d; } /** * The invariant life time of this particle. */ Time lifeTime() const { return lifeLength().m(); } //@} /** @name Functions for (Lorentz) transformations. */ //@{ /** * Do Lorentz transformations on this particle. */ void transform(const LorentzRotation & r); /** * Do Lorentz transformations on this particle. \a bx, \a by and \a * bz are the boost vector components. */ void boost(double bx, double by, double bz) { transform(LorentzRotation(Boost(bx, by, bz))); } /** * Do Lorentz transformations on this particle. \a b is the boost * vector. */ void boost(const Boost & b) { transform(LorentzRotation(b)); } /** * Rotate around the x-axis. */ void rotateX(double a); /** * Rotate around the y-axis. */ void rotateY(double a); /** * Rotate around the z-axis. */ void rotateZ(double a); /** * Rotate around the given \a axis. */ void rotate(double a, const Axis & axis); /** * Mirror in the xy-plane. */ void mirror() { theMomentum.setZ(-theMomentum.z()); } /** * Do Lorentz transformations on this particle and its decendants. */ void deepTransform(const LorentzRotation & r); /** * Do Lorentz transformations on this particle and its * decendants. \a bx, \a by and \a bz are the boost vector * components. */ void deepBoost(double bx, double by, double bz) { deepTransform(LorentzRotation(Boost(bx, by, bz))); } /** * Do Lorentz transformations on this particle and its * decendants. \a b is the boost vector. */ void deepBoost(const Boost & b) { deepTransform(LorentzRotation(b)); } /** * Rotate this particle and its decendants around the x-axis. */ void deepRotateX(double a); /** * Rotate this particle and its decendants around the y-axis. */ void deepRotateY(double a); /** * Rotate this particle and its decendants around the z-axis. */ void deepRotateZ(double a); /** * Rotate this particle and its decendants around the given \a axis. */ void deepRotate(double a, const Axis & axis); //@} /** @name Functions controlling possible mass/momentum inconsistencies. */ //@{ /** * Return the relative inconsistency in the mass component. */ double massError() const { return theMomentum.massError(); } /** * Return the relative inconsistency in the energy component. */ double energyError() const { return theMomentum.energyError(); } /** * Return the relative inconsistency in the spatial components. */ double rhoError() const { return theMomentum.rhoError(); } /** * Rescale energy, so that the invariant length/mass of the * LorentzVector agrees with the current one. */ void rescaleEnergy() { theMomentum.rescaleEnergy(); } /** * Rescale spatial component, so that the invariant length/mass of * the LorentzVector agrees with the current one. */ void rescaleRho() { theMomentum.rescaleRho(); } /** * Set the invariant length/mass member, so that it agrees with the * invariant length/mass of the LorentzVector. */ void rescaleMass() { theMomentum.rescaleMass(); } //@} /** @name Acces incormation about colour connections */ //@{ /** * True if this particle has colour information. To determine if * this particle is actually coloured, the coloured(), hasColour() or * hasAntiColour() methods should be used instead. */ bool hasColourInfo() const { return hasRep() && rep().theColourInfo; } /** * Return the colour lines to which this particles anti-colour is * connected. */ tColinePtr antiColourLine() const { return hasColourInfo() ? colourInfo()->antiColourLine() : tColinePtr(); } /** * Return the colour lines to which this particles (\a anti-)colour * is connected. */ tColinePtr colourLine(bool anti = false) const { if ( anti ) return antiColourLine(); return hasColourInfo() ? colourInfo()->colourLine() : tColinePtr(); } /** * Return true if the particle is connected to the given (\a anti-) * colour \a line. */ bool hasColourLine(tcColinePtr line, bool anti = false) const { return hasColourInfo() ? colourInfo()->hasColourLine(line, anti) : false; } /** * Return true if the particle is connected to the given anti-colour * \a line. */ bool hasAntiColourLine(tcColinePtr line) const { return hasColourLine(line, true); } /** * True if this particle type is not a colour singlet. */ bool coloured() const { return data().coloured(); } /** * True if this particle type carries (\a anti-)colour. */ bool hasColour(bool anti = false) const { return data().hasColour(anti); } /** * True if this particle type carries anti-colour. */ bool hasAntiColour() const { return data().hasAntiColour(); } /** * Get the ColourBase object. */ tcCBPtr colourInfo() const { return hasRep() ? rep().theColourInfo : CBPtr(); } /** * Get the ColourBase object. */ tCBPtr colourInfo() { if ( !rep().theColourInfo ) { switch(theData->iColour()) { case PDT::Colour6: case PDT::Colour6bar: rep().theColourInfo = new_ptr(MultiColour()); break; default: rep().theColourInfo = new_ptr(ColourBase()); } } return rep().theColourInfo; } /** * Set the ColourBase object. */ void colourInfo(tCBPtr c) { rep().theColourInfo = c; } /** * Get a pointer to the colour neighbor. Returns a particle in the * range \a first to \a last which colour is connected to the same * line as this particles anti-colour. If \a anti is true return * antiColourNeighbour(). */ template typename std::iterator_traits::value_type colourNeighbour(Iterator first, Iterator last, bool anti = false) const; /** * Get a pointer to the anti-colour neighbor. Returns a particle in * the range \a first to \a last which anti-colour is * connected to the same line as this particles colour. */ template typename std::iterator_traits::value_type antiColourNeighbour(Iterator first, Iterator last) const { return colourNeighbour(first, last, true); } /** * Set the colour neighbor. Connects the given particles colour to * the same colour line as this particles anti-colour. If \a anti is * true call antiColourNeighbour(tPPtr). */ void colourNeighbour(tPPtr, bool anti = false); /** * Set the anti-colour neighbor. Connects the given particles * anti-colour to the same colour line as this particles colour. */ void antiColourNeighbour(tPPtr p) { colourNeighbour(p, true); } /** * Connect colour. Create a colour line connecting to it this * particles colour and the given particles anti-colour. */ void antiColourConnect(tPPtr neighbour) { colourConnect(neighbour, true); } /** * Connect colour. Create a colour line connecting to it this * particles anti-colour and the given particles colour. If \a anti * is true call antiColourConnect(tPPtr). */ void colourConnect(tPPtr neighbour, bool anti = false) { colourNeighbour(neighbour, anti); } /** * Incoming colour. Return the parent particle which colour is * connected to the same colour line as this particle. If \a anti is * true return incomingAntiColour(). */ tPPtr incomingColour(bool anti = false) const; /** * Incoming anti-colour. Return the parent particle which * anti-colour is connected to the same colour line as this * particle. */ tPPtr incomingAntiColour() const { return incomingColour(true); } /** * Set incoming colour. Connect this particles colour to the same * colour line as the given particle. If \a anti * is true call incomingAntiColour(tPPtr). */ void incomingColour(tPPtr p, bool anti = false) { p->outgoingColour(this, anti); } /** * Set incoming anti-colour. Connect this particles anti colour to * the same colour line as the given particle. */ void incomingAntiColour(tPPtr p) { p->outgoingColour(this, true); } /** * Outgoing colour. Return the daughter particle which colour is * connected to the same colour line as this particle. If \a anti is * true return outgoingAntiColour(). */ tPPtr outgoingColour(bool anti = false) const; /** * Outgoing anti-colour. Return the daughter particle which * anti-colour is connected to the same colour line as this * particle. */ tPPtr outgoingAntiColour() const { return outgoingColour(true); } /** * Set outgoing colour. Connect this particles colour to the same * colour line as the given particle. If \a anti * is true call outgoingAntiColour(tPPtr). */ void outgoingColour(tPPtr, bool anti = false); /** * Set outgoing anti-colour. Connect this particles anti-colour to * the same colour line as the given particle. */ void outgoingAntiColour(tPPtr p) { outgoingColour(p, true); } /** * Specify colour flow. Calls outgoingColour(tPPtr,bool). */ void colourFlow(tPPtr child, bool anti = false) { outgoingColour(child, anti); } /** * Specify anticolour flow. Calls outgoingAntiColour(tPPtr,bool). */ void antiColourFlow(tPPtr child) { colourFlow(child, true); } /** * Remove all colour information; */ void resetColour() { if ( hasColourInfo() ) rep().theColourInfo = CBPtr(); } //@} /** @name Functions to access spin. */ //@{ /** * Return the Spin object. */ tcSpinPtr spinInfo() const { return hasRep() ? rep().theSpinInfo : SpinPtr(); } /** * Return the Spin object. */ tSpinPtr spinInfo() { return hasRep() ? rep().theSpinInfo : SpinPtr(); } /** * Set the Spin object. */ void spinInfo(tSpinPtr s) { rep().theSpinInfo = s; } //@} /** @name Accessing user-defined information. */ //@{ /** * Access user-defined information as a vector of EventInfoBase pointers. */ const EIVector & getInfo() const { static const EIVector null; return hasRep() ? rep().theExtraInfo : null; } /** * Access user-defined information as a vector of EventInfoBase pointers. */ EIVector & getInfo() { return rep().theExtraInfo; } //@} public: /** @name Accessing user-defined information. */ //@{ /** * True if this particle has instantiated the object with * information other than type and momentum. */ bool hasRep() const { return theRep; } /** * If this particle has only a type and momentum, instantiate the * rest of the information. */ void initFull(); //@} public: /** @name Input and output functions. */ //@{ /** * Standard function for writing to a persistent stream. */ void persistentOutput(PersistentOStream &) const; /** * Standard function for reading from a persistent stream. */ void persistentInput(PersistentIStream &, int); //@} /** * Print particle info to a stream \a os. The \a step is used to * access information about colour neighbors and other struff. */ ostream & print(ostream & os, tcStepPtr step = tcStepPtr()) const; /** * Print a range of particles. */ template static void PrintParticles(ostream & os, Iterator first, Iterator last, tcStepPtr step = tcStepPtr()); /** * Print a container of particles. */ template static inline void PrintParticles(ostream & os, const Cont & c, tcStepPtr step = tcStepPtr()) { PrintParticles(os, c.begin(), c.end(), step); } /** * Standard Init function. @see Base::Init(). */ static void Init(); /** * Specify how to print particles. The format string is analogous to * the one used by eg. the unix 'date' command as described above. */ static string outputFormat; private: /** * Standard clone function. */ virtual PPtr clone() const; /** * Rebind to cloned objects. When an Event is cloned, a shallow * copy is done first, then all Particles etc, are * cloned, and finally this method is used to see to that the * pointers in the cloned Particle points to the cloned objects. */ virtual void rebind(const EventTranslationMap &); /** * Set the order-number for this particle in the current event. */ void number(int n) { rep().theNumber = n; } /** * Remove the given particle from the list of children. */ void removeChild(tPPtr c) { if ( hasRep() ) rep().theChildren.erase(remove(rep().theChildren.begin(), rep().theChildren.end(), c), rep().theChildren.end()); } /** * Remove the given particle from the list of parents. */ void removeParent(tPPtr p) { if ( hasRep() ) rep().theParents.erase(remove(rep().theParents.begin(), rep().theParents.end(), p), rep().theParents.end()); } /** * Set the mass of this particle. */ void mass(Energy m) { theMomentum.setMass(m); } /** * Set the invaiant life time of this particle. */ void lifeTime(Length t) { rep().theLifeLength.setTau(t); } /** * Return a reference to the bulk information of this particle. if * no ParticleRep object exists, one is created. */ ParticleRep & rep() { if ( !hasRep() ) initFull(); return *theRep; } /** * Return a reference to the bulk information of this particle. if * no ParticleRep object exists, we return the default values. */ const ParticleRep & rep() const { static const ParticleRep null; return hasRep() ? *theRep : null; } /** * The pointer to the ParticleData object */ cEventPDPtr theData; /** * The momentum. */ Lorentz5Momentum theMomentum; /** * The rest of the information in this particle is only instantiated * if needed. */ ParticleRep * theRep; + /** + * The status code of the particle + */ + int theStatus; + public: /** * This class is used internally in the Particle class to represent * information besides momentum and type. A corresponding object * will only be instantiated if needed to save memory and time when * temporarily creating particles. */ struct ParticleRep { /** * Default constructor. */ ParticleRep() : theScale(-1.0*GeV2), theVetoScale(-1.0*GeV2), theNumber(0) {} /** * Copy constructor. */ ParticleRep(const ParticleRep &); /** * The pointers to the parents. */ tParticleVector theParents; /** * The pointers to the children. */ ParticleVector theChildren; /** * The pointer to the previous instance. */ tPPtr thePrevious; /** * The pointer to the next instance. */ PPtr theNext; /** * If this particle has decayed this is the pointer to the * corresponding decay mode. */ tDMPtr theDecayMode; /** * The pointer to the first step where this particle occurred. */ tStepPtr theBirthStep; /** * The creation point. */ LorentzPoint theVertex; /** * The life time/length. */ Lorentz5Distance theLifeLength; /** * the resolution scale. */ Energy2 theScale; /** * the veto scale. */ Energy2 theVetoScale; /** * The order-number for this particle in the current event. */ int theNumber; /** * A pointer to the colour information object. */ CBPtr theColourInfo; /** * Spin information */ SpinPtr theSpinInfo; /** * Additional used-defined information. */ EIVector theExtraInfo; }; public: /** * Print out debugging information for this object on std::cerr. To * be called from within a debugger via the debug() function. */ virtual void debugme() const; protected: /** * Private default constructor must only be used by the * PersistentIStream class via the ClassTraits class. */ - Particle() : theRep(0) {} + Particle() : theRep(0), theStatus(0) {} /** * The ClassTraits class must be a friend to be able to * use the private default constructor. */ friend struct ClassTraits; private: /** * Private and non-existent assignment. */ Particle & operator=(const Particle &) = delete; /** * Describe concrete class with persistent data. */ static ClassDescription initParticle; }; /** * Write a Particle object to a stream. */ ostream & operator<<(ostream &, const Particle &); /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base class of Particle. */ template <> struct BaseClassTrait: public ClassTraitsType { /** Typedef of the first base class of Collision. */ typedef EventRecordBase NthBase; }; /** This template specialization informs ThePEG about the name of * the Particle class and how to create it. */ template <> struct ClassTraits: public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "ThePEG::Particle"; } /** Create a Particle object. */ static TPtr create() { return TPtr::Create(Particle()); } }; /** @endcond */ } #ifndef ThePEG_TEMPLATES_IN_CC_FILE #include "Particle.tcc" #endif #endif /* ThePEG_Particle_H */ diff --git a/PDF/PartonExtractor.cc b/PDF/PartonExtractor.cc --- a/PDF/PartonExtractor.cc +++ b/PDF/PartonExtractor.cc @@ -1,653 +1,653 @@ // -*- C++ -*- // // PartonExtractor.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the PartonExtractor class. // #include "PartonExtractor.h" #include "ThePEG/Handlers/XComb.h" #include "ThePEG/Cuts/Cuts.h" #include "ThePEG/PDF/NoPDF.h" #include "ThePEG/PDF/RemnantHandler.h" #include "ThePEG/PDF/BeamParticleData.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/EventRecord/Step.h" #include "ThePEG/EventRecord/SubProcess.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" #include "ThePEG/Utilities/UtilityBase.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/PDT/EnumParticles.h" using namespace ThePEG; PartonExtractor::PartonExtractor() : theMaxTries(100), flatSHatY(false) {} PartonExtractor::~PartonExtractor() {} IBPtr PartonExtractor::clone() const { return new_ptr(*this); } IBPtr PartonExtractor::fullclone() const { return new_ptr(*this); } PartonPairVec PartonExtractor:: getPartons(Energy maxEnergy, const cPDPair & incoming, const Cuts & kc) const { PartonPairVec result; PartonVector first; PDFCuts cuts1(kc, true, maxEnergy); PBPtr p1 = new_ptr(PartonBin(PDPtr(), PBPtr(), incoming.first, PDFPtr(), cuts1)); addPartons(p1, cuts1, theFirstPDF, first); PartonVector second; PDFCuts cuts2(kc, false, maxEnergy); PBPtr p2 = new_ptr(PartonBin(PDPtr(), PBPtr(), incoming.second, PDFPtr(), cuts2)); addPartons(p2, cuts2, theSecondPDF, second); for ( PartonVector::iterator it1 = first.begin(); it1 != first.end(); ++it1 ) for ( PartonVector::iterator it2 = second.begin(); it2 != second.end(); ++it2 ) result.push_back(PBPair(*it1, *it2)); // We add the original parton bins as well to avoid them being // deleted. result.push_back(PBPair(p1, p2)); return result; } void PartonExtractor:: addPartons(tPBPtr incoming, const PDFCuts & cuts, tcPDFPtr pdf, PartonVector & pbins) const { if(!pdf) pdf = getPDF(incoming->parton()); if ( dynamic_ptr_cast::tcp>(pdf) || incoming->parton() == incoming->particle() ) { pbins.push_back(incoming); return; } cPDVector partons = pdf->partons(incoming->parton()); for ( int i = 0, N = partons.size(); i < N; ++i ) { PBPtr pb = new_ptr(PartonBin(incoming->parton(), incoming, partons[i], pdf, cuts)); incoming->addOutgoing(pb); addPartons(pb, cuts, PDFPtr(), pbins); } } tcPDFPtr PartonExtractor::getPDF(tcPDPtr particle) const { for ( vector::const_iterator it = theSpecialDensities.begin(); it != theSpecialDensities.end(); ++it ) if ( (**it).canHandle(particle) ) return *it; Ptr::tcp p = dynamic_ptr_cast::tcp>(particle); if ( !p || !p->pdf() ) return noPDF(); return p->pdf(); } void PartonExtractor::select(tXCombPtr newXComb) { theLastXComb = newXComb; } tPBIPtr PartonExtractor::partonBinInstance(tcPPtr p) const { PartonBinInstanceMap::const_iterator it = partonBinInstances().find(p); return it == partonBinInstances().end()? PBIPtr(): it->second; } void PartonExtractor:: colourConnect(tPPtr particle, tPPtr parton, const tPVector & remnants) const { // Sorry cannot handle coloured resolved particles. if ( particle->coloured() ) throw RemColException(*this); // First connect the loose colour line from the extacted parton. if ( parton->hasColour() ) findConnect(parton->colourLine(), parton, true, remnants.rbegin(), remnants.rend()); // First connect the loose anti-colour line from the extacted parton. if ( parton->hasAntiColour() ) findConnect(parton->antiColourLine(), parton, false, remnants.begin(), remnants.end()); // Go through the rest of the remnants and create new colour lines // if needed. Go through it forwards and backwards to catch possible // inconsistencies. for ( tPVector::const_iterator it = remnants.begin(); it != remnants.end(); ++it ) { if ( (**it).hasAntiColour() && !(**it).antiColourLine() ) findConnect(ColourLine::create(*it, true), *it, false, it + 1, remnants.end()); } for ( tPVector::const_reverse_iterator it = remnants.rbegin(); it != remnants.rend(); ++it ) { if ( (**it).hasColour() && !(**it).colourLine() ) findConnect(ColourLine::create(*it), *it, true, it + 1, remnants.rend()); } } Energy2 PartonExtractor::newScale() { return lastScale(); } pair PartonExtractor::nDims(const PBPair & pbins) { // if photon from a lepton or proton generate scale bool genscale[2]={false,false}; for(unsigned int ix=0;ix<2;++ix) { PBPtr bin = ix==0 ? pbins.first : pbins.second; if (!bin || !bin->particle() || !bin->parton()) continue; if(bin->pdf()->partons(bin->particle()).size()==1 && bin->particle()->id()!=bin->parton()->id()) genscale[ix]=true; } return make_pair(pbins.first ->nDim(genscale[0]), pbins.second->nDim(genscale[1])); } void PartonExtractor::prepare(const PBIPair & pbins) { partonBinInstances().clear(); pbins.first->prepare(); pbins.second->prepare(); } void PartonExtractor::updatePartonBinInstances(const PBIPair & pbins) { partonBinInstances().clear(); tPBIPtr current = pbins.first; while ( current->incoming() ) { partonBinInstances()[current->parton()] = current; current = current->incoming(); } current = pbins.second; while ( current->incoming() ) { partonBinInstances()[current->parton()] = current; current = current->incoming(); } } bool PartonExtractor:: generateL(const PBIPair & pbins, const double * r1, const double * r2) { Direction<0> dir(true); generateL(*pbins.first, r1); dir.reverse(); generateL(*pbins.second, r2); if ( !flatSHatY || pbins.first->hasPoleIn1() || pbins.second->hasPoleIn1() ) return true; Energy2 shmax = lastCuts().sHatMax(); Energy2 shmin = lastCuts().sHatMin(); Energy2 sh = shmin*pow(shmax/shmin, *r1); double ymax = lastCuts().yHatMax(); double ymin = lastCuts().yHatMin(); double km = log(shmax/shmin); ymax = min(ymax, log(lastCuts().x1Max()*sqrt(lastS()/sh))); ymin = max(ymin, -log(lastCuts().x2Max()*sqrt(lastS()/sh))); double y = ymin + (*r2)*(ymax - ymin); double l1 = 0.5*log(lastS()/sh) - y; double l2 = 0.5*log(lastS()/sh) + y; pbins.first->li(l1 - pbins.first->l() + pbins.first->li()); pbins.first->l(l1); pbins.first->jacobian(km*(ymax - ymin)); pbins.second->li(l2 - pbins.second->l() + pbins.second->li()); pbins.second->l(l2); pbins.second->jacobian(1.0); return ( pbins.first->li() >= 0.0 && pbins.second->li() >= 0.0 ); } Energy2 PartonExtractor:: generateSHat(Energy2, const PBIPair & pbins, const double * r1, const double * r2, bool haveMEPartons) { Direction<0> dir(true); if(pbins.first->bin()->pdfDim()<=1) pbins.first->scale(-lastScale()); if ( !generate(*pbins.first, r1, lastSHat(), pbins.first->getFirst()->parton()->momentum(), haveMEPartons) ) return -1.0*GeV2; dir.reverse(); if(pbins.second->bin()->pdfDim()<=1) pbins.second->scale(-lastScale()); if ( !generate(*pbins.second, r2, lastSHat(), pbins.second->getFirst()->parton()->momentum(), haveMEPartons) ) return -1.0*GeV2; return (pbins.first->parton()->momentum() + pbins.second->parton()->momentum()).m2(); } void PartonExtractor:: generateL(PartonBinInstance & pb, const double * r) { if ( !pb.incoming() ) return; pb.parton(pb.partonData()->produceParticle(Lorentz5Momentum())); generateL(*pb.incoming(), r + pb.bin()->pdfDim() + pb.bin()->remDim()); pb.particle(pb.incoming()->parton()); if ( pb.li() >= 0 ) return; double jac = 1.0; if ( pb.bin()->pdfDim() ) pb.li(pb.pdf()->flattenL(pb.particleData(), pb.partonData(), pb.bin()->cuts(), *r++, jac)); pb.scale(-1.0*GeV2); if ( pb.bin()->pdfDim() > 1 ) pb.scale(pb.pdf()->flattenScale(pb.particleData(), pb.partonData(), pb.bin()->cuts(), pb.li(), *r++, jac) *pb.bin()->cuts().scaleMaxL(pb.li())); pb.jacobian(jac); pb.l(pb.incoming()->l() + pb.li()); } bool PartonExtractor:: generate(PartonBinInstance & pb, const double * r, Energy2 shat, const Lorentz5Momentum & first, bool haveMEPartons) { if ( !pb.incoming() ) return true; if ( !generate(*pb.incoming(), r + pb.bin()->pdfDim() + pb.bin()->remDim(), shat/pb.xi(), first) ) return false; pb.remnantWeight(1.0); - pb.parton()->setMomentum + pb.parton()->set5Momentum (pb.remnantHandler()->generate(pb, r + pb.bin()->pdfDim(), pb.scale(), shat, pb.particle()->momentum(),haveMEPartons)); if ( pb.remnantWeight() <= 0.0 ) return false; partonBinInstances()[pb.parton()] = &pb; return true; } void PartonExtractor:: constructRemnants(const PBIPair & pbins, tSubProPtr sub, tStepPtr step) const { partonBinInstances().clear(); LorentzMomentum k1 = pbins.first->parton()->momentum(); LorentzMomentum k2 = pbins.second->parton()->momentum(); LorentzMomentum Ph = k1 + k2; LorentzMomentum Phold = Ph; LorentzRotation Rh = Utilities::getBoostToCM(make_pair(k1, k2)); bool pickside = rndbool(); if ( pickside && pbins.first->incoming() ) { Direction<0> dir(true); constructRemnants(*pbins.first, Ph, k2); construct(*pbins.first, step, false); } if ( pbins.second->incoming() ) { Direction<0> dir(false); constructRemnants(*pbins.second, Ph, pbins.first->parton()->momentum()); construct(*pbins.second, step, false); } if ( (!pickside) && pbins.first->incoming() ) { Direction<0> dir(true); constructRemnants(*pbins.first, Ph, pbins.second->parton()->momentum()); construct(*pbins.first, step, false); } // LorentzRotation rot = Utilities::transformToMomentum(Phold, Ph); k1 = pbins.first->parton()->momentum(); k2 = pbins.second->parton()->momentum(); LorentzRotation rot = Utilities::getBoostFromCM(make_pair(k1, k2))*Rh; Utilities::transform(sub->outgoing(), rot); Utilities::transform(sub->intermediates(), rot); Ph = k1 + k2; if ( abs(Ph.m2() - Phold.m2())/Phold.m2() > 0.000001 ) cerr << Ph.m2()/GeV2 << " was (" << Phold.m2()/GeV2 << ")" << endl; } void PartonExtractor:: constructRemnants(PartonBinInstance & pb, LorentzMomentum & Ph, const LorentzMomentum & k) const { LorentzMomentum P = pb.particle()->momentum(); DVector r = UseRandom::rndvec(pb.bin()->remDim()); if ( r.empty() ) r.push_back(0.0); pb.parton()->setMomentum(pb.remnantHandler()-> generate(pb, &r[0], pb.scale(), Ph.m2(), P)); if ( pb.remnantWeight() <= 0.0 ) throw Veto(); pb.remnantHandler()->boostRemnants(pb); LorentzMomentum Pr = Utilities::sumMomentum(pb.remnants()); transformRemnants(Ph, Pr, k, pb.particle()->momentum()); pb.parton()->setMomentum(pb.particle()->momentum() - Pr); try { Utilities::setMomentum(pb.remnants().begin(), pb.remnants().end(), static_cast(Pr)); } catch ( ThePEG::Exception & e) { throw e; } catch ( ThePEG::Veto ) { throw; } catch ( std::exception & e ) { throw Exception() << "Caught non-ThePEG exception " << e.what() << "in " << "PartonExtractor::constructRemnants" << Exception::eventerror; } partonBinInstances()[pb.parton()] = &pb; if ( !pb.incoming()->incoming() ) return; // We get here if we need to construct remnants recursively. LorentzMomentum Phnew = Ph + Pr; constructRemnants(*pb.incoming(), Phnew, k); LorentzRotation rot = Utilities::transformToMomentum(Ph + Pr, Phnew); Utilities::transform(pb.remnants(), rot); Ph.transform(rot); } LorentzRotation PartonExtractor:: boostRemnants(PBIPair & bins, LorentzMomentum k1, LorentzMomentum k2, bool side1, bool side2) const { if ( !side1 && !side2 ) return LorentzRotation(); LorentzMomentum P1 = bins.first? LorentzMomentum(bins.first->parton()->momentum()): k1; LorentzMomentum Pr1; if ( side1 ) { P1 = bins.first->particle()->momentum(); Pr1 = Utilities::sumMomentum(bins.first->remnants()); } LorentzMomentum P2 = bins.second? LorentzMomentum(bins.second->parton()->momentum()): k2; LorentzMomentum Pr2; if ( side2 ) { P2 = bins.second->particle()->momentum(); Pr2 = Utilities::sumMomentum(bins.second->remnants()); } LorentzRotation Rh = Utilities::getBoostToCM(make_pair(k1, k2)); LorentzMomentum Ph = k1 + k2; // LorentzMomentum Phold = Ph; bool otherside = rndbool(); if ( otherside && side2 ){ Direction<0> dir(false); transformRemnants(Ph, Pr2, k1, P2); k2 = P2 - Pr2; } if ( side1 ){ Direction<0> dir(true); transformRemnants(Ph, Pr1, k2, P1); k1 = P1 - Pr1; } if ( side2 && !otherside ) { Direction<0> dir(false); transformRemnants(Ph, Pr2, k1, P2); k2 = P2 - Pr2; } if ( bins.first ) { if ( bins.first->remnants().size() == 1 ) bins.first->remnants()[0]->setMomentum(Pr1); else Utilities::setMomentum(bins.first->remnants().begin(), bins.first->remnants().end(), static_cast(Pr1)); bins.first->parton()->setMomentum(k1); } if ( bins.second ) { if ( bins.second->remnants().size() == 1 ) bins.second->remnants()[0]->setMomentum(Pr2); else Utilities::setMomentum(bins.second->remnants().begin(), bins.second->remnants().end(), static_cast(Pr2)); bins.second->parton()->setMomentum(k2); } Rh.transform(Utilities::getBoostFromCM(make_pair(k1, k2))); // LorentzMomentum phh = Rh*Phold; return Rh; // return Utilities::transformToMomentum(Phold, Ph); } void PartonExtractor:: transformRemnants(LorentzMomentum & Ph, LorentzMomentum & Pr, const LorentzMomentum & k, const LorentzMomentum & P) const { // don't do this for very soft remnants, as // we may run into numerical troubles; threshold // needs to become a parameter at some point if ( Pr.vect().mag2()/k.vect().mag2() < 1e-10 && sqr(Pr.e()/k.e()) < 1e-10 ) return; TransverseMomentum pt = Pr; try { if ( Direction<0>::pos() ) SimplePhaseSpace::CMS(Pr, Ph, (P + k).m2(), 1.0, 0.0); else SimplePhaseSpace::CMS(Ph, Pr, (k + P).m2(), 1.0, 0.0); LorentzRotation rpt; if ( sqr(Pr.z()) > ZERO ) rpt.rotateY(asin(pt.pt()/Pr.z())); rpt.rotateZ(pt.phi()); rpt = Direction<0>::pos()? Utilities::getBoostFromCM(make_pair(P, k))*rpt: Utilities::getBoostFromCM(make_pair(k, P))*rpt; Ph.transform(rpt); Pr.transform(rpt); } catch ( ImpossibleKinematics & e ) {} } double PartonExtractor::fullFn(const PBIPair & pbins, Energy2 scale, pair noLastPDF) { if(pbins.first->bin()->pdfDim()<=1) pbins.first->scale(scale); if(pbins.second->bin()->pdfDim()<=1) pbins.second->scale(scale); return fullFn(*pbins.first,noLastPDF.first)*fullFn(*pbins.second,noLastPDF.second); } double PartonExtractor::fullFn(const PartonBinInstance & pb, bool noLastPDF) { if ( !pb.incoming() ) return 1.0; if (noLastPDF) return fullFn(*pb.incoming(),false) * pb.jacobian() * pb.remnantWeight() * exp(-pb.li()); return fullFn(*pb.incoming(),false) * pb.jacobian() * pb.remnantWeight() * pb.pdf()->xfl(pb.particleData(), pb.partonData(), pb.scale(), pb.li(), pb.incoming()->scale()); } void PartonExtractor:: construct(const PBIPair & pbins, tStepPtr step) const { // if a long chain we need to break some mother/child relationships if(pbins.first->incoming()) { if(pbins.first->incoming()->incoming()) { if(!pbins.first->parton()->parents().empty()) { tParticleVector parents=pbins.first->parton()->parents(); tPPtr parton = pbins.first->parton(); for(unsigned int ix=0;ixabandonChild(parton); } } } if(pbins.second->incoming()) { if(pbins.second->incoming()->incoming()) { if(!pbins.second->parton()->parents().empty()) { tParticleVector parents=pbins.second->parton()->parents(); tPPtr parton = pbins.second->parton(); for(unsigned int ix=0;ixabandonChild(parton); } } } Direction<0> dir(true); construct(*pbins.first, step); dir.reverse(); construct(*pbins.second, step); } void PartonExtractor:: construct(PartonBinInstance & pb, tStepPtr step, bool boost) const { if ( !pb.incoming() ) return; if ( boost ) pb.remnantHandler()->boostRemnants(pb); if ( pb.incoming()->incoming() ) { step->insertIntermediate(pb.particle(),pb.incoming()->particle(),pb.parton()); } tPVector rem(pb.remnants().begin(), pb.remnants().end()); if ( !step->addDecayProduct(pb.particle(), rem.begin(), rem.end(), false) ) {} colourConnect(pb.particle(), pb.parton(), rem); construct(*pb.incoming(), step); } PBIPair PartonExtractor::newRemnants(tPPair oldp, tPPair newp, tStepPtr step) { PBIPair pb; Direction<0> dir(true); pb.first = newRemnants(partonBinInstance(oldp.first), newp.first, newp.second->momentum()); dir.reverse(); pb.second = newRemnants(partonBinInstance(oldp.second), newp.second, newp.first->momentum()); addNewRemnants(partonBinInstance(oldp.first), pb.first, step); addNewRemnants(partonBinInstance(oldp.second), pb.second, step); return pb; } PBIPtr PartonExtractor:: newRemnants(tPBIPtr oldpb, tPPtr newp, const LorentzMomentum & k) { if ( ! oldpb || !oldpb->incoming() ) return oldpb; Energy2 shat = (k + newp->momentum()).m2(); // Loop over all possible PartonBin sisters to find the one // corresponding to the newly extracted parton. const PartonBin::PBVector & sisters = oldpb->incoming()->bin()->outgoing(); for ( int i = 0, N = sisters.size(); i < N; ++i ) if ( sisters[i]->parton() == newp->dataPtr() ) { // Setup necessary info in new PartonBinInstance object. PBIPtr newpb = new_ptr(PartonBinInstance(sisters[i], oldpb->incoming())); newpb->particle(oldpb->particle()); newpb->parton(newp); newpb->li(log(oldpb->particle()->momentum().dirPlus()/ newp->momentum().dirPlus())); newpb->l(oldpb->l() - oldpb->li() + newpb->li()); Energy2 sc = -newp->scale(); newpb->scale(newp->scale()); if ( oldpb->incoming()->incoming() ) sc = -newpb->particle()->momentum().m2(); // Now we can construct the new remnants. newpb->remnantWeight(1.0); if ( !newpb->remnantHandler()-> recreateRemnants(*newpb, oldpb->parton(), newp, newpb->li(), sc, shat, newpb->particle()->momentum()) ) throw Veto(); if ( newpb->remnantWeight() <= 0.0 ) throw Veto(); return newpb; } throw Veto(); } void PartonExtractor:: addNewRemnants(tPBIPtr oldpb, tPBIPtr newpb, tStepPtr step) { if ( oldpb == newpb ) return; if ( oldpb->parton() != newpb->parton() ) { step->removeDecayProduct(newpb->particle(), oldpb->parton()); if ( !step->addDecayProduct(newpb->particle(), newpb->parton()) ) throw Veto(); } tPVector rem(newpb->remnants().begin(), newpb->remnants().end()); colourConnect(newpb->particle(), newpb->parton(), rem); partonBinInstances()[newpb->parton()] = newpb; if ( !step->addDecayProduct(oldpb->remnants().begin(), oldpb->remnants().end(), rem.begin(), rem.end()) ) throw Veto(); } void PartonExtractor::persistentOutput(PersistentOStream & os) const { os << theLastXComb << theSpecialDensities << theNoPDF << theMaxTries << flatSHatY << theFirstPDF << theSecondPDF; } void PartonExtractor::persistentInput(PersistentIStream & is, int) { is >> theLastXComb >> theSpecialDensities >> theNoPDF >> theMaxTries >> flatSHatY >> theFirstPDF >> theSecondPDF; } ClassDescription PartonExtractor::initPartonExtractor; void PartonExtractor::Init() { static ClassDocumentation documentation ("There is no documentation for the ThePEG::PartonExtractor class"); static RefVector interfaceSpecialDensities ("SpecialDensities", "A list of parton density objects to be used for incoming particles " "overriding possible densities given for particles of the " "BeamParticleData class.", &PartonExtractor::theSpecialDensities, 0, false, false, true, false); static Reference interfaceNoPDF ("NoPDF", "A fixed reference to a NoPDF object to be used for particles without " "substructure.", &PartonExtractor::theNoPDF, true, true, true, false); static Parameter interfaceMaxTries ("MaxTries", "The maximum number of attempts allowed when trying to generate " "remnants.", &PartonExtractor::theMaxTries, 100, 1, 1000, false, false, true); static Switch interfaceFlatSHatY ("FlatSHatY", "The possibility to override the l-generation in the PDFs and generate " "a flat distribution in \\f$\\log(\\hat{s})\\f$ and \\f$y\\f$. This only " "applies if the parton densities do not have poles in \\f$x=1\\f$.", &PartonExtractor::flatSHatY, false, false, false); static SwitchOption interfaceFlatSHatY0 (interfaceFlatSHatY, "Off", "Use the l-generation defined by the PDFs", false); static SwitchOption interfaceFlatSHatY1 (interfaceFlatSHatY, "On", "Generate flat rapidity and \\f$\\log(\\hat{s})\\f$", true); static SwitchOption interfaceFlatSHatNo (interfaceFlatSHatY, "No", "Use the l-generation defined by the PDFs", false); static SwitchOption interfaceFlatSHatYes (interfaceFlatSHatY, "Yes", "Generate flat rapidity and \\f$\\log(\\hat{s})\\f$", true); static Reference interfaceFirstPDF ("FirstPDF", "PDF to override the default PDF for the first beam particle", &PartonExtractor::theFirstPDF, false, false, true, true, false); static Reference interfaceSecondPDF ("SecondPDF", "PDF to override the default PDF for the second beam particle", &PartonExtractor::theSecondPDF, false, false, true, true, false); } RemColException::RemColException(const PartonExtractor & pe) { theMessage << "Parton extractor '" << pe.name() << "' failed to connect " << "the colours of the outgoing partons and the remnants."; severity(maybeabort); } void PartonExtractor::dofinish() { partonBinInstances().clear(); HandlerBase::dofinish(); } diff --git a/Utilities/SimplePhaseSpace.cc b/Utilities/SimplePhaseSpace.cc --- a/Utilities/SimplePhaseSpace.cc +++ b/Utilities/SimplePhaseSpace.cc @@ -1,112 +1,116 @@ // -*- C++ -*- // // SimplePhaseSpace.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #include "SimplePhaseSpace.h" #ifdef ThePEG_TEMPLATES_IN_CC_FILE #include "SimplePhaseSpace.tcc" #endif using namespace ThePEG; Energy SimplePhaseSpace::getMagnitude(Energy2 s, Energy m1, Energy m2) { const Energy2 eps = 10.0*s*Constants::epsilon; if ( m1 < ZERO && sqr(m1) < eps ) m1 = ZERO; if ( m2 < ZERO && sqr(m2) < eps ) m2 = ZERO; if ( m1 >= ZERO && m2 >= ZERO ) { Energy2 aa = s - sqr(m1+m2); if ( aa < ZERO && aa > -eps ) return ZERO; if ( aa < ZERO ) throw ImpossibleKinematics(); return 0.5*sqrt(aa*(s-sqr(m1-m2))/s); } - Energy2 m12 = m1 < ZERO? -sqr(m1): sqr(m1); - Energy2 m22 = m2 < ZERO? -sqr(m2): sqr(m2); - Energy2 r2 = 0.25*(sqr(m12) + sqr(m22 - s) -2.0*m12*(m22 + s))/s; - if ( r2 < ZERO || r2 + m12 < ZERO || r2 + m22 < ZERO ) - throw ImpossibleKinematics(); - return sqrt(r2); + if ( m1 >= ZERO && m2 < ZERO ) { + return sqrt(sqr(m2)+sqr(s-sqr(m1)-sqr(m2))/(4.*s)); + } + if ( m1 < ZERO && m2 >= ZERO ) { + return sqrt(sqr(m1)+sqr(s-sqr(m1)-sqr(m2))/(4.*s)); + } + if ( m1 < ZERO && m2 < ZERO ) { + return sqrt(sqr(m1)+sqr(s-sqr(m1)+sqr(m2))/(4.*s)); + } + return ZERO; } // Energy SimplePhaseSpace::checkMagnitude(Energy2 s, Energy m1, Energy m2) // { // if ( s < ZERO ) return -1.0*GeV; // const Energy2 eps = 10.0*s*Constants::epsilon; // if ( m1 < ZERO && sqr(m1) < eps ) m1 = ZERO; // if ( m2 < ZERO && sqr(m2) < eps ) m2 = ZERO; // if ( m1 >= ZERO && m2 >= ZERO ) { // Energy2 aa = s - sqr(m1+m2); // if ( aa < ZERO && aa > -eps ) return ZERO; // if ( aa < ZERO ) return -1.0*GeV; // return 0.5*sqrt(aa*(s-sqr(m1-m2))/s); // } // Energy2 m12 = m1 < ZERO? -sqr(m1): sqr(m1); // Energy2 m22 = m2 < ZERO? -sqr(m2): sqr(m2); // Energy2 r2 = 0.25*(sqr(m12) + sqr(m22 - s) -2.0*m12*(m22 + s))/s; // if ( r2 < ZERO || r2 + m12 < ZERO || r2 + m22 < ZERO ) // return -1.0*GeV; // return sqrt(r2); // } vector SimplePhaseSpace:: CMSn(Energy m0, const vector & m) { using Constants::pi; // Setup constants. int Np = m.size(); vector ret(Np); Energy summ = std::accumulate(m.begin(), m.end(), Energy()); if ( summ >= m0 ) throw ImpossibleKinematics(); while ( true ) { // First get an ordered list of random numbers. vector rndv(Np); rndv[0] = 1.0; rndv.back() = 0.0; for ( int i = 1; i < Np - 1; ++i ) rndv[i] = UseRandom::rnd(); std::sort(rndv.begin() + 1, rndv.end() - 1, std::greater()); // Now setup masses of subsystems. vector sm(Np); Energy tmass = m0 - summ; Energy tmp = summ; for ( int i = 0; i < Np; ++i ) { sm[i] = rndv[i]*tmass + tmp; tmp -= m[i]; } // Now the magnitude of all the momenta can be calculated. This // gives the weight. double weight = 1.0; vector p(Np); p[Np - 1] = getMagnitude(sqr(sm[Np - 2]), m[Np -2], sm[Np - 1]); for ( int i = Np - 2; i >= 0; --i ) weight *= (p[i] = getMagnitude(sqr(sm[i]), m[i], sm[i + 1]))/sm[i]; if ( weight > UseRandom::rnd() ) continue; // Now we just have to generate the angles. ret[Np - 1] = LorentzMomentum(ZERO, ZERO, ZERO, m[Np - 1]); for ( int i = Np - 2; i >= 0; --i ) { Momentum3 p3 = polar3Vector(p[i], 2.0*UseRandom::rnd() - 1.0, 2.0*pi*UseRandom::rnd()); ret[i] = LorentzMomentum(-p3, sqrt(sqr(p[i]) + sqr(m[i]))); if ( i == Np -2 ) { ret[Np - 1] = LorentzMomentum(p3, sqrt(sqr(m[Np - 1]) + p3.mag2())); } else { Boost bv = p3*(1.0/sqrt(sqr(p[i]) + sqr(sm[i + 1]))); if ( bv.mag2() >= 1.0 ) throw ImpossibleKinematics(); LorentzRotation r(bv); for ( int j = i + 1; j < Np; ++j ) ret[j]*=r.one(); } } return ret; } } diff --git a/Utilities/SimplePhaseSpace.tcc b/Utilities/SimplePhaseSpace.tcc --- a/Utilities/SimplePhaseSpace.tcc +++ b/Utilities/SimplePhaseSpace.tcc @@ -1,136 +1,141 @@ // -*- C++ -*- // // SimplePhaseSpace.tcc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined templated member // functions of the SimplePhaseSpace class. // namespace ThePEG { template void SimplePhaseSpace::CMS(PType & p1, PType & p2, Energy2 s) { typedef ParticleTraits Traits; - Energy z = getMagnitude(s, Traits::mass(p1), Traits::mass(p2)); - Traits::set3Momentum(p1, Momentum3(ZERO, ZERO, z)); - Traits::set3Momentum(p2, Momentum3(ZERO, ZERO, -z)); + Energy m1 = Traits::mass(p1); Energy m2 = Traits::mass(p2); + Energy z = getMagnitude(s, m1, m2); + Energy2 m12 = m1 >= ZERO ? sqr(m1) : -sqr(m1); + Energy2 m22 = m2 >= ZERO ? sqr(m2) : -sqr(m2); + Energy2 c1 = (s+m12-m22); + Energy2 c2 = (s-m12+m22); + Traits::set5Momentum(p1, Lorentz5Momentum(ZERO, ZERO, z, (c1 > ZERO ? 1. : -1.) * 0.5*sqrt(sqr(c1)/s), m1)); + Traits::set5Momentum(p2, Lorentz5Momentum(ZERO, ZERO, -z, (c2 > ZERO ? 1. : -1.) * 0.5*sqrt(sqr(c2)/s), m2)); } template void SimplePhaseSpace::CMS(Energy2 s, PType & p1, PType & p2) { CMS(p1, p2, s, 2.0*UseRandom::rnd() - 1.0, Constants::twopi*UseRandom::rnd()); } template void SimplePhaseSpace::CMS(PType & p1, PType & p2, Energy2 s, double cthe, double phi) { typedef ParticleTraits Traits; Energy r = getMagnitude(s, Traits::mass(p1), Traits::mass(p2)); double sthe = sqrt(1.0-sqr(cthe)); Momentum3 p(r*sthe*cos(phi), r*sthe*sin(phi), r*cthe); Traits::set3Momentum(p1, p); Traits::set3Momentum(p2, -p); } template void SimplePhaseSpace:: CMS(PType & p1, PType & p2, PType & p3, Energy2 s, double x1, double x3) { CMS(p1, p2, p3, s, x1, x3, Constants::twopi*UseRandom::rnd(), acos(2.0*UseRandom::rnd() - 1.0), Constants::twopi*UseRandom::rnd()); } template void SimplePhaseSpace:: CMS(PType & p1, PType & p2, PType & p3, Energy2 s, double x1, double x3, double phii, double the, double phi) { typedef ParticleTraits Traits; Energy Etot = sqrt(s); Energy m1 = Traits::mass(p1); Energy m2 = Traits::mass(p2); Energy m3 = Traits::mass(p3); Energy e1 = 0.5*x1*Etot; Energy e3 = 0.5*x3*Etot; Energy e2 = Etot - e1 - e3; if ( e1 < m1 || e2 < m2 || e3 < m3 ) throw ImpossibleKinematics(); Energy r1 = sqrt(sqr(e1)-sqr(m1)); Energy r2 = sqrt(sqr(e2)-sqr(m2)); Energy r3 = sqrt(sqr(e3)-sqr(m3)); Traits::set3Momentum(p1, Momentum3(ZERO, ZERO, r1)); double cthe2 = (sqr(r3)-sqr(r2)-sqr(r1))/(2.0*r2*r1); double cthe3 = (sqr(r2)-sqr(r3)-sqr(r1))/(2.0*r3*r1); if ( abs(cthe2) > 1.0 || abs(cthe3) > 1.0 ) throw ImpossibleKinematics(); double sthe2 = sqrt(1.0-sqr(cthe2)); Energy px = r2*sthe2*cos(phii); Energy py = r2*sthe2*sin(phii); Traits::set3Momentum(p2, Momentum3(px, py, r2*cthe2)); Traits::set3Momentum(p3, Momentum3(-px, -py, r3*cthe3)); if ( the == 0.0 && phi == 0.0 ) return; LorentzRotation r; r.rotateZ(phi); r.rotateX(the); Traits::transform(p1, r); Traits::transform(p2, r); Traits::transform(p3, r); } template void SimplePhaseSpace:: CMS(PType & p1, PType & p2, Energy2 s, Energy2 t, double phi, const PType & p0) { typedef ParticleTraits Traits; Energy r = getMagnitude(s, Traits::mass(p1), Traits::mass(p2)); Energy e = sqrt(sqr(r) + sqr(Traits::mass(p1))); Energy r0 = Traits::momentum(p0).rho(); Energy e0 = Traits::momentum(p0).e(); double cthe = (t + sqr(e - e0) + sqr(r) + sqr(r0))/(2.0*r*r0); if ( abs(cthe) > 1.0 ) throw ImpossibleKinematics(); double sthe = sqrt(1.0-sqr(cthe)); Momentum3 p(r*sthe*cos(phi), r*sthe*sin(phi), r*cthe); Traits::set3Momentum(p1, p); Traits::set3Momentum(p2, -p); if ( Traits::momentum(p0).perp2() > ZERO ) { LorentzRotation r; r.rotateX(Traits::momentum(p0).theta()); r.rotateZ(Traits::momentum(p0).phi()); Traits::transform(p1, r); Traits::transform(p2, r); } } template void SimplePhaseSpace::CMSn(Container & particles, Energy m0) { typedef typename Container::value_type PType; typedef typename Container::iterator Iterator; if ( particles.size() == 2 ) { Iterator it = particles.begin(); PType & p1 = *it++; PType & p2 = *it; CMS(sqr(m0), p1, p2); return; } typedef ParticleTraits Traits; vector masses(particles.size()); int j = 0; for ( Iterator i = particles.begin();i != particles.end(); ++i, ++j ) masses[j] = Traits::mass(*i); vector p = CMSn(m0, masses); j = 0; for ( Iterator i = particles.begin();i != particles.end(); ++i, ++j ) Traits::set5Momentum(*i, p[j]); } } diff --git a/Vectors/HepMCConverter.tcc b/Vectors/HepMCConverter.tcc --- a/Vectors/HepMCConverter.tcc +++ b/Vectors/HepMCConverter.tcc @@ -1,337 +1,337 @@ // -*- C++ -*- // // HepMCConverter.tcc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the HepMCConverter class. // #include "ThePEG/StandardModel/StandardModelBase.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/EventRecord/StandardSelectors.h" #include "ThePEG/EventRecord/Collision.h" #include "ThePEG/EventRecord/Step.h" #include "ThePEG/EventRecord/SubProcess.h" #include "ThePEG/Handlers/XComb.h" #include "ThePEG/Handlers/EventHandler.h" #include "ThePEG/PDF/PartonExtractor.h" #include "ThePEG/PDF/PDF.h" #include "ThePEG/PDT/StandardMatchers.h" #include "ThePEG/Utilities/Throw.h" namespace ThePEG { template typename HepMCConverter::GenEvent * HepMCConverter:: convert(const Event & ev, bool nocopies, Energy eunit, Length lunit) { HepMCConverter converter(ev, nocopies, eunit, lunit); return converter.geneve; } template void HepMCConverter:: convert(const Event & ev, GenEvent & gev, bool nocopies) { HepMCConverter converter(ev, gev, nocopies, Traits::momentumUnit(gev), Traits::lengthUnit(gev)); } template void HepMCConverter:: convert(const Event & ev, GenEvent & gev, bool nocopies, Energy eunit, Length lunit) { HepMCConverter converter(ev, gev, nocopies, eunit, lunit); } template HepMCConverter:: HepMCConverter(const Event & ev, bool nocopies, Energy eunit, Length lunit) : energyUnit(eunit), lengthUnit(lunit) { geneve = Traits::newEvent(ev.number(), ev.weight(), ev.optionalWeights()); init(ev, nocopies); } template HepMCConverter:: HepMCConverter(const Event & ev, GenEvent & gev, bool nocopies, Energy eunit, Length lunit) : energyUnit(eunit), lengthUnit(lunit) { geneve = &gev; Traits::resetEvent(geneve, ev.number(), ev.weight(), ev.optionalWeights()); init(ev, nocopies); } struct ParticleOrderNumberCmp { bool operator()(tcPPtr a, tcPPtr b) const { return a->number() < b->number(); } }; template void HepMCConverter::init(const Event & ev, bool nocopies) { if ( lengthUnit != millimeter && lengthUnit != centimeter ) Throw() << "Length unit used for HepMC::GenEvent was not MM nor CM." << Exception::runerror; if ( energyUnit != GeV && energyUnit != MeV ) Throw() << "Momentum unit used for HepMC::GenEvent was not GEV nor MEV." << Exception::runerror; Traits::setUnits(*geneve, energyUnit, lengthUnit); tcEHPtr eh; if ( ev.primaryCollision() && ( eh = dynamic_ptr_cast(ev.primaryCollision()->handler()) ) ) { // Get general event info if present. Traits::setScaleAndAlphas(*geneve, eh->lastScale(), eh->lastAlphaS(),eh->lastAlphaEM(), energyUnit); } // Extract all particles and order them. tcPVector all; ev.select(back_inserter(all), SelectAll()); stable_sort(all.begin(), all.end(), ParticleOrderNumberCmp()); vertices.reserve(all.size()*2); // Create GenParticle's and map them to the ThePEG particles. for ( int i = 0, N = all.size(); i < N; ++i ) { tcPPtr p = all[i]; if ( nocopies && p->next() ) continue; if ( pmap.find(p) != pmap.end() ) continue; pmap[p] = createParticle(p); if ( !p->children().empty() || p->next() ) { // If the particle has children it should have a decay vertex: vertices.push_back(Vertex()); decv[p] = &vertices.back(); vertices.back().in.insert(p); } if ( !p->parents().empty() || p->previous() || (p->children().empty() && !p->next()) ) { // If the particle has parents it should have a production // vertex. If neither parents or children it should still have a // dummy production vertex. vertices.push_back(Vertex()); prov[p] = &vertices.back(); vertices.back().out.insert(p); } } // Now go through the the particles again, and join the vertices. for ( int i = 0, N = all.size(); i < N; ++i ) { tcPPtr p = all[i]; if ( nocopies ) { if ( p->next() ) continue; for ( int i = 0, N = p->children().size(); i < N; ++i ) join(p, p->children()[i]->final()); tcPPtr pp = p; while ( pp->parents().empty() && pp->previous() ) pp = pp->previous(); for ( int i = 0, N = pp->parents().size(); i < N; ++i ) join(pp->parents()[i]->final(), p); } else { for ( int i = 0, N = p->children().size(); i < N; ++i ) join(p, p->children()[i]); if ( p->next() ) join(p, p->next()); for ( int i = 0, N = p->parents().size(); i < N; ++i ) join(p->parents()[i], p); if ( p->previous() ) join(p->previous(), p); } } // Time to create the GenVertex's for ( typename VertexMap::iterator it = prov.begin(); it != prov.end(); ++it ) if ( !member(vmap, it->second) ) vmap[it->second] = createVertex(it->second); for ( typename VertexMap::iterator it = decv.begin(); it != decv.end(); ++it ) if ( !member(vmap, it->second) ) vmap[it->second] = createVertex(it->second); // Now find the primary signal process vertex defined to be the // decay vertex of the first parton coming into the primary hard // sub-collision. tSubProPtr sub = ev.primarySubProcess(); if ( sub && sub->incoming().first ) { const Vertex * prim = decv[sub->incoming().first]; Traits::setSignalProcessVertex(*geneve, vmap[prim]); vmap.erase(prim); } // Then add the rest of the vertices. for ( typename GenVertexMap::iterator it = vmap.begin(); it != vmap.end(); ++it ) Traits::addVertex(*geneve, it->second); // and the incoming beam particles Traits::setBeamParticles(*geneve,pmap[ev.incoming().first], pmap[ev.incoming().second]); // and the PDF info setPdfInfo(ev); // and the cross section info Traits::setCrossSection(*geneve, eh->integratedXSec()/picobarn, eh->integratedXSecErr()/picobarn); for ( int i = 0, N = all.size(); i < N; ++i ) { tcPPtr p = all[i]; if ( pmap.find(p) == pmap.end() ) continue; GenParticlePtrT gp = pmap[p]; if ( p->hasColourInfo() ) { // Check if the particle is connected to colour lines, in which // case the lines are mapped to an integer and set in the // GenParticle's Flow info. tcColinePtr l; if ( (l = p->colourLine()) ) { if ( !member(flowmap, l) ) flowmap[l] = flowmap.size() + 500; Traits::setColourLine(*gp, 1, flowmap[l]); } if ( (l = p->antiColourLine()) ) { if ( !member(flowmap, l) ) flowmap[l] = flowmap.size() + 500; Traits::setColourLine(*gp, 2, flowmap[l]); } } if ( p->spinInfo() && p->spinInfo()->hasPolarization() ) { DPair pol = p->spinInfo()->polarization(); Traits::setPolarization(*gp, pol.first, pol.second); } } } template typename HepMCConverter::GenParticlePtrT HepMCConverter::createParticle(tcPPtr p) const { int status = 1; size_t nChildren = p->children().size(); if ( nChildren > 0 || p->next() ) status = 11; if ( nChildren > 1 ) { long id = p->data().id(); if ( BaryonMatcher::Check(id) || MesonMatcher::Check(id) || id == ParticleID::muminus || id == ParticleID::muplus || id == ParticleID::tauminus || id == ParticleID::tauplus ) { bool child = false; for(unsigned int ix=0;ixchildren()[ix]->id()==id) { child = true; break; } } if ( !child ) { if(p->data().widthCut()!=ZERO) { if(p->mass() <= p->data().massMax() && p->mass() >= p->data().massMin() ) status = 2; } else { status = 2; } } } } GenParticlePtrT gp = - Traits::newParticle(p->momentum(), p->id(), status, energyUnit); + Traits::newParticle(p->momentum(), p->id(), p->status() ? p->status() : status, energyUnit); if ( p->spinInfo() && p->spinInfo()->hasPolarization() ) { DPair pol = p->spinInfo()->polarization(); Traits::setPolarization(*gp, pol.first, pol.second); } return gp; } template void HepMCConverter::join(tcPPtr parent, tcPPtr child) { Vertex * dec = decv[parent]; Vertex * pro = prov[child]; if ( !pro || !dec ) Throw() << "Found a reference to a ThePEG::Particle which was not in the Event." << Exception::eventerror; if ( pro == dec ) return; while ( !pro->in.empty() ) { dec->in.insert(*(pro->in.begin())); decv[*(pro->in.begin())] = dec; pro->in.erase(pro->in.begin()); } while ( !pro->out.empty() ) { dec->out.insert(*(pro->out.begin())); prov[*(pro->out.begin())] = dec; pro->out.erase(pro->out.begin()); } } template typename HepMCConverter::GenVertexPtrT HepMCConverter::createVertex(Vertex * v) { if ( !v ) Throw() << "Found internal null Vertex." << Exception::abortnow; GenVertexPtrT gv = Traits::newVertex(); // We assume that the vertex position is the average of the decay // vertices of all incoming and the creation vertices of all // outgoing particles in the lab. Note that this will probably not // be useful information for very small distances. LorentzPoint p; for ( tcParticleSet::iterator it = v->in.begin(); it != v->in.end(); ++it ) { p += (**it).labDecayVertex(); Traits::addIncoming(*gv, pmap[*it]); } for ( tcParticleSet::iterator it = v->out.begin(); it != v->out.end(); ++it ) { p += (**it).labVertex(); Traits::addOutgoing(*gv, pmap[*it]); } p /= double(v->in.size() + v->out.size()); Traits::setPosition(*gv, p, lengthUnit); return gv; } template void HepMCConverter::setPdfInfo(const Event & e) { // ids of the partons going into the primary sub process tSubProPtr sub = e.primarySubProcess(); int id1 = sub->incoming().first ->id(); int id2 = sub->incoming().second->id(); // get the event handler tcEHPtr eh = dynamic_ptr_cast(e.handler()); // get the values of x double x1 = eh->lastX1(); double x2 = eh->lastX2(); // get the pdfs pair pdfs; pdfs.first = eh->pdf(sub->incoming().first ); pdfs.second = eh->pdf(sub->incoming().second); // get the scale Energy2 scale = eh->lastScale(); // get the values of the pdfs double xf1 = pdfs.first.xfx(sub->incoming().first->dataPtr(), scale, x1); double xf2 = pdfs.second.xfx(sub->incoming().second->dataPtr(), scale, x2); Traits::setPdfInfo(*geneve, id1, id2, x1, x2, sqrt(scale/GeV2), xf1, xf2); } }