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,211 +1,213 @@ // -*- C++ -*- // // FS_QTildeShowerKinematics1to2.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the 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; + if(ShowerHandler::currentHandlerIsSet() && + ! dynamic_ptr_cast(ShowerHandler::currentHandler())->correlations()) return; SpinPtr pspin(parent->spinInfo()); - if(!pspin || !dynamic_ptr_cast(ShowerHandler::currentHandler())->spinCorrelations() ) return; + if(!pspin || (ShowerHandler::currentHandlerIsSet() && + !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()) + if(!(mass > ZERO) && (!ShowerHandler::currentHandlerIsSet() ||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; 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/KinematicsReconstructor.cc b/Shower/QTilde/Kinematics/KinematicsReconstructor.cc --- a/Shower/QTilde/Kinematics/KinematicsReconstructor.cc +++ b/Shower/QTilde/Kinematics/KinematicsReconstructor.cc @@ -1,2832 +1,2826 @@ // -*- C++ -*- // // KinematicsReconstructor.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the KinematicsReconstructor class. // #include "KinematicsReconstructor.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/EventRecord/Event.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Deleted.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/RefVector.h" #include "Herwig/Shower/QTilde/Base/PartnerFinder.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Shower/QTilde/SplittingFunctions/SplittingFunction.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/EventRecord/ColourLine.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Shower/QTilde/QTildeShowerHandler.h" #include #include "KinematicsReconstructor.tcc" using namespace Herwig; DescribeClass describeKinematicsReconstructor("Herwig::KinematicsReconstructor", "HwShower.so"); namespace { /** * Struct to order the jets in off-shellness */ struct JetOrdering { bool operator() (const JetKinStruct & j1, const JetKinStruct & j2) { Energy diff1 = j1.q.m()-j1.p.m(); Energy diff2 = j2.q.m()-j2.p.m(); if(diff1!=diff2) { return diff1>diff2; } else if( j1.q.e() != j2.q.e() ) return j1.q.e()>j2.q.e(); else return j1.parent->uniqueId>j2.parent->uniqueId; } }; } void KinematicsReconstructor::persistentOutput(PersistentOStream & os) const { os << _reconopt << _initialBoost << ounit(_minQ,GeV) << _noRescale << _noRescaleVector << _initialStateReconOption << _finalFinalWeight; } void KinematicsReconstructor::persistentInput(PersistentIStream & is, int) { is >> _reconopt >> _initialBoost >> iunit(_minQ,GeV) >> _noRescale >> _noRescaleVector >> _initialStateReconOption >> _finalFinalWeight; } void KinematicsReconstructor::Init() { static ClassDocumentation documentation ( "This class is responsible for the kinematics reconstruction of the showering,", " including the kinematics reshuffling necessary to compensate for the recoil" "of the emissions." ); static Switch interfaceReconstructionOption ("ReconstructionOption", "Option for the kinematics reconstruction", &KinematicsReconstructor::_reconopt, 0, false, false); static SwitchOption interfaceReconstructionOptionGeneral (interfaceReconstructionOption, "General", "Use the general solution which ignores the colour structure for all processes", 0); static SwitchOption interfaceReconstructionOptionColour (interfaceReconstructionOption, "Colour", "Use the colour structure of the process to determine the reconstruction procedure.", 1); static SwitchOption interfaceReconstructionOptionColour2 (interfaceReconstructionOption, "Colour2", "Make the most use possible of the colour structure of the process to determine the reconstruction procedure. " "Start with FF, then IF then II colour connections", 2); static SwitchOption interfaceReconstructionOptionColour3 (interfaceReconstructionOption, "Colour3", "Make the most use possible of the colour structure of the process to determine the reconstruction procedure. " "Do the colour connections in order of the pT's emitted in the shower starting with the hardest." " The colour partner is fully reconstructed at the same time.", 3); static SwitchOption interfaceReconstructionOptionColour4 (interfaceReconstructionOption, "Colour4", "Make the most use possible of the colour structure of the process to determine the reconstruction procedure. " "Do the colour connections in order of the pT's emitted in the shower starting with the hardest, while leaving" " the colour partner on mass-shell", 4); static Parameter interfaceMinimumQ2 ("MinimumQ2", "The minimum Q2 for the reconstruction of initial-final systems", &KinematicsReconstructor::_minQ, GeV, 0.001*GeV, 1e-6*GeV, 10.0*GeV, false, false, Interface::limited); static RefVector interfaceNoRescale ("NoRescale", "Particles which shouldn't be rescaled to be on shell by the shower", &KinematicsReconstructor::_noRescaleVector, -1, false, false, true, false, false); static Switch interfaceInitialInitialBoostOption ("InitialInitialBoostOption", "Option for how the boost from the system before ISR to that after ISR is applied.", &KinematicsReconstructor::_initialBoost, 0, false, false); static SwitchOption interfaceInitialInitialBoostOptionOneBoost (interfaceInitialInitialBoostOption, "OneBoost", "Apply one boost from old CMS to new CMS", 0); static SwitchOption interfaceInitialInitialBoostOptionLongTransBoost (interfaceInitialInitialBoostOption, "LongTransBoost", "First apply a longitudinal and then a transverse boost", 1); static Deleted delFinalStateReconOption ("FinalStateReconOption", "The old default (0) is now the only choice"); static Switch interfaceInitialStateReconOption ("InitialStateReconOption", "Option for the reconstruction of initial state radiation", &KinematicsReconstructor::_initialStateReconOption, 0, false, false); static SwitchOption interfaceInitialStateReconOptionRapidity (interfaceInitialStateReconOption, "Rapidity", "Preserve shat and rapidity", 0); static SwitchOption interfaceInitialStateReconOptionLongitudinal (interfaceInitialStateReconOption, "Longitudinal", "Preserve longitudinal momentum", 1); static SwitchOption interfaceInitialStateReconOptionSofterFraction (interfaceInitialStateReconOption, "SofterFraction", "Preserve the momentum fraction of the parton which has emitted softer.", 2); static Switch interfaceFinalFinalWeight ("FinalFinalWeight", "Apply kinematic rejection weight for final-states", &KinematicsReconstructor::_finalFinalWeight, false, false, false); static SwitchOption interfaceFinalFinalWeightNo (interfaceFinalFinalWeight, "No", "Don't apply the weight", false); static SwitchOption interfaceFinalFinalWeightYes (interfaceFinalFinalWeight, "Yes", "Apply the weight", true); } void KinematicsReconstructor::doinit() { Interfaced::doinit(); _noRescale = set(_noRescaleVector.begin(),_noRescaleVector.end()); } bool KinematicsReconstructor:: -reconstructTimeLikeJet(const tShowerParticlePtr particleJetParent) const { +reconstructTimeLikeJet(const tShowerParticlePtr particleJetParent, + const tShowerParticlePtr progenitor) const { assert(particleJetParent); bool emitted=true; // if this is not a fixed point in the reconstruction if( !particleJetParent->children().empty() ) { // if not a reconstruction fixpoint, dig deeper for all children: for ( ParticleVector::const_iterator cit = particleJetParent->children().begin(); cit != particleJetParent->children().end(); ++cit ) - reconstructTimeLikeJet(dynamic_ptr_cast(*cit)); + reconstructTimeLikeJet(dynamic_ptr_cast(*cit),progenitor); } // it is a reconstruction fixpoint, ie kinematical data has to be available else { // check if the parent was part of the shower ShowerParticlePtr jetGrandParent; if(!particleJetParent->parents().empty()) jetGrandParent= dynamic_ptr_cast (particleJetParent->parents()[0]); // update if so if (jetGrandParent) { if (jetGrandParent->showerKinematics()) { - if(particleJetParent->id()==_progenitor->id()&& - !_progenitor->data().stable()&&abs(_progenitor->data().id())!=ParticleID::tauminus) { + if(particleJetParent->id()==progenitor->id()&& + !progenitor->data().stable()&&abs(progenitor->data().id())!=ParticleID::tauminus) { jetGrandParent->showerKinematics()->reconstructLast(particleJetParent, - _progenitor->mass()); + progenitor->mass()); } else { jetGrandParent->showerKinematics()->reconstructLast(particleJetParent); } } } // otherwise else { - Energy dm = ShowerHandler::currentHandler()->retConstituentMasses()? + Energy dm = !ShowerHandler::currentHandlerIsSet() || + ShowerHandler::currentHandler()->retConstituentMasses()? particleJetParent->data().constituentMass(): particleJetParent->data().mass(); if (abs(dm-particleJetParent->momentum().m())>0.001*MeV &&(particleJetParent->dataPtr()->stable() || abs(particleJetParent->id())==ParticleID::tauminus) &&particleJetParent->id()!=ParticleID::gamma &&_noRescale.find(particleJetParent->dataPtr())==_noRescale.end()) { Lorentz5Momentum dum = particleJetParent->momentum(); dum.setMass(dm); dum.rescaleEnergy(); if(abs(particleJetParent->id())==15&&particleJetParent->spinInfo()) { if(particleJetParent->spinInfo()->isNear(particleJetParent->momentum())) { particleJetParent->spinInfo()->SpinInfo::transform(dum,LorentzRotation()); } } particleJetParent->set5Momentum(dum); } else { emitted=false; } } } // recursion has reached an endpoint once, ie we can reconstruct the // kinematics from the children. if( !particleJetParent->children().empty() ) particleJetParent->showerKinematics() ->reconstructParent( particleJetParent, particleJetParent->children() ); return emitted; } bool KinematicsReconstructor:: reconstructHardJets(ShowerTreePtr hard, const map > & intrinsic, ShowerInteraction type, bool switchRecon) const { _currentTree = hard; _intrinsic=intrinsic; // extract the particles from the ShowerTree vector ShowerHardJets=hard->extractProgenitors(); for(unsigned int ix=0;ixprogenitor()] = vector(); } for(map >::const_iterator tit = _currentTree->treelinks().begin(); tit != _currentTree->treelinks().end();++tit) { _treeBoosts[tit->first] = vector(); } try { // old recon method, using new member functions if(_reconopt == 0 || switchRecon ) { reconstructGeneralSystem(ShowerHardJets); } // reconstruction based on coloured systems else if( _reconopt == 1) { reconstructColourSinglets(ShowerHardJets,type); } // reconstruction of FF, then IF, then II else if( _reconopt == 2) { reconstructFinalFirst(ShowerHardJets); } // reconstruction based on coloured systems else if( _reconopt == 3 || _reconopt == 4) { reconstructColourPartner(ShowerHardJets); } else assert(false); } catch(KinematicsReconstructionVeto) { - _progenitor=tShowerParticlePtr(); _intrinsic.clear(); for(map >::const_iterator bit=_boosts.begin();bit!=_boosts.end();++bit) { for(vector::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) { LorentzRotation rot = rit->inverse(); bit->first->transform(rot); } } _boosts.clear(); for(map >::const_iterator bit=_treeBoosts.begin();bit!=_treeBoosts.end();++bit) { for(vector::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) { LorentzRotation rot = rit->inverse(); bit->first->transform(rot,false); } } _currentTree = tShowerTreePtr(); _treeBoosts.clear(); return false; } catch (Exception & ex) { - _progenitor=tShowerParticlePtr(); _intrinsic.clear(); _currentTree = tShowerTreePtr(); _boosts.clear(); _treeBoosts.clear(); throw ex; } - _progenitor=tShowerParticlePtr(); _intrinsic.clear(); // ensure x<1 for(map::const_iterator cit=hard->incomingLines().begin();cit!=hard->incomingLines().end();++cit) { tPPtr parent = cit->first->progenitor(); while (!parent->parents().empty()) { parent = parent->parents()[0]; } tPPtr hadron; if ( cit->first->original()->parents().empty() ) { hadron = cit->first->original(); } else { hadron = cit->first->original()->parents()[0]; } if( ! (hadron->id() == parent->id() && hadron->children().size() <= 1) && parent->momentum().rho() > hadron->momentum().rho()) { - _progenitor=tShowerParticlePtr(); _intrinsic.clear(); for(map >::const_iterator bit=_boosts.begin();bit!=_boosts.end();++bit) { for(vector::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) { LorentzRotation rot = rit->inverse(); bit->first->transform(rot); } } _boosts.clear(); for(map >::const_iterator bit=_treeBoosts.begin();bit!=_treeBoosts.end();++bit) { for(vector::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) { LorentzRotation rot = rit->inverse(); bit->first->transform(rot,false); } } _currentTree = tShowerTreePtr(); _treeBoosts.clear(); return false; } } _boosts.clear(); _treeBoosts.clear(); _currentTree = tShowerTreePtr(); return true; } double KinematicsReconstructor::solveKfactor(const Energy & root_s, const JetKinVect & jets) const { Energy2 s = sqr(root_s); // must be at least two jets if ( jets.size() < 2) throw KinematicsReconstructionVeto(); // sum of jet masses must be less than roots if(momConsEq( 0.0, root_s, jets )>ZERO) throw KinematicsReconstructionVeto(); // if two jets simple solution if ( jets.size() == 2 ) { static const Energy2 eps = 1.0e-4 * MeV2; if ( sqr(jets[0].p.x()+jets[1].p.x()) < eps && sqr(jets[0].p.y()+jets[1].p.y()) < eps && sqr(jets[0].p.z()+jets[1].p.z()) < eps ) { Energy test = (jets[0].p+jets[1].p).vect().mag(); if(test > 1.0e-4 * MeV) throw KinematicsReconstructionVeto(); if ( jets[0].p.vect().mag2() < eps ) throw KinematicsReconstructionVeto(); Energy2 m1sq(jets[0].q.m2()),m2sq(jets[1].q.m2()); return sqrt( ( sqr(s - m1sq - m2sq) - 4.*m1sq*m2sq ) /(4.*s*jets[0].p.vect().mag2()) ); } else throw KinematicsReconstructionVeto(); } // i.e. jets.size() > 2, numerically // check convergence, if it's a problem maybe use Newton iteration? else { double k1 = 0.,k2 = 1.,k = 0.; if ( momConsEq( k1, root_s, jets ) < ZERO ) { while ( momConsEq( k2, root_s, jets ) < ZERO ) { k1 = k2; k2 *= 2; } while ( fabs( (k1 - k2)/(k1 + k2) ) > 1.e-10 ) { if( momConsEq( k2, root_s, jets ) == ZERO ) { return k2; } else { k = (k1+k2)/2.; if ( momConsEq( k, root_s, jets ) > ZERO ) { k2 = k; } else { k1 = k; } } } return k1; } else throw KinematicsReconstructionVeto(); } throw KinematicsReconstructionVeto(); } bool KinematicsReconstructor:: reconstructSpaceLikeJet( const tShowerParticlePtr p) const { bool emitted = true; tShowerParticlePtr child; tShowerParticlePtr parent; if(!p->parents().empty()) parent = dynamic_ptr_cast(p->parents()[0]); if(parent) { emitted=true; reconstructSpaceLikeJet(parent); } // if branching reconstruct time-like child if(p->children().size()==2) child = dynamic_ptr_cast(p->children()[1]); if(p->perturbative()==0 && child) { dynamic_ptr_cast(p->children()[0])-> showerKinematics()->reconstructParent(p,p->children()); if(!child->children().empty()) { - _progenitor=child; - reconstructTimeLikeJet(child); + reconstructTimeLikeJet(child,child); // calculate the momentum of the particle Lorentz5Momentum pnew=p->momentum()-child->momentum(); pnew.rescaleMass(); p->children()[0]->set5Momentum(pnew); } } return emitted; } Boost KinematicsReconstructor:: solveBoostBeta( const double k, const Lorentz5Momentum & newq, const Lorentz5Momentum & oldp ) { // try something different, purely numerical first: // a) boost to rest frame of newq, b) boost with kp/E Energy q = newq.vect().mag(); Energy2 qs = sqr(q); Energy2 Q2 = newq.m2(); Energy kp = k*(oldp.vect().mag()); Energy2 kps = sqr(kp); // usually we take the minus sign, since this boost will be smaller. // we only require |k \vec p| = |\vec q'| which leaves the sign of // the boost open but the 'minus' solution gives a smaller boost // parameter, i.e. the result should be closest to the previous // result. this is to be changed if we would get many momentum // conservation violations at the end of the shower from a hard // process. double betam = (q*sqrt(qs + Q2) - kp*sqrt(kps + Q2))/(kps + qs + Q2); // move directly to 'return' Boost beta = -betam*(k/kp)*oldp.vect(); // note that (k/kp)*oldp.vect() = oldp.vect()/oldp.vect().mag() but cheaper. // leave this out if it's running properly! if ( betam >= 0 ) return beta; else return Boost(0., 0., 0.); } bool KinematicsReconstructor:: reconstructDecayJets(ShowerTreePtr decay, ShowerInteraction) const { _currentTree = decay; // extract the particles from the ShowerTree vector ShowerHardJets=decay->extractProgenitors(); for(unsigned int ix=0;ixprogenitor()] = vector(); } for(map >::const_iterator tit = _currentTree->treelinks().begin(); tit != _currentTree->treelinks().end();++tit) { _treeBoosts[tit->first] = vector(); } try { bool radiated[2]={false,false}; // find the decaying particle and check if particles radiated ShowerProgenitorPtr initial; for(unsigned int ix=0;ixprogenitor()->isFinalState()) { radiated[1] |=ShowerHardJets[ix]->hasEmitted(); } else { initial=ShowerHardJets[ix]; radiated[0]|=ShowerHardJets[ix]->hasEmitted(); } } // find boost to the rest frame if needed Boost boosttorest=-initial->progenitor()->momentum().boostVector(); double gammarest = initial->progenitor()->momentum().e()/ initial->progenitor()->momentum().mass(); // check if need to boost to rest frame bool gottaBoost = (boosttorest.mag() > 1e-12); // if initial state radiation reconstruct the jet and set up the basis vectors Lorentz5Momentum pjet; Lorentz5Momentum nvect; // find the partner ShowerParticlePtr partner = initial->progenitor()->partner(); Lorentz5Momentum ppartner[2]; if(partner) ppartner[0]=partner->momentum(); // get the n reference vector if(partner) { if(initial->progenitor()->showerKinematics()) { nvect = initial->progenitor()->showerBasis()->getBasis()[1]; } else { Lorentz5Momentum ppartner=initial->progenitor()->partner()->momentum(); if(gottaBoost) ppartner.boost(boosttorest,gammarest); nvect = Lorentz5Momentum( ZERO,0.5*initial->progenitor()->mass()* ppartner.vect().unit()); nvect.boost(-boosttorest,gammarest); } } // if ISR if(radiated[0]) { // reconstruct the decay jet reconstructDecayJet(initial->progenitor()); // momentum of decaying particle after ISR pjet=initial->progenitor()->momentum() -decay->incomingLines().begin()->second->momentum(); pjet.rescaleMass(); } // boost initial state jet and basis vector if needed if(gottaBoost) { pjet.boost(boosttorest,gammarest); nvect.boost(boosttorest,gammarest); ppartner[0].boost(boosttorest,gammarest); } // loop over the final-state particles and do the reconstruction JetKinVect possiblepartners; JetKinVect jetKinematics; bool atLeastOnce = radiated[0]; LorentzRotation restboost(boosttorest,gammarest); Energy inmass(ZERO); for(unsigned int ix=0;ixprogenitor()->isFinalState()) { inmass=ShowerHardJets[ix]->progenitor()->mass(); continue; } // do the reconstruction JetKinStruct tempJetKin; tempJetKin.parent = ShowerHardJets[ix]->progenitor(); if(ShowerHardJets.size()==2) { Lorentz5Momentum dum=ShowerHardJets[ix]->progenitor()->momentum(); dum.setMass(inmass); dum.rescaleRho(); tempJetKin.parent->set5Momentum(dum); } tempJetKin.p = ShowerHardJets[ix]->progenitor()->momentum(); if(gottaBoost) tempJetKin.p.boost(boosttorest,gammarest); - _progenitor=tempJetKin.parent; if(ShowerHardJets[ix]->reconstructed()==ShowerProgenitor::notReconstructed) { - atLeastOnce |= reconstructTimeLikeJet(tempJetKin.parent); + atLeastOnce |= reconstructTimeLikeJet(tempJetKin.parent,tempJetKin.parent); ShowerHardJets[ix]->reconstructed(ShowerProgenitor::done); } if(gottaBoost) deepTransform(tempJetKin.parent,restboost); tempJetKin.q = ShowerHardJets[ix]->progenitor()->momentum(); jetKinematics.push_back(tempJetKin); } if(partner) ppartner[1]=partner->momentum(); // calculate the rescaling parameters double k1,k2; Lorentz5Momentum qt; if(!solveDecayKFactor(initial->progenitor()->mass(),nvect,pjet, jetKinematics,partner,ppartner,k1,k2,qt)) { for(map >::const_iterator bit=_boosts.begin();bit!=_boosts.end();++bit) { for(vector::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) { LorentzRotation rot = rit->inverse(); bit->first->transform(rot); } } _boosts.clear(); for(map >::const_iterator bit=_treeBoosts.begin();bit!=_treeBoosts.end();++bit) { for(vector::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) { LorentzRotation rot = rit->inverse(); bit->first->transform(rot,false); } } _treeBoosts.clear(); _currentTree = tShowerTreePtr(); return false; } // apply boosts and rescalings to final-state jets for(JetKinVect::iterator it = jetKinematics.begin(); it != jetKinematics.end(); ++it) { LorentzRotation Trafo = LorentzRotation(); if(it->parent!=partner) { // boost for rescaling if(atLeastOnce) { map >::const_iterator tit; for(tit = _currentTree->treelinks().begin(); tit != _currentTree->treelinks().end();++tit) { if(tit->second.first && tit->second.second==it->parent) break; } if(it->parent->children().empty()&&!it->parent->spinInfo() && tit==_currentTree->treelinks().end()) { Lorentz5Momentum pnew(k2*it->p.vect(), sqrt(sqr(k2*it->p.vect().mag())+it->q.mass2()), it->q.mass()); it->parent->set5Momentum(pnew); } else { // rescaling boost can't ever work in this case if(k2<0. && it->q.mass()==ZERO) throw KinematicsReconstructionVeto(); Trafo = solveBoost(k2, it->q, it->p); } } if(gottaBoost) Trafo.boost(-boosttorest,gammarest); if(atLeastOnce || gottaBoost) deepTransform(it->parent,Trafo); } else { Lorentz5Momentum pnew=ppartner[0]; pnew *=k1; pnew-=qt; pnew.setMass(ppartner[1].mass()); pnew.rescaleEnergy(); LorentzRotation Trafo=solveBoost(1.,ppartner[1],pnew); if(gottaBoost) Trafo.boost(-boosttorest,gammarest); deepTransform(partner,Trafo); } } } catch(KinematicsReconstructionVeto) { for(map >::const_iterator bit=_boosts.begin();bit!=_boosts.end();++bit) { for(vector::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) { LorentzRotation rot = rit->inverse(); bit->first->transform(rot); } } _boosts.clear(); for(map >::const_iterator bit=_treeBoosts.begin();bit!=_treeBoosts.end();++bit) { for(vector::const_reverse_iterator rit=bit->second.rbegin();rit!=bit->second.rend();++rit) { LorentzRotation rot = rit->inverse(); bit->first->transform(rot,false); } } _treeBoosts.clear(); _currentTree = tShowerTreePtr(); return false; } catch (Exception & ex) { _currentTree = tShowerTreePtr(); _boosts.clear(); _treeBoosts.clear(); throw ex; } _boosts.clear(); _treeBoosts.clear(); _currentTree = tShowerTreePtr(); return true; } bool KinematicsReconstructor:: reconstructDecayJet( const tShowerParticlePtr p) const { if(p->children().empty()) return false; tShowerParticlePtr child; // if branching reconstruct time-like child child = dynamic_ptr_cast(p->children()[1]); if(child) { - _progenitor=child; - reconstructTimeLikeJet(child); + reconstructTimeLikeJet(child,child); // calculate the momentum of the particle Lorentz5Momentum pnew=p->momentum()-child->momentum(); pnew.rescaleMass(); p->children()[0]->set5Momentum(pnew); child=dynamic_ptr_cast(p->children()[0]); reconstructDecayJet(child); return true; } return false; } bool KinematicsReconstructor:: solveDecayKFactor(Energy mb, const Lorentz5Momentum & n, const Lorentz5Momentum & pjet, const JetKinVect & jetKinematics, ShowerParticlePtr partner, Lorentz5Momentum ppartner[2], double & k1, double & k2, Lorentz5Momentum & qt) const { Energy2 pjn = partner ? pjet.vect()*n.vect() : ZERO; Energy2 pcn = partner ? ppartner[0].vect()*n.vect() : 1.*MeV2; Energy2 nmag = n.vect().mag2(); Lorentz5Momentum pn = partner ? (pjn/nmag)*n : Lorentz5Momentum(); qt=pjet-pn; qt.setE(ZERO); Energy2 pt2=qt.vect().mag2(); Energy Ejet = pjet.e(); // magnitudes of the momenta for fast access vector pmag; Energy total(Ejet); for(unsigned int ix=0;ixmb) return false; Energy2 pcmag=ppartner[0].vect().mag2(); // used newton-raphson to get the rescaling static const Energy eps=1e-8*GeV; long double d1(1.),d2(1.); Energy roots, ea, ec, ds; unsigned int ix=0; do { ++ix; d2 = d1 + pjn/pcn; roots = Ejet; ds = ZERO; for(unsigned int iy=0;iyeps && ix<100); k1=d1; k2=d2; // return true if N-R succeed, otherwise false return ix<100; } bool KinematicsReconstructor:: deconstructDecayJets(HardTreePtr decay,ShowerInteraction) const { // extract the momenta of the particles vector pin; vector pout; // on-shell masses of the decay products vector mon; Energy mbar(-GeV); // the hard branchings of the particles set::iterator cit; set branchings=decay->branchings(); // properties of the incoming particle bool ISR = false; HardBranchingPtr initial; Lorentz5Momentum qisr; // find the incoming particle, both before and after // any ISR for(cit=branchings.begin();cit!=branchings.end();++cit){ if((*cit)->status()==HardBranching::Incoming|| (*cit)->status()==HardBranching::Decay) { // search back up isr if needed HardBranchingPtr branch = *cit; while(branch->parent()) branch=branch->parent(); initial=branch; // momentum or original parent pin.push_back(branch->branchingParticle()->momentum()); // ISR? ISR = !branch->branchingParticle()->children().empty(); // ISR momentum qisr = pin.back()-(**cit).branchingParticle()->momentum(); qisr.rescaleMass(); } } assert(pin.size()==1); // compute boost to rest frame Boost boostv=-pin[0].boostVector(); // partner for ISR ShowerParticlePtr partner; Lorentz5Momentum ppartner; if(initial->branchingParticle()->partner()) { partner=initial->branchingParticle()->partner(); ppartner=partner->momentum(); } // momentum of the decay products for(cit=branchings.begin();cit!=branchings.end();++cit) { if((*cit)->status()!=HardBranching::Outgoing) continue; // find the mass of the particle // including special treatment for off-shell resonances // to preserve off-shell mass Energy mass; if(!(**cit).branchingParticle()->dataPtr()->stable()) { HardBranchingPtr branch=*cit; while(!branch->children().empty()) { for(unsigned int ix=0;ixchildren().size();++ix) { if(branch->children()[ix]->branchingParticle()->id()== (**cit).branchingParticle()->id()) { branch = branch->children()[ix]; continue; } } }; mass = branch->branchingParticle()->mass(); } else { mass = (**cit).branchingParticle()->dataPtr()->mass(); } // if not evolution partner of decaying particle if((*cit)->branchingParticle()!=partner) { pout.push_back((*cit)->branchingParticle()->momentum()); mon.push_back(mass); } // evolution partner of decaying particle else { mbar = mass; } } // boost all the momenta to the rest frame of the decaying particle for(unsigned int ix=0;ixbranchingParticle()->partner()) { ppartner.boost(boostv); qisr.boost(boostv); } // compute the rescaling factors double k1,k2; if(!ISR) { if(partner) { pout.push_back(ppartner); mon.push_back(mbar); } k1=k2=inverseRescalingFactor(pout,mon,pin[0].mass()); if(partner) { pout.pop_back(); mon.pop_back(); } } else { if(!inverseDecayRescalingFactor(pout,mon,pin[0].mass(), ppartner,mbar,k1,k2)) return false; } // now calculate the p reference vectors unsigned int ifinal=0; for(cit=branchings.begin();cit!=branchings.end();++cit) { if((**cit).status()!=HardBranching::Outgoing) continue; // for partners other than colour partner of decaying particle if((*cit)->branchingParticle()!=partner) { Lorentz5Momentum pvect = (*cit)->branchingParticle()->momentum(); pvect.boost(boostv); pvect /= k1; pvect.setMass(mon[ifinal]); ++ifinal; pvect.rescaleEnergy(); pvect.boost(-boostv); (*cit)->pVector(pvect); (*cit)->showerMomentum(pvect); } // for colour partner of decaying particle else { Lorentz5Momentum pvect = (*cit)->branchingParticle()->momentum(); pvect.boost(boostv); Lorentz5Momentum qtotal; for(unsigned int ix=0;ixpVector(pvect); (*cit)->showerMomentum(pvect); } } // For initial-state if needed if(initial) { tShowerParticlePtr newPartner=initial->branchingParticle()->partner(); if(newPartner) { tHardBranchingPtr branch; for( set::iterator clt = branchings.begin(); clt != branchings.end(); ++clt ) { if((**clt).branchingParticle()==newPartner) { initial->colourPartner(*clt); branch=*clt; break; } } Lorentz5Momentum pvect = initial->branchingParticle()->momentum(); initial->pVector(pvect); Lorentz5Momentum ptemp = branch->pVector(); ptemp.boost(boostv); Lorentz5Momentum nvect = Lorentz5Momentum( ZERO, 0.5*initial->branchingParticle()->mass()* ptemp.vect().unit()); nvect.boost(-boostv); initial->nVector(nvect); } } // calculate the reference vectors, then for outgoing particles for(cit=branchings.begin();cit!=branchings.end();++cit){ if((**cit).status()!=HardBranching::Outgoing) continue; // find the partner branchings tShowerParticlePtr newPartner=(*cit)->branchingParticle()->partner(); if(!newPartner) continue; tHardBranchingPtr branch; for( set::iterator clt = branchings.begin(); clt != branchings.end(); ++clt ) { if(cit==clt) continue; if((**clt).branchingParticle()==newPartner) { (**cit).colourPartner(*clt); branch=*clt; break; } } if((**decay->incoming().begin()).branchingParticle()==newPartner) { (**cit).colourPartner(*decay->incoming().begin()); branch = *decay->incoming().begin(); } // final-state colour partner if(branch->status()==HardBranching::Outgoing) { Boost boost=((*cit)->pVector()+branch->pVector()).findBoostToCM(); Lorentz5Momentum pcm = branch->pVector(); pcm.boost(boost); Lorentz5Momentum nvect = Lorentz5Momentum(ZERO,pcm.vect()); nvect.boost( -boost); (*cit)->nVector(nvect); } // initial-state colour partner else { Boost boost=branch->pVector().findBoostToCM(); Lorentz5Momentum pcm = (*cit)->pVector(); pcm.boost(boost); Lorentz5Momentum nvect = Lorentz5Momentum( ZERO, -pcm.vect()); nvect.boost( -boost); (*cit)->nVector(nvect); } } // now compute the new momenta // and calculate the shower variables for(cit=branchings.begin();cit!=branchings.end();++cit) { if((**cit).status()!=HardBranching::Outgoing) continue; LorentzRotation B=LorentzRotation(-boostv); LorentzRotation A=LorentzRotation(boostv),R; if((*cit)->branchingParticle()==partner) { Lorentz5Momentum qnew; Energy2 dot=(*cit)->pVector()*(*cit)->nVector(); double beta = 0.5*((*cit)->branchingParticle()->momentum().m2() -sqr((*cit)->pVector().mass()))/dot; qnew=(*cit)->pVector()+beta*(*cit)->nVector(); qnew.rescaleMass(); // compute the boost R=B*solveBoost(A*qnew,A*(*cit)->branchingParticle()->momentum())*A; } else { Lorentz5Momentum qnew; if((*cit)->branchingParticle()->partner()) { Energy2 dot=(*cit)->pVector()*(*cit)->nVector(); double beta = 0.5*((*cit)->branchingParticle()->momentum().m2() -sqr((*cit)->pVector().mass()))/dot; qnew=(*cit)->pVector()+beta*(*cit)->nVector(); qnew.rescaleMass(); } else { qnew = (*cit)->pVector(); } // compute the boost R=B*solveBoost(A*qnew,A*(*cit)->branchingParticle()->momentum())*A; } // reconstruct the momenta (*cit)->setMomenta(R,1.0,Lorentz5Momentum()); } if(initial) { initial->setMomenta(LorentzRotation(),1.0,Lorentz5Momentum()); } return true; } double KinematicsReconstructor:: inverseRescalingFactor(vector pout, vector mon, Energy roots) const { double lambda=1.; if(pout.size()==2) { double mu_q1(pout[0].m()/roots), mu_q2(pout[1].m()/roots); double mu_p1(mon[0]/roots) , mu_p2(mon[1]/roots); lambda = ((1.+mu_q1+mu_q2)*(1.-mu_q1-mu_q2)*(mu_q1-1.-mu_q2)*(mu_q2-1.-mu_q1))/ ((1.+mu_p1+mu_p2)*(1.-mu_p1-mu_p2)*(mu_p1-1.-mu_p2)*(mu_p2-1.-mu_p1)); if(lambda<0.) throw Exception() << "Rescaling factor is imaginary in KinematicsReconstructor::" << "inverseRescalingFactor lambda^2= " << lambda << Exception::eventerror; lambda = sqrt(lambda); } else { unsigned int ntry=0; // compute magnitudes once for speed vector pmag; for(unsigned int ix=0;ix root(pout.size()); do { // compute new energies Energy sum(ZERO); for(unsigned int ix=0;ix::const_iterator it=tree->branchings().begin(); it!=tree->branchings().end();++it) { if((**it).status()==HardBranching::Incoming) in .jets.push_back(*it); else out.jets.push_back(*it); } LorentzRotation toRest,fromRest; bool applyBoost(false); // do the initial-state reconstruction deconstructInitialInitialSystem(applyBoost,toRest,fromRest, tree,in.jets,type); // do the final-state reconstruction deconstructFinalStateSystem(toRest,fromRest,tree, out.jets,type); // only at this point that we can be sure all the reference vectors // are correct for(set::const_iterator it=tree->branchings().begin(); it!=tree->branchings().end();++it) { if((**it).status()==HardBranching::Incoming) continue; if((**it).branchingParticle()->coloured()) (**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false); } for(set::const_iterator it=tree->incoming().begin(); it!=tree->incoming().end();++it) { (**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false); } return true; } bool KinematicsReconstructor::deconstructHardJets(HardTreePtr tree, ShowerInteraction type) const { // inverse of old recon method if(_reconopt == 0) { return deconstructGeneralSystem(tree,type); } else if(_reconopt == 1) { return deconstructColourSinglets(tree,type); } else if(_reconopt == 2) { throw Exception() << "Inverse reconstruction is not currently supported for ReconstructionOption Colour2 " << "in KinematicsReconstructor::deconstructHardJets(). Please use one of the other options\n" << Exception::runerror; } else if(_reconopt == 3 || _reconopt == 4 ) { return deconstructColourPartner(tree,type); } else { assert(false); return false; } } bool KinematicsReconstructor:: deconstructColourSinglets(HardTreePtr tree, ShowerInteraction type) const { // identify the colour singlet systems unsigned int nnun(0),nnii(0),nnif(0),nnf(0),nni(0); vector systems(identifySystems(tree->branchings(),nnun,nnii,nnif,nnf,nni)); // now decide what to do LorentzRotation toRest,fromRest; bool applyBoost(false); bool general(false); // initial-initial connection and final-state colour singlet systems // Drell-Yan type if(nnun==0&&nnii==1&&nnif==0&&nnf>0&&nni==0) { // reconstruct initial-initial system for(unsigned int ix=0;ix0&&nni==1)|| (nnif==2&& nni==0))) { for(unsigned int ix=0;ix0&&nni==2) { // only FS needed // but need to boost to rest frame if QED ISR Lorentz5Momentum ptotal; for(unsigned int ix=0;ixbranchingParticle()->momentum(); } toRest = LorentzRotation(ptotal.findBoostToCM()); fromRest = toRest; fromRest.invert(); if(type!=ShowerInteraction::QCD) { combineFinalState(systems); general=false; } } // general type else { general = true; } // final-state systems except for general recon if(!general) { for(unsigned int ix=0;ix::const_iterator it=tree->branchings().begin(); it!=tree->branchings().end();++it) { if((**it).status()==HardBranching::Incoming) continue; if((**it).branchingParticle()->coloured()) (**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false); } for(set::const_iterator it=tree->incoming().begin(); it!=tree->incoming().end();++it) { (**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false); } return true; } else { return deconstructGeneralSystem(tree,type); } return true; } bool KinematicsReconstructor:: deconstructColourPartner(HardTreePtr tree, ShowerInteraction type) const { Lorentz5Momentum ptotal; HardBranchingPtr emitter; ColourSingletShower incomingShower,outgoingShower; for(set::const_iterator it=tree->branchings().begin(); it!=tree->branchings().end();++it) { if((**it).status()==HardBranching::Incoming) { incomingShower.jets.push_back(*it); ptotal += (*it)->branchingParticle()->momentum(); // check for emitting particle if((**it).parent() ) { if(!emitter) emitter = *it; else throw Exception() << "Only one emitting particle allowed in " << "KinematicsReconstructor::deconstructColourPartner()" << Exception::runerror; } } else if ((**it).status()==HardBranching::Outgoing) { outgoingShower.jets.push_back(*it); // check for emitting particle if(!(**it).children().empty() ) { if(!emitter) emitter = *it; else throw Exception() << "Only one emitting particle allowed in " << "KinematicsReconstructor::deconstructColourPartner()" << Exception::runerror; } } } assert(emitter); assert(emitter->colourPartner()); ColourSingletShower system; system.jets.push_back(emitter); system.jets.push_back(emitter->colourPartner()); LorentzRotation toRest,fromRest; bool applyBoost(false); // identify the colour singlet system if(emitter->status() == HardBranching::Outgoing && emitter->colourPartner()->status() == HardBranching::Outgoing ) { system.type=F; // need to boost to rest frame if QED ISR if( !incomingShower.jets[0]->branchingParticle()->coloured() && !incomingShower.jets[1]->branchingParticle()->coloured() ) { Boost boost = ptotal.findBoostToCM(); toRest = LorentzRotation( boost); fromRest = LorentzRotation(-boost); } else findInitialBoost(ptotal,ptotal,toRest,fromRest); deconstructFinalStateSystem(toRest,fromRest,tree, system.jets,type); } else if (emitter->status() == HardBranching::Incoming && emitter->colourPartner()->status() == HardBranching::Incoming) { system.type=II; deconstructInitialInitialSystem(applyBoost,toRest,fromRest,tree,system.jets,type); // make sure the recoil gets applied deconstructFinalStateSystem(toRest,fromRest,tree, outgoingShower.jets,type); } else if ((emitter->status() == HardBranching::Outgoing && emitter->colourPartner()->status() == HardBranching::Incoming ) || (emitter->status() == HardBranching::Incoming && emitter->colourPartner()->status() == HardBranching::Outgoing)) { system.type=IF; // enusre incoming first if(system.jets[0]->status() == HardBranching::Outgoing) swap(system.jets[0],system.jets[1]); deconstructInitialFinalSystem(tree,system.jets,type); } else { throw Exception() << "Unknown type of system in " << "KinematicsReconstructor::deconstructColourPartner()" << Exception::runerror; } // only at this point that we can be sure all the reference vectors // are correct for(set::const_iterator it=tree->branchings().begin(); it!=tree->branchings().end();++it) { if((**it).status()==HardBranching::Incoming) continue; if((**it).branchingParticle()->coloured()) (**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false); } for(set::const_iterator it=tree->incoming().begin(); it!=tree->incoming().end();++it) { (**it).setMomenta(LorentzRotation(),1.,Lorentz5Momentum(),false); } for(set::const_iterator it=tree->branchings().begin(); it!=tree->branchings().end();++it) { if((**it).status()!=HardBranching::Incoming) continue; if(*it==system.jets[0] || *it==system.jets[1]) continue; if((**it).branchingParticle()->momentum().z()>ZERO) { (**it).z((**it).branchingParticle()->momentum().plus()/(**it).beam()->momentum().plus()); } else { (**it).z((**it).branchingParticle()->momentum().minus()/(**it).beam()->momentum().minus()); } } return true; } void KinematicsReconstructor:: reconstructInitialFinalSystem(vector jets) const { Lorentz5Momentum pin[2],pout[2],pbeam; for(unsigned int ix=0;ixprogenitor()->isFinalState()) { pout[0] +=jets[ix]->progenitor()->momentum(); - _progenitor = jets[ix]->progenitor(); if(jets[ix]->reconstructed()==ShowerProgenitor::notReconstructed) { - reconstructTimeLikeJet(jets[ix]->progenitor()); + reconstructTimeLikeJet(jets[ix]->progenitor(),jets[ix]->progenitor()); jets[ix]->reconstructed(ShowerProgenitor::done); } } // initial-state parton else { pin[0] +=jets[ix]->progenitor()->momentum(); if(jets[ix]->progenitor()->showerKinematics()) { pbeam = jets[ix]->progenitor()->showerBasis()->getBasis()[0]; } else { if ( jets[ix]->original()->parents().empty() ) { pbeam = jets[ix]->progenitor()->momentum(); } else { pbeam = jets[ix]->original()->parents()[0]->momentum(); } } if(jets[ix]->reconstructed()==ShowerProgenitor::notReconstructed) { reconstructSpaceLikeJet(jets[ix]->progenitor()); jets[ix]->reconstructed(ShowerProgenitor::done); } assert(!jets[ix]->original()->parents().empty()); } } // add intrinsic pt if needed addIntrinsicPt(jets); // momenta after showering for(unsigned int ix=0;ixprogenitor()->isFinalState()) pout[1] += jets[ix]->progenitor()->momentum(); else pin[1] += jets[ix]->progenitor()->momentum(); } // work out the boost to the Breit frame Lorentz5Momentum pa = pout[0]-pin[0]; Axis axis(pa.vect().unit()); LorentzRotation rot; double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); if ( sinth > 1.e-9 ) rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); rot.rotateX(Constants::pi); rot.boostZ( pa.e()/pa.vect().mag()); Lorentz5Momentum ptemp=rot*pbeam; Boost trans = -1./ptemp.e()*ptemp.vect(); trans.setZ(0.); if ( trans.mag2() - 1. >= 0. ) throw KinematicsReconstructionVeto(); rot.boost(trans); pa *=rot; // project and calculate rescaling // reference vectors Lorentz5Momentum n1(ZERO,ZERO,-pa.z(),-pa.z()); Lorentz5Momentum n2(ZERO,ZERO, pa.z(),-pa.z()); Energy2 n1n2 = n1*n2; // decompose the momenta Lorentz5Momentum qbp=rot*pin[1],qcp=rot*pout[1]; qbp.rescaleMass(); qcp.rescaleMass(); double a[2],b[2]; a[0] = n2*qbp/n1n2; b[0] = n1*qbp/n1n2; Lorentz5Momentum qperp = qbp-a[0]*n1-b[0]*n2; b[1] = 0.5; a[1] = 0.5*(qcp.m2()-qperp.m2())/n1n2/b[1]; double kb; if(a[0]!=0.) { double A(0.5*a[0]),B(b[0]*a[0]-a[1]*b[1]-0.25),C(-0.5*b[0]); if(sqr(B)-4.*A*C<0.) throw KinematicsReconstructionVeto(); kb = 0.5*(-B+sqrt(sqr(B)-4.*A*C))/A; } else { kb = 0.5*b[0]/(b[0]*a[0]-a[1]*b[1]-0.25); } // changed to improve stability if(kb==0.) throw KinematicsReconstructionVeto(); if ( a[1]>b[1] && abs(a[1]) < 1e-12 ) throw KinematicsReconstructionVeto(); if ( a[1]<=b[1] && abs(0.5+b[0]/kb) < 1e-12 ) throw KinematicsReconstructionVeto(); double kc = (a[1]>b[1]) ? (a[0]*kb-0.5)/a[1] : b[1]/(0.5+b[0]/kb); if(kc==0.) throw KinematicsReconstructionVeto(); Lorentz5Momentum pnew[2] = { a[0]*kb*n1+b[0]/kb*n2+qperp, a[1]*kc*n1+b[1]/kc*n2+qperp}; LorentzRotation rotinv=rot.inverse(); for(unsigned int ix=0;ixprogenitor()->isFinalState()) { deepTransform(jets[ix]->progenitor(),rot); deepTransform(jets[ix]->progenitor(),solveBoost(pnew[1],qcp)); Energy delta = jets[ix]->progenitor()->momentum().m()-jets[ix]->progenitor()->momentum().mass(); if ( abs(delta) > MeV ) throw KinematicsReconstructionVeto(); deepTransform(jets[ix]->progenitor(),rotinv); } else { tPPtr parent; boostChain(jets[ix]->progenitor(),rot,parent); boostChain(jets[ix]->progenitor(),solveBoostZ(pnew[0],qbp),parent); // check the first boost worked, and if not apply small correction to // fix energy/momentum conservation // this is a kludge but it reduces momentum non-conservation dramatically Lorentz5Momentum pdiff = pnew[0]-jets[ix]->progenitor()->momentum(); Energy2 delta = sqr(pdiff.x())+sqr(pdiff.y())+sqr(pdiff.z())+sqr(pdiff.t()); unsigned int ntry=0; while(delta>1e-6*GeV2 && ntry<5 ) { ntry +=1; boostChain(jets[ix]->progenitor(),solveBoostZ(pnew[0],jets[ix]->progenitor()->momentum()),parent); pdiff = pnew[0]-jets[ix]->progenitor()->momentum(); delta = sqr(pdiff.x())+sqr(pdiff.y())+sqr(pdiff.z())+sqr(pdiff.t()); } // apply test in breit-frame Lorentz5Momentum ptest1 = parent->momentum(); Lorentz5Momentum ptest2 = rot*pbeam; if(ptest1.z()/ptest2.z()<0. || ptest1.z()/ptest2.z()>1.) throw KinematicsReconstructionVeto(); boostChain(jets[ix]->progenitor(),rotinv,parent); } } } bool KinematicsReconstructor::addIntrinsicPt(vector jets) const { bool added=false; // add the intrinsic pt if needed for(unsigned int ix=0;ixprogenitor()->isFinalState()|| jets[ix]->hasEmitted()|| jets[ix]->reconstructed()==ShowerProgenitor::dontReconstruct) continue; if(_intrinsic.find(jets[ix])==_intrinsic.end()) continue; pair pt=_intrinsic[jets[ix]]; Energy etemp = jets[ix]->original()->parents()[0]->momentum().z(); Lorentz5Momentum p_basis(ZERO, ZERO, etemp, abs(etemp)), n_basis(ZERO, ZERO,-etemp, abs(etemp)); double alpha = jets[ix]->progenitor()->x(); double beta = 0.5*(sqr(jets[ix]->progenitor()->data().mass())+ sqr(pt.first))/alpha/(p_basis*n_basis); Lorentz5Momentum pnew=alpha*p_basis+beta*n_basis; pnew.setX(pt.first*cos(pt.second)); pnew.setY(pt.first*sin(pt.second)); pnew.rescaleMass(); jets[ix]->progenitor()->set5Momentum(pnew); added = true; } return added; } namespace { double defaultSolveBoostGamma(const double & betam,const Energy2 & kps, const Energy2 & qs, const Energy2 & Q2, const Energy & kp, const Energy & q, const Energy & qE) { if(betam<0.5) { return 1./sqrt(1.-sqr(betam)); } else { return ( kps+ qs + Q2)/ sqrt(2.*kps*qs + kps*Q2 + qs*Q2 + sqr(Q2) + 2.*q*qE*kp*sqrt(kps + Q2)); } } } LorentzRotation KinematicsReconstructor:: solveBoost(const double k, const Lorentz5Momentum & newq, const Lorentz5Momentum & oldp ) const { Energy q = newq.vect().mag(); Energy2 qs = sqr(q); Energy2 Q2 = newq.mass2(); Energy kp = k*(oldp.vect().mag()); Energy2 kps = sqr(kp); double betam = (q*newq.e() - kp*sqrt(kps + Q2))/(kps + qs + Q2); if ( abs(betam) - 1. >= 0. ) throw KinematicsReconstructionVeto(); Boost beta = -betam*(k/kp)*oldp.vect(); double gamma = 0.; if(Q2/sqr(oldp.e())>1e-4) { gamma = defaultSolveBoostGamma(betam,kps,qs,Q2,kp,q,newq.e()); } else { if(k>0) { gamma = 4.*kps*qs/sqr(kps +qs) + 2.*sqr(kps-qs)*Q2/pow<3,1>(kps +qs) - 0.25*( sqr(kps) + 14.*kps*qs + sqr(qs))*sqr(kps-qs)/(pow<4,1>(kps +qs)*kps*qs)*sqr(Q2); } else { gamma = 0.25*sqr(Q2)/(kps*qs)*(1. - 0.5*(kps+qs)/(kps*qs)*Q2); } if(gamma<=0.) throw KinematicsReconstructionVeto(); gamma = 1./sqrt(gamma); if(gamma>2.) gamma = defaultSolveBoostGamma(betam,kps,qs,Q2,kp,q,newq.e()); } // note that (k/kp)*oldp.vect() = oldp.vect()/oldp.vect().mag() but cheaper. ThreeVector ax = newq.vect().cross( oldp.vect() ); double delta; if (newq.x()*oldp.x()+newq.y()*oldp.y()+newq.z()*oldp.z()< 1e-16*GeV2) { throw KinematicsReconstructionVeto(); }else{ delta = newq.vect().angle( oldp.vect() ); } LorentzRotation R; using Constants::pi; Energy2 scale1 = sqr(newq.x())+ sqr(newq.y())+sqr(newq.z()); Energy2 scale2 = sqr(oldp.x())+ sqr(oldp.y())+sqr(oldp.z()); if ( ax.mag2()/scale1/scale2 > 1e-28 ) { R.rotate( delta, unitVector(ax) ).boost( beta , gamma ); } else if(abs(delta-pi)/pi < 0.001) { double phi=2.*pi*UseRandom::rnd(); Axis axis(cos(phi),sin(phi),0.); axis.rotateUz(newq.vect().unit()); R.rotate(delta,axis).boost( beta , gamma ); } else { R.boost( beta , gamma ); } return R; } LorentzRotation KinematicsReconstructor::solveBoost(const Lorentz5Momentum & q, const Lorentz5Momentum & p ) const { Energy modp = p.vect().mag(); Energy modq = q.vect().mag(); double betam = (p.e()*modp-q.e()*modq)/(sqr(modq)+sqr(modp)+p.mass2()); if ( abs(betam)-1. >= 0. ) throw KinematicsReconstructionVeto(); Boost beta = -betam*q.vect().unit(); ThreeVector ax = p.vect().cross( q.vect() ); double delta = p.vect().angle( q.vect() ); LorentzRotation R; using Constants::pi; if ( beta.mag2() - 1. >= 0. ) throw KinematicsReconstructionVeto(); if ( ax.mag2()/GeV2/MeV2 > 1e-16 ) { R.rotate( delta, unitVector(ax) ).boost( beta ); } else { R.boost( beta ); } return R; } LorentzRotation KinematicsReconstructor::solveBoostZ(const Lorentz5Momentum & q, const Lorentz5Momentum & p ) const { static const double eps = 1e-6; LorentzRotation R; double beta; Energy2 mt2 = p.mass()eps) { double erat = (q.t()+q.z())/(p.t()+p.z()); Energy2 den = mt2*(erat+1./erat); Energy2 num = (q.z()-p.z())*(q.t()+p.t()) + (p.z()+q.z())*(p.t()-q.t()); beta = num/den; if ( abs(beta) - 1. >= 0. ) throw KinematicsReconstructionVeto(); R.boostZ(beta); } else { double er = sqr(p.t()/q.t()); double x = ratio+0.125*(er+10.+1./er)*sqr(ratio); beta = -(p.t()-q.t())*(p.t()+q.t())/(sqr(p.t())+sqr(q.t()))*(1.+x); double gamma = (4.*sqr(p.t()*q.t()) +sqr(p.t()-q.t())*sqr(p.t()+q.t())* (-2.*x+sqr(x)))/sqr(sqr(p.t())+sqr(q.t())); if ( abs(beta) - 1. >= 0. ) throw KinematicsReconstructionVeto(); gamma = 1./sqrt(gamma); R.boost(0.,0.,beta,gamma); } Lorentz5Momentum ptest = R*p; if(ptest.z()/q.z() < 0. || ptest.t()/q.t() < 0. ) { throw KinematicsReconstructionVeto(); } return R; } void KinematicsReconstructor:: reconstructFinalStateSystem(bool applyBoost, const LorentzRotation & toRest, const LorentzRotation & fromRest, vector jets) const { LorentzRotation trans = applyBoost? toRest : LorentzRotation(); // special for case of individual particle if(jets.size()==1) { deepTransform(jets[0]->progenitor(),trans); deepTransform(jets[0]->progenitor(),fromRest); return; } bool radiated(false); // find the hard process centre-of-mass energy Lorentz5Momentum pcm; // check if radiated and calculate total momentum for(unsigned int ix=0;ixhasEmitted(); pcm += jets[ix]->progenitor()->momentum(); } if(applyBoost) pcm *= trans; // check if in CMF frame Boost beta_cm = pcm.findBoostToCM(); bool gottaBoost(false); if(beta_cm.mag() > 1e-12) { gottaBoost = true; trans.boost(beta_cm); } // collection of pointers to initial hard particle and jet momenta // for final boosts JetKinVect jetKinematics; vector::const_iterator cit; for(cit = jets.begin(); cit != jets.end(); cit++) { JetKinStruct tempJetKin; tempJetKin.parent = (*cit)->progenitor(); if(applyBoost || gottaBoost) { deepTransform(tempJetKin.parent,trans); } tempJetKin.p = (*cit)->progenitor()->momentum(); - _progenitor=tempJetKin.parent; if((**cit).reconstructed()==ShowerProgenitor::notReconstructed) { - radiated |= reconstructTimeLikeJet((*cit)->progenitor()); + radiated |= reconstructTimeLikeJet((*cit)->progenitor(),tempJetKin.parent); (**cit).reconstructed(ShowerProgenitor::done); } else { radiated |= !(*cit)->progenitor()->children().empty(); } tempJetKin.q = (*cit)->progenitor()->momentum(); jetKinematics.push_back(tempJetKin); } if(_finalFinalWeight && jetKinematics.size()==2) { Energy m1 = jetKinematics[0].q.m(); Energy m2 = jetKinematics[1].q.m(); Energy m0 = pcm.m(); if(m0sqrt(lambdaNew/lambdaOld)) throw KinematicsReconstructionVeto(); } // default option rescale everything with the same factor // find the rescaling factor double k = 0.0; if(radiated) { k = solveKfactor(pcm.m(), jetKinematics); // perform the rescaling and boosts for(JetKinVect::iterator it = jetKinematics.begin(); it != jetKinematics.end(); ++it) { LorentzRotation Trafo = solveBoost(k, it->q, it->p); deepTransform(it->parent,Trafo); } } // apply the final boosts if(gottaBoost || applyBoost) { LorentzRotation finalBoosts; if(gottaBoost) finalBoosts.boost(-beta_cm); if(applyBoost) finalBoosts.transform(fromRest); for(JetKinVect::iterator it = jetKinematics.begin(); it != jetKinematics.end(); ++it) { deepTransform(it->parent,finalBoosts); } } } void KinematicsReconstructor:: reconstructInitialInitialSystem(bool & applyBoost, LorentzRotation & toRest, LorentzRotation & fromRest, vector jets) const { bool radiated = false; Lorentz5Momentum pcm; // check whether particles radiated and calculate total momentum for( unsigned int ix = 0; ix < jets.size(); ++ix ) { radiated |= jets[ix]->hasEmitted(); pcm += jets[ix]->progenitor()->momentum(); if(jets[ix]->original()->parents().empty()) return; } pcm.rescaleMass(); // check if intrinsic pt to be added radiated |= !_intrinsic.empty(); // if no radiation return if(!radiated) { for(unsigned int ix=0;ixreconstructed()==ShowerProgenitor::notReconstructed) jets[ix]->reconstructed(ShowerProgenitor::done); } return; } // initial state shuffling applyBoost=false; vector p, pq, p_in; vector pts; for(unsigned int ix=0;ixprogenitor()->momentum()); // reconstruct the jet if(jets[ix]->reconstructed()==ShowerProgenitor::notReconstructed) { radiated |= reconstructSpaceLikeJet(jets[ix]->progenitor()); jets[ix]->reconstructed(ShowerProgenitor::done); } assert(!jets[ix]->original()->parents().empty()); Energy etemp = jets[ix]->original()->parents()[0]->momentum().z(); Lorentz5Momentum ptemp = Lorentz5Momentum(ZERO, ZERO, etemp, abs(etemp)); pq.push_back(ptemp); pts.push_back(jets[ix]->highestpT()); } // add the intrinsic pt if needed radiated |=addIntrinsicPt(jets); for(unsigned int ix=0;ixprogenitor()->momentum()); } double x1 = p_in[0].z()/pq[0].z(); double x2 = p_in[1].z()/pq[1].z(); vector beta=initialStateRescaling(x1,x2,p_in[0]+p_in[1],p,pq,pts); // if not need don't apply boosts if(!(radiated && p.size() == 2 && pq.size() == 2)) return; applyBoost=true; // apply the boosts Lorentz5Momentum newcmf; for(unsigned int ix=0;ixprogenitor(); Boost betaboost(0, 0, beta[ix]); tPPtr parent; boostChain(toBoost, LorentzRotation(0.,0.,beta[ix]),parent); if(parent->momentum().e()/pq[ix].e()>1.|| parent->momentum().z()/pq[ix].z()>1.) throw KinematicsReconstructionVeto(); newcmf+=toBoost->momentum(); } if(newcmf.m() jets, ShowerInteraction) const { assert(jets.size()==2); // put beam with +z first if(jets[0]->beam()->momentum().z() pin,pq; for(unsigned int ix=0;ixbranchingParticle()->momentum()); Energy etemp = jets[ix]->beam()->momentum().z(); pq.push_back(Lorentz5Momentum(ZERO, ZERO,etemp, abs(etemp))); } // calculate the rescaling double x[2]; Lorentz5Momentum pcm=pin[0]+pin[1]; assert(pcm.mass2()>ZERO); pcm.rescaleMass(); vector boost = inverseInitialStateRescaling(x[0],x[1],pcm,pin,pq); set::const_iterator cjt=tree->incoming().begin(); HardBranchingPtr incoming[2]; incoming[0] = *cjt; ++cjt; incoming[1] = *cjt; if((*tree->incoming().begin())->beam()->momentum().z()/pq[0].z()<0.) swap(incoming[0],incoming[1]); // apply the boost the the particles unsigned int iswap[2]={1,0}; for(unsigned int ix=0;ix<2;++ix) { LorentzRotation R(0.,0.,-boost[ix]); incoming[ix]->pVector(pq[ix]); incoming[ix]->nVector(pq[iswap[ix]]); incoming[ix]->setMomenta(R,1.,Lorentz5Momentum()); jets[ix]->showerMomentum(x[ix]*jets[ix]->pVector()); } // and calculate the boosts applyBoost=true; // do one boost if(_initialBoost==0) { toRest = LorentzRotation(-pcm.boostVector()); } else if(_initialBoost==1) { // first the transverse boost Energy pT = sqrt(sqr(pcm.x())+sqr(pcm.y())); double beta = -pT/pcm.t(); toRest=LorentzRotation(Boost(beta*pcm.x()/pT,beta*pcm.y()/pT,0.)); // the longitudinal beta = pcm.z()/sqrt(pcm.m2()+sqr(pcm.z())); toRest.boost(Boost(0.,0.,-beta)); } else assert(false); fromRest = LorentzRotation((jets[0]->showerMomentum()+ jets[1]->showerMomentum()).boostVector()); } void KinematicsReconstructor:: deconstructFinalStateSystem(const LorentzRotation & toRest, const LorentzRotation & fromRest, HardTreePtr tree, vector jets, ShowerInteraction type) const { LorentzRotation trans = toRest; if(jets.size()==1) { Lorentz5Momentum pnew = toRest*(jets[0]->branchingParticle()->momentum()); pnew *= fromRest; jets[0]-> original(pnew); jets[0]->showerMomentum(pnew); // find the colour partners ShowerParticleVector particles; vector ptemp; set::const_iterator cjt; for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) { ptemp.push_back((**cjt).branchingParticle()->momentum()); (**cjt).branchingParticle()->set5Momentum((**cjt).showerMomentum()); particles.push_back((**cjt).branchingParticle()); } dynamic_ptr_cast(ShowerHandler::currentHandler())->partnerFinder() ->setInitialEvolutionScales(particles,false,type,false); // calculate the reference vectors unsigned int iloc(0); set::iterator clt; for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) { // reset the momentum (**cjt).branchingParticle()->set5Momentum(ptemp[iloc]); ++iloc; // sort out the partners tShowerParticlePtr partner = (*cjt)->branchingParticle()->partner(); if(!partner) continue; for(clt=tree->branchings().begin();clt!=tree->branchings().end();++clt) { if((**clt).branchingParticle()==partner) { (**cjt).colourPartner(*clt); break; } } tHardBranchingPtr branch; for(clt=tree->branchings().begin();clt!=tree->branchings().end();++clt) { if(clt==cjt) continue; if((*clt)->branchingParticle()==partner) { branch=*clt; break; } } } return; } vector::iterator cit; vector pout; vector mon; Lorentz5Momentum pin; for(cit=jets.begin();cit!=jets.end();++cit) { pout.push_back((*cit)->branchingParticle()->momentum()); mon.push_back(findMass(*cit)); pin+=pout.back(); } // boost all the momenta to the rest frame of the decaying particle pin.rescaleMass(); pin *=trans; Boost beta_cm = pin.findBoostToCM(); bool gottaBoost(false); if(beta_cm.mag() > 1e-12) { gottaBoost = true; trans.boost(beta_cm); pin.boost(beta_cm); } for(unsigned int ix=0;ixbranchingParticle()->momentum(); pvect.transform(trans); pvect /= lambda; pvect.setMass(mon[ix]); pvect.rescaleEnergy(); if(gottaBoost) pvect.boost(-beta_cm); pvect.transform(fromRest); jets[ix]->pVector(pvect); jets[ix]->showerMomentum(pvect); } // find the colour partners ShowerParticleVector particles; vector ptemp; set::const_iterator cjt; for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) { ptemp.push_back((**cjt).branchingParticle()->momentum()); (**cjt).branchingParticle()->set5Momentum((**cjt).showerMomentum()); particles.push_back((**cjt).branchingParticle()); } dynamic_ptr_cast(ShowerHandler::currentHandler())->partnerFinder() ->setInitialEvolutionScales(particles,false,type,false); // calculate the reference vectors unsigned int iloc(0); set::iterator clt; for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) { // reset the momentum (**cjt).branchingParticle()->set5Momentum(ptemp[iloc]); ++iloc; } for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) { // sort out the partners tShowerParticlePtr partner = (*cjt)->branchingParticle()->partner(); if(!partner) continue; for(clt=tree->branchings().begin();clt!=tree->branchings().end();++clt) { if((**clt).branchingParticle()==partner) { (**cjt).colourPartner(*clt); break; } } tHardBranchingPtr branch; for(clt=tree->branchings().begin();clt!=tree->branchings().end();++clt) { if(clt==cjt) continue; if((*clt)->branchingParticle()==partner) { branch=*clt; break; } } // compute the reference vectors // both incoming, should all ready be done if((**cjt).status()==HardBranching::Incoming && (**clt).status()==HardBranching::Incoming) { continue; } // both outgoing else if((**cjt).status()!=HardBranching::Incoming&& branch->status()==HardBranching::Outgoing) { Boost boost=((*cjt)->pVector()+branch->pVector()).findBoostToCM(); Lorentz5Momentum pcm = branch->pVector(); pcm.boost(boost); Lorentz5Momentum nvect = Lorentz5Momentum(ZERO,pcm.vect()); nvect.boost( -boost); (**cjt).nVector(nvect); } else if((**cjt).status()==HardBranching::Incoming) { Lorentz5Momentum pa = -(**cjt).showerMomentum()+branch->showerMomentum(); Lorentz5Momentum pb = (**cjt).showerMomentum(); Axis axis(pa.vect().unit()); LorentzRotation rot; double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); rot.rotateX(Constants::pi); rot.boostZ( pa.e()/pa.vect().mag()); pb*=rot; Boost trans = -1./pb.e()*pb.vect(); trans.setZ(0.); rot.boost(trans); Energy scale=(**cjt).beam()->momentum().e(); Lorentz5Momentum pbasis(ZERO,(**cjt).beam()->momentum().vect().unit()*scale); Lorentz5Momentum pcm = rot*pbasis; rot.invert(); (**cjt).nVector(rot*Lorentz5Momentum(ZERO,-pcm.vect())); tHardBranchingPtr branch2 = *cjt;; while (branch2->parent()) { branch2=branch2->parent(); branch2->nVector(rot*Lorentz5Momentum(ZERO,-pcm.vect())); } } else if(branch->status()==HardBranching::Incoming) { (**cjt).nVector(Lorentz5Momentum(ZERO,branch->showerMomentum().vect())); } } // now compute the new momenta for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) { if(!(*cjt)->branchingParticle()->isFinalState()) continue; Lorentz5Momentum qnew; if((*cjt)->branchingParticle()->partner()) { Energy2 dot=(*cjt)->pVector()*(*cjt)->nVector(); double beta = 0.5*((*cjt)->branchingParticle()->momentum().m2() -sqr((*cjt)->pVector().mass()))/dot; qnew=(*cjt)->pVector()+beta*(*cjt)->nVector(); qnew.rescaleMass(); } else { qnew = (*cjt)->pVector(); } // qnew is the unshuffled momentum in the rest frame of the p basis vectors, // for the simple case Z->q qbar g this was checked against analytic formulae. // compute the boost LorentzRotation R=solveBoost(qnew, toRest*(*cjt)->branchingParticle()->momentum())*toRest; (*cjt)->setMomenta(R,1.0,Lorentz5Momentum()); } } Energy KinematicsReconstructor::momConsEq(double k, const Energy & root_s, const JetKinVect & jets) const { static const Energy2 eps=1e-8*GeV2; Energy dum = ZERO; for(JetKinVect::const_iterator it = jets.begin(); it != jets.end(); ++it) { Energy2 dum2 = (it->q).m2() + sqr(k)*(it->p).vect().mag2(); if(dum2 < ZERO) { if(dum2 < -eps) throw KinematicsReconstructionVeto(); dum2 = ZERO; } dum += sqrt(dum2); } return dum - root_s; } void KinematicsReconstructor::boostChain(tPPtr p, const LorentzRotation &bv, tPPtr & parent) const { if(!p->parents().empty()) boostChain(p->parents()[0], bv,parent); else parent=p; p->transform(bv); if(p->children().size()==2) { if(dynamic_ptr_cast(p->children()[1])) deepTransform(p->children()[1],bv); } } namespace { bool sortJets(ShowerProgenitorPtr j1, ShowerProgenitorPtr j2) { return j1->highestpT()>j2->highestpT(); } } void KinematicsReconstructor:: reconstructGeneralSystem(vector & ShowerHardJets) const { // find initial- and final-state systems ColourSingletSystem in,out; for(unsigned int ix=0;ixprogenitor()->isFinalState()) out.jets.push_back(ShowerHardJets[ix]); else in.jets.push_back(ShowerHardJets[ix]); } // reconstruct initial-initial system LorentzRotation toRest,fromRest; bool applyBoost(false); // reconstruct initial-initial system reconstructInitialInitialSystem(applyBoost,toRest,fromRest,in.jets); // reconstruct the final-state systems reconstructFinalStateSystem(applyBoost,toRest,fromRest,out.jets); } void KinematicsReconstructor:: reconstructFinalFirst(vector & ShowerHardJets) const { static const Energy2 minQ2 = 1e-4*GeV2; map used; for(unsigned int ix=0;ix outgoing; // first find any particles with final state partners for(unsigned int ix=0;ixprogenitor()->isFinalState()&& ShowerHardJets[ix]->progenitor()->partner()&& ShowerHardJets[ix]->progenitor()->partner()->isFinalState()) outgoing.insert(ShowerHardJets[ix]); } // then find the colour partners if(!outgoing.empty()) { set partners; for(set::const_iterator it=outgoing.begin();it!=outgoing.end();++it) { for(unsigned int ix=0;ixpartner()==ShowerHardJets[ix]->progenitor()) { partners.insert(ShowerHardJets[ix]); break; } } } outgoing.insert(partners.begin(),partners.end()); } // do the final-state reconstruction if needed if(!outgoing.empty()) { assert(outgoing.size()!=1); LorentzRotation toRest,fromRest; vector outgoingJets(outgoing.begin(),outgoing.end()); reconstructFinalStateSystem(false,toRest,fromRest,outgoingJets); } // Now do any initial-final systems which are needed vector IFSystems; // find the systems N.B. can have duplicates // find initial-state with FS partners or FS with IS partners for(unsigned int ix=0;ixprogenitor()->isFinalState()&& ShowerHardJets[ix]->progenitor()->partner()&& ShowerHardJets[ix]->progenitor()->partner()->isFinalState()) { IFSystems.push_back(ColourSingletSystem(IF,ShowerHardJets[ix])); } else if(ShowerHardJets[ix]->progenitor()->isFinalState()&& ShowerHardJets[ix]->progenitor()->partner()&& !ShowerHardJets[ix]->progenitor()->partner()->isFinalState()) { IFSystems.push_back(ColourSingletSystem(IF,ShowerHardJets[ix])); } } // then add the partners for(unsigned int is=0;isprogenitor()->partner()==ShowerHardJets[ix]->progenitor()) { IFSystems[is].jets.push_back(ShowerHardJets[ix]); } } // ensure incoming first if(IFSystems[is].jets[0]->progenitor()->isFinalState()) swap(IFSystems[is].jets[0],IFSystems[is].jets[1]); } if(!IFSystems.empty()) { unsigned int istart = UseRandom::irnd(IFSystems.size()); unsigned int istop=IFSystems.size(); for(unsigned int is=istart;is<=istop;++is) { if(is==IFSystems.size()) { if(istart!=0) { istop = istart-1; is=0; } else break; } // skip duplicates if(used[IFSystems[is].jets[0]] && used[IFSystems[is].jets[1]] ) continue; if(IFSystems[is].jets[0]->original()&&IFSystems[is].jets[0]->original()->parents().empty()) continue; Lorentz5Momentum psum; for(unsigned int ix=0;ixprogenitor()->isFinalState()) psum += IFSystems[is].jets[ix]->progenitor()->momentum(); else psum -= IFSystems[is].jets[ix]->progenitor()->momentum(); } if(-psum.m2()>minQ2) { reconstructInitialFinalSystem(IFSystems[is].jets); for(unsigned int ix=0;ixprogenitor()->isFinalState()) out.jets.push_back(ShowerHardJets[ix]); else in.jets.push_back(ShowerHardJets[ix]); } // reconstruct initial-initial system bool doRecon = false; for(unsigned int ix=0;ix & ShowerHardJets) const { static const Energy2 minQ2 = 1e-4*GeV2; // sort the vector by hardness of emission std::sort(ShowerHardJets.begin(),ShowerHardJets.end(),sortJets); // map between particles and progenitors for easy lookup map progenitorMap; for(unsigned int ix=0;ixprogenitor()] = ShowerHardJets[ix]; } // check that the IF systems can be reconstructed bool canReconstruct = true; for(unsigned int ix=0;ixprogenitor(); tShowerParticlePtr partner = progenitor->partner(); if(!partner) continue; else if((progenitor->isFinalState() && !partner->isFinalState()) || (!progenitor->isFinalState() && partner->isFinalState()) ) { vector jets(2); jets[0] = ShowerHardJets[ix]; jets[1] = progenitorMap[partner]; Lorentz5Momentum psum; for(unsigned int iy=0;iyprogenitor()->isFinalState()) psum += jets[iy]->progenitor()->momentum(); else psum -= jets[iy]->progenitor()->momentum(); } if(-psum.m2() used; for(unsigned int ix=0;ixreconstructed()==ShowerProgenitor::done) continue; // already reconstructed if(used[ShowerHardJets[ix]]) continue; // no partner continue tShowerParticlePtr progenitor = ShowerHardJets[ix]->progenitor(); tShowerParticlePtr partner = progenitor->partner(); if(!partner) { // check if there's a daughter tree which also needs boosting Lorentz5Momentum porig = progenitor->momentum(); map >::const_iterator tit; for(tit = _currentTree->treelinks().begin(); tit != _currentTree->treelinks().end();++tit) { // if there is, boost it if(tit->second.first && tit->second.second==progenitor) { Lorentz5Momentum pnew = tit->first->incomingLines().begin() ->first->progenitor()->momentum(); pnew *= tit->first->transform(); Lorentz5Momentum pdiff = porig-pnew; Energy2 test = sqr(pdiff.x()) + sqr(pdiff.y()) + sqr(pdiff.z()) + sqr(pdiff.t()); LorentzRotation rot; if(test>1e-6*GeV2) rot = solveBoost(porig,pnew); tit->first->transform(rot,false); _treeBoosts[tit->first].push_back(rot); } } ShowerHardJets[ix]->reconstructed(ShowerProgenitor::done); continue; } // do the reconstruction // final-final if(progenitor->isFinalState() && partner->isFinalState() ) { LorentzRotation toRest,fromRest; vector jets(2); jets[0] = ShowerHardJets[ix]; jets[1] = progenitorMap[partner]; if(_reconopt==4 && jets[1]->reconstructed()==ShowerProgenitor::notReconstructed) jets[1]->reconstructed(ShowerProgenitor::dontReconstruct); reconstructFinalStateSystem(false,toRest,fromRest,jets); if(_reconopt==4 && jets[1]->reconstructed()==ShowerProgenitor::dontReconstruct) jets[1]->reconstructed(ShowerProgenitor::notReconstructed); used[jets[0]] = true; if(_reconopt==3) used[jets[1]] = true; } // initial-final else if((progenitor->isFinalState() && !partner->isFinalState()) || (!progenitor->isFinalState() && partner->isFinalState()) ) { vector jets(2); jets[0] = ShowerHardJets[ix]; jets[1] = progenitorMap[partner]; if(jets[0]->progenitor()->isFinalState()) swap(jets[0],jets[1]); if(jets[0]->original()&&jets[0]->original()->parents().empty()) continue; Lorentz5Momentum psum; for(unsigned int iy=0;iyprogenitor()->isFinalState()) psum += jets[iy]->progenitor()->momentum(); else psum -= jets[iy]->progenitor()->momentum(); } if(_reconopt==4 && progenitorMap[partner]->reconstructed()==ShowerProgenitor::notReconstructed) progenitorMap[partner]->reconstructed(ShowerProgenitor::dontReconstruct); reconstructInitialFinalSystem(jets); if(_reconopt==4 && progenitorMap[partner]->reconstructed()==ShowerProgenitor::dontReconstruct) progenitorMap[partner]->reconstructed(ShowerProgenitor::notReconstructed); used[ShowerHardJets[ix]] = true; if(_reconopt==3) used[progenitorMap[partner]] = true; } // initial-initial else if(!progenitor->isFinalState() && !partner->isFinalState() ) { ColourSingletSystem in,out; in.jets.push_back(ShowerHardJets[ix]); in.jets.push_back(progenitorMap[partner]); for(unsigned int iy=0;iyprogenitor()->isFinalState()) out.jets.push_back(ShowerHardJets[iy]); } LorentzRotation toRest,fromRest; bool applyBoost(false); if(_reconopt==4 && in.jets[1]->reconstructed()==ShowerProgenitor::notReconstructed) in.jets[1]->reconstructed(ShowerProgenitor::dontReconstruct); reconstructInitialInitialSystem(applyBoost,toRest,fromRest,in.jets); if(_reconopt==4 && in.jets[1]->reconstructed()==ShowerProgenitor::dontReconstruct) in.jets[1]->reconstructed(ShowerProgenitor::notReconstructed); used[in.jets[0]] = true; if(_reconopt==3) used[in.jets[1]] = true; for(unsigned int iy=0;iyreconstructed()==ShowerProgenitor::notReconstructed) out.jets[iy]->reconstructed(ShowerProgenitor::dontReconstruct); } // reconstruct the final-state systems LorentzRotation finalBoosts; finalBoosts.transform( toRest); finalBoosts.transform(fromRest); for(unsigned int iy=0;iyprogenitor(),finalBoosts); } for(unsigned int iy=0;iyreconstructed()==ShowerProgenitor::dontReconstruct) out.jets[iy]->reconstructed(ShowerProgenitor::notReconstructed); } } } } bool KinematicsReconstructor:: inverseDecayRescalingFactor(vector pout, vector mon,Energy roots, Lorentz5Momentum ppartner, Energy mbar, double & k1, double & k2) const { ThreeVector qtotal; vector pmag; for(unsigned int ix=0;ix1e10) return false; } while (abs(numer)>eps&&itry<100); k1 = abs(k1); k2 = a*k1; return itry<100; } void KinematicsReconstructor:: deconstructInitialFinalSystem(HardTreePtr tree,vector jets, ShowerInteraction type) const { HardBranchingPtr incoming; Lorentz5Momentum pin[2],pout[2],pbeam; HardBranchingPtr initial; Energy mc(ZERO); for(unsigned int ix=0;ixstatus()==HardBranching::Outgoing) { pout[0] += jets[ix]->branchingParticle()->momentum(); mc = jets[ix]->branchingParticle()->thePEGBase() ? jets[ix]->branchingParticle()->thePEGBase()->mass() : jets[ix]->branchingParticle()->dataPtr()->mass(); } // initial-state parton else { pin[0] += jets[ix]->branchingParticle()->momentum(); initial = jets[ix]; pbeam = jets[ix]->beam()->momentum(); Energy scale=pbeam.t(); pbeam = Lorentz5Momentum(ZERO,pbeam.vect().unit()*scale); incoming = jets[ix]; while(incoming->parent()) incoming = incoming->parent(); } } if(jets.size()>2) { pout[0].rescaleMass(); mc = pout[0].mass(); } // work out the boost to the Breit frame Lorentz5Momentum pa = pout[0]-pin[0]; Axis axis(pa.vect().unit()); LorentzRotation rot; double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); if(axis.perp2()>0.) { rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); rot.rotateX(Constants::pi); rot.boostZ( pa.e()/pa.vect().mag()); } // transverse part Lorentz5Momentum paxis=rot*pbeam; Boost trans = -1./paxis.e()*paxis.vect(); trans.setZ(0.); rot.boost(trans); pa *= rot; // reference vectors Lorentz5Momentum n1(ZERO,ZERO,-pa.z(),-pa.z()); Lorentz5Momentum n2(ZERO,ZERO, pa.z(),-pa.z()); Energy2 n1n2 = n1*n2; // decompose the momenta Lorentz5Momentum qbp=rot*pin[0],qcp= rot*pout[0]; double a[2],b[2]; a[0] = n2*qbp/n1n2; b[0] = n1*qbp/n1n2; a[1] = n2*qcp/n1n2; b[1] = n1*qcp/n1n2; Lorentz5Momentum qperp = qbp-a[0]*n1-b[0]*n2; // before reshuffling Energy Q = abs(pa.z()); double c = sqr(mc/Q); Lorentz5Momentum pb(ZERO,ZERO,0.5*Q*(1.+c),0.5*Q*(1.+c)); Lorentz5Momentum pc(ZERO,ZERO,0.5*Q*(c-1.),0.5*Q*(1.+c)); double anew[2],bnew[2]; anew[0] = pb*n2/n1n2; bnew[0] = 0.5*(qbp.m2()-qperp.m2())/n1n2/anew[0]; bnew[1] = pc*n1/n1n2; anew[1] = 0.5*qcp.m2()/bnew[1]/n1n2; Lorentz5Momentum qnewb = (anew[0]*n1+bnew[0]*n2+qperp); Lorentz5Momentum qnewc = (anew[1]*n1+bnew[1]*n2); // initial-state boost LorentzRotation rotinv=rot.inverse(); LorentzRotation transb=rotinv*solveBoostZ(qnewb,qbp)*rot; // final-state boost LorentzRotation transc=rotinv*solveBoost(qnewc,qcp)*rot; // this will need changing for more than one outgoing particle // set the pvectors for(unsigned int ix=0;ixstatus()==HardBranching::Incoming) { jets[ix]->pVector(pbeam); jets[ix]->showerMomentum(rotinv*pb); incoming->pVector(jets[ix]->pVector()); } else { jets[ix]->pVector(rotinv*pc); jets[ix]->showerMomentum(jets[ix]->pVector()); } } // find the colour partners ShowerParticleVector particles; vector ptemp; set::const_iterator cjt; for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) { ptemp.push_back((**cjt).branchingParticle()->momentum()); (**cjt).branchingParticle()->set5Momentum((**cjt).showerMomentum()); particles.push_back((**cjt).branchingParticle()); } dynamic_ptr_cast(ShowerHandler::currentHandler())->partnerFinder() ->setInitialEvolutionScales(particles,false,type,false); unsigned int iloc(0); for(cjt=tree->branchings().begin();cjt!=tree->branchings().end();++cjt) { // reset the momentum (**cjt).branchingParticle()->set5Momentum(ptemp[iloc]); ++iloc; } for(vector::const_iterator cjt=jets.begin(); cjt!=jets.end();++cjt) { // sort out the partners tShowerParticlePtr partner = (*cjt)->branchingParticle()->partner(); if(!partner) continue; tHardBranchingPtr branch; for(set::const_iterator clt=tree->branchings().begin();clt!=tree->branchings().end();++clt) { if((**clt).branchingParticle()==partner) { (**cjt).colourPartner(*clt); branch=*clt; break; } } // compute the reference vectors // both incoming, should all ready be done if((**cjt).status()==HardBranching::Incoming && branch->status()==HardBranching::Incoming) { Energy etemp = (*cjt)->beam()->momentum().z(); Lorentz5Momentum nvect(ZERO, ZERO,-etemp, abs(etemp)); tHardBranchingPtr branch2 = *cjt; (**cjt).nVector(nvect); while (branch2->parent()) { branch2=branch2->parent(); branch2->nVector(nvect); } } // both outgoing else if((**cjt).status()==HardBranching::Outgoing&& branch->status()==HardBranching::Outgoing) { Boost boost=((*cjt)->pVector()+branch->pVector()).findBoostToCM(); Lorentz5Momentum pcm = branch->pVector(); pcm.boost(boost); Lorentz5Momentum nvect = Lorentz5Momentum(ZERO,pcm.vect()); nvect.boost( -boost); (**cjt).nVector(nvect); } else if((**cjt).status()==HardBranching::Incoming) { Lorentz5Momentum pa = -(**cjt).showerMomentum()+branch->showerMomentum(); Lorentz5Momentum pb = (**cjt).showerMomentum(); Axis axis(pa.vect().unit()); LorentzRotation rot; double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); if(axis.perp2()>1e-20) { rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); rot.rotateX(Constants::pi); } if(abs(1.-pa.e()/pa.vect().mag())>1e-6) rot.boostZ( pa.e()/pa.vect().mag()); pb*=rot; Boost trans = -1./pb.e()*pb.vect(); trans.setZ(0.); rot.boost(trans); Energy scale=(**cjt).beam()->momentum().t(); Lorentz5Momentum pbasis(ZERO,(**cjt).beam()->momentum().vect().unit()*scale); Lorentz5Momentum pcm = rot*pbasis; rot.invert(); Lorentz5Momentum nvect = rot*Lorentz5Momentum(ZERO,-pcm.vect()); (**cjt).nVector(nvect); tHardBranchingPtr branch2 = *cjt; while (branch2->parent()) { branch2=branch2->parent(); branch2->nVector(nvect); } } else if(branch->status()==HardBranching::Incoming) { Lorentz5Momentum nvect=Lorentz5Momentum(ZERO,branch->showerMomentum().vect()); (**cjt).nVector(nvect); } } // now compute the new momenta for(vector::const_iterator cjt=jets.begin(); cjt!=jets.end();++cjt) { if((**cjt).status()==HardBranching::Outgoing) { (**cjt).setMomenta(transc,1.,Lorentz5Momentum()); } } incoming->setMomenta(transb,1.,Lorentz5Momentum()); } void KinematicsReconstructor::deepTransform(PPtr particle, const LorentzRotation & r, bool match, PPtr original) const { if(_boosts.find(particle)!=_boosts.end()) { _boosts[particle].push_back(r); } Lorentz5Momentum porig = particle->momentum(); if(!original) original = particle; for ( int i = 0, N = particle->children().size(); i < N; ++i ) { deepTransform(particle->children()[i],r, particle->children()[i]->id()==original->id()&&match,original); } particle->transform(r); // transform the p and n vectors ShowerParticlePtr sparticle = dynamic_ptr_cast(particle); if(sparticle && sparticle->showerBasis()) { sparticle->showerBasis()->transform(r); } if ( particle->next() ) deepTransform(particle->next(),r,match,original); if(!match) return; if(!particle->children().empty()) return; // force the mass shell if(particle->dataPtr()->stable()) { Lorentz5Momentum ptemp = particle->momentum(); ptemp.rescaleEnergy(); particle->set5Momentum(ptemp); } // check if there's a daughter tree which also needs boosting + if(!_currentTree) return; map >::const_iterator tit; for(tit = _currentTree->treelinks().begin(); tit != _currentTree->treelinks().end();++tit) { // if there is, boost it if(tit->second.first && tit->second.second==original) { Lorentz5Momentum pnew = tit->first->incomingLines().begin() ->first->progenitor()->momentum(); pnew *= tit->first->transform(); Lorentz5Momentum pdiff = porig-pnew; Energy2 test = sqr(pdiff.x()) + sqr(pdiff.y()) + sqr(pdiff.z()) + sqr(pdiff.t()); LorentzRotation rot; if(test>1e-6*GeV2) rot = solveBoost(porig,pnew); tit->first->transform(r*rot,false); _treeBoosts[tit->first].push_back(r*rot); } } } Energy KinematicsReconstructor::findMass(HardBranchingPtr branch) const { // KH - 230909 - If the particle has no children then it will // not have showered and so it should be "on-shell" so we can // get it's mass from it's momentum. This means that the // inverseRescalingFactor doesn't give any nans or do things // it shouldn't if it gets e.g. two Z bosons generated with // off-shell masses. This is for sure not the best solution. // PR 1/1/10 modification to previous soln // PR 28/8/14 change to procedure and factorize into a function if(branch->children().empty()) { return branch->branchingParticle()->mass(); } else if(!branch->children().empty() && !branch->branchingParticle()->dataPtr()->stable() ) { for(unsigned int ix=0;ixchildren().size();++ix) { if(branch->branchingParticle()->id()== branch->children()[ix]->branchingParticle()->id()) return findMass(branch->children()[ix]); } } return branch->branchingParticle()->dataPtr()->mass(); } vector KinematicsReconstructor::inverseInitialStateRescaling(double & x1, double & x2, const Lorentz5Momentum & pold, const vector & p, const vector & pq) const { // hadronic CMS Energy2 s = (pq[0] +pq[1] ).m2(); // partonic CMS Energy MDY = pold.m(); // find alpha, beta and pt Energy2 p12=pq[0]*pq[1]; double a[2],b[2]; Lorentz5Momentum pt[2]; for(unsigned int ix=0;ix<2;++ix) { a[ix] = p[ix]*pq[1]/p12; b [ix] = p[ix]*pq[0]/p12; pt[ix] = p[ix]-a[ix]*pq[0]-b[ix]*pq[1]; } // compute kappa // we always want to preserve the mass of the system double k1(1.),k2(1.); if(_initialStateReconOption==0) { double rap=pold.rapidity(); x2 = MDY/sqrt(s*exp(2.*rap)); x1 = sqr(MDY)/s/x2; k1=a[0]/x1; k2=b[1]/x2; } // longitudinal momentum else if(_initialStateReconOption==1) { double A = 1.; double C = -sqr(MDY)/s; double B = 2.*pold.z()/sqrt(s); if(abs(B)>1e-10) { double discrim = 1.-4.*A*C/sqr(B); if(discrim < 0.) throw KinematicsReconstructionVeto(); x1 = B>0. ? 0.5*B/A*(1.+sqrt(discrim)) : 0.5*B/A*(1.-sqrt(discrim)); } else { x1 = -C/A; if( x1 <= 0.) throw KinematicsReconstructionVeto(); x1 = sqrt(x1); } x2 = sqr(MDY)/s/x1; k1=a[0]/x1; k2=b[1]/x2; } // preserve mass and don't scale the softer system // to reproduce the dipole kinematics else if(_initialStateReconOption==2) { // in this case kp = k1 or k2 depending on who's the harder guy k1 = a[0]*b[1]*s/sqr(MDY); if ( pt[0].perp2() < pt[1].perp2() ) swap(k1,k2); x1 = a[0]/k1; x2 = b[1]/k2; } else assert(false); // decompose the momenta double anew[2] = {a[0]/k1,a[1]*k2}; double bnew[2] = {b[0]*k1,b[1]/k2}; vector boost(2); for(unsigned int ix=0;ix<2;++ix) { boost[ix] = getBeta(a [ix]+b [ix], a[ix] -b [ix], anew[ix]+bnew[ix], anew[ix]-bnew[ix]); } return boost; } vector KinematicsReconstructor::initialStateRescaling(double x1, double x2, const Lorentz5Momentum & pold, const vector & p, const vector & pq, const vector& highestpts) const { Energy2 S = (pq[0]+pq[1]).m2(); // find alphas and betas in terms of desired basis Energy2 p12 = pq[0]*pq[1]; double a[2] = {p[0]*pq[1]/p12,p[1]*pq[1]/p12}; double b[2] = {p[0]*pq[0]/p12,p[1]*pq[0]/p12}; Lorentz5Momentum p1p = p[0] - a[0]*pq[0] - b[0]*pq[1]; Lorentz5Momentum p2p = p[1] - a[1]*pq[0] - b[1]*pq[1]; // compute kappa // we always want to preserve the mass of the system Energy MDY = pold.m(); Energy2 A = a[0]*b[1]*S; Energy2 B = Energy2(sqr(MDY)) - (a[0]*b[0]+a[1]*b[1])*S - (p1p+p2p).m2(); Energy2 C = a[1]*b[0]*S; double rad = 1.-4.*A*C/sqr(B); if(rad < 0.) throw KinematicsReconstructionVeto(); double kp = B/(2.*A)*(1.+sqrt(rad)); // now compute k1 // conserve rapidity double k1(0.); double k2(0.); if(_initialStateReconOption==0) { rad = kp*(b[0]+kp*b[1])/(kp*a[0]+a[1]); rad *= pq[0].z()1e-10) { double discrim = 1.-4.*a2*c2/sqr(b2); if(discrim < 0.) throw KinematicsReconstructionVeto(); k1 = b2>0. ? 0.5*b2/a2*(1.+sqrt(discrim)) : 0.5*b2/a2*(1.-sqrt(discrim)); } else { k1 = -c2/a2; if( k1 <= 0.) throw KinematicsReconstructionVeto(); k1 = sqrt(k1); } k2 = kp/k1; } // preserve mass and don't scale the softer system // to reproduce the dipole kinematics else if(_initialStateReconOption==2) { // in this case kp = k1 or k2 depending on who's the harder guy k1 = kp; k2 = 1.; if ( highestpts[0] < highestpts[1] ) swap(k1,k2); } else assert(false); // calculate the boosts vector beta(2); beta[0] = getBeta((a[0]+b[0]), (a[0]-b[0]), (k1*a[0]+b[0]/k1), (k1*a[0]-b[0]/k1)); beta[1] = getBeta((a[1]+b[1]), (a[1]-b[1]), (a[1]/k2+k2*b[1]), (a[1]/k2-k2*b[1])); if (pq[0].z() > ZERO) { beta[0] = -beta[0]; beta[1] = -beta[1]; } return beta; } void KinematicsReconstructor:: reconstructColourSinglets(vector & ShowerHardJets, ShowerInteraction type) const { // identify and catagorize the colour singlet systems unsigned int nnun(0),nnii(0),nnif(0),nnf(0),nni(0); vector systems(identifySystems(set(ShowerHardJets.begin(),ShowerHardJets.end()), nnun,nnii,nnif,nnf,nni)); // now decide what to do // initial-initial connection and final-state colour singlet systems LorentzRotation toRest,fromRest; bool applyBoost(false),general(false); // Drell-Yan type if(nnun==0&&nnii==1&&nnif==0&&nnf>0&&nni==0) { // reconstruct initial-initial system for(unsigned int ix=0;ix0&&nni==1)|| (nnif==2&& nni==0))) { // check these systems can be reconstructed for(unsigned int ix=0;ixprogenitor()->isFinalState()) q += systems[ix].jets[iy]->progenitor()->momentum(); else q -= systems[ix].jets[iy]->progenitor()->momentum(); } q.rescaleMass(); // check above cut if(abs(q.m())>=_minQ) continue; if(nnif==1&&nni==1) { throw KinematicsReconstructionVeto(); } else { general = true; break; } } if(!general) { for(unsigned int ix=0;ix0&&nni==2) { general = type!=ShowerInteraction::QCD; } // general type else { general = true; } // final-state systems except for general recon if(!general) { for(unsigned int ix=0;ix namespace Herwig { using namespace ThePEG; /**\ingroup Shower * Exception class * used to communicate failure of kinematics * reconstruction. */ struct KinematicsReconstructionVeto {}; /** \ingroup Shower * A simple struct to store the information we need on the * showering */ struct JetKinStruct { /** * Parent particle of the jet */ tShowerParticlePtr parent; /** * Momentum of the particle before reconstruction */ Lorentz5Momentum p; /** * Momentum of the particle after reconstruction */ Lorentz5Momentum q; }; /** * typedef for a vector of JetKinStruct */ typedef vector JetKinVect; /** * Enum to identify types of colour singlet systems */ enum SystemType { UNDEFINED=-1, II, IF, F ,I }; /** * Struct to store colour singlets */ template struct ColourSinglet { typedef vector > VecType; ColourSinglet() : type(UNDEFINED) {} ColourSinglet(SystemType intype,Value inpart) : type(intype),jets(1,inpart) {} /** * The type of system */ SystemType type; /** * The jets in the system */ vector jets; }; /** * Struct to store a colour singlet system of particles */ typedef ColourSinglet ColourSingletSystem; /** * Struct to store a colour singlet shower */ typedef ColourSinglet ColourSingletShower; /** \ingroup Shower * * This class is responsible for the kinematical reconstruction * after each showering step, and also for the necessary Lorentz boosts * in order to preserve energy-momentum conservation in the overall collision, * and also the invariant mass and the rapidity of the hard subprocess system. * In the case of multi-step showering, there will be not unnecessary * kinematical reconstructions. * * There is also the option of taking a set of momenta for the particles * and inverting the reconstruction to give the evolution variables for the * shower. * * Notice: * - although we often use the term "jet" in either methods or variables names, * or in comments, which could appear applicable only for QCD showering, * there is indeed no "dynamics" represented in this class: only kinematics * is involved, as the name of this class remainds. Therefore it can be used * for any kind of showers (QCD-,QED-,EWK-,... bremsstrahlung). * * @see ShowerParticle * @see ShowerKinematics * @see \ref KinematicsReconstructorInterfaces "The interfaces" * defined for KinematicsReconstructor. */ class KinematicsReconstructor: public Interfaced { public: /** * Default constructor */ KinematicsReconstructor() : _reconopt(0), _initialBoost(0), _finalStateReconOption(0), _initialStateReconOption(0), _finalFinalWeight(false), _minQ(MeV) {}; /** * Methods to reconstruct the kinematics of a scattering or decay process */ //@{ /** * Given in input a vector of the particles which initiated the showers * the method does the reconstruction of such jets, * including the appropriate boosts (kinematics reshufflings) * needed to conserve the total energy-momentum of the collision * and preserving the invariant mass and the rapidity of the * hard subprocess system. */ virtual bool reconstructHardJets(ShowerTreePtr hard, const map > & pt, ShowerInteraction type, bool switchRecon) const; /** * Given in input a vector of the particles which initiated the showers * the method does the reconstruction of such jets, * including the appropriate boosts (kinematics reshufflings) * needed to conserve the total energy-momentum of the collision * and preserving the invariant mass and the rapidity of the * hard subprocess system. */ virtual bool reconstructDecayJets(ShowerTreePtr decay, ShowerInteraction type) const; //@} /** * Methods to invert the reconstruction of the shower for * a scattering or decay process and calculate * the variables used to generate the * shower given the particles produced. * This is needed for the CKKW and POWHEG approaches */ //@{ /** * Given the particles, with a history which we wish to interpret * as a shower reconstruct the variables used to generate the * shower */ virtual bool deconstructDecayJets(HardTreePtr,ShowerInteraction) const; /** * Given the particles, with a history which we wish to interpret * as a shower reconstruct the variables used to generate the shower * for a hard process */ virtual bool deconstructHardJets(HardTreePtr,ShowerInteraction) const; //@} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); -protected: - - /** - * Methods to reconstruct the kinematics of individual jets - */ - //@{ +public: + /** * Given the particle (ShowerParticle object) that * originates a forward (time-like) jet, this method reconstructs the kinematics * of the jet. That is, by starting from the final grand-children (which * originates directly or indirectly from particleJetParent, * and which don't have children), and moving "backwards" (in a physical * time picture), towards the particleJetParent, the * ShowerKinematics objects associated with the various particles, * which have been created during the showering, are now completed. * In particular, at the end, we get the mass of the jet, which is the * main information we want. * This methods returns false if there was no radiation or rescaling required */ - virtual bool reconstructTimeLikeJet(const tShowerParticlePtr particleJetParent) const; + virtual bool reconstructTimeLikeJet(const tShowerParticlePtr particleJetParent, + const tShowerParticlePtr progenitor) const; + + /** + * Given a vector of 5-momenta of jets, where the 3-momenta are the initial + * ones before showering and the masses are reconstructed after the showering, + * this method returns the overall scaling factor for the 3-momenta of the + * vector of particles, vec{P}_i -> k * vec{P}_i, such to preserve energy- + * momentum conservation, i.e. after the rescaling the center of mass 5-momentum + * is equal to the one specified in input, cmMomentum. + * The method returns 0 if such factor cannot be found. + * @param root_s Centre-of-mass energy + * @param jets The jets + */ + double solveKfactor( const Energy & root_s, const JetKinVect & jets ) const; + + /** + * Compute the boost to get from the the old momentum to the new + */ + LorentzRotation solveBoost(const double k, + const Lorentz5Momentum & newq, + const Lorentz5Momentum & oldp) const; + + /** + * Apply a transform to the particle and any child, including child ShowerTree + * objects + * @param particle The particle + * @param r The Lorentz transformation + * @param match Whether or not to look at children etc + * @param original The original particle + */ + void deepTransform(PPtr particle,const LorentzRotation & r, + bool match=true,PPtr original=PPtr()) const; + +protected: + + /** + * Methods to reconstruct the kinematics of individual jets + */ + //@{ /** * Exactly similar to the previous one, but for a space-like jet. * Also in this case we start from the final grand-children (which * are childless) of the particle which originates the jet, but in * this case we proceed "forward" (in the physical time picture) * towards the particleJetParent. * This methods returns false if there was no radiation or rescaling required */ bool reconstructSpaceLikeJet(const tShowerParticlePtr particleJetParent) const; /** * Exactly similar to the previous one, but for a decay jet * This methods returns false if there was no radiation or rescaling required */ bool reconstructDecayJet(const tShowerParticlePtr particleJetParent) const; //@} /** * Methods to perform the reconstruction of various types of colour * singlet systems */ //@{ /** * Perform the reconstruction of a system with one incoming and at least one * outgoing particle */ void reconstructInitialFinalSystem(vector) const; /** * Perform the reconstruction of a system with only final-state * particles */ void reconstructFinalStateSystem(bool applyBoost, const LorentzRotation & toRest, const LorentzRotation & fromRest, vector) const; /** * Reconstruction of a general coloured system */ void reconstructGeneralSystem(vector & ShowerHardJets) const; /** * Reconstruction of a general coloured system doing * final-final, then initial-final and then initial-initial */ void reconstructFinalFirst(vector & ShowerHardJets) const; /** * Reconstruction of a general coloured system doing * colour parners */ void reconstructColourPartner(vector & ShowerHardJets) const; /** * Reconstruction based on colour singlet systems */ void reconstructColourSinglets(vector & ShowerHardJets, ShowerInteraction type) const; /** * Perform the reconstruction of a system with only final-state * particles */ void reconstructInitialInitialSystem(bool & applyBoost, LorentzRotation & toRest, LorentzRotation & fromRest, vector) const; //@} /** * Methods to perform the inverse reconstruction of various types of * colour singlet systems */ //@{ /** * Perform the inverse reconstruction of a system with only final-state * particles */ void deconstructFinalStateSystem(const LorentzRotation & toRest, const LorentzRotation & fromRest, HardTreePtr, vector, ShowerInteraction) const; /** * Perform the inverse reconstruction of a system with only initial-state * particles */ void deconstructInitialInitialSystem(bool & applyBoost, LorentzRotation & toRest, LorentzRotation & fromRest, HardTreePtr, vector, ShowerInteraction ) const; /** * Perform the inverse reconstruction of a system with only initial-state * particles */ void deconstructInitialFinalSystem(HardTreePtr, vector, ShowerInteraction ) const; bool deconstructGeneralSystem(HardTreePtr, ShowerInteraction) const; bool deconstructColourSinglets(HardTreePtr, ShowerInteraction) const; bool deconstructColourPartner(HardTreePtr, ShowerInteraction) const; //@} /** * Recursively treat the most off-shell paricle seperately * for final-final reconstruction */ void reconstructFinalFinalOffShell(JetKinVect orderedJets, Energy2 s, bool recursive) const; /** * Various methods for the Lorentz transforms needed to do the * rescalings */ //@{ /** * Compute the boost to get from the the old momentum to the new */ - LorentzRotation solveBoost(const double k, - const Lorentz5Momentum & newq, - const Lorentz5Momentum & oldp) const; - - /** - * Compute the boost to get from the the old momentum to the new - */ LorentzRotation solveBoost(const Lorentz5Momentum & newq, const Lorentz5Momentum & oldq) const; /** * Compute the boost to get from the the old momentum to the new */ LorentzRotation solveBoostZ(const Lorentz5Momentum & newq, const Lorentz5Momentum & oldq) const; /** * Recursively boost the initial-state shower * @param p The particle * @param bv The boost * @param parent The parent of the chain */ void boostChain(tPPtr p, const LorentzRotation & bv, tPPtr & parent) const; /** * Given a 5-momentum and a scale factor, the method returns the * Lorentz boost that transforms the 3-vector vec{momentum} ---> * k*vec{momentum}. The method returns the null boost in the case no * solution exists. This will only work in the case where the * outgoing jet-momenta are parallel to the momenta of the particles * leaving the hard subprocess. */ Boost solveBoostBeta( const double k, const Lorentz5Momentum & newq, const Lorentz5Momentum & oldp); /** * Compute boost parameter along z axis to get (Ep, any perp, qp) * from (E, same perp, q). */ double getBeta(const double E, const double q, const double Ep, const double qp) const {return (q*E-qp*Ep)/(sqr(qp)+sqr(E));} //@} /** * Methods to calculate the various scaling factors */ //@{ - /** - * Given a vector of 5-momenta of jets, where the 3-momenta are the initial - * ones before showering and the masses are reconstructed after the showering, - * this method returns the overall scaling factor for the 3-momenta of the - * vector of particles, vec{P}_i -> k * vec{P}_i, such to preserve energy- - * momentum conservation, i.e. after the rescaling the center of mass 5-momentum - * is equal to the one specified in input, cmMomentum. - * The method returns 0 if such factor cannot be found. - * @param root_s Centre-of-mass energy - * @param jets The jets - */ - double solveKfactor( const Energy & root_s, const JetKinVect & jets ) const; /** * Calculate the rescaling factors for the jets in a particle decay where * there was initial-state radiation * @param mb The mass of the decaying particle * @param n The reference vector for the initial state radiation * @param pjet The momentum of the initial-state jet * @param jetKinematics The JetKinStruct objects for the jets * @param partner The colour partner * @param ppartner The momentum of the colour partner of the decaying particle * before and after radiation * @param k1 The rescaling parameter for the partner * @param k2 The rescaling parameter for the outgoing singlet * @param qt The transverse momentum vector */ bool solveDecayKFactor(Energy mb, const Lorentz5Momentum & n, const Lorentz5Momentum & pjet, const JetKinVect & jetKinematics, ShowerParticlePtr partner, Lorentz5Momentum ppartner[2], double & k1, double & k2, Lorentz5Momentum & qt) const; /** * Compute the momentum rescaling factor needed to invert the shower * @param pout The momenta of the outgoing particles * @param mon The on-shell masses * @param roots The mass of the decaying particle */ double inverseRescalingFactor(vector pout, vector mon,Energy roots) const; /** * Compute the momentum rescaling factor needed to invert the shower * @param pout The momenta of the outgoing particles * @param mon The on-shell masses * @param roots The mass of the decaying particle * @param ppartner The momentum of the colour partner * @param mbar The mass of the decaying particle * @param k1 The first scaling factor * @param k2 The second scaling factor */ bool inverseDecayRescalingFactor(vector pout, vector mon,Energy roots, Lorentz5Momentum ppartner, Energy mbar, double & k1, double & k2) const; /** * Check the rescaling conserves momentum * @param k The rescaling * @param root_s The centre-of-mass energy * @param jets The jets */ Energy momConsEq(double k, const Energy & root_s, const JetKinVect & jets) const; void findInitialBoost(const Lorentz5Momentum & pold, const Lorentz5Momentum & pnew, LorentzRotation & toRest, LorentzRotation & fromRest) const; //@} /** * Find the colour partners of a particle to identify the colour singlet * systems for the reconstruction. */ template void findPartners(Value branch,set & done, const set & branchings, vector & jets) const; /** * Add the intrinsic \f$p_T\f$ to the system if needed */ bool addIntrinsicPt(vector) const; /** - * Apply a transform to the particle and any child, including child ShowerTree - * objects - * @param particle The particle - * @param r The Lorentz transformation - * @param match Whether or not to look at children etc - * @param original The original particle - */ - void deepTransform(PPtr particle,const LorentzRotation & r, - bool match=true,PPtr original=PPtr()) const; - - /** * Find the mass of a particle in the hard branching */ Energy findMass(HardBranchingPtr) const; /** * Calculate the initial-state rescaling factors */ vector initialStateRescaling(double x1, double x2, const Lorentz5Momentum & pold, const vector & p, const vector & pq, const vector& highespts) const; /** * Calculate the inverse of the initial-state rescaling factor */ vector inverseInitialStateRescaling(double & x1, double & x2, const Lorentz5Momentum & pold, const vector & p, const vector & pq) const; /** * Find the colour singlet systems */ template typename ColourSinglet::VecType identifySystems(set jets, unsigned int & nnun,unsigned int & nnii, unsigned int & nnif,unsigned int & nnf, unsigned int & nni) const; /** * Combine final-state colour systems */ template void combineFinalState(vector > & systems) const; protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const {return new_ptr(*this);} /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const {return new_ptr(*this);} //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ KinematicsReconstructor & operator=(const KinematicsReconstructor &) = delete; private: /** * Option for handling the reconstruction */ unsigned int _reconopt; /** * Option for the boost for initial-initial reconstruction */ unsigned int _initialBoost; /** * Option for the reconstruction of final state systems */ unsigned int _finalStateReconOption; /** * Option for the initial state reconstruction */ unsigned int _initialStateReconOption; /** * Option for FF kinematic factor */ bool _finalFinalWeight; /** * Minimum invariant mass for initial-final dipoles to allow the * reconstruction */ Energy _minQ; /** - * The progenitor of the jet currently being reconstructed - */ - mutable tShowerParticlePtr _progenitor; - - /** * Storage of the intrinsic \f$p_T\f$ */ mutable map > _intrinsic; /** * Current ShowerTree */ mutable tShowerTreePtr _currentTree; /** * Particles which shouldn't have their masses rescaled as * vector for the interface */ PDVector _noRescaleVector; /** * Particles which shouldn't have their masses rescaled as * set for quick access */ set _noRescale; /** * Storage of the boosts applied to enable resetting after failure */ mutable map > _boosts; /** * Storage of the boosts applied to enable resetting after failure */ mutable map > _treeBoosts; }; } #endif /* HERWIG_KinematicsReconstructor_H */ diff --git a/Shower/QTilde/SplittingFunctions/SudakovFormFactor.cc b/Shower/QTilde/SplittingFunctions/SudakovFormFactor.cc --- a/Shower/QTilde/SplittingFunctions/SudakovFormFactor.cc +++ b/Shower/QTilde/SplittingFunctions/SudakovFormFactor.cc @@ -1,1225 +1,1232 @@ // -*- C++ -*- // // SudakovFormFactor.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the SudakovFormFactor class. // #include "SudakovFormFactor.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "Herwig/Shower/QTilde/Kinematics/ShowerKinematics.h" #include "Herwig/Shower/QTilde/Base/ShowerParticle.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Shower/QTilde/QTildeShowerHandler.h" #include "Herwig/Shower/QTilde/Kinematics/FS_QTildeShowerKinematics1to2.h" #include "Herwig/Shower/QTilde/Kinematics/IS_QTildeShowerKinematics1to2.h" #include "Herwig/Shower/QTilde/Kinematics/Decay_QTildeShowerKinematics1to2.h" #include "Herwig/Shower/QTilde/Kinematics/KinematicHelpers.h" #include "SudakovCutOff.h" #include using std::array; using namespace Herwig; DescribeClass describeSudakovFormFactor ("Herwig::SudakovFormFactor",""); void SudakovFormFactor::persistentOutput(PersistentOStream & os) const { os << splittingFn_ << alpha_ << pdfmax_ << particles_ << pdffactor_ << cutoff_; } void SudakovFormFactor::persistentInput(PersistentIStream & is, int) { is >> splittingFn_ >> alpha_ >> pdfmax_ >> particles_ >> pdffactor_ >> cutoff_; } void SudakovFormFactor::Init() { static ClassDocumentation documentation ("The SudakovFormFactor class is the base class for the implementation of Sudakov" " form factors in Herwig"); static Reference interfaceSplittingFunction("SplittingFunction", "A reference to the SplittingFunction object", &Herwig::SudakovFormFactor::splittingFn_, false, false, true, false); static Reference interfaceAlpha("Alpha", "A reference to the Alpha object", &Herwig::SudakovFormFactor::alpha_, false, false, true, false); static Reference interfaceCutoff("Cutoff", "A reference to the SudakovCutOff object", &Herwig::SudakovFormFactor::cutoff_, false, false, true, false); static Parameter interfacePDFmax ("PDFmax", "Maximum value of PDF weight. ", &SudakovFormFactor::pdfmax_, 35.0, 1.0, 1000000.0, false, false, Interface::limited); static Switch interfacePDFFactor ("PDFFactor", "Include additional factors in the overestimate for the PDFs", &SudakovFormFactor::pdffactor_, 0, false, false); static SwitchOption interfacePDFFactorNo (interfacePDFFactor, "No", "Don't include any factors", 0); static SwitchOption interfacePDFFactorOverZ (interfacePDFFactor, "OverZ", "Include an additional factor of 1/z", 1); static SwitchOption interfacePDFFactorOverOneMinusZ (interfacePDFFactor, "OverOneMinusZ", "Include an additional factor of 1/(1-z)", 2); static SwitchOption interfacePDFFactorOverZOneMinusZ (interfacePDFFactor, "OverZOneMinusZ", "Include an additional factor of 1/z/(1-z)", 3); static SwitchOption interfacePDFFactorOverRootZ (interfacePDFFactor, "OverRootZ", "Include an additional factor of 1/sqrt(z)", 4); static SwitchOption interfacePDFFactorRootZ (interfacePDFFactor, "RootZ", "Include an additional factor of sqrt(z)", 5); } bool SudakovFormFactor::alphaSVeto(Energy2 pt2) const { double ratio=alphaSVetoRatio(pt2,1.); return UseRandom::rnd() > ratio; } double SudakovFormFactor::alphaSVetoRatio(Energy2 pt2, double factor) const { - factor *= ShowerHandler::currentHandler()->renormalizationScaleFactor(); + if(ShowerHandler::currentHandlerIsSet()) + factor *= ShowerHandler::currentHandler()->renormalizationScaleFactor(); return alpha_->ratio(pt2, factor); } bool SudakovFormFactor::PDFVeto(const Energy2 t, const double x, const tcPDPtr parton0, const tcPDPtr parton1, Ptr::transient_const_pointer beam) const { double ratio=PDFVetoRatio(t,x,parton0,parton1,beam,1.); return UseRandom::rnd() > ratio; } double SudakovFormFactor::PDFVetoRatio(const Energy2 t, const double x, const tcPDPtr parton0, const tcPDPtr parton1, Ptr::transient_const_pointer beam,double factor) const { assert(pdf_); Energy2 theScale = t * sqr(ShowerHandler::currentHandler()->factorizationScaleFactor()*factor); if (theScale < sqr(freeze_)) theScale = sqr(freeze_); const double newpdf=pdf_->xfx(beam,parton0,theScale,x/z()); if(newpdf<=0.) return 0.; const double oldpdf=pdf_->xfx(beam,parton1,theScale,x); if(oldpdf<=0.) return 1.; const double ratio = newpdf/oldpdf; double maxpdf = pdfmax_; switch (pdffactor_) { case 0: break; case 1: maxpdf /= z(); break; case 2: maxpdf /= 1.-z(); break; case 3: maxpdf /= (z()*(1.-z())); break; case 4: maxpdf /= sqrt(z()); break; case 5: maxpdf *= sqrt(z()); break; default : throw Exception() << "SudakovFormFactor::PDFVetoRatio invalid PDFfactor = " << pdffactor_ << Exception::runerror; } if (ratio > maxpdf) { generator()->log() << "PDFVeto warning: Ratio > " << name() << ":PDFmax (by a factor of " << ratio/maxpdf <<") for " << parton0->PDGName() << " to " << parton1->PDGName() << "\n"; } return ratio/maxpdf ; } void SudakovFormFactor::addSplitting(const IdList & in) { bool add=true; for(unsigned int ix=0;ix::iterator it=particles_.begin(); it!=particles_.end();++it) { if(it->size()==in.size()) { bool match=true; for(unsigned int iy=0;iy::iterator itemp=it; --itemp; particles_.erase(it); it = itemp; } } } } void SudakovFormFactor::guesstz(Energy2 t1,unsigned int iopt, const IdList &ids, double enhance,bool ident, double detune, Energy2 &t_main, double &z_main){ unsigned int pdfopt = iopt!=1 ? 0 : pdffactor_; double lower = splittingFn_->integOverP(zlimits_.first ,ids,pdfopt); double upper = splittingFn_->integOverP(zlimits_.second,ids,pdfopt); double c = 1./((upper - lower) * alpha_->overestimateValue()/Constants::twopi*enhance*detune); double r = UseRandom::rnd(); assert(iopt<=2); if(iopt==1) { c/=pdfmax_; //symmetry of FS gluon splitting if(ident) c*=0.5; } else if(iopt==2) c*=-1.; // guessing t if(iopt!=2 || c*log(r)invIntegOverP(lower + UseRandom::rnd() *(upper - lower),ids,pdfopt); } bool SudakovFormFactor::guessTimeLike(Energy2 &t,Energy2 tmin,double enhance, double detune) { Energy2 told = t; // calculate limits on z and if lower>upper return if(!computeTimeLikeLimits(t)) return false; // guess values of t and z guesstz(told,0,ids_,enhance,ids_[1]==ids_[2],detune,t,z_); // actual values for z-limits if(!computeTimeLikeLimits(t)) return false; if(tupper return if(!computeSpaceLikeLimits(t,x)) return false; // guess values of t and z guesstz(told,1,ids_,enhance,ids_[1]==ids_[2],detune,t,z_); // actual values for z-limits if(!computeSpaceLikeLimits(t,x)) return false; if(t zlimits_.second) return true; Energy2 m02 = (ids_[0]->id()!=ParticleID::g && ids_[0]->id()!=ParticleID::gamma) ? masssquared_[0] : Energy2(); Energy2 pt2 = QTildeKinematics::pT2_FSR(t,z(),m02,masssquared_[1],masssquared_[2], masssquared_[1],masssquared_[2]); // if pt2<0 veto if (pt2pT2min()) return true; // otherwise calculate pt and return pT_ = sqrt(pt2); return false; } ShoKinPtr SudakovFormFactor::generateNextTimeBranching(const Energy startingScale, const IdList &ids, const RhoDMatrix & rho, double enhance, double detuning) { // First reset the internal kinematics variables that can // have been eventually set in the previous call to the method. q_ = ZERO; z_ = 0.; phi_ = 0.; // perform initialization Energy2 tmax(sqr(startingScale)),tmin; initialize(ids,tmin); // check max > min if(tmax<=tmin) return ShoKinPtr(); // calculate next value of t using veto algorithm Energy2 t(tmax); // no shower variations to calculate - if(ShowerHandler::currentHandler()->showerVariations().empty()){ + if(!ShowerHandler::currentHandlerIsSet() || + ShowerHandler::currentHandler()->showerVariations().empty()) { // Without variations do the usual Veto algorithm // No need for more if-statements in this loop. do { if(!guessTimeLike(t,tmin,enhance,detuning)) break; } while(PSVeto(t) || SplittingFnVeto(z()*(1.-z())*t,ids,true,rho,detuning) || alphaSVeto(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t)); } else { bool alphaRew(true),PSRew(true),SplitRew(true); do { if(!guessTimeLike(t,tmin,enhance,detuning)) break; PSRew=PSVeto(t); if (PSRew) continue; SplitRew=SplittingFnVeto(z()*(1.-z())*t,ids,true,rho,detuning); alphaRew=alphaSVeto(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t); double factor=alphaSVetoRatio(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t,1.)* SplittingFnVetoRatio(z()*(1.-z())*t,ids,true,rho,detuning); tShowerHandlerPtr ch = ShowerHandler::currentHandler(); if( !(SplitRew || alphaRew) ) { //Emission q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV; if (q_ <= ZERO) break; } for ( map::const_iterator var = ch->showerVariations().begin(); var != ch->showerVariations().end(); ++var ) { if ( ( ch->firstInteraction() && var->second.firstInteraction ) || ( !ch->firstInteraction() && var->second.secondaryInteractions ) ) { double newfactor = alphaSVetoRatio(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t,var->second.renormalizationScaleFactor) * SplittingFnVetoRatio(z()*(1.-z())*t,ids,true,rho,detuning); double varied; if ( SplitRew || alphaRew ) { // No Emission varied = (1. - newfactor) / (1. - factor); } else { // Emission varied = newfactor / factor; } map::iterator wi = ch->currentWeights().find(var->first); if ( wi != ch->currentWeights().end() ) wi->second *= varied; else { assert(false); //ch->currentWeights()[var->first] = varied; } } } } while(PSRew || SplitRew || alphaRew); } q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV; if(q_ < ZERO) return ShoKinPtr(); // return the ShowerKinematics object return new_ptr(FS_QTildeShowerKinematics1to2(q_,z(),phi(),pT(),this)); } ShoKinPtr SudakovFormFactor:: generateNextSpaceBranching(const Energy startingQ, const IdList &ids, double x, const RhoDMatrix & rho, double enhance, Ptr::transient_const_pointer beam, double detuning) { // First reset the internal kinematics variables that can // have been eventually set in the previous call to the method. q_ = ZERO; z_ = 0.; phi_ = 0.; // perform the initialization Energy2 tmax(sqr(startingQ)),tmin; initialize(ids,tmin); // check max > min if(tmax<=tmin) return ShoKinPtr(); // calculate next value of t using veto algorithm Energy2 t(tmax),pt2(ZERO); // no shower variations if(ShowerHandler::currentHandler()->showerVariations().empty()){ // Without variations do the usual Veto algorithm // No need for more if-statements in this loop. do { if(!guessSpaceLike(t,tmin,x,enhance,detuning)) break; pt2 = QTildeKinematics::pT2_ISR(t,z(),masssquared_[2]); } while(pt2 < cutoff_->pT2min()|| z() > zlimits_.second|| SplittingFnVeto((1.-z())*t/z(),ids,false,rho,detuning)|| alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t)|| PDFVeto(t,x,ids[0],ids[1],beam)); } // shower variations else { bool alphaRew(true),PDFRew(true),ptRew(true),zRew(true),SplitRew(true); do { if(!guessSpaceLike(t,tmin,x,enhance,detuning)) break; pt2 = QTildeKinematics::pT2_ISR(t,z(),masssquared_[2]); ptRew=pt2 < cutoff_->pT2min(); zRew=z() > zlimits_.second; if (ptRew||zRew) continue; SplitRew=SplittingFnVeto((1.-z())*t/z(),ids,false,rho,detuning); alphaRew=alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t); PDFRew=PDFVeto(t,x,ids[0],ids[1],beam); double factor=PDFVetoRatio(t,x,ids[0],ids[1],beam,1.)* alphaSVetoRatio(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t,1.)* SplittingFnVetoRatio((1.-z())*t/z(),ids,false,rho,detuning); tShowerHandlerPtr ch = ShowerHandler::currentHandler(); if( !(PDFRew || SplitRew || alphaRew) ) { //Emission q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV; if (q_ <= ZERO) break; } for ( map::const_iterator var = ch->showerVariations().begin(); var != ch->showerVariations().end(); ++var ) { if ( ( ch->firstInteraction() && var->second.firstInteraction ) || ( !ch->firstInteraction() && var->second.secondaryInteractions ) ) { double newfactor = PDFVetoRatio(t,x,ids[0],ids[1],beam,var->second.factorizationScaleFactor)* alphaSVetoRatio(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t,var->second.renormalizationScaleFactor) *SplittingFnVetoRatio((1.-z())*t/z(),ids,false,rho,detuning); double varied; if( PDFRew || SplitRew || alphaRew) { // No Emission varied = (1. - newfactor) / (1. - factor); } else { // Emission varied = newfactor / factor; } map::iterator wi = ch->currentWeights().find(var->first); if ( wi != ch->currentWeights().end() ) wi->second *= varied; else { assert(false); //ch->currentWeights()[var->first] = varied; } } } } while( PDFRew || SplitRew || alphaRew); } if(t > ZERO && zlimits_.first < zlimits_.second) q_ = sqrt(t); else return ShoKinPtr(); pT_ = sqrt(pt2); // create the ShowerKinematics and return it return new_ptr(IS_QTildeShowerKinematics1to2(q_,z(),phi(),pT(),this)); } void SudakovFormFactor::initialize(const IdList & ids, Energy2 & tmin) { ids_=ids; tmin = 4.*cutoff_->pT2min(); masses_ = cutoff_->virtualMasses(ids); masssquared_.clear(); for(unsigned int ix=0;ix0) tmin=max(masssquared_[ix],tmin); } } ShoKinPtr SudakovFormFactor::generateNextDecayBranching(const Energy startingScale, const Energy stoppingScale, const Energy minmass, const IdList &ids, const RhoDMatrix & rho, double enhance, double detuning) { // First reset the internal kinematics variables that can // have been eventually set in the previous call to this method. q_ = Constants::MaxEnergy; z_ = 0.; phi_ = 0.; // perform initialisation Energy2 tmax(sqr(stoppingScale)),tmin; initialize(ids,tmin); tmin=sqr(startingScale); // check some branching possible if(tmax<=tmin) return ShoKinPtr(); // perform the evolution Energy2 t(tmin),pt2(-MeV2); do { if(!guessDecay(t,tmax,minmass,enhance,detuning)) break; pt2 = QTildeKinematics::pT2_Decay(t,z(),masssquared_[0],masssquared_[2]); } while(SplittingFnVeto((1.-z())*t/z(),ids,true,rho,detuning)|| alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t ) || pt2pT2min() || t*(1.-z())>masssquared_[0]-sqr(minmass)); if(t > ZERO) { q_ = sqrt(t); pT_ = sqrt(pt2); } else return ShoKinPtr(); phi_ = 0.; // create the ShowerKinematics object return new_ptr(Decay_QTildeShowerKinematics1to2(q_,z(),phi(),pT(),this)); } bool SudakovFormFactor::guessDecay(Energy2 &t,Energy2 tmax, Energy minmass, double enhance, double detune) { minmass = max(minmass,GeV); // previous scale Energy2 told = t; // overestimated limits on z if(tmaxpT2min()+ 0.25*sqr(masssquared_[2])/tm2)/tm +0.5*masssquared_[2]/tm2); if(zlimits_.secondpT2min()+ 0.25*sqr(masssquared_[2])/tm2)/tm +0.5*masssquared_[2]/tm2); if(t>tmax||zlimits_.secondpT2min(); // special case for gluon radiating if(ids_[0]->id()==ParticleID::g||ids_[0]->id()==ParticleID::gamma) { // no emission possible if(t<16.*(masssquared_[1]+pT2min)) { t=-1.*GeV2; return false; } // overestimate of the limits zlimits_.first = 0.5*(1.-sqrt(1.-4.*sqrt((masssquared_[1]+pT2min)/t))); zlimits_.second = 1.-zlimits_.first; } // special case for radiated particle is gluon else if(ids_[2]->id()==ParticleID::g||ids_[2]->id()==ParticleID::gamma) { zlimits_.first = sqrt((masssquared_[1]+pT2min)/t); zlimits_.second = 1.-sqrt((masssquared_[2]+pT2min)/t); } else if(ids_[1]->id()==ParticleID::g||ids_[1]->id()==ParticleID::gamma) { zlimits_.second = sqrt((masssquared_[2]+pT2min)/t); zlimits_.first = 1.-sqrt((masssquared_[1]+pT2min)/t); } else { zlimits_.first = (masssquared_[1]+pT2min)/t; zlimits_.second = 1.-(masssquared_[2]+pT2min)/t; } if(zlimits_.first>=zlimits_.second) { t=-1.*GeV2; return false; } return true; } bool SudakovFormFactor::computeSpaceLikeLimits(Energy2 & t, double x) { if (t < 1e-20 * GeV2) { t=-1.*GeV2; return false; } // compute the limits zlimits_.first = x; double yy = 1.+0.5*masssquared_[2]/t; zlimits_.second = yy - sqrt(sqr(yy)-1.+cutoff_->pT2min()/t); // return false if lower>upper if(zlimits_.second(particle.parents()[0]) : tShowerParticlePtr(); } else { mother = particle.children().size()==2 ? dynamic_ptr_cast(&particle) : tShowerParticlePtr(); } tShowerParticlePtr partner; while(mother) { tPPtr otherChild; if(forward) { for (unsigned int ix=0;ixchildren().size();++ix) { if(mother->children()[ix]!=child) { otherChild = mother->children()[ix]; break; } } } else { otherChild = mother->children()[1]; } tShowerParticlePtr other = dynamic_ptr_cast(otherChild); if((inter==ShowerInteraction::QCD && otherChild->dataPtr()->coloured()) || (inter==ShowerInteraction::QED && otherChild->dataPtr()->charged())) { partner = other; break; } if(forward && !other->isFinalState()) { partner = dynamic_ptr_cast(mother); break; } child = mother; if(forward) { mother = ! mother->parents().empty() ? dynamic_ptr_cast(mother->parents()[0]) : tShowerParticlePtr(); } else { if(mother->children()[0]->children().size()!=2) break; tShowerParticlePtr mtemp = dynamic_ptr_cast(mother->children()[0]); if(!mtemp) break; else mother=mtemp; } } if(!partner) { if(forward) { partner = dynamic_ptr_cast( child)->partner(); } else { if(mother) { tShowerParticlePtr parent; if(!mother->children().empty()) { parent = dynamic_ptr_cast(mother->children()[0]); } if(!parent) { parent = dynamic_ptr_cast(mother); } partner = parent->partner(); } else { partner = dynamic_ptr_cast(&particle)->partner(); } } } return partner; } pair softPhiMin(double phi0, double phi1, double A, double B, double C, double D) { double c01 = cos(phi0 - phi1); double s01 = sin(phi0 - phi1); double s012(sqr(s01)), c012(sqr(c01)); double A2(A*A), B2(B*B), C2(C*C), D2(D*D); if(abs(B/A)<1e-10 && abs(D/C)<1e-10) return make_pair(phi0,phi0+Constants::pi); double root = sqr(B2)*C2*D2*sqr(s012) + 2.*A*B2*B*C2*C*D*c01*s012 + 2.*A*B2*B*C*D2*D*c01*s012 + 4.*A2*B2*C2*D2*c012 - A2*B2*C2*D2*s012 - A2*B2*sqr(D2)*s012 - sqr(B2)*sqr(C2)*s012 - sqr(B2)*C2*D2*s012 - 4.*A2*A*B*C*D2*D*c01 - 4.*A*B2*B*C2*C*D*c01 + sqr(A2)*sqr(D2) + 2.*A2*B2*C2*D2 + sqr(B2)*sqr(C2); if(root<0.) return make_pair(phi0,phi0+Constants::pi); root = sqrt(root); double denom = (-2.*A*B*C*D*c01 + A2*D2 + B2*C2); double denom2 = (-B*C*c01 + A*D); double num = B2*C*D*s012; double y1 = B*s01*(-C*(num + root) + D*denom) / denom2; double y2 = B*s01*(-C*(num - root) + D*denom) / denom2; double x1 = -(num + root ); double x2 = -(num - root ); if(denom<0.) { y1*=-1.; y2*=-1.; x1*=-1.; x2*=-1.; } return make_pair(atan2(y1,x1) + phi0,atan2(y2,x2) + phi0); } } double SudakovFormFactor::generatePhiForward(ShowerParticle & particle, const IdList & ids, ShoKinPtr kinematics, const RhoDMatrix & rho) { // no correlations, return flat phi - if(! dynamic_ptr_cast(ShowerHandler::currentHandler())->correlations()) + if(ShowerHandler::currentHandlerIsSet() && + ! dynamic_ptr_cast(ShowerHandler::currentHandler())->correlations()) return Constants::twopi*UseRandom::rnd(); // get the kinematic variables double z = kinematics->z(); Energy2 t = z*(1.-z)*sqr(kinematics->scale()); Energy pT = kinematics->pT(); // if soft correlations Energy2 pipj,pik; bool canBeSoft[2] = {ids[1]->id()==ParticleID::g || ids[1]->id()==ParticleID::gamma, ids[2]->id()==ParticleID::g || ids[2]->id()==ParticleID::gamma }; array pjk; array Ek; Energy Ei,Ej; Energy2 m12(ZERO),m22(ZERO); InvEnergy2 aziMax(ZERO); - bool softAllowed = dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()&& + bool softAllowed = (!ShowerHandler::currentHandlerIsSet() || + dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()) && (canBeSoft[0] || canBeSoft[1]); if(softAllowed) { // find the partner for the soft correlations tShowerParticlePtr partner=findCorrelationPartner(particle,true,splittingFn()->interactionType()); // remember we want the softer gluon bool swapOrder = !canBeSoft[1] || (canBeSoft[0] && canBeSoft[1] && z < 0.5); double zFact = !swapOrder ? (1.-z) : z; // compute the transforms to the shower reference frame // first the boost Lorentz5Momentum pVect = particle.showerBasis()->pVector(); Lorentz5Momentum nVect = particle.showerBasis()->nVector(); Boost beta_bb; if(particle.showerBasis()->frame()==ShowerBasis::BackToBack) { beta_bb = -(pVect + nVect).boostVector(); } else if(particle.showerBasis()->frame()==ShowerBasis::Rest) { beta_bb = -pVect.boostVector(); } else assert(false); pVect.boost(beta_bb); nVect.boost(beta_bb); Axis axis; if(particle.showerBasis()->frame()==ShowerBasis::BackToBack) { axis = pVect.vect().unit(); } else if(particle.showerBasis()->frame()==ShowerBasis::Rest) { axis = nVect.vect().unit(); } else assert(false); // and then the rotation LorentzRotation rot; if(axis.perp2()>0.) { double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); } else if(axis.z()<0.) { rot.rotate(Constants::pi,Axis(1.,0.,0.)); } rot.invert(); pVect *= rot; nVect *= rot; // shower parameters Energy2 pn = pVect*nVect, m2 = pVect.m2(); double alpha0 = particle.showerParameters().alpha; double beta0 = 0.5/alpha0/pn* (sqr(particle.dataPtr()->mass())-sqr(alpha0)*m2+sqr(particle.showerParameters().pt)); Lorentz5Momentum qperp0(particle.showerParameters().ptx, particle.showerParameters().pty,ZERO,ZERO); assert(partner); Lorentz5Momentum pj = partner->momentum(); pj.boost(beta_bb); pj *= rot; // compute the two phi independent dot products pik = 0.5*zFact*(sqr(alpha0)*m2 - sqr(particle.showerParameters().pt) + 2.*alpha0*beta0*pn ) +0.5*sqr(pT)/zFact; Energy2 dot1 = pj*pVect; Energy2 dot2 = pj*nVect; Energy2 dot3 = pj*qperp0; pipj = alpha0*dot1+beta0*dot2+dot3; // compute the constants for the phi dependent dot product pjk[0] = zFact*(alpha0*dot1+dot3-0.5*dot2/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) +0.5*sqr(pT)*dot2/pn/zFact/alpha0; pjk[1] = (pj.x() - dot2/alpha0/pn*qperp0.x())*pT; pjk[2] = (pj.y() - dot2/alpha0/pn*qperp0.y())*pT; m12 = sqr(particle.dataPtr()->mass()); m22 = sqr(partner->dataPtr()->mass()); if(swapOrder) { pjk[1] *= -1.; pjk[2] *= -1.; } Ek[0] = zFact*(alpha0*pVect.t()-0.5*nVect.t()/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) +0.5*sqr(pT)*nVect.t()/pn/zFact/alpha0; Ek[1] = -nVect.t()/alpha0/pn*qperp0.x()*pT; Ek[2] = -nVect.t()/alpha0/pn*qperp0.y()*pT; if(swapOrder) { Ek[1] *= -1.; Ek[2] *= -1.; } Energy mag2=sqrt(sqr(Ek[1])+sqr(Ek[2])); Ei = alpha0*pVect.t()+beta0*nVect.t(); Ej = pj.t(); double phi0 = atan2(-pjk[2],-pjk[1]); if(phi0<0.) phi0 += Constants::twopi; double phi1 = atan2(-Ek[2],-Ek[1]); if(phi1<0.) phi1 += Constants::twopi; double xi_min = pik/Ei/(Ek[0]+mag2), xi_max = pik/Ei/(Ek[0]-mag2), xi_ij = pipj/Ei/Ej; if(xi_min>xi_max) swap(xi_min,xi_max); if(xi_min>xi_ij) softAllowed = false; Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2])); - if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { - aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); - } - else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { + if(!ShowerHandler::currentHandlerIsSet() || + dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { double A = (pipj*Ek[0]- Ej*pik)/Ej/sqr(Ej); double B = -sqrt(sqr(pipj)*(sqr(Ek[1])+sqr(Ek[2])))/Ej/sqr(Ej); double C = pjk[0]/sqr(Ej); double D = -sqrt(sqr(pjk[1])+sqr(pjk[2]))/sqr(Ej); pair minima = softPhiMin(phi0,phi1,A,B,C,D); aziMax = 0.5/pik/(Ek[0]-mag2)*(Ei-m12*(Ek[0]-mag2)/pik + max(Ej*(A+B*cos(minima.first -phi1))/(C+D*cos(minima.first -phi0)), Ej*(A+B*cos(minima.second-phi1))/(C+D*cos(minima.second-phi0)))); } + else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { + aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); + } else assert(false); } // if spin correlations vector > wgts; - if(dynamic_ptr_cast(ShowerHandler::currentHandler())->spinCorrelations()) { + if(!ShowerHandler::currentHandlerIsSet() || + dynamic_ptr_cast(ShowerHandler::currentHandler())->spinCorrelations()) { // calculate the weights wgts = splittingFn()->generatePhiForward(z,t,ids,rho); } else { wgts = {{ {0, 1.} }}; } // generate the azimuthal angle double phi,wgt; static const Complex ii(0.,1.); unsigned int ntry(0); double phiMax(0.),wgtMax(0.); do { phi = Constants::twopi*UseRandom::rnd(); // first the spin correlations bit (gives 1 if correlations off) Complex spinWgt = 0.; for(unsigned int ix=0;ix1e-10) { generator()->log() << "Forward spin weight problem " << wgt << " " << wgt-1. << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n"; generator()->log() << "Weights \n"; for(unsigned int ix=0;ixlog() << wgts[ix].first << " " << wgts[ix].second << "\n"; } // soft correlations bit double aziWgt = 1.; if(softAllowed) { Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi); Energy Eg = Ek[0]+Ek[1]*cos(phi)+Ek[2]*sin(phi); if(pipj*Eg>pik*Ej) { - if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { + if(!ShowerHandler::currentHandlerIsSet() || + dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { + aziWgt = max(ZERO,0.5/pik/Eg*(Ei-m12*Eg/pik + (pipj*Eg - Ej*pik)/dot)/aziMax); + } + else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; } - else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { - aziWgt = max(ZERO,0.5/pik/Eg*(Ei-m12*Eg/pik + (pipj*Eg - Ej*pik)/dot)/aziMax); - } if(aziWgt-1.>1e-10||aziWgt<-1e-10) { generator()->log() << "Forward soft weight problem " << aziWgt << " " << aziWgt-1. << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n"; } } else { aziWgt = 0.; } } wgt *= aziWgt; if(wgt>wgtMax) { phiMax = phi; wgtMax = wgt; } ++ntry; } while(wgtlog() << "Too many tries to generate phi in forward evolution\n"; phi = phiMax; } // return the azimuthal angle return phi; } double SudakovFormFactor::generatePhiBackward(ShowerParticle & particle, const IdList & ids, ShoKinPtr kinematics, const RhoDMatrix & rho) { // no correlations, return flat phi if(! dynamic_ptr_cast(ShowerHandler::currentHandler())->correlations()) return Constants::twopi*UseRandom::rnd(); // get the kinematic variables double z = kinematics->z(); Energy2 t = (1.-z)*sqr(kinematics->scale())/z; Energy pT = kinematics->pT(); // if soft correlations bool softAllowed = dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations() && (ids[2]->id()==ParticleID::g || ids[2]->id()==ParticleID::gamma); Energy2 pipj,pik,m12(ZERO),m22(ZERO); array pjk; Energy Ei,Ej,Ek; InvEnergy2 aziMax(ZERO); if(softAllowed) { // find the partner for the soft correlations tShowerParticlePtr partner=findCorrelationPartner(particle,false,splittingFn()->interactionType()); double zFact = (1.-z); // compute the transforms to the shower reference frame // first the boost Lorentz5Momentum pVect = particle.showerBasis()->pVector(); Lorentz5Momentum nVect = particle.showerBasis()->nVector(); assert(particle.showerBasis()->frame()==ShowerBasis::BackToBack); Boost beta_bb = -(pVect + nVect).boostVector(); pVect.boost(beta_bb); nVect.boost(beta_bb); Axis axis = pVect.vect().unit(); // and then the rotation LorentzRotation rot; if(axis.perp2()>0.) { double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); } else if(axis.z()<0.) { rot.rotate(Constants::pi,Axis(1.,0.,0.)); } rot.invert(); pVect *= rot; nVect *= rot; // shower parameters Energy2 pn = pVect*nVect; Energy2 m2 = pVect.m2(); double alpha0 = particle.x(); double beta0 = -0.5/alpha0/pn*sqr(alpha0)*m2; Lorentz5Momentum pj = partner->momentum(); pj.boost(beta_bb); pj *= rot; double beta2 = 0.5*(1.-zFact)*(sqr(alpha0*zFact/(1.-zFact))*m2+sqr(pT))/alpha0/zFact/pn; // compute the two phi independent dot products Energy2 dot1 = pj*pVect; Energy2 dot2 = pj*nVect; pipj = alpha0*dot1+beta0*dot2; pik = alpha0*(alpha0*zFact/(1.-zFact)*m2+pn*(beta2+zFact/(1.-zFact)*beta0)); // compute the constants for the phi dependent dot product pjk[0] = alpha0*zFact/(1.-zFact)*dot1+beta2*dot2; pjk[1] = pj.x()*pT; pjk[2] = pj.y()*pT; m12 = ZERO; m22 = sqr(partner->dataPtr()->mass()); Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2])); if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { Ek = alpha0*zFact/(1.-zFact)*pVect.t()+beta2*nVect.t(); Ei = alpha0*pVect.t()+beta0*nVect.t(); Ej = pj.t(); if(pipj*Ek> Ej*pik) { aziMax = 0.5/pik/Ek*(Ei-m12*Ek/pik + (pipj*Ek- Ej*pik)/(pjk[0]-mag)); } else { aziMax = 0.5/pik/Ek*(Ei-m12*Ek/pik); } } else { assert(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==0); } } // if spin correlations vector > wgts; if(dynamic_ptr_cast(ShowerHandler::currentHandler())->spinCorrelations()) { // get the weights wgts = splittingFn()->generatePhiBackward(z,t,ids,rho); } else { wgts = {{ {0, 1.} }}; } // generate the azimuthal angle double phi,wgt; static const Complex ii(0.,1.); unsigned int ntry(0); double phiMax(0.),wgtMax(0.); do { phi = Constants::twopi*UseRandom::rnd(); Complex spinWgt = 0.; for(unsigned int ix=0;ix1e-10) { generator()->log() << "Backward weight problem " << wgt << " " << wgt-1. << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << z << " " << phi << "\n"; generator()->log() << "Weights \n"; for(unsigned int ix=0;ixlog() << wgts[ix].first << " " << wgts[ix].second << "\n"; } // soft correlations bit double aziWgt = 1.; if(softAllowed) { Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi); if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { aziWgt = max(ZERO,0.5/pik/Ek*(Ei-m12*Ek/pik + pipj*Ek/dot - Ej*pik/dot)/aziMax); } if(aziWgt-1.>1e-10||aziWgt<-1e-10) { generator()->log() << "Backward soft weight problem " << aziWgt << " " << aziWgt-1. << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n"; } } wgt *= aziWgt; if(wgt>wgtMax) { phiMax = phi; wgtMax = wgt; } ++ntry; } while(wgtlog() << "Too many tries to generate phi in backward evolution\n"; phi = phiMax; } // return the azimuthal angle return phi; } double SudakovFormFactor::generatePhiDecay(ShowerParticle & particle, const IdList & ids, ShoKinPtr kinematics, const RhoDMatrix &) { // only soft correlations in this case // no correlations, return flat phi if( !(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations() && (ids[2]->id()==ParticleID::g || ids[2]->id()==ParticleID::gamma ))) return Constants::twopi*UseRandom::rnd(); // get the kinematic variables double z = kinematics->z(); Energy pT = kinematics->pT(); // if soft correlations // find the partner for the soft correlations tShowerParticlePtr partner = findCorrelationPartner(particle,true,splittingFn()->interactionType()); double zFact(1.-z); // compute the transforms to the shower reference frame // first the boost Lorentz5Momentum pVect = particle.showerBasis()->pVector(); Lorentz5Momentum nVect = particle.showerBasis()->nVector(); assert(particle.showerBasis()->frame()==ShowerBasis::Rest); Boost beta_bb = -pVect.boostVector(); pVect.boost(beta_bb); nVect.boost(beta_bb); Axis axis = nVect.vect().unit(); // and then the rotation LorentzRotation rot; if(axis.perp2()>0.) { double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); } else if(axis.z()<0.) { rot.rotate(Constants::pi,Axis(1.,0.,0.)); } rot.invert(); pVect *= rot; nVect *= rot; // shower parameters Energy2 pn = pVect*nVect; Energy2 m2 = pVect.m2(); double alpha0 = particle.showerParameters().alpha; double beta0 = 0.5/alpha0/pn* (sqr(particle.dataPtr()->mass())-sqr(alpha0)*m2+sqr(particle.showerParameters().pt)); Lorentz5Momentum qperp0(particle.showerParameters().ptx, particle.showerParameters().pty,ZERO,ZERO); Lorentz5Momentum pj = partner->momentum(); pj.boost(beta_bb); pj *= rot; // compute the two phi independent dot products Energy2 pik = 0.5*zFact*(sqr(alpha0)*m2 - sqr(particle.showerParameters().pt) + 2.*alpha0*beta0*pn ) +0.5*sqr(pT)/zFact; Energy2 dot1 = pj*pVect; Energy2 dot2 = pj*nVect; Energy2 dot3 = pj*qperp0; Energy2 pipj = alpha0*dot1+beta0*dot2+dot3; // compute the constants for the phi dependent dot product array pjk; pjk[0] = zFact*(alpha0*dot1+dot3-0.5*dot2/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) +0.5*sqr(pT)*dot2/pn/zFact/alpha0; pjk[1] = (pj.x() - dot2/alpha0/pn*qperp0.x())*pT; pjk[2] = (pj.y() - dot2/alpha0/pn*qperp0.y())*pT; Energy2 m12 = sqr(particle.dataPtr()->mass()); Energy2 m22 = sqr(partner->dataPtr()->mass()); Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2])); InvEnergy2 aziMax; array Ek; Energy Ei,Ej; if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { Ek[0] = zFact*(alpha0*pVect.t()+-0.5*nVect.t()/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) +0.5*sqr(pT)*nVect.t()/pn/zFact/alpha0; Ek[1] = -nVect.t()/alpha0/pn*qperp0.x()*pT; Ek[2] = -nVect.t()/alpha0/pn*qperp0.y()*pT; Energy mag2=sqrt(sqr(Ek[1])+sqr(Ek[2])); Ei = alpha0*pVect.t()+beta0*nVect.t(); Ej = pj.t(); aziMax = 0.5/pik/(Ek[0]-mag2)*(Ei-m12*(Ek[0]-mag2)/pik + pipj*(Ek[0]+mag2)/(pjk[0]-mag) - Ej*pik/(pjk[0]-mag) ); } else assert(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==0); // generate the azimuthal angle double phi,wgt(0.); unsigned int ntry(0); double phiMax(0.),wgtMax(0.); do { phi = Constants::twopi*UseRandom::rnd(); Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi); if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { wgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { if(qperp0.m2()==ZERO) { wgt = 1.; } else { Energy Eg = Ek[0]+Ek[1]*cos(phi)+Ek[2]*sin(phi); wgt = max(ZERO,0.5/pik/Eg*(Ei-m12*Eg/pik + (pipj*Eg - Ej*pik)/dot)/aziMax); } } if(wgt-1.>1e-10||wgt<-1e-10) { generator()->log() << "Decay soft weight problem " << wgt << " " << wgt-1. << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n"; } if(wgt>wgtMax) { phiMax = phi; wgtMax = wgt; } ++ntry; } while(wgtlog() << "Too many tries to generate phi\n"; } // return the azimuthal angle return phi; } Energy SudakovFormFactor::calculateScale(double zin, Energy pt, IdList ids, unsigned int iopt) { Energy2 tmin; initialize(ids,tmin); // final-state branching if(iopt==0) { Energy2 scale=(sqr(pt)+masssquared_[1]*(1.-zin)+masssquared_[2]*zin); if(ids[0]->id()!=ParticleID::g) scale -= zin*(1.-zin)*masssquared_[0]; scale /= sqr(zin*(1-zin)); return scale<=ZERO ? sqrt(tmin) : sqrt(scale); } else if(iopt==1) { Energy2 scale=(sqr(pt)+zin*masssquared_[2])/sqr(1.-zin); return scale<=ZERO ? sqrt(tmin) : sqrt(scale); } else if(iopt==2) { Energy2 scale = (sqr(pt)+zin*masssquared_[2])/sqr(1.-zin)+masssquared_[0]; return scale<=ZERO ? sqrt(tmin) : sqrt(scale); } else { throw Exception() << "Unknown option in SudakovFormFactor::calculateScale() " << "iopt = " << iopt << Exception::runerror; } }