diff --git a/Shower/QTilde/Kinematics/Decay_QTildeShowerKinematics1to2.cc b/Shower/QTilde/Kinematics/Decay_QTildeShowerKinematics1to2.cc --- a/Shower/QTilde/Kinematics/Decay_QTildeShowerKinematics1to2.cc +++ b/Shower/QTilde/Kinematics/Decay_QTildeShowerKinematics1to2.cc @@ -1,114 +1,115 @@ // -*- C++ -*- // // Decay_QTildeShowerKinematics1to2.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 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 Decay_QTildeShowerKinematics1to2 class. // #include "Decay_QTildeShowerKinematics1to2.h" #include "ThePEG/PDT/EnumParticles.h" #include "Herwig/Shower/QTilde/SplittingFunctions/SplittingFunction.h" #include "Herwig/Shower/QTilde/Base/ShowerParticle.h" #include #include "Herwig/Shower/ShowerHandler.h" #include "Herwig/Shower/QTilde/Base/ShowerVertex.h" using namespace Herwig; void Decay_QTildeShowerKinematics1to2:: updateChildren(const tShowerParticlePtr parent, const ShowerParticleVector & children, ShowerPartnerType partnerType) const { assert(children.size() == 2); // calculate the scales splittingFn()->evaluateDecayScales(partnerType,scale(),z(),parent, children[0],children[1]); // set the maximum virtual masses IdList ids(3); ids[0] = parent->dataPtr(); ids[1] = children[0]->dataPtr(); ids[2] = children[1]->dataPtr(); const vector & virtualMasses = SudakovFormFactor()->virtualMasses(ids); Energy2 q2 = sqr(virtualMasses[0])-(1.-z())*sqr(scale()); children[0]->virtualMass(sqrt(q2)); // determine alphas of children according to interpretation of z const ShowerParticle::Parameters & params = parent->showerParameters(); ShowerParticle::Parameters & child0 = children[0]->showerParameters(); ShowerParticle::Parameters & child1 = children[1]->showerParameters(); child0.alpha = z() * params.alpha; child1.alpha = (1.-z()) * params.alpha; child0.ptx = pT() * cos(phi()) + z()* params.ptx; child0.pty = pT() * sin(phi()) + z()* params.pty; child0.pt = sqrt( sqr(child0.ptx) + sqr(child0.pty) ); child1.ptx = -pT() * cos(phi()) + (1.-z()) * params.ptx; child1.pty = -pT() * sin(phi()) + (1.-z()) * params.pty; child1.pt = sqrt( sqr(child1.ptx) + sqr(child1.pty) ); // set up the colour connections splittingFn()->colourConnection(parent,children[0],children[1],partnerType,false); // make the products children of the parent parent->addChild(children[0]); parent->addChild(children[1]); // set the momenta of the children for(ShowerParticleVector::const_iterator pit=children.begin(); pit!=children.end();++pit) { (**pit).showerBasis(parent->showerBasis(),true); (**pit).setShowerMomentum(true); } } void Decay_QTildeShowerKinematics1to2:: reconstructParent( const tShowerParticlePtr, const ParticleVector &) const { throw Exception() << "Decay_QTildeShowerKinematics1to2::reconstructParent not implemented" << Exception::abortnow; } void Decay_QTildeShowerKinematics1to2:: reconstructLast(const tShowerParticlePtr last, Energy mass) const { // set beta component and consequently all missing data from that, // using the nominal (i.e. PDT) mass. Energy theMass =ZERO; if(!(mass > ZERO) && ShowerHandler::currentHandler()->retConstituentMasses()) theMass = last->data().constituentMass(); else theMass = mass > ZERO ? mass : last->data().mass(); last->showerParameters().beta= (sqr(theMass) + sqr(last->showerParameters().pt) - sqr( last->showerParameters().alpha )*last->showerBasis()->pVector().m2()) / ( 2.*last->showerParameters().alpha*last->showerBasis()->p_dot_n() ); // set that new momentum last->set5Momentum( last->showerBasis()->sudakov2Momentum( last->showerParameters().alpha, last->showerParameters().beta, last->showerParameters().ptx, last->showerParameters().pty) ); } void Decay_QTildeShowerKinematics1to2::updateParent(const tShowerParticlePtr parent, const ShowerParticleVector & children, + unsigned int , ShowerPartnerType) const { IdList ids(3); ids[0] = parent->dataPtr(); ids[1] = children[0]->dataPtr(); ids[2] = children[1]->dataPtr(); const vector & virtualMasses = SudakovFormFactor()->virtualMasses(ids); children[0]->virtualMass(sqrt(sqr(virtualMasses[0])-(1.-z())*sqr(scale()))); if(children[1]->children().empty()) children[1]->virtualMass(virtualMasses[2]); // compute the new pT of the branching Energy2 pt2=(1.-z())*(z()*sqr(virtualMasses[0])-sqr(children[0]->virtualMass())) -z()*sqr(children[1]->virtualMass()); if(pt2>ZERO) { pT(sqrt(pt2)); } else { parent->virtualMass(ZERO); } } diff --git a/Shower/QTilde/Kinematics/Decay_QTildeShowerKinematics1to2.h b/Shower/QTilde/Kinematics/Decay_QTildeShowerKinematics1to2.h --- a/Shower/QTilde/Kinematics/Decay_QTildeShowerKinematics1to2.h +++ b/Shower/QTilde/Kinematics/Decay_QTildeShowerKinematics1to2.h @@ -1,113 +1,114 @@ // -*- C++ -*- // // Decay_QTildeShowerKinematics1to2.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 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. // #ifndef HERWIG_Decay_QTildeShowerKinematics1to2_H #define HERWIG_Decay_QTildeShowerKinematics1to2_H // // This is the declaration of the Decay_QTildeShowerKinematics1to2 class. // #include "Herwig/Shower/QTilde/Kinematics/ShowerKinematics.h" namespace Herwig { using namespace ThePEG; /** \ingroup Shower * * This (concrete) class provides the specific decay shower * kinematics information. * * @see ShowerKinematics * @see IS_QTildeShowerKinematics1to2 * @see FS_QTildeShowerKinematics1to2 * @see KinematicsReconstructor * */ class Decay_QTildeShowerKinematics1to2: public ShowerKinematics { public: /** * Default constructor */ Decay_QTildeShowerKinematics1to2() = default; /** * The constructor. */ Decay_QTildeShowerKinematics1to2(Energy scale, double z, double phi, Energy pt, tSudakovPtr sud) : ShowerKinematics(scale,z,phi,pt,sud) {} public: /** * The updateChildren, updateParent and updateLast * members to update the values of the \f$\alpha\f$ and * \f$p_\perp\f$ variables during the shower evolution. */ //@{ /** * Along with the showering evolution --- going forward for * time-like (forward) evolution, and going backward for space-like * (backward) evolution --- the kinematical variables of the * branching products are calculated and updated from the knowledge * of the parent kinematics. This method is used by the * ForwardShowerEvolver. * @param parent The branching particle * @param children The particles produced in the branching * @param partnerType The type of evolution partner */ virtual void updateChildren( const tShowerParticlePtr parent, const ShowerParticleVector & children, ShowerPartnerType partnerType) const; /** * Update the parent Kinematics from the knowledge of the kinematics * of the children. This method will be used by the * KinematicsReconstructor. */ virtual void reconstructParent( const tShowerParticlePtr parent, const ParticleVector & children ) const; /** * Update the kinematical data of a particle when a reconstruction * fixpoint was found. This will highly depend on the kind of * kinematics chosen and will be defined in the inherited concrete * classes. This method will be used by the KinematicsReconstructor. * @param last The particle to update * @param mass The mass to be used, if less than zero on-shell */ virtual void reconstructLast(const tShowerParticlePtr last, Energy mass=-1.*GeV) const; /** * Update the parent Kinematics from the knowledge of the kinematics * of the children. This method will be used by the KinematicsReconstructor. * @param parent The parent * @param children The children * @param partnerType The type of evolution partner */ virtual void updateParent(const tShowerParticlePtr parent, const ShowerParticleVector & children, + unsigned int pTscheme, ShowerPartnerType partnerType) const; //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ Decay_QTildeShowerKinematics1to2 & operator=(const Decay_QTildeShowerKinematics1to2 &) = delete; }; } #endif /* HERWIG_Decay_QTildeShowerKinematics1to2_H */ diff --git a/Shower/QTilde/Kinematics/FS_QTildeShowerKinematics1to2.cc b/Shower/QTilde/Kinematics/FS_QTildeShowerKinematics1to2.cc --- a/Shower/QTilde/Kinematics/FS_QTildeShowerKinematics1to2.cc +++ b/Shower/QTilde/Kinematics/FS_QTildeShowerKinematics1to2.cc @@ -1,196 +1,211 @@ // -*- C++ -*- // // FS_QTildeShowerKinematics1to2.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 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 FS_QTildeShowerKinematics1to2 class. // #include "FS_QTildeShowerKinematics1to2.h" #include "ThePEG/PDT/EnumParticles.h" #include "Herwig/Shower/QTilde/SplittingFunctions/SplittingFunction.h" #include "Herwig/Shower/QTilde/Base/ShowerParticle.h" #include "ThePEG/Utilities/Debug.h" #include "Herwig/Shower/QTilde/QTildeShowerHandler.h" #include "Herwig/Shower/QTilde/Base/PartnerFinder.h" #include "Herwig/Shower/QTilde/Kinematics/KinematicsReconstructor.h" #include "Herwig/Shower/QTilde/Kinematics/KinematicHelpers.h" #include "Herwig/Shower/QTilde/Base/ShowerVertex.h" using namespace Herwig; void FS_QTildeShowerKinematics1to2:: updateParameters(tShowerParticlePtr theParent, tShowerParticlePtr theChild0, tShowerParticlePtr theChild1, bool setAlpha) const { const ShowerParticle::Parameters & parent = theParent->showerParameters(); ShowerParticle::Parameters & child0 = theChild0->showerParameters(); ShowerParticle::Parameters & child1 = theChild1->showerParameters(); // determine alphas of children according to interpretation of z if ( setAlpha ) { child0.alpha = z() * parent.alpha; child1.alpha = (1.-z()) * parent.alpha; } // set the values double cphi = cos(phi()); double sphi = sin(phi()); child0.ptx = pT() * cphi + z() * parent.ptx; child0.pty = pT() * sphi + z() * parent.pty; child0.pt = sqrt( sqr(child0.ptx) + sqr(child0.pty) ); child1.ptx = -pT() * cphi + (1.-z())* parent.ptx; child1.pty = -pT() * sphi + (1.-z())* parent.pty; child1.pt = sqrt( sqr(child1.ptx) + sqr(child1.pty) ); } void FS_QTildeShowerKinematics1to2:: updateChildren(const tShowerParticlePtr parent, const ShowerParticleVector & children, ShowerPartnerType partnerType) const { assert(children.size()==2); // calculate the scales splittingFn()->evaluateFinalStateScales(partnerType,scale(),z(),parent, children[0],children[1]); // update the parameters updateParameters(parent, children[0], children[1], true); // set up the colour connections splittingFn()->colourConnection(parent,children[0],children[1],partnerType,false); // make the products children of the parent parent->addChild(children[0]); parent->addChild(children[1]); // set the momenta of the children for(ShowerParticleVector::const_iterator pit=children.begin(); pit!=children.end();++pit) { (**pit).showerBasis(parent->showerBasis(),true); (**pit).setShowerMomentum(true); } // sort out the helicity stuff if(! dynamic_ptr_cast(ShowerHandler::currentHandler())->correlations()) return; SpinPtr pspin(parent->spinInfo()); if(!pspin || !dynamic_ptr_cast(ShowerHandler::currentHandler())->spinCorrelations() ) return; Energy2 t = sqr(scale())*z()*(1.-z()); IdList ids; ids.push_back(parent->dataPtr()); ids.push_back(children[0]->dataPtr()); ids.push_back(children[1]->dataPtr()); // create the vertex SVertexPtr vertex(new_ptr(ShowerVertex())); // set the matrix element vertex->ME(splittingFn()->matrixElement(z(),t,ids,phi(),true)); RhoDMatrix mapping; SpinPtr inspin; bool needMapping = parent->getMapping(inspin,mapping); if(needMapping) vertex->incomingBasisTransform(mapping); // set the incoming particle for the vertex parent->spinInfo()->decayVertex(vertex); for(ShowerParticleVector::const_iterator pit=children.begin(); pit!=children.end();++pit) { // construct the spin info for the children (**pit).constructSpinInfo(true); // connect the spinInfo object to the vertex (*pit)->spinInfo()->productionVertex(vertex); } } void FS_QTildeShowerKinematics1to2:: reconstructParent(const tShowerParticlePtr parent, const ParticleVector & children ) const { assert(children.size() == 2); ShowerParticlePtr c1 = dynamic_ptr_cast(children[0]); ShowerParticlePtr c2 = dynamic_ptr_cast(children[1]); parent->showerParameters().beta= c1->showerParameters().beta + c2->showerParameters().beta; Lorentz5Momentum pnew = c1->momentum() + c2->momentum(); Energy2 m2 = sqr(pT())/z()/(1.-z()) + sqr(c1->mass())/z() + sqr(c2->mass())/(1.-z()); pnew.setMass(sqrt(m2)); parent->set5Momentum( pnew ); } void FS_QTildeShowerKinematics1to2::reconstructLast(const tShowerParticlePtr last, Energy mass) const { // set beta component and consequently all missing data from that, // using the nominal (i.e. PDT) mass. Energy theMass =ZERO; if(!(mass > ZERO) && ShowerHandler::currentHandler()->retConstituentMasses()) theMass = last->data().constituentMass(); else theMass = mass > ZERO ? mass : last->data().mass(); Lorentz5Momentum pVector = last->showerBasis()->pVector(); ShowerParticle::Parameters & lastParam = last->showerParameters(); Energy2 denom = 2. * lastParam.alpha * last->showerBasis()->p_dot_n(); if(abs(denom)/(sqr(pVector.e())+pVector.rho2())<1e-10) { throw KinematicsReconstructionVeto(); } lastParam.beta = ( sqr(theMass) + sqr(lastParam.pt) - sqr(lastParam.alpha) * pVector.m2() ) / denom; // set that new momentum Lorentz5Momentum newMomentum = last->showerBasis()-> sudakov2Momentum( lastParam.alpha, lastParam.beta, lastParam.ptx , lastParam.pty); newMomentum.setMass(theMass); newMomentum.rescaleEnergy(); if(last->data().stable()) { last->set5Momentum( newMomentum ); } else { last->boost(last->momentum().findBoostToCM()); last->boost(newMomentum.boostVector()); } } void FS_QTildeShowerKinematics1to2::updateParent(const tShowerParticlePtr parent, const ShowerParticleVector & children, + unsigned int pTscheme, ShowerPartnerType) const { IdList ids(3); ids[0] = parent->dataPtr(); ids[1] = children[0]->dataPtr(); ids[2] = children[1]->dataPtr(); const vector & virtualMasses = SudakovFormFactor()->virtualMasses(ids); if(children[0]->children().empty()) children[0]->virtualMass(virtualMasses[1]); if(children[1]->children().empty()) children[1]->virtualMass(virtualMasses[2]); // compute the new pT of the branching Energy2 m02 = (ids[0]->id()!=ParticleID::g && ids[0]->id()!=ParticleID::gamma) ? sqr(virtualMasses[0]) : Energy2(); - - Energy2 pt2 = QTildeKinematics::pT2_FSR(sqr(scale()), z(), m02, - sqr(children[0]->virtualMass()), sqr(children[1]->virtualMass()), - sqr(children[0]->virtualMass()), sqr(children[1]->virtualMass())); - + Energy2 pt2; + if(pTscheme==0) { + pt2 = QTildeKinematics::pT2_FSR(sqr(scale()), z(), m02, + sqr(virtualMasses[1]) ,sqr(virtualMasses[2]), + sqr(virtualMasses[1]) ,sqr(virtualMasses[2])); + } + else if(pTscheme==1) { + pt2 = QTildeKinematics::pT2_FSR(sqr(scale()), z(), m02, + sqr(virtualMasses[1]) ,sqr(virtualMasses[2]), + sqr(children[0]->virtualMass()), sqr(children[1]->virtualMass())); + } + else if(pTscheme==2) { + pt2 = QTildeKinematics::pT2_FSR(sqr(scale()), z(), m02, + sqr(children[0]->virtualMass()), sqr(children[1]->virtualMass()), + sqr(children[0]->virtualMass()), sqr(children[1]->virtualMass())); + } + else + assert(false); + if(pt2>ZERO) { pT(sqrt(pt2)); } else { pt2=ZERO; pT(ZERO); } Energy2 q2 = QTildeKinematics::q2_FSR( pt2, z(), sqr(children[0]->virtualMass()), sqr(children[1]->virtualMass()) ); parent->virtualMass(sqrt(q2)); } void FS_QTildeShowerKinematics1to2:: resetChildren(const tShowerParticlePtr parent, const ShowerParticleVector & children) const { updateParameters(parent, children[0], children[1], false); for(unsigned int ix=0;ixchildren().empty()) continue; ShowerParticleVector newChildren; for(unsigned int iy=0;iychildren().size();++iy) newChildren.push_back(dynamic_ptr_cast (children[ix]->children()[iy])); children[ix]->showerKinematics()->resetChildren(children[ix],newChildren); } } diff --git a/Shower/QTilde/Kinematics/FS_QTildeShowerKinematics1to2.h b/Shower/QTilde/Kinematics/FS_QTildeShowerKinematics1to2.h --- a/Shower/QTilde/Kinematics/FS_QTildeShowerKinematics1to2.h +++ b/Shower/QTilde/Kinematics/FS_QTildeShowerKinematics1to2.h @@ -1,123 +1,124 @@ // -*- C++ -*- // // FS_QTildeShowerKinematics1to2.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 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. // #ifndef HERWIG_FS_QTildeShowerKinematics1to2_H #define HERWIG_FS_QTildeShowerKinematics1to2_H // // This is the declaration of the FS_QTildeShowerKinematics1to2 class. // #include "Herwig/Shower/QTilde/Kinematics/ShowerKinematics.h" namespace Herwig { using namespace ThePEG; /** \ingroup Shower * * This (concrete) class provides the specific Final State shower * kinematics information. * * @see ShowerKinematics * @see IS_QTildeShowerKinematics1to2 * @see Decay_QTildeShowerKinematics1to2 * @see KinematicsReconstructor */ class FS_QTildeShowerKinematics1to2: public ShowerKinematics { public: /** * Default constructor */ FS_QTildeShowerKinematics1to2() = default; /** * The constructor. */ FS_QTildeShowerKinematics1to2(Energy scale, double z, double phi, Energy pt, tSudakovPtr sud) : ShowerKinematics(scale,z,phi,pt,sud) {} /** * The updateChildren, updateParent and updateLast * members to update the values of the \f$\alpha\f$ and * \f$p_\perp\f$ variables during the shower evolution. */ //@{ /** * Along with the showering evolution --- going forward for * time-like (forward) evolution, and going backward for space-like * (backward) evolution --- the kinematical variables of the * branching products are calculated and updated from the knowledge * of the parent kinematics. This method is used by the * ForwardShowerEvolver. * @param parent The branching particle * @param children The particles produced in the branching * @param partnerType The type of evolution partner */ private: void updateParameters(tShowerParticlePtr theParent, tShowerParticlePtr theChild0, tShowerParticlePtr theChild1, bool setAlpha) const; public: virtual void updateChildren( const tShowerParticlePtr parent, const ShowerParticleVector & children, ShowerPartnerType partnerType) const; virtual void resetChildren( const tShowerParticlePtr parent, const ShowerParticleVector & children) const; /** * Update the parent Kinematics from the knowledge of the kinematics * of the children. This method will be used by the KinematicsReconstructor. * @param parent The parent * @param children The children * @param partnerType The type of evolution partner */ virtual void updateParent(const tShowerParticlePtr parent, const ShowerParticleVector & children, + unsigned int pTscheme, ShowerPartnerType partnerType) const; /** * Update the parent Kinematics from the knowledge of the kinematics * of the children. This method will be used by the * KinematicsReconstructor. */ virtual void reconstructParent( const tShowerParticlePtr parent, const ParticleVector & children ) const; /** * Update the kinematical data of a particle when a reconstruction * fixpoint was found. This will highly depend on the kind of * kinematics chosen and will be defined in the inherited concrete * classes. This method will be used by the KinematicsReconstructor. * @param last The particle to update * @param mass The mass to be used, if less than zero on-shell */ virtual void reconstructLast(const tShowerParticlePtr last, Energy mass=-1.*GeV) const; //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ FS_QTildeShowerKinematics1to2 & operator=(const FS_QTildeShowerKinematics1to2 &) = delete; }; } #endif /* HERWIG_FS_QTildeShowerKinematics1to2_H */ diff --git a/Shower/QTilde/Kinematics/IS_QTildeShowerKinematics1to2.cc b/Shower/QTilde/Kinematics/IS_QTildeShowerKinematics1to2.cc --- a/Shower/QTilde/Kinematics/IS_QTildeShowerKinematics1to2.cc +++ b/Shower/QTilde/Kinematics/IS_QTildeShowerKinematics1to2.cc @@ -1,149 +1,150 @@ // -*- C++ -*- // // IS_QTildeShowerKinematics1to2.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 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 IS_QTildeShowerKinematics1to2 class. // #include "IS_QTildeShowerKinematics1to2.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "Herwig/Shower/QTilde/Base/ShowerParticle.h" #include "ThePEG/Utilities/Debug.h" #include "Herwig/Shower/QTilde/QTildeShowerHandler.h" #include "Herwig/Shower/QTilde/Base/PartnerFinder.h" #include "Herwig/Shower/QTilde/Kinematics/KinematicsReconstructor.h" #include "Herwig/Shower/QTilde/Base/ShowerVertex.h" #include using namespace Herwig; void IS_QTildeShowerKinematics1to2:: updateChildren( const tShowerParticlePtr theParent, const ShowerParticleVector & children, ShowerPartnerType) const { const ShowerParticle::Parameters & parent = theParent->showerParameters(); ShowerParticle::Parameters & child0 = children[0]->showerParameters(); ShowerParticle::Parameters & child1 = children[1]->showerParameters(); double cphi = cos(phi()); double sphi = sin(phi()); child1.alpha = (1.-z()) * parent.alpha; child1.ptx = (1.-z()) * parent.ptx - cphi * pT(); child1.pty = (1.-z()) * parent.pty - sphi * pT(); child1.pt = sqrt( sqr(child1.ptx) + sqr(child1.pty) ); // space-like child child0.alpha = parent.alpha - child1.alpha; child0.beta = parent.beta - child1.beta; child0.ptx = parent.ptx - child1.ptx; child0.pty = parent.pty - child1.pty; } void IS_QTildeShowerKinematics1to2:: updateParent(const tShowerParticlePtr parent, const ShowerParticleVector & children, + unsigned int , ShowerPartnerType partnerType) const { // calculate the scales splittingFn()->evaluateInitialStateScales(partnerType,scale(),z(),parent, children[0],children[1]); // set proper colour connections splittingFn()->colourConnection(parent,children[0],children[1], partnerType,true); // set proper parent/child relationships parent->addChild(children[0]); parent->addChild(children[1]); parent->x(children[0]->x()/z()); // sort out the helicity stuff // construct the spin info for parent and timelike child // temporary assignment of shower parameters to calculate correlations parent->showerParameters().alpha = parent->x(); children[1]->showerParameters().alpha = (1.-z()) * parent->x(); children[1]->showerParameters().ptx = - cos(phi()) * pT(); children[1]->showerParameters().pty = - sin(phi()) * pT(); children[1]->showerParameters().pt = pT(); parent ->showerBasis(children[0]->showerBasis(),true); children[1]->showerBasis(children[0]->showerBasis(),true); parent ->setShowerMomentum(false); children[1]->setShowerMomentum(true); if(! dynamic_ptr_cast(ShowerHandler::currentHandler())->correlations()) return; SpinPtr pspin(children[0]->spinInfo()); if(!pspin || !dynamic_ptr_cast(ShowerHandler::currentHandler())->spinCorrelations() ) return; // compute the matrix element for spin correlations IdList ids; ids.push_back(parent->dataPtr()); ids.push_back(children[0]->dataPtr()); ids.push_back(children[1]->dataPtr()); Energy2 t = (1.-z())*sqr(scale())/z(); // create the vertex SVertexPtr vertex(new_ptr(ShowerVertex())); // set the matrix element vertex->ME(splittingFn()->matrixElement(z(),t,ids,phi(),false)); // set the incoming particle for the vertex // (in reality the first child as going backwards) pspin->decayVertex(vertex); // construct the spin infos parent ->constructSpinInfo(false); children[1]->constructSpinInfo(true); // connect the spinInfo objects to the vertex parent ->spinInfo()->productionVertex(vertex); children[1]->spinInfo()->productionVertex(vertex); } void IS_QTildeShowerKinematics1to2:: reconstructParent(const tShowerParticlePtr parent, const ParticleVector & children ) const { PPtr c1 = children[0]; ShowerParticlePtr c2 = dynamic_ptr_cast(children[1]); ShowerParticle::Parameters & c2param = c2->showerParameters(); // get shower variables from 1st child in order to keep notation // parent->(c1, c2) clean even though the splitting was initiated // from c1. The name updateParent is still referring to the // timelike branching though. // on-shell child auto m= ShowerHandler::currentHandler()->retConstituentMasses()? c2->data().constituentMass(): c2->data().mass(); c2param.beta = 0.5*( sqr(m) + sqr(c2param.pt) ) / ( c2param.alpha * parent->showerBasis()->p_dot_n() ); Lorentz5Momentum pnew = parent->showerBasis()-> sudakov2Momentum(c2param.alpha, c2param.beta, c2param.ptx , c2param.pty); pnew.setMass(m); pnew.rescaleEnergy(); c2->set5Momentum( pnew ); // spacelike child Lorentz5Momentum pc1(parent->momentum() - c2->momentum()); pc1.rescaleMass(); c1->set5Momentum(pc1); } void IS_QTildeShowerKinematics1to2:: updateLast( const tShowerParticlePtr theLast,Energy px,Energy py) const { if(theLast->isFinalState()) return; ShowerParticle::Parameters & last = theLast->showerParameters(); Lorentz5Momentum pVector = theLast->showerBasis()->pVector(); Energy2 pt2 = sqr(px) + sqr(py); last.alpha = theLast->x(); last.beta = 0.5 * pt2 / last.alpha / theLast->showerBasis()->p_dot_n(); last.ptx = ZERO; last.pty = ZERO; last.pt = ZERO; // momentum Lorentz5Momentum ntemp = Lorentz5Momentum(ZERO,-pVector.vect()); double beta = 0.5 * pt2 / last.alpha / ( pVector * ntemp); Lorentz5Momentum plast = Lorentz5Momentum( (pVector.z()>ZERO ? px : -px), py, ZERO, ZERO) + theLast->x() * pVector + beta * ntemp; plast.rescaleMass(); theLast->set5Momentum(plast); } diff --git a/Shower/QTilde/Kinematics/IS_QTildeShowerKinematics1to2.h b/Shower/QTilde/Kinematics/IS_QTildeShowerKinematics1to2.h --- a/Shower/QTilde/Kinematics/IS_QTildeShowerKinematics1to2.h +++ b/Shower/QTilde/Kinematics/IS_QTildeShowerKinematics1to2.h @@ -1,117 +1,118 @@ // -*- C++ -*- // // IS_QTildeShowerKinematics1to2.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 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. // #ifndef HERWIG_IS_QTildeShowerKinematics1to2_H #define HERWIG_IS_QTildeShowerKinematics1to2_H // // This is the declaration of the IS_QTildeShowerKinematics1to2 class. // #include "Herwig/Shower/QTilde/Kinematics/ShowerKinematics.h" namespace Herwig { using namespace ThePEG; /** \ingroup Shower * * This (concrete) class provides the specific Intial State shower * kinematics information. * * @see ShowerKinematics * @see FS_QTildeShowerKinematics1to2 * @see Decay_QTildeShowerKinematics1to2 * @see KinematicsReconstructor */ class IS_QTildeShowerKinematics1to2: public ShowerKinematics { public: /** @name Standard constructors and destructors. */ //@{ /** * Construct in terms of the basis states */ IS_QTildeShowerKinematics1to2()= default; /** * The default constructor. */ IS_QTildeShowerKinematics1to2(Energy scale, double z, double phi, Energy pt, tSudakovPtr sud) : ShowerKinematics(scale,z,phi,pt,sud) {} //@} public: /** * The updateChildren, updateParent and updateLast * members to update the values of the \f$\alpha\f$ and * \f$p_\perp\f$ variables during the shower evolution. */ //@{ /** * Along with the showering evolution --- going forward for * time-like (forward) evolution, and going backward for space-like * (backward) evolution --- the kinematical variables of the * branching products are calculated and updated from the knowledge * of the parent kinematics. This method is used by the * ForwardShowerEvolver. * @param parent The branching particle * @param children The particles produced in the branching * @param partnerType The type of evolution partner */ virtual void updateChildren( const tShowerParticlePtr parent, const ShowerParticleVector & children, ShowerPartnerType partnerType) const; /** * Update the parent Kinematics from the knowledge of the kinematics * of the children. This method will be used by the * KinematicsReconstructor. * @param parent The branching particle * @param children The particles produced in the branching * @param partnerType The type of evolution partner */ virtual void updateParent( const tShowerParticlePtr parent, const ShowerParticleVector & children, + unsigned int pTscheme, ShowerPartnerType partnerType) const; /** * Update the parent Kinematics from the knowledge of the kinematics * of the children. This method will be used by the * KinematicsReconstructor. */ virtual void reconstructParent( const tShowerParticlePtr parent, const ParticleVector & children ) const; /** * Update the kinematical data of a particle when a reconstruction * fixpoint was found. This will highly depend on the kind of * kinematics chosen and will be defined in the inherited concrete * classes. This method will be used by the KinematicsReconstructor. * @param theLast The particle. * @param px The \f$x\f$ component of the \f$p_T\f$. * @param py The \f$y\f$ component of the \f$p_T\f$. */ virtual void updateLast(const tShowerParticlePtr theLast, Energy px, Energy py) const; //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ IS_QTildeShowerKinematics1to2 & operator=(const IS_QTildeShowerKinematics1to2 &) = delete; }; } #endif /* HERWIG_IS_QTildeShowerKinematics1to2_H */ 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,69 +1,70 @@ // -*- C++ -*- // // ShowerKinematics.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 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 &, 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/Kinematics/ShowerKinematics.h b/Shower/QTilde/Kinematics/ShowerKinematics.h --- a/Shower/QTilde/Kinematics/ShowerKinematics.h +++ b/Shower/QTilde/Kinematics/ShowerKinematics.h @@ -1,261 +1,262 @@ // -*- C++ -*- // // ShowerKinematics.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 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. // #ifndef HERWIG_ShowerKinematics_H #define HERWIG_ShowerKinematics_H // // This is the declaration of the ShowerKinematics class. // #include "Herwig/Shower/QTilde/ShowerConfig.h" #include "ThePEG/Config/ThePEG.h" #include "Herwig/Shower/QTilde/SplittingFunctions/SudakovFormFactor.h" #include "Herwig/Shower/QTilde/Kinematics/ShowerKinematics.fh" namespace Herwig { using namespace ThePEG; /**\ingroup Shower * * This is the abstract base class from which all other shower * kinematics classes derive. The main purpose of the * shower kinematics classes is to allow the reconstruction * of jet masses, at the end of the showering (indeed, for * multi-scale showering, at the end of each scale-range evolution). * This is necessary for the kinematics reshuffling * in order to compensate the recoil of the emissions. * The KinematicsReconstructor class is in * charge of this job, and which is the main "user" of * ShowerKinematics and its derived classes. * How this is done depends on the choice of kinematics variables * and whether the jet is time-like (forward evolved) or * space-like (backward evolved), whereas the class ShowerKinematics * describes only the common features which are independent by them. * * In general there are a number of methods specific to a shower approach * * @see KinematicsReconstructor */ class ShowerKinematics: public Base { public: /** * The default constructor. */ ShowerKinematics() : Base(), _scale(), _z( 0.0 ), _phi( 0.0 ), _pt(), _sudakov() {} /** * The default constructor. */ ShowerKinematics(Energy scale, double z, double phi, Energy pt, tSudakovPtr sud) : Base(), _scale(scale), _z(z), _phi(phi), _pt(pt), _sudakov(sud) {} /** * The updateChildren and updateParent * members to update the values of the \f$\alpha\f$ and * \f$p_\perp\f$ variables during the shower evolution. */ //@{ /** * Along with the showering evolution --- going forward for * time-like (forward) evolution, and going backward for space-like * (backward) evolution --- the kinematical variables of the * branching products are calculated and updated from the knowledge * of the parent kinematics. * @param parent The parent * @param children The children * @param partnerType The type of evolution partner */ virtual void updateChildren(const tShowerParticlePtr parent, const ShowerParticleVector & children, ShowerPartnerType partnerType) const; virtual void resetChildren( const tShowerParticlePtr parent, const ShowerParticleVector & children) const; /** * Update the parent Kinematics from the knowledge of the kinematics * of the children. This method will be used by the KinematicsReconstructor. * @param parent The parent * @param children The children * @param partnerType The type of evolution partner */ virtual void updateParent(const tShowerParticlePtr parent, const ShowerParticleVector & children, + unsigned int pTscheme, ShowerPartnerType partnerType) const; /** * Update the kinematical data of a particle when a reconstruction * fixpoint was found. This will highly depend on the kind of * kinematics chosen and will be defined in the inherited concrete * classes. This method will be used by the KinematicsReconstructor. * @param last The particle. * @param px The \f$x\f$ component of the \f$p_T\f$. * @param py The \f$y\f$ component of the \f$p_T\f$. */ virtual void updateLast(const tShowerParticlePtr last, Energy px, Energy py) const; //@} /** * The reconstructLast, reconstructChildren and reconstructParent members * are used during the reconstruction */ //@{ /** * Along with the showering evolution --- going forward for * time-like (forward) evolution, and going backward for space-like * (backward) evolution --- the kinematical variables of the * branching products are calculated and updated from the knowledge * of the parent kinematics. * @param parent The parent * @param children The children */ virtual void reconstructChildren(const tShowerParticlePtr parent, const ShowerParticleVector & children) const; /** * Reconstruct the parent Kinematics from the knowledge of the kinematics * of the children. This method will be used by the KinematicsReconstructor. * @param parent The parent * @param children The children */ virtual void reconstructParent(const tShowerParticlePtr parent, const ParticleVector & children) const; /** * Update the kinematical data of a particle when a reconstruction * fixpoint was found. This will highly depend on the kind of * kinematics chosen and will be defined in the inherited concrete * classes. This method will be used by the KinematicsReconstructor. * @param last The particle. * @param mass The mass to be used, if less than zero on-shell */ virtual void reconstructLast(const tShowerParticlePtr last, Energy mass=-1.*GeV) const; //@} public: /** * Set/Get methods for the kinematic variables */ //@{ /** * Access the scale of the splitting. */ Energy scale() const { return _scale; } /** * Set the scale of the splitting. */ void scale(const Energy in) { _scale=in; } /** * Access the energy fraction, \f$z\f$. */ double z() const { return _z; } /** * Set the energy fraction, \f$z\f$. */ void z(const double in) { _z=in; } /** * Access the azimuthal angle, \f$\phi\f$. */ double phi() const { return _phi; } /** * Set the azimuthal angle, \f$\phi\f$. */ void phi(const double in) { _phi=in; } /** * Access the relative \f$p_T\f$ for the branching */ Energy pT() const { return _pt; } /** * Set the relative \f$p_T\f$ for the branching */ void pT(const Energy in) const { _pt=in; } //@} /** * Set and get methods for the SplittingFunction object */ //@{ /** * Access the SplittingFunction object responsible of the * eventual branching of this particle. */ tSplittingFnPtr splittingFn() const { return _sudakov-> splittingFn(); } //@} /** * Set and get methods for the SudakovFormFactor object */ /** * Access the SudakovFormFactor object responsible of the * eventual branching of this particle. */ tSudakovPtr SudakovFormFactor() const { return _sudakov; } /** * Set the SudakovFormFactor object responsible of the * eventual branching of this particle. */ void SudakovFormFactor(const tSudakovPtr sud) { _sudakov=sud; } //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ ShowerKinematics & operator=(const ShowerKinematics &) = delete; private: /** * The \f$\tilde{q}\f$ evolution variable. */ Energy _scale; /** * The energy fraction, \f$z\f$ */ double _z; /** * The azimuthal angle, \f$\phi\f$. */ double _phi; /** * The relative \f$p_T\f$ */ mutable Energy _pt; /** * The splitting function for the branching of the particle */ tSudakovPtr _sudakov; }; } #endif /* HERWIG_ShowerKinematics_H */ diff --git a/Shower/QTilde/QTildeShowerHandler.cc b/Shower/QTilde/QTildeShowerHandler.cc --- a/Shower/QTilde/QTildeShowerHandler.cc +++ b/Shower/QTilde/QTildeShowerHandler.cc @@ -1,3208 +1,3208 @@ // -*- C++ -*- // // QTildeShowerHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 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/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" using namespace Herwig; bool QTildeShowerHandler::_hardEmissionWarn = true; bool QTildeShowerHandler::_missingTruncWarn = true; QTildeShowerHandler::QTildeShowerHandler() : _maxtry(100), _meCorrMode(1), _reconOpt(0), _hardVetoReadOption(false), _iptrms(ZERO), _beta(0.), _gamma(ZERO), _iptmax(), _limitEmissions(0), _initialenhance(1.), _finalenhance(1.), _nReWeight(100), _reWeight(false), interaction_(ShowerInteraction::Both), _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 << _reconOpt << 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 >> _reconOpt >> 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{Marchesini:1983bm,Marchesini:1987cf,Gieseke:2003rz,Bahr:2008pv}.", "%\\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" ); 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::Both, false, false); static SwitchOption interfaceInteractionsQCD (interfaceInteractions, "QCD", "Only QCD radiation", ShowerInteraction::QCD); static SwitchOption interfaceInteractionsQED (interfaceInteractions, "QED", "Only QEd radiation", ShowerInteraction::QED); static SwitchOption interfaceInteractionsQCDandQED (interfaceInteractions, "QCDandQED", "Both QED and QCD radiation", ShowerInteraction::Both); static Switch interfaceReconstructionOption ("ReconstructionOption", "Treatment of the reconstruction of the transverse momentum of " "a branching from the evolution scale.", &QTildeShowerHandler::_reconOpt, 0, false, false); static SwitchOption interfaceReconstructionOptionCutOff (interfaceReconstructionOption, "CutOff", "Use the cut-off masses in the calculation", 0); static SwitchOption interfaceReconstructionOptionOffShell5 (interfaceReconstructionOption, "OffShell5", "Try and preserve q2 but if pt negative just zero it", 5); 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) { _intrinsic.clear(); 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,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(); if(_reconOpt>=1) - particle->showerKinematics()->updateParent(particle, children,fb.type); + particle->showerKinematics()->updateParent(particle, children,2,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(firstPDF().particle() == _beam) pdf = firstPDF().pdf(); if(secondPDF().particle() == _beam) 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,bb.type); + updateParent(newParent, theChildren,2,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,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); updateHistory(otherChild); if(theChildren[1]->spinInfo()) theChildren[1]->spinInfo()->develop(); // return the emitted if(particle->spinInfo()) particle->spinInfo()->develop(); 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, 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_noAOchildren().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,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,fb.type); + particle->showerKinematics()->updateParent(particle, children,2,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, branch->type()); + updateParent( newParent, theChildren,2, 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,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, bb.type); + updateParent( newParent, theChildren,2, 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, 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, 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 if(hard) { 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 { // 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_, switchRecon && ntry>maximumTries()/2) : 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(); }