diff --git a/Decay/Dalitz/DalitzResonance.h b/Decay/Dalitz/DalitzResonance.h --- a/Decay/Dalitz/DalitzResonance.h +++ b/Decay/Dalitz/DalitzResonance.h @@ -1,156 +1,158 @@ // -*- C++ -*- #ifndef Herwig_DalitzResonance_H #define Herwig_DalitzResonance_H // // This is the declaration of the DalitzResonance class. // #include "ThePEG/Config/ThePEG.h" #include "DalitzResonance.fh" #include namespace Herwig { using namespace ThePEG; namespace ResonanceType { /** * Enum for the type of resonace */ enum Type {NonResonant=0, Spin0=1,Spin1=3,Spin2=5, Spin0E691=11,Spin1E691=13,Spin2E691=15, BABARf0=21, Spin0Gauss=31, Flattef0=41, Spin0Complex=51, Flattea0=61, FlatteKstar0=71, Sigma=81, Spin1GS=23, Spin0NonResonant=101, Spin1NonResonant=103, Spin2NonResonant=105, Spin0MIPWA=-1, PiPiI2=-11, KMatrix=-21, LASS=-31}; } /** * The DalitzResonance class provides a container class for * information on resonances in multi-body dalitz decays. */ class DalitzResonance: public Base { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ - DalitzResonance() {} + DalitzResonance() : id(0), type(ResonanceType::NonResonant),mass(ZERO),width(ZERO), + daughter1(0),daughter2(0),spectator(0),amp(0.),R(ZERO) + {} /** * Constructor specifiying the parameters */ DalitzResonance(long pid, ResonanceType::Type rtype, Energy m, Energy w, unsigned int d1, unsigned int d2, unsigned int s, double mag, double phi, InvEnergy rr) : id(pid), type(rtype), mass(m),width(w), daughter1(d1),daughter2(d2),spectator(s), amp(mag*exp(Complex(0,phi))), R(rr) {} //@} public: /** * Return the Breit-Wigner times the form factor */ virtual Complex BreitWigner(const Energy & mAB, const Energy & mA, const Energy & mB) const; /** * Output the parameters */ virtual void dataBaseOutput(ofstream & output); /** * Read the parameters for a Dalitz resonance */ static DalitzResonancePtr readResonance(string arg, string & error); public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ DalitzResonance & operator=(const DalitzResonance &) = delete; public: /** * PID of resonant particle */ long id; /** * Type of the resonance */ ResonanceType::Type type; /** * Mass of the resonance */ Energy mass; /** * Width of the resonance */ Energy width; /** * The children */ unsigned int daughter1,daughter2; /** * The spectactor */ unsigned int spectator; /** * The amplitude */ Complex amp; /** * Radius for the Ballt-Weisskopf formfactor */ InvEnergy R; }; } #endif /* Herwig_DalitzResonance_H */ diff --git a/MatrixElement/Gamma/MEGammaGamma2PiPi.cc b/MatrixElement/Gamma/MEGammaGamma2PiPi.cc --- a/MatrixElement/Gamma/MEGammaGamma2PiPi.cc +++ b/MatrixElement/Gamma/MEGammaGamma2PiPi.cc @@ -1,182 +1,176 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the MEGammaGamma2PiPi class. // #include "MEGammaGamma2PiPi.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/MatrixElement/Tree2toNDiagram.h" #include "ThePEG/StandardModel/StandardModelBase.h" #include "Herwig/MatrixElement/HardVertex.h" using namespace Herwig; using ThePEG::Helicity::incoming; MEGammaGamma2PiPi::MEGammaGamma2PiPi() { massOption(vector(2,1)); } void MEGammaGamma2PiPi::getDiagrams() const { tcPDPtr gamma = getParticleData(ParticleID::gamma); tcPDPtr pip = getParticleData(ParticleID::piplus); tcPDPtr pim = pip->CC(); // first t-channel add(new_ptr((Tree2toNDiagram(3),gamma,pip,gamma,1,pip, 2,pim,-1))); // interchange add(new_ptr((Tree2toNDiagram(3),gamma,pip,gamma,2,pip, 1,pim,-2))); } Energy2 MEGammaGamma2PiPi::scale() const { return 2.*sHat()*tHat()*uHat()/(sqr(sHat())+sqr(tHat())+sqr(uHat())); } double MEGammaGamma2PiPi::me2() const { VectorWaveFunction p1w(rescaledMomenta()[0],mePartonData()[0],incoming); VectorWaveFunction p2w(rescaledMomenta()[1],mePartonData()[1],incoming); vector p1,p2; for(unsigned int ix=0;ix<2;++ix) { p1w.reset(2*ix);p1.push_back(p1w); p2w.reset(2*ix);p2.push_back(p2w); } // calculate the matrix element return helicityME(p1,p2,rescaledMomenta(),false); } unsigned int MEGammaGamma2PiPi::orderInAlphaS() const { return 0; } unsigned int MEGammaGamma2PiPi::orderInAlphaEW() const { return 2; } Selector MEGammaGamma2PiPi::diagrams(const DiagramVector & diags) const { Selector sel; for ( DiagramIndex i = 0; i < diags.size(); ++i ) if ( diags[i]->id() == -1 ) sel.insert(meInfo()[0], i); else if ( diags[i]->id() == -2 ) sel.insert(meInfo()[1], i); return sel; } Selector MEGammaGamma2PiPi::colourGeometries(tcDiagPtr ) const { static ColourLines c1(""); Selector sel; sel.insert(1.0, &c1); return sel; } IBPtr MEGammaGamma2PiPi::clone() const { return new_ptr(*this); } IBPtr MEGammaGamma2PiPi::fullclone() const { return new_ptr(*this); } -void MEGammaGamma2PiPi::persistentOutput(PersistentOStream & os) const { -} - -void MEGammaGamma2PiPi::persistentInput(PersistentIStream & is, int) { -} - // The following static variable is needed for the type // description system in ThePEG. -DescribeClass +DescribeNoPIOClass describeHerwigMEGammaGamma2PiPi("Herwig::MEGammaGamma2PiPi", "HwMEGammaGamma.so"); void MEGammaGamma2PiPi::Init() { static ClassDocumentation documentation ("The MEGammaGamma2PiPi class provides a simple matrix element for gamma gamma -> pi+pi-"); } void MEGammaGamma2PiPi::constructVertex(tSubProPtr sub) { // extract the particles in the hard process ParticleVector hard; hard.push_back(sub->incoming().first); hard.push_back(sub->incoming().second); hard.push_back(sub->outgoing()[0]); hard.push_back(sub->outgoing()[1]); // order of particles unsigned int order[4]={0,1,2,3}; // identify the process and calculate matrix element vector p1,p2; if(hard[order[2]]->id()<0) swap(order[2],order[3]); VectorWaveFunction (p1 ,hard[order[0]],incoming,false,true); VectorWaveFunction (p2 ,hard[order[1]],incoming,false,true); p1[1]=p1[2]; p2[1]=p2[2]; vector momenta; for(unsigned int ix=0;ix<4;++ix) momenta.push_back(hard[order[ix]]->momentum()); helicityME(p1,p2,momenta,true); // construct the vertex HardVertexPtr hardvertex=new_ptr(HardVertex()); // set the matrix element for the vertex hardvertex->ME(me_); // set the pointers and to and from the vertex for(unsigned int ix=0;ix<4;++ix) { if(hard[order[ix]]->spinInfo()) hard[order[ix]]->spinInfo()->productionVertex(hardvertex); } } double MEGammaGamma2PiPi::helicityME(vector &p1, vector &p2, const vector & momenta, bool calc) const { // for(unsigned int ix=0;ix<4;++ix) // cerr << mePartonData()[ix]->PDGName() <<" " << meMomenta()[ix]/GeV << " " << meMomenta()[ix].mass()/GeV << "\n"; // double beta = meMomenta()[2].vect().mag()/meMomenta()[2].t(); // double ct = meMomenta()[2].cosTheta(); // double phi = meMomenta()[2].phi(); // scale (external photons so scale in couplings is 0) Energy2 mt(0.*GeV2); // matrix element to be stored if(calc) me_.reset(ProductionMatrixElement(PDT::Spin1,PDT::Spin1, PDT::Spin0,PDT::Spin0)); // calculate the matrix element double output(0.),sumdiag[2]={0.,0.}; Complex diag[2]; Energy2 d1 = momenta[0]*momenta[2], d2 = momenta[1]*momenta[2]; for(unsigned int ihel1=0;ihel1<2;++ihel1) { complex a1 = momenta[2]*p1[ihel1].wave(); complex a2 = momenta[3]*p1[ihel1].wave(); for(unsigned int ihel2=0;ihel2<2;++ihel2) { complex b1 = momenta[2]*p2[ihel2].wave(); complex b2 = momenta[3]*p2[ihel2].wave(); // t-channel diagram diag[0] = a1*b2/d1; // u-channel diagram diag[1] = a2*b1/d2; sumdiag[0] += norm(diag[0]); sumdiag[1] += norm(diag[1]); Complex amp = p1[ihel1].wave()*p2[ihel2].wave()-diag[0]-diag[1]; output += norm(amp); // store the me if needed if(calc) me_(2*ihel1,2*ihel2,0,0) = amp; } } // diagrams DVector save; save.push_back(sumdiag[0]); save.push_back(sumdiag[1]); meInfo(save); // code to test vs the analytic result // Energy2 m2 = sqr(momenta[2].mass()); // double test = 2. + ((d1 + d2)*m2*(-2*d1*d2 + (d1 + d2)*m2))/(sqr(d1)*sqr(d2)); // cerr << "testing ME " << (output-test)/(output+test) << " " << output << " " << test << " " << output/test << "\n"; // coupling factors return output*sqr(SM().alphaEM()*4.*Constants::pi); } diff --git a/MatrixElement/Gamma/MEGammaGamma2PiPi.h b/MatrixElement/Gamma/MEGammaGamma2PiPi.h --- a/MatrixElement/Gamma/MEGammaGamma2PiPi.h +++ b/MatrixElement/Gamma/MEGammaGamma2PiPi.h @@ -1,168 +1,152 @@ // -*- C++ -*- #ifndef Herwig_MEGammaGamma2PiPi_H #define Herwig_MEGammaGamma2PiPi_H // // This is the declaration of the MEGammaGamma2PiPi class. // #include "Herwig/MatrixElement/HwMEBase.h" #include "Herwig/MatrixElement/ProductionMatrixElement.h" #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h" namespace Herwig { using namespace ThePEG; using ThePEG::Helicity::VectorWaveFunction; /** * The MEGammaGamma2PiPi class provides a smiple matrix element for \f$\gamma\gamma\to \pi^+\pi^-$ using scalar QED * * @see \ref MEGammaGamma2PiPiInterfaces "The interfaces" * defined for MEGammaGamma2PiPi. */ class MEGammaGamma2PiPi: public HwMEBase { public: /** * The default constructor. */ MEGammaGamma2PiPi(); public: /** @name Virtual functions required by the MEBase class. */ //@{ /** * Return the order in \f$\alpha_S\f$ in which this matrix * element is given. */ virtual unsigned int orderInAlphaS() const; /** * Return the order in \f$\alpha_{EW}\f$ in which this matrix * element is given. */ virtual unsigned int orderInAlphaEW() const; /** * The matrix element for the kinematical configuration * previously provided by the last call to setKinematics(), suitably * scaled by sHat() to give a dimension-less number. * @return the matrix element scaled with sHat() to give a * dimensionless number. */ virtual double me2() const; /** * Return the scale associated with the last set phase space point. */ virtual Energy2 scale() const; /** * Add all possible diagrams with the add() function. */ virtual void getDiagrams() const; /** * Get diagram selector. With the information previously supplied with the * setKinematics method, a derived class may optionally * override this method to weight the given diagrams with their * (although certainly not physical) relative probabilities. * @param dv the diagrams to be weighted. * @return a Selector relating the given diagrams to their weights. */ virtual Selector diagrams(const DiagramVector & dv) const; /** * Return a Selector with possible colour geometries for the selected * diagram weighted by their relative probabilities. * @param diag the diagram chosen. * @return the possible colour geometries weighted by their * relative probabilities. */ virtual Selector colourGeometries(tcDiagPtr diag) const; //@} /** * Construct the vertex of spin correlations. */ virtual void constructVertex(tSubProPtr); public: - /** @name Functions used by the persistent I/O system. */ - //@{ - /** - * Function used to write out object persistently. - * @param os the persistent output stream written to. - */ - void persistentOutput(PersistentOStream & os) const; - - /** - * Function used to read in object persistently. - * @param is the persistent input stream read from. - * @param version the version number of the object when written. - */ - void persistentInput(PersistentIStream & is, int version); - //@} - /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** * Matrix element for \f$\gamma\gamma\to q\bar{q}\f$ * @param p1 The wavefunctions for the first incoming photon * @param p2 The wavefunctions for the second incoming photon * @param momenta The momenta * @param calc Whether or not to calculate the matrix element */ double helicityME(vector &p1,vector &p2, const vector & momenta, bool calc) const; protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ MEGammaGamma2PiPi & operator=(const MEGammaGamma2PiPi &) = delete; private: /** * Matrix element */ mutable ProductionMatrixElement me_; }; } #endif /* Herwig_MEGammaGamma2PiPi_H */ diff --git a/Shower/Dipole/Merging/Node.cc b/Shower/Dipole/Merging/Node.cc --- a/Shower/Dipole/Merging/Node.cc +++ b/Shower/Dipole/Merging/Node.cc @@ -1,458 +1,456 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the Node class. // #include "Node.h" #include "MergingFactory.h" #include "Merger.h" using namespace Herwig; Node::Node(MatchboxMEBasePtr nodeME, int cutstage, MergerPtr mh) :Interfaced(), thenodeMEPtr(nodeME), thedipol(), theparent(), theCutStage(cutstage), theSubtractedReal(false), theVirtualContribution(false), theMergingHelper(mh) { nodeME->maxMultCKKW(1); nodeME->minMultCKKW(0); } Node::Node(NodePtr deephead, NodePtr head, SubtractionDipolePtr dipol, MatchboxMEBasePtr nodeME, int cutstage) :Interfaced(), thenodeMEPtr(nodeME), thedipol(dipol), theparent(head), theDeepHead(deephead), theCutStage(cutstage), theSubtractedReal(false), theVirtualContribution(false), theMergingHelper() //The subnodes have no merging helper { } -Node::~Node() { } - SubtractionDipolePtr Node::dipole() const { return thedipol; } /** returns the matrix element pointer */ const MatchboxMEBasePtr Node::nodeME() const { return thenodeMEPtr; } /** access the matrix element pointer */ MatchboxMEBasePtr Node::nodeME() { return thenodeMEPtr; } pair Node::getInOut( ){ PVector in; const auto me= nodeME(); const auto pd=me->mePartonData(); for( auto i : {0 , 1} ) in.push_back(pd[i]->produceParticle( me->lastMEMomenta()[i] ) ); PVector out; for ( size_t i = 2;i< pd.size();i++ ){ PPtr p = pd[i]->produceParticle( me->lastMEMomenta()[i] ); out.push_back( p ); } return { in , out }; } int Node::legsize() const {return nodeME()->legsize();} NodePtr Node::randomChild() { return thechildren[UseRandom::irnd(thechildren.size())]; } bool Node::allAbove(Energy pt) { for (NodePtr child : thechildren) if ( child->pT() < pt ) return false; return true; } Energy Node::maxChildPt(){ Energy maxi=-1*GeV; for (NodePtr child : thechildren)maxi=max(child->pT(),maxi); return maxi; } bool Node::isInHistoryOf(NodePtr other) { while (other->parent()) { if (other == this) return true; other = other->parent(); } return false; } void Node::flushCaches() { if (didflush) return; didflush=true; for ( auto const & ch: thechildren) { ch->xcomb()->clean(); ch->nodeME()->flushCaches(); ch->flushCaches(); } } void Node::setKinematics() { for (auto const & ch: thechildren) { ch->dipole()->setXComb(ch->xcomb()); ch->dipole()->setKinematics(); ch->nodeME()->setKinematics(); ch->setKinematics(); } } void Node::clearKinematics() { for (auto const & ch: thechildren) { ch->dipole()->setXComb(ch->xcomb()); ch->nodeME()->clearKinematics(); ch->dipole()->clearKinematics(); ch->clearKinematics(); } } bool Node::generateKinematics(const double *r, bool directCut) { didflush=false; // If there are no children to the child process we are done. if(children().empty()) return true; assert(parent()); if ( ! directCut && pT() < deepHead()->MH()->mergePt()) { // Real emission: // If there are children to the child process, // we now require that all subsequent children // with pt < merging scale are in their ME region. // Since the possible children of the real emission // contribution are now in their ME region, // it is clear that a second clustering is possible // -- modulo phase space restrictions. // Therefore the real emission contribution // are unitarised and the cross section is // hardly modified. auto inOutPair = getInOut(); NodePtr rc = randomChild(); rc->dipole()->setXComb(rc->xcomb()); if(!rc->dipole()->generateKinematics(r))assert(false); // If not in ME -> return false if(!deepHead()->MH()->matrixElementRegion( inOutPair.first , inOutPair.second , rc->pT() , deepHead()->MH()->mergePt() ) )return false; } for (auto & ch : children() ) { ch->dipole()->setXComb(ch->xcomb()); if ( !ch->dipole()->generateKinematics(r) ) { assert(false); } ch->generateKinematics( r, true); } return true; } bool Node::firstgenerateKinematics(const double *r, bool directCut) { didflush=false; // This is called form the merging helper for the first node. So: assert(!parent()); assert(xcomb()); if(MH()->treefactory()->nonQCDCuts()){ tcPDVector outdata(xcomb()->mePartonData().begin()+2, xcomb()->mePartonData().end()); vector outmomenta(xcomb()->meMomenta().begin()+2, xcomb()->meMomenta().end()); if ( !MH()->treefactory()->nonQCDCuts()->passCuts(outdata,outmomenta, xcomb()->mePartonData()[0], xcomb()->mePartonData()[1]) ) return false; } ///// This should not be needed!!! ///// ( Warning inMerger::matrixElementRegion gets triggered.) flushCaches(); //Set here the new merge Pt for the next phase space point.( Smearing!!!) MH()->smearMergePt(); // If there are no children to this node, we are done here: if (children().empty()) return true; // directCut is for born and for virtual contributions. // if directCut is true, then cut on the first ME region. // call recursiv generate kinematics for subsequent nodes. if ( directCut ){ auto inOutPair = getInOut(); NodePtr rc = randomChild(); rc->dipole()->setXComb(rc->xcomb()); if ( !rc->dipole()->generateKinematics(r) ) { return false; } rc->nodeME()->setXComb(rc->xcomb()); if(!MH()->matrixElementRegion( inOutPair.first , inOutPair.second , rc->pT() , MH()->mergePt() ) ){ return false; } } for (auto const & ch: thechildren) { ch->dipole()->setXComb(ch->xcomb()); if ( !ch->dipole()->generateKinematics(r) ) { cout<<"\nCould not generate dipole kinematics";;return false; } if( ! ch->generateKinematics(r,directCut) )return false; } return true; } StdXCombPtr Node::xcomb() const { assert(thexcomb); return thexcomb; } StdXCombPtr Node::xcomb(){ if(thexcomb)return thexcomb; assert(parent()); thexcomb=dipole()->makeBornXComb(parent()->xcomb()); xcomb()->head(parent()->xcomb()); dipole()->setXComb(thexcomb); return thexcomb; } void Node::setXComb(tStdXCombPtr xc) { assert ( !parent() ); thexcomb=xc; assert(thexcomb->lastParticles().first); } #include "Herwig/MatrixElement/Matchbox/Base/DipoleRepository.h" void Node::birth(const vector & vec) { // produce the children vector dipoles = nodeME()->getDipoles(DipoleRepository::dipoles( nodeME()->factory()->dipoleSet()), vec, true); for ( auto const & dip : dipoles ) { dip->doSubtraction(); NodePtr node = new_ptr(Node(theDeepHead, this, dip, dip->underlyingBornME(), theDeepHead->cutStage())); thechildren.push_back(node); } } vector Node::getNextOrderedNodes(bool normal, double hardScaleFactor) const { vector temp = children(); vector res; for (NodePtr const & child : children()) { if(deepHead()->MH()->mergePt()>child->pT()) { res.clear(); return res; } } for (NodePtr const & child: children()) { if (parent()&& normal) { if ( child->pT() < pT() ) { continue; } } if ( child->children().size() != 0 ) { for (NodePtr itChild: child->children()) { if( itChild->pT() > child->pT()&&child->inShowerPS(itChild->pT()) ) { res.push_back(child); break; } } } else { const auto sc=child->nodeME()->factory()->scaleChoice(); sc->setXComb(child->xcomb()); if ( sqr(hardScaleFactor)* sc->renormalizationScale() >= sqr(child->pT()) && child->inShowerPS(hardScaleFactor*sqrt(sc->renormalizationScale()))) { res.push_back(child); } } } return res; } bool Node::inShowerPS(Energy hardpT)const { // Here we decide if the current phase space // point can be reached from the underlying Node. // Full phase space available -> Tilde Kinematic is always fine. if(deepHead()->MH()->openZBoundaries()==1) return true; double z_ = dipole()->lastZ(); // restrict according to hard scale if(deepHead()->MH()->openZBoundaries()==0){ pair zbounds = dipole()->tildeKinematics()->zBounds(pT(), hardpT); return (zbounds.first temp = getNextOrderedNodes(normal, hardScaleFactor); Energy minpt = Constants::MaxEnergy; Selector subprosel; while (temp.size() != 0) { minpt = Constants::MaxEnergy; subprosel.clear(); for (NodePtr const & child : temp) { // colour basis is not strictly needed for evaluating the colour correlator //assert(deepHead()->MH()->largeNBasis()); if( child->dipole()->underlyingBornME()->largeNColourCorrelatedME2( {child->dipole()->bornEmitter(), child->dipole()->bornSpectator()}, deepHead()->MH()->largeNBasis()) != 0. ) { double weight = 1.; if ( deepHead()->MH()->chooseHistory() == 0 ) weight = abs(child->dipole()->dSigHatDR()/nanobarn); else if ( deepHead()->MH()->chooseHistory() == 1 ) weight = abs(child->dipole()->dSigHatDR()/child->nodeME()->dSigHatDRB()); else if ( deepHead()->MH()->chooseHistory() == 2 ) weight = 1.; else if ( deepHead()->MH()->chooseHistory() == 3 ) weight = 1_GeV/child->pT(); else assert(false); if(weight != 0.) { subprosel.insert(weight , child); minpt = min(minpt, child->pT()); } } } if (subprosel.empty()) return res; res = subprosel.select(UseRandom::rnd()); temp = res->getNextOrderedNodes(true, hardScaleFactor); } return res; } pair Node::calcDipandPS(Energy scale)const { return dipole()->dipandPs(sqr(scale), deepHead()->MH()->largeNBasis()); } CrossSection Node::calcPs(Energy scale)const { return dipole()->ps(sqr(scale), deepHead()->MH()->largeNBasis()); } CrossSection Node::calcDip(Energy scale)const { return dipole()->dip(sqr(scale)); } IBPtr Node::clone() const { return new_ptr(*this); } IBPtr Node::fullclone() const { return new_ptr(*this); } #include "ThePEG/Persistency/PersistentOStream.h" void Node::persistentOutput(PersistentOStream & os) const { os << thexcomb<< thenodeMEPtr<< thedipol<< thechildren<< theparent<< theProjector<< theDeepHead<< theCutStage<< ounit(theRunningPt, GeV)<< theSubtractedReal<< theVirtualContribution<< theMergingHelper; } #include "ThePEG/Persistency/PersistentIStream.h" void Node::persistentInput(PersistentIStream & is, int) { is >> thexcomb>> thenodeMEPtr>> thedipol>> thechildren>> theparent>> theProjector>> theDeepHead>> theCutStage>> iunit(theRunningPt, GeV)>> theSubtractedReal>> theVirtualContribution>> theMergingHelper; } // *** Attention *** The following static variable is needed for the type // description system in ThePEG. Please check that the template arguments // are correct (the class and its base class), and that the constructor // arguments are correct (the class name and the name of the dynamically // loadable library where the class implementation can be found). #include "ThePEG/Utilities/DescribeClass.h" DescribeClass describeHerwigNode("Herwig::Node", "HwDipoleShower.so"); void Node::Init() { static ClassDocumentation documentation("There is no documentation for the Node class"); } diff --git a/Shower/Dipole/Merging/Node.h b/Shower/Dipole/Merging/Node.h --- a/Shower/Dipole/Merging/Node.h +++ b/Shower/Dipole/Merging/Node.h @@ -1,242 +1,240 @@ // -*- C++ -*- #ifndef Herwig_Node_H #define Herwig_Node_H // // This is the declaration of the Node class. // #include "Node.fh" #include "MergingFactory.fh" #include "Merger.h" #include "ThePEG/Config/ThePEG.h" #include "ThePEG/Config/std.h" #include "ThePEG/Interface/Interfaced.h" #include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.fh" #include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh" #include "Herwig/Shower/Dipole/Base/DipoleEventRecord.h" #include "ThePEG/MatrixElement/MEBase.h" #include namespace Herwig { using namespace ThePEG; /** * The Node class represents a part of the shower history * in the merging procedure. * Each node can search for its "childen" in the "birth" step. * The children processes that can be created by clustering two legs * with a given spectator -- just like in the shower. * To perform the "birth" a vector of processes is given to the node * and the possible dipoles are given by the factory object. * * * @see \ref NodeInterfaces "The interfaces" * defined for Node. */ class Node : public Interfaced { public: /** @name Standard constructors and destructors. */ //@{ Node (){}; // constructor for first nodes Node(MatchboxMEBasePtr nodeME , int cutstage , MergerPtr mh ); // another constructor for underlying nodes Node(NodePtr deephead, NodePtr head, SubtractionDipolePtr dipol, MatchboxMEBasePtr nodeME, int cutstage); - /// The destructor. - virtual ~Node(); //@} public: // get children from vector void birth(const vector & vec); /// recursive setXComb. proStage is the number of clusterings /// before the projectors get filled. void setXComb(tStdXCombPtr xc); /// calculate the dipole and ps approximation pair calcDipandPS(Energy scale)const; /// calculate the ps approximation CrossSection calcPs(Energy scale)const; /// calculate the dipole CrossSection calcDip(Energy scale)const; /// recursive flush caches and clean up XCombs. void flushCaches(); /// recursive clearKinematics void clearKinematics(); /// recursive setKinematics void setKinematics(); /// recursive generateKinematics using tilde kinematics of the dipoles bool generateKinematics(const double *r, bool directCut); /// generate the kinamatics of the first node bool firstgenerateKinematics(const double *r, bool directCut); //return the ME const MatchboxMEBasePtr nodeME() const; //return the node ME MatchboxMEBasePtr nodeME(); //return the parent Node NodePtr parent() const {return theparent;} /// vector of children nodes created in birth vector< NodePtr > children() const {return thechildren;} //pick a random child (flat) NodePtr randomChild(); /// true if all children show scales above pt bool allAbove(Energy pt); /// return maximum of all child pts. Energy maxChildPt(); /// true if the node is in the history of other. bool isInHistoryOf(NodePtr other); /// legsize of the node ME int legsize() const; /// set the first node (first men). only use in factory void deepHead(NodePtr deephead) {theDeepHead = deephead;} /// return the first node NodePtr deepHead() const {return theDeepHead;} /// returns the dipol of the node. SubtractionDipolePtr dipole() const; /// return the xcomb StdXCombPtr xcomb() const; /// return the xcomb (if not created, create one from head) StdXCombPtr xcomb() ; /// return the current running pt Energy runningPt() const { return theRunningPt; } /// set the current running pt void runningPt(Energy x) { theRunningPt=x; } /// return the cut stage to cut on merging pt in generate kinematics int cutStage() const { return theCutStage; } /// get a vector of the next nodes, ordered in pt (and in parton shower phace space) vector getNextOrderedNodes(bool normal=true, double hardscalefactor=1.) const; //true if the node is in shower history for a given pt bool inShowerPS(Energy hardpt)const; //get the history NodePtr getHistory(bool normal=true, double hardscalefactor=1.); //true if node correspond to a subtracted real. bool subtractedReal() const {return theSubtractedReal;} /// set if node correspont to a subtracted real. void subtractedReal(bool x) { theSubtractedReal = x;} //true if node correspond to a virtual contribution. bool virtualContribution() const { return theVirtualContribution ;} /// set if node correspont to a virtual contribution. void virtualContribution(bool x) {theVirtualContribution = x;} //pointer to the merging helper MergerPtr MH()const{return theMergingHelper;} /// set the merging helper void MH(MergerPtr a){theMergingHelper=a;} /// pT of the dipole Energy pT()const{return dipole()->lastPt();} /// get incoming and outgoing particles (TODO: expensive) pair getInOut(); private: /// the Matrixelement representing this node. MatchboxMEBasePtr thenodeMEPtr; /// the dipol used to substract /// and generate kinematics using tilde kinematics SubtractionDipolePtr thedipol; /// the parent node NodePtr theparent; /// The godfather node of whole tree.(Firstnode) NodePtr theDeepHead; /** * The CutStage is number of clusterings which are possible without * introducing a merging scale to cut away singularities. * -> subtracted MEs have the CutStage 1. * -> virtual and normal tree level ME get 0. */ int theCutStage; /// flag to tell if node is subtracted real bool theSubtractedReal; /// flag to tell if node is virtual contribution bool theVirtualContribution; /// the merging helper (should be static) MergerPtr theMergingHelper; //the xcomb of the node StdXCombPtr thexcomb; /// vector of the children node vector< NodePtr > thechildren; /// the current running pt Energy theRunningPt; /// The nodes of the projection stage. NodePtr theProjector; /// flag not to enter infinite loop. (There should be a better solution...) bool didflush=false; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ Node & operator=(const Node &) = delete; }; } #endif /* Herwig_Node_H */ diff --git a/Shower/QTilde/Base/PartnerFinder.cc b/Shower/QTilde/Base/PartnerFinder.cc --- a/Shower/QTilde/Base/PartnerFinder.cc +++ b/Shower/QTilde/Base/PartnerFinder.cc @@ -1,699 +1,699 @@ // -*- C++ -*- // // PartnerFinder.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig 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 PartnerFinder class. // #include "PartnerFinder.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Shower/QTilde/Base/ShowerParticle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Utilities/Debug.h" #include "ThePEG/Utilities/DescribeClass.h" using namespace Herwig; DescribeClass describePartnerFinder ("Herwig::PartnerFinder","HwShower.so"); // some useful functions to avoid using #define namespace { // return bool if final-state particle inline bool FS(const tShowerParticlePtr a) { return a->isFinalState(); } // return colour line pointer inline Ptr::transient_pointer CL(const tShowerParticlePtr a, unsigned int index=0) { return a->colourInfo()->colourLines().empty() ? ThePEG::tColinePtr() : const_ptr_cast(a->colourInfo()->colourLines()[index]); } // return colour line size inline size_t CLSIZE(const tShowerParticlePtr a) { return a->colourInfo()->colourLines().size(); } inline Ptr::transient_pointer ACL(const tShowerParticlePtr a, unsigned int index=0) { return a->colourInfo()->antiColourLines().empty() ? ThePEG::tColinePtr() : const_ptr_cast(a->colourInfo()->antiColourLines()[index]); } inline size_t ACLSIZE(const tShowerParticlePtr a) { return a->colourInfo()->antiColourLines().size(); } } void PartnerFinder::persistentOutput(PersistentOStream & os) const { os << partnerMethod_ << QEDPartner_ << scaleChoice_; } void PartnerFinder::persistentInput(PersistentIStream & is, int) { is >> partnerMethod_ >> QEDPartner_ >> scaleChoice_; } void PartnerFinder::Init() { static ClassDocumentation documentation ("This class is responsible for finding the partners for each interaction types ", "and within the evolution scale range specified by the ShowerVariables ", "then to determine the initial evolution scales for each pair of partners."); static Switch interfacePartnerMethod ("PartnerMethod", "Choice of partner finding method for gluon evolution.", &PartnerFinder::partnerMethod_, 0, false, false); static SwitchOption interfacePartnerMethodRandom (interfacePartnerMethod, "Random", "Choose partners of a gluon randomly.", 0); static SwitchOption interfacePartnerMethodMaximum (interfacePartnerMethod, "Maximum", "Choose partner of gluon with largest angle.", 1); static Switch interfaceQEDPartner ("QEDPartner", "Control of which particles to use as the partner for QED radiation", &PartnerFinder::QEDPartner_, 0, false, false); static SwitchOption interfaceQEDPartnerAll (interfaceQEDPartner, "All", "Consider all possible choices which give a positive contribution" " in the soft limit.", 0); static SwitchOption interfaceQEDPartnerIIandFF (interfaceQEDPartner, "IIandFF", "Only allow initial-initial or final-final combinations", 1); static SwitchOption interfaceQEDPartnerIF (interfaceQEDPartner, "IF", "Only allow initial-final combinations", 2); static Switch interfaceScaleChoice ("ScaleChoice", "The choice of the evolution scales", &PartnerFinder::scaleChoice_, 0, false, false); static SwitchOption interfaceScaleChoicePartner (interfaceScaleChoice, "Partner", "Scale of all interactions is that of the evolution partner", 0); static SwitchOption interfaceScaleChoiceDifferent (interfaceScaleChoice, "Different", "Allow each interaction to have different scales", 1); } void PartnerFinder::setInitialEvolutionScales(const ShowerParticleVector &particles, const bool isDecayCase, ShowerInteraction type, const bool setPartners) { // clear the existing partners for(ShowerParticleVector::const_iterator cit = particles.begin(); cit != particles.end(); ++cit) (*cit)->clearPartners(); // set them if(type==ShowerInteraction::QCD) { setInitialQCDEvolutionScales(particles,isDecayCase,setPartners); } else if(type==ShowerInteraction::QED) { setInitialQEDEvolutionScales(particles,isDecayCase,setPartners); } else if(type==ShowerInteraction::EW) { setInitialEWEvolutionScales(particles,isDecayCase,false); } else if(type==ShowerInteraction::QEDQCD) { setInitialQCDEvolutionScales(particles,isDecayCase,setPartners); setInitialQEDEvolutionScales(particles,isDecayCase,false); } else if(type==ShowerInteraction::ALL) { setInitialQCDEvolutionScales(particles,isDecayCase,setPartners); setInitialQEDEvolutionScales(particles,isDecayCase,false); setInitialEWEvolutionScales(particles,isDecayCase,false); } else assert(false); // \todo EW scales here // print out for debugging if(Debug::level>=10) { for(ShowerParticleVector::const_iterator cit = particles.begin(); cit != particles.end(); ++cit) { generator()->log() << "Particle: " << **cit << "\n"; if(!(**cit).partner()) continue; generator()->log() << "Primary partner: " << *(**cit).partner() << "\n"; for(vector::const_iterator it= (**cit).partners().begin(); it!=(**cit).partners().end();++it) { generator()->log() << static_cast(it->type) << " " << it->weight << " " << it->scale/GeV << " " << *(it->partner) << "\n"; } } generator()->log() << flush; } } void PartnerFinder::setInitialQCDEvolutionScales(const ShowerParticleVector &particles, const bool isDecayCase, const bool setPartners) { // Loop over particles and consider only coloured particles which don't // have already their colour partner fixed and that don't have children // (the latter requirement is relaxed in the case isDecayCase is true). // Build a map which has as key one of these particles (i.e. a pointer // to a ShowerParticle object) and as a corresponding value the vector // of all its possible *normal* candidate colour partners, defined as follows: // --- have colour, and no children (this is not required in the case // isDecayCase is true); // --- if both are initial (incoming) state particles, then the (non-null) colourLine() // of one of them must match the (non-null) antiColourLine() of the other. // --- if one is an initial (incoming) state particle and the other is // a final (outgoing) state particle, then both must have the // same (non-null) colourLine() or the same (non-null) antiColourLine(); // Notice that this definition exclude the special case of baryon-violating // processes (as in R-parity Susy), which will show up as particles // without candidate colour partners, and that we will be treated a part later // (this means that no modifications of the following loop is needed! for ( const auto & sp : particles ) { // Skip colourless particles if(!sp->data().coloured()) continue; // find the partners auto partners = findQCDPartners(sp,particles); // must have a partner if(partners.empty()) { throw Exception() << "`Failed to make colour connections in " << "PartnerFinder::setQCDInitialEvolutionScales" << *sp << Exception::eventerror; } // Calculate the evolution scales for all possible pairs of of particles vector > scales; int position = -1; for(size_t ix=0; ix< partners.size(); ++ix) { scales.push_back(calculateInitialEvolutionScales(ShowerPPair(sp, partners[ix].second),isDecayCase)); if (!setPartners && partners[ix].second) position = ix; } assert(setPartners || position >= 0); // set partners if required if (setPartners) { // In the case of more than one candidate colour partners, // there are now two approaches to choosing the partner. The // first method is based on two assumptions: // 1) the choice of which is the colour partner is done // *randomly* between the available candidates. // 2) the choice of which is the colour partner is done // *independently* from each particle: in other words, // if for a particle "i" its selected colour partner is // the particle "j", then the colour partner of "j" // does not have to be necessarily "i". // The second method always chooses the furthest partner // for hard gluons and gluinos. // random choice if( partnerMethod_ == 0 ) { // random choice of partner position = UseRandom::irnd(partners.size()); } // take the one with largest angle else if (partnerMethod_ == 1 ) { if (sp->perturbative() == 1 && sp->dataPtr()->iColour()==PDT::Colour8 ) { assert(partners.size()==2); // Determine largest angle double maxAngle(0.); for(unsigned int ix=0;ixmomentum().vect(). angle(partners[ix].second->momentum().vect()); if(angle>maxAngle) { maxAngle = angle; position = ix; } } } else position = UseRandom::irnd(partners.size()); } else assert(false); // set the evolution partner sp->partner(partners[position].second); } // primary partner set, set the others and do the scale for(size_t ix=0; ixaddPartner(ShowerParticle::EvolutionPartner(partners[ix].second,1.,partners[ix].first, scales[ix].first)); } // set scales for all interactions to that of the partner, default Energy scale = scales[position].first; for(unsigned int ix=0;ixscales().QCD_c = sp->scales().QCD_c_noAO = (scaleChoice_==0 ? scale : scales[ix].first); } else if(partners[ix].first==ShowerPartnerType::QCDAntiColourLine) { sp->scales().QCD_ac = sp->scales().QCD_ac_noAO = (scaleChoice_==0 ? scale : scales[ix].first); } else assert(false); } } } void PartnerFinder::setInitialQEDEvolutionScales(const ShowerParticleVector &particles, const bool isDecayCase, const bool setPartners) { // loop over all the particles for(const auto & sp : particles) { // not charged or photon continue if(!sp->dataPtr()->charged() && sp->dataPtr()->id()!=ParticleID::gamma) continue; // find the potential partners vector > partners = findQEDPartners(sp,particles,isDecayCase); if(partners.empty()) { throw Exception() << "Failed to find partner in " << "PartnerFinder::setQEDInitialEvolutionScales" << *sp << Exception::eventerror; } // calculate the probabilities double prob(0.); for(unsigned int ix=0;ixpartner()) { for(unsigned int ix=0;ixpartner()==partners[ix].second) { position = ix; break; } } } // set the partner if(setPartners||!sp->partner()||position<0) { prob = UseRandom::rnd(); for(unsigned int ix=0;ixprob) { position = ix; break; } prob -= partners[ix].first; } if(position>=0&&(setPartners||!sp->partner())) { sp->partner(partners[position].second); } } // must have a partner if(position<0) throw Exception() << "Failed to find partner in " << "PartnerFinder::setQEDInitialEvolutionScales" << *sp << Exception::eventerror; // Calculate the evolution scales for all possible pairs of of particles vector > scales; for(unsigned int ix=0;ix< partners.size();++ix) { scales.push_back(calculateInitialEvolutionScales(ShowerPPair(sp,partners[ix].second), isDecayCase)); } // store all the possible partners for(unsigned int ix=0;ixaddPartner(ShowerParticle::EvolutionPartner(partners[ix].second, partners[ix].first, ShowerPartnerType::QED, scales[ix].first)); } // set scales sp->scales().QED = scales[position].first; sp->scales().QED_noAO = scales[position].first; } } pair PartnerFinder:: calculateInitialEvolutionScales(const ShowerPPair &particlePair, const bool isDecayCase, int key) { bool FS1=FS(particlePair.first),FS2= FS(particlePair.second); if(FS1 && FS2){ return calculateFinalFinalScales(particlePair.first->momentum(),particlePair.second->momentum(), key); } else if(FS1 && !FS2) { pair rval = calculateInitialFinalScales(particlePair.second->momentum(), particlePair.first->momentum(), isDecayCase); return { rval.second, rval.first }; } else if(!FS1 &&FS2) return calculateInitialFinalScales(particlePair.first->momentum(),particlePair.second->momentum(),isDecayCase); else return calculateInitialInitialScales(particlePair.first->momentum(),particlePair.second->momentum()); } vector< pair > PartnerFinder::findQCDPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles) { vector< pair > partners; for(const auto & sp : particles) { if(!sp->data().coloured() || particle==sp) continue; // one initial-state and one final-state particle if(FS(particle) != FS(sp)) { // loop over all the colours of both particles for(size_t ix=0; ixsourceNeighbours().first) { tColinePair cpair = col->sourceNeighbours(); for(const auto & sp : particles) { if(( FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second))|| (!FS(sp) && (ACL(sp) == cpair.first || ACL(sp) == cpair.second ))) { partners.push_back({ ShowerPartnerType:: QCDColourLine, sp }); } } } else if(col&&col->sinkNeighbours().first) { tColinePair cpair = col->sinkNeighbours(); for(const auto & sp : particles) { if(( FS(sp) && (ACL(sp) == cpair.first || ACL(sp) == cpair.second))|| (!FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second))) { partners.push_back({ ShowerPartnerType:: QCDColourLine, sp }); } } } col = ACL(particle); if(FS(particle)&&col&&col->sinkNeighbours().first) { tColinePair cpair = col->sinkNeighbours(); for(const auto & sp : particles) { if(( FS(sp) && (ACL(sp) == cpair.first || ACL(sp) == cpair.second))|| (!FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second ))) { partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp }); } } } else if(col&&col->sourceNeighbours().first) { tColinePair cpair = col->sourceNeighbours(); for(const auto & sp : particles) { if(( FS(sp) && ( CL(sp) == cpair.first || CL(sp) == cpair.second))|| (!FS(sp) && (ACL(sp) == cpair.first ||ACL(sp) == cpair.second))) { partners.push_back({ ShowerPartnerType::QCDAntiColourLine, sp }); } } } } // return the partners return partners; } vector< pair > PartnerFinder::findQEDPartners(tShowerParticlePtr particle, const ShowerParticleVector & particles, const bool isDecayCase) { vector< pair > partners; const double pcharge = particle->id()==ParticleID::gamma ? 1 : double(particle->data().iCharge()); vector< pair > photons; for (const auto & sp : particles) { if (particle == sp) continue; if (sp->id()==ParticleID::gamma) photons.push_back(make_pair(1.,sp)); if (!sp->data().charged() ) continue; double charge = pcharge*double((sp)->data().iCharge()); if ( FS(particle) != FS(sp) ) charge *=-1.; if ( QEDPartner_ != 0 && !isDecayCase ) { // only include II and FF as requested if ( QEDPartner_ == 1 && FS(particle) != FS(sp) ) continue; // only include IF is requested else if (QEDPartner_ == 2 && FS(particle) == FS(sp) ) continue; } if (particle->id()==ParticleID::gamma) charge = -abs(charge); // only keep positive dipoles if (charge<0.) partners.push_back({ -charge, sp }); } if (particle->id()==ParticleID::gamma && partners.empty()) return photons; return partners; } void PartnerFinder::setInitialEWEvolutionScales(const ShowerParticleVector &particles, const bool isDecayCase, const bool setPartners) { // loop over all the particles for(ShowerParticleVector::const_iterator cit = particles.begin(); cit != particles.end(); ++cit) { // if not weakly interacting continue if( !weaklyInteracting( (**cit).dataPtr())) continue; // find the potential partners vector > partners = findEWPartners(*cit,particles,isDecayCase); if(partners.empty()) { throw Exception() << "Failed to find partner in " << "PartnerFinder::setEWInitialEvolutionScales" << (**cit) << Exception::eventerror; } // calculate the probabilities double prob(0.); for(unsigned int ix=0;ixpartner()) { for(unsigned int ix=0;ixpartner()==partners[ix].second) { position = ix; break; } } } // set the partner if(setPartners||!(*cit)->partner()||position<0) { prob = UseRandom::rnd(); for(unsigned int ix=0;ixprob) { position = ix; break; } prob -= partners[ix].first; } if(position>=0&&(setPartners||!(*cit)->partner())) { (*cit)->partner(partners[position].second); } } // must have a partner if(position<0) throw Exception() << "Failed to find partner in " << "PartnerFinder::setEWInitialEvolutionScales" << (**cit) << Exception::eventerror; // Calculate the evolution scales for all possible pairs of of particles vector > scales; for(unsigned int ix=0;ix< partners.size();++ix) { scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second), isDecayCase, 0)); } // store all the possible partners for(unsigned int ix=0;ix > PartnerFinder::findEWPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles, - const bool isDecayCase ) { + const bool ) { vector< pair > partners; for(ShowerParticlePtr partner : particles) { if(particle==partner || !weaklyInteracting(partner->dataPtr())) continue; partners.push_back(make_pair(1.,partner)); } return partners; } pair PartnerFinder::calculateFinalFinalScales( const Lorentz5Momentum & p1, const Lorentz5Momentum & p2, int key) { static const double eps=1e-7; // Using JHEP 12(2003)045 we find that we need ktilde = 1/2(1+b-c+lambda) // ktilde = qtilde^2/Q^2 therefore qtilde = sqrt(ktilde*Q^2) // find momenta in rest frame of system // calculate quantities for the scales Energy2 Q2 = (p1+p2).m2(); double b = p1.mass2()/Q2; double c = p2.mass2()/Q2; if(b<0.) { if(b<-eps) { throw Exception() << "Negative Mass squared b = " << b << "in PartnerFinder::calculateFinalFinalScales()" << Exception::eventerror; } b = 0.; } if(c<0.) { if(c<-eps) { throw Exception() << "Negative Mass squared c = " << c << "in PartnerFinder::calculateFinalFinalScales()" << Exception::eventerror; } c = 0.; } // KMH & PR - 16 May 2008 - swapped lambda calculation from // double lam=2.*p1.vect().mag()/Q; to sqrt(kallen(1,b,c)), // which should be identical for p1 & p2 onshell in their COM // but in the inverse construction for the Nason method, this // was not the case, leading to misuse. const double lam=sqrt((1.+sqrt(b)+sqrt(c))*(1.-sqrt(b)-sqrt(c)) *(sqrt(b)-1.-sqrt(c))*(sqrt(c)-1.-sqrt(b))); Energy firstQ, secondQ; double kappab(0.), kappac(0.); //key = 0; // symmetric case pre-selection switch(key) { case 0: // symmetric case firstQ = sqrt(0.5*Q2*(1.+b-c+lam)); secondQ = sqrt(0.5*Q2*(1.-b+c+lam)); break; case 1: // maximum emission from both legs kappab=4.*(1.-2.*sqrt(c)-b+c); kappac=4.*(1.-2.*sqrt(b)-c+b); firstQ = sqrt(Q2*kappab); secondQ = sqrt(Q2*kappac); break; default: assert(false); } // calculate the scales return pair(firstQ, secondQ); } pair PartnerFinder::calculateInitialFinalScales(const Lorentz5Momentum& pb, const Lorentz5Momentum& pc, const bool isDecayCase) { if(!isDecayCase) { // In this case from JHEP 12(2003)045 we find the conditions // ktilde_b = (1+c) and ktilde_c = (1+2c) // We also find that c = m_c^2/Q^2. The process is a+b->c where // particle a is not colour connected (considered as a colour singlet). // Therefore we simply find that q_b = sqrt(Q^2+m_c^2) and // q_c = sqrt(Q^2+2 m_c^2) // We also assume that the first particle in the pair is the initial // state particle and the second is the final state one c const Energy2 mc2 = sqr(pc.mass()); const Energy2 Q2 = -(pb-pc).m2(); return { sqrt(Q2+mc2), sqrt(Q2+2*mc2) }; } else { // In this case from JHEP 12(2003)045 we find, for the decay // process b->c+a(neutral), the condition // (ktilde_b-1)*(ktilde_c-c)=(1/4)*sqr(1-a+c+lambda). // We also assume that the first particle in the pair is the initial // state particle (b) and the second is the final state one (c). // - We find maximal phase space coverage through emissions from // c if we set ktilde_c = 4.*(sqr(1.-sqrt(a))-c) // - We find the most 'symmetric' way to populate the phase space // occurs for (ktilde_b-1)=(ktilde_c-c)=(1/2)*(1-a+c+lambda) // - We find the most 'smooth' way to populate the phase space // occurs for... Energy2 mb2(sqr(pb.mass())); double a=(pb-pc).m2()/mb2; double c=sqr(pc.mass())/mb2; double lambda = 1. + a*a + c*c - 2.*a - 2.*c - 2.*a*c; lambda = sqrt(max(lambda,0.)); const double PROD = 0.25*sqr(1. - a + c + lambda); int key = 0; double ktilde_b, ktilde_c, cosi = 0.; switch (key) { case 0: // the 'symmetric' choice ktilde_c = 0.5*(1-a+c+lambda) + c ; ktilde_b = 1.+PROD/(ktilde_c-c) ; break; case 1: // the 'maximal' choice ktilde_c = 4.0*(sqr(1.-sqrt(a))-c); ktilde_b = 1.+PROD/(ktilde_c-c) ; break; case 2: // the 'smooth' choice c = max(c,1.*GeV2/mb2); cosi = (sqr(1-sqrt(c))-a)/lambda; ktilde_b = 2.0/(1.0-cosi); ktilde_c = (1.0-a+c+lambda)*(1.0+c-a-lambda*cosi)/(2.0*(1.0+cosi)); break; } return { sqrt(mb2*ktilde_b), sqrt(mb2*ktilde_c) }; } } pair PartnerFinder::calculateInitialInitialScales(const Lorentz5Momentum& p1, const Lorentz5Momentum& p2) { // This case is quite simple. From JHEP 12(2003)045 we find the condition // that ktilde_b = ktilde_c = 1. In this case we have the process // b+c->a so we need merely boost to the CM frame of the two incoming // particles and then qtilde is equal to the energy in that frame const Energy Q = sqrt((p1+p2).m2()); return {Q,Q}; } diff --git a/Shower/QTilde/Kinematics/ShowerKinematics.cc b/Shower/QTilde/Kinematics/ShowerKinematics.cc --- a/Shower/QTilde/Kinematics/ShowerKinematics.cc +++ b/Shower/QTilde/Kinematics/ShowerKinematics.cc @@ -1,71 +1,71 @@ // -*- C++ -*- // // ShowerKinematics.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig 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 ShowerKinematics class. // #include "ShowerKinematics.h" #include "ThePEG/Utilities/DescribeClass.h" using namespace Herwig; DescribeAbstractNoPIOClass describeShowerKinematics("Herwig::ShowerKinematics","Herwig.so"); void ShowerKinematics::updateChildren(const tShowerParticlePtr, const ShowerParticleVector &, - unsigned int pTscheme, + unsigned int , ShowerPartnerType) const { throw Exception() << "Base class ShowerKinematics::updateChildren called," << " should have been overriden in an inheriting class" << Exception::runerror; } void ShowerKinematics::resetChildren(const tShowerParticlePtr, const ShowerParticleVector &) const { throw Exception() << "Base class ShowerKinematics::resetChildren called," << " should have been overriden in an inheriting class" << Exception::runerror; } void ShowerKinematics::updateParent(const tShowerParticlePtr, const ShowerParticleVector &, unsigned int , ShowerPartnerType) const { throw Exception() << "Base class ShowerKinematics::updateParent called," << " should have been overriden in an inheriting class" << Exception::runerror; } void ShowerKinematics::reconstructChildren(const tShowerParticlePtr, const ShowerParticleVector &) const { throw Exception() << "Base class ShowerKinematics::reconstructChildren called," << " should have been overriden in an inheriting class" << Exception::runerror; } void ShowerKinematics::reconstructParent(const tShowerParticlePtr, const ParticleVector &) const { throw Exception() << "Base class ShowerKinematics::reconstructParent called," << " should have been overriden in an inheriting class" << Exception::runerror; } void ShowerKinematics::reconstructLast(const tShowerParticlePtr, Energy) const { throw Exception() << "Base class ShowerKinematics::reconstructLast called," << " should have been overriden in an inheriting class" << Exception::runerror; } void ShowerKinematics::updateLast(const tShowerParticlePtr, Energy,Energy) const { throw Exception() << "Base class ShowerKinematics::updatetLast called," << " should have been overriden in an inheriting class" << Exception::runerror; } diff --git a/Shower/QTilde/QTildeShowerHandler.cc b/Shower/QTilde/QTildeShowerHandler.cc --- a/Shower/QTilde/QTildeShowerHandler.cc +++ b/Shower/QTilde/QTildeShowerHandler.cc @@ -1,3253 +1,3251 @@ // -*- C++ -*- // // QTildeShowerHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig 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 QTildeShowerHandler class. // #include "QTildeShowerHandler.h" #include "ThePEG/Interface/Deleted.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Utilities/EnumIO.h" #include "Herwig/Shower/QTilde/Base/ShowerParticle.h" #include "Herwig/PDF/MPIPDF.h" #include "Herwig/PDF/MinBiasPDF.h" #include "Herwig/Shower/QTilde/Base/ShowerTree.h" #include "Herwig/Shower/QTilde/Base/HardTree.h" #include "Herwig/Shower/QTilde/Kinematics/KinematicsReconstructor.h" #include "Herwig/Shower/QTilde/Base/PartnerFinder.h" #include "Herwig/PDF/HwRemDecayer.h" #include "Herwig/Shower/QTilde/Base/ShowerVertex.h" #include "ThePEG/Repository/CurrentGenerator.h" #include "Herwig/MatrixElement/Matchbox/Base/SubtractedME.h" #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h" #include "ThePEG/PDF/PartonExtractor.h" #include "Herwig/Shower/RealEmissionProcess.h" #include "Herwig/Shower/QTilde/Kinematics/FS_QTildeShowerKinematics1to2.h" #include "Herwig/Shower/QTilde/Kinematics/IS_QTildeShowerKinematics1to2.h" #include "Herwig/Shower/QTilde/Kinematics/Decay_QTildeShowerKinematics1to2.h" #include "Herwig/MatrixElement/MEMinBias.h" using namespace Herwig; bool QTildeShowerHandler::_hardEmissionWarn = true; bool QTildeShowerHandler::_missingTruncWarn = true; QTildeShowerHandler::QTildeShowerHandler() : _maxtry(100), _meCorrMode(1), _evolutionScheme(1), _hardVetoReadOption(false), _iptrms(ZERO), _beta(0.), _gamma(ZERO), _iptmax(), _limitEmissions(0), _initialenhance(1.), _finalenhance(1.), _nReWeight(100), _reWeight(false), interaction_(ShowerInteraction::ALL), _trunc_Mode(true), _hardEmission(1), _softOpt(2), _hardPOWHEG(false), muPt(ZERO) {} -QTildeShowerHandler::~QTildeShowerHandler() {} - IBPtr QTildeShowerHandler::clone() const { return new_ptr(*this); } IBPtr QTildeShowerHandler::fullclone() const { return new_ptr(*this); } void QTildeShowerHandler::persistentOutput(PersistentOStream & os) const { os << _splittingGenerator << _maxtry << _meCorrMode << _hardVetoReadOption << _limitEmissions << _softOpt << _hardPOWHEG << ounit(_iptrms,GeV) << _beta << ounit(_gamma,GeV) << ounit(_iptmax,GeV) << _vetoes << _fullShowerVetoes << _nReWeight << _reWeight << _trunc_Mode << _hardEmission << _evolutionScheme << ounit(muPt,GeV) << oenum(interaction_) << _reconstructor << _partnerfinder; } void QTildeShowerHandler::persistentInput(PersistentIStream & is, int) { is >> _splittingGenerator >> _maxtry >> _meCorrMode >> _hardVetoReadOption >> _limitEmissions >> _softOpt >> _hardPOWHEG >> iunit(_iptrms,GeV) >> _beta >> iunit(_gamma,GeV) >> iunit(_iptmax,GeV) >> _vetoes >> _fullShowerVetoes >> _nReWeight >> _reWeight >> _trunc_Mode >> _hardEmission >> _evolutionScheme >> iunit(muPt,GeV) >> ienum(interaction_) >> _reconstructor >> _partnerfinder; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigQTildeShowerHandler("Herwig::QTildeShowerHandler", "HwShower.so"); void QTildeShowerHandler::Init() { static ClassDocumentation documentation ("TheQTildeShowerHandler class is the main class" " for the angular-ordered parton shower", "The Shower evolution was performed using an algorithm described in " "\\cite{Bahr:2008pv,Marchesini:1983bm,Marchesini:1987cf,Gieseke:2003rz,Masouminia:2021kne}.", // "%\\cite{Marchesini:1983bm}\n" "\\bibitem{Marchesini:1983bm}\n" " G.~Marchesini and B.~R.~Webber,\n" " ``Simulation Of QCD Jets Including Soft Gluon Interference,''\n" " Nucl.\\ Phys.\\ B {\\bf 238}, 1 (1984).\n" " %%CITATION = NUPHA,B238,1;%%\n" // "%\\cite{Marchesini:1987cf}\n" "\\bibitem{Marchesini:1987cf}\n" " G.~Marchesini and B.~R.~Webber,\n" " ``Monte Carlo Simulation of General Hard Processes with Coherent QCD\n" " Radiation,''\n" " Nucl.\\ Phys.\\ B {\\bf 310}, 461 (1988).\n" " %%CITATION = NUPHA,B310,461;%%\n" // "%\\cite{Gieseke:2003rz}\n" "\\bibitem{Gieseke:2003rz}\n" " S.~Gieseke, P.~Stephens and B.~Webber,\n" " ``New formalism for QCD parton showers,''\n" " JHEP {\\bf 0312}, 045 (2003)\n" " [arXiv:hep-ph/0310083].\n" " %%CITATION = JHEPA,0312,045;%%\n" // "%\\cite{Masouminia:2021kne}\n" "\\bibitem{Masouminia:2021kne}\n" " M.~R.~Masouminia and P.~Richardson,\n" " ``Implementation of angularly ordered electroweak parton shower in Herwig 7,''\n" " JHEP {\\bf 04}, 112 (2022)\n" " [arXiv:2108.10817 [hep-ph].\n" " %%CITATION = JHEP,04,112;%%\n" ); static Reference interfaceSplitGen("SplittingGenerator", "A reference to the SplittingGenerator object", &Herwig::QTildeShowerHandler::_splittingGenerator, false, false, true, false); static Parameter interfaceMaxTry ("MaxTry", "The maximum number of attempts to generate the shower from a" " particular ShowerTree", &QTildeShowerHandler::_maxtry, 100, 1, 100000, false, false, Interface::limited); static Parameter interfaceNReWeight ("NReWeight", "The number of attempts for the shower when reweighting", &QTildeShowerHandler::_nReWeight, 100, 10, 10000, false, false, Interface::limited); static Switch ifaceMECorrMode ("MECorrMode", "Choice of the ME Correction Mode", &QTildeShowerHandler::_meCorrMode, 1, false, false); static SwitchOption on (ifaceMECorrMode,"HardPlusSoft","hard+soft on", 1); static SwitchOption hard (ifaceMECorrMode,"Hard","only hard on", 2); static SwitchOption soft (ifaceMECorrMode,"Soft","only soft on", 3); static Switch ifaceHardVetoReadOption ("HardVetoReadOption", "Apply read-in scale veto to all collisions or just the primary one?", &QTildeShowerHandler::_hardVetoReadOption, false, false, false); static SwitchOption AllCollisions (ifaceHardVetoReadOption, "AllCollisions", "Read-in pT veto applied to primary and secondary collisions.", false); static SwitchOption PrimaryCollision (ifaceHardVetoReadOption, "PrimaryCollision", "Read-in pT veto applied to primary but not secondary collisions.", true); static Parameter ifaceiptrms ("IntrinsicPtGaussian", "RMS of intrinsic pT of Gaussian distribution:\n" "2*(1-Beta)*exp(-sqr(intrinsicpT/RMS))/sqr(RMS)", &QTildeShowerHandler::_iptrms, GeV, ZERO, ZERO, 1000000.0*GeV, false, false, Interface::limited); static Parameter ifacebeta ("IntrinsicPtBeta", "Proportion of inverse quadratic distribution in generating intrinsic pT.\n" "(1-Beta) is the proportion of Gaussian distribution", &QTildeShowerHandler::_beta, 0, 0, 1, false, false, Interface::limited); static Parameter ifacegamma ("IntrinsicPtGamma", "Parameter for inverse quadratic:\n" "2*Beta*Gamma/(sqr(Gamma)+sqr(intrinsicpT))", &QTildeShowerHandler::_gamma,GeV, ZERO, ZERO, 100000.0*GeV, false, false, Interface::limited); static Parameter ifaceiptmax ("IntrinsicPtIptmax", "Upper bound on intrinsic pT for inverse quadratic", &QTildeShowerHandler::_iptmax,GeV, ZERO, ZERO, 100000.0*GeV, false, false, Interface::limited); static RefVector ifaceVetoes ("Vetoes", "The vetoes to be checked during showering", &QTildeShowerHandler::_vetoes, -1, false,false,true,true,false); static RefVector interfaceFullShowerVetoes ("FullShowerVetoes", "The vetos to be appliede on the full final state of the shower", &QTildeShowerHandler::_fullShowerVetoes, -1, false, false, true, false, false); static Switch interfaceLimitEmissions ("LimitEmissions", "Limit the number and type of emissions for testing", &QTildeShowerHandler::_limitEmissions, 0, false, false); static SwitchOption interfaceLimitEmissionsNoLimit (interfaceLimitEmissions, "NoLimit", "Allow an arbitrary number of emissions", 0); static SwitchOption interfaceLimitEmissionsOneInitialStateEmission (interfaceLimitEmissions, "OneInitialStateEmission", "Allow one emission in the initial state and none in the final state", 1); static SwitchOption interfaceLimitEmissionsOneFinalStateEmission (interfaceLimitEmissions, "OneFinalStateEmission", "Allow one emission in the final state and none in the initial state", 2); static SwitchOption interfaceLimitEmissionsHardOnly (interfaceLimitEmissions, "HardOnly", "Only allow radiation from the hard ME correction", 3); static SwitchOption interfaceLimitEmissionsOneEmission (interfaceLimitEmissions, "OneEmission", "Allow one emission in either the final state or initial state, but not both", 4); static Switch interfaceTruncMode ("TruncatedShower", "Include the truncated shower?", &QTildeShowerHandler::_trunc_Mode, 1, false, false); static SwitchOption interfaceTruncMode0 (interfaceTruncMode,"No","Truncated Shower is OFF", 0); static SwitchOption interfaceTruncMode1 (interfaceTruncMode,"Yes","Truncated Shower is ON", 1); static Switch interfaceHardEmission ("HardEmission", "Whether to use ME corrections or POWHEG for the hardest emission", &QTildeShowerHandler::_hardEmission, 0, false, false); static SwitchOption interfaceHardEmissionNone (interfaceHardEmission, "None", "No Corrections", 0); static SwitchOption interfaceHardEmissionMECorrection (interfaceHardEmission, "MECorrection", "Old fashioned ME correction", 1); static SwitchOption interfaceHardEmissionPOWHEG (interfaceHardEmission, "POWHEG", "Powheg style hard emission", 2); static Switch interfaceInteractions ("Interactions", "The interactions to be used in the shower", &QTildeShowerHandler::interaction_, ShowerInteraction::ALL, false, false); static SwitchOption interfaceInteractionsQCD (interfaceInteractions, "QCD", "Only QCD radiation", ShowerInteraction::QCD); static SwitchOption interfaceInteractionsQED (interfaceInteractions, "QED", "Only QEd radiation", ShowerInteraction::QED); static SwitchOption interfaceInteractionEWOnly (interfaceInteractions, "EWOnly", "Only EW", ShowerInteraction::EW); static SwitchOption interfaceInteractionQEDQCD (interfaceInteractions, "QEDQCD", "QED and QCD", ShowerInteraction::QEDQCD); static SwitchOption interfaceInteractionALL (interfaceInteractions, "ALL", "QED, QCD and EW", ShowerInteraction::ALL); static Deleted delReconstructionOption ("ReconstructionOption", "The old reconstruction option switch has been replaced with" " the new EvolutionScheme switch, see arXiv:1904.11866 for details"); static Switch interfaceEvolutionScheme ("EvolutionScheme", "The scheme to interpret the evolution variable in the case of multple emission.", &QTildeShowerHandler::_evolutionScheme, 1, false, false); static SwitchOption interfaceEvolutionSchemepT (interfaceEvolutionScheme, "pT", "pT scheme", 0); static SwitchOption interfaceEvolutionSchemeDotProduct (interfaceEvolutionScheme, "DotProduct", "Dot-product scheme", 1); static SwitchOption interfaceEvolutionSchemeQ2 (interfaceEvolutionScheme, "Q2", "Q2 scheme", 2); static Switch interfaceSoftCorrelations ("SoftCorrelations", "Option for the treatment of soft correlations in the parton shower", &QTildeShowerHandler::_softOpt, 2, false, false); static SwitchOption interfaceSoftCorrelationsNone (interfaceSoftCorrelations, "No", "No soft correlations", 0); static SwitchOption interfaceSoftCorrelationsFull (interfaceSoftCorrelations, "Full", "Use the full eikonal", 1); static SwitchOption interfaceSoftCorrelationsSingular (interfaceSoftCorrelations, "Singular", "Use original Webber-Marchisini form", 2); static Switch interfaceHardPOWHEG ("HardPOWHEG", "Treatment of powheg emissions which are too hard to have a shower interpretation", &QTildeShowerHandler::_hardPOWHEG, false, false, false); static SwitchOption interfaceHardPOWHEGAsShower (interfaceHardPOWHEG, "AsShower", "Still interpret as shower emissions", false); static SwitchOption interfaceHardPOWHEGRealEmission (interfaceHardPOWHEG, "RealEmission", "Generate shower from the real emmission configuration", true); static Reference interfaceKinematicsReconstructor ("KinematicsReconstructor", "Reference to the KinematicsReconstructor object", &QTildeShowerHandler::_reconstructor, false, false, true, false, false); static Reference interfacePartnerFinder ("PartnerFinder", "Reference to the PartnerFinder object", &QTildeShowerHandler::_partnerfinder, false, false, true, false, false); } tPPair QTildeShowerHandler::cascade(tSubProPtr sub, XCPtr xcomb) { // use me for reference in tex file etc useMe(); prepareCascade(sub); // set things up in the base class resetWeights(); hard_=ShowerTreePtr(); decay_.clear(); done_.clear(); // start of the try block for the whole showering process unsigned int countFailures=0; while (countFailuresoutgoing().begin(), currentSubProcess()->outgoing().end()), hard,decay); ShowerTree::constructTrees(hard_,decay_,hard,decay); // if no hard process if(!hard_) throw Exception() << "Shower starting with a decay" << "is not implemented" << Exception::runerror; // perform the shower for the hard process showerHardProcess(hard_,xcomb); done_.push_back(hard_); hard_->updateAfterShower(decay_); // if no decaying particles to shower break out of the loop if(decay_.empty()) break; // shower the decay products while(!decay_.empty()) { // find particle whose production process has been showered ShowerDecayMap::iterator dit = decay_.begin(); while(!dit->second->parent()->hasShowered() && dit!=decay_.end()) ++dit; assert(dit!=decay_.end()); // get the particle ShowerTreePtr decayingTree = dit->second; // remove it from the multimap decay_.erase(dit); // make sure the particle has been decayed QTildeShowerHandler::decay(decayingTree,decay_); // now shower the decay showerDecay(decayingTree); done_.push_back(decayingTree); decayingTree->updateAfterShower(decay_); } // suceeded break out of the loop break; } catch (KinematicsReconstructionVeto) { resetWeights(); ++countFailures; } catch ( ... ) { hard_=ShowerTreePtr(); decay_.clear(); done_.clear(); throw; } } // if loop exited because of too many tries, throw event away if (countFailures >= maxtry()) { resetWeights(); hard_=ShowerTreePtr(); decay_.clear(); done_.clear(); throw Exception() << "Too many tries for main while loop " << "in QTildeShowerHandler::cascade()." << Exception::eventerror; } //enter the particles in the event record fillEventRecord(); // clear storage hard_=ShowerTreePtr(); decay_.clear(); done_.clear(); // non hadronic case return if (!isResolvedHadron(incomingBeams().first ) && !isResolvedHadron(incomingBeams().second) ) return incomingBeams(); // remake the remnants (needs to be after the colours are sorted // out in the insertion into the event record) if ( firstInteraction() ) return remakeRemnant(sub->incoming()); //Return the new pair of incoming partons. remakeRemnant is not //necessary here, because the secondary interactions are not yet //connected to the remnants. return make_pair(findFirstParton(sub->incoming().first ), findFirstParton(sub->incoming().second)); } void QTildeShowerHandler::fillEventRecord() { // create a new step StepPtr pstep = newStep(); assert(!done_.empty()); assert(done_[0]->isHard()); // insert the steps for(unsigned int ix=0;ixfillEventRecord(pstep,doISR(),doFSR()); } } HardTreePtr QTildeShowerHandler::generateCKKW(ShowerTreePtr ) const { return HardTreePtr(); } void QTildeShowerHandler::doinit() { ShowerHandler::doinit(); // check on the reweighting for(unsigned int ix=0;ix<_fullShowerVetoes.size();++ix) { if(_fullShowerVetoes[ix]->behaviour()==1) { _reWeight = true; break; } } if(_reWeight && maximumTries()<_nReWeight) { throw Exception() << "Reweight being performed in the shower but the number of attempts for the" << "shower is less than that for the reweighting.\n" << "Maximum number of attempt for the shower " << fullName() << ":MaxTry is " << maximumTries() << "\nand for reweighting is " << fullName() << ":NReWeight is " << _nReWeight << "\n" << "we recommend the number of attempts is 10 times the number for reweighting\n" << Exception::runerror; } ShowerTree::_vmin2 = vMin(); ShowerTree::_spaceTime = includeSpaceTime(); } void QTildeShowerHandler::doinitrun() { ShowerHandler::doinitrun(); ShowerTree::_vmin2 = vMin(); ShowerTree::_spaceTime = includeSpaceTime(); } void QTildeShowerHandler::generateIntrinsicpT(vector particlesToShower) { if ( !ipTon() || !doISR() ) return; // don't do anything for the moment for secondary scatters if( !firstInteraction() ) return; // generate intrinsic pT for(unsigned int ix=0;ixprogenitor()->isFinalState()) continue; if(!particlesToShower[ix]->progenitor()->dataPtr()->coloured()) continue; Energy ipt; if(UseRandom::rnd() > _beta) { ipt=_iptrms*sqrt(-log(UseRandom::rnd())); } else { ipt=_gamma*sqrt(pow(1.+sqr(_iptmax/_gamma), UseRandom::rnd())-1.); } pair pt = make_pair(ipt,UseRandom::rnd(Constants::twopi)); _intrinsic[particlesToShower[ix]] = pt; } } void QTildeShowerHandler::setupMaximumScales(const vector & p, XCPtr xcomb) { // let POWHEG events radiate freely if(_hardEmission==2&&hardTree()) { vector::const_iterator ckt = p.begin(); for (; ckt != p.end(); ckt++) (*ckt)->maxHardPt(Constants::MaxEnergy); return; } // return if no vetos if (!restrictPhasespace()) return; // find out if hard partonic subprocess. bool isPartonic(false); map::const_iterator cit = _currenttree->incomingLines().begin(); Lorentz5Momentum pcm; for(; cit!=currentTree()->incomingLines().end(); ++cit) { pcm += cit->first->progenitor()->momentum(); isPartonic |= cit->first->progenitor()->coloured(); } // find minimum pt from hard process, the maximum pt from all outgoing // coloured lines (this is simpler and more general than // 2stu/(s^2+t^2+u^2)). Maximum scale for scattering processes will // be transverse mass. Energy ptmax = generator()->maximumCMEnergy(); // general case calculate the scale if ( !hardScaleIsMuF() || (hardVetoReadOption()&&!firstInteraction()) ) { // scattering process if(currentTree()->isHard()) { assert(xcomb); // coloured incoming particles if (isPartonic) { map::const_iterator cjt = currentTree()->outgoingLines().begin(); for(; cjt!=currentTree()->outgoingLines().end(); ++cjt) { if (cjt->first->progenitor()->coloured()) ptmax = min(ptmax,cjt->first->progenitor()->momentum().mt()); } } if (ptmax == generator()->maximumCMEnergy() ) ptmax = pcm.m(); if(hardScaleIsMuF()&&hardVetoReadOption()&& !firstInteraction()) { ptmax=min(ptmax,sqrt(xcomb->lastShowerScale())); } } // decay, incoming() is the decaying particle. else { ptmax = currentTree()->incomingLines().begin()->first ->progenitor()->momentum().mass(); } } // hepeup.SCALUP is written into the lastXComb by the // LesHouchesReader itself - use this by user's choice. // Can be more general than this. else { if(currentTree()->isHard()) { assert(xcomb); ptmax = sqrt( xcomb->lastShowerScale() ); } else { ptmax = currentTree()->incomingLines().begin()->first ->progenitor()->momentum().mass(); } } ptmax *= hardScaleFactor(); // set maxHardPt for all progenitors. For partonic processes this // is now the max pt in the FS, for non-partonic processes or // processes with no coloured FS the invariant mass of the IS vector::const_iterator ckt = p.begin(); for (; ckt != p.end(); ckt++) (*ckt)->maxHardPt(ptmax); } void QTildeShowerHandler::setupHardScales(const vector & p, XCPtr xcomb) { if ( hardScaleIsMuF() && (!hardVetoReadOption() || firstInteraction()) ) { Energy hardScale = ZERO; if(currentTree()->isHard()) { assert(xcomb); hardScale = sqrt( xcomb->lastShowerScale() ); } else { hardScale = currentTree()->incomingLines().begin()->first ->progenitor()->momentum().mass(); } hardScale *= hardScaleFactor(); vector::const_iterator ckt = p.begin(); for (; ckt != p.end(); ckt++) (*ckt)->hardScale(hardScale); muPt = hardScale; } } void QTildeShowerHandler::showerHardProcess(ShowerTreePtr hard, XCPtr xcomb) { _hardme = HwMEBasePtr(); // extract the matrix element tStdXCombPtr lastXC = dynamic_ptr_cast(xcomb); if(lastXC) { _hardme = dynamic_ptr_cast(lastXC->matrixElement()); } _decayme = HwDecayerBasePtr(); // set the current tree currentTree(hard); hardTree(HardTreePtr()); // work out the type of event currentTree()->xcombPtr(dynamic_ptr_cast(xcomb)); currentTree()->identifyEventType(); checkFlags(); // generate the showering doShowering(true,xcomb); } RealEmissionProcessPtr QTildeShowerHandler::hardMatrixElementCorrection(bool hard) { // set the initial enhancement factors for the soft correction _initialenhance = 1.; _finalenhance = 1.; // see if we can get the correction from the matrix element // or decayer RealEmissionProcessPtr real; if(hard) { if(_hardme&&_hardme->hasMECorrection()) { _hardme->initializeMECorrection(_currenttree->perturbativeProcess(), _initialenhance,_finalenhance); if(hardMEC()) real = _hardme->applyHardMatrixElementCorrection(_currenttree->perturbativeProcess()); } } else { if(_decayme&&_decayme->hasMECorrection()) { _decayme->initializeMECorrection(_currenttree->perturbativeProcess(), _initialenhance,_finalenhance); if(hardMEC()) real = _decayme->applyHardMatrixElementCorrection(_currenttree->perturbativeProcess()); } } return real; } ShowerParticleVector QTildeShowerHandler::createTimeLikeChildren(tShowerParticlePtr, IdList ids) { // Create the ShowerParticle objects for the two children of // the emitting particle; set the parent/child relationship // if same as definition create particles, otherwise create cc ShowerParticleVector children; for(unsigned int ix=0;ix<2;++ix) { children.push_back(new_ptr(ShowerParticle(ids[ix+1],true))); if(children[ix]->id()==_progenitor->id()&&!ids[ix+1]->stable()&&abs(ids[ix+1]->id())!=ParticleID::tauminus) children[ix]->set5Momentum(Lorentz5Momentum(_progenitor->progenitor()->mass())); else children[ix]->set5Momentum(Lorentz5Momentum(ids[ix+1]->mass())); } return children; } bool QTildeShowerHandler::timeLikeShower(tShowerParticlePtr particle, ShowerInteraction type, Branching fb, bool first) { // don't do anything if not needed if(_limitEmissions == 1 || hardOnly() || ( _limitEmissions == 2 && _nfs != 0) || ( _limitEmissions == 4 && _nfs + _nis != 0) ) { if(particle->spinInfo()) particle->spinInfo()->develop(); return false; } // generate the emission ShowerParticleVector children; // generate the emission if(!fb.kinematics) fb = selectTimeLikeBranching(particle,type,HardBranchingPtr()); // no emission, return if(!fb.kinematics) { if(particle->spinInfo()) particle->spinInfo()->develop(); return false; } Branching fc[2] = {Branching(),Branching()}; assert(fb.kinematics); // has emitted // Assign the shower kinematics to the emitting particle. particle->showerKinematics(fb.kinematics); // check highest pT if(fb.kinematics->pT()>progenitor()->highestpT()) progenitor()->highestpT(fb.kinematics->pT()); // create the children children = createTimeLikeChildren(particle,fb.ids); // update the children particle->showerKinematics()-> updateChildren(particle, children,_evolutionScheme,fb.type); // update number of emissions ++_nfs; if(_limitEmissions!=0) { if(children[0]->spinInfo()) children[0]->spinInfo()->develop(); if(children[1]->spinInfo()) children[1]->spinInfo()->develop(); if(particle->spinInfo()) particle->spinInfo()->develop(); return true; } // select branchings for children fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr()); fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr()); // shower the first particle if(fc[0].kinematics) timeLikeShower(children[0],type,fc[0],false); if(children[0]->spinInfo()) children[0]->spinInfo()->develop(); // shower the second particle if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],false); if(children[1]->spinInfo()) children[1]->spinInfo()->develop(); particle->showerKinematics()->updateParent(particle, children,_evolutionScheme,fb.type); // branching has happened if(first&&!children.empty()) particle->showerKinematics()->resetChildren(particle,children); if(particle->spinInfo()) particle->spinInfo()->develop(); return true; } bool QTildeShowerHandler::spaceLikeShower(tShowerParticlePtr particle, PPtr beam, ShowerInteraction type) { //using the pdf's associated with the ShowerHandler assures, that //modified pdf's are used for the secondary interactions via //CascadeHandler::resetPDFs(...) tcPDFPtr pdf; if(beam == incomingBeams().first) pdf = firstPDF().pdf(); if(beam == incomingBeams().second) pdf = secondPDF().pdf(); Energy freeze = pdfFreezingScale(); // don't do anything if not needed if(_limitEmissions == 2 || hardOnly() || ( _limitEmissions == 1 && _nis != 0 ) || ( _limitEmissions == 4 && _nis + _nfs != 0 ) ) { if(particle->spinInfo()) particle->spinInfo()->develop(); return false; } Branching bb; // generate branching while (true) { bb=_splittingGenerator->chooseBackwardBranching(*particle,beam, _initialenhance, _beam,type, pdf,freeze); // return if no emission if(!bb.kinematics) { if(particle->spinInfo()) particle->spinInfo()->develop(); return false; } // if not vetoed break if(!spaceLikeVetoed(bb,particle)) break; // otherwise reset scale and continue particle->vetoEmission(bb.type,bb.kinematics->scale()); if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr()); } // assign the splitting function and shower kinematics particle->showerKinematics(bb.kinematics); if(bb.kinematics->pT()>progenitor()->highestpT()) progenitor()->highestpT(bb.kinematics->pT()); // For the time being we are considering only 1->2 branching // particles as in Sudakov form factor tcPDPtr part[2]={bb.ids[0],bb.ids[2]}; // Now create the actual particles, make the otherChild a final state // particle, while the newParent is not ShowerParticlePtr newParent = new_ptr(ShowerParticle(part[0],false)); ShowerParticlePtr otherChild = new_ptr(ShowerParticle(part[1],true,true)); ShowerParticleVector theChildren; theChildren.push_back(particle); theChildren.push_back(otherChild); //this updates the evolution scale particle->showerKinematics()-> updateParent(newParent, theChildren,_evolutionScheme,bb.type); // update the history if needed _currenttree->updateInitialStateShowerProduct(_progenitor,newParent); _currenttree->addInitialStateBranching(particle,newParent,otherChild); // for the reconstruction of kinematics, parent/child // relationships are according to the branching process: // now continue the shower ++_nis; bool emitted = _limitEmissions==0 ? spaceLikeShower(newParent,beam,type) : false; if(newParent->spinInfo()) newParent->spinInfo()->develop(); // now reconstruct the momentum if(!emitted) { if(_intrinsic.find(_progenitor)==_intrinsic.end()) { bb.kinematics->updateLast(newParent,ZERO,ZERO); } else { pair kt=_intrinsic[_progenitor]; bb.kinematics->updateLast(newParent, kt.first*cos(kt.second), kt.first*sin(kt.second)); } } particle->showerKinematics()-> updateChildren(newParent, theChildren,_evolutionScheme,bb.type); if(_limitEmissions!=0) { if(particle->spinInfo()) particle->spinInfo()->develop(); return true; } // perform the shower of the final-state particle timeLikeShower(otherChild,type,Branching(),true); particle->showerKinematics()-> updateChildren(newParent, theChildren,_evolutionScheme,bb.type); updateHistory(otherChild); if(theChildren[1]->spinInfo()) theChildren[1]->spinInfo()->develop(); // return the emitted if(particle->spinInfo()) particle->spinInfo()->develop(); if(!theChildren.empty()){ particle->showerKinematics()->resetChildren(newParent,theChildren); } return true; } void QTildeShowerHandler::showerDecay(ShowerTreePtr decay) { // work out the type of event currentTree()->xcombPtr(StdXCombPtr()); currentTree()->identifyEventType(); _decayme = HwDecayerBasePtr(); _hardme = HwMEBasePtr(); // find the decayer // try the normal way if possible tDMPtr dm = decay->incomingLines().begin()->first->original() ->decayMode(); if(!dm) dm = decay->incomingLines().begin()->first->copy() ->decayMode(); if(!dm) dm = decay->incomingLines().begin()->first->progenitor()->decayMode(); // otherwise make a string and look it up if(!dm) { string tag = decay->incomingLines().begin()->first->original()->dataPtr()->name() + "->"; OrderedParticles outgoing; for(map::const_iterator it=decay->outgoingLines().begin();it!=decay->outgoingLines().end();++it) { if(abs(decay->incomingLines().begin()->first->original()->id()) == ParticleID::t && abs(it->first->original()->id())==ParticleID::Wplus && decay->treelinks().size() == 1) { ShowerTreePtr Wtree = decay->treelinks().begin()->first; for(map::const_iterator it2=Wtree->outgoingLines().begin();it2!=Wtree->outgoingLines().end();++it2) { outgoing.insert(it2->first->original()->dataPtr()); } } else { outgoing.insert(it->first->original()->dataPtr()); } } for(OrderedParticles::const_iterator it=outgoing.begin(); it!=outgoing.end();++it) { if(it!=outgoing.begin()) tag += ","; tag +=(**it).name(); } tag += ";"; dm = findDecayMode(tag); } if(dm) _decayme = dynamic_ptr_cast(dm->decayer()); // set the ShowerTree to be showered currentTree(decay); decay->applyTransforms(); hardTree(HardTreePtr()); // generate the showering doShowering(false,XCPtr()); // if no vetos // force calculation of spin correlations SpinPtr spInfo = decay->incomingLines().begin()->first->progenitor()->spinInfo(); if(spInfo) { if(!spInfo->developed()) spInfo->needsUpdate(); spInfo->develop(); } } bool QTildeShowerHandler::spaceLikeDecayShower(tShowerParticlePtr particle, const ShowerParticle::EvolutionScales & maxScales, Energy minmass,ShowerInteraction type, Branching fb) { // don't do anything if not needed if(_limitEmissions == 1 || hardOnly() || ( _limitEmissions == 3 && _nis != 0) || ( _limitEmissions == 4 && _nfs + _nis != 0) ) { return false; } // generate the emission ShowerParticleVector children; // generate the emission if(!fb.kinematics) fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type, HardBranchingPtr()); // no emission, return if(!fb.kinematics) return false; Branching fc[2]; if(particle->virtualMass()==ZERO) particle->virtualMass(_progenitor->progenitor()->mass()); fc[0] = Branching(); fc[1] = Branching(); assert(fb.kinematics); // has emitted // Assign the shower kinematics to the emitting particle. particle->showerKinematics(fb.kinematics); if(fb.kinematics->pT()>progenitor()->highestpT()) progenitor()->highestpT(fb.kinematics->pT()); // create the ShowerParticle objects for the two children children = createTimeLikeChildren(particle,fb.ids); // updateChildren the children particle->showerKinematics()-> updateChildren(particle, children,_evolutionScheme,fb.type); // select branchings for children fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass, type,HardBranchingPtr()); fc[1] = selectTimeLikeBranching (children[1],type,HardBranchingPtr()); // old default ++_nis; // shower the first particle _currenttree->updateInitialStateShowerProduct(_progenitor,children[0]); _currenttree->addInitialStateBranching(particle,children[0],children[1]); if(fc[0].kinematics) spaceLikeDecayShower(children[0],maxScales,minmass,type,Branching()); // shower the second particle if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],true); updateHistory(children[1]); // TODO NEED AN UPDATE HERE FOR RECONOPT!=0 // branching has happened return true; } vector QTildeShowerHandler::setupShower(bool hard) { RealEmissionProcessPtr real; // generate hard me if needed if(_hardEmission==1) { real = hardMatrixElementCorrection(hard); if(real&&!real->outgoing().empty()) setupMECorrection(real); } // generate POWHEG hard emission if needed else if(_hardEmission==2) hardestEmission(hard); // set the initial colour partners setEvolutionPartners(hard,interaction_,false); // get the particles to be showered vector particlesToShower = currentTree()->extractProgenitors(); // return the answer return particlesToShower; } void QTildeShowerHandler::setEvolutionPartners(bool hard,ShowerInteraction type, bool clear) { // match the particles in the ShowerTree and hardTree if(hardTree() && !hardTree()->connect(currentTree())) throw Exception() << "Can't match trees in " << "QTildeShowerHandler::setEvolutionPartners()" << Exception::eventerror; // extract the progenitors vector particles = currentTree()->extractProgenitorParticles(); // clear the partners if needed if(clear) { for(unsigned int ix=0;ixpartner(ShowerParticlePtr()); particles[ix]->clearPartners(); } } // sort out the colour partners if(hardTree()) { // find the partner for(unsigned int ix=0;ixparticles()[particles[ix]]->branchingParticle()->partner(); if(!partner) continue; for(map::const_iterator it=hardTree()->particles().begin(); it!=hardTree()->particles().end();++it) { if(it->second->branchingParticle()==partner) { particles[ix]->partner(it->first); break; } } if(!particles[ix]->partner()) throw Exception() << "Can't match partners in " << "QTildeShowerHandler::setEvolutionPartners()" << Exception::eventerror; } } // Set the initial evolution scales partnerFinder()-> setInitialEvolutionScales(particles,!hard,interaction_,!_hardtree); if(hardTree() && _hardPOWHEG) { bool tooHard=false; map::const_iterator eit=hardTree()->particles().end(); for(unsigned int ix=0;ix::const_iterator mit = hardTree()->particles().find(particles[ix]); Energy hardScale(ZERO); ShowerPartnerType type(ShowerPartnerType::Undefined); // final-state if(particles[ix]->isFinalState()) { if(mit!= eit && !mit->second->children().empty()) { hardScale = mit->second->scale(); type = mit->second->type(); } } // initial-state else { if(mit!= eit && mit->second->parent()) { hardScale = mit->second->parent()->scale(); type = mit->second->parent()->type(); } } if(type!=ShowerPartnerType::Undefined) { if(type==ShowerPartnerType::QED) { tooHard |= particles[ix]->scales().QED_noAOscales().QCD_c_noAOscales().QCD_ac_noAOscales().EWchildren().empty()) { ShowerParticleVector theChildren; for(unsigned int ix=0;ixchildren().size();++ix) { ShowerParticlePtr part = dynamic_ptr_cast (particle->children()[ix]); theChildren.push_back(part); } // update the history if needed if(particle==_currenttree->getFinalStateShowerProduct(_progenitor)) _currenttree->updateFinalStateShowerProduct(_progenitor, particle,theChildren); _currenttree->addFinalStateBranching(particle,theChildren); for(unsigned int ix=0;ixprogenitor()->partner()) return false; progenitor()->progenitor()->initializeFinalState(); if(hardTree()) { map::const_iterator eit=hardTree()->particles().end(), mit = hardTree()->particles().find(progenitor()->progenitor()); if( mit != eit && !mit->second->children().empty() ) { bool output=truncatedTimeLikeShower(progenitor()->progenitor(), mit->second ,type,Branching(),true); if(output) updateHistory(progenitor()->progenitor()); return output; } } // do the shower bool output = hardOnly() ? false : timeLikeShower(progenitor()->progenitor() ,type,Branching(),true) ; if(output) updateHistory(progenitor()->progenitor()); return output; } bool QTildeShowerHandler::startSpaceLikeShower(PPtr parent, ShowerInteraction type) { // initialise the basis vectors if(!progenitor()->progenitor()->partner()) return false; progenitor()->progenitor()->initializeInitialState(parent); if(hardTree()) { map::const_iterator eit =hardTree()->particles().end(), mit = hardTree()->particles().find(progenitor()->progenitor()); if( mit != eit && mit->second->parent() ) { return truncatedSpaceLikeShower( progenitor()->progenitor(), parent, mit->second->parent(), type ); } } // perform the shower return hardOnly() ? false : spaceLikeShower(progenitor()->progenitor(),parent,type); } bool QTildeShowerHandler:: startSpaceLikeDecayShower(const ShowerParticle::EvolutionScales & maxScales, Energy minimumMass,ShowerInteraction type) { // set up the particle basis vectors if(!progenitor()->progenitor()->partner()) return false; progenitor()->progenitor()->initializeDecay(); if(hardTree()) { map::const_iterator eit =hardTree()->particles().end(), mit = hardTree()->particles().find(progenitor()->progenitor()); if( mit != eit && mit->second->parent() ) { HardBranchingPtr branch=mit->second; while(branch->parent()) branch=branch->parent(); return truncatedSpaceLikeDecayShower(progenitor()->progenitor(),maxScales, minimumMass, branch ,type, Branching()); } } // perform the shower return hardOnly() ? false : spaceLikeDecayShower(progenitor()->progenitor(),maxScales,minimumMass,type,Branching()); } bool QTildeShowerHandler::timeLikeVetoed(const Branching & fb, ShowerParticlePtr particle) { // work out type of interaction ShowerInteraction type = convertInteraction(fb.type); // check whether emission was harder than largest pt of hard subprocess if ( restrictPhasespace() && fb.kinematics->pT() > _progenitor->maxHardPt() ) return true; // soft matrix element correction veto if( softMEC()) { if(_hardme && _hardme->hasMECorrection()) { if(_hardme->softMatrixElementVeto(particle, _progenitor->progenitor(), particle->isFinalState(), _progenitor->highestpT(), fb.ids, fb.kinematics->z(), fb.kinematics->scale(), fb.kinematics->pT())) return true; } else if(_decayme && _decayme->hasMECorrection()) { if(_decayme->softMatrixElementVeto(particle, _progenitor->progenitor(), particle->isFinalState(), _progenitor->highestpT(), fb.ids, fb.kinematics->z(), fb.kinematics->scale(), fb.kinematics->pT())) return true; } } // veto on maximum pt if(fb.kinematics->pT()>_progenitor->maximumpT(type)) return true; // general vetos if (fb.kinematics && !_vetoes.empty()) { bool vetoed=false; for (vector::iterator v = _vetoes.begin(); v != _vetoes.end(); ++v) { bool test = (**v).vetoTimeLike(_progenitor,particle,fb,currentTree()); switch((**v).vetoType()) { case ShowerVeto::Emission: vetoed |= test; break; case ShowerVeto::Shower: if(test) throw VetoShower(); break; case ShowerVeto::Event: if(test) throw Veto(); break; } } if(vetoed) return true; } if ( firstInteraction() && profileScales() ) { double weight = profileScales()-> hardScaleProfile(_progenitor->hardScale(),fb.kinematics->pT()); if ( UseRandom::rnd() > weight ) return true; } return false; } bool QTildeShowerHandler::spaceLikeVetoed(const Branching & bb, ShowerParticlePtr particle) { // work out type of interaction ShowerInteraction type = convertInteraction(bb.type); // check whether emission was harder than largest pt of hard subprocess if (restrictPhasespace() && bb.kinematics->pT() > _progenitor->maxHardPt()) return true; // apply the soft correction if( softMEC() && _hardme && _hardme->hasMECorrection() ) { if(_hardme->softMatrixElementVeto(particle, _progenitor->progenitor(), particle->isFinalState(), _progenitor->highestpT(), bb.ids, bb.kinematics->z(), bb.kinematics->scale(), bb.kinematics->pT())) return true; } // the more general vetos // check vs max pt for the shower if(bb.kinematics->pT()>_progenitor->maximumpT(type)) return true; if (!_vetoes.empty()) { bool vetoed=false; for (vector::iterator v = _vetoes.begin(); v != _vetoes.end(); ++v) { bool test = (**v).vetoSpaceLike(_progenitor,particle,bb,currentTree()); switch ((**v).vetoType()) { case ShowerVeto::Emission: vetoed |= test; break; case ShowerVeto::Shower: if(test) throw VetoShower(); break; case ShowerVeto::Event: if(test) throw Veto(); break; } } if (vetoed) return true; } if ( firstInteraction() && profileScales() ) { double weight = profileScales()-> hardScaleProfile(_progenitor->hardScale(),bb.kinematics->pT()); if ( UseRandom::rnd() > weight ) return true; } return false; } bool QTildeShowerHandler::spaceLikeDecayVetoed( const Branching & fb, ShowerParticlePtr particle) { // work out type of interaction ShowerInteraction type = convertInteraction(fb.type); // apply the soft correction if( softMEC() && _decayme && _decayme->hasMECorrection() ) { if(_decayme->softMatrixElementVeto(particle, _progenitor->progenitor(), particle->isFinalState(), _progenitor->highestpT(), fb.ids, fb.kinematics->z(), fb.kinematics->scale(), fb.kinematics->pT())) return true; } // veto on hardest pt in the shower if(fb.kinematics->pT()> _progenitor->maximumpT(type)) return true; // general vetos if (!_vetoes.empty()) { bool vetoed=false; for (vector::iterator v = _vetoes.begin(); v != _vetoes.end(); ++v) { bool test = (**v).vetoSpaceLike(_progenitor,particle,fb,currentTree()); switch((**v).vetoType()) { case ShowerVeto::Emission: vetoed |= test; break; case ShowerVeto::Shower: if(test) throw VetoShower(); break; case ShowerVeto::Event: if(test) throw Veto(); break; } if (vetoed) return true; } } return false; } void QTildeShowerHandler::hardestEmission(bool hard) { HardTreePtr ISRTree; // internal POWHEG in production or decay if( (( _hardme && _hardme->hasPOWHEGCorrection()!=0 ) || ( _decayme && _decayme->hasPOWHEGCorrection()!=0 ) ) ) { RealEmissionProcessPtr real; unsigned int type(0); // production if(_hardme) { assert(hard); real = _hardme->generateHardest( currentTree()->perturbativeProcess(), interaction_); type = _hardme->hasPOWHEGCorrection(); } // decay else { assert(!hard); real = _decayme->generateHardest( currentTree()->perturbativeProcess() ); type = _decayme->hasPOWHEGCorrection(); } if(real) { // set up ther hard tree if(!real->outgoing().empty()) _hardtree = new_ptr(HardTree(real)); // set up the vetos currentTree()->setVetoes(real->pT(),type); } // store initial state POWHEG radiation if(_hardtree && _hardme && _hardme->hasPOWHEGCorrection()==1) ISRTree = _hardtree; } else if (hard) { // Get minimum pT cutoff used in shower approximation Energy maxpt = 1.*GeV; if ( currentTree()->showerApproximation() ) { int colouredIn = 0; int colouredOut = 0; for( map< ShowerProgenitorPtr, tShowerParticlePtr >::iterator it = currentTree()->outgoingLines().begin(); it != currentTree()->outgoingLines().end(); ++it ) { if( it->second->coloured() ) ++colouredOut; } for( map< ShowerProgenitorPtr, ShowerParticlePtr >::iterator it = currentTree()->incomingLines().begin(); it != currentTree()->incomingLines().end(); ++it ) { if( it->second->coloured() ) ++colouredIn; } if ( currentTree()->showerApproximation()->ffPtCut() == currentTree()->showerApproximation()->fiPtCut() && currentTree()->showerApproximation()->ffPtCut() == currentTree()->showerApproximation()->iiPtCut() ) maxpt = currentTree()->showerApproximation()->ffPtCut(); else if ( colouredIn == 2 && colouredOut == 0 ) maxpt = currentTree()->showerApproximation()->iiPtCut(); else if ( colouredIn == 0 && colouredOut > 1 ) maxpt = currentTree()->showerApproximation()->ffPtCut(); else if ( colouredIn == 2 && colouredOut == 1 ) maxpt = min(currentTree()->showerApproximation()->iiPtCut(), currentTree()->showerApproximation()->fiPtCut()); else if ( colouredIn == 1 && colouredOut > 1 ) maxpt = min(currentTree()->showerApproximation()->ffPtCut(), currentTree()->showerApproximation()->fiPtCut()); else maxpt = min(min(currentTree()->showerApproximation()->iiPtCut(), currentTree()->showerApproximation()->fiPtCut()), currentTree()->showerApproximation()->ffPtCut()); } // Generate hardtree from born and real emission subprocesses _hardtree = generateCKKW(currentTree()); // Find transverse momentum of hardest emission if (_hardtree){ for(set::iterator it=_hardtree->branchings().begin(); it!=_hardtree->branchings().end();++it) { if ((*it)->parent() && (*it)->status()==HardBranching::Incoming) maxpt=(*it)->branchingParticle()->momentum().perp(); if ((*it)->children().size()==2 && (*it)->status()==HardBranching::Outgoing){ if ((*it)->branchingParticle()->id()!=21 && abs((*it)->branchingParticle()->id())>5 ){ if ((*it)->children()[0]->branchingParticle()->id()==21 || abs((*it)->children()[0]->branchingParticle()->id())<6) maxpt=(*it)->children()[0]->branchingParticle()->momentum().perp(); else if ((*it)->children()[1]->branchingParticle()->id()==21 || abs((*it)->children()[1]->branchingParticle()->id())<6) maxpt=(*it)->children()[1]->branchingParticle()->momentum().perp(); } else { if ( abs((*it)->branchingParticle()->id())<6){ if (abs((*it)->children()[0]->branchingParticle()->id())<6) maxpt = (*it)->children()[1]->branchingParticle()->momentum().perp(); else maxpt = (*it)->children()[0]->branchingParticle()->momentum().perp(); } else maxpt = (*it)->children()[1]->branchingParticle()->momentum().perp(); } } } } // Hardest (pt) emission should be the first powheg emission. maxpt=min(sqrt(lastXCombPtr()->lastShowerScale()),maxpt); // set maximum pT for subsequent emissions from S events if ( currentTree()->isPowhegSEvent() ) { for( map< ShowerProgenitorPtr, tShowerParticlePtr >::iterator it = currentTree()->outgoingLines().begin(); it != currentTree()->outgoingLines().end(); ++it ) { if( ! it->second->coloured() ) continue; it->first->maximumpT(maxpt, ShowerInteraction::QCD ); } for( map< ShowerProgenitorPtr, ShowerParticlePtr >::iterator it = currentTree()->incomingLines().begin(); it != currentTree()->incomingLines().end(); ++it ) { if( ! it->second->coloured() ) continue; it->first->maximumpT(maxpt, ShowerInteraction::QCD ); } } } else _hardtree = generateCKKW(currentTree()); // if hard me doesn't have a FSR powheg // correction use decay powheg correction if (_hardme && _hardme->hasPOWHEGCorrection()<2) { addFSRUsingDecayPOWHEG(ISRTree); } // connect the trees if(_hardtree) { connectTrees(currentTree(),_hardtree,hard); } } void QTildeShowerHandler::addFSRUsingDecayPOWHEG(HardTreePtr ISRTree) { // check for intermediate colour singlet resonance const ParticleVector inter = _hardme->subProcess()->intermediates(); if (inter.size()!=1 || inter[0]->momentum().m2()/GeV2 < 0 || inter[0]->dataPtr()->iColour()!=PDT::Colour0) { return; } // ignore cases where outgoing particles are not coloured map out = currentTree()->outgoingLines(); if (out.size() != 2 || out. begin()->second->dataPtr()->iColour()==PDT::Colour0 || out.rbegin()->second->dataPtr()->iColour()==PDT::Colour0) { return; } // look up decay mode tDMPtr dm; string tag; string inParticle = inter[0]->dataPtr()->name() + "->"; vector outParticles; outParticles.push_back(out.begin ()->first->progenitor()->dataPtr()->name()); outParticles.push_back(out.rbegin()->first->progenitor()->dataPtr()->name()); for (int it=0; it<2; ++it){ tag = inParticle + outParticles[it] + "," + outParticles[(it+1)%2] + ";"; dm = generator()->findDecayMode(tag); if(dm) break; } // get the decayer HwDecayerBasePtr decayer; if(dm) decayer = dynamic_ptr_cast(dm->decayer()); // check if decayer has a FSR POWHEG correction if (!decayer || decayer->hasPOWHEGCorrection()<2) { return; } // generate the hardest emission // create RealEmissionProcess PPtr in = new_ptr(*inter[0]); RealEmissionProcessPtr newProcess(new_ptr(RealEmissionProcess())); newProcess->bornIncoming().push_back(in); newProcess->bornOutgoing().push_back(out.begin ()->first->progenitor()); newProcess->bornOutgoing().push_back(out.rbegin()->first->progenitor()); // generate the FSR newProcess = decayer->generateHardest(newProcess); HardTreePtr FSRTree; if(newProcess) { // set up ther hard tree if(!newProcess->outgoing().empty()) FSRTree = new_ptr(HardTree(newProcess)); // set up the vetos currentTree()->setVetoes(newProcess->pT(),2); } if(!FSRTree) return; // if there is no ISRTree make _hardtree from FSRTree if (!ISRTree){ vector inBranch,hardBranch; for(map::const_iterator cit =currentTree()->incomingLines().begin(); cit!=currentTree()->incomingLines().end();++cit ) { inBranch.push_back(new_ptr(HardBranching(cit->second,SudakovPtr(), HardBranchingPtr(), HardBranching::Incoming))); inBranch.back()->beam(cit->first->original()->parents()[0]); hardBranch.push_back(inBranch.back()); } if(inBranch[0]->branchingParticle()->dataPtr()->coloured()) { inBranch[0]->colourPartner(inBranch[1]); inBranch[1]->colourPartner(inBranch[0]); } for(set::iterator it=FSRTree->branchings().begin(); it!=FSRTree->branchings().end();++it) { if((**it).branchingParticle()->id()!=in->id()) hardBranch.push_back(*it); } hardBranch[2]->colourPartner(hardBranch[3]); hardBranch[3]->colourPartner(hardBranch[2]); HardTreePtr newTree = new_ptr(HardTree(hardBranch,inBranch, ShowerInteraction::QCD)); _hardtree = newTree; } // Otherwise modify the ISRTree to include the emission in FSRTree else { vector FSROut, ISROut; set::iterator itFSR, itISR; // get outgoing particles for(itFSR =FSRTree->branchings().begin(); itFSR!=FSRTree->branchings().end();++itFSR){ if ((**itFSR).status()==HardBranching::Outgoing) FSROut.push_back((*itFSR)->branchingParticle()); } for(itISR =ISRTree->branchings().begin(); itISR!=ISRTree->branchings().end();++itISR){ if ((**itISR).status()==HardBranching::Outgoing) ISROut.push_back((*itISR)->branchingParticle()); } // find COM frame formed by outgoing particles LorentzRotation eventFrameFSR, eventFrameISR; eventFrameFSR = ((FSROut[0]->momentum()+FSROut[1]->momentum()).findBoostToCM()); eventFrameISR = ((ISROut[0]->momentum()+ISROut[1]->momentum()).findBoostToCM()); // find rotation between ISR and FSR frames int j=0; if (ISROut[0]->id()!=FSROut[0]->id()) j=1; eventFrameISR.rotateZ( (eventFrameFSR*FSROut[0]->momentum()).phi()- (eventFrameISR*ISROut[j]->momentum()).phi() ); eventFrameISR.rotateY( (eventFrameFSR*FSROut[0]->momentum()).theta()- (eventFrameISR*ISROut[j]->momentum()).theta() ); eventFrameISR.invert(); for (itFSR=FSRTree->branchings().begin(); itFSR!=FSRTree->branchings().end();++itFSR){ if ((**itFSR).branchingParticle()->id()==in->id()) continue; for (itISR =ISRTree->branchings().begin(); itISR!=ISRTree->branchings().end();++itISR){ if ((**itISR).status()==HardBranching::Incoming) continue; if ((**itFSR).branchingParticle()->id()== (**itISR).branchingParticle()->id()){ // rotate FSRTree particle to ISRTree event frame (**itISR).branchingParticle()->setMomentum(eventFrameISR* eventFrameFSR* (**itFSR).branchingParticle()->momentum()); (**itISR).branchingParticle()->rescaleMass(); // add the children of the FSRTree particles to the ISRTree if(!(**itFSR).children().empty()){ (**itISR).addChild((**itFSR).children()[0]); (**itISR).addChild((**itFSR).children()[1]); // rotate momenta to ISRTree event frame (**itISR).children()[0]->branchingParticle()->setMomentum(eventFrameISR* eventFrameFSR* (**itFSR).children()[0]->branchingParticle()->momentum()); (**itISR).children()[1]->branchingParticle()->setMomentum(eventFrameISR* eventFrameFSR* (**itFSR).children()[1]->branchingParticle()->momentum()); } } } } _hardtree = ISRTree; } } bool QTildeShowerHandler::truncatedTimeLikeShower(tShowerParticlePtr particle, HardBranchingPtr branch, ShowerInteraction type, Branching fb, bool first) { // select a branching if we don't have one if(!fb.kinematics) fb = selectTimeLikeBranching(particle,type,branch); // must be an emission, the forced one it not a truncated one assert(fb.kinematics); ShowerParticleVector children; Branching fc[2] = {Branching(),Branching()}; // Assign the shower kinematics to the emitting particle. particle->showerKinematics(fb.kinematics); if(fb.kinematics->pT()>progenitor()->highestpT()) progenitor()->highestpT(fb.kinematics->pT()); // create the children children = createTimeLikeChildren(particle,fb.ids); // update the children particle->showerKinematics()-> updateChildren(particle, children,_evolutionScheme,fb.type); // select branchings for children if(!fc[0].kinematics) { // select branching for first particle if(!fb.hard && fb.iout ==1 ) fc[0] = selectTimeLikeBranching(children[0],type,branch); else if(fb.hard && !branch->children()[0]->children().empty() ) fc[0] = selectTimeLikeBranching(children[0],type,branch->children()[0]); else fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr()); } // select branching for the second particle if(!fc[1].kinematics) { // select branching for first particle if(!fb.hard && fb.iout ==2 ) fc[1] = selectTimeLikeBranching(children[1],type,branch); else if(fb.hard && !branch->children()[1]->children().empty() ) fc[1] = selectTimeLikeBranching(children[1],type,branch->children()[1]); else fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr()); } // shower the first particle if(fc[0].kinematics) { // the parent has truncated emission and following line if(!fb.hard && fb.iout == 1) truncatedTimeLikeShower(children[0],branch,type,fc[0],false); // hard emission and subsquent hard emissions else if(fb.hard && !branch->children()[0]->children().empty() ) truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false); // normal shower else timeLikeShower(children[0],type,fc[0],false); } if(children[0]->spinInfo()) children[0]->spinInfo()->develop(); // shower the second particle if(fc[1].kinematics) { // the parent has truncated emission and following line if(!fb.hard && fb.iout == 2) truncatedTimeLikeShower(children[1],branch,type,fc[1],false); // hard emission and subsquent hard emissions else if(fb.hard && !branch->children()[1]->children().empty() ) truncatedTimeLikeShower(children[1],branch->children()[1],type,fc[1],false); else timeLikeShower(children[1],type,fc[1],false); } if(children[1]->spinInfo()) children[1]->spinInfo()->develop(); // branching has happened particle->showerKinematics()->updateParent(particle, children,_evolutionScheme,fb.type); if(first&&!children.empty()) particle->showerKinematics()->resetChildren(particle,children); if(particle->spinInfo()) particle->spinInfo()->develop(); // TODO NEED AN UPDATE HERE FOR RECONOPT!=0 ? return true; } bool QTildeShowerHandler::truncatedSpaceLikeShower(tShowerParticlePtr particle, PPtr beam, HardBranchingPtr branch, ShowerInteraction type) { tcPDFPtr pdf; if(firstPDF().particle() == beamParticle()) pdf = firstPDF().pdf(); if(secondPDF().particle() == beamParticle()) pdf = secondPDF().pdf(); Energy freeze = pdfFreezingScale(); Branching bb; // parameters of the force branching double z(0.); HardBranchingPtr timelike; for( unsigned int ix = 0; ix < branch->children().size(); ++ix ) { if( branch->children()[ix]->status() ==HardBranching::Outgoing) { timelike = branch->children()[ix]; } if( branch->children()[ix]->status() ==HardBranching::Incoming ) z = branch->children()[ix]->z(); } // generate truncated branching tcPDPtr part[2]; if(z>=0.&&z<=1.) { while (true) { if( !isTruncatedShowerON() || hardOnly() ) break; bb = splittingGenerator()->chooseBackwardBranching( *particle, beam, 1., beamParticle(), type , pdf,freeze); if( !bb.kinematics || bb.kinematics->scale() < branch->scale() ) { bb = Branching(); break; } // particles as in Sudakov form factor part[0] = bb.ids[0]; part[1] = bb.ids[2]; double zsplit = bb.kinematics->z(); // apply the vetos for the truncated shower // if doesn't carry most of momentum ShowerInteraction type2 = convertInteraction(bb.type); if(type2==branch->sudakov()->interactionType() && zsplit < 0.5) { particle->vetoEmission(bb.type,bb.kinematics->scale()); continue; } // others if( part[0]->id() != particle->id() || // if particle changes type bb.kinematics->pT() > progenitor()->maximumpT(type2) || // pt veto bb.kinematics->scale() < branch->scale()) { // angular ordering veto particle->vetoEmission(bb.type,bb.kinematics->scale()); continue; } // and those from the base class if(spaceLikeVetoed(bb,particle)) { particle->vetoEmission(bb.type,bb.kinematics->scale()); continue; } break; } } if( !bb.kinematics ) { //do the hard emission ShoKinPtr kinematics = new_ptr(IS_QTildeShowerKinematics1to2( branch->scale(), z, branch->phi(), branch->children()[0]->pT(), branch->sudakov() )); // assign the splitting function and shower kinematics particle->showerKinematics( kinematics ); if(kinematics->pT()>progenitor()->highestpT()) progenitor()->highestpT(kinematics->pT()); // For the time being we are considering only 1->2 branching // Now create the actual particles, make the otherChild a final state // particle, while the newParent is not ShowerParticlePtr newParent = new_ptr( ShowerParticle( branch->branchingParticle()->dataPtr(), false ) ); ShowerParticlePtr otherChild = new_ptr( ShowerParticle( timelike->branchingParticle()->dataPtr(), true, true ) ); ShowerParticleVector theChildren; theChildren.push_back( particle ); theChildren.push_back( otherChild ); particle->showerKinematics()-> updateParent( newParent, theChildren,_evolutionScheme, branch->type()); // update the history if needed currentTree()->updateInitialStateShowerProduct( progenitor(), newParent ); currentTree()->addInitialStateBranching( particle, newParent, otherChild ); // for the reconstruction of kinematics, parent/child // relationships are according to the branching process: // now continue the shower bool emitted=false; if(!hardOnly()) { if( branch->parent() ) { emitted = truncatedSpaceLikeShower( newParent, beam, branch->parent() , type); } else { emitted = spaceLikeShower( newParent, beam , type); } } if( !emitted ) { if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) { kinematics->updateLast( newParent, ZERO, ZERO ); } else { pair kt = intrinsicpT()[progenitor()]; kinematics->updateLast( newParent, kt.first*cos( kt.second ), kt.first*sin( kt.second ) ); } } particle->showerKinematics()-> updateChildren( newParent, theChildren,_evolutionScheme,bb.type); if(hardOnly()) return true; // perform the shower of the final-state particle if( timelike->children().empty() ) { timeLikeShower( otherChild , type,Branching(),true); } else { truncatedTimeLikeShower( otherChild, timelike , type,Branching(), true); } updateHistory(otherChild); // return the emitted return true; } // assign the splitting function and shower kinematics particle->showerKinematics( bb.kinematics ); if(bb.kinematics->pT()>progenitor()->highestpT()) progenitor()->highestpT(bb.kinematics->pT()); // For the time being we are considering only 1->2 branching // Now create the actual particles, make the otherChild a final state // particle, while the newParent is not ShowerParticlePtr newParent = new_ptr( ShowerParticle( part[0], false ) ); ShowerParticlePtr otherChild = new_ptr( ShowerParticle( part[1], true, true ) ); ShowerParticleVector theChildren; theChildren.push_back( particle ); theChildren.push_back( otherChild ); particle->showerKinematics()-> updateParent( newParent, theChildren,_evolutionScheme, bb.type); // update the history if needed currentTree()->updateInitialStateShowerProduct( progenitor(), newParent ); currentTree()->addInitialStateBranching( particle, newParent, otherChild ); // for the reconstruction of kinematics, parent/child // relationships are according to the branching process: // now continue the shower bool emitted = truncatedSpaceLikeShower( newParent, beam, branch,type); // now reconstruct the momentum if( !emitted ) { if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) { bb.kinematics->updateLast( newParent, ZERO, ZERO ); } else { pair kt = intrinsicpT()[ progenitor() ]; bb.kinematics->updateLast( newParent, kt.first*cos( kt.second ), kt.first*sin( kt.second ) ); } } particle->showerKinematics()-> updateChildren( newParent, theChildren,_evolutionScheme, bb.type); // perform the shower of the final-state particle timeLikeShower( otherChild , type,Branching(),true); updateHistory(otherChild); // return the emitted return true; } bool QTildeShowerHandler:: truncatedSpaceLikeDecayShower(tShowerParticlePtr particle, const ShowerParticle::EvolutionScales & maxScales, Energy minmass, HardBranchingPtr branch, ShowerInteraction type, Branching fb) { // select a branching if we don't have one if(!fb.kinematics) fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,branch); // must be an emission, the forced one it not a truncated one assert(fb.kinematics); ShowerParticleVector children; Branching fc[2]={Branching(),Branching()}; // Assign the shower kinematics to the emitting particle. particle->showerKinematics(fb.kinematics); if(fb.kinematics->pT()>progenitor()->highestpT()) progenitor()->highestpT(fb.kinematics->pT()); // create the ShowerParticle objects for the two children children = createTimeLikeChildren(particle,fb.ids); // updateChildren the children particle->showerKinematics()-> updateChildren(particle, children,_evolutionScheme, fb.type); // select branchings for children if(!fc[0].kinematics) { if(children[0]->id()==particle->id()) { // select branching for first particle if(!fb.hard) fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type,branch); else if(fb.hard && ! branch->children()[0]->children().empty() ) fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type, branch->children()[0]); else fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type, HardBranchingPtr()); } else { // select branching for first particle if(fb.hard && !branch->children()[0]->children().empty() ) fc[0] = selectTimeLikeBranching(children[0],type,branch->children()[0]); else fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr()); } } // select branching for the second particle if(!fc[1].kinematics) { if(children[1]->id()==particle->id()) { // select branching for first particle if(!fb.hard) fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type,branch); else if(fb.hard && ! branch->children()[1]->children().empty() ) fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type, branch->children()[1]); else fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type, HardBranchingPtr()); } else { if(fb.hard && !branch->children()[1]->children().empty() ) fc[1] = selectTimeLikeBranching(children[1],type,branch->children()[1]); else fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr()); } } // old default // update the history if needed currentTree()->updateInitialStateShowerProduct(progenitor(),children[0]); currentTree()->addInitialStateBranching(particle,children[0],children[1]); // shower the first particle if(fc[0].kinematics) { if(children[0]->id()==particle->id()) { if(!fb.hard) truncatedSpaceLikeDecayShower( children[0],maxScales,minmass, branch,type,fc[0]); else if(fb.hard && ! branch->children()[0]->children().empty() ) truncatedSpaceLikeDecayShower( children[0],maxScales,minmass, branch->children()[0],type,fc[0]); else spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[0]); } else { if(fb.hard && !branch->children()[0]->children().empty() ) truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false); // normal shower else timeLikeShower(children[0],type,fc[0],false); } } // shower the second particle if(fc[1].kinematics) { if(children[0]->id()==particle->id()) { if(!fb.hard) truncatedSpaceLikeDecayShower( children[0],maxScales,minmass, branch,type,fc[1]); else if(fb.hard && ! branch->children()[0]->children().empty() ) truncatedSpaceLikeDecayShower( children[0],maxScales,minmass, branch->children()[0],type,fc[1]); else spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[1]); } else { if(fb.hard && !branch->children()[0]->children().empty() ) truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[1],false); // normal shower else timeLikeShower(children[0],type,fc[1],false); } } updateHistory(children[1]); // TODO DO WE NEED A CHECK IF RECONPT !=0 return true; } void QTildeShowerHandler::connectTrees(ShowerTreePtr showerTree, HardTreePtr hardTree, bool hard ) { ShowerParticleVector particles; // find the Sudakovs for(set::iterator cit=hardTree->branchings().begin(); cit!=hardTree->branchings().end();++cit) { // Sudakovs for ISR if((**cit).parent()&&(**cit).status()==HardBranching::Incoming) { ++_nis; array br; br[0] = (**cit).parent()->branchingParticle()->id(); br[1] = (**cit). branchingParticle()->id(); br[2] = (**cit).parent()->children()[0]==*cit ? (**cit).parent()->children()[1]->branchingParticle()->id() : (**cit).parent()->children()[0]->branchingParticle()->id(); BranchingList branchings = splittingGenerator()->initialStateBranchings(); if(br[1]<0&&br[0]==br[1]) { br[0] = abs(br[0]); br[1] = abs(br[1]); } else if(br[1]<0) { br[1] = -br[1]; br[2] = -br[2]; } long index = abs(br[1]); SudakovPtr sudakov; for(BranchingList::const_iterator cjt = branchings.lower_bound(index); cjt != branchings.upper_bound(index); ++cjt ) { IdList ids = cjt->second.particles; if(ids[0]->id()==br[0]&&ids[1]->id()==br[1]&&ids[2]->id()==br[2]) { sudakov=cjt->second.sudakov; break; } } if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in " << "QTildeShowerHandler::connectTrees() for ISR" << Exception::runerror; (**cit).parent()->sudakov(sudakov); } // Sudakovs for FSR else if(!(**cit).children().empty()) { ++_nfs; array br; br[0] = (**cit) .branchingParticle()->id(); br[1] = (**cit).children()[0]->branchingParticle()->id(); br[2] = (**cit).children()[1]->branchingParticle()->id(); BranchingList branchings = splittingGenerator()->finalStateBranchings(); if(br[0]<0) { br[0] = abs(br[0]); br[1] = abs(br[1]); br[2] = abs(br[2]); } long index = br[0]; SudakovPtr sudakov; for(BranchingList::const_iterator cjt = branchings.lower_bound(index); cjt != branchings.upper_bound(index); ++cjt ) { IdList ids = cjt->second.particles; if(ids[0]->id()==br[0]&&ids[1]->id()==br[1]&&ids[2]->id()==br[2]) { sudakov=cjt->second.sudakov; break; } } if(!sudakov) { throw Exception() << "Can't find Sudakov for the hard emission in " << "QTildeShowerHandler::connectTrees()" << Exception::runerror; } (**cit).sudakov(sudakov); } } // calculate the evolution scale for(set::iterator cit=hardTree->branchings().begin(); cit!=hardTree->branchings().end();++cit) { particles.push_back((*cit)->branchingParticle()); } partnerFinder()-> setInitialEvolutionScales(particles,!hard,interaction_,true); hardTree->partnersSet(true); // inverse reconstruction if(hard) { kinematicsReconstructor()-> deconstructHardJets(hardTree,interaction_); } else kinematicsReconstructor()-> deconstructDecayJets(hardTree,interaction_); // now reset the momenta of the showering particles vector particlesToShower=showerTree->extractProgenitors(); // match them map partners; for(set::const_iterator bit=hardTree->branchings().begin(); bit!=hardTree->branchings().end();++bit) { Energy2 dmin( 1e30*GeV2 ); ShowerProgenitorPtr partner; for(vector::const_iterator pit=particlesToShower.begin(); pit!=particlesToShower.end();++pit) { if(partners.find(*pit)!=partners.end()) continue; if( (**bit).branchingParticle()->id() != (**pit).progenitor()->id() ) continue; if( (**bit).branchingParticle()->isFinalState() != (**pit).progenitor()->isFinalState() ) continue; if( (**pit).progenitor()->isFinalState() ) { Energy2 dtest = sqr( (**pit).progenitor()->momentum().x() - (**bit).showerMomentum().x() ) + sqr( (**pit).progenitor()->momentum().y() - (**bit).showerMomentum().y() ) + sqr( (**pit).progenitor()->momentum().z() - (**bit).showerMomentum().z() ) + sqr( (**pit).progenitor()->momentum().t() - (**bit).showerMomentum().t() ); // add mass difference for identical particles (e.g. Z0 Z0 production) dtest += 1e10*sqr((**pit).progenitor()->momentum().m()-(**bit).showerMomentum().m()); if( dtest < dmin ) { partner = *pit; dmin = dtest; } } else { // ensure directions are right if((**pit).progenitor()->momentum().z()/(**bit).showerMomentum().z()>ZERO) { partner = *pit; break; } } } if(!partner) throw Exception() << "Failed to match shower and hard trees in QTildeShowerHandler::hardestEmission" << Exception::eventerror; partners[partner] = *bit; } for(vector::const_iterator pit=particlesToShower.begin(); pit!=particlesToShower.end();++pit) { HardBranchingPtr partner = partners[*pit]; if((**pit).progenitor()->dataPtr()->stable()) { (**pit).progenitor()->set5Momentum(partner->showerMomentum()); (**pit).copy()->set5Momentum(partner->showerMomentum()); } else { Lorentz5Momentum oldMomentum = (**pit).progenitor()->momentum(); Lorentz5Momentum newMomentum = partner->showerMomentum(); LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass()); (**pit).progenitor()->transform(boost); (**pit).copy() ->transform(boost); boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass()); (**pit).progenitor()->transform(boost); (**pit).copy() ->transform(boost); } } // correction boosts for daughter trees for(map >::const_iterator tit = showerTree->treelinks().begin(); tit != showerTree->treelinks().end();++tit) { ShowerTreePtr decayTree = tit->first; map::const_iterator cit = decayTree->incomingLines().begin(); // reset the momentum of the decay particle Lorentz5Momentum oldMomentum = cit->first->progenitor()->momentum(); Lorentz5Momentum newMomentum = tit->second.second->momentum(); LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass()); decayTree->transform(boost,true); boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass()); decayTree->transform(boost,true); } } void QTildeShowerHandler::doShowering(bool hard,XCPtr xcomb) { // zero number of emissions _nis = _nfs = 0; // if MC@NLO H event and limited emissions // indicate both final and initial state emission if ( currentTree()->isMCatNLOHEvent() && _limitEmissions != 0 ) { _nis = _nfs = 1; } // extract particles to shower vector particlesToShower(setupShower(hard)); // check if we should shower bool colCharge = false; for(unsigned int ix=0;ixprogenitor()->dataPtr()->coloured() || particlesToShower[ix]->progenitor()->dataPtr()->charged()) { colCharge = true; break; } } if(!colCharge) { _currenttree->hasShowered(true); return; } // setup the maximum scales for the shower if (restrictPhasespace()) setupMaximumScales(particlesToShower,xcomb); // set the hard scales for the profiles setupHardScales(particlesToShower,xcomb); // specific stuff for hard processes and decays Energy minmass(ZERO), mIn(ZERO); // hard process generate the intrinsic p_T once and for all bool minBias = dynamic_ptr_cast::ptr>(_hardme); if(hard) { _intrinsic.clear(); if(!minBias) generateIntrinsicpT(particlesToShower); } // decay compute the minimum mass of the final-state else { for(unsigned int ix=0;ixprogenitor()->isFinalState()) { if(particlesToShower[ix]->progenitor()->dataPtr()->stable()) { auto dm= ShowerHandler::currentHandler()->retConstituentMasses()? particlesToShower[ix]->progenitor()->dataPtr()->constituentMass(): particlesToShower[ix]->progenitor()->dataPtr()->mass(); minmass += dm; } else minmass += particlesToShower[ix]->progenitor()->mass(); } else { mIn = particlesToShower[ix]->progenitor()->mass(); } } // throw exception if decay can't happen if ( minmass > mIn ) { throw Exception() << "QTildeShowerHandler.cc: Mass of decaying particle is " << "below constituent masses of decay products." << Exception::eventerror; } } // setup for reweighted bool reWeighting = _reWeight && hard && ShowerHandler::currentHandler()->firstInteraction(); double eventWeight=0.; unsigned int nTryReWeight(0); // create random particle vector (only need to do once) vector tmp; unsigned int nColouredIncoming = 0; while(particlesToShower.size()>0){ unsigned int xx=UseRandom::irnd(particlesToShower.size()); tmp.push_back(particlesToShower[xx]); particlesToShower.erase(particlesToShower.begin()+xx); } particlesToShower=tmp; for(unsigned int ix=0;ixprogenitor()->isFinalState() && particlesToShower[ix]->progenitor()->coloured()) ++nColouredIncoming; } bool switchRecon = hard && nColouredIncoming !=1; // main shower loop unsigned int ntry(0); bool reconstructed = false; do { // type of recon to use bool generalRecon = minBias >> (switchRecon && ntry>maximumTries()/2); // clear results of last attempt if needed if(ntry!=0) { currentTree()->clear(); setEvolutionPartners(hard,interaction_,true); _nis = _nfs = 0; // if MC@NLO H event and limited emissions // indicate both final and initial state emission if ( currentTree()->isMCatNLOHEvent() && _limitEmissions != 0 ) { _nis = _nfs = 1; } for(unsigned int ix=0; ixprogenitor()->spinInfo(); if(spin && spin->decayVertex() && dynamic_ptr_cast(spin->decayVertex())) { spin->decayVertex(VertexPtr()); } } for(unsigned int ix=0;ixprogenitor()->isFinalState() || (hard && !particlesToShower[ix]->progenitor()->isFinalState())) { if(particlesToShower[ix]->progenitor()->spinInfo()) particlesToShower[ix]->progenitor()->spinInfo()->reset(); } } } // loop over particles for(unsigned int ix=0;ixprogenitor()->isFinalState()) { if(!doFSR()) continue; // perform shower progenitor()->hasEmitted(startTimeLikeShower(interaction_)); } // initial-state radiation else { if(!doISR()) continue; // hard process if(hard) { // get the PDF setBeamParticle(_progenitor->beam()); if(!beamParticle()) { throw Exception() << "Incorrect type of beam particle in " << "QTildeShowerHandler::doShowering(). " << "This should not happen for conventional choices but may happen if you have used a" << " non-default choice and have not changed the create ParticleData line in the input files" << " for this particle to create BeamParticleData." << Exception::runerror; } // perform the shower // set the beam particle tPPtr beamparticle=progenitor()->original(); if(!beamparticle->parents().empty()) beamparticle=beamparticle->parents()[0]; // generate the shower progenitor()->hasEmitted(startSpaceLikeShower(beamparticle, interaction_)); } // decay else { // skip colour and electrically neutral particles if(!progenitor()->progenitor()->dataPtr()->coloured() && !progenitor()->progenitor()->dataPtr()->charged()) { progenitor()->hasEmitted(false); continue; } // perform shower // set the scales correctly. The current scale is the maximum scale for // emission not the starting scale ShowerParticle::EvolutionScales maxScales(progenitor()->progenitor()->scales()); progenitor()->progenitor()->scales() = ShowerParticle::EvolutionScales(); if(progenitor()->progenitor()->dataPtr()->charged()) { progenitor()->progenitor()->scales().QED = progenitor()->progenitor()->mass(); progenitor()->progenitor()->scales().QED_noAO = progenitor()->progenitor()->mass(); } if(progenitor()->progenitor()->hasColour()) { progenitor()->progenitor()->scales().QCD_c = progenitor()->progenitor()->mass(); progenitor()->progenitor()->scales().QCD_c_noAO = progenitor()->progenitor()->mass(); } if(progenitor()->progenitor()->hasAntiColour()) { progenitor()->progenitor()->scales().QCD_ac = progenitor()->progenitor()->mass(); progenitor()->progenitor()->scales().QCD_ac_noAO = progenitor()->progenitor()->mass(); } // perform the shower progenitor()->hasEmitted(startSpaceLikeDecayShower(maxScales,minmass, interaction_)); } } } // do the kinematic reconstruction, checking if it worked reconstructed = hard ? kinematicsReconstructor()-> reconstructHardJets (currentTree(),intrinsicpT(),interaction_,generalRecon) : kinematicsReconstructor()-> reconstructDecayJets(currentTree(),interaction_); if(!reconstructed) continue; // apply vetos on the full shower for(vector::const_iterator it=_fullShowerVetoes.begin(); it!=_fullShowerVetoes.end();++it) { int veto = (**it).applyVeto(currentTree()); if(veto<0) continue; // veto the shower if(veto==0) { reconstructed = false; break; } // veto the shower and reweight else if(veto==1) { reconstructed = false; break; } // veto the event else if(veto==2) { throw Veto(); } } if(reWeighting) { if(reconstructed) eventWeight += 1.; reconstructed=false; ++nTryReWeight; if(nTryReWeight==_nReWeight) { reWeighting = false; if(eventWeight==0.) throw Veto(); } } } while(!reconstructed&&maximumTries()>++ntry); // check if failed to generate the shower if(ntry==maximumTries()) { if(hard) throw ShowerHandler::ShowerTriesVeto(ntry); else throw Exception() << "Failed to generate the shower after " << ntry << " attempts in QTildeShowerHandler::showerDecay()" << Exception::eventerror; } // handle the weights and apply any reweighting required if(nTryReWeight>0) { tStdEHPtr seh = dynamic_ptr_cast(generator()->currentEventHandler()); static bool first = true; if(seh) { seh->reweight(eventWeight/double(nTryReWeight)); } else if(first) { generator()->log() << "Reweighting the shower only works with internal Herwig7 processes" << "Presumably you are showering Les Houches Events. These will not be" << "reweighted\n"; first = false; } } // tree has now showered _currenttree->hasShowered(true); hardTree(HardTreePtr()); } void QTildeShowerHandler:: convertHardTree(bool hard,ShowerInteraction type) { map cmap; // incoming particles for(map::const_iterator cit=currentTree()->incomingLines().begin();cit!=currentTree()->incomingLines().end();++cit) { map::const_iterator mit = hardTree()->particles().find(cit->first->progenitor()); // put the colour lines in the map ShowerParticlePtr oldParticle = cit->first->progenitor(); ShowerParticlePtr newParticle = mit->second->branchingParticle(); ColinePtr cLine = oldParticle-> colourLine(); ColinePtr aLine = oldParticle->antiColourLine(); if(newParticle->colourLine() && cmap.find(newParticle-> colourLine())==cmap.end()) cmap[newParticle-> colourLine()] = cLine; if(newParticle->antiColourLine() && cmap.find(newParticle->antiColourLine())==cmap.end()) cmap[newParticle->antiColourLine()] = aLine; // check whether or not particle emits bool emission = mit->second->parent(); if(emission) { if(newParticle->colourLine()) { ColinePtr ctemp = newParticle-> colourLine(); ctemp->removeColoured(newParticle); } if(newParticle->antiColourLine()) { ColinePtr ctemp = newParticle->antiColourLine(); ctemp->removeAntiColoured(newParticle); } newParticle = mit->second->parent()->branchingParticle(); } // get the new colour lines ColinePtr newCLine,newALine; // sort out colour lines if(newParticle->colourLine()) { ColinePtr ctemp = newParticle-> colourLine(); ctemp->removeColoured(newParticle); if(cmap.find(ctemp)!=cmap.end()) { newCLine = cmap[ctemp]; } else { newCLine = new_ptr(ColourLine()); cmap[ctemp] = newCLine; } } // and anticolour lines if(newParticle->antiColourLine()) { ColinePtr ctemp = newParticle->antiColourLine(); ctemp->removeAntiColoured(newParticle); if(cmap.find(ctemp)!=cmap.end()) { newALine = cmap[ctemp]; } else { newALine = new_ptr(ColourLine()); cmap[ctemp] = newALine; } } // remove colour lines from old particle if(aLine) { aLine->removeAntiColoured(cit->first->copy()); aLine->removeAntiColoured(cit->first->progenitor()); } if(cLine) { cLine->removeColoured(cit->first->copy()); cLine->removeColoured(cit->first->progenitor()); } // add particle to colour lines if(newCLine) newCLine->addColoured (newParticle); if(newALine) newALine->addAntiColoured(newParticle); // insert new particles cit->first->copy(newParticle); ShowerParticlePtr sp(new_ptr(ShowerParticle(*newParticle,1,false))); cit->first->progenitor(sp); currentTree()->incomingLines()[cit->first]=sp; cit->first->perturbative(!emission); // and the emitted particle if needed if(emission) { ShowerParticlePtr newOut = mit->second->parent()->children()[1]->branchingParticle(); if(newOut->colourLine()) { ColinePtr ctemp = newOut-> colourLine(); ctemp->removeColoured(newOut); assert(cmap.find(ctemp)!=cmap.end()); cmap[ctemp]->addColoured (newOut); } if(newOut->antiColourLine()) { ColinePtr ctemp = newOut->antiColourLine(); ctemp->removeAntiColoured(newOut); assert(cmap.find(ctemp)!=cmap.end()); cmap[ctemp]->addAntiColoured(newOut); } ShowerParticlePtr sout=new_ptr(ShowerParticle(*newOut,1,true)); ShowerProgenitorPtr out=new_ptr(ShowerProgenitor(cit->first->original(),newOut,sout)); out->perturbative(false); currentTree()->outgoingLines().insert(make_pair(out,sout)); } if(hard) { // sort out the value of x if(mit->second->beam()->momentum().z()>ZERO) { sp->x(newParticle->momentum(). plus()/mit->second->beam()->momentum(). plus()); } else { sp->x(newParticle->momentum().minus()/mit->second->beam()->momentum().minus()); } } } // outgoing particles for(map::const_iterator cit=currentTree()->outgoingLines().begin();cit!=currentTree()->outgoingLines().end();++cit) { map >::const_iterator tit; for(tit = currentTree()->treelinks().begin(); tit != currentTree()->treelinks().end();++tit) { if(tit->second.first && tit->second.second==cit->first->progenitor()) break; } map::const_iterator mit = hardTree()->particles().find(cit->first->progenitor()); if(mit==hardTree()->particles().end()) continue; // put the colour lines in the map ShowerParticlePtr oldParticle = cit->first->progenitor(); ShowerParticlePtr newParticle = mit->second->branchingParticle(); ShowerParticlePtr newOut; ColinePtr cLine = oldParticle-> colourLine(); ColinePtr aLine = oldParticle->antiColourLine(); if(newParticle->colourLine() && cmap.find(newParticle-> colourLine())==cmap.end()) cmap[newParticle-> colourLine()] = cLine; if(newParticle->antiColourLine() && cmap.find(newParticle->antiColourLine())==cmap.end()) cmap[newParticle->antiColourLine()] = aLine; // check whether or not particle emits bool emission = !mit->second->children().empty(); if(emission) { if(newParticle->colourLine()) { ColinePtr ctemp = newParticle-> colourLine(); ctemp->removeColoured(newParticle); } if(newParticle->antiColourLine()) { ColinePtr ctemp = newParticle->antiColourLine(); ctemp->removeAntiColoured(newParticle); } newParticle = mit->second->children()[0]->branchingParticle(); newOut = mit->second->children()[1]->branchingParticle(); if(newParticle->id()!=oldParticle->id()&&newParticle->id()==newOut->id()) swap(newParticle,newOut); } // get the new colour lines ColinePtr newCLine,newALine; // sort out colour lines if(newParticle->colourLine()) { ColinePtr ctemp = newParticle-> colourLine(); ctemp->removeColoured(newParticle); if(cmap.find(ctemp)!=cmap.end()) { newCLine = cmap[ctemp]; } else { newCLine = new_ptr(ColourLine()); cmap[ctemp] = newCLine; } } // and anticolour lines if(newParticle->antiColourLine()) { ColinePtr ctemp = newParticle->antiColourLine(); ctemp->removeAntiColoured(newParticle); if(cmap.find(ctemp)!=cmap.end()) { newALine = cmap[ctemp]; } else { newALine = new_ptr(ColourLine()); cmap[ctemp] = newALine; } } // remove colour lines from old particle if(aLine) { aLine->removeAntiColoured(cit->first->copy()); aLine->removeAntiColoured(cit->first->progenitor()); } if(cLine) { cLine->removeColoured(cit->first->copy()); cLine->removeColoured(cit->first->progenitor()); } // special for unstable particles if(newParticle->id()==oldParticle->id() && (tit!=currentTree()->treelinks().end()||!oldParticle->dataPtr()->stable())) { Lorentz5Momentum oldMomentum = oldParticle->momentum(); Lorentz5Momentum newMomentum = newParticle->momentum(); LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass()); if(tit!=currentTree()->treelinks().end()) tit->first->transform(boost,false); oldParticle->transform(boost); boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass()); oldParticle->transform(boost); if(tit!=currentTree()->treelinks().end()) tit->first->transform(boost,false); newParticle=oldParticle; } // add particle to colour lines if(newCLine) newCLine->addColoured (newParticle); if(newALine) newALine->addAntiColoured(newParticle); // insert new particles cit->first->copy(newParticle); ShowerParticlePtr sp(new_ptr(ShowerParticle(*newParticle,1,true))); cit->first->progenitor(sp); currentTree()->outgoingLines()[cit->first]=sp; cit->first->perturbative(!emission); // and the emitted particle if needed if(emission) { if(newOut->colourLine()) { ColinePtr ctemp = newOut-> colourLine(); ctemp->removeColoured(newOut); assert(cmap.find(ctemp)!=cmap.end()); cmap[ctemp]->addColoured (newOut); } if(newOut->antiColourLine()) { ColinePtr ctemp = newOut->antiColourLine(); ctemp->removeAntiColoured(newOut); assert(cmap.find(ctemp)!=cmap.end()); cmap[ctemp]->addAntiColoured(newOut); } ShowerParticlePtr sout=new_ptr(ShowerParticle(*newOut,1,true)); ShowerProgenitorPtr out=new_ptr(ShowerProgenitor(cit->first->original(),newOut,sout)); out->perturbative(false); currentTree()->outgoingLines().insert(make_pair(out,sout)); } // update any decay products if(tit!=currentTree()->treelinks().end()) currentTree()->updateLink(tit->first,make_pair(cit->first,sp)); } // reset the tree currentTree()->resetShowerProducts(); // reextract the particles and set the colour partners vector particles = currentTree()->extractProgenitorParticles(); // clear the partners for(unsigned int ix=0;ixpartner(ShowerParticlePtr()); particles[ix]->clearPartners(); } // clear the tree hardTree(HardTreePtr()); // Set the initial evolution scales partnerFinder()-> setInitialEvolutionScales(particles,!hard,type,!_hardtree); } Branching QTildeShowerHandler::selectTimeLikeBranching(tShowerParticlePtr particle, ShowerInteraction type, HardBranchingPtr branch) { Branching fb; unsigned int iout=0; while (true) { // break if doing truncated shower and no truncated shower needed if(branch && (!isTruncatedShowerON()||hardOnly())) break; fb=_splittingGenerator->chooseForwardBranching(*particle,_finalenhance,type); // no emission break if(!fb.kinematics) break; // special for truncated shower if(branch) { // check haven't evolved too far if(fb.kinematics->scale() < branch->scale()) { fb=Branching(); break; } // find the truncated line iout=0; if(fb.ids[1]->id()!=fb.ids[2]->id()) { if(fb.ids[1]->id()==particle->id()) iout=1; else if (fb.ids[2]->id()==particle->id()) iout=2; } else if(fb.ids[1]->id()==particle->id()) { if(fb.kinematics->z()>0.5) iout=1; else iout=2; } // apply the vetos for the truncated shower // no flavour changing branchings if(iout==0) { particle->vetoEmission(fb.type,fb.kinematics->scale()); continue; } double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z(); // only if same interaction for forced branching ShowerInteraction type2 = convertInteraction(fb.type); // and evolution if(type2==branch->sudakov()->interactionType()) { if(zsplit < 0.5 || // hardest line veto fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto particle->vetoEmission(fb.type,fb.kinematics->scale()); continue; } } // pt veto if(fb.kinematics->pT() > progenitor()->maximumpT(type2)) { particle->vetoEmission(fb.type,fb.kinematics->scale()); continue; } } // standard vetos for all emissions if(timeLikeVetoed(fb,particle)) { particle->vetoEmission(fb.type,fb.kinematics->scale()); if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr()); continue; } // special for already decayed particles // don't allow flavour changing branchings bool vetoDecay = false; for(map >::const_iterator tit = currentTree()->treelinks().begin(); tit != currentTree()->treelinks().end();++tit) { if(tit->second.first == progenitor()) { map::const_iterator it = currentTree()->outgoingLines().find(progenitor()); if(it!=currentTree()->outgoingLines().end() && particle == it->second && fb.ids[0]!=fb.ids[1] && fb.ids[1]!=fb.ids[2]) { vetoDecay = true; break; } } } if(vetoDecay) { particle->vetoEmission(fb.type,fb.kinematics->scale()); if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr()); continue; } break; } // normal case if(!branch) { if(fb.kinematics) fb.hard = false; return fb; } // truncated emission if(fb.kinematics) { fb.hard = false; fb.iout = iout; return fb; } // otherwise need to return the hard emission // construct the kinematics for the hard emission ShoKinPtr showerKin = new_ptr(FS_QTildeShowerKinematics1to2( branch->scale(), branch->children()[0]->z(), branch->phi(), branch->children()[0]->pT(), branch->sudakov() )); IdList idlist(3); idlist[0] = particle->dataPtr(); idlist[1] = branch->children()[0]->branchingParticle()->dataPtr(); idlist[2] = branch->children()[1]->branchingParticle()->dataPtr(); fb = Branching( showerKin, idlist, branch->sudakov(),branch->type() ); fb.hard = true; fb.iout=0; // return it return fb; } Branching QTildeShowerHandler::selectSpaceLikeDecayBranching(tShowerParticlePtr particle, const ShowerParticle::EvolutionScales & maxScales, Energy minmass,ShowerInteraction type, HardBranchingPtr branch) { Branching fb; unsigned int iout=0; while (true) { // break if doing truncated shower and no truncated shower needed if(branch && (!isTruncatedShowerON()||hardOnly())) break; // select branching fb=_splittingGenerator->chooseDecayBranching(*particle,maxScales,minmass, _initialenhance,type); // return if no radiation if(!fb.kinematics) break; // special for truncated shower if(branch) { // check haven't evolved too far if(fb.kinematics->scale() < branch->scale()) { fb=Branching(); break; } // find the truncated line iout=0; if(fb.ids[1]->id()!=fb.ids[2]->id()) { if(fb.ids[1]->id()==particle->id()) iout=1; else if (fb.ids[2]->id()==particle->id()) iout=2; } else if(fb.ids[1]->id()==particle->id()) { if(fb.kinematics->z()>0.5) iout=1; else iout=2; } // apply the vetos for the truncated shower // no flavour changing branchings if(iout==0) { particle->vetoEmission(fb.type,fb.kinematics->scale()); continue; } ShowerInteraction type2 = convertInteraction(fb.type); double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z(); if(type2==branch->sudakov()->interactionType()) { if(zsplit < 0.5 || // hardest line veto fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto particle->vetoEmission(fb.type,fb.kinematics->scale()); continue; } } // pt veto if(fb.kinematics->pT() > progenitor()->maximumpT(type2)) { particle->vetoEmission(fb.type,fb.kinematics->scale()); continue; } } // if not vetoed break if(spaceLikeDecayVetoed(fb,particle)) { // otherwise reset scale and continue particle->vetoEmission(fb.type,fb.kinematics->scale()); continue; } break; } // normal case if(!branch) { if(fb.kinematics) fb.hard = false; return fb; } // truncated emission if(fb.kinematics) { fb.hard = false; fb.iout = iout; return fb; } // otherwise need to return the hard emission // construct the kinematics for the hard emission ShoKinPtr showerKin = new_ptr(Decay_QTildeShowerKinematics1to2( branch->scale(), branch->children()[0]->z(), branch->phi(), branch->children()[0]->pT(), branch->sudakov())); IdList idlist(3); idlist[0] = particle->dataPtr(); idlist[1] = branch->children()[0]->branchingParticle()->dataPtr(); idlist[2] = branch->children()[1]->branchingParticle()->dataPtr(); // create the branching fb = Branching( showerKin, idlist, branch->sudakov(),ShowerPartnerType::QCDColourLine ); fb.hard=true; fb.iout=0; // return it return fb; } void QTildeShowerHandler::checkFlags() { string error = "Inconsistent hard emission set-up in QTildeShowerHandler::showerHardProcess(). "; if ( ( currentTree()->isMCatNLOSEvent() || currentTree()->isMCatNLOHEvent() ) ) { if (_hardEmission ==2 ) throw Exception() << error << "Cannot generate POWHEG matching with MC@NLO shower " << "approximation. Add 'set QTildeShowerHandler:HardEmission 0' to input file." << Exception::runerror; if ( canHandleMatchboxTrunc() ) throw Exception() << error << "Cannot use truncated qtilde shower with MC@NLO shower " << "approximation. Set LHCGenerator:EventHandler" << ":CascadeHandler to '/Herwig/Shower/ShowerHandler' or " << "'/Herwig/Shower/Dipole/DipoleShowerHandler'." << Exception::runerror; } else if ( ((currentTree()->isPowhegSEvent() || currentTree()->isPowhegHEvent()) ) && _hardEmission != 2){ if ( canHandleMatchboxTrunc()) throw Exception() << error << "Unmatched events requested for POWHEG shower " << "approximation. Set QTildeShowerHandler:HardEmission to " << "'POWHEG'." << Exception::runerror; else if (_hardEmissionWarn) { _hardEmissionWarn = false; _hardEmission=2; throw Exception() << error << "Unmatched events requested for POWHEG shower " << "approximation. Changing QTildeShowerHandler:HardEmission from " << _hardEmission << " to 2" << Exception::warning; } } if ( currentTree()->isPowhegSEvent() || currentTree()->isPowhegHEvent()) { if (currentTree()->showerApproximation()->needsTruncatedShower() && !canHandleMatchboxTrunc() ) throw Exception() << error << "Current shower handler cannot generate truncated shower. " << "Set Generator:EventHandler:CascadeHandler to " << "'/Herwig/Shower/PowhegShowerHandler'." << Exception::runerror; } else if ( currentTree()->truncatedShower() && _missingTruncWarn) { _missingTruncWarn=false; throw Exception() << "Warning: POWHEG shower approximation used without " << "truncated shower. Set Generator:EventHandler:" << "CascadeHandler to '/Herwig/Shower/PowhegShowerHandler' and " << "'MEMatching:TruncatedShower Yes'." << Exception::warning; } // else if ( !dipme && _hardEmissionMode > 1 && // firstInteraction()) // throw Exception() << error // << "POWHEG matching requested for LO events. Include " // << "'set Factory:ShowerApproximation MEMatching' in input file." // << Exception::runerror; } tPPair QTildeShowerHandler::remakeRemnant(tPPair oldp){ // get the parton extractor PartonExtractor & pex = *lastExtractor(); // get the new partons tPPair newp = make_pair(findFirstParton(oldp.first ), findFirstParton(oldp.second)); // if the same do nothing if(newp == oldp) return oldp; // Creates the new remnants and returns the new PartonBinInstances // ATTENTION Broken here for very strange configuration PBIPair newbins = pex.newRemnants(oldp, newp, newStep()); newStep()->addIntermediate(newp.first); newStep()->addIntermediate(newp.second); // return the new partons return newp; } PPtr QTildeShowerHandler::findFirstParton(tPPtr seed) const{ if(seed->parents().empty()) return seed; tPPtr parent = seed->parents()[0]; //if no parent there this is a loose end which will //be connected to the remnant soon. if(!parent || parent == incomingBeams().first || parent == incomingBeams().second ) return seed; else return findFirstParton(parent); } void QTildeShowerHandler::decay(ShowerTreePtr tree, ShowerDecayMap & decay) { // must be one incoming particle assert(tree->incomingLines().size()==1); // apply any transforms tree->applyTransforms(); // if already decayed return if(!tree->outgoingLines().empty()) return; // now we need to replace the particle with a new copy after the shower // find particle after the shower map >::const_iterator tit = tree->parent()->treelinks().find(tree); assert(tit!=tree->parent()->treelinks().end()); ShowerParticlePtr newparent=tit->second.second; PerturbativeProcessPtr newProcess = new_ptr(PerturbativeProcess()); newProcess->incoming().push_back(make_pair(newparent,PerturbativeProcessPtr())); DecayProcessMap decayMap; ShowerHandler::decay(newProcess,decayMap); ShowerTree::constructTrees(tree,decay,newProcess,decayMap); } namespace { ShowerProgenitorPtr findFinalStateLine(ShowerTreePtr tree, long id, Lorentz5Momentum momentum) { map::iterator partner; Energy2 dmin(1e30*GeV2); for(map::iterator cit =tree->outgoingLines().begin(); cit!=tree->outgoingLines().end(); ++cit) { if(cit->second->id()!=id) continue; Energy2 test = sqr(cit->second->momentum().x()-momentum.x())+ sqr(cit->second->momentum().y()-momentum.y())+ sqr(cit->second->momentum().z()-momentum.z())+ sqr(cit->second->momentum().t()-momentum.t()); if(testfirst; } ShowerProgenitorPtr findInitialStateLine(ShowerTreePtr tree, long id, Lorentz5Momentum momentum) { map::iterator partner; Energy2 dmin(1e30*GeV2); for(map::iterator cit =tree->incomingLines().begin(); cit!=tree->incomingLines().end(); ++cit) { if(cit->second->id()!=id) continue; Energy2 test = sqr(cit->second->momentum().x()-momentum.x())+ sqr(cit->second->momentum().y()-momentum.y())+ sqr(cit->second->momentum().z()-momentum.z())+ sqr(cit->second->momentum().t()-momentum.t()); if(testfirst; } void fixSpectatorColours(PPtr newSpect,ShowerProgenitorPtr oldSpect, ColinePair & cline,ColinePair & aline, bool reconnect) { cline.first = oldSpect->progenitor()->colourLine(); cline.second = newSpect->colourLine(); aline.first = oldSpect->progenitor()->antiColourLine(); aline.second = newSpect->antiColourLine(); if(!reconnect) return; if(cline.first) { cline.first ->removeColoured(oldSpect->copy()); cline.first ->removeColoured(oldSpect->progenitor()); cline.second->removeColoured(newSpect); cline.first ->addColoured(newSpect); } if(aline.first) { aline.first ->removeAntiColoured(oldSpect->copy()); aline.first ->removeAntiColoured(oldSpect->progenitor()); aline.second->removeAntiColoured(newSpect); aline.first ->addAntiColoured(newSpect); } } void fixInitialStateEmitter(ShowerTreePtr tree, PPtr newEmit,PPtr emitted, ShowerProgenitorPtr emitter, ColinePair cline,ColinePair aline,double x) { // sort out the colours if(emitted->dataPtr()->iColour()==PDT::Colour8) { // emitter if(cline.first && cline.first == emitter->progenitor()->antiColourLine() && cline.second !=newEmit->antiColourLine()) { // sort out not radiating line ColinePtr col = emitter->progenitor()->colourLine(); if(col) { col->removeColoured(emitter->copy()); col->removeColoured(emitter->progenitor()); newEmit->colourLine()->removeColoured(newEmit); col->addColoured(newEmit); } } else if(aline.first && aline.first == emitter->progenitor()->colourLine() && aline.second !=newEmit->colourLine()) { // sort out not radiating line ColinePtr anti = emitter->progenitor()->antiColourLine(); if(anti) { anti->removeAntiColoured(emitter->copy()); anti->removeAntiColoured(emitter->progenitor()); newEmit->colourLine()->removeAntiColoured(newEmit); anti->addAntiColoured(newEmit); } } else assert(false); // emitted if(cline.first && cline.second==emitted->colourLine()) { cline.second->removeColoured(emitted); cline.first->addColoured(emitted); } else if(aline.first && aline.second==emitted->antiColourLine()) { aline.second->removeAntiColoured(emitted); aline.first->addAntiColoured(emitted); } else assert(false); } else { if(emitter->progenitor()->antiColourLine() ) { ColinePtr col = emitter->progenitor()->antiColourLine(); col->removeAntiColoured(emitter->copy()); col->removeAntiColoured(emitter->progenitor()); if(newEmit->antiColourLine()) { newEmit->antiColourLine()->removeAntiColoured(newEmit); col->addAntiColoured(newEmit); } else if (emitted->colourLine()) { emitted->colourLine()->removeColoured(emitted); col->addColoured(emitted); } else assert(false); } if(emitter->progenitor()->colourLine() ) { ColinePtr col = emitter->progenitor()->colourLine(); col->removeColoured(emitter->copy()); col->removeColoured(emitter->progenitor()); if(newEmit->colourLine()) { newEmit->colourLine()->removeColoured(newEmit); col->addColoured(newEmit); } else if (emitted->antiColourLine()) { emitted->antiColourLine()->removeAntiColoured(emitted); col->addAntiColoured(emitted); } else assert(false); } } // update the emitter emitter->copy(newEmit); ShowerParticlePtr sp = new_ptr(ShowerParticle(*newEmit,1,false)); sp->x(x); emitter->progenitor(sp); tree->incomingLines()[emitter]=sp; emitter->perturbative(false); // add emitted sp=new_ptr(ShowerParticle(*emitted,1,true)); ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(emitter->original(),emitted,sp)); gluon->perturbative(false); tree->outgoingLines().insert(make_pair(gluon,sp)); } void fixFinalStateEmitter(ShowerTreePtr tree, PPtr newEmit,PPtr emitted, ShowerProgenitorPtr emitter, ColinePair cline,ColinePair aline) { map >::const_iterator tit; // special case if decayed for(tit = tree->treelinks().begin(); tit != tree->treelinks().end();++tit) { if(tit->second.first && tit->second.second==emitter->progenitor()) break; } // sort out the colour lines if(cline.first && cline.first == emitter->progenitor()->antiColourLine() && cline.second !=newEmit->antiColourLine()) { // sort out not radiating line ColinePtr col = emitter->progenitor()->colourLine(); if(col) { col->removeColoured(emitter->copy()); col->removeColoured(emitter->progenitor()); newEmit->colourLine()->removeColoured(newEmit); col->addColoured(newEmit); } } else if(aline.first && aline.first == emitter->progenitor()->colourLine() && aline.second !=newEmit->colourLine()) { // sort out not radiating line ColinePtr anti = emitter->progenitor()->antiColourLine(); if(anti) { anti->removeAntiColoured(emitter->copy()); anti->removeAntiColoured(emitter->progenitor()); newEmit->colourLine()->removeAntiColoured(newEmit); anti->addAntiColoured(newEmit); } } else assert(false); // update the emitter emitter->copy(newEmit); ShowerParticlePtr sp = new_ptr(ShowerParticle(*newEmit,1,true)); emitter->progenitor(sp); tree->outgoingLines()[emitter]=sp; emitter->perturbative(false); // update for decaying particles if(tit!=tree->treelinks().end()) tree->updateLink(tit->first,make_pair(emitter,sp)); // add the emitted particle // sort out the colour if(cline.first && cline.second==emitted->antiColourLine()) { cline.second->removeAntiColoured(emitted); cline.first->addAntiColoured(emitted); } else if(aline.first && aline.second==emitted->colourLine()) { aline.second->removeColoured(emitted); aline.first->addColoured(emitted); } else assert(false); sp=new_ptr(ShowerParticle(*emitted,1,true)); ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(emitter->original(), emitted,sp)); gluon->perturbative(false); tree->outgoingLines().insert(make_pair(gluon,sp)); } } void QTildeShowerHandler::setupMECorrection(RealEmissionProcessPtr real) { assert(real); // II emission if(real->emitter() < real->incoming().size() && real->spectator() < real->incoming().size()) { // recoiling system for( map::const_iterator cjt= currentTree()->outgoingLines().begin(); cjt != currentTree()->outgoingLines().end();++cjt ) { cjt->first->progenitor()->transform(real->transformation()); cjt->first->copy()->transform(real->transformation()); } // the the radiating system ShowerProgenitorPtr emitter,spectator; unsigned int iemit = real->emitter(); unsigned int ispect = real->spectator(); int ig = int(real->emitted())-int(real->incoming().size()); emitter = findInitialStateLine(currentTree(), real->bornIncoming()[iemit]->id(), real->bornIncoming()[iemit]->momentum()); spectator = findInitialStateLine(currentTree(), real->bornIncoming()[ispect]->id(), real->bornIncoming()[ispect]->momentum()); // sort out the colours ColinePair cline,aline; fixSpectatorColours(real->incoming()[ispect],spectator,cline,aline,true); // update the spectator spectator->copy(real->incoming()[ispect]); ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->incoming()[ispect],1,false))); sp->x(ispect ==0 ? real->x().first :real->x().second); spectator->progenitor(sp); currentTree()->incomingLines()[spectator]=sp; spectator->perturbative(true); // now for the emitter fixInitialStateEmitter(currentTree(),real->incoming()[iemit],real->outgoing()[ig], emitter,cline,aline,iemit ==0 ? real->x().first :real->x().second); } // FF emission else if(real->emitter() >= real->incoming().size() && real->spectator() >= real->incoming().size()) { assert(real->outgoing()[real->emitted()-real->incoming().size()]->id()==ParticleID::g); // find the emitter and spectator in the shower tree ShowerProgenitorPtr emitter,spectator; int iemit = int(real->emitter())-int(real->incoming().size()); emitter = findFinalStateLine(currentTree(), real->bornOutgoing()[iemit]->id(), real->bornOutgoing()[iemit]->momentum()); int ispect = int(real->spectator())-int(real->incoming().size()); spectator = findFinalStateLine(currentTree(), real->bornOutgoing()[ispect]->id(), real->bornOutgoing()[ispect]->momentum()); map >::const_iterator tit; // first the spectator // special case if decayed for(tit = currentTree()->treelinks().begin(); tit != currentTree()->treelinks().end();++tit) { if(tit->second.first && tit->second.second==spectator->progenitor()) break; } // sort out the colours ColinePair cline,aline; fixSpectatorColours(real->outgoing()[ispect],spectator,cline,aline,true); // update the spectator spectator->copy(real->outgoing()[ispect]); ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->outgoing()[ispect],1,true))); spectator->progenitor(sp); currentTree()->outgoingLines()[spectator]=sp; spectator->perturbative(true); // update for decaying particles if(tit!=currentTree()->treelinks().end()) currentTree()->updateLink(tit->first,make_pair(spectator,sp)); // now the emitting particle int ig = int(real->emitted())-int(real->incoming().size()); fixFinalStateEmitter(currentTree(),real->outgoing()[iemit], real->outgoing()[ig], emitter,cline,aline); } // IF emission else { // scattering process if(real->incoming().size()==2) { ShowerProgenitorPtr emitter,spectator; unsigned int iemit = real->emitter(); unsigned int ispect = real->spectator(); int ig = int(real->emitted())-int(real->incoming().size()); ColinePair cline,aline; // incoming spectator if(ispect<2) { spectator = findInitialStateLine(currentTree(), real->bornIncoming()[ispect]->id(), real->bornIncoming()[ispect]->momentum()); fixSpectatorColours(real->incoming()[ispect],spectator,cline,aline,true); // update the spectator spectator->copy(real->incoming()[ispect]); ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->incoming()[ispect],1,false))); sp->x(ispect ==0 ? real->x().first :real->x().second); spectator->progenitor(sp); currentTree()->incomingLines()[spectator]=sp; spectator->perturbative(true); } // outgoing spectator else { spectator = findFinalStateLine(currentTree(), real->bornOutgoing()[ispect-real->incoming().size()]->id(), real->bornOutgoing()[ispect-real->incoming().size()]->momentum()); // special case if decayed map >::const_iterator tit; for(tit = currentTree()->treelinks().begin(); tit != currentTree()->treelinks().end();++tit) { if(tit->second.first && tit->second.second==spectator->progenitor()) break; } fixSpectatorColours(real->outgoing()[ispect-real->incoming().size()],spectator,cline,aline,true); // update the spectator spectator->copy(real->outgoing()[ispect-real->incoming().size()]); ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->outgoing()[ispect-real->incoming().size()],1,true))); spectator->progenitor(sp); currentTree()->outgoingLines()[spectator]=sp; spectator->perturbative(true); // update for decaying particles if(tit!=currentTree()->treelinks().end()) currentTree()->updateLink(tit->first,make_pair(spectator,sp)); } // incoming emitter if(iemit<2) { emitter = findInitialStateLine(currentTree(), real->bornIncoming()[iemit]->id(), real->bornIncoming()[iemit]->momentum()); fixInitialStateEmitter(currentTree(),real->incoming()[iemit],real->outgoing()[ig], emitter,aline,cline,iemit ==0 ? real->x().first :real->x().second); } // outgoing emitter else { emitter = findFinalStateLine(currentTree(), real->bornOutgoing()[iemit-real->incoming().size()]->id(), real->bornOutgoing()[iemit-real->incoming().size()]->momentum()); fixFinalStateEmitter(currentTree(),real->outgoing()[iemit-real->incoming().size()], real->outgoing()[ig],emitter,aline,cline); } } // decay process else { assert(real->spectator()==0); unsigned int iemit = real->emitter()-real->incoming().size(); int ig = int(real->emitted())-int(real->incoming().size()); ColinePair cline,aline; // incoming spectator ShowerProgenitorPtr spectator = findInitialStateLine(currentTree(), real->bornIncoming()[0]->id(), real->bornIncoming()[0]->momentum()); fixSpectatorColours(real->incoming()[0],spectator,cline,aline,false); // find the emitter ShowerProgenitorPtr emitter = findFinalStateLine(currentTree(), real->bornOutgoing()[iemit]->id(), real->bornOutgoing()[iemit]->momentum()); // recoiling system for( map::const_iterator cjt= currentTree()->outgoingLines().begin(); cjt != currentTree()->outgoingLines().end();++cjt ) { if(cjt->first==emitter) continue; cjt->first->progenitor()->transform(real->transformation()); cjt->first->copy()->transform(real->transformation()); } // sort out the emitter fixFinalStateEmitter(currentTree(),real->outgoing()[iemit], real->outgoing()[ig],emitter,aline,cline); } } // clean up the shower tree _currenttree->resetShowerProducts(); } diff --git a/Shower/QTilde/QTildeShowerHandler.h b/Shower/QTilde/QTildeShowerHandler.h --- a/Shower/QTilde/QTildeShowerHandler.h +++ b/Shower/QTilde/QTildeShowerHandler.h @@ -1,841 +1,833 @@ // -*- C++ -*- #ifndef Herwig_QTildeShowerHandler_H #define Herwig_QTildeShowerHandler_H // // This is the declaration of the QTildeShowerHandler class. // #include "QTildeShowerHandler.fh" #include "Herwig/Shower/ShowerHandler.h" #include "Herwig/Shower/QTilde/SplittingFunctions/SplittingGenerator.h" #include "Herwig/Shower/QTilde/Base/ShowerTree.h" #include "Herwig/Shower/QTilde/Base/ShowerProgenitor.fh" #include "Herwig/Shower/QTilde/Base/HardTree.h" #include "Herwig/Shower/QTilde/Base/Branching.h" #include "Herwig/Shower/QTilde/Base/ShowerVeto.h" #include "Herwig/Shower/QTilde/Base/FullShowerVeto.h" #include "Herwig/Shower/QTilde/Kinematics/KinematicsReconstructor.fh" #include "Herwig/Shower/QTilde/Base/PartnerFinder.fh" #include "Herwig/Shower/QTilde/SplittingFunctions/SudakovFormFactor.fh" #include "Herwig/MatrixElement/HwMEBase.h" #include "Herwig/Decay/HwDecayerBase.h" #include "Herwig/MatrixElement/Matchbox/Matching/ShowerApproximation.h" #include "Herwig/Shower/RealEmissionProcess.h" #include "Herwig/Utilities/Statistic.h" namespace Herwig { using namespace ThePEG; /** * The QTildeShowerHandler class. * * @see \ref QTildeShowerHandlerInterfaces "The interfaces" * defined for QTildeShowerHandler. */ class QTildeShowerHandler: public ShowerHandler { public: /** * Pointer to an XComb object */ typedef Ptr::pointer XCPtr; public: - /** @name Standard constructors and destructors. */ - //@{ /** * The default constructor. */ QTildeShowerHandler(); - /** - * The destructor. - */ - virtual ~QTildeShowerHandler(); - //@} - public: /** * At the end of the Showering, transform ShowerParticle objects * into ThePEG particles and fill the event record with them. * Notice that the parent/child relationships and the * transformation from ShowerColourLine objects into ThePEG * ColourLine ones must be properly handled. */ void fillEventRecord(); /** * Return the relevant hard scale to be used in the profile scales */ virtual Energy hardScale() const { return muPt; } /** * Hook to allow vetoing of event after showering hard sub-process * as in e.g. MLM merging. */ virtual bool showerHardProcessVeto() const { return false; } /** * Generate hard emissions for CKKW etc */ virtual HardTreePtr generateCKKW(ShowerTreePtr tree) const; /** * Members to perform the shower */ //@{ /** * Perform the shower of the hard process */ virtual void showerHardProcess(ShowerTreePtr,XCPtr); /** * Perform the shower of a decay */ virtual void showerDecay(ShowerTreePtr); //@} /** * Access to the flags and shower variables */ //@{ /** * Get the SplittingGenerator */ tSplittingGeneratorPtr splittingGenerator() const { return _splittingGenerator; } /** * Mode for hard emissions */ int hardEmission() const {return _hardEmission;} //@} /** * Connect the Hard and Shower trees */ virtual void connectTrees(ShowerTreePtr showerTree, HardTreePtr hardTree, bool hard ); /** * Access to switches for spin correlations */ //@{ /** * Soft correlations */ unsigned int softCorrelations() const { return _softOpt; } /** * Any correlations */ virtual bool correlations() const { return spinCorrelations()!=0||_softOpt!=0; } //@} public: /** * Access methods to access the objects */ //@{ /** * Access to the KinematicsReconstructor object */ tKinematicsReconstructorPtr kinematicsReconstructor() const { return _reconstructor; } /** * Access to the PartnerFinder object */ tPartnerFinderPtr partnerFinder() const { return _partnerfinder; } //@} protected: /** * Perform the shower */ void doShowering(bool hard,XCPtr); /** * Generate the hard matrix element correction */ virtual RealEmissionProcessPtr hardMatrixElementCorrection(bool); /** * Generate the hardest emission */ virtual void hardestEmission(bool hard); /** * Set up for applying a matrix element correction */ void setupMECorrection(RealEmissionProcessPtr real); /** * Extract the particles to be showered, set the evolution scales * and apply the hard matrix element correction * @param hard Whether this is a hard process or decay * @return The particles to be showered */ virtual vector setupShower(bool hard); /** * set the colour partners */ virtual void setEvolutionPartners(bool hard,ShowerInteraction, bool clear); /** * Methods to perform the evolution of an individual particle, including * recursive calling on the products */ //@{ /** * It does the forward evolution of the time-like input particle * (and recursively for all its radiation products). * accepting only emissions which conforms to the showerVariables * and soft matrix element correction. * If at least one emission has occurred then the method returns true. * @param particle The particle to be showered */ virtual bool timeLikeShower(tShowerParticlePtr particle, ShowerInteraction, Branching fb, bool first); /** * It does the backward evolution of the space-like input particle * (and recursively for all its time-like radiation products). * accepting only emissions which conforms to the showerVariables. * If at least one emission has occurred then the method returns true * @param particle The particle to be showered * @param beam The beam particle */ virtual bool spaceLikeShower(tShowerParticlePtr particle,PPtr beam, ShowerInteraction); /** * If does the forward evolution of the input on-shell particle * involved in a decay * (and recursively for all its time-like radiation products). * accepting only emissions which conforms to the showerVariables. * @param particle The particle to be showered * @param maxscale The maximum scale for the shower. * @param minimumMass The minimum mass of the final-state system */ virtual bool spaceLikeDecayShower(tShowerParticlePtr particle, const ShowerParticle::EvolutionScales & maxScales, Energy minimumMass,ShowerInteraction, Branching fb); /** * Truncated shower from a time-like particle */ virtual bool truncatedTimeLikeShower(tShowerParticlePtr particle, HardBranchingPtr branch, ShowerInteraction type, Branching fb, bool first); /** * Truncated shower from a space-like particle */ virtual bool truncatedSpaceLikeShower(tShowerParticlePtr particle,PPtr beam, HardBranchingPtr branch, ShowerInteraction type); /** * Truncated shower from a time-like particle */ virtual bool truncatedSpaceLikeDecayShower(tShowerParticlePtr particle, const ShowerParticle::EvolutionScales & maxScales, Energy minimumMass, HardBranchingPtr branch, ShowerInteraction type, Branching fb); //@} /** * Switches for matrix element corrections */ //@{ /** * Any ME correction? */ bool MECOn() const { return _hardEmission == 1; } /** * Any hard ME correction? */ bool hardMEC() const { return _hardEmission == 1 && (_meCorrMode == 1 || _meCorrMode == 2); } /** * Any soft ME correction? */ bool softMEC() const { return _hardEmission == 1 && (_meCorrMode == 1 || _meCorrMode > 2); } //@} /** * Is the truncated shower on? */ bool isTruncatedShowerON() const {return _trunc_Mode;} /** * Switch for intrinsic pT */ //@{ /** * Any intrinsic pT? */ bool ipTon() const { return _iptrms != ZERO || ( _beta == 1.0 && _gamma != ZERO && _iptmax !=ZERO ); } //@} /**@name Additional shower vetoes */ //@{ /** * Insert a veto. */ void addVeto (ShowerVetoPtr v) { _vetoes.push_back(v); } /** * Remove a veto. */ void removeVeto (ShowerVetoPtr v) { vector::iterator vit = find(_vetoes.begin(),_vetoes.end(),v); if (vit != _vetoes.end()) _vetoes.erase(vit); } //@} /** * Switches for vetoing hard emissions */ //@{ /** * Returns true if the hard veto read-in is to be applied to only * the primary collision and false otherwise. */ bool hardVetoReadOption() const {return _hardVetoReadOption;} //@} /** * Enhancement factors for radiation needed to generate the soft matrix * element correction. */ //@{ /** * Access the enhancement factor for initial-state radiation */ double initialStateRadiationEnhancementFactor() const { return _initialenhance; } /** * Access the enhancement factor for final-state radiation */ double finalStateRadiationEnhancementFactor() const { return _finalenhance; } /** * Set the enhancement factor for initial-state radiation */ void initialStateRadiationEnhancementFactor(double in) { _initialenhance=in; } /** * Set the enhancement factor for final-state radiation */ void finalStateRadiationEnhancementFactor(double in) { _finalenhance=in; } //@} /** * Access to set/get the HardTree currently beinging showered */ //@{ /** * The HardTree currently being showered */ tHardTreePtr hardTree() {return _hardtree;} /** * The HardTree currently being showered */ void hardTree(tHardTreePtr in) {_hardtree = in;} //@} /** * Access/set the beam particle for the current initial-state shower */ //@{ /** * Get the beam particle data */ Ptr::const_pointer beamParticle() const { return _beam; } /** * Set the beam particle data */ void setBeamParticle(Ptr::const_pointer in) { _beam=in; } //@} /** * Set/Get the current tree being evolver for inheriting classes */ //@{ /** * Get the tree */ tShowerTreePtr currentTree() { return _currenttree; } /** * Set the tree */ void currentTree(tShowerTreePtr tree) { _currenttree=tree; } //@} /** * Access the maximum number of attempts to generate the shower */ unsigned int maximumTries() const { return _maxtry; } /** * Set/Get the ShowerProgenitor for the current shower */ //@{ /** * Access the progenitor */ ShowerProgenitorPtr progenitor() { return _progenitor; } /** * Set the progenitor */ void progenitor(ShowerProgenitorPtr in) { _progenitor=in; } //@} /** * Calculate the intrinsic \f$p_T\f$. */ virtual void generateIntrinsicpT(vector); /** * Access to the intrinsic \f$p_T\f$ for inheriting classes */ map > & intrinsicpT() { return _intrinsic; } /** * find the maximally allowed pt acc to the hard process. */ void setupMaximumScales(const vector &,XCPtr); /** * find the relevant hard scales for profile scales. */ void setupHardScales(const vector &,XCPtr); /** * Convert the HardTree into an extra shower emission */ void convertHardTree(bool hard,ShowerInteraction type); protected: /** * Find the parton extracted from the incoming particle after ISR */ PPtr findFirstParton(tPPtr seed) const; /** * Fix Remnant connections after ISR */ tPPair remakeRemnant(tPPair oldp); protected: /** * Start the shower of a timelike particle */ virtual bool startTimeLikeShower(ShowerInteraction); /** * Update of the time-like stuff */ void updateHistory(tShowerParticlePtr particle); /** * Start the shower of a spacelike particle */ virtual bool startSpaceLikeShower(PPtr,ShowerInteraction); /** * Start the shower of a spacelike particle */ virtual bool startSpaceLikeDecayShower(const ShowerParticle::EvolutionScales & maxScales, Energy minimumMass,ShowerInteraction); /** * Select the branching for the next time-like emission */ Branching selectTimeLikeBranching(tShowerParticlePtr particle, ShowerInteraction type, HardBranchingPtr branch); /** * Select the branching for the next space-like emission in a decay */ Branching selectSpaceLikeDecayBranching(tShowerParticlePtr particle, const ShowerParticle::EvolutionScales & maxScales, Energy minmass,ShowerInteraction type, HardBranchingPtr branch); /** * Create the timelike child of a branching */ ShowerParticleVector createTimeLikeChildren(tShowerParticlePtr particle, IdList ids); /** * Vetos for the timelike shower */ virtual bool timeLikeVetoed(const Branching &,ShowerParticlePtr); /** * Vetos for the spacelike shower */ virtual bool spaceLikeVetoed(const Branching &,ShowerParticlePtr); /** * Vetos for the spacelike shower */ virtual bool spaceLikeDecayVetoed(const Branching &,ShowerParticlePtr); /** * Only generate the hard emission, for testing only. */ bool hardOnly() const {return _limitEmissions==3;} /** * Check the flags */ void checkFlags(); /** * */ void addFSRUsingDecayPOWHEG(HardTreePtr ISRTree); public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** * The main method which manages the showering of a subprocess. */ virtual tPPair cascade(tSubProPtr sub, XCPtr xcomb); /** * Decay a ShowerTree */ void decay(ShowerTreePtr tree, ShowerDecayMap & decay); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ QTildeShowerHandler & operator=(const QTildeShowerHandler &) = delete; private: /** * Stuff from the ShowerHandler */ //@{ /** * The ShowerTree for the hard process */ ShowerTreePtr hard_; /** * The ShowerTree for the decays */ ShowerDecayMap decay_; /** * The ShowerTrees for which the initial shower */ vector done_; //@} private : /** * Pointer to the splitting generator */ SplittingGeneratorPtr _splittingGenerator; /** * Maximum number of tries to generate the shower of a particular tree */ unsigned int _maxtry; /** * Matrix element correction switch */ unsigned int _meCorrMode; /** * Control of the reconstruction option */ unsigned int _evolutionScheme; /** * If hard veto pT scale is being read-in this determines * whether the read-in value is applied to primary and * secondary (MPI) scatters or just the primary one, with * the usual computation of the veto being performed for * the secondary (MPI) scatters. */ bool _hardVetoReadOption; /** * rms intrinsic pT of Gaussian distribution */ Energy _iptrms; /** * Proportion of inverse quadratic intrinsic pT distribution */ double _beta; /** * Parameter for inverse quadratic: 2*Beta*Gamma/(sqr(Gamma)+sqr(intrinsicpT)) */ Energy _gamma; /** * Upper bound on intrinsic pT for inverse quadratic */ Energy _iptmax; /** * Limit the number of emissions for testing */ unsigned int _limitEmissions; /** * The progenitor of the current shower */ ShowerProgenitorPtr _progenitor; /** * Matrix element */ HwMEBasePtr _hardme; /** * Decayer */ HwDecayerBasePtr _decayme; /** * The ShowerTree currently being showered */ ShowerTreePtr _currenttree; /** * The HardTree currently being showered */ HardTreePtr _hardtree; /** * Radiation enhancement factors for use with the veto algorithm * if needed by the soft matrix element correction */ //@{ /** * Enhancement factor for initial-state radiation */ double _initialenhance; /** * Enhancement factor for final-state radiation */ double _finalenhance; //@} /** * The beam particle data for the current initial-state shower */ Ptr::const_pointer _beam; /** * Storage of the intrinsic \f$p_t\f$ of the particles */ map > _intrinsic; /** * Vetoes */ vector _vetoes; /** * Full Shower Vetoes */ vector _fullShowerVetoes; /** * Number of iterations for reweighting */ unsigned int _nReWeight; /** * Whether or not we are reweighting */ bool _reWeight; /** * number of IS emissions */ unsigned int _nis; /** * Number of FS emissions */ unsigned int _nfs; /** * The option for wqhich interactions to use */ ShowerInteraction interaction_; /** * Truncated shower switch */ bool _trunc_Mode; /** * Mode for the hard emissions */ int _hardEmission; /** * Option for the kernal for soft correlations */ unsigned int _softOpt; /** * Option for hard radiation in POWHEG events */ bool _hardPOWHEG; /** * True if no warnings about incorrect hard emission * mode setting have been issued yet */ static bool _hardEmissionWarn; /** * True if no warnings about missing truncated shower * have been issued yet */ static bool _missingTruncWarn; /** * The relevant hard scale to be used in the profile scales */ Energy muPt; private: /** * Pointer to the various objects */ //@{ /** * Pointer to the KinematicsReconstructor object */ KinematicsReconstructorPtr _reconstructor; /** * Pointer to the PartnerFinder object */ PartnerFinderPtr _partnerfinder; //@} }; } #endif /* HERWIG_QTildeShowerHandler_H */ diff --git a/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc b/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc --- a/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc +++ b/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc @@ -1,243 +1,243 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the HalfHalfZeroEWSplitFn class. // #include "HalfHalfZeroEWSplitFn.h" #include "ThePEG/StandardModel/StandardModelBase.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/ParticleData.h" #include "Herwig/Decay/TwoBodyDecayMatrixElement.h" #include "Herwig/Models/StandardModel/SMFFHVertex.h" using namespace Herwig; IBPtr HalfHalfZeroEWSplitFn::clone() const { return new_ptr(*this); } IBPtr HalfHalfZeroEWSplitFn::fullclone() const { return new_ptr(*this); } void HalfHalfZeroEWSplitFn::persistentOutput(PersistentOStream & os) const { os << ghqq_ << _theSM; } void HalfHalfZeroEWSplitFn::persistentInput(PersistentIStream & is, int) { is >> ghqq_ >> _theSM; } // The following static variable is needed for the type description system in ThePEG. DescribeClass describeHerwigHalfHalfZeroEWSplitFn("Herwig::HalfHalfZeroEWSplitFn", "HwShower.so"); void HalfHalfZeroEWSplitFn::Init() { static ClassDocumentation documentation ("The HalfHalfZeroEWSplitFn class implements the splitting q->qH"); } void HalfHalfZeroEWSplitFn::doinit() { SplittingFunction::doinit(); tcSMPtr sm = generator()->standardModel(); double sw2 = sm->sin2ThetaW(); ghqq_ = 1./sqrt(4.*sw2); _theSM = dynamic_ptr_cast(generator()->standardModel()); } void HalfHalfZeroEWSplitFn::getCouplings(double & gH, const IdList & ids) const { if(abs(ids[2]->id())==ParticleID::h0) { //get quark masses Energy mq; if(abs(ids[0]->id())==ParticleID::c) mq = getParticleData(ParticleID::c)->mass(); else if(abs(ids[0]->id())==ParticleID::b) mq = getParticleData(ParticleID::b)->mass(); else if(abs(ids[0]->id())==ParticleID::t) mq = getParticleData(ParticleID::t)->mass(); Energy mW = getParticleData(ParticleID::Wplus)->mass(); gH = ghqq_*(mq/mW); } else assert(false); } void HalfHalfZeroEWSplitFn::getCouplings(double & gH, const IdList & ids, const Energy2 t) const { if(abs(ids[2]->id())==ParticleID::h0) { //get quark masses Energy mq; if(abs(ids[0]->id())==ParticleID::c) mq = _theSM->mass(t,getParticleData(ParticleID::c)); else if(abs(ids[0]->id())==ParticleID::b) mq = _theSM->mass(t,getParticleData(ParticleID::b)); else if(abs(ids[0]->id())==ParticleID::t) mq = _theSM->mass(t,getParticleData(ParticleID::t)); Energy mW = getParticleData(ParticleID::Wplus)->mass(); //Energy mW = _theSM->mass(t,getParticleData(ParticleID::Wplus)); gH = ghqq_*(mq/mW); } else assert(false); } double HalfHalfZeroEWSplitFn::P(const double z, const Energy2 t, - const IdList &ids, const bool mass, const RhoDMatrix & rho) const { + const IdList &ids, const bool mass, const RhoDMatrix &) const { double gH(0.); getCouplings(gH,ids,t); double val = (1.-z); Energy mq, mH; //get masses if(mass) { mq = ids[0]->mass(); mH = ids[2]->mass(); } else { // to assure the particle mass in non-zero if(abs(ids[0]->id())==ParticleID::c) mq = getParticleData(ParticleID::c)->mass(); else if(abs(ids[0]->id())==ParticleID::b) mq = getParticleData(ParticleID::b)->mass(); else if(abs(ids[0]->id())==ParticleID::t) mq = getParticleData(ParticleID::t)->mass(); mH = getParticleData(ParticleID::h0)->mass(); } val += (4.*sqr(mq) - sqr(mH))/(t*(1. - z)*z); val *= sqr(gH); return colourFactor(ids)*val; } double HalfHalfZeroEWSplitFn::overestimateP(const double z, const IdList & ids) const { double gH(0.); getCouplings(gH,ids); return sqr(gH)*colourFactor(ids)*(1.-z); } double HalfHalfZeroEWSplitFn::ratioP(const double z, const Energy2 t, const IdList & ids, const bool mass, - const RhoDMatrix & rho) const { + const RhoDMatrix & ) const { double gH(0.); getCouplings(gH,ids,t); double val = 1.; Energy mq, mH; if(mass) { mq = ids[0]->mass(); mH = ids[2]->mass(); } else { // to assure the particle mass in non-zero if(abs(ids[0]->id())==ParticleID::c) mq = getParticleData(ParticleID::c)->mass(); else if(abs(ids[0]->id())==ParticleID::b) mq = getParticleData(ParticleID::b)->mass(); else if(abs(ids[0]->id())==ParticleID::t) mq = getParticleData(ParticleID::t)->mass(); mH = getParticleData(ParticleID::h0)->mass(); } val += (4.*sqr(mq) - sqr(mH))/(t*(1. - z)*z); return val; } double HalfHalfZeroEWSplitFn::integOverP(const double z, const IdList & ids, unsigned int PDFfactor) const { double gH(0.); getCouplings(gH,ids); double pre = colourFactor(ids)*sqr(gH); switch (PDFfactor) { case 0: //OverP return pre*(z-sqr(z)/2.); case 1: //OverP/z return pre*(log(z)-z); case 2: //OverP/(1-z) return pre*z; case 3: //OverP/[z(1-z)] return pre*log(z); default: throw Exception() << "HalfHalfZeroEWSplitFn::integOverP() invalid PDFfactor = " << PDFfactor << Exception::runerror; } } double HalfHalfZeroEWSplitFn::invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor) const { double gH(0.); getCouplings(gH,ids); double pre = colourFactor(ids)*sqr(gH); switch (PDFfactor) { case 0: return 1. - sqrt(1. - 2.*r/pre); case 1: //OverP/z case 2: //OverP/(1-z) return r/pre; case 3: //OverP/[z(1-z)] return exp(r/pre); default: throw Exception() << "HalfHalfZeroEWSplitFn::invIntegOverP() invalid PDFfactor = " << PDFfactor << Exception::runerror; } } bool HalfHalfZeroEWSplitFn::accept(const IdList &ids) const { if(ids.size()!=3) return false; if(ids[2]->id()==ParticleID::h0) { if(ids[0]->id()==ids[1]->id() && (ids[0]->id()==4 || ids[0]->id()==5 || ids[0]->id()==6)) return true; } return false; } vector > HalfHalfZeroEWSplitFn::generatePhiForward(const double, const Energy2, const IdList & , const RhoDMatrix &) { // no dependence on the spin density matrix, dependence on off-diagonal terms cancels // and rest = splitting function for Tr(rho)=1 as required by defn return vector >(1,make_pair(0,1.)); } vector > HalfHalfZeroEWSplitFn::generatePhiBackward(const double, const Energy2, const IdList & , const RhoDMatrix &) { // no dependence on the spin density matrix, dependence on off-diagonal terms cancels // and rest = splitting function for Tr(rho)=1 as required by defn return vector >(1,make_pair(0,1.)); } DecayMEPtr HalfHalfZeroEWSplitFn::matrixElement(const double z, const Energy2 t, const IdList & ids, const double phi, bool) { // calculate the kernal DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin0))); //get masses Energy mq, mH; if(abs(ids[0]->id())==ParticleID::c) mq = getParticleData(ParticleID::c)->mass(); else if(abs(ids[0]->id())==ParticleID::b) mq = getParticleData(ParticleID::b)->mass(); else if(abs(ids[0]->id())==ParticleID::t) mq = getParticleData(ParticleID::t)->mass(); mH = getParticleData(ParticleID::h0)->mass(); double gH(0.); Energy2 tC = t/(z*(1-z)); getCouplings(gH,ids,tC); double mqt = mq/sqrt(tC); double mHt = mH/sqrt(tC); double num1 = gH*(1.+z)*mqt; double num2 = gH*sqrt(-sqr(mqt)*(1.-z) - sqr(mHt)*z + z*(1.-z)*(sqr(mqt)+z*(1.-z))); //watch this double dnum = sqrt(2.)*sqrt((1.-z)*sqr(z)); Complex phase = exp(Complex(0.,1.)*phi); Complex cphase = conj(phase); (*kernal)(0,0,0) = num1/dnum; (*kernal)(0,1,0) = cphase*num2/dnum; (*kernal)(1,0,0) = -phase*num2/dnum; (*kernal)(1,1,0) = num1/dnum; // return the answer return kernal; }