diff --git a/Shower/QTilde/Base/KinematicsReconstructor.cc b/Shower/QTilde/Base/KinematicsReconstructor.cc --- a/Shower/QTilde/Base/KinematicsReconstructor.cc +++ b/Shower/QTilde/Base/KinematicsReconstructor.cc @@ -1,30 +1,2942 @@ // -*- C++ -*- // // KinematicsReconstructor.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the 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/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 using namespace Herwig; -DescribeAbstractNoPIOClass -describeKinematicsReconstructor ("Herwig::KinematicsReconstructor","HwShower.so"); +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 << _finalStateReconOption + << _initialStateReconOption; +} + +void KinematicsReconstructor::persistentInput(PersistentIStream & is, int) { + is >> _reconopt >> _initialBoost >> iunit(_minQ,GeV) >> _noRescale + >> _noRescaleVector >> _finalStateReconOption + >> _initialStateReconOption; +} 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 Switch interfaceFinalStateReconOption + ("FinalStateReconOption", + "Option for how to reconstruct the momenta of the final-state system", + &KinematicsReconstructor::_finalStateReconOption, 0, false, false); + static SwitchOption interfaceFinalStateReconOptionDefault + (interfaceFinalStateReconOption, + "Default", + "All the momenta are rescaled in the rest frame", + 0); + static SwitchOption interfaceFinalStateReconOptionMostOffShell + (interfaceFinalStateReconOption, + "MostOffShell", + "All particles put on the new-mass shell and then the most off-shell and" + " recoiling system are rescaled to ensure 4-momentum is conserved.", + 1); + static SwitchOption interfaceFinalStateReconOptionRecursive + (interfaceFinalStateReconOption, + "Recursive", + "Recursively put on shell by putting the most off-shell particle which" + " hasn't been rescaled on-shell by rescaling the particles and the recoiling system. ", + 2); + static SwitchOption interfaceFinalStateReconOptionRestMostOffShell + (interfaceFinalStateReconOption, + "RestMostOffShell", + "The most off-shell is put on shell by rescaling it and the recoiling system," + " the recoiling system is then put on-shell in its rest frame.", + 3); + static SwitchOption interfaceFinalStateReconOptionRestRecursive + (interfaceFinalStateReconOption, + "RestRecursive", + "As 3 but recursive treated the currently most-off shell," + " only makes a difference if more than 3 partons.", + 4); + + 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); + } + +void KinematicsReconstructor::doinit() { + Interfaced::doinit(); + _noRescale = set(_noRescaleVector.begin(),_noRescaleVector.end()); +} + +bool KinematicsReconstructor:: +reconstructTimeLikeJet(const tShowerParticlePtr particleJetParent) 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)); + } + // 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()) { + jetGrandParent->showerKinematics()->reconstructLast(particleJetParent, + _progenitor->mass()); + } + else { + jetGrandParent->showerKinematics()->reconstructLast(particleJetParent); + } + } + } + // otherwise + else { + Energy dm = particleJetParent->data().constituentMass(); + if (abs(dm-particleJetParent->momentum().m())>0.001*MeV + &&particleJetParent->dataPtr()->stable() + &&particleJetParent->id()!=ParticleID::gamma + &&_noRescale.find(particleJetParent->dataPtr())==_noRescale.end()) { + Lorentz5Momentum dum = particleJetParent->momentum(); + dum.setMass(dm); + dum.rescaleEnergy(); + 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); + // 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); + 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); + // 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); +} + +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()); + 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()); + (**cit).reconstructed(ShowerProgenitor::done); + } + else { + radiated |= !(*cit)->progenitor()->children().empty(); + } + tempJetKin.q = (*cit)->progenitor()->momentum(); + jetKinematics.push_back(tempJetKin); + } + // default option rescale everything with the same factor + if( _finalStateReconOption == 0 || jetKinematics.size() <= 2 ) { + // 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); + } + } + } + // different treatment of most off-shell + else if ( _finalStateReconOption <= 4 ) { + // sort the jets by virtuality + std::sort(jetKinematics.begin(),jetKinematics.end(),JetOrdering()); + // Bryan's procedures from FORTRAN + if( _finalStateReconOption <=2 ) { + // loop over the off-shell partons, _finalStateReconOption==1 only first ==2 all + JetKinVect::const_iterator jend = _finalStateReconOption==1 ? jetKinematics.begin()+1 : jetKinematics.end(); + for(JetKinVect::const_iterator jit=jetKinematics.begin(); jit!=jend;++jit) { + // calculate the 4-momentum of the recoiling system + Lorentz5Momentum psum; + bool done = true; + for(JetKinVect::const_iterator it=jetKinematics.begin();it!=jetKinematics.end();++it) { + if(it==jit) { + done = false; + continue; + } + // first option put on-shell and sum 4-momenta + if( _finalStateReconOption == 1 ) { + LorentzRotation Trafo = solveBoost(1., it->q, it->p); + deepTransform(it->parent,Trafo); + psum += it->parent->momentum(); + } + // second option, sum momenta + else { + // already rescaled + if(done) psum += it->parent->momentum(); + // still needs to be rescaled + else psum += it->p; + } + } + // set the mass + psum.rescaleMass(); + // calculate the 3-momentum rescaling factor + Energy2 s(pcm.m2()); + Energy2 m1sq(jit->q.m2()),m2sq(psum.m2()); + Energy4 num = sqr(s - m1sq - m2sq) - 4.*m1sq*m2sq; + if(nump.vect().mag2()) ); + // boost the off-shell parton + LorentzRotation B1 = solveBoost(k, jit->q, jit->p); + deepTransform(jit->parent,B1); + // boost everything else to rescale + LorentzRotation B2 = solveBoost(k, psum, psum); + for(JetKinVect::iterator it=jetKinematics.begin();it!=jetKinematics.end();++it) { + if(it==jit) continue; + deepTransform(it->parent,B2); + it->p *= B2; + it->q *= B2; + } + } + } + // Peter's C++ procedures + else { + reconstructFinalFinalOffShell(jetKinematics,pcm.m2(), _finalStateReconOption == 4); + } + } + else + assert(false); + // 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 + 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); + } + } +} + +void KinematicsReconstructor::reconstructFinalFinalOffShell(JetKinVect orderedJets, + Energy2 s, + bool recursive) const { + JetKinVect::iterator jit; + jit = orderedJets.begin(); ++jit; + // 4-momentum of recoiling system + Lorentz5Momentum psum; + for( ; jit!=orderedJets.end(); ++jit) psum += jit->p; + psum.rescaleMass(); + // calculate the 3-momentum rescaling factor + Energy2 m1sq(orderedJets.begin()->q.m2()),m2sq(psum.m2()); + Energy4 num = sqr(s - m1sq - m2sq) - 4.*m1sq*m2sq; + if(nump.vect().mag2()) ); + // boost the most off-shell + LorentzRotation B1 = solveBoost(k, orderedJets.begin()->q, orderedJets.begin()->p); + deepTransform(orderedJets.begin()->parent,B1); + // boost everything else + // first to rescale + LorentzRotation B2 = solveBoost(k, psum, psum); + // and then to rest frame of new system + Lorentz5Momentum pnew = B2*psum; + pnew.rescaleMass(); + B2.transform(pnew.findBoostToCM()); + // apply transform (calling routine ensures at least 3 elements) + jit = orderedJets.begin(); ++jit; + for(;jit!=orderedJets.end();++jit) { + deepTransform(jit->parent,B2); + jit->p *= B2; + jit->q *= B2; + } + JetKinVect newJets(orderedJets.begin()+1,orderedJets.end()); + // final reconstruction + if(newJets.size()==2 || !recursive ) { + // rescaling factor + double k = solveKfactor(psum.m(), newJets); + // rescale jets in the new CMF + for(JetKinVect::iterator it = newJets.begin(); it != newJets.end(); ++it) { + LorentzRotation Trafo = solveBoost(k, it->q, it->p); + deepTransform(it->parent,Trafo); + } + } + // recursive + else { + std::sort(newJets.begin(),newJets.end(),JetOrdering()); + reconstructFinalFinalOffShell(newJets,psum.m2(),recursive); + } + // finally boost back from new CMF + LorentzRotation back(-pnew.findBoostToCM()); + for(JetKinVect::iterator it = newJets.begin(); it != newJets.end(); ++it) { + deepTransform(it->parent,back); + } +} + +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), _minQ(MeV) {}; + + /** * Methods to reconstruct the kinematics of a scattering or decay process */ //@{ /** - * Given the ShowerTree for the shower from a hard process - * the method does the reconstruction of the jets, + * 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=0; + bool switchRecon) const; /** - * Given the ShowerTree for a decay shower - * the method does the reconstruction of the jets, + * 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=0; + 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 for a decay process + * shower */ - virtual bool deconstructDecayJets(HardTreePtr decay,ShowerInteraction) const=0; + 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 hard,ShowerInteraction) const=0; + 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 + */ + //@{ + /** + * 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; + + /** + * 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 &); + 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; + + /** + * 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; }; } +#include "KinematicsReconstructor.tcc" #endif /* HERWIG_KinematicsReconstructor_H */ diff --git a/Shower/QTilde/Default/QTildeReconstructor.tcc b/Shower/QTilde/Base/KinematicsReconstructor.tcc rename from Shower/QTilde/Default/QTildeReconstructor.tcc rename to Shower/QTilde/Base/KinematicsReconstructor.tcc --- a/Shower/QTilde/Default/QTildeReconstructor.tcc +++ b/Shower/QTilde/Base/KinematicsReconstructor.tcc @@ -1,247 +1,247 @@ // -*- C++ -*- // -// QTildeReconstructor.tcc is a part of Herwig - A multi-purpose Monte Carlo event generator +// KinematicsReconstructor.tcc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined templated member -// functions of the QTildeReconstructor class. +// functions of the KinematicsReconstructor class. // using namespace Herwig; using namespace ThePEG; namespace { /** * find showering particle for hard branchings */ tShowerParticlePtr SHOWERINGPARTICLE(HardBranchingPtr a) { return a->branchingParticle(); } /** * find showering particle for progenitors */ tShowerParticlePtr SHOWERINGPARTICLE(ShowerProgenitorPtr a) { return a->progenitor(); } /** * Return colour line progenitor pointer for ShowerProgenitor */ template Ptr::transient_pointer CL(Value a, unsigned int index=0) { return const_ptr_cast(SHOWERINGPARTICLE(a)->colourInfo()->colourLines()[index]); } /** * Return progenitor colour line size for ShowerProgenitor */ template unsigned int CLSIZE(Value a) { return SHOWERINGPARTICLE(a)->colourInfo()->colourLines().size(); } /** * Return anti-colour line progenitor pointer for ShowerProgenitor */ template Ptr::transient_pointer ACL(Value a, unsigned int index=0) { return const_ptr_cast(SHOWERINGPARTICLE(a)->colourInfo()->antiColourLines()[index]); } /** * Return progenitor anti-colour line size for ShowerProgenitor */ template unsigned int ACLSIZE(Value a) { return SHOWERINGPARTICLE(a)->colourInfo()->antiColourLines().size(); } } -template void QTildeReconstructor:: +template void KinematicsReconstructor:: findPartners(Value jet,set & done, const set & jets, vector & system) const { tShowerParticlePtr part=SHOWERINGPARTICLE(jet); unsigned int partNumColourLines = part->colourInfo()-> colourLines().size(); unsigned int partNumAColourLines = part->colourInfo()->antiColourLines().size(); for(typename set::const_iterator cit=jets.begin();cit!=jets.end();++cit) { if(done.find(*cit)!=done.end()||!SHOWERINGPARTICLE(*cit)->coloured()) continue; bool isPartner = false; // one initial one final if(part->isFinalState()!=SHOWERINGPARTICLE(*cit)->isFinalState()) { //loop over all the colours of both for(unsigned int ix=0; ixcolourLine()) { for(unsigned int ix=0; ixantiColourLine()&&!isPartner) { for(unsigned int ix=0; ixcolourLine()) { if(part->colourLine()->sourceNeighbours().first) { tColinePair lines = part->colourLine()->sourceNeighbours(); if(lines.first == CL(*cit) || lines.first == ACL(*cit) || lines.second == CL(*cit) || lines.second == ACL(*cit) ) isPartner = true; } if(part->colourLine()->sinkNeighbours().first) { tColinePair lines = part->colourLine()->sinkNeighbours(); if(lines.first == CL(*cit) || lines.first == ACL(*cit) || lines.second == CL(*cit) || lines.second == ACL(*cit) ) isPartner = true; } } if(part->antiColourLine()) { if(part->antiColourLine()->sourceNeighbours().first) { tColinePair lines = part->antiColourLine()->sourceNeighbours(); if(lines.first == CL(*cit) || lines.first == ACL(*cit) || lines.second == CL(*cit) || lines.second == ACL(*cit) ) isPartner = true; } if(part->antiColourLine()->sinkNeighbours().first) { tColinePair lines = part->antiColourLine()->sinkNeighbours(); if(lines.first == CL(*cit) || lines.first == ACL(*cit) || lines.second == CL(*cit) || lines.second == ACL(*cit) ) isPartner = true; } } if(isPartner) { system.push_back(*cit); done.insert(*cit); findPartners(*cit,done,jets,system); } } } template -typename Herwig::ColourSinglet::VecType QTildeReconstructor:: +typename Herwig::ColourSinglet::VecType KinematicsReconstructor:: identifySystems(set jets, unsigned int & nnun,unsigned int & nnii,unsigned int & nnif, unsigned int & nnf ,unsigned int & nni ) const { vector > systems; set done; for(typename set::const_iterator it=jets.begin();it!=jets.end();++it) { // if not treated create new system if(done.find(*it)!=done.end()) continue; done.insert(*it); systems.push_back(ColourSinglet (UNDEFINED,*it)); if(!SHOWERINGPARTICLE(*it)->coloured()) continue; findPartners(*it,done,jets,systems.back().jets); } for(unsigned int ix=0;ixisFinalState()) ++nf; else ++ni; } // type // initial-initial if(ni==2&&nf==0) { systems[ix].type = II; ++nnii; } // initial only else if(ni==1&&nf==0) { systems[ix].type = I; ++nni; } // initial-final else if(ni==1&&nf>0) { systems[ix].type = IF; ++nnif; } // final only else if(ni==0&&nf>0) { systems[ix].type = F; ++nnf; } // otherwise unknown else { systems[ix].type = UNDEFINED; ++nnun; } } return systems; } template -void QTildeReconstructor::combineFinalState(vector > & systems) const { +void KinematicsReconstructor::combineFinalState(vector > & systems) const { // check that 1 particle final-state systems which can be combine bool canCombine(true); for(unsigned int ix=0;ix > oldsystems=systems; systems.clear(); ColourSinglet finalState; finalState.type = F; for(unsigned int ix=0;ix - -using namespace Herwig; - -DescribeClass -describeQTildeReconstructor("Herwig::QTildeReconstructor", "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 QTildeReconstructor::persistentOutput(PersistentOStream & os) const { - os << _reconopt << _initialBoost << ounit(_minQ,GeV) << _noRescale - << _noRescaleVector << _finalStateReconOption - << _initialStateReconOption; -} - -void QTildeReconstructor::persistentInput(PersistentIStream & is, int) { - is >> _reconopt >> _initialBoost >> iunit(_minQ,GeV) >> _noRescale - >> _noRescaleVector >> _finalStateReconOption - >> _initialStateReconOption; -} - -void QTildeReconstructor::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", - &QTildeReconstructor::_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", - &QTildeReconstructor::_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", - &QTildeReconstructor::_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.", - &QTildeReconstructor::_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 Switch interfaceFinalStateReconOption - ("FinalStateReconOption", - "Option for how to reconstruct the momenta of the final-state system", - &QTildeReconstructor::_finalStateReconOption, 0, false, false); - static SwitchOption interfaceFinalStateReconOptionDefault - (interfaceFinalStateReconOption, - "Default", - "All the momenta are rescaled in the rest frame", - 0); - static SwitchOption interfaceFinalStateReconOptionMostOffShell - (interfaceFinalStateReconOption, - "MostOffShell", - "All particles put on the new-mass shell and then the most off-shell and" - " recoiling system are rescaled to ensure 4-momentum is conserved.", - 1); - static SwitchOption interfaceFinalStateReconOptionRecursive - (interfaceFinalStateReconOption, - "Recursive", - "Recursively put on shell by putting the most off-shell particle which" - " hasn't been rescaled on-shell by rescaling the particles and the recoiling system. ", - 2); - static SwitchOption interfaceFinalStateReconOptionRestMostOffShell - (interfaceFinalStateReconOption, - "RestMostOffShell", - "The most off-shell is put on shell by rescaling it and the recoiling system," - " the recoiling system is then put on-shell in its rest frame.", - 3); - static SwitchOption interfaceFinalStateReconOptionRestRecursive - (interfaceFinalStateReconOption, - "RestRecursive", - "As 3 but recursive treated the currently most-off shell," - " only makes a difference if more than 3 partons.", - 4); - - static Switch interfaceInitialStateReconOption - ("InitialStateReconOption", - "Option for the reconstruction of initial state radiation", - &QTildeReconstructor::_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); - -} - -void QTildeReconstructor::doinit() { - KinematicsReconstructor::doinit(); - _noRescale = set(_noRescaleVector.begin(),_noRescaleVector.end()); -} - -bool QTildeReconstructor:: -reconstructTimeLikeJet(const tShowerParticlePtr particleJetParent) 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)); - } - // 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()) { - jetGrandParent->showerKinematics()->reconstructLast(particleJetParent, - _progenitor->mass()); - } - else { - jetGrandParent->showerKinematics()->reconstructLast(particleJetParent); - } - } - } - // otherwise - else { - Energy dm = particleJetParent->data().constituentMass(); - if (abs(dm-particleJetParent->momentum().m())>0.001*MeV - &&particleJetParent->dataPtr()->stable() - &&particleJetParent->id()!=ParticleID::gamma - &&_noRescale.find(particleJetParent->dataPtr())==_noRescale.end()) { - Lorentz5Momentum dum = particleJetParent->momentum(); - dum.setMass(dm); - dum.rescaleEnergy(); - 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 QTildeReconstructor:: -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 -QTildeReconstructor::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 QTildeReconstructor:: -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); - // calculate the momentum of the particle - Lorentz5Momentum pnew=p->momentum()-child->momentum(); - pnew.rescaleMass(); - p->children()[0]->set5Momentum(pnew); - } - } - return emitted; -} - -Boost QTildeReconstructor:: -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 QTildeReconstructor:: -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); - 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 QTildeReconstructor:: -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); - // 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 QTildeReconstructor:: -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 QTildeReconstructor:: -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 QTildeReconstructor:: -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 QTildeReconstructor::" - << "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 QTildeReconstructor::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 QTildeReconstructor::deconstructHardJets(). Please use one of the other options\n" - << Exception::runerror; - } - else if(_reconopt == 3 || _reconopt == 4 ) { - return deconstructColourPartner(tree,type); - } - else - assert(false); -} - -bool QTildeReconstructor:: -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 QTildeReconstructor:: -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 " - << "QTildeReconstructor::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 " - << "QTildeReconstructor::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 " - << "QTildeReconstructor::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 QTildeReconstructor:: -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()); - 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 QTildeReconstructor::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 QTildeReconstructor:: -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 QTildeReconstructor::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 QTildeReconstructor::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 QTildeReconstructor:: -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()); - (**cit).reconstructed(ShowerProgenitor::done); - } - else { - radiated |= !(*cit)->progenitor()->children().empty(); - } - tempJetKin.q = (*cit)->progenitor()->momentum(); - jetKinematics.push_back(tempJetKin); - } - // default option rescale everything with the same factor - if( _finalStateReconOption == 0 || jetKinematics.size() <= 2 ) { - // 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); - } - } - } - // different treatment of most off-shell - else if ( _finalStateReconOption <= 4 ) { - // sort the jets by virtuality - std::sort(jetKinematics.begin(),jetKinematics.end(),JetOrdering()); - // Bryan's procedures from FORTRAN - if( _finalStateReconOption <=2 ) { - // loop over the off-shell partons, _finalStateReconOption==1 only first ==2 all - JetKinVect::const_iterator jend = _finalStateReconOption==1 ? jetKinematics.begin()+1 : jetKinematics.end(); - for(JetKinVect::const_iterator jit=jetKinematics.begin(); jit!=jend;++jit) { - // calculate the 4-momentum of the recoiling system - Lorentz5Momentum psum; - bool done = true; - for(JetKinVect::const_iterator it=jetKinematics.begin();it!=jetKinematics.end();++it) { - if(it==jit) { - done = false; - continue; - } - // first option put on-shell and sum 4-momenta - if( _finalStateReconOption == 1 ) { - LorentzRotation Trafo = solveBoost(1., it->q, it->p); - deepTransform(it->parent,Trafo); - psum += it->parent->momentum(); - } - // second option, sum momenta - else { - // already rescaled - if(done) psum += it->parent->momentum(); - // still needs to be rescaled - else psum += it->p; - } - } - // set the mass - psum.rescaleMass(); - // calculate the 3-momentum rescaling factor - Energy2 s(pcm.m2()); - Energy2 m1sq(jit->q.m2()),m2sq(psum.m2()); - Energy4 num = sqr(s - m1sq - m2sq) - 4.*m1sq*m2sq; - if(nump.vect().mag2()) ); - // boost the off-shell parton - LorentzRotation B1 = solveBoost(k, jit->q, jit->p); - deepTransform(jit->parent,B1); - // boost everything else to rescale - LorentzRotation B2 = solveBoost(k, psum, psum); - for(JetKinVect::iterator it=jetKinematics.begin();it!=jetKinematics.end();++it) { - if(it==jit) continue; - deepTransform(it->parent,B2); - it->p *= B2; - it->q *= B2; - } - } - } - // Peter's C++ procedures - else { - reconstructFinalFinalOffShell(jetKinematics,pcm.m2(), _finalStateReconOption == 4); - } - } - else - assert(false); - // 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 QTildeReconstructor:: -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 QTildeReconstructor:: -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 QTildeReconstructor::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 QTildeReconstructor::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 QTildeReconstructor:: -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 QTildeReconstructor:: -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 QTildeReconstructor:: -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 QTildeReconstructor:: -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 QTildeReconstructor::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 - 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); - } - } -} - -void QTildeReconstructor::reconstructFinalFinalOffShell(JetKinVect orderedJets, - Energy2 s, - bool recursive) const { - JetKinVect::iterator jit; - jit = orderedJets.begin(); ++jit; - // 4-momentum of recoiling system - Lorentz5Momentum psum; - for( ; jit!=orderedJets.end(); ++jit) psum += jit->p; - psum.rescaleMass(); - // calculate the 3-momentum rescaling factor - Energy2 m1sq(orderedJets.begin()->q.m2()),m2sq(psum.m2()); - Energy4 num = sqr(s - m1sq - m2sq) - 4.*m1sq*m2sq; - if(nump.vect().mag2()) ); - // boost the most off-shell - LorentzRotation B1 = solveBoost(k, orderedJets.begin()->q, orderedJets.begin()->p); - deepTransform(orderedJets.begin()->parent,B1); - // boost everything else - // first to rescale - LorentzRotation B2 = solveBoost(k, psum, psum); - // and then to rest frame of new system - Lorentz5Momentum pnew = B2*psum; - pnew.rescaleMass(); - B2.transform(pnew.findBoostToCM()); - // apply transform (calling routine ensures at least 3 elements) - jit = orderedJets.begin(); ++jit; - for(;jit!=orderedJets.end();++jit) { - deepTransform(jit->parent,B2); - jit->p *= B2; - jit->q *= B2; - } - JetKinVect newJets(orderedJets.begin()+1,orderedJets.end()); - // final reconstruction - if(newJets.size()==2 || !recursive ) { - // rescaling factor - double k = solveKfactor(psum.m(), newJets); - // rescale jets in the new CMF - for(JetKinVect::iterator it = newJets.begin(); it != newJets.end(); ++it) { - LorentzRotation Trafo = solveBoost(k, it->q, it->p); - deepTransform(it->parent,Trafo); - } - } - // recursive - else { - std::sort(newJets.begin(),newJets.end(),JetOrdering()); - reconstructFinalFinalOffShell(newJets,psum.m2(),recursive); - } - // finally boost back from new CMF - LorentzRotation back(-pnew.findBoostToCM()); - for(JetKinVect::iterator it = newJets.begin(); it != newJets.end(); ++it) { - deepTransform(it->parent,back); - } -} - -Energy QTildeReconstructor::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 -QTildeReconstructor::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 -QTildeReconstructor::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 QTildeReconstructor:: -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 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 QTildeReconstructorInterfaces "The interfaces" - * defined for QTildeReconstructor. - */ -class QTildeReconstructor: public KinematicsReconstructor { - -public: - - /** - * Default constructor - */ - QTildeReconstructor() : _reconopt(0), _initialBoost(0), - _finalStateReconOption(0), - _initialStateReconOption(0), _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 - */ - //@{ - /** - * 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; - - /** - * 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. - */ - QTildeReconstructor & operator=(const QTildeReconstructor &); - -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; - - /** - * 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; -}; - -} - -#include "QTildeReconstructor.tcc" -#endif /* HERWIG_QTildeReconstructor_H */ diff --git a/Shower/QTilde/Makefile.am b/Shower/QTilde/Makefile.am --- a/Shower/QTilde/Makefile.am +++ b/Shower/QTilde/Makefile.am @@ -1,42 +1,42 @@ SUBDIRS = Matching pkglib_LTLIBRARIES = HwShower.la HwShower_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 25:0:0 HwShower_la_SOURCES = \ Couplings/ShowerAlphaQCD.h Couplings/ShowerAlphaQCD.cc \ Couplings/ShowerAlphaQED.h Couplings/ShowerAlphaQED.cc\ QTildeShowerHandler.h QTildeShowerHandler.fh QTildeShowerHandler.cc \ SplittingFunctions/HalfHalfOneSplitFn.h SplittingFunctions/HalfHalfOneSplitFn.cc\ SplittingFunctions/OneOneOneSplitFn.h SplittingFunctions/OneOneOneSplitFn.cc\ SplittingFunctions/OneOneOneMassiveSplitFn.h SplittingFunctions/OneOneOneMassiveSplitFn.cc\ SplittingFunctions/ZeroZeroOneSplitFn.h SplittingFunctions/ZeroZeroOneSplitFn.cc\ SplittingFunctions/OneHalfHalfSplitFn.h SplittingFunctions/OneHalfHalfSplitFn.cc\ SplittingFunctions/HalfOneHalfSplitFn.h SplittingFunctions/HalfOneHalfSplitFn.cc\ SplittingFunctions/CMWOneOneOneSplitFn.h SplittingFunctions/CMWOneOneOneSplitFn.cc\ SplittingFunctions/CMWHalfHalfOneSplitFn.h SplittingFunctions/CMWHalfHalfOneSplitFn.cc\ Default/QTildeSudakov.cc Default/QTildeSudakov.h\ Default/Decay_QTildeShowerKinematics1to2.cc \ Default/Decay_QTildeShowerKinematics1to2.h \ Default/IS_QTildeShowerKinematics1to2.cc Default/IS_QTildeShowerKinematics1to2.h \ Default/FS_QTildeShowerKinematics1to2.cc Default/FS_QTildeShowerKinematics1to2.h \ -Default/QTildeReconstructor.cc Default/QTildeReconstructor.h Default/QTildeReconstructor.tcc \ Base/KinematicsReconstructor.cc \ +Base/KinematicsReconstructor.tcc \ Base/KinematicsReconstructor.h \ Base/KinematicsReconstructor.fh \ Base/HardTree.cc Base/HardTree.h Base/HardTree.fh \ Base/HardBranching.h Base/HardBranching.fh Base/HardBranching.cc\ Base/PartnerFinder.h Base/PartnerFinder.fh Base/PartnerFinder.cc \ Base/ShowerVeto.h Base/ShowerVeto.fh Base/ShowerVeto.cc \ Base/FullShowerVeto.h Base/FullShowerVeto.fh Base/FullShowerVeto.cc \ SplittingFunctions/SplittingGenerator.cc SplittingFunctions/SplittingGenerator.h\ SplittingFunctions/SplittingGenerator.fh \ Base/ShowerTree.h Base/ShowerTree.fh Base/ShowerTree.cc \ ShowerConfig.h ShowerConfig.cc \ Base/Branching.h \ Base/ShowerParticle.cc Base/ShowerParticle.fh Base/ShowerParticle.h \ Base/ShowerKinematics.fh Base/ShowerKinematics.h Base/ShowerKinematics.cc \ Base/ShowerBasis.fh Base/ShowerBasis.h Base/ShowerBasis.cc \ Base/ShowerProgenitor.fh Base/ShowerProgenitor.h \ Base/SudakovFormFactor.cc Base/SudakovFormFactor.h Base/SudakovFormFactor.fh \ SplittingFunctions/SplittingFunction.h SplittingFunctions/SplittingFunction.fh \ SplittingFunctions/SplittingFunction.cc \ Base/ShowerVertex.cc Base/ShowerVertex.fh Base/ShowerVertex.h diff --git a/src/defaults/Shower.in b/src/defaults/Shower.in --- a/src/defaults/Shower.in +++ b/src/defaults/Shower.in @@ -1,310 +1,310 @@ # -*- ThePEG-repository -*- ############################################################ # Setup of default parton shower # # Useful switches for users are marked near the top of # this file. # # Don't edit this file directly, but reset the switches # in your own input files! ############################################################ library HwMPI.so library HwShower.so library HwMatching.so mkdir /Herwig/Shower cd /Herwig/Shower create Herwig::QTildeShowerHandler ShowerHandler newdef ShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler newdef ShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer # use LO PDFs for Shower, can be changed later newdef ShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF newdef ShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF newdef ShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF newdef ShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF ##################################### # initial setup, don't change these! ##################################### create Herwig::SplittingGenerator SplittingGenerator create Herwig::ShowerAlphaQCD AlphaQCD create Herwig::ShowerAlphaQED AlphaQED create Herwig::PartnerFinder PartnerFinder newdef PartnerFinder:PartnerMethod 1 newdef PartnerFinder:ScaleChoice 1 -create Herwig::QTildeReconstructor KinematicsReconstructor +create Herwig::KinematicsReconstructor KinematicsReconstructor newdef KinematicsReconstructor:ReconstructionOption Colour3 newdef KinematicsReconstructor:InitialStateReconOption SofterFraction newdef KinematicsReconstructor:InitialInitialBoostOption LongTransBoost newdef /Herwig/Partons/RemnantDecayer:AlphaS AlphaQCD newdef /Herwig/Partons/RemnantDecayer:AlphaEM AlphaQED newdef ShowerHandler:PartnerFinder PartnerFinder newdef ShowerHandler:KinematicsReconstructor KinematicsReconstructor newdef ShowerHandler:SplittingGenerator SplittingGenerator newdef ShowerHandler:SpinCorrelations Yes newdef ShowerHandler:SoftCorrelations Singular ################################################################## # Intrinsic pT # # Recommended: # 1.9 GeV for Tevatron W/Z production. # 2.1 GeV for LHC W/Z production at 10 TeV # 2.2 GeV for LHC W/Z production at 14 TeV # # Set all parameters to 0 to disable ################################################################## newdef ShowerHandler:IntrinsicPtGaussian 1.3*GeV newdef ShowerHandler:IntrinsicPtBeta 0 newdef ShowerHandler:IntrinsicPtGamma 0*GeV newdef ShowerHandler:IntrinsicPtIptmax 0*GeV ############################################################# # Set up truncated shower handler. ############################################################# create Herwig::PowhegShowerHandler PowhegShowerHandler set PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler set PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF newdef PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler newdef PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF newdef PowhegShowerHandler:PartnerFinder PartnerFinder newdef PowhegShowerHandler:KinematicsReconstructor KinematicsReconstructor newdef PowhegShowerHandler:SplittingGenerator SplittingGenerator newdef PowhegShowerHandler:Interactions QCDandQED newdef PowhegShowerHandler:SpinCorrelations Yes newdef PowhegShowerHandler:SoftCorrelations Singular newdef PowhegShowerHandler:IntrinsicPtGaussian 1.3*GeV newdef PowhegShowerHandler:IntrinsicPtBeta 0 newdef PowhegShowerHandler:IntrinsicPtGamma 0*GeV newdef PowhegShowerHandler:IntrinsicPtIptmax 0*GeV newdef PowhegShowerHandler:ReconstructionOption OffShell5 ############################################################# # End of interesting user servicable section. # # Anything that follows below should only be touched if you # know what you're doing. # # Really. ############################################################# # # a few default values newdef ShowerHandler:MECorrMode 1 newdef ShowerHandler:ReconstructionOption OffShell5 newdef AlphaQCD:ScaleFactor 1.0 newdef AlphaQCD:NPAlphaS 2 newdef AlphaQCD:Qmin 0.935 newdef AlphaQCD:NumberOfLoops 2 newdef AlphaQCD:InputOption 1 newdef AlphaQCD:AlphaMZ 0.126234 # # # Lets set up all the splittings create Herwig::HalfHalfOneSplitFn QtoQGammaSplitFn set QtoQGammaSplitFn:InteractionType QED set QtoQGammaSplitFn:ColourStructure ChargedChargedNeutral set QtoQGammaSplitFn:AngularOrdered Yes create Herwig::HalfHalfOneSplitFn QtoQGSplitFn newdef QtoQGSplitFn:InteractionType QCD newdef QtoQGSplitFn:ColourStructure TripletTripletOctet set QtoQGSplitFn:AngularOrdered Yes create Herwig::OneOneOneSplitFn GtoGGSplitFn newdef GtoGGSplitFn:InteractionType QCD newdef GtoGGSplitFn:ColourStructure OctetOctetOctet set GtoGGSplitFn:AngularOrdered Yes create Herwig::OneOneOneMassiveSplitFn WtoWGammaSplitFn newdef WtoWGammaSplitFn:InteractionType QED newdef WtoWGammaSplitFn:ColourStructure ChargedChargedNeutral set WtoWGammaSplitFn:AngularOrdered Yes create Herwig::OneHalfHalfSplitFn GtoQQbarSplitFn newdef GtoQQbarSplitFn:InteractionType QCD newdef GtoQQbarSplitFn:ColourStructure OctetTripletTriplet set GtoQQbarSplitFn:AngularOrdered Yes create Herwig::OneHalfHalfSplitFn GammatoQQbarSplitFn newdef GammatoQQbarSplitFn:InteractionType QED newdef GammatoQQbarSplitFn:ColourStructure NeutralChargedCharged set GammatoQQbarSplitFn:AngularOrdered Yes create Herwig::HalfOneHalfSplitFn QtoGQSplitFn newdef QtoGQSplitFn:InteractionType QCD newdef QtoGQSplitFn:ColourStructure TripletOctetTriplet set QtoGQSplitFn:AngularOrdered Yes create Herwig::HalfOneHalfSplitFn QtoGammaQSplitFn newdef QtoGammaQSplitFn:InteractionType QED newdef QtoGammaQSplitFn:ColourStructure ChargedNeutralCharged set QtoGammaQSplitFn:AngularOrdered Yes # # Now the Sudakovs create Herwig::QTildeSudakov SudakovCommon newdef SudakovCommon:Alpha AlphaQCD newdef SudakovCommon:cutoffKinScale 0.0*GeV newdef SudakovCommon:PDFmax 1.0 newdef SudakovCommon:CutOffOption pT newdef SudakovCommon:pTmin 1.222798*GeV cp SudakovCommon QtoQGSudakov newdef QtoQGSudakov:SplittingFunction QtoQGSplitFn newdef QtoQGSudakov:PDFmax 1.9 cp SudakovCommon QtoQGammaSudakov set QtoQGammaSudakov:SplittingFunction QtoQGammaSplitFn set QtoQGammaSudakov:Alpha AlphaQED set QtoQGammaSudakov:PDFmax 1.9 cp QtoQGammaSudakov LtoLGammaSudakov # Technical parameter to stop evolution. set LtoLGammaSudakov:pTmin 0.000001 cp SudakovCommon GtoGGSudakov newdef GtoGGSudakov:SplittingFunction GtoGGSplitFn newdef GtoGGSudakov:PDFmax 2.0 cp SudakovCommon WtoWGammaSudakov newdef WtoWGammaSudakov:SplittingFunction WtoWGammaSplitFn set WtoWGammaSudakov:Alpha AlphaQED cp SudakovCommon GtoQQbarSudakov newdef GtoQQbarSudakov:SplittingFunction GtoQQbarSplitFn newdef GtoQQbarSudakov:PDFmax 120.0 cp SudakovCommon GammatoQQbarSudakov newdef GammatoQQbarSudakov:SplittingFunction GammatoQQbarSplitFn set GammatoQQbarSudakov:Alpha AlphaQED newdef GammatoQQbarSudakov:PDFmax 120.0 cp SudakovCommon GtobbbarSudakov newdef GtobbbarSudakov:SplittingFunction GtoQQbarSplitFn newdef GtobbbarSudakov:PDFmax 40000.0 cp SudakovCommon GtoccbarSudakov newdef GtoccbarSudakov:SplittingFunction GtoQQbarSplitFn newdef GtoccbarSudakov:PDFmax 2000.0 cp SudakovCommon QtoGQSudakov newdef QtoGQSudakov:SplittingFunction QtoGQSplitFn cp SudakovCommon QtoGammaQSudakov newdef QtoGammaQSudakov:SplittingFunction QtoGammaQSplitFn set QtoGammaQSudakov:Alpha AlphaQED cp SudakovCommon utoGuSudakov newdef utoGuSudakov:SplittingFunction QtoGQSplitFn newdef utoGuSudakov:PDFFactor OverOneMinusZ newdef utoGuSudakov:PDFmax 5.0 cp SudakovCommon dtoGdSudakov newdef dtoGdSudakov:SplittingFunction QtoGQSplitFn newdef dtoGdSudakov:PDFFactor OverOneMinusZ # # Now add the final splittings # do SplittingGenerator:AddFinalSplitting u->u,g; QtoQGSudakov do SplittingGenerator:AddFinalSplitting d->d,g; QtoQGSudakov do SplittingGenerator:AddFinalSplitting s->s,g; QtoQGSudakov do SplittingGenerator:AddFinalSplitting c->c,g; QtoQGSudakov do SplittingGenerator:AddFinalSplitting b->b,g; QtoQGSudakov do SplittingGenerator:AddFinalSplitting t->t,g; QtoQGSudakov # do SplittingGenerator:AddFinalSplitting g->g,g; GtoGGSudakov # do SplittingGenerator:AddFinalSplitting g->u,ubar; GtoQQbarSudakov do SplittingGenerator:AddFinalSplitting g->d,dbar; GtoQQbarSudakov do SplittingGenerator:AddFinalSplitting g->s,sbar; GtoQQbarSudakov do SplittingGenerator:AddFinalSplitting g->c,cbar; GtoccbarSudakov do SplittingGenerator:AddFinalSplitting g->b,bbar; GtobbbarSudakov do SplittingGenerator:AddFinalSplitting g->t,tbar; GtoQQbarSudakov # do SplittingGenerator:AddFinalSplitting gamma->u,ubar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->d,dbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->s,sbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->c,cbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->b,bbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->t,tbar; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->e-,e+; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->mu-,mu+; GammatoQQbarSudakov do SplittingGenerator:AddFinalSplitting gamma->tau-,tau+; GammatoQQbarSudakov # do SplittingGenerator:AddFinalSplitting u->u,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting d->d,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting s->s,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting c->c,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting b->b,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting t->t,gamma; QtoQGammaSudakov do SplittingGenerator:AddFinalSplitting e-->e-,gamma; LtoLGammaSudakov do SplittingGenerator:AddFinalSplitting mu-->mu-,gamma; LtoLGammaSudakov do SplittingGenerator:AddFinalSplitting tau-->tau-,gamma; LtoLGammaSudakov do SplittingGenerator:AddFinalSplitting W+->W+,gamma; WtoWGammaSudakov # # Now lets add the initial splittings. Remember the form a->b,c; means # that the current particle b is given and we backward branch to new # particle a which is initial state and new particle c which is final state # do SplittingGenerator:AddInitialSplitting u->u,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting d->d,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting s->s,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting c->c,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting b->b,g; QtoQGSudakov do SplittingGenerator:AddInitialSplitting u->u,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting d->d,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting s->s,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting c->c,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting b->b,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting t->t,gamma; QtoQGammaSudakov do SplittingGenerator:AddInitialSplitting g->g,g; GtoGGSudakov # do SplittingGenerator:AddInitialSplitting g->d,dbar; GtoQQbarSudakov do SplittingGenerator:AddInitialSplitting g->u,ubar; GtoQQbarSudakov do SplittingGenerator:AddInitialSplitting g->s,sbar; GtoQQbarSudakov do SplittingGenerator:AddInitialSplitting g->c,cbar; GtoccbarSudakov do SplittingGenerator:AddInitialSplitting g->b,bbar; GtobbbarSudakov # do SplittingGenerator:AddInitialSplitting gamma->d,dbar; GammatoQQbarSudakov do SplittingGenerator:AddInitialSplitting gamma->u,ubar; GammatoQQbarSudakov do SplittingGenerator:AddInitialSplitting gamma->s,sbar; GammatoQQbarSudakov do SplittingGenerator:AddInitialSplitting gamma->c,cbar; GammatoQQbarSudakov do SplittingGenerator:AddInitialSplitting gamma->b,bbar; GammatoQQbarSudakov # do SplittingGenerator:AddInitialSplitting d->g,d; dtoGdSudakov do SplittingGenerator:AddInitialSplitting u->g,u; utoGuSudakov do SplittingGenerator:AddInitialSplitting s->g,s; QtoGQSudakov do SplittingGenerator:AddInitialSplitting c->g,c; QtoGQSudakov do SplittingGenerator:AddInitialSplitting b->g,b; QtoGQSudakov do SplittingGenerator:AddInitialSplitting dbar->g,dbar; dtoGdSudakov do SplittingGenerator:AddInitialSplitting ubar->g,ubar; utoGuSudakov do SplittingGenerator:AddInitialSplitting sbar->g,sbar; QtoGQSudakov do SplittingGenerator:AddInitialSplitting cbar->g,cbar; QtoGQSudakov do SplittingGenerator:AddInitialSplitting bbar->g,bbar; QtoGQSudakov # do SplittingGenerator:AddInitialSplitting d->gamma,d; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting u->gamma,u; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting s->gamma,s; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting c->gamma,c; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting b->gamma,b; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting dbar->gamma,dbar; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting ubar->gamma,ubar; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting sbar->gamma,sbar; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting cbar->gamma,cbar; QtoGammaQSudakov do SplittingGenerator:AddInitialSplitting bbar->gamma,bbar; QtoGammaQSudakov