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]); } }