diff --git a/Decay/Partonic/WeakPartonicDecayer.cc b/Decay/Partonic/WeakPartonicDecayer.cc --- a/Decay/Partonic/WeakPartonicDecayer.cc +++ b/Decay/Partonic/WeakPartonicDecayer.cc @@ -1,637 +1,640 @@ // -*- C++ -*- // // WeakPartonicDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the WeakPartonicDecayer class. // #include "WeakPartonicDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Utilities/Kinematics.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/PDT/ConstituentParticleData.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h" using namespace Herwig; using namespace ThePEG::Helicity; -WeakPartonicDecayer::WeakPartonicDecayer() : MECode(0), _radprob(0.0), _maxtry(300), - _threemax(3.), _fourmax(3.) +WeakPartonicDecayer::WeakPartonicDecayer() : MECode(0), _radprob(0.0), _maxtry(300), + _threemax(3.), _fourmax(3.) {} IBPtr WeakPartonicDecayer::clone() const { return new_ptr(*this); } IBPtr WeakPartonicDecayer::fullclone() const { return new_ptr(*this); } bool WeakPartonicDecayer::accept(tcPDPtr parent, const tPDVector & prod) const { // check we can find the flavours of the quarks in the decaying meson long id = parent->id(); int flav1, flav2; if((id / 1000)%10) { flav1 = (id/1000)%10; flav2 = (id/10)%100; - } + } else { - flav1 = id/100; + flav1 = id/100; flav2 = (id/10)%10; } if(!flav1 || !flav2) return false; // if two decay products one must be in triplet and one antitriplet if(prod.size()==2) { if((prod[0]->iColour()==PDT::Colour3&&prod[1]->iColour()==PDT::Colour3bar)|| (prod[0]->iColour()==PDT::Colour3bar&&prod[1]->iColour()==PDT::Colour3)) return true; } else if(prod.size()==3) { if(((prod[0]->iColour()==PDT::Colour3 &&prod[2]->iColour()==PDT::Colour3bar)|| (prod[0]->iColour()==PDT::Colour3bar&&prod[2]->iColour()==PDT::Colour3)) &&prod[1]->iColour()==PDT::Colour8) return true; } else if(prod.size()==4) { // first two particles should be leptons or q qbar if((prod[0]->id()>=11&&prod[0]->id()<=16&&prod[1]->id()<=-11&&prod[1]->id()>=-16)|| (prod[1]->id()>=11&&prod[1]->id()<=16&&prod[0]->id()<=-11&&prod[0]->id()>=-16)|| (prod[0]->iColour()==PDT::Colour3 &&prod[1]->iColour()==PDT::Colour3bar )|| (prod[1]->iColour()==PDT::Colour3 &&prod[0]->iColour()==PDT::Colour3bar )) { // third particle quark and fourth colour anti-triplet or // thrid particle antiquark and fourth colour triplet if((prod[2]->iColour()==PDT::Colour3bar&&prod[3]->iColour()==PDT::Colour3 )|| (prod[2]->iColour()==PDT::Colour3 &&prod[3]->iColour()==PDT::Colour3bar)) return true; } else return false; } return false; } ParticleVector WeakPartonicDecayer::decay(const Particle & parent, const tPDVector & children) const { static tcPDPtr gluon=getParticleData(ParticleID::g); // make the particles ParticleVector partons; for(unsigned int ix=0;ixproduceParticle()); // these products have the mass but should have constituent mass partons[ix]->set5Momentum(Lorentz5Momentum(children[ix]->constituentMass())); } // 2-body decays if(partons.size()==2) { // no gluon if not select based on probability or if three body not allowed if(UseRandom::rnd()>_radprob|| parent.mass()constituentMass()+partons[0]->mass()+partons[1]->mass()) { double ctheta,phi; Lorentz5Momentum pout[2]; for(unsigned int ix=0;ix<2;++ix) pout[ix].setMass(partons[ix]->mass()); Kinematics::generateAngles(ctheta,phi); Kinematics::twoBodyDecay(parent.momentum(),pout[0].mass(),pout[1].mass(), ctheta,phi,pout[0],pout[1]); for(unsigned int ix=0; ix<2;++ix) partons[ix]->setMomentum(pout[ix]); if(partons[0]->dataPtr()->iColour()==PDT::Colour3) { partons[0]->antiColourNeighbour(partons[1]); } - else { + else { partons[0]-> colourNeighbour(partons[1]); } } else { Lorentz5Momentum pout[3]; for(unsigned int ix=0;ix<2;++ix) pout[ix].setMass(partons[ix]->mass()); // add gluon partons.push_back(gluon->produceParticle()); partons.back()->set5Momentum(gluon->constituentMass()); // momentum of gluon pout[2] = Lorentz5Momentum(gluon->constituentMass()); Kinematics::threeBodyDecay(parent.momentum(),pout[1],pout[0],pout[2]); for(unsigned int ix=0; ix<3;++ix) partons[ix]->setMomentum(pout[ix]); if(partons[0]->dataPtr()->iColour()==PDT::Colour3) { partons[0]->antiColourNeighbour(partons[2]); partons[1]-> colourNeighbour(partons[2]); } - else { + else { partons[0]-> colourNeighbour(partons[2]); partons[1]->antiColourNeighbour(partons[2]); } } } // 3-body decays else if(partons.size()==3) { // set masses of products Lorentz5Momentum pout[3],pin(parent.momentum()); for(unsigned int ix=0;ix<3;++ix) pout[ix].setMass(partons[ix]->mass()); double xs(partons[2]->mass()/pin.mass()),xb(1.-xs); pout[2]=xs*pin; // Get the particle quark that is decaying long idQ, idSpec; idSpec = partons[2]->id(); idQ = (parent.id()/1000)%10; if(!idQ) idQ = (parent.id()/100)%10; // Now the odd case of a B_c where the c decays, not the b if(idSpec == idQ) idQ = (parent.id()/10)%10; // momentum of the decaying quark PPtr inter = getParticleData(idQ)->produceParticle(parent.momentum()*xb); // two body decay of heavy quark double ctheta,phi; Kinematics::generateAngles(ctheta,phi); Kinematics::twoBodyDecay(inter->momentum(),pout[0].mass(),pout[1].mass(), ctheta,phi,pout[0],pout[1]); // set the momenta of the decay products for(unsigned int ix=0; ix<3;++ix) partons[ix]->setMomentum(pout[ix]); // make the colour connections // quark first if(partons[0]->data().iColour()==PDT::Colour3) { partons[0]->antiColourNeighbour(partons[1]); partons[1]->colourNeighbour(partons[0]); partons[1]->antiColourNeighbour(partons[2]); partons[2]->colourNeighbour(partons[1]); } // antiquark first else { partons[0]->colourNeighbour(partons[1]); partons[1]->antiColourNeighbour(partons[0]); partons[1]->colourNeighbour(partons[2]); partons[2]->antiColourNeighbour(partons[1]); } } // 4-body decays else if(partons.size()==4) { // swap 0 and 1 if needed if(partons[1]->dataPtr()->iColour()!=PDT::Colour0&& partons[1]->dataPtr()->iColour()!=partons[2]->dataPtr()->iColour()) swap(partons[0],partons[1]); // get the momenta of the decaying quark and the spectator Lorentz5Momentum pin(parent.momentum()); double xs(partons[3]->mass()/pin.mass()),xb(1.-xs); Lorentz5Momentum pspect(xs*pin),pdec(xb*pin); pspect.setMass(partons[3]->mass()); pdec.rescaleMass(); // Get the particle quark that is decaying long idQ, idSpec; idSpec = partons[3]->id(); idQ = (abs(parent.id())/1000)%10; if(!idQ) idQ = (abs(parent.id())/100)%10; // Now the odd case of a B_c where the c decays, not the b if(abs(idSpec) == idQ) idQ = (abs(parent.id())/10)%10; // change sign if spectator quark or antidiquark if((idSpec>0&&idSpec<6)||idSpec<-6) idQ = -idQ; // check if W products coloured bool Wcol = partons[0]->coloured(); // particle data object tcPDPtr dec = getParticleData(idQ); // spin density matrix for the decaying quark RhoDMatrix rhoin(PDT::Spin1Half,true); if(parent.dataPtr()->iSpin()!=PDT::Spin0 && parent.spinInfo()) { parent.spinInfo()->decay(); RhoDMatrix rhoHadron = parent.spinInfo()->rhoMatrix(); // particles with spin 0 diquark if(abs(parent.id())==5122 || abs(parent.id())==4122 || abs(parent.id())==5122 || abs(parent.id())==4122 || abs(parent.id())==5132 || abs(parent.id())==4132 || abs(parent.id())==5232 || abs(parent.id())==4232) { for(unsigned int ix=0;ix<2;++ix) rhoin(ix,ix) = rhoHadron(ix,ix); } // particles with spin 1 diquark else if(abs(parent.id())==5332 || abs(parent.id())==4332) { rhoin(0,0) = 2./3.*rhoHadron(1,1)+1./3.*rhoHadron(0,0); rhoin(1,1) = 2./3.*rhoHadron(0,0)+1./3.*rhoHadron(1,1); } } // momenta of the decay products vector pout(3,Lorentz5Momentum()); for(unsigned int ix=0;ix<3;++ix) pout[ix].setMass(partons[ix]->mass()); // charges of the exchanged boson and check if colour rearranged int c1 = dec ->iCharge()-partons[2]->dataPtr()->iCharge(); int c2 = partons[0]->dataPtr()->iCharge()+partons[1]->dataPtr()->iCharge(); bool rearranged = !(c1==c2&&abs(c1)==3); if(MECode==0) rearranged=false; if(rearranged) { int c3 = dec ->iCharge()-partons[1]->dataPtr()->iCharge(); int c4 = partons[0]->dataPtr()->iCharge()+partons[2]->dataPtr()->iCharge(); if(!(c3==c4&&abs(c3)==3)) { generator()->log() << "Unknown order for colour rearranged decay" << " in WeakPartonicDecayer::decay()\n"; generator()->log() << c1 << " " << c2 << " " << c3 << " " << c4 << "\n"; generator()->log() << parent << "\n" << dec->PDGName() << "\n"; for(unsigned int ix=0;ix<4;++ix) generator()->log() << *partons[ix] << "\n"; throw Exception() << "Unknown order for colour rearranged decay" << " in WeakPartonicDecayer::decay() " << Exception::runerror; } swap(pout[1] ,pout[2] ); swap(partons[1],partons[2]); } // decide if three or four body using prob bool threeBody = UseRandom::rnd() > _radprob; // if four body not kinematically possible must be three body if(pdec.mass()constituentMass()+pout[0].mass()+ pout[1].mass()+pout[2].mass()) threeBody=true; // if code ==0 always three body if(MECode==0) threeBody=true; // three body decay if( threeBody ) { if(MECode==0) { - Kinematics::threeBodyDecay(pdec,pout[1],pout[0],pout[2]); + Kinematics::threeBodyDecay(pdec,pout[1],pout[0],pout[2]); + // set momenta of particles + for(unsigned int ix=0;ixsetMomentum(pout[ix]); } else { // generate the kinematics double wgt(0.); Energy2 mb2max = sqr(pdec.mass() - pout[2].mass()); Energy2 mb2min = sqr(pout[0].mass() + pout[1].mass()); unsigned int ntry = 0; do { ++ntry; Energy2 mb2 = (mb2max-mb2min)*UseRandom::rnd()+mb2min; double CosAngle, AzmAngle; // perform first decay Lorentz5Momentum p01; p01.setMass(sqrt(mb2)); Kinematics::generateAngles(CosAngle,AzmAngle); Kinematics::twoBodyDecay(pdec,pout[2].mass(),p01.mass(), CosAngle,AzmAngle,pout[2],p01); // perform second decay Kinematics::generateAngles(CosAngle,AzmAngle); Kinematics::twoBodyDecay(p01,pout[0].mass(),pout[1].mass(), CosAngle,AzmAngle,pout[0],pout[1]); // kinematic piece of the weight - wgt = + wgt = Kinematics::pstarTwoBodyDecay(pdec.mass(),p01 .mass(),pout[2].mass())/pdec.mass()* Kinematics::pstarTwoBodyDecay(p01 .mass(),pout[0].mass(),pout[1].mass())/p01.mass(); // piece to improve weight variation (not kinematics dependent) // and integration over m23^2 (N.B. m23^2 fac on bottom for efficiency) wgt *= pdec.mass()/Kinematics::pstarTwoBodyDecay(pdec.mass(),sqrt(mb2min),pout[2].mass()) /(mb2max-mb2min)*sqr(pdec.mass()); // set momenta of particles for(unsigned int ix=0;ixsetMomentum(pout[ix]); // matrix element piece wgt *= threeBodyMatrixElement(dec,rhoin,pdec,partons); // check doesn't violate max if(wgt>_threemax) { ostringstream message; message << "Maximum weight for three-body decay " << "violated in WeakPartonicDecayer::decay()" << "Maximum = " << _threemax << " weight = " << wgt; generator()->logWarning( Exception(message.str(),Exception::warning) ); } } while( wgt < _threemax*UseRandom::rnd() && ntry < _maxtry ); - if(ntry==_maxtry) throw Exception() + if(ntry==_maxtry) throw Exception() << "Too many attempts to generate three body kinematics in " << "WeakPartonicDecayer::decay()" << Exception::eventerror; } partons[3]->setMomentum(pspect); // set up the colour connections if(rearranged) swap(partons[1],partons[2]); if(Wcol) { if(partons[0]->data().iColour()==PDT::Colour3) partons[0]->antiColourNeighbour(partons[1]); else partons[0]-> colourNeighbour(partons[1]); } if(partons[2]->data().iColour()==PDT::Colour3) { partons[2]->antiColourNeighbour(partons[3]); } else { partons[2]-> colourNeighbour(partons[3]); } - } + } // four body decay else { // generate the extra gluon partons.push_back(gluon->produceParticle()); partons.back()->set5Momentum(gluon->constituentMass()); // momentum of gluon pout.push_back(Lorentz5Momentum(gluon->constituentMass())); // generate the kinematics Energy2 ms2min(sqr(pout[0].mass()+pout[1].mass()+pout[2].mass())); Energy2 ms2max(sqr(pdec.mass()-pout[3].mass())); double wgt(0.); unsigned int ntry=0; bool initial = true; do { ++ntry; Energy2 ms2 = ms2min+UseRandom::rnd()*(ms2max-ms2min); Energy ms = sqrt(ms2); // and the W Energy2 mb2max = sqr(ms -pout[2].mass()); Energy2 mb2min = sqr(pout[0].mass()+pout[1].mass()); Energy2 mb2 = (mb2max-mb2min)*UseRandom::rnd()+mb2min; wgt = (mb2max-mb2min)/(ms2max-mb2min); // perform first decay Lorentz5Momentum ps; double CosAngle,AzmAngle; Kinematics::generateAngles(CosAngle,AzmAngle); Kinematics::twoBodyDecay(pdec,pout[3].mass(),ms,CosAngle,AzmAngle,pout[3],ps); // generate the kinematics // perform second decay Kinematics::generateAngles(CosAngle,AzmAngle); Lorentz5Momentum p01; p01.setMass(sqrt(mb2)); Kinematics::twoBodyDecay(ps,pout[2].mass(),p01.mass(), CosAngle,AzmAngle,pout[2],p01); // perform third decay Kinematics::generateAngles(CosAngle,AzmAngle); Kinematics::twoBodyDecay(p01,pout[0].mass(),pout[1].mass(), CosAngle,AzmAngle,pout[0],pout[1]); // kinematic piece of the weight wgt *= 16.* Kinematics::pstarTwoBodyDecay(pdec.mass(),pout[3].mass(),ms )/pdec.mass()* Kinematics::pstarTwoBodyDecay(ms ,p01 .mass(),pout[2].mass())/ms* Kinematics::pstarTwoBodyDecay(p01 .mass(),pout[0].mass(),pout[1].mass())/p01.mass(); wgt *= fourBodyMatrixElement(pdec,pout[2],pout[0],pout[1],pout[3],Wcol,initial); // check doesn't violate max if(wgt>_threemax) { ostringstream message; message << "Maximum weight for four-body decay " << "violated in WeakPartonicDecayer::decay()" << "Maximum = " << _fourmax << " weight = " << wgt; generator()->logWarning( Exception(message.str(),Exception::warning) ); } } while ( wgt < _fourmax*UseRandom::rnd() && ntry < _maxtry ); - if(ntry==_maxtry) throw Exception() + if(ntry==_maxtry) throw Exception() << "Too many attempts to generate four body kinematics in " << "WeakPartonicDecayer::decay()" << Exception::eventerror; // set momenta of particles for(unsigned int ix=0;ix<3;++ix) partons[ix]->setMomentum(pout[ix]); partons[3]->setMomentum(pspect); partons[4]->setMomentum(pout[3]); // special for tau leptons to get correlations threeBodyMatrixElement(dec,rhoin,pdec,partons); // set up the colour connections if(rearranged) swap(partons[1],partons[2]); // radiation from initial-state if(initial) { if(Wcol) { if(partons[0]->data().iColour()==PDT::Colour3) partons[0]->antiColourNeighbour(partons[1]); else partons[0]-> colourNeighbour(partons[1]); } if(partons[2]->data().iColour()==PDT::Colour3) { partons[2]->antiColourNeighbour(partons[4]); partons[3]-> colourNeighbour(partons[4]); } else { partons[2]-> colourNeighbour(partons[4]); partons[3]->antiColourNeighbour(partons[4]); } } // radiation from final-state else { if(partons[0]->data().iColour()==PDT::Colour3) { partons[0]->antiColourNeighbour(partons[4]); partons[1]-> colourNeighbour(partons[4]); } else { partons[0]-> colourNeighbour(partons[4]); partons[1]->antiColourNeighbour(partons[4]); } if(partons[2]->data().iColour()==PDT::Colour3) { partons[2]->antiColourNeighbour(partons[3]); } else { partons[2]-> colourNeighbour(partons[3]); } } } } return partons; } - + void WeakPartonicDecayer::persistentOutput(PersistentOStream & os) const { os << MECode << _radprob << _maxtry << _threemax << _fourmax; } void WeakPartonicDecayer::persistentInput(PersistentIStream & is, int) { is >> MECode >> _radprob >> _maxtry >> _threemax >> _fourmax; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigWeakPartonicDecayer("Herwig::WeakPartonicDecayer", "HwPartonicDecay.so"); void WeakPartonicDecayer::Init() { static ClassDocumentation documentation ("The WeakPartonicDecayer class performs partonic decays of hadrons containing a " "heavy quark."); static Switch interfaceMECode ("MECode", "The code for the type of matrix element to be used.", &WeakPartonicDecayer::MECode, 0, false, false); static SwitchOption interfaceMECodePhaseSpace (interfaceMECode, "PhaseSpace", "Phase space decays", 0); static SwitchOption interfaceMECodeWeak (interfaceMECode, "Weak", "Weak matrix element", 100); static Parameter interfaceRadiationProbability ("RadiationProbability", "The probability that QCD radiation produces an extra q qbar pair", &WeakPartonicDecayer::_radprob, 0., 0.0, 1., false, false, Interface::limited); static Parameter interfaceMaxTry ("MaxTry", "The maximum number of attempts to generate the kinematics", &WeakPartonicDecayer::_maxtry, 300, 10, 1000, false, false, Interface::limited); static Parameter interfaceThreeMax ("ThreeMax", "Maximum weight for sampling of three-body decays", &WeakPartonicDecayer::_threemax, 3.0, 1.0, 1000.0, false, false, Interface::limited); static Parameter interfaceFourMax ("FourMax", "Maximum weight for sampling of four-body decays", &WeakPartonicDecayer::_fourmax, 3.0, 1.0, 1000.0, false, false, Interface::limited); } double WeakPartonicDecayer::VAWt(Energy2 t0, Energy2 t1, Energy2 t2, InvEnergy4 t3) { - return (t1-t0)*(t0-t2)*t3; + return (t1-t0)*(t0-t2)*t3; } void WeakPartonicDecayer::dataBaseOutput(ofstream & output, bool header) const { if(header) output << "update decayers set parameters=\""; // parameters for the PartonicDecayerBase base class PartonicDecayerBase::dataBaseOutput(output,false); output << "newdef " << name() << ":MECode " << MECode << " \n"; - if(header) output << "\n\" where BINARY ThePEGName=\"" + if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } double WeakPartonicDecayer:: threeBodyMatrixElement(tcPDPtr dec, const RhoDMatrix & rhoin, Lorentz5Momentum & pdec, ParticleVector & partons) const { // spinors LorentzSpinor w0[2],w2[2]; LorentzSpinorBar w1[2],w3[2]; // spinors for the decaying particle and first product if(dec->id()>0) { SpinorWaveFunction win(pdec,dec,0,incoming); w0[0] = win.dimensionedWave(); win.reset(1); w0[1] = win.dimensionedWave(); SpinorBarWaveFunction wout(partons[2]->momentum(), partons[2]->dataPtr(),0,outgoing); w1[0] = wout.dimensionedWave(); wout.reset(1); w1[1] = wout.dimensionedWave(); } else { SpinorBarWaveFunction win(pdec,dec,0,incoming); w1[0] = win.dimensionedWave(); win.reset(1); w1[1] = win.dimensionedWave(); SpinorWaveFunction wout(partons[2]->momentum(), partons[2]->dataPtr(),0,outgoing); w0[0] = wout.dimensionedWave(); wout.reset(1); w0[1] = wout.dimensionedWave(); } // spinors for the W decay products bool lorder = true; if(partons[0]->id()<0) { SpinorWaveFunction wout2(partons[0]->momentum(), partons[0]->dataPtr(),0,outgoing); SpinorBarWaveFunction wout3(partons[1]->momentum(), partons[1]->dataPtr(),0,outgoing); lorder = partons[0]->dataPtr()->charged(); w2[0] = wout2.dimensionedWave(); w3[0] = wout3.dimensionedWave(); wout2.reset(1); wout3.reset(1); w2[1] = wout2.dimensionedWave(); w3[1] = wout3.dimensionedWave(); } else { SpinorWaveFunction wout2(partons[1]->momentum(), partons[1]->dataPtr(),0,outgoing); SpinorBarWaveFunction wout3(partons[0]->momentum(), partons[0]->dataPtr(),0,outgoing); lorder = partons[1]->dataPtr()->charged(); w2[0] = wout2.dimensionedWave(); w3[0] = wout3.dimensionedWave(); wout2.reset(1); wout3.reset(1); w2[1] = wout2.dimensionedWave(); w3[1] = wout3.dimensionedWave(); } bool tau = abs(partons[0]->id())==ParticleID::tauminus || abs(partons[1]->id())==ParticleID::tauminus; // calculate the currents LorentzPolarizationVectorE Jbc[2][2],Jdec[2][2]; for(unsigned int ix=0;ix<2;++ix) { for(unsigned int iy=0;iy<2;++iy) { if(dec->id()>0) Jbc [ix][iy] = w0[ix].leftCurrent(w1[iy]); else Jbc [ix][iy] = w0[iy].leftCurrent(w1[ix]); if(lorder) Jdec[ix][iy] = w2[ix].leftCurrent(w3[iy]); - else + else Jdec[ix][iy] = w2[iy].leftCurrent(w3[ix]); } } // compute the matrix element Complex me[2][2][2][2]; double total=0.; for(unsigned int i0=0;i0<2;++i0) { for(unsigned int i1=0;i1<2;++i1) { for(unsigned int i2=0;i2<2;++i2) { for(unsigned int i3=0;i3<2;++i3) { me[i0][i1][i2][i3] = Jbc[i0][i1].dot(Jdec[i2][i3])/sqr(pdec.mass()); total += rhoin(i0,i0).real()*norm(me[i0][i1][i2][i3]); } } } } total *=2.; if(tau) { RhoDMatrix rho(PDT::Spin1Half); for(unsigned int it1=0;it1<2;++it1) { for(unsigned int it2=0;it2<2;++it2) { for(unsigned int i0=0;i0<2;++i0) { for(unsigned int i1=0;i1<2;++i1) { for(unsigned int i2=0;i2<2;++i2) { rho(it1,it2) += me[i0][i1][it1][i2 ]*conj(me[i0][i1][it2][i2 ]); } } } } } // normalize matrix to unit trace rho.normalize(); for(unsigned int ix=0;ix<2;++ix) { if(abs(partons[ix]->id())!=ParticleID::tauminus) continue; bool loc = partons[ix]->id() < 0; // create the spin info object FermionSpinPtr spin = new_ptr(FermionSpinInfo(partons[ix]->momentum(),true)); // assign spinors for(unsigned int iy=0;iy<2;++iy) { spin->setBasisState(iy, loc ? w2[iy] : w3[iy].bar()); } // assign rho spin->rhoMatrix() = rho; // assign spin info partons[ix]->spinInfo(spin); } } return total; } double WeakPartonicDecayer:: fourBodyMatrixElement(Lorentz5Momentum & p0,Lorentz5Momentum & p1, Lorentz5Momentum & p2,Lorentz5Momentum & p3, Lorentz5Momentum & pg, bool Wcol, bool & initial) const { Energy2 d01(p0*p1),d02(p0*p2),d03(p0*p3),d0g(p0*pg); Energy2 d12(p1*p2),d13(p1*p3),d1g(p1*pg); Energy2 d23(p2*p3),d2g(p2*pg),d3g(p3*pg); Energy2 m02(sqr(p0.mass())),m12(sqr(p1.mass())),m22(sqr(p2.mass())), m32(sqr(p3.mass())); - Energy2 mei = + Energy2 mei = +1./d0g/d1g *( -d01*d12*d3g+d01*d03*d2g+2*d01*d03*d12 ) +1./d0g *( d12*d3g-d03*d12-d02*d03 ) +1./d1g *( d12*d13+d03*d2g+d03*d12 ) +m12/sqr(d1g)*( -d03*d2g-d03*d12 ) +m02/sqr(d0g)*( d12*d3g-d03*d12 ); - Energy2 mef = !Wcol ? ZERO : + Energy2 mef = !Wcol ? ZERO : +1./d2g/d3g *( d0g*d12*d23+d03*d1g*d23+2*d03*d12*d23 ) +1./d2g *( d03*d1g+d03*d12-d02*d12 ) +1./d3g *( d0g*d12-d03*d13+d03*d12 ) +m32/sqr(d3g)*( -d0g*d12-d03*d12 ) +m22/sqr(d2g)*( -d03*d1g-d03*d12 ); initial = mef/(mei+mef) #include #include #include #include #include #include #include "Herwig/Utilities/Kinematics.h" -#include "CheckId.h" #include "Cluster.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include using namespace Herwig; DescribeClass describeClusterFissioner("Herwig::ClusterFissioner","Herwig.so"); ClusterFissioner::ClusterFissioner() : _clMaxLight(3.35*GeV), _clMaxBottom(3.35*GeV), _clMaxCharm(3.35*GeV), _clMaxExotic(3.35*GeV), _clPowLight(2.0), _clPowBottom(2.0), _clPowCharm(2.0), _clPowExotic(2.0), _pSplitLight(1.0), _pSplitBottom(1.0), _pSplitCharm(1.0), _pSplitExotic(1.0), _fissionPwtUquark(1), _fissionPwtDquark(1), _fissionPwtSquark(0.5), _fissionCluster(0), _btClM(1.0*GeV), _iopRem(1), _kappa(1.0e15*GeV/meter), _enhanceSProb(0), _m0Fission(2.*GeV), _massMeasure(0) {} IBPtr ClusterFissioner::clone() const { return new_ptr(*this); } IBPtr ClusterFissioner::fullclone() const { return new_ptr(*this); } void ClusterFissioner::persistentOutput(PersistentOStream & os) const { os << _hadronsSelector << ounit(_clMaxLight,GeV) << ounit(_clMaxBottom,GeV) << ounit(_clMaxCharm,GeV) << ounit(_clMaxExotic,GeV) << _clPowLight << _clPowBottom << _clPowCharm << _clPowExotic << _pSplitLight << _pSplitBottom << _pSplitCharm << _pSplitExotic << _fissionCluster << _fissionPwtUquark << _fissionPwtDquark << _fissionPwtSquark << ounit(_btClM,GeV) << _iopRem << ounit(_kappa, GeV/meter) << _enhanceSProb << ounit(_m0Fission,GeV) << _massMeasure; } void ClusterFissioner::persistentInput(PersistentIStream & is, int) { is >> _hadronsSelector >> iunit(_clMaxLight,GeV) >> iunit(_clMaxBottom,GeV) >> iunit(_clMaxCharm,GeV) >> iunit(_clMaxExotic,GeV) >> _clPowLight >> _clPowBottom >> _clPowCharm >> _clPowExotic >> _pSplitLight >> _pSplitBottom >> _pSplitCharm >> _pSplitExotic >> _fissionCluster >> _fissionPwtUquark >> _fissionPwtDquark >> _fissionPwtSquark >> iunit(_btClM,GeV) >> _iopRem >> iunit(_kappa, GeV/meter) >> _enhanceSProb >> iunit(_m0Fission,GeV) >> _massMeasure; } void ClusterFissioner::Init() { static ClassDocumentation documentation ("Class responsibles for chopping up the clusters"); static Reference interfaceHadronSelector("HadronSelector", "A reference to the HadronSelector object", &Herwig::ClusterFissioner::_hadronsSelector, false, false, true, false); // ClMax for light, Bottom, Charm and exotic (e.g. Susy) quarks static Parameter interfaceClMaxLight ("ClMaxLight","cluster max mass for light quarks (unit [GeV])", &ClusterFissioner::_clMaxLight, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); static Parameter interfaceClMaxBottom ("ClMaxBottom","cluster max mass for b quarks (unit [GeV])", &ClusterFissioner::_clMaxBottom, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); static Parameter interfaceClMaxCharm ("ClMaxCharm","cluster max mass for c quarks (unit [GeV])", &ClusterFissioner::_clMaxCharm, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); static Parameter interfaceClMaxExotic ("ClMaxExotic","cluster max mass for exotic quarks (unit [GeV])", &ClusterFissioner::_clMaxExotic, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); // ClPow for light, Bottom, Charm and exotic (e.g. Susy) quarks static Parameter interfaceClPowLight ("ClPowLight","cluster mass exponent for light quarks", &ClusterFissioner::_clPowLight, 0, 2.0, 0.0, 10.0,false,false,false); static Parameter interfaceClPowBottom ("ClPowBottom","cluster mass exponent for b quarks", &ClusterFissioner::_clPowBottom, 0, 2.0, 0.0, 10.0,false,false,false); static Parameter interfaceClPowCharm ("ClPowCharm","cluster mass exponent for c quarks", &ClusterFissioner::_clPowCharm, 0, 2.0, 0.0, 10.0,false,false,false); static Parameter interfaceClPowExotic ("ClPowExotic","cluster mass exponent for exotic quarks", &ClusterFissioner::_clPowExotic, 0, 2.0, 0.0, 10.0,false,false,false); // PSplit for light, Bottom, Charm and exotic (e.g. Susy) quarks static Parameter interfacePSplitLight ("PSplitLight","cluster mass splitting param for light quarks", &ClusterFissioner::_pSplitLight, 0, 1.0, 0.0, 10.0,false,false,false); static Parameter interfacePSplitBottom ("PSplitBottom","cluster mass splitting param for b quarks", &ClusterFissioner::_pSplitBottom, 0, 1.0, 0.0, 10.0,false,false,false); static Parameter interfacePSplitCharm ("PSplitCharm","cluster mass splitting param for c quarks", &ClusterFissioner::_pSplitCharm, 0, 1.0, 0.0, 10.0,false,false,false); static Parameter interfacePSplitExotic ("PSplitExotic","cluster mass splitting param for exotic quarks", &ClusterFissioner::_pSplitExotic, 0, 1.0, 0.0, 10.0,false,false,false); static Switch interfaceFission ("Fission", "Option for different Fission options", &ClusterFissioner::_fissionCluster, 1, false, false); static SwitchOption interfaceFissionDefault (interfaceFission, "default", "Normal cluster fission which depends on the hadron selector class.", 0); static SwitchOption interfaceFissionNew (interfaceFission, "new", "Alternative cluster fission which does not depend on the hadron selector class", 1); static Parameter interfaceFissionPwtUquark ("FissionPwtUquark", "Weight for fission in U quarks", &ClusterFissioner::_fissionPwtUquark, 1, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceFissionPwtDquark ("FissionPwtDquark", "Weight for fission in D quarks", &ClusterFissioner::_fissionPwtDquark, 1, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceFissionPwtSquark ("FissionPwtSquark", "Weight for fission in S quarks", &ClusterFissioner::_fissionPwtSquark, 0.5, 0.0, 1.0, false, false, Interface::limited); static Switch interfaceRemnantOption ("RemnantOption", "Option for the treatment of remnant clusters", &ClusterFissioner::_iopRem, 1, false, false); static SwitchOption interfaceRemnantOptionSoft (interfaceRemnantOption, "Soft", "Both clusters produced in the fission of the beam cluster" " are treated as soft clusters.", 0); static SwitchOption interfaceRemnantOptionHard (interfaceRemnantOption, "Hard", "Only the cluster containing the remnant is treated as a soft cluster.", 1); static SwitchOption interfaceRemnantOptionVeryHard (interfaceRemnantOption, "VeryHard", "Even remnant clusters are treated as hard, i.e. all clusters the same", 2); static Parameter interfaceBTCLM ("SoftClusterFactor", "Parameter for the mass spectrum of remnant clusters", &ClusterFissioner::_btClM, GeV, 1.*GeV, 0.1*GeV, 10.0*GeV, false, false, Interface::limited); static Parameter interfaceStringTension ("StringTension", "String tension used in vertex displacement calculation", &ClusterFissioner::_kappa, GeV/meter, 1.0e15*GeV/meter, ZERO, ZERO, false, false, Interface::lowerlim); static Switch interfaceEnhanceSProb ("EnhanceSProb", "Option for enhancing strangeness", &ClusterFissioner::_enhanceSProb, 0, false, false); static SwitchOption interfaceEnhanceSProbNo (interfaceEnhanceSProb, "No", "No strangeness enhancement.", 0); static SwitchOption interfaceEnhanceSProbScaled (interfaceEnhanceSProb, "Scaled", "Scaled strangeness enhancement", 1); static SwitchOption interfaceEnhanceSProbExponential (interfaceEnhanceSProb, "Exponential", "Exponential strangeness enhancement", 2); static Switch interfaceMassMeasure ("MassMeasure", "Option to use different mass measures", &ClusterFissioner::_massMeasure,0,false,false); static SwitchOption interfaceMassMeasureMass (interfaceMassMeasure, "Mass", "Mass Measure", 0); static SwitchOption interfaceMassMeasureLambda (interfaceMassMeasure, "Lambda", "Lambda Measure", 1); static Parameter interfaceFissionMassScale ("FissionMassScale", "Cluster fission mass scale", &ClusterFissioner::_m0Fission, GeV, 2.0*GeV, 0.1*GeV, 50.*GeV, false, false, Interface::limited); } tPVector ClusterFissioner::fission(ClusterVector & clusters, bool softUEisOn) { // return if no clusters if (clusters.empty()) return tPVector(); /***************** * Loop over the (input) collection of cluster pointers, and store in * the vector splitClusters all the clusters that need to be split * (these are beam clusters, if soft underlying event is off, and * heavy non-beam clusters). ********************/ stack splitClusters; for(ClusterVector::iterator it = clusters.begin() ; it != clusters.end() ; ++it) { /************** * Skip 3-component clusters that have been redefined (as 2-component * clusters) or not available clusters. The latter check is indeed * redundant now, but it is used for possible future extensions in which, * for some reasons, some of the clusters found by ClusterFinder are tagged * straight away as not available. **************/ if((*it)->isRedefined() || !(*it)->isAvailable()) continue; // if the cluster is a beam cluster add it to the vector of clusters // to be split or if it is heavy if((*it)->isBeamCluster() || isHeavy(*it)) splitClusters.push(*it); } tPVector finalhadrons; cut(splitClusters, clusters, finalhadrons, softUEisOn); return finalhadrons; } void ClusterFissioner::cut(stack & clusterStack, ClusterVector &clusters, tPVector & finalhadrons, bool softUEisOn) { /************************************************** * This method does the splitting of the cluster pointed by cluPtr * and "recursively" by all of its cluster children, if heavy. All of these * new children clusters are added (indeed the pointers to them) to the * collection of cluster pointers collecCluPtr. The method works as follows. * Initially the vector vecCluPtr contains just the input pointer to the * cluster to be split. Then it will be filled "recursively" by all * of the cluster's children that are heavy enough to require, in their turn, * to be split. In each loop, the last element of the vector vecCluPtr is * considered (only once because it is then removed from the vector). * This approach is conceptually recursive, but avoid the overhead of * a concrete recursive function. Furthermore it requires minimal changes * in the case that the fission of an heavy cluster could produce more * than two cluster children as assumed now. * * Draw the masses: for normal, non-beam clusters a power-like mass dist * is used, whereas for beam clusters a fast-decreasing exponential mass * dist is used instead (to avoid many iterative splitting which could * produce an unphysical large transverse energy from a supposed soft beam * remnant process). ****************************************/ // Here we recursively loop over clusters in the stack and cut them while (!clusterStack.empty()) { // take the last element of the vector ClusterPtr iCluster = clusterStack.top(); clusterStack.pop(); // split it cutType ct = iCluster->numComponents() == 2 ? cutTwo(iCluster, finalhadrons, softUEisOn) : cutThree(iCluster, finalhadrons, softUEisOn); // There are cases when we don't want to split, even if it fails mass test if(!ct.first.first || !ct.second.first) { // if an unsplit beam cluster leave if for the underlying event if(iCluster->isBeamCluster() && softUEisOn) iCluster->isAvailable(false); continue; } // check if clusters ClusterPtr one = dynamic_ptr_cast(ct.first.first); ClusterPtr two = dynamic_ptr_cast(ct.second.first); // is a beam cluster must be split into two clusters if(iCluster->isBeamCluster() && (!one||!two) && softUEisOn) { iCluster->isAvailable(false); continue; } // There should always be a intermediate quark(s) from the splitting assert(ct.first.second && ct.second.second); /// \todo sort out motherless quark pairs here. Watch out for 'quark in final state' errors iCluster->addChild(ct.first.first); // iCluster->addChild(ct.first.second); // ct.first.second->addChild(ct.first.first); iCluster->addChild(ct.second.first); // iCluster->addChild(ct.second.second); // ct.second.second->addChild(ct.second.first); // Sometimes the clusters decay C -> H + C' rather then C -> C' + C'' if(one) { clusters.push_back(one); if(one->isBeamCluster() && softUEisOn) one->isAvailable(false); if(isHeavy(one) && one->isAvailable()) clusterStack.push(one); } if(two) { clusters.push_back(two); if(two->isBeamCluster() && softUEisOn) two->isAvailable(false); if(isHeavy(two) && two->isAvailable()) clusterStack.push(two); } } } namespace { /** * Check if can't make a hadron from the partons */ bool cantMakeHadron(tcPPtr p1, tcPPtr p2) { return ! CheckId::canBeHadron(p1->dataPtr(), p2->dataPtr()); } /** * Check if can't make a diquark from the partons */ bool cantMakeDiQuark(tcPPtr p1, tcPPtr p2) { long id1 = p1->id(), id2 = p2->id(); return ! (QuarkMatcher::Check(id1) && QuarkMatcher::Check(id2) && id1*id2>0); } } ClusterFissioner::cutType ClusterFissioner::cutTwo(ClusterPtr & cluster, tPVector & finalhadrons, bool softUEisOn) { // need to make sure only 2-cpt clusters get here assert(cluster->numComponents() == 2); tPPtr ptrQ1 = cluster->particle(0); tPPtr ptrQ2 = cluster->particle(1); Energy Mc = cluster->mass(); assert(ptrQ1); assert(ptrQ2); // And check if those particles are from a beam remnant bool rem1 = cluster->isBeamRemnant(0); bool rem2 = cluster->isBeamRemnant(1); // workout which distribution to use bool soft1(false),soft2(false); switch (_iopRem) { case 0: soft1 = rem1 || rem2; soft2 = rem2 || rem1; break; case 1: soft1 = rem1; soft2 = rem2; break; } // Initialization for the exponential ("soft") mass distribution. static const int max_loop = 1000; int counter = 0; Energy Mc1 = ZERO, Mc2 = ZERO,m1=ZERO,m2=ZERO,m=ZERO; tcPDPtr toHadron1, toHadron2; PPtr newPtr1 = PPtr (); PPtr newPtr2 = PPtr (); bool succeeded = false; + Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2; + // pClu1(Mc1), pClu2(Mc2), pQ1(m1), pQone(m), pQtwo(m), pQ2(m2); do { succeeded = false; ++counter; - if (_enhanceSProb == 0){ - drawNewFlavour(newPtr1,newPtr2); - } - else { - Energy2 mass2 = clustermass(cluster); - drawNewFlavourEnhanced(newPtr1,newPtr2,mass2); - } + + // get a flavour for the qqbar pair + drawNewFlavour(newPtr1,newPtr2,cluster); + // check for right ordering assert (ptrQ2); assert (newPtr2); assert (ptrQ2->dataPtr()); assert (newPtr2->dataPtr()); if(cantMakeHadron(ptrQ1, newPtr1) || cantMakeHadron(ptrQ2, newPtr2)) { swap(newPtr1, newPtr2); // check again if(cantMakeHadron(ptrQ1, newPtr1) || cantMakeHadron(ptrQ2, newPtr2)) { throw Exception() << "ClusterFissioner cannot split the cluster (" << ptrQ1->PDGName() << ' ' << ptrQ2->PDGName() << ") into hadrons.\n" << Exception::runerror; } } // Check that new clusters can produce particles and there is enough // phase space to choose the drawn flavour m1 = ptrQ1->data().constituentMass(); m2 = ptrQ2->data().constituentMass(); m = newPtr1->data().constituentMass(); // Do not split in the case there is no phase space available if(Mc < m1+m + m2+m) continue; - // power for splitting - double exp1=_pSplitLight; - double exp2=_pSplitLight; - if (CheckId::isExotic(ptrQ1->dataPtr())) exp1 = _pSplitExotic; - else if(CheckId::hasBottom(ptrQ1->dataPtr()))exp1 = _pSplitBottom; - else if(CheckId::hasCharm(ptrQ1->dataPtr())) exp1 = _pSplitCharm; + pQ1.setMass(m1); + pQone.setMass(m); + pQtwo.setMass(m); + pQ2.setMass(m2); - if (CheckId::isExotic(ptrQ2->dataPtr())) exp2 = _pSplitExotic; - else if(CheckId::hasBottom(ptrQ2->dataPtr())) exp2 = _pSplitBottom; - else if(CheckId::hasCharm(ptrQ2->dataPtr())) exp2 = _pSplitCharm; + pair res = drawNewMasses(Mc, soft1, soft2, pClu1, pClu2, + ptrQ1, pQ1, newPtr1, pQone, + newPtr2, pQtwo, ptrQ2, pQ2); - // If, during the drawing of candidate masses, too many attempts fail - // (because the phase space available is tiny) - /// \todo run separate loop here? - Mc1 = drawChildMass(Mc,m1,m2,m,exp1,soft1); - Mc2 = drawChildMass(Mc,m2,m1,m,exp2,soft2); + Mc1 = res.first; Mc2 = res.second; + if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > Mc) continue; /************************** * New (not present in Fortran Herwig): * check whether the fragment masses Mc1 and Mc2 are above the * threshold for the production of the lightest pair of hadrons with the * right flavours. If not, then set by hand the mass to the lightest * single hadron with the right flavours, in order to solve correctly * the kinematics, and (later in this method) create directly such hadron * and add it to the children hadrons of the cluster that undergoes the * fission (i.e. the one pointed by iCluPtr). Notice that in this special * case, the heavy cluster that undergoes the fission has one single * cluster child and one single hadron child. We prefer this approach, * rather than to create a light cluster, with the mass set equal to * the lightest hadron, and let then the class LightClusterDecayer to do * the job to decay it to that single hadron, for two reasons: * First, because the sum of the masses of the two constituents can be, * in this case, greater than the mass of that hadron, hence it would * be impossible to solve the kinematics for such two components, and * therefore we would have a cluster whose components are undefined. * Second, the algorithm is faster, because it avoids the reshuffling * procedure that would be necessary if we used LightClusterDecayer * to decay the light cluster to the lightest hadron. ****************************/ + + // override chosen masses if needed toHadron1 = _hadronsSelector->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc1); - if(toHadron1) Mc1 = toHadron1->mass(); + if(toHadron1) { Mc1 = toHadron1->mass(); pClu1.setMass(Mc1); } toHadron2 = _hadronsSelector->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc2); - if(toHadron2) Mc2 = toHadron2->mass(); + if(toHadron2) { Mc2 = toHadron2->mass(); pClu2.setMass(Mc2); } // if a beam cluster not allowed to decay to hadrons if(cluster->isBeamCluster() && (toHadron1||toHadron2) && softUEisOn) continue; // Check if the decay kinematics is still possible: if not then // force the one-hadron decay for the other cluster as well. if(Mc1 + Mc2 > Mc) { if(!toHadron1) { toHadron1 = _hadronsSelector->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc-Mc2); - if(toHadron1) Mc1 = toHadron1->mass(); + if(toHadron1) { Mc1 = toHadron1->mass(); pClu1.setMass(Mc1); } } else if(!toHadron2) { toHadron2 = _hadronsSelector->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc-Mc1); - if(toHadron2) Mc2 = toHadron2->mass(); + if(toHadron2) { Mc2 = toHadron2->mass(); pClu2.setMass(Mc2); } } } succeeded = (Mc >= Mc1+Mc2); } while (!succeeded && counter < max_loop); if(counter >= max_loop) { static const PPtr null = PPtr(); return cutType(PPair(null,null),PPair(null,null)); } // Determined the (5-components) momenta (all in the LAB frame) Lorentz5Momentum pClu = cluster->momentum(); // known Lorentz5Momentum p0Q1 = ptrQ1->momentum(); // known (mom Q1 before fission) - Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2; //unknown - pClu1.setMass(Mc1); - pClu2.setMass(Mc2); - pQ1.setMass(m1); - pQ2.setMass(m2); - pQone.setMass(m); - pQtwo.setMass(m); calculateKinematics(pClu,p0Q1,toHadron1,toHadron2, - pClu1,pClu2,pQ1,pQone,pQtwo,pQ2); // out + pClu1,pClu2,pQ1,pQone,pQtwo,pQ2); /****************** * The previous methods have determined the kinematics and positions * of C -> C1 + C2. * In the case that one of the two product is light, that means either * decayOneHadronClu1 or decayOneHadronClu2 is true, then the momenta * of the components of that light product have not been determined, * and a (light) cluster will not be created: the heavy father cluster * decays, in this case, into a single (not-light) cluster and a * single hadron. In the other, "normal", cases the father cluster * decays into two clusters, each of which has well defined components. * Notice that, in the case of components which point to particles, the * momenta of the components is properly set to the new values, whereas * we do not change the momenta of the pointed particles, because we * want to keep all of the information (that is the new momentum of a * component after the splitting, which is contained in the _momentum * member of the Component class, and the (old) momentum of that component * before the splitting, which is contained in the momentum of the * pointed particle). Please not make confusion of this only apparent * inconsistency! ********************/ LorentzPoint posC,pos1,pos2; posC = cluster->vertex(); calculatePositions(pClu, posC, pClu1, pClu2, pos1, pos2); cutType rval; if(toHadron1) { rval.first = produceHadron(toHadron1, newPtr1, pClu1, pos1); finalhadrons.push_back(rval.first.first); } else { rval.first = produceCluster(ptrQ1, newPtr1, pClu1, pos1, pQ1, pQone, rem1); } if(toHadron2) { rval.second = produceHadron(toHadron2, newPtr2, pClu2, pos2); finalhadrons.push_back(rval.second.first); } else { rval.second = produceCluster(ptrQ2, newPtr2, pClu2, pos2, pQ2, pQtwo, rem2); } return rval; } ClusterFissioner::cutType ClusterFissioner::cutThree(ClusterPtr & cluster, tPVector & finalhadrons, bool softUEisOn) { // need to make sure only 3-cpt clusters get here assert(cluster->numComponents() == 3); // extract quarks tPPtr ptrQ[3] = {cluster->particle(0),cluster->particle(1),cluster->particle(2)}; assert( ptrQ[0] && ptrQ[1] && ptrQ[2] ); // find maximum mass pair Energy mmax(ZERO); Lorentz5Momentum pDiQuark; int iq1(-1),iq2(-1); Lorentz5Momentum psum; for(int q1=0;q1<3;++q1) { psum+= ptrQ[q1]->momentum(); for(int q2=q1+1;q2<3;++q2) { Lorentz5Momentum ptest = ptrQ[q1]->momentum()+ptrQ[q2]->momentum(); ptest.rescaleMass(); Energy mass = ptest.m(); if(mass>mmax) { mmax = mass; pDiQuark = ptest; iq1 = q1; iq2 = q2; } } } // and the spectators int iother(-1); for(int ix=0;ix<3;++ix) if(ix!=iq1&&ix!=iq2) iother=ix; assert(iq1>=0&&iq2>=0&&iother>=0); // And check if those particles are from a beam remnant bool rem1 = cluster->isBeamRemnant(iq1); bool rem2 = cluster->isBeamRemnant(iq2); // workout which distribution to use bool soft1(false),soft2(false); switch (_iopRem) { case 0: soft1 = rem1 || rem2; soft2 = rem2 || rem1; break; case 1: soft1 = rem1; soft2 = rem2; break; } // Initialization for the exponential ("soft") mass distribution. static const int max_loop = 1000; int counter = 0; Energy Mc1 = ZERO, Mc2 = ZERO, m1=ZERO, m2=ZERO, m=ZERO; tcPDPtr toHadron; bool toDiQuark(false); PPtr newPtr1 = PPtr(),newPtr2 = PPtr(); PDPtr diquark; bool succeeded = false; + Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2; do { succeeded = false; ++counter; - if (_enhanceSProb == 0) { - drawNewFlavour(newPtr1,newPtr2); - } - else { - Energy2 mass2 = clustermass(cluster); - drawNewFlavourEnhanced(newPtr1,newPtr2, mass2); - } + + // get a flavour for the qqbar pair + drawNewFlavour(newPtr1,newPtr2,cluster); + // randomly pick which will be (anti)diquark and which a mesonic cluster if(UseRandom::rndbool()) { swap(iq1,iq2); swap(rem1,rem2); } // check first order if(cantMakeHadron(ptrQ[iq1], newPtr1) || cantMakeDiQuark(ptrQ[iq2], newPtr2)) { swap(newPtr1,newPtr2); } // check again if(cantMakeHadron(ptrQ[iq1], newPtr1) || cantMakeDiQuark(ptrQ[iq2], newPtr2)) { throw Exception() << "ClusterFissioner cannot split the cluster (" << ptrQ[iq1]->PDGName() << ' ' << ptrQ[iq2]->PDGName() << ") into a hadron and diquark.\n" << Exception::runerror; } // Check that new clusters can produce particles and there is enough // phase space to choose the drawn flavour m1 = ptrQ[iq1]->data().constituentMass(); m2 = ptrQ[iq2]->data().constituentMass(); m = newPtr1->data().constituentMass(); // Do not split in the case there is no phase space available if(mmax < m1+m + m2+m) continue; - // power for splitting - double exp1(_pSplitLight),exp2(_pSplitLight); - if (CheckId::isExotic (ptrQ[iq1]->dataPtr())) exp1 = _pSplitExotic; - else if(CheckId::hasBottom(ptrQ[iq1]->dataPtr())) exp1 = _pSplitBottom; - else if(CheckId::hasCharm (ptrQ[iq1]->dataPtr())) exp1 = _pSplitCharm; - if (CheckId::isExotic (ptrQ[iq2]->dataPtr())) exp2 = _pSplitExotic; - else if(CheckId::hasBottom(ptrQ[iq2]->dataPtr())) exp2 = _pSplitBottom; - else if(CheckId::hasCharm (ptrQ[iq2]->dataPtr())) exp2 = _pSplitCharm; + pQ1.setMass(m1); + pQone.setMass(m); + pQtwo.setMass(m); + pQ2.setMass(m2); - // If, during the drawing of candidate masses, too many attempts fail - // (because the phase space available is tiny) - /// \todo run separate loop here? - Mc1 = drawChildMass(mmax,m1,m2,m,exp1,soft1); - Mc2 = drawChildMass(mmax,m2,m1,m,exp2,soft2); + pair res = drawNewMasses(mmax, soft1, soft2, pClu1, pClu2, + ptrQ[iq1], pQ1, newPtr1, pQone, + newPtr2, pQtwo, ptrQ[iq1], pQ2); + + Mc1 = res.first; Mc2 = res.second; + if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > mmax) continue; + // check if need to force meson clster to hadron toHadron = _hadronsSelector->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),Mc1); - if(toHadron) Mc1 = toHadron->mass(); + if(toHadron) { Mc1 = toHadron->mass(); pClu1.setMass(Mc1); } // check if need to force diquark cluster to be on-shell toDiQuark = false; diquark = CheckId::makeDiquark(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr()); if(Mc2 < diquark->constituentMass()) { - Mc2 = diquark->constituentMass(); + Mc2 = diquark->constituentMass(); pClu2.setMass(Mc2); toDiQuark = true; } // if a beam cluster not allowed to decay to hadrons if(cluster->isBeamCluster() && toHadron && softUEisOn) continue; // Check if the decay kinematics is still possible: if not then // force the one-hadron decay for the other cluster as well. if(Mc1 + Mc2 > mmax) { if(!toHadron) { toHadron = _hadronsSelector->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),mmax-Mc2); - if(toHadron) Mc1 = toHadron->mass(); + if(toHadron) { Mc1 = toHadron->mass(); pClu1.setMass(Mc1); } } else if(!toDiQuark) { - Mc2 = _hadronsSelector->massLightestHadron(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr()); + Mc2 = _hadronsSelector->massLightestHadron(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr()); pClu2.setMass(Mc2); toDiQuark = true; } } succeeded = (mmax >= Mc1+Mc2); } while (!succeeded && counter < max_loop); // check no of tries if(counter >= max_loop) return cutType(); // Determine the (5-components) momenta (all in the LAB frame) Lorentz5Momentum p0Q1 = ptrQ[iq1]->momentum(); - // to be determined - Lorentz5Momentum pClu1(Mc1), pClu2(Mc2), pQ1(m1), pQone(m), pQtwo(m), pQ2(m2); calculateKinematics(pDiQuark,p0Q1,toHadron,toDiQuark, pClu1,pClu2,pQ1,pQone,pQtwo,pQ2); // positions of the new clusters LorentzPoint pos1,pos2; Lorentz5Momentum pBaryon = pClu2+ptrQ[iother]->momentum(); calculatePositions(cluster->momentum(), cluster->vertex(), pClu1, pBaryon, pos1, pos2); // first the mesonic cluster/meson cutType rval; if(toHadron) { rval.first = produceHadron(toHadron, newPtr1, pClu1, pos1); finalhadrons.push_back(rval.first.first); } else { rval.first = produceCluster(ptrQ[iq1], newPtr1, pClu1, pos1, pQ1, pQone, rem1); } if(toDiQuark) { rem2 |= cluster->isBeamRemnant(iother); PPtr newDiQuark = diquark->produceParticle(pClu2); rval.second = produceCluster(newDiQuark, ptrQ[iother], pBaryon, pos2, pClu2, ptrQ[iother]->momentum(), rem2); } else { rval.second = produceCluster(ptrQ[iq2], newPtr2, pBaryon, pos2, pQ2, pQtwo, rem2, ptrQ[iother],cluster->isBeamRemnant(iother)); } cluster->isAvailable(false); return rval; } ClusterFissioner::PPair ClusterFissioner::produceHadron(tcPDPtr hadron, tPPtr newPtr, const Lorentz5Momentum &a, const LorentzPoint &b) const { PPair rval; if(hadron->coloured()) { rval.first = (_hadronsSelector->lightestHadron(hadron,newPtr->dataPtr()))->produceParticle(); } else rval.first = hadron->produceParticle(); rval.second = newPtr; rval.first->set5Momentum(a); rval.first->setVertex(b); return rval; } ClusterFissioner::PPair ClusterFissioner::produceCluster(tPPtr ptrQ, tPPtr newPtr, const Lorentz5Momentum & a, const LorentzPoint & b, const Lorentz5Momentum & c, const Lorentz5Momentum & d, bool isRem, tPPtr spect, bool remSpect) const { PPair rval; rval.second = newPtr; ClusterPtr cluster = !spect ? new_ptr(Cluster(ptrQ,rval.second)) : new_ptr(Cluster(ptrQ,rval.second,spect)); rval.first = cluster; cluster->set5Momentum(a); cluster->setVertex(b); assert(cluster->particle(0)->id() == ptrQ->id()); cluster->particle(0)->set5Momentum(c); cluster->particle(1)->set5Momentum(d); cluster->setBeamRemnant(0,isRem); if(remSpect) cluster->setBeamRemnant(2,remSpect); return rval; } void ClusterFissioner::drawNewFlavour(PPtr& newPtrPos,PPtr& newPtrNeg) const { // Flavour is assumed to be only u, d, s, with weights // (which are not normalized probabilities) given // by the same weights as used in HadronsSelector for // the decay of clusters into two hadrons. double prob_d; double prob_u; double prob_s; switch(_fissionCluster){ case 0: prob_d = _hadronsSelector->pwtDquark(); prob_u = _hadronsSelector->pwtUquark(); prob_s = _hadronsSelector->pwtSquark(); break; case 1: prob_d = _fissionPwtDquark; prob_u = _fissionPwtUquark; prob_s = _fissionPwtSquark; break; default : assert(false); } int choice = UseRandom::rnd3(prob_u, prob_d, prob_s); long idNew = 0; switch (choice) { case 0: idNew = ThePEG::ParticleID::u; break; case 1: idNew = ThePEG::ParticleID::d; break; case 2: idNew = ThePEG::ParticleID::s; break; } newPtrPos = getParticle(idNew); newPtrNeg = getParticle(-idNew); assert (newPtrPos); assert(newPtrNeg); assert (newPtrPos->dataPtr()); assert(newPtrNeg->dataPtr()); } void ClusterFissioner::drawNewFlavourEnhanced(PPtr& newPtrPos,PPtr& newPtrNeg, Energy2 mass2) const { // Flavour is assumed to be only u, d, s, with weights // (which are not normalized probabilities) given // by the same weights as used in HadronsSelector for // the decay of clusters into two hadrons. double prob_d; double prob_u; double prob_s = 0.; double scale = abs(double(sqr(_m0Fission)/mass2)); // Choose which splitting weights you wish to use switch(_fissionCluster){ // 0: ClusterFissioner and ClusterDecayer use the same weights case 0: prob_d = _hadronsSelector->pwtDquark(); prob_u = _hadronsSelector->pwtUquark(); /* Strangeness enhancement: Case 1: probability scaling Case 2: Exponential scaling */ if (_enhanceSProb == 1) prob_s = (_maxScale < scale) ? 0. : pow(_hadronsSelector->pwtSquark(),scale); else if (_enhanceSProb == 2) prob_s = (_maxScale < scale) ? 0. : exp(-scale); break; /* 1: ClusterFissioner uses its own unique set of weights, i.e. decoupled from ClusterDecayer */ case 1: prob_d = _fissionPwtDquark; prob_u = _fissionPwtUquark; if (_enhanceSProb == 1) prob_s = (_maxScale < scale) ? 0. : pow(_fissionPwtSquark,scale); else if (_enhanceSProb == 2) prob_s = (_maxScale < scale) ? 0. : exp(-scale); break; } int choice = UseRandom::rnd3(prob_u, prob_d, prob_s); long idNew = 0; switch (choice) { case 0: idNew = ThePEG::ParticleID::u; break; case 1: idNew = ThePEG::ParticleID::d; break; case 2: idNew = ThePEG::ParticleID::s; break; } newPtrPos = getParticle(idNew); newPtrNeg = getParticle(-idNew); assert (newPtrPos); assert(newPtrNeg); assert (newPtrPos->dataPtr()); assert(newPtrNeg->dataPtr()); } -Energy2 ClusterFissioner::clustermass(const ClusterPtr & cluster){ +Energy2 ClusterFissioner::clustermass(const ClusterPtr & cluster) const{ Lorentz5Momentum pIn = cluster->momentum(); Energy2 endpointmass2 = sqr(cluster->particle(0)->mass() + cluster->particle(1)->mass()); Energy2 singletm2 = pIn.m2(); // Return either the cluster mass, or the lambda measure return (_massMeasure == 0) ? singletm2 : singletm2 - endpointmass2; } Energy ClusterFissioner::drawChildMass(const Energy M, const Energy m1, const Energy m2, const Energy m, const double expt, const bool soft) const { /*************************** * This method, given in input the cluster mass Mclu of an heavy cluster C, * made of consituents of masses m1 and m2, draws the masses Mclu1 and Mclu2 * of, respectively, the children cluster C1, made of constituent masses m1 * and m, and cluster C2, of mass Mclu2 and made of constituent masses m2 * and m. The mass is extracted from one of the two following mass * distributions: * --- power-like ("normal" distribution) * d(Prob) / d(M^exponent) = const * where the exponent can be different from the two children C1 (exp1) * and C2 (exponent2). * --- exponential ("soft" distribution) * d(Prob) / d(M^2) = exp(-b*M) * where b = 2.0 / average. * Such distributions are limited below by the masses of * the constituents quarks, and above from the mass of decaying cluster C. * The choice of which of the two mass distributions to use for each of the * two cluster children is dictated by iRemnant (see below). * If the number of attempts to extract a pair of mass values that are * kinematically acceptable is above some fixed number (max_loop, see below) * the method gives up and returns false; otherwise, when it succeeds, it * returns true. * * These distributions have been modified from HERWIG: * Before these were: * Mclu1 = m1 + (Mclu - m1 - m2)*pow( rnd(), 1.0/exponent1 ); * The new one coded here is a more efficient version, same density * but taking into account 'in phase space from' beforehand ***************************/ // hard cluster if(!soft) { return pow(UseRandom::rnd(pow((M-m1-m2-m)*UnitRemoval::InvE, expt), pow(m*UnitRemoval::InvE, expt)), 1./expt )*UnitRemoval::E + m1; } // Otherwise it uses a soft mass distribution else { static const InvEnergy b = 2.0 / _btClM; Energy max = M-m1-m2-2.0*m; double rmin = b*max; rmin = ( rmin < 50 ) ? exp(-rmin) : 0.; double r1; do { r1 = UseRandom::rnd(rmin, 1.0) * UseRandom::rnd(rmin, 1.0); } while (r1 < rmin); return m1 + m - log(r1)/b; } } void ClusterFissioner::calculateKinematics(const Lorentz5Momentum & pClu, const Lorentz5Momentum & p0Q1, const bool toHadron1, const bool toHadron2, Lorentz5Momentum & pClu1, Lorentz5Momentum & pClu2, Lorentz5Momentum & pQ1, Lorentz5Momentum & pQbar, Lorentz5Momentum & pQ, Lorentz5Momentum & pQ2bar) const { /****************** * This method solves the kinematics of the two body cluster decay: * C (Q1 Q2bar) ---> C1 (Q1 Qbar) + C2 (Q Q2bar) * In input we receive the momentum of C, pClu, and the momentum * of the quark Q1 (constituent of C), p0Q1, both in the LAB frame. * Furthermore, two boolean variables inform whether the two fission * products (C1, C2) decay immediately into a single hadron (in which * case the cluster itself is identify with that hadron) and we do * not have to solve the kinematics of the components (Q1,Qbar) for * C1 and (Q,Q2bar) for C2. * The output is given by the following momenta (all 5-components, * and all in the LAB frame): * pClu1 , pClu2 respectively of C1 , C2 * pQ1 , pQbar respectively of Q1 , Qbar in C1 * pQ , pQ2bar respectively of Q , Q2 in C2 * The assumption, suggested from the string model, is that, in C frame, * C1 and its constituents Q1 and Qbar are collinear, and collinear to * the direction of Q1 in C (that is before cluster decay); similarly, * (always in the C frame) C2 and its constituents Q and Q2bar are * collinear (and therefore anti-collinear with C1,Q1,Qbar). * The solution is then obtained by using Lorentz boosts, as follows. * The kinematics of C1 and C2 is solved in their parent C frame, * and then boosted back in the LAB. The kinematics of Q1 and Qbar * is solved in their parent C1 frame and then boosted back in the LAB; * similarly, the kinematics of Q and Q2bar is solved in their parent * C2 frame and then boosted back in the LAB. In each of the three * "two-body decay"-like cases, we use the fact that the direction * of the motion of the decay products is known in the rest frame of * their parent. This is obvious for the first case in which the * parent rest frame is C; but it is also true in the other two cases * where the rest frames are C1 and C2. This is because C1 and C2 * are boosted w.r.t. C in the same direction where their components, * respectively (Q1,Qbar) and (Q,Q2bar) move in C1 and C2 rest frame * respectively. * Of course, although the notation used assumed that C = (Q1 Q2bar) * where Q1 is a quark and Q2bar an antiquark, indeed everything remain * unchanged also in all following cases: * Q1 quark, Q2bar antiquark; --> Q quark; * Q1 antiquark , Q2bar quark; --> Q antiquark; * Q1 quark, Q2bar diquark; --> Q quark * Q1 antiquark, Q2bar anti-diquark; --> Q antiquark * Q1 diquark, Q2bar quark --> Q antiquark * Q1 anti-diquark, Q2bar antiquark; --> Q quark **************************/ // Calculate the unit three-vector, in the C frame, along which // all of the constituents and children clusters move. Lorentz5Momentum u(p0Q1); u.boost( -pClu.boostVector() ); // boost from LAB to C // the unit three-vector is then u.vect().unit() // Calculate the momenta of C1 and C2 in the (parent) C frame first, // where the direction of C1 is u.vect().unit(), and then boost back in the // LAB frame. if (pClu.m() < pClu1.mass() + pClu2.mass() ) { throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (A)" << Exception::eventerror; } Kinematics::twoBodyDecay(pClu, pClu1.mass(), pClu2.mass(), u.vect().unit(), pClu1, pClu2); // In the case that cluster1 does not decay immediately into a single hadron, // calculate the momenta of Q1 (as constituent of C1) and Qbar in the // (parent) C1 frame first, where the direction of Q1 is u.vect().unit(), // and then boost back in the LAB frame. if(!toHadron1) { if (pClu1.m() < pQ1.mass() + pQbar.mass() ) { throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (B)" << Exception::eventerror; } Kinematics::twoBodyDecay(pClu1, pQ1.mass(), pQbar.mass(), u.vect().unit(), pQ1, pQbar); } // In the case that cluster2 does not decay immediately into a single hadron, // Calculate the momenta of Q and Q2bar (as constituent of C2) in the // (parent) C2 frame first, where the direction of Q is u.vect().unit(), // and then boost back in the LAB frame. if(!toHadron2) { if (pClu2.m() < pQ.mass() + pQ2bar.mass() ) { throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (C)" << Exception::eventerror; } Kinematics::twoBodyDecay(pClu2, pQ.mass(), pQ2bar.mass(), u.vect().unit(), pQ, pQ2bar); } } void ClusterFissioner::calculatePositions(const Lorentz5Momentum & pClu, const LorentzPoint & positionClu, const Lorentz5Momentum & pClu1, const Lorentz5Momentum & pClu2, LorentzPoint & positionClu1, LorentzPoint & positionClu2) const { // Determine positions of cluster children. // See Marc Smith's thesis, page 127, formulas (4.122) and (4.123). Energy Mclu = pClu.m(); Energy Mclu1 = pClu1.m(); Energy Mclu2 = pClu2.m(); // Calculate the unit three-vector, in the C frame, along which // children clusters move. Lorentz5Momentum u(pClu1); u.boost( -pClu.boostVector() ); // boost from LAB to C frame // the unit three-vector is then u.vect().unit() Energy pstarChild = Kinematics::pstarTwoBodyDecay(Mclu,Mclu1,Mclu2); // First, determine the relative positions of the children clusters // in the parent cluster reference frame. Length x1 = ( 0.25*Mclu + 0.5*( pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa; Length t1 = Mclu/_kappa - x1; LorentzDistance distanceClu1( x1 * u.vect().unit(), t1 ); Length x2 = (-0.25*Mclu + 0.5*(-pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa; Length t2 = Mclu/_kappa + x2; LorentzDistance distanceClu2( x2 * u.vect().unit(), t2 ); // Then, transform such relative positions from the parent cluster // reference frame to the Lab frame. distanceClu1.boost( pClu.boostVector() ); distanceClu2.boost( pClu.boostVector() ); // Finally, determine the absolute positions in the Lab frame. positionClu1 = positionClu + distanceClu1; positionClu2 = positionClu + distanceClu2; } bool ClusterFissioner::isHeavy(tcClusterPtr clu) { // default double clpow = _clPowLight; Energy clmax = _clMaxLight; // particle data for constituents tcPDPtr cptr[3]={tcPDPtr(),tcPDPtr(),tcPDPtr()}; for(int ix=0;ixnumComponents(),3);++ix) { cptr[ix]=clu->particle(ix)->dataPtr(); } // different parameters for exotic, bottom and charm clusters if(CheckId::isExotic(cptr[0],cptr[1],cptr[1])) { clpow = _clPowExotic; clmax = _clMaxExotic; } else if(CheckId::hasBottom(cptr[0],cptr[1],cptr[1])) { clpow = _clPowBottom; clmax = _clMaxBottom; } else if(CheckId::hasCharm(cptr[0],cptr[1],cptr[1])) { clpow = _clPowCharm; clmax = _clMaxCharm; } bool aboveCutoff = ( pow(clu->mass()*UnitRemoval::InvE , clpow) > pow(clmax*UnitRemoval::InvE, clpow) + pow(clu->sumConstituentMasses()*UnitRemoval::InvE, clpow) ); // required test for SUSY clusters, since aboveCutoff alone // cannot guarantee (Mc > m1 + m2 + 2*m) in cut() static const Energy minmass = getParticleData(ParticleID::d)->constituentMass(); bool canSplitMinimally = clu->mass() > clu->sumConstituentMasses() + 2.0 * minmass; return aboveCutoff && canSplitMinimally; } diff --git a/Hadronization/ClusterFissioner.h b/Hadronization/ClusterFissioner.h --- a/Hadronization/ClusterFissioner.h +++ b/Hadronization/ClusterFissioner.h @@ -1,439 +1,490 @@ // -*- C++ -*- // // ClusterFissioner.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_ClusterFissioner_H #define HERWIG_ClusterFissioner_H #include #include "CluHadConfig.h" #include "HadronSelector.h" #include "ClusterFissioner.fh" +#include "CheckId.h" namespace Herwig { using namespace ThePEG; //class Cluster; // forward declaration /** \ingroup Hadronization * \class ClusterFissioner * \brief This class handles clusters which are too heavy. * \author Philip Stephens * \author Alberto Ribon * \author Stefan Gieseke * * This class does the job of chopping up either heavy clusters or beam * clusters in two lighter ones. The procedure is repeated recursively until * all of the cluster children have masses below some threshold values. * * For the beam remnant clusters, at the moment what is done is the following. * In the case that the soft underlying event is switched on, the * beam remnant clusters are tagged as not available, * therefore they will not be treated at all during the hadronization. * In the case instead that the soft underlying event is switched off, * then the beam remnant clusters are treated exactly as "normal" clusters, * with the only exception of the mass spectrum used to generate the * cluster children masses. For non-beam clusters, the masses of the cluster * children are draw from a power-like mass distribution; for beam clusters, * according to the value of the flag _IOpRem, either both * children masses are draw from a fast-decreasing exponential mass * distribution (case _IOpRem == 0, or, indendently by * _IOpRem, in the special case that the beam cluster contains two * beam remnants), or one mass from the exponential distribution (corresponding * of the cluster child with the beam remnant) and the other with the usual * power-like distribution (case _IOpRem == 1, which is the * default one, as in Herwig 6.3). * * The reason behind the use of a fast-decreasing exponential distribution * is that to avoid a large transverse energy from the many sequential * fissions that would otherwise occur due to the typical large cluster * mass of beam clusters. Using instead an exponential distribution * the masses of the two cluster children will be very small (order of * GeV). * * The rationale behind the implementation of the splitting of clusters * has been to preserve *all* of the information about such splitting * process. More explicitly a ThePEG::Step class is passed in and the * new clusters are added to the step as the decay products of the * heavy cluster. This approach has the twofold * advantage to provide all of the information that could be needed * (expecially in future developments), without any information loss, * and furthermore it allows a better debugging. * * @see HadronSelector * @see \ref ClusterFissionerInterfaces "The interfaces" * defined for ClusterFissioner. */ class ClusterFissioner: public Interfaced { public: /** @name Standard constructors and destructors. */ //@{ /** * Default constructor. */ ClusterFissioner(); //@} /** Splits the clusters which are too heavy. * * Split either heavy clusters or beam clusters recursively until all * children have mass below some threshold. Heavy clusters are those that * satisfy the condition * \f[ M^P > C^P + S^P \f] * where \f$ M \f$ is the clusters mass, \f$ P \f$ is the parameter * ClPow, \f$ C \f$ is the parameter ClMax and \f$ S \f$ is the * sum of the clusters constituent partons. * For beam clusters, they are split only if the soft underlying event * is switched off, otherwise these clusters will be tagged as unavailable * and they will not be treated by the hadronization altogether. * In the case beam clusters will be split, the procedure is exactly * the same as for normal non-beam clusters, with the only exception * of the mass spectrum from which to draw the masses of the two * cluster children (see method drawChildrenMasses for details). */ tPVector fission(ClusterVector & clusters, bool softUEisOn); 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); //@} /** * Standard Init function used to initialize the interfaces. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} private: /** * Private and non-existent assignment operator. */ ClusterFissioner & operator=(const ClusterFissioner &) = delete; /** * This method directs the splitting of the heavy clusters * * This method does the splitting of the clusters and all of its cluster * children, if heavy. All of these new children clusters are added to the * collection of clusters. The method works as follows. * Initially the vector contains just the stack of input pointers to the * clusters to be split. Then it will be filled recursively by all * of the cluster's children that are heavy enough to require * to be split. In each loop, the last element of the vector is * considered (only once because it is then removed from the vector). * * \todo is the following still true? * For normal, non-beam clusters, a power-like mass distribution * is used, whereas for beam clusters a fast-decreasing exponential mass * distribution is used instead. This avoids many iterative splitting which * could produce an unphysical large transverse energy from a supposed * soft beam remnant process. */ void cut(stack &, ClusterVector&, tPVector & finalhadrons, bool softUEisOn); public: /** * Definition for easy passing of two particles. */ typedef pair PPair; /** * Definition for use in the cut function. */ typedef pair cutType; /** * Splits the input cluster. * * Split the input cluster (which can be either an heavy non-beam * cluster or a beam cluster). The result is two pairs of particles. The * first element of each pair is new cluster/hadron, while the second * element of each pair is the particle drawn from the vacuum to create * the new cluster/hadron. * Notice that this method treats also beam clusters by using a different * mass spectrum used to generate the cluster child masses (see method * drawChildMass). */ //@{ /** * Split two-component cluster */ virtual cutType cutTwo(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn); /** * Split three-component cluster */ virtual cutType cutThree(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn); //@} public: /** * Produces a hadron and returns the flavour drawn from the vacuum. * * This routine produces a new hadron. It * also sets the momentum and vertex to the values given. */ PPair produceHadron(tcPDPtr hadron, tPPtr newPtr, const Lorentz5Momentum &a, const LorentzPoint &b) const; + +protected: + + /** + * Function that returns either the cluster mass or the lambda measure + */ + Energy2 clustermass(const ClusterPtr & cluster) const; + + /** + * Draw a new flavour for the given cluster; currently defaults to + * the default model + */ + virtual void drawNewFlavour(PPtr& newPtr1, PPtr& newPtr2, const ClusterPtr & cluster) const { + if (_enhanceSProb == 0){ + drawNewFlavour(newPtr1,newPtr2); + } + else { + drawNewFlavourEnhanced(newPtr1,newPtr2,clustermass(cluster)); + } + } + + /** + * Calculate the masses and possibly kinematics of the cluster + * fission at hand; if claculateKineamtics is perfomring non-trivial + * steps kinematics claulcated here will be overriden. Currentl;y resorts to the default + */ + virtual pair drawNewMasses(Energy Mc, bool soft1, bool soft2, + Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2, + tPPtr ptrQ1, Lorentz5Momentum& pQ1, + tPPtr, Lorentz5Momentum& pQone, + tPPtr, Lorentz5Momentum& pQtwo, + tPPtr ptrQ2, Lorentz5Momentum& pQ2) const { + + pair result; + + double exp1=_pSplitLight; + double exp2=_pSplitLight; + + if (CheckId::isExotic(ptrQ1->dataPtr())) exp1 = _pSplitExotic; + else if(CheckId::hasBottom(ptrQ1->dataPtr()))exp1 = _pSplitBottom; + else if(CheckId::hasCharm(ptrQ1->dataPtr())) exp1 = _pSplitCharm; + + if (CheckId::isExotic(ptrQ2->dataPtr())) exp2 = _pSplitExotic; + else if(CheckId::hasBottom(ptrQ2->dataPtr())) exp2 = _pSplitBottom; + else if(CheckId::hasCharm(ptrQ2->dataPtr())) exp2 = _pSplitCharm; + + result.first = drawChildMass(Mc,pQ1.mass(),pQ2.mass(),pQone.mass(),exp1,soft1); + result.second = drawChildMass(Mc,pQ2.mass(),pQ1.mass(),pQtwo.mass(),exp2,soft2); + + pClu1.setMass(result.first); + pClu2.setMass(result.second); + + return result; + + } + + /** + * Calculate the final kinematics of a heavy cluster decay C->C1 + + * C2, if not already performed by drawNewMasses + */ + virtual void calculateKinematics(const Lorentz5Momentum &pClu, + const Lorentz5Momentum &p0Q1, + const bool toHadron1, const bool toHadron2, + Lorentz5Momentum &pClu1, Lorentz5Momentum &pClu2, + Lorentz5Momentum &pQ1, Lorentz5Momentum &pQb, + Lorentz5Momentum &pQ2, Lorentz5Momentum &pQ2b) const; + protected: /** * Produces a cluster from the flavours passed in. * * This routine produces a new cluster with the flavours given by ptrQ and newPtr. * The new 5 momentum is a and the parent momentum are c and d. C is for the * ptrQ and d is for the new particle newPtr. rem specifies whether the existing * particle is a beam remnant or not. */ PPair produceCluster(tPPtr ptrQ, tPPtr newPtr, const Lorentz5Momentum &a, const LorentzPoint &b, const Lorentz5Momentum &c, const Lorentz5Momentum &d, const bool rem, tPPtr spect=tPPtr(), bool remSpect=false) const; /** * Returns the new quark-antiquark pair * needed for fission of a heavy cluster. Equal probabilities * are assumed for producing u, d, or s pairs. */ - void drawNewFlavour(PPtr& newPtrPos,PPtr& newPtrNeg) const; + void drawNewFlavour(PPtr& newPtrPos, PPtr& newPtrNeg) const; /** * Returns the new quark-antiquark pair * needed for fission of a heavy cluster. Equal probabilities * are assumed for producing u, d, or s pairs. * Extra argument is used when performing strangeness enhancement */ void drawNewFlavourEnhanced(PPtr& newPtrPos,PPtr& newPtrNeg, Energy2 mass2) const; /** * Produces the mass of a child cluster. * * Draw the masses \f$M'\f$ of the the cluster child produced * by the fission of an heavy cluster (of mass M). m1, m2 are the masses * of the constituents of the cluster; m is the mass of the quark extract * from the vacuum (together with its antiparticle). The algorithm produces * the mass of the cluster formed with consituent m1. * Two mass distributions can be used for the child cluster mass: * -# power-like mass distribution ("normal" mass) with power exp * \f[ M' = {\rm rnd}((M-m_1-m_2-m)^P, m^p)^{1/P} + m_1 \f] * where \f$ P \f$ is a parameter of the model and \f$ \rm{rnd} \f$ is * the function: * \f[ \rm{rnd}(a,b) = (1-r)a + r b \f] * and here \f$ r \f$ is a random number [0,1]. * -# fast-decreasing exponential mass distribution ("soft" mass) with * rmin. rmin is given by * \f[ r_{\rm min} = \exp(-b (M - m_1 - m_2 - 2 m)) \f] * where \f$ b \f$ is a parameter of the model. The generated mass is * given by * \f[ M' = m_1 + m - \frac{\log\left( * {\rm rnd}(r_{\rm min}, 1-r_{\rm min})\right)}{b} \f]. * * The choice of which mass distribution should be used for each of the two * cluster children is dictated by the parameter soft. */ Energy drawChildMass(const Energy M, const Energy m1, const Energy m2, const Energy m, const double exp, const bool soft) const; /** - * Determines the kinematics of a heavy cluster decay C->C1 + C2 - */ - void calculateKinematics(const Lorentz5Momentum &pClu, - const Lorentz5Momentum &p0Q1, - const bool toHadron1, const bool toHadron2, - Lorentz5Momentum &pClu1, Lorentz5Momentum &pClu2, - Lorentz5Momentum &pQ1, Lorentz5Momentum &pQb, - Lorentz5Momentum &pQ2, Lorentz5Momentum &pQ2b) const; - - /** * Determine the positions of the two children clusters. * * This routine generates the momentum of the decay products. It also * generates the momentum in the lab frame of the partons drawn out of * the vacuum. */ void calculatePositions(const Lorentz5Momentum &pClu, const LorentzPoint & positionClu, const Lorentz5Momentum & pClu1, const Lorentz5Momentum & pClu2, LorentzPoint & positionClu1, LorentzPoint & positionClu2 ) const; - - - /** - * Function that returns either the cluster mass or the lambda measure - */ - Energy2 clustermass(const ClusterPtr & cluster); - protected: /** @name Access members for child classes. */ //@{ /** * Access to the hadron selector */ HadronSelectorPtr hadronsSelector() const {return _hadronsSelector;} /** * Access to soft-cluster parameter */ Energy btClM() const {return _btClM;} /** * Cluster splitting paramater for light quarks */ double pSplitLight() const {return _pSplitLight;} /** * Cluster splitting paramater for bottom quarks */ double pSplitBottom() const {return _pSplitBottom;} /** * Cluster splitting paramater for charm quarks */ double pSplitCharm() const {return _pSplitCharm;} /** * Cluster splitting paramater for exotic particles */ double pSplitExotic() const {return _pSplitExotic;} //@} private: /** * Check if a cluster is heavy enough to split again */ bool isHeavy(tcClusterPtr ); /** * A pointer to a Herwig::HadronSelector object for generating hadrons. */ HadronSelectorPtr _hadronsSelector; /** * @name The Cluster max mass,dependant on which quarks are involved, used to determine when * fission will occur. */ //@{ Energy _clMaxLight; Energy _clMaxBottom; Energy _clMaxCharm; Energy _clMaxExotic; //@} /** * @name The power used to determine when cluster fission will occur. */ //@{ double _clPowLight; double _clPowBottom; double _clPowCharm; double _clPowExotic; //@} /** * @name The power, dependant on whic quarks are involved, used in the cluster mass generation. */ //@{ double _pSplitLight; double _pSplitBottom; double _pSplitCharm; double _pSplitExotic; // weights for alternaive cluster fission double _fissionPwtUquark; double _fissionPwtDquark; double _fissionPwtSquark; /** * Flag used to determine between normal cluster fission and alternative cluster fission */ int _fissionCluster; //@} /** * Parameter used (2/b) for the beam cluster mass generation. * Currently hard coded value. */ Energy _btClM; /** * Flag used to determine what distributions to use for the cluster masses. */ int _iopRem; /** * The string constant */ Tension _kappa; /** * Flag that switches between no strangeness enhancement, scaling enhancement, * and exponential enhancement (in numerical order) */ int _enhanceSProb; /** * Parameter that governs the strangeness enhancement scaling */ Energy _m0Fission; /** * Flag that switches between mass measures used in strangeness enhancement: * cluster mass, or the lambda measure - ( m_{clu}^2 - (m_q + m_{qbar})^2 ) */ int _massMeasure; /** * Constant variable which stops the scale from being to large, and not worth * calculating */ const double _maxScale = 20.; }; } #endif /* HERWIG_ClusterFissioner_H */ diff --git a/Hadronization/ClusterHadronizationHandler.cc b/Hadronization/ClusterHadronizationHandler.cc --- a/Hadronization/ClusterHadronizationHandler.cc +++ b/Hadronization/ClusterHadronizationHandler.cc @@ -1,320 +1,463 @@ // -*- C++ -*- // // ClusterHadronizationHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ClusterHadronizationHandler class. // #include "ClusterHadronizationHandler.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "Herwig/Utilities/EnumParticles.h" #include "CluHadConfig.h" #include "Cluster.h" #include using namespace Herwig; ClusterHadronizationHandler * ClusterHadronizationHandler::currentHandler_ = 0; DescribeClass describeClusterHadronizationHandler("Herwig::ClusterHadronizationHandler","Herwig.so"); IBPtr ClusterHadronizationHandler::clone() const { return new_ptr(*this); } IBPtr ClusterHadronizationHandler::fullclone() const { return new_ptr(*this); } void ClusterHadronizationHandler::persistentOutput(PersistentOStream & os) const { os << _partonSplitter << _clusterFinder << _colourReconnector << _clusterFissioner << _lightClusterDecayer << _clusterDecayer + << reshuffle_ << reshuffleMode_ << gluonMassGenerator_ << ounit(_minVirtuality2,GeV2) << ounit(_maxDisplacement,mm) << _underlyingEventHandler << _reduceToTwoComponents; } void ClusterHadronizationHandler::persistentInput(PersistentIStream & is, int) { is >> _partonSplitter >> _clusterFinder >> _colourReconnector >> _clusterFissioner >> _lightClusterDecayer >> _clusterDecayer + >> reshuffle_ >> reshuffleMode_ >> gluonMassGenerator_ >> iunit(_minVirtuality2,GeV2) >> iunit(_maxDisplacement,mm) >> _underlyingEventHandler >> _reduceToTwoComponents; } void ClusterHadronizationHandler::Init() { static ClassDocumentation documentation ("This is the main handler class for the Cluster Hadronization", "The hadronization was performed using the cluster model of \\cite{Webber:1983if}.", "%\\cite{Webber:1983if}\n" "\\bibitem{Webber:1983if}\n" " B.~R.~Webber,\n" " ``A QCD Model For Jet Fragmentation Including Soft Gluon Interference,''\n" " Nucl.\\ Phys.\\ B {\\bf 238}, 492 (1984).\n" " %%CITATION = NUPHA,B238,492;%%\n" // main manual ); static Reference interfacePartonSplitter("PartonSplitter", "A reference to the PartonSplitter object", &Herwig::ClusterHadronizationHandler::_partonSplitter, false, false, true, false); static Reference interfaceClusterFinder("ClusterFinder", "A reference to the ClusterFinder object", &Herwig::ClusterHadronizationHandler::_clusterFinder, false, false, true, false); static Reference interfaceColourReconnector("ColourReconnector", "A reference to the ColourReconnector object", &Herwig::ClusterHadronizationHandler::_colourReconnector, false, false, true, false); static Reference interfaceClusterFissioner("ClusterFissioner", "A reference to the ClusterFissioner object", &Herwig::ClusterHadronizationHandler::_clusterFissioner, false, false, true, false); static Reference interfaceLightClusterDecayer("LightClusterDecayer", "A reference to the LightClusterDecayer object", &Herwig::ClusterHadronizationHandler::_lightClusterDecayer, false, false, true, false); static Reference interfaceClusterDecayer("ClusterDecayer", "A reference to the ClusterDecayer object", &Herwig::ClusterHadronizationHandler::_clusterDecayer, false, false, true, false); + static Reference interfaceGluonMassGenerator + ("GluonMassGenerator", + "Set a reference to a gluon mass generator.", + &ClusterHadronizationHandler::gluonMassGenerator_, false, false, true, true, false); + + static Switch interfaceReshuffle + ("Reshuffle", + "Perform reshuffling if constituent masses have not yet been included by the shower", + &ClusterHadronizationHandler::reshuffle_, false, false, false); + static SwitchOption interfaceReshuffleYes + (interfaceReshuffle, + "Global", + "Do reshuffle.", + true); + static SwitchOption interfaceReshuffleNo + (interfaceReshuffle, + "No", + "Do not reshuffle.", + false); + + static Switch interfaceReshuffleMode + ("ReshuffleMode", + "Which mode is used for the reshuffling to constituent masses", + &ClusterHadronizationHandler::reshuffleMode_, 0, false, false); + static SwitchOption interfaceReshuffleModeGlobal + (interfaceReshuffleMode, + "Global", + "Global reshuffling on all final state partons", + 0); + static SwitchOption interfaceReshuffleModeColourConnected + (interfaceReshuffleMode, + "ColourConnected", + "Separate reshuffling for colour connected partons", + 1); + static Parameter interfaceMinVirtuality2 ("MinVirtuality2", "Minimum virtuality^2 of partons to use in calculating distances (unit [GeV2]).", &ClusterHadronizationHandler::_minVirtuality2, GeV2, 0.1*GeV2, ZERO, 10.0*GeV2,false,false,false); static Parameter interfaceMaxDisplacement ("MaxDisplacement", "Maximum displacement that is allowed for a particle (unit [millimeter]).", &ClusterHadronizationHandler::_maxDisplacement, mm, 1.0e-10*mm, 0.0*mm, 1.0e-9*mm,false,false,false); static Reference interfaceUnderlyingEventHandler ("UnderlyingEventHandler", "Pointer to the handler for the Underlying Event. " "Set to NULL to disable.", &ClusterHadronizationHandler::_underlyingEventHandler, false, false, true, true, false); static Switch interfaceReduceToTwoComponents ("ReduceToTwoComponents", "Whether or not to reduce three component baryon-number violating clusters to two components before cluster splitting or leave" " this till after the cluster splitting", &ClusterHadronizationHandler::_reduceToTwoComponents, true, false, false); static SwitchOption interfaceReduceToTwoComponentsYes (interfaceReduceToTwoComponents, "BeforeSplitting", "Reduce to two components", true); static SwitchOption interfaceReduceToTwoComponentsNo (interfaceReduceToTwoComponents, "AfterSplitting", "Treat as three components", false); } namespace { void extractChildren(tPPtr p, set & all) { if (p->children().empty()) return; for (PVector::const_iterator child = p->children().begin(); child != p->children().end(); ++child) { all.insert(*child); extractChildren(*child, all); } } } void ClusterHadronizationHandler:: handle(EventHandler & ch, const tPVector & tagged, const Hint &) { useMe(); currentHandler_ = this; - PVector currentlist(tagged.begin(),tagged.end()); + + PVector theList(tagged.begin(),tagged.end()); + + if ( reshuffle_ ) { + + vector reshufflelists; + + if (reshuffleMode_==0){ // global reshuffling + reshufflelists.push_back(theList); + } + else if (reshuffleMode_==1){// colour connected reshuffling + splitIntoColourSinglets(theList, reshufflelists); + } + + for (auto currentlist : reshufflelists){ + // get available energy and energy needed for constituent mass shells + LorentzMomentum totalQ; + Energy needQ = ZERO; + size_t nGluons = 0; // number of gluons for which a mass need be generated + for ( auto p : currentlist ) { + totalQ += p->momentum(); + if ( p->id() == ParticleID::g && gluonMassGenerator() ) { + ++nGluons; + continue; + } + needQ += p->dataPtr()->constituentMass(); + } + Energy Q = totalQ.m(); + if ( needQ > Q ) + throw Exception() << "cannot reshuffle to constituent mass shells" << Exception::eventerror; + + // generate gluon masses if needed + list gluonMasses; + if ( nGluons && gluonMassGenerator() ) + gluonMasses = gluonMassGenerator()->generateMany(nGluons,Q-needQ); + + // set masses for inidividual particles + vector masses; + for ( auto p : currentlist ) { + if ( p->id() == ParticleID::g && gluonMassGenerator() ) { + list::const_iterator it = gluonMasses.begin(); + advance(it,UseRandom::irnd(gluonMasses.size())); + masses.push_back(*it); + gluonMasses.erase(it); + } + else { + masses.push_back(p->dataPtr()->constituentMass()); + } + } + + // reshuffle to new masses + reshuffle(currentlist,masses); + + } + + } + // set the scale for coloured particles to just above the gluon mass squared // if less than this so they are classed as perturbative Energy2 Q02 = 1.01*sqr(getParticleData(ParticleID::g)->constituentMass()); - for(unsigned int ix=0;ixscale()scale(Q02); + for(unsigned int ix=0;ixscale()scale(Q02); } // split the gluons - _partonSplitter->split(currentlist); + _partonSplitter->split(theList); // form the clusters ClusterVector clusters = - _clusterFinder->formClusters(currentlist); + _clusterFinder->formClusters(theList); // reduce BV clusters to two components now if needed if(_reduceToTwoComponents) _clusterFinder->reduceToTwoComponents(clusters); // perform colour reconnection if needed and then // decay the clusters into one hadron bool lightOK = false; short tried = 0; const ClusterVector savedclusters = clusters; tPVector finalHadrons; // only needed for partonic decayer while (!lightOK && tried++ < 10) { // no colour reconnection with baryon-number-violating (BV) clusters ClusterVector CRclusters, BVclusters; CRclusters.reserve( clusters.size() ); BVclusters.reserve( clusters.size() ); for (size_t ic = 0; ic < clusters.size(); ++ic) { ClusterPtr cl = clusters.at(ic); bool hasClusterParent = false; for (unsigned int ix=0; ix < cl->parents().size(); ++ix) { if (cl->parents()[ix]->id() == ParticleID::Cluster) { hasClusterParent = true; break; } } if (cl->numComponents() > 2 || hasClusterParent) BVclusters.push_back(cl); else CRclusters.push_back(cl); } // colour reconnection _colourReconnector->rearrange(CRclusters); // tag new clusters as children of the partons to hadronize _setChildren(CRclusters); // forms diquarks _clusterFinder->reduceToTwoComponents(CRclusters); // recombine vectors of (possibly) reconnected and BV clusters clusters.clear(); clusters.insert( clusters.end(), CRclusters.begin(), CRclusters.end() ); clusters.insert( clusters.end(), BVclusters.begin(), BVclusters.end() ); // fission of heavy clusters // NB: during cluster fission, light hadrons might be produced straight away finalHadrons = _clusterFissioner->fission(clusters,isSoftUnderlyingEventON()); // if clusters not previously reduced to two components do it now if(!_reduceToTwoComponents) _clusterFinder->reduceToTwoComponents(clusters); lightOK = _lightClusterDecayer->decay(clusters,finalHadrons); // if the decay of the light clusters was not successful, undo the cluster // fission and decay steps and revert to the original state of the event // record if (!lightOK) { clusters = savedclusters; for_each(clusters.begin(), clusters.end(), std::mem_fn(&Particle::undecay)); } } if (!lightOK) { throw Exception( "CluHad::handle(): tried LightClusterDecayer 10 times!", Exception::eventerror); } // decay the remaining clusters _clusterDecayer->decay(clusters,finalHadrons); // ***************************************** // ***************************************** // ***************************************** bool finalStateCluster=false; StepPtr pstep = newStep(); set allDecendants; for (tPVector::const_iterator it = tagged.begin(); it != tagged.end(); ++it) { extractChildren(*it, allDecendants); } for(set::const_iterator it = allDecendants.begin(); it != allDecendants.end(); ++it) { // this is a workaround because the set sometimes // re-orders parents after their children if ((*it)->children().empty()){ // If there is a cluster in the final state throw an event error if((*it)->id()==81) { finalStateCluster=true; } pstep->addDecayProduct(*it); } else { pstep->addDecayProduct(*it); pstep->addIntermediate(*it); } } // For very small center of mass energies it might happen that baryonic clusters cannot decay into hadrons if (finalStateCluster){ throw Exception( "CluHad::Handle(): Cluster in the final state", Exception::eventerror); } // ***************************************** // ***************************************** // ***************************************** // soft underlying event if needed if (isSoftUnderlyingEventON()) { assert(_underlyingEventHandler); ch.performStep(_underlyingEventHandler,Hint::Default()); } } // Sets parent child relationship of all clusters with two components // Relationships for clusters with more than two components are set elsewhere in the Colour Reconnector void ClusterHadronizationHandler::_setChildren(const ClusterVector & clusters) const { // erase existing information about the partons' children tPVector partons; for ( const auto & cl : clusters ) { if ( cl->numComponents() > 2 ) continue; partons.push_back( cl->colParticle() ); partons.push_back( cl->antiColParticle() ); } // erase all previous information about parent child relationship for_each(partons.begin(), partons.end(), std::mem_fn(&Particle::undecay)); // give new parents to the clusters: their constituents for ( const auto & cl : clusters ) { if ( cl->numComponents() > 2 ) continue; cl->colParticle()->addChild(cl); cl->antiColParticle()->addChild(cl); } } + +void ClusterHadronizationHandler::splitIntoColourSinglets(PVector copylist, + vector& reshufflelists){ + + PVector currentlist; + bool gluonloop; + PPtr firstparticle, temp; + reshufflelists.clear(); + + while (copylist.size()>0){ + gluonloop=false; + currentlist.clear(); + + firstparticle=copylist.back(); + copylist.pop_back(); + + if (!firstparticle->coloured()){ + continue; //non-coloured particles are not included + } + + currentlist.push_back(firstparticle); + + //go up the anitColourLine and check if we are in a gluon loop + temp=firstparticle; + while( temp->hasAntiColour()){ + temp = temp->antiColourLine()->endParticle(); + if(temp==firstparticle){ + gluonloop=true; + break; + } + else{ + currentlist.push_back(temp); + copylist.erase(remove(copylist.begin(),copylist.end(), temp), copylist.end()); + } + } + + //if not a gluon loop, go up the ColourLine + if(!gluonloop){ + temp=firstparticle; + while( temp->hasColour()){ + temp=temp->colourLine()->startParticle(); + currentlist.push_back(temp); + copylist.erase(remove(copylist.begin(),copylist.end(), temp), copylist.end()); + } + } + + reshufflelists.push_back(currentlist); + } + +} diff --git a/Hadronization/ClusterHadronizationHandler.h b/Hadronization/ClusterHadronizationHandler.h --- a/Hadronization/ClusterHadronizationHandler.h +++ b/Hadronization/ClusterHadronizationHandler.h @@ -1,214 +1,247 @@ // -*- C++ -*- // // ClusterHadronizationHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_ClusterHadronizationHandler_H #define HERWIG_ClusterHadronizationHandler_H #include #include "PartonSplitter.h" #include "ClusterFinder.h" #include "ColourReconnector.h" #include "ClusterFissioner.h" #include "LightClusterDecayer.h" #include "ClusterDecayer.h" #include "ClusterHadronizationHandler.fh" +#include "Herwig/Utilities/Reshuffler.h" +#include "GluonMassGenerator.h" namespace Herwig { using namespace ThePEG; /** \ingroup Hadronization * \class ClusterHadronizationHandler * \brief Class that controls the cluster hadronization algorithm. * \author Philip Stephens // cerr << *ch.currentEvent() << '\n'; cerr << finalHadrons.size() << '\n'; cerr << "Finished hadronizing \n"; * \author Alberto Ribon * * This class is the main driver of the cluster hadronization: it is * responsible for the proper handling of all other specific collaborating * classes PartonSplitter, ClusterFinder, ColourReconnector, ClusterFissioner, * LightClusterDecayer, ClusterDecayer; * and for the storing of the produced particles in the Event record. * * @see PartonSplitter * @see ClusterFinder * @see ColourReconnector * @see ClusterFissioner * @see LightClusterDecayer * @see ClusterDecayer * @see Cluster * @see \ref ClusterHadronizationHandlerInterfaces "The interfaces" * defined for ClusterHadronizationHandler. */ -class ClusterHadronizationHandler: public HadronizationHandler { +class ClusterHadronizationHandler: + public HadronizationHandler, public Reshuffler { public: /** * The main method which manages the all cluster hadronization. * * This routine directs "traffic". It determines which function is called * and on which particles/clusters. This function also handles the * situation of vetos on the hadronization. */ virtual void handle(EventHandler & ch, const tPVector & tagged, const Hint & hint); /** * It returns minimum virtuality^2 of partons to use in calculating * distances. It is used both in the Showering and Hadronization. */ Energy2 minVirtuality2() const { return _minVirtuality2; } /** * It returns the maximum displacement that is allowed for a particle * (used to determine the position of a cluster with two components). */ Length maxDisplacement() const { return _maxDisplacement; } /** * It returns true/false according if the soft underlying model * is switched on/off. */ bool isSoftUnderlyingEventON() const { return _underlyingEventHandler; } /** * pointer to "this", the current HadronizationHandler. */ static const ClusterHadronizationHandler * currentHandler() { if(!currentHandler_){ cerr<< " \nCreating new ClusterHadronizationHandler without input from infiles."; cerr<< " \nWhen using for example the string model "; cerr<< " hadronic decays are still treated by the Cluster model\n"; currentHandler_=new ClusterHadronizationHandler();; } return currentHandler_; } + /** + * A pointer to a gluon mass generator for the reshuffling + */ + Ptr::tptr gluonMassGenerator() const { + return gluonMassGenerator_; + } + 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); //@} /** * Standard Init function used to initialize the interfaces. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} private: /** * Private and non-existent assignment operator. */ ClusterHadronizationHandler & operator=(const ClusterHadronizationHandler &) = delete; /** * This is a pointer to a Herwig::PartonSplitter object. */ PartonSplitterPtr _partonSplitter; /** * This is a pointer to a Herwig::ClusterFinder object. */ ClusterFinderPtr _clusterFinder; /** * This is a pointer to a Herwig::ColourReconnector object. */ ColourReconnectorPtr _colourReconnector; /** * This is a pointer to a Herwig::ClusterFissioner object. */ ClusterFissionerPtr _clusterFissioner; /** * This is a pointer to a Herwig::LightClusterDecayer object. */ LightClusterDecayerPtr _lightClusterDecayer; /** * This is a pointer to a Herwig::ClusterDecayer object. */ ClusterDecayerPtr _clusterDecayer; /** + * Perform reshuffling to constituent masses. + */ + bool reshuffle_ = false; + + /** + * Which type of reshuffling (global (default) or colour connected) is used + */ + int reshuffleMode_ = 0; + + /** + * A pointer to a gluon mass generator for the reshuffling + */ + Ptr::ptr gluonMassGenerator_; + + /** * The minimum virtuality^2 of partons to use in calculating * distances. */ Energy2 _minVirtuality2 = 0.1_GeV2; /** * The maximum displacement that is allowed for a particle * (used to determine the position of a cluster with two components). */ Length _maxDisplacement = 1.0e-10_mm; /** * The pointer to the Underlying Event handler. */ StepHdlPtr _underlyingEventHandler; /** * How to handle baryon-number clusters */ bool _reduceToTwoComponents = true; /** * Tag the constituents of the clusters as their parents */ void _setChildren(const ClusterVector & clusters) const; + + + /** + * Split the list of partons into colour connected sub-lists before reshuffling + */ + void splitIntoColourSinglets(PVector thelist, + vector& reshufflelists); + /** * pointer to "this", the current HadronizationHandler. */ static ClusterHadronizationHandler * currentHandler_; }; } #endif /* HERWIG_ClusterHadronizationHandler_H */ diff --git a/Hadronization/GluonMassGenerator.cc b/Hadronization/GluonMassGenerator.cc new file mode 100644 --- /dev/null +++ b/Hadronization/GluonMassGenerator.cc @@ -0,0 +1,63 @@ +// -*- C++ -*- +// +// GluonMassGenerator.cc is a part of Herwig - A multi-purpose Monte Carlo event generator +// Copyright (C) 2002-2019 The Herwig Collaboration +// +// Herwig is licenced under version 3 of the GPL, see COPYING for details. +// Please respect the MCnet academic guidelines, see GUIDELINES for details. +// +// +// This is the implementation of the non-inlined, non-templated member +// functions of the GluonMassGenerator class. +// + +#include "GluonMassGenerator.h" +#include "ThePEG/Interface/ClassDocumentation.h" +#include "ThePEG/EventRecord/Particle.h" +#include "ThePEG/Repository/UseRandom.h" +#include "ThePEG/Repository/EventGenerator.h" +#include "ThePEG/Utilities/DescribeClass.h" +#include "ThePEG/StandardModel/StandardModelBase.h" +#include "ClusterHadronizationHandler.h" + +#include "ThePEG/Persistency/PersistentOStream.h" +#include "ThePEG/Persistency/PersistentIStream.h" + +using namespace Herwig; + +GluonMassGenerator::GluonMassGenerator() {} + +GluonMassGenerator::~GluonMassGenerator() {} + +IBPtr GluonMassGenerator::clone() const { + return new_ptr(*this); +} + +IBPtr GluonMassGenerator::fullclone() const { + return new_ptr(*this); +} + +// If needed, insert default implementations of virtual function defined +// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). + + +void GluonMassGenerator::persistentOutput(PersistentOStream &) const {} + +void GluonMassGenerator::persistentInput(PersistentIStream &, int) {} + + +// *** Attention *** The following static variable is needed for the type +// description system in ThePEG. Please check that the template arguments +// are correct (the class and its base class), and that the constructor +// arguments are correct (the class name and the name of the dynamically +// loadable library where the class implementation can be found). +DescribeClass + describeHerwigGluonMassGenerator("Herwig::GluonMassGenerator", ""); + +void GluonMassGenerator::Init() { + + static ClassDocumentation documentation + ("Dynamic gluon mass generation"); + +} + diff --git a/Hadronization/GluonMassGenerator.h b/Hadronization/GluonMassGenerator.h new file mode 100644 --- /dev/null +++ b/Hadronization/GluonMassGenerator.h @@ -0,0 +1,171 @@ +// -*- C++ -*- +// +// GluonMassGenerator.h is a part of Herwig - A multi-purpose Monte Carlo event generator +// Copyright (C) 2002-2019 The Herwig Collaboration +// +// Herwig is licenced under version 3 of the GPL, see COPYING for details. +// Please respect the MCnet academic guidelines, see GUIDELINES for details. +// +#ifndef Herwig_GluonMassGenerator_H +#define Herwig_GluonMassGenerator_H +// +// This is the declaration of the GluonMassGenerator class. +// + +#include "ThePEG/Handlers/HandlerBase.h" +#include "ThePEG/EventRecord/Particle.h" + +namespace Herwig { + +using namespace ThePEG; + +/** + * \ingroup Hadronization + * \brief Dynamic gluon mass generator; the default returns a constant mass. + * + * @see \ref GluonMassGeneratorInterfaces "The interfaces" + * defined for GluonMassGenerator. + */ +class GluonMassGenerator: public HandlerBase { + +public: + + /** @name Standard constructors and destructors. */ + //@{ + /** + * The default constructor. + */ + GluonMassGenerator(); + + /** + * The destructor. + */ + virtual ~GluonMassGenerator(); + //@} + +public: + + /** + * Generate a single gluon mass with possible reference to a hard + * scale Q and up to a maximum value + */ + virtual Energy generate(Energy, Energy) const { + return generate(); + } + + /** + * Generate a single gluon mass with possible reference to a hard + * scale Q + */ + virtual Energy generate(Energy) const { + return generate(); + } + + /** + * Generate a single gluon mass without further constraints + */ + virtual Energy generate() const { + return getParticleData(ThePEG::ParticleID::g)->constituentMass(); + } + + /** + * Generate a list of n gluon masses, with a maximum available energy + */ + list generateMany(size_t n, Energy QMax) const { + list res; + Energy m0, mu, md, ms, mg, mgmax, summg; + + mu=getParticleData(ThePEG::ParticleID::u)->constituentMass(); + md=getParticleData(ThePEG::ParticleID::d)->constituentMass(); + ms=getParticleData(ThePEG::ParticleID::s)->constituentMass(); + + m0=md; + if(mu QMax - 2.0*m0*(n-k-1) ){ + repeat=true; + break; + } + } + } + + return res; + + } + +public: + + /** @name Functions used by the persistent I/O system. */ + //@{ + /** + * Function used to write out object persistently. + * @param os the persistent output stream written to. + */ + void persistentOutput(PersistentOStream & os) const; + + /** + * Function used to read in object persistently. + * @param is the persistent input stream read from. + * @param version the version number of the object when written. + */ + void persistentInput(PersistentIStream & is, int version); + //@} + + /** + * The standard Init function used to initialize the interfaces. + * Called exactly once for each class by the class description system + * before the main function starts or + * when this class is dynamically loaded. + */ + static void Init(); + +protected: + + /** @name Clone Methods. */ + //@{ + /** + * Make a simple clone of this object. + * @return a pointer to the new object. + */ + virtual IBPtr clone() const; + + /** Make a clone of this object, possibly modifying the cloned object + * to make it sane. + * @return a pointer to the new object. + */ + virtual IBPtr fullclone() const; + //@} + + +// If needed, insert declarations of virtual function defined in the +// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). + + +private: + + /** + * The assignment operator is private and must never be called. + * In fact, it should not even be implemented. + */ + GluonMassGenerator & operator=(const GluonMassGenerator &); + +}; + +} + +#endif /* Herwig_GluonMassGenerator_H */ diff --git a/Hadronization/Makefile.am b/Hadronization/Makefile.am --- a/Hadronization/Makefile.am +++ b/Hadronization/Makefile.am @@ -1,17 +1,18 @@ noinst_LTLIBRARIES = libHwHadronization.la libHwHadronization_la_SOURCES = \ CheckId.cc CheckId.h \ CluHadConfig.h \ Cluster.h Cluster.cc Cluster.fh \ ClusterDecayer.cc ClusterDecayer.h ClusterDecayer.fh \ ClusterFinder.cc ClusterFinder.h ClusterFinder.fh \ ClusterFissioner.cc ClusterFissioner.h ClusterFissioner.fh \ ClusterHadronizationHandler.cc ClusterHadronizationHandler.h \ ClusterHadronizationHandler.fh \ ColourReconnector.cc ColourReconnector.h ColourReconnector.fh\ +GluonMassGenerator.h GluonMassGenerator.cc \ HadronSelector.cc HadronSelector.h HadronSelector.fh\ Hw64Selector.cc Hw64Selector.h Hw64Selector.fh\ HwppSelector.cc HwppSelector.h HwppSelector.fh\ LightClusterDecayer.cc LightClusterDecayer.h LightClusterDecayer.fh \ PartonSplitter.cc PartonSplitter.h PartonSplitter.fh \ SpinHadronizer.h SpinHadronizer.cc diff --git a/Shower/Dipole/Utility/ConstituentReshuffler.cc b/Shower/Dipole/Utility/ConstituentReshuffler.cc --- a/Shower/Dipole/Utility/ConstituentReshuffler.cc +++ b/Shower/Dipole/Utility/ConstituentReshuffler.cc @@ -1,663 +1,622 @@ // -*- C++ -*- // // ConstituentReshuffler.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ConstituentReshuffler class. // #include #include "ConstituentReshuffler.h" #include "ThePEG/Interface/ClassDocumentation.h" #include #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "DipolePartonSplitter.h" #include "Herwig/Utilities/GSLBisection.h" #include "Herwig/Shower/Dipole/DipoleShowerHandler.h" #include "Herwig/Shower/ShowerHandler.h" using namespace Herwig; ConstituentReshuffler::ConstituentReshuffler() : HandlerBase() {} ConstituentReshuffler::~ConstituentReshuffler() {} IBPtr ConstituentReshuffler::clone() const { return new_ptr(*this); } IBPtr ConstituentReshuffler::fullclone() const { return new_ptr(*this); } -double ConstituentReshuffler::ReshuffleEquation::aUnit() { - return 1.; -} - -double ConstituentReshuffler::ReshuffleEquation::vUnit() { - return 1.; -} - -double ConstituentReshuffler::DecayReshuffleEquation::aUnit() { - return 1.; -} - -double ConstituentReshuffler::DecayReshuffleEquation::vUnit() { - return 1.; -} - - -double ConstituentReshuffler::ReshuffleEquation::operator() (double xi) const { - - double r = - w/GeV; - - for (PList::iterator p = p_begin; p != p_end; ++p) { - r += sqrt(sqr((**p).dataPtr()->constituentMass()) + - xi*xi*(sqr((**p).momentum().t())-sqr((**p).dataPtr()->mass()))) / GeV; - } - return r; - -} - -double ConstituentReshuffler::DecayReshuffleEquation::operator() (double xi) const { - double r = - w/GeV; - - for (PList::iterator pIt = p_begin; pIt != p_end; ++pIt) { - r += sqrt(sqr((**pIt).dataPtr()->constituentMass()) + - xi*xi*(sqr((**pIt).momentum().t())-sqr((**pIt).dataPtr()->mass()))) / GeV; - } - - for (PList::iterator rIt = r_begin; rIt != r_end; ++rIt) { - r += sqrt(sqr((**rIt).momentum().m()) + - xi*xi*(sqr((**rIt).momentum().t())-sqr((**rIt).momentum().m()))) / GeV; - } - - return r; - -} - void ConstituentReshuffler::reshuffle(PList& out, PPair& in, PList& intermediates, const bool decay, PList& decayPartons, PList& decayRecoilers) { assert(ShowerHandler::currentHandler()->retConstituentMasses()); if ( !decay ) { if (out.size() == 0) return; if (out.size() == 1) { PPtr recoiler; PPtr parton = out.front(); if (DipolePartonSplitter::colourConnected(parton,in.first) && DipolePartonSplitter::colourConnected(parton,in.second)) { if (UseRandom::rnd() < .5) recoiler = in.first; else recoiler = in.second; } else if (DipolePartonSplitter::colourConnected(parton,in.first)) { recoiler = in.first; } else if (DipolePartonSplitter::colourConnected(parton,in.second)) { recoiler = in.second; } else assert(false); assert(abs(recoiler->momentum().vect().perp2()/GeV2) < 1e-6); double sign = recoiler->momentum().z() < 0.*GeV ? -1. : 1.; Energy2 qperp2 = parton->momentum().perp2(); if (qperp2/GeV2 < Constants::epsilon) { // no emission off a 2 -> singlet process which // needed a single forced splitting: should never happen (?) assert(false); throw Veto(); } Energy2 m2 = sqr(parton->dataPtr()->constituentMass()); Energy abs_q = parton->momentum().vect().mag(); Energy qz = parton->momentum().z(); Energy abs_pz = recoiler->momentum().t(); assert(abs_pz > 0.*GeV); Energy xi_pz = sign*(2.*qperp2*abs_pz + m2*(abs_q + sign*qz))/(2.*qperp2); Energy x_qz = (2.*qperp2*qz + m2*(qz+sign*abs_q))/(2.*qperp2); Lorentz5Momentum recoiler_momentum (0.*GeV,0.*GeV,xi_pz,xi_pz < 0.*GeV ? - xi_pz : xi_pz); recoiler_momentum.rescaleMass(); Lorentz5Momentum parton_momentum (parton->momentum().x(),parton->momentum().y(),x_qz,sqrt(m2+qperp2+x_qz*x_qz)); parton_momentum.rescaleMass(); PPtr n_parton = new_ptr(Particle(parton->dataPtr())); n_parton->set5Momentum(parton_momentum); DipolePartonSplitter::change(parton,n_parton,false); out.pop_front(); intermediates.push_back(parton); out.push_back(n_parton); PPtr n_recoiler = new_ptr(Particle(recoiler->dataPtr())); n_recoiler->set5Momentum(recoiler_momentum); DipolePartonSplitter::change(recoiler,n_recoiler,true); intermediates.push_back(recoiler); if (recoiler == in.first) { in.first = n_recoiler; } if (recoiler == in.second) { in.second = n_recoiler; } return; } } Energy zero (0.*GeV); Lorentz5Momentum Q (zero,zero,zero,zero); for (PList::iterator p = out.begin(); p != out.end(); ++p) { Q += (**p).momentum(); } Boost beta = Q.findBoostToCM(); list mbackup; bool need_boost = (beta.mag2() > Constants::epsilon); if (need_boost) { for (PList::iterator p = out.begin(); p != out.end(); ++p) { Lorentz5Momentum mom = (**p).momentum(); mbackup.push_back(mom); (**p).set5Momentum(mom.boost(beta)); } } double xi; // Only partons if ( decayRecoilers.size()==0 ) { - ReshuffleEquation solve (Q.m(),out.begin(),out.end()); + list masses; + for ( auto p : out ) + masses.push_back(p->dataPtr()->constituentMass()); + ReshuffleEquation::const_iterator> + solve (Q.m(),out.begin(),out.end(), + masses.begin(),masses.end()); GSLBisection solver(1e-10,1e-8,10000); try { xi = solver.value(solve,0.0,1.1); } catch (GSLBisection::GSLerror) { throw DipoleShowerHandler::RedoShower(); } catch (GSLBisection::IntervalError) { throw DipoleShowerHandler::RedoShower(); } } // Partons and decaying recoilers else { DecayReshuffleEquation solve (Q.m(),decayPartons.begin(),decayPartons.end(),decayRecoilers.begin(),decayRecoilers.end()); GSLBisection solver(1e-10,1e-8,10000); try { xi = solver.value(solve,0.0,1.1); } catch (GSLBisection::GSLerror) { throw DipoleShowerHandler::RedoShower(); } catch (GSLBisection::IntervalError) { throw DipoleShowerHandler::RedoShower(); } } PList reshuffled; list::const_iterator backup_it; if (need_boost) backup_it = mbackup.begin(); // Reshuffling of non-decaying partons only if ( decayRecoilers.size()==0 ) { for (PList::iterator p = out.begin(); p != out.end(); ++p) { PPtr rp = new_ptr(Particle((**p).dataPtr())); DipolePartonSplitter::change(*p,rp,false); Lorentz5Momentum rm; rm = Lorentz5Momentum (xi*(**p).momentum().x(), xi*(**p).momentum().y(), xi*(**p).momentum().z(), sqrt(sqr((**p).dataPtr()->constituentMass()) + xi*xi*(sqr((**p).momentum().t())-sqr((**p).dataPtr()->mass())))); rm.rescaleMass(); if (need_boost) { (**p).set5Momentum(*backup_it); ++backup_it; rm.boost(-beta); } rp->set5Momentum(rm); intermediates.push_back(*p); reshuffled.push_back(rp); } } // For the case of a decay process with non-partonic recoilers else { assert ( decay ); for (PList::iterator p = out.begin(); p != out.end(); ++p) { // Flag to update spinInfo bool updateSpin = false; PPtr rp = new_ptr(Particle((**p).dataPtr())); DipolePartonSplitter::change(*p,rp,false); Lorentz5Momentum rm; // If the particle is a parton and not a recoiler if ( find( decayRecoilers.begin(), decayRecoilers.end(), *p ) == decayRecoilers.end() ) { rm = Lorentz5Momentum (xi*(**p).momentum().x(), xi*(**p).momentum().y(), xi*(**p).momentum().z(), sqrt(sqr((**p).dataPtr()->constituentMass()) + xi*xi*(sqr((**p).momentum().t())-sqr((**p).dataPtr()->mass())))); } // Otherwise the parton is a recoiler // and its invariant mass must be preserved else { if ( (*p)-> spinInfo() ) updateSpin = true; rm = Lorentz5Momentum (xi*(**p).momentum().x(), xi*(**p).momentum().y(), xi*(**p).momentum().z(), sqrt(sqr((**p).momentum().m()) + xi*xi*(sqr((**p).momentum().t())-sqr((**p).momentum().m())))); } rm.rescaleMass(); if (need_boost) { (**p).set5Momentum(*backup_it); ++backup_it; rm.boost(-beta); } rp->set5Momentum(rm); // Update SpinInfo if required if ( updateSpin ) updateSpinInfo(*p, rp); intermediates.push_back(*p); reshuffled.push_back(rp); } } out.clear(); out.splice(out.end(),reshuffled); } void ConstituentReshuffler::hardProcDecayReshuffle(PList& decaying, PList& eventOutgoing, PList& eventHard, PPair& eventIncoming, PList& eventIntermediates) { // Note, when this function is called, the particle pointers // in theDecays/decaying are those prior to the showering. // Here we find the newest pointers in the outgoing. // The update of the PPtrs in theDecays is done in DipoleShowerHandler::constituentReshuffle() // as this needs to be done if ConstituentReshuffling is switched off. //Make sure the shower should return constituent masses: assert(ShowerHandler::currentHandler()->retConstituentMasses()); // Find the outgoing decaying particles PList recoilers; for ( PList::iterator decIt = decaying.begin(); decIt != decaying.end(); ++decIt) { // First find the particles in the intermediates PList::iterator pos = find(eventIntermediates.begin(),eventIntermediates.end(), *decIt); // Colourless particle or coloured particle that did not radiate. if(pos==eventIntermediates.end()) { // Check that this is not a particle from a subsequent decay. // e.g. the W from a top decay from an LHE file. if ( find( eventHard.begin(), eventHard.end(), *decIt ) == eventHard.end() && find( eventOutgoing.begin(), eventOutgoing.end(), *decIt ) == eventOutgoing.end() ) continue; else recoilers.push_back( *decIt ); } // Coloured decaying particle that radiated else { PPtr unstable = *pos; while(!unstable->children().empty()) { unstable = unstable->children()[0]; } assert( find( eventOutgoing.begin(),eventOutgoing.end(), unstable ) != eventOutgoing.end() ); recoilers.push_back( unstable ); } } // Make a list of partons PList partons; for ( PList::iterator outPos = eventOutgoing.begin(); outPos != eventOutgoing.end(); ++outPos ) { if ( find (recoilers.begin(), recoilers.end(), *outPos ) == recoilers.end() ) { partons.push_back( *outPos ); } } // If no outgoing partons, do nothing if ( partons.size() == 0 ){ return; } // Otherwise reshuffling needs to be done. // If there is only one parton, attempt to reshuffle with // the incoming to be consistent with the reshuffle for a // hard process with no decays. else if ( partons.size() == 1 && ( DipolePartonSplitter::colourConnected(partons.front(),eventIncoming.first) || DipolePartonSplitter::colourConnected(partons.front(),eventIncoming.second) ) ) { // Erase the parton from the event outgoing eventOutgoing.erase( find( eventOutgoing.begin(), eventOutgoing.end(), partons.front() ) ); // Perform the reshuffle, this update the intermediates and the incoming reshuffle(partons, eventIncoming, eventIntermediates); // Update the outgoing eventOutgoing.push_back(partons.front()); return; } // If reshuffling amongst the incoming is not possible // or if we have multiple outgoing partons. else { // Create a complete list of the outgoing from the process PList out; // Make an empty list for storing the new intermediates PList intermediates; // Empty incoming particles pair PPair in; // A single parton which cannot be reshuffled // with the incoming. if ( partons.size() == 1 ) { // Populate the out for the reshuffling out.insert(out.end(),partons.begin(),partons.end()); out.insert(out.end(),recoilers.begin(),recoilers.end()); assert( out.size() > 1 ); // Perform the reshuffle with the temporary particle lists reshuffle(out, in, intermediates, true, partons, recoilers); } // If there is more than one parton, reshuffle only // amongst the partons else { assert(partons.size() > 1); // Populate the out for the reshuffling out.insert(out.end(),partons.begin(),partons.end()); assert( out.size() > 1 ); // Perform the reshuffle with the temporary particle lists reshuffle(out, in, intermediates, true); } // Update the dipole event record updateEvent(intermediates, eventIntermediates, out, eventOutgoing, eventHard ); return; } } void ConstituentReshuffler::decayReshuffle(PerturbativeProcessPtr& decayProc, PList& eventOutgoing, PList& eventHard, PList& eventIntermediates ) { // Separate particles into those to be assigned constituent masses // i.e. non-decaying coloured partons // and those which must only absorb recoil // i.e. non-coloured and decaying particles PList partons; PList recoilers; //Make sure the shower should return constituent masses: assert(ShowerHandler::currentHandler()->retConstituentMasses()); // Populate the particle lists from the outgoing of the decay process for( unsigned int ix = 0; ixoutgoing().size(); ++ix) { // Identify recoilers if ( !decayProc->outgoing()[ix].first->coloured() || ShowerHandler::currentHandler()->decaysInShower(decayProc->outgoing()[ix].first->id() ) ) recoilers.push_back(decayProc->outgoing()[ix].first); else partons.push_back(decayProc->outgoing()[ix].first); } // If there are no outgoing partons, then no reshuffling // needs to be done if ( partons.size() == 0 ) return; // Reshuffling needs to be done: else { // Create a complete list of the outgoing from the process PList out; // Make an empty list for storing the new intermediates PList intermediates; // Empty incoming particles pair PPair in; // SW - 15/06/2018, 31/01/2019 - Always include 'recoilers' in // reshuffling, regardless of the number of partons to be put on their // constituent mass shell. This is because reshuffling between 2 partons // frequently leads to a redoShower exception. This treatment is // consistent with the AO shower // Populate the out for the reshuffling out.insert(out.end(),partons.begin(),partons.end()); out.insert(out.end(),recoilers.begin(),recoilers.end()); assert( out.size() > 1 ); // Perform the reshuffle with the temporary particle lists reshuffle(out, in, intermediates, true, partons, recoilers); // Update the dipole event record and the decay process updateEvent(intermediates, eventIntermediates, out, eventOutgoing, eventHard, decayProc ); return; } } void ConstituentReshuffler::updateEvent( PList& intermediates, PList& eventIntermediates, #ifndef NDEBUG PList& out, #else PList&, #endif PList& eventOutgoing, PList& eventHard, PerturbativeProcessPtr decayProc ) { // Loop over the new intermediates following the reshuffling for (PList::iterator p = intermediates.begin(); p != intermediates.end(); ++p) { // Update the event record intermediates eventIntermediates.push_back(*p); // Identify the reshuffled particle assert( (*p)->children().size()==1 ); PPtr reshuffled = (*p)->children()[0]; assert( find(out.begin(), out.end(), reshuffled) != out.end() ); // Update the event record outgoing PList::iterator posOut = find(eventOutgoing.begin(), eventOutgoing.end(), *p); if ( posOut != eventOutgoing.end() ) { eventOutgoing.erase(posOut); eventOutgoing.push_back(reshuffled); } else { PList::iterator posHard = find(eventHard.begin(), eventHard.end(), *p); assert( posHard != eventHard.end() ); eventHard.erase(posHard); eventHard.push_back(reshuffled); } // Replace the particle in the the decay process outgoing if ( decayProc ) { vector >::iterator decayOutIt = decayProc->outgoing().end(); for ( decayOutIt = decayProc->outgoing().begin(); decayOutIt!= decayProc->outgoing().end(); ++decayOutIt ) { if ( decayOutIt->first == *p ){ break; } } assert( decayOutIt != decayProc->outgoing().end() ); decayOutIt->first = reshuffled; } } } void ConstituentReshuffler::updateSpinInfo( PPtr& oldPart, PPtr& newPart ) { const Lorentz5Momentum& oldMom = oldPart->momentum(); const Lorentz5Momentum& newMom = newPart->momentum(); // Rotation from old momentum to +ve z-axis LorentzRotation oldToZAxis; Axis axisOld(oldMom.vect().unit()); if( axisOld.perp2() > 1e-12 ) { double sinth(sqrt(1.-sqr(axisOld.z()))); oldToZAxis.rotate( -acos(axisOld.z()),Axis(-axisOld.y()/sinth,axisOld.x()/sinth,0.)); } // Rotation from new momentum to +ve z-axis LorentzRotation newToZAxis; Axis axisNew(newMom.vect().unit()); if( axisNew.perp2() > 1e-12 ) { double sinth(sqrt(1.-sqr(axisNew.z()))); newToZAxis.rotate( -acos(axisNew.z()),Axis(-axisNew.y()/sinth,axisNew.x()/sinth,0.)); } // Boost from old momentum to new momentum along z-axis Lorentz5Momentum momOldRotated = oldToZAxis*Lorentz5Momentum(oldMom); Lorentz5Momentum momNewRotated = newToZAxis*Lorentz5Momentum(newMom); Energy2 a = sqr(momOldRotated.z()) + sqr(momNewRotated.t()); Energy2 b = 2.*momOldRotated.t()*momOldRotated.z(); Energy2 c = sqr(momOldRotated.t()) - sqr(momNewRotated.t()); double beta; // The rotated momentum should always lie along the +ve z-axis if ( momOldRotated.z() > ZERO ) beta = (-b + sqrt(sqr(b)-4.*a*c)) / 2. / a; else beta = (-b - sqrt(sqr(b)-4.*a*c)) / 2. / a; LorentzRotation boostOldToNew(0., 0., beta); // Total transform LorentzRotation transform = (newToZAxis.inverse())*boostOldToNew*oldToZAxis; // Assign the same spin info to the old and new particles newPart->spinInfo(oldPart->spinInfo()); newPart->spinInfo()->transform(oldMom, transform); } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void ConstituentReshuffler::persistentOutput(PersistentOStream &) const { } void ConstituentReshuffler::persistentInput(PersistentIStream &, int) { } ClassDescription ConstituentReshuffler::initConstituentReshuffler; // Definition of the static class description member. void ConstituentReshuffler::Init() { static ClassDocumentation documentation ("The ConstituentReshuffler class implements reshuffling " "of partons on their nominal mass shell to their constituent " "mass shells."); } diff --git a/Shower/Dipole/Utility/ConstituentReshuffler.h b/Shower/Dipole/Utility/ConstituentReshuffler.h --- a/Shower/Dipole/Utility/ConstituentReshuffler.h +++ b/Shower/Dipole/Utility/ConstituentReshuffler.h @@ -1,299 +1,274 @@ // -*- C++ -*- // // ConstituentReshuffler.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_ConstituentReshuffler_H #define HERWIG_ConstituentReshuffler_H // // This is the declaration of the ConstituentReshuffler class. // #include "ThePEG/Handlers/HandlerBase.h" #include "ThePEG/Utilities/Exception.h" #include "Herwig/Shower/PerturbativeProcess.h" +#include "Herwig/Utilities/Reshuffler.h" namespace Herwig { using namespace ThePEG; /** * \ingroup DipoleShower * \author Simon Platzer, Stephen Webster * * \brief The ConstituentReshuffler class implements reshuffling * of partons on their nominal mass shell to their constituent * mass shells. * */ -class ConstituentReshuffler: public HandlerBase { +class ConstituentReshuffler: public HandlerBase, public Reshuffler { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ ConstituentReshuffler(); /** * The destructor. */ virtual ~ConstituentReshuffler(); //@} public: /** * Reshuffle the outgoing partons to constituent * masses. Optionally, incoming partons are given * to absorb recoils. Add the non-reshuffled partons * to the intermediates list. Throw ConstituentReshufflerProblem * if a numerical problem prevents the solution of * the reshuffling equation. */ void reshuffle(PList& out, PPair& in, PList& intermediates, const bool decay, PList& decayPartons, PList& decayRecoilers); /** * Reshuffle the outgoing partons to constituent * masses. Optionally, incoming partons are given * to absorb recoils. Add the non-reshuffled partons * to the intermediates list. Throw ConstituentReshufflerProblem * if a numerical problem prevents the solution of * the reshuffling equation. */ void reshuffle(PList& out, PPair& in, PList& intermediates, const bool decay=false) { PList decayPartons; PList decayRecoilers; reshuffle(out, in, intermediates, decay, decayPartons, decayRecoilers); } /** * Reshuffle the outgoing partons following the showering * of the initial hard interaction to constituent masses, * for the case of outgoing decaying particles. * Throw ConstituentReshufflerProblem * if a numerical problem prevents the solution of * the reshuffling equation. */ void hardProcDecayReshuffle(PList& decaying, PList& eventOutgoing, PList& eventHard, PPair& eventIncoming, PList& eventIntermediates) ; /** * Reshuffle the outgoing partons following the showering * of a particle decay to constituent masses. * Throw ConstituentReshufflerProblem * if a numerical problem prevents the solution of * the reshuffling equation. */ void decayReshuffle(PerturbativeProcessPtr& decayProc, PList& eventOutgoing, PList& eventHard, PList& eventIntermediates) ; /** * Update the dipole event record and, if appropriate, * the relevant decay process. **/ void updateEvent( PList& intermediates, PList& eventIntermediates, PList& out, PList& eventOutgoing, PList& eventHard, PerturbativeProcessPtr decayProc = PerturbativeProcessPtr() ) ; /** * Update the spinInfo of a particle following reshuffling * to take account of the change in momentum. * Used only for unstable particles that need to be dealt with. **/ void updateSpinInfo( PPtr& oldPart, PPtr& newPart ) ; protected: /** * The function object defining the equation - * to be solved. - */ - struct ReshuffleEquation { - - ReshuffleEquation (Energy q, - PList::iterator m_begin, - PList::iterator m_end) - : w(q), p_begin(m_begin), p_end(m_end) {} - - typedef double ArgType; - typedef double ValType; - - static double aUnit(); - static double vUnit(); - - double operator() (double xi) const; - - Energy w; - - PList::iterator p_begin; - PList::iterator p_end; - - }; - - /** - * The function object defining the equation * to be solved in the case of separate recoilers * TODO - refine the whole implementation of separate partons and recoilers */ struct DecayReshuffleEquation { DecayReshuffleEquation (Energy q, PList::iterator m_begin, PList::iterator m_end, PList::iterator n_begin, PList::iterator n_end) : w(q), p_begin(m_begin), p_end(m_end), r_begin(n_begin), r_end(n_end) {} typedef double ArgType; typedef double ValType; static double aUnit(); static double vUnit(); double operator() (double xi) const; Energy w; PList::iterator p_begin; PList::iterator p_end; PList::iterator r_begin; PList::iterator r_end; }; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** * The static object used to initialize the description of this class. * Indicates that this is a concrete class with persistent data. */ static ClassDescription initConstituentReshuffler; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ ConstituentReshuffler & operator=(const ConstituentReshuffler &) = delete; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of ConstituentReshuffler. */ template <> struct BaseClassTrait { /** Typedef of the first base class of ConstituentReshuffler. */ typedef HandlerBase NthBase; }; /** This template specialization informs ThePEG about the name of * the ConstituentReshuffler class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "Herwig::ConstituentReshuffler"; } /** * The name of a file containing the dynamic library where the class * ConstituentReshuffler is implemented. It may also include several, space-separated, * libraries if the class ConstituentReshuffler depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "HwDipoleShower.so"; } }; /** @endcond */ } #endif /* HERWIG_ConstituentReshuffler_H */ diff --git a/Shower/QTilde/SplittingFunctions/SudakovFormFactor.cc b/Shower/QTilde/SplittingFunctions/SudakovFormFactor.cc --- a/Shower/QTilde/SplittingFunctions/SudakovFormFactor.cc +++ b/Shower/QTilde/SplittingFunctions/SudakovFormFactor.cc @@ -1,1225 +1,1225 @@ // -*- C++ -*- // // SudakovFormFactor.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the SudakovFormFactor class. // #include "SudakovFormFactor.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "Herwig/Shower/QTilde/Kinematics/ShowerKinematics.h" #include "Herwig/Shower/QTilde/Base/ShowerParticle.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Shower/QTilde/QTildeShowerHandler.h" #include "Herwig/Shower/QTilde/Kinematics/FS_QTildeShowerKinematics1to2.h" #include "Herwig/Shower/QTilde/Kinematics/IS_QTildeShowerKinematics1to2.h" #include "Herwig/Shower/QTilde/Kinematics/Decay_QTildeShowerKinematics1to2.h" #include "Herwig/Shower/QTilde/Kinematics/KinematicHelpers.h" #include "SudakovCutOff.h" #include using std::array; using namespace Herwig; DescribeClass describeSudakovFormFactor ("Herwig::SudakovFormFactor",""); void SudakovFormFactor::persistentOutput(PersistentOStream & os) const { os << splittingFn_ << alpha_ << pdfmax_ << particles_ << pdffactor_ << cutoff_; } void SudakovFormFactor::persistentInput(PersistentIStream & is, int) { is >> splittingFn_ >> alpha_ >> pdfmax_ >> particles_ >> pdffactor_ >> cutoff_; } void SudakovFormFactor::Init() { static ClassDocumentation documentation ("The SudakovFormFactor class is the base class for the implementation of Sudakov" " form factors in Herwig"); static Reference interfaceSplittingFunction("SplittingFunction", "A reference to the SplittingFunction object", &Herwig::SudakovFormFactor::splittingFn_, false, false, true, false); static Reference interfaceAlpha("Alpha", "A reference to the Alpha object", &Herwig::SudakovFormFactor::alpha_, false, false, true, false); static Reference interfaceCutoff("Cutoff", "A reference to the SudakovCutOff object", &Herwig::SudakovFormFactor::cutoff_, false, false, true, false); static Parameter interfacePDFmax ("PDFmax", "Maximum value of PDF weight. ", &SudakovFormFactor::pdfmax_, 35.0, 1.0, 1000000.0, false, false, Interface::limited); static Switch interfacePDFFactor ("PDFFactor", "Include additional factors in the overestimate for the PDFs", &SudakovFormFactor::pdffactor_, 0, false, false); static SwitchOption interfacePDFFactorNo (interfacePDFFactor, "No", "Don't include any factors", 0); static SwitchOption interfacePDFFactorOverZ (interfacePDFFactor, "OverZ", "Include an additional factor of 1/z", 1); static SwitchOption interfacePDFFactorOverOneMinusZ (interfacePDFFactor, "OverOneMinusZ", "Include an additional factor of 1/(1-z)", 2); static SwitchOption interfacePDFFactorOverZOneMinusZ (interfacePDFFactor, "OverZOneMinusZ", "Include an additional factor of 1/z/(1-z)", 3); static SwitchOption interfacePDFFactorOverRootZ (interfacePDFFactor, "OverRootZ", "Include an additional factor of 1/sqrt(z)", 4); static SwitchOption interfacePDFFactorRootZ (interfacePDFFactor, "RootZ", "Include an additional factor of sqrt(z)", 5); } bool SudakovFormFactor::alphaSVeto(Energy2 pt2) const { double ratio=alphaSVetoRatio(pt2,1.); return UseRandom::rnd() > ratio; } double SudakovFormFactor::alphaSVetoRatio(Energy2 pt2, double factor) const { factor *= ShowerHandler::currentHandler()->renormalizationScaleFactor(); - return alpha_->ratio(pt2, factor); + return alpha_->showerRatio(pt2, factor); } bool SudakovFormFactor::PDFVeto(const Energy2 t, const double x, const tcPDPtr parton0, const tcPDPtr parton1, Ptr::transient_const_pointer beam) const { double ratio=PDFVetoRatio(t,x,parton0,parton1,beam,1.); return UseRandom::rnd() > ratio; } double SudakovFormFactor::PDFVetoRatio(const Energy2 t, const double x, const tcPDPtr parton0, const tcPDPtr parton1, Ptr::transient_const_pointer beam,double factor) const { assert(pdf_); Energy2 theScale = t * sqr(ShowerHandler::currentHandler()->factorizationScaleFactor()*factor); if (theScale < sqr(freeze_)) theScale = sqr(freeze_); const double newpdf=pdf_->xfx(beam,parton0,theScale,x/z()); if(newpdf<=0.) return 0.; const double oldpdf=pdf_->xfx(beam,parton1,theScale,x); if(oldpdf<=0.) return 1.; const double ratio = newpdf/oldpdf; double maxpdf = pdfmax_; switch (pdffactor_) { case 0: break; case 1: maxpdf /= z(); break; case 2: maxpdf /= 1.-z(); break; case 3: maxpdf /= (z()*(1.-z())); break; case 4: maxpdf /= sqrt(z()); break; case 5: maxpdf *= sqrt(z()); break; default : throw Exception() << "SudakovFormFactor::PDFVetoRatio invalid PDFfactor = " << pdffactor_ << Exception::runerror; } if (ratio > maxpdf) { generator()->log() << "PDFVeto warning: Ratio > " << name() << ":PDFmax (by a factor of " << ratio/maxpdf <<") for " << parton0->PDGName() << " to " << parton1->PDGName() << "\n"; } return ratio/maxpdf ; } void SudakovFormFactor::addSplitting(const IdList & in) { bool add=true; for(unsigned int ix=0;ix::iterator it=particles_.begin(); it!=particles_.end();++it) { if(it->size()==in.size()) { bool match=true; for(unsigned int iy=0;iy::iterator itemp=it; --itemp; particles_.erase(it); it = itemp; } } } } void SudakovFormFactor::guesstz(Energy2 t1,unsigned int iopt, const IdList &ids, double enhance,bool ident, double detune, Energy2 &t_main, double &z_main){ unsigned int pdfopt = iopt!=1 ? 0 : pdffactor_; double lower = splittingFn_->integOverP(zlimits_.first ,ids,pdfopt); double upper = splittingFn_->integOverP(zlimits_.second,ids,pdfopt); double c = 1./((upper - lower) - * alpha_->overestimateValue()/Constants::twopi*enhance*detune); + * alpha_->showerOverestimateValue()/Constants::twopi*enhance*detune); double r = UseRandom::rnd(); assert(iopt<=2); if(iopt==1) { c/=pdfmax_; //symmetry of FS gluon splitting if(ident) c*=0.5; } else if(iopt==2) c*=-1.; // guessing t if(iopt!=2 || c*log(r)invIntegOverP(lower + UseRandom::rnd() *(upper - lower),ids,pdfopt); } bool SudakovFormFactor::guessTimeLike(Energy2 &t,Energy2 tmin,double enhance, double detune) { Energy2 told = t; // calculate limits on z and if lower>upper return if(!computeTimeLikeLimits(t)) return false; // guess values of t and z guesstz(told,0,ids_,enhance,ids_[1]==ids_[2],detune,t,z_); // actual values for z-limits if(!computeTimeLikeLimits(t)) return false; if(tupper return if(!computeSpaceLikeLimits(t,x)) return false; // guess values of t and z guesstz(told,1,ids_,enhance,ids_[1]==ids_[2],detune,t,z_); // actual values for z-limits if(!computeSpaceLikeLimits(t,x)) return false; if(t zlimits_.second) return true; Energy2 m02 = (ids_[0]->id()!=ParticleID::g && ids_[0]->id()!=ParticleID::gamma) ? masssquared_[0] : Energy2(); Energy2 pt2 = QTildeKinematics::pT2_FSR(t,z(),m02,masssquared_[1],masssquared_[2], masssquared_[1],masssquared_[2]); // if pt2<0 veto if (pt2pT2min()) return true; // otherwise calculate pt and return pT_ = sqrt(pt2); return false; } ShoKinPtr SudakovFormFactor::generateNextTimeBranching(const Energy startingScale, const IdList &ids, const RhoDMatrix & rho, double enhance, double detuning) { // First reset the internal kinematics variables that can // have been eventually set in the previous call to the method. q_ = ZERO; z_ = 0.; phi_ = 0.; // perform initialization Energy2 tmax(sqr(startingScale)),tmin; initialize(ids,tmin); // check max > min if(tmax<=tmin) return ShoKinPtr(); // calculate next value of t using veto algorithm Energy2 t(tmax); // no shower variations to calculate if(ShowerHandler::currentHandler()->showerVariations().empty()){ // Without variations do the usual Veto algorithm // No need for more if-statements in this loop. do { if(!guessTimeLike(t,tmin,enhance,detuning)) break; } while(PSVeto(t) || SplittingFnVeto(z()*(1.-z())*t,ids,true,rho,detuning) || alphaSVeto(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t)); } else { bool alphaRew(true),PSRew(true),SplitRew(true); do { if(!guessTimeLike(t,tmin,enhance,detuning)) break; PSRew=PSVeto(t); if (PSRew) continue; SplitRew=SplittingFnVeto(z()*(1.-z())*t,ids,true,rho,detuning); alphaRew=alphaSVeto(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t); double factor=alphaSVetoRatio(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t,1.)* SplittingFnVetoRatio(z()*(1.-z())*t,ids,true,rho,detuning); tShowerHandlerPtr ch = ShowerHandler::currentHandler(); if( !(SplitRew || alphaRew) ) { //Emission q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV; if (q_ <= ZERO) break; } for ( map::const_iterator var = ch->showerVariations().begin(); var != ch->showerVariations().end(); ++var ) { if ( ( ch->firstInteraction() && var->second.firstInteraction ) || ( !ch->firstInteraction() && var->second.secondaryInteractions ) ) { double newfactor = alphaSVetoRatio(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t,var->second.renormalizationScaleFactor) * SplittingFnVetoRatio(z()*(1.-z())*t,ids,true,rho,detuning); double varied; if ( SplitRew || alphaRew ) { // No Emission varied = (1. - newfactor) / (1. - factor); } else { // Emission varied = newfactor / factor; } map::iterator wi = ch->currentWeights().find(var->first); if ( wi != ch->currentWeights().end() ) wi->second *= varied; else { assert(false); //ch->currentWeights()[var->first] = varied; } } } } while(PSRew || SplitRew || alphaRew); } q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV; if(q_ < ZERO) return ShoKinPtr(); // return the ShowerKinematics object return new_ptr(FS_QTildeShowerKinematics1to2(q_,z(),phi(),pT(),this)); } ShoKinPtr SudakovFormFactor:: generateNextSpaceBranching(const Energy startingQ, const IdList &ids, double x, const RhoDMatrix & rho, double enhance, Ptr::transient_const_pointer beam, double detuning) { // First reset the internal kinematics variables that can // have been eventually set in the previous call to the method. q_ = ZERO; z_ = 0.; phi_ = 0.; // perform the initialization Energy2 tmax(sqr(startingQ)),tmin; initialize(ids,tmin); // check max > min if(tmax<=tmin) return ShoKinPtr(); // calculate next value of t using veto algorithm Energy2 t(tmax),pt2(ZERO); // no shower variations if(ShowerHandler::currentHandler()->showerVariations().empty()){ // Without variations do the usual Veto algorithm // No need for more if-statements in this loop. do { if(!guessSpaceLike(t,tmin,x,enhance,detuning)) break; pt2 = QTildeKinematics::pT2_ISR(t,z(),masssquared_[2]); } while(pt2 < cutoff_->pT2min()|| z() > zlimits_.second|| SplittingFnVeto((1.-z())*t/z(),ids,false,rho,detuning)|| alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t)|| PDFVeto(t,x,ids[0],ids[1],beam)); } // shower variations else { bool alphaRew(true),PDFRew(true),ptRew(true),zRew(true),SplitRew(true); do { if(!guessSpaceLike(t,tmin,x,enhance,detuning)) break; pt2 = QTildeKinematics::pT2_ISR(t,z(),masssquared_[2]); ptRew=pt2 < cutoff_->pT2min(); zRew=z() > zlimits_.second; if (ptRew||zRew) continue; SplitRew=SplittingFnVeto((1.-z())*t/z(),ids,false,rho,detuning); alphaRew=alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t); PDFRew=PDFVeto(t,x,ids[0],ids[1],beam); double factor=PDFVetoRatio(t,x,ids[0],ids[1],beam,1.)* alphaSVetoRatio(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t,1.)* SplittingFnVetoRatio((1.-z())*t/z(),ids,false,rho,detuning); tShowerHandlerPtr ch = ShowerHandler::currentHandler(); if( !(PDFRew || SplitRew || alphaRew) ) { //Emission q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV; if (q_ <= ZERO) break; } for ( map::const_iterator var = ch->showerVariations().begin(); var != ch->showerVariations().end(); ++var ) { if ( ( ch->firstInteraction() && var->second.firstInteraction ) || ( !ch->firstInteraction() && var->second.secondaryInteractions ) ) { double newfactor = PDFVetoRatio(t,x,ids[0],ids[1],beam,var->second.factorizationScaleFactor)* alphaSVetoRatio(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t,var->second.renormalizationScaleFactor) *SplittingFnVetoRatio((1.-z())*t/z(),ids,false,rho,detuning); double varied; if( PDFRew || SplitRew || alphaRew) { // No Emission varied = (1. - newfactor) / (1. - factor); } else { // Emission varied = newfactor / factor; } map::iterator wi = ch->currentWeights().find(var->first); if ( wi != ch->currentWeights().end() ) wi->second *= varied; else { assert(false); //ch->currentWeights()[var->first] = varied; } } } } while( PDFRew || SplitRew || alphaRew); } if(t > ZERO && zlimits_.first < zlimits_.second) q_ = sqrt(t); else return ShoKinPtr(); pT_ = sqrt(pt2); // create the ShowerKinematics and return it return new_ptr(IS_QTildeShowerKinematics1to2(q_,z(),phi(),pT(),this)); } void SudakovFormFactor::initialize(const IdList & ids, Energy2 & tmin) { ids_=ids; tmin = 4.*cutoff_->pT2min(); masses_ = cutoff_->virtualMasses(ids); masssquared_.clear(); for(unsigned int ix=0;ix0) tmin=max(masssquared_[ix],tmin); } } ShoKinPtr SudakovFormFactor::generateNextDecayBranching(const Energy startingScale, const Energy stoppingScale, const Energy minmass, const IdList &ids, const RhoDMatrix & rho, double enhance, double detuning) { // First reset the internal kinematics variables that can // have been eventually set in the previous call to this method. q_ = Constants::MaxEnergy; z_ = 0.; phi_ = 0.; // perform initialisation Energy2 tmax(sqr(stoppingScale)),tmin; initialize(ids,tmin); tmin=sqr(startingScale); // check some branching possible if(tmax<=tmin) return ShoKinPtr(); // perform the evolution Energy2 t(tmin),pt2(-MeV2); do { if(!guessDecay(t,tmax,minmass,enhance,detuning)) break; pt2 = QTildeKinematics::pT2_Decay(t,z(),masssquared_[0],masssquared_[2]); } while(SplittingFnVeto((1.-z())*t/z(),ids,true,rho,detuning)|| alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t ) || pt2pT2min() || t*(1.-z())>masssquared_[0]-sqr(minmass)); if(t > ZERO) { q_ = sqrt(t); pT_ = sqrt(pt2); } else return ShoKinPtr(); phi_ = 0.; // create the ShowerKinematics object return new_ptr(Decay_QTildeShowerKinematics1to2(q_,z(),phi(),pT(),this)); } bool SudakovFormFactor::guessDecay(Energy2 &t,Energy2 tmax, Energy minmass, double enhance, double detune) { minmass = max(minmass,GeV); // previous scale Energy2 told = t; // overestimated limits on z if(tmaxpT2min()+ 0.25*sqr(masssquared_[2])/tm2)/tm +0.5*masssquared_[2]/tm2); if(zlimits_.secondpT2min()+ 0.25*sqr(masssquared_[2])/tm2)/tm +0.5*masssquared_[2]/tm2); if(t>tmax||zlimits_.secondpT2min(); // special case for gluon radiating if(ids_[0]->id()==ParticleID::g||ids_[0]->id()==ParticleID::gamma) { // no emission possible if(t<16.*(masssquared_[1]+pT2min)) { t=-1.*GeV2; return false; } // overestimate of the limits zlimits_.first = 0.5*(1.-sqrt(1.-4.*sqrt((masssquared_[1]+pT2min)/t))); zlimits_.second = 1.-zlimits_.first; } // special case for radiated particle is gluon else if(ids_[2]->id()==ParticleID::g||ids_[2]->id()==ParticleID::gamma) { zlimits_.first = sqrt((masssquared_[1]+pT2min)/t); zlimits_.second = 1.-sqrt((masssquared_[2]+pT2min)/t); } else if(ids_[1]->id()==ParticleID::g||ids_[1]->id()==ParticleID::gamma) { zlimits_.second = sqrt((masssquared_[2]+pT2min)/t); zlimits_.first = 1.-sqrt((masssquared_[1]+pT2min)/t); } else { zlimits_.first = (masssquared_[1]+pT2min)/t; zlimits_.second = 1.-(masssquared_[2]+pT2min)/t; } if(zlimits_.first>=zlimits_.second) { t=-1.*GeV2; return false; } return true; } bool SudakovFormFactor::computeSpaceLikeLimits(Energy2 & t, double x) { if (t < 1e-20 * GeV2) { t=-1.*GeV2; return false; } // compute the limits zlimits_.first = x; double yy = 1.+0.5*masssquared_[2]/t; zlimits_.second = yy - sqrt(sqr(yy)-1.+cutoff_->pT2min()/t); // return false if lower>upper if(zlimits_.second(particle.parents()[0]) : tShowerParticlePtr(); } else { mother = particle.children().size()==2 ? dynamic_ptr_cast(&particle) : tShowerParticlePtr(); } tShowerParticlePtr partner; while(mother) { tPPtr otherChild; if(forward) { for (unsigned int ix=0;ixchildren().size();++ix) { if(mother->children()[ix]!=child) { otherChild = mother->children()[ix]; break; } } } else { otherChild = mother->children()[1]; } tShowerParticlePtr other = dynamic_ptr_cast(otherChild); if((inter==ShowerInteraction::QCD && otherChild->dataPtr()->coloured()) || (inter==ShowerInteraction::QED && otherChild->dataPtr()->charged())) { partner = other; break; } if(forward && !other->isFinalState()) { partner = dynamic_ptr_cast(mother); break; } child = mother; if(forward) { mother = ! mother->parents().empty() ? dynamic_ptr_cast(mother->parents()[0]) : tShowerParticlePtr(); } else { if(mother->children()[0]->children().size()!=2) break; tShowerParticlePtr mtemp = dynamic_ptr_cast(mother->children()[0]); if(!mtemp) break; else mother=mtemp; } } if(!partner) { if(forward) { partner = dynamic_ptr_cast( child)->partner(); } else { if(mother) { tShowerParticlePtr parent; if(!mother->children().empty()) { parent = dynamic_ptr_cast(mother->children()[0]); } if(!parent) { parent = dynamic_ptr_cast(mother); } partner = parent->partner(); } else { partner = dynamic_ptr_cast(&particle)->partner(); } } } return partner; } pair softPhiMin(double phi0, double phi1, double A, double B, double C, double D) { double c01 = cos(phi0 - phi1); double s01 = sin(phi0 - phi1); double s012(sqr(s01)), c012(sqr(c01)); double A2(A*A), B2(B*B), C2(C*C), D2(D*D); if(abs(B/A)<1e-10 && abs(D/C)<1e-10) return make_pair(phi0,phi0+Constants::pi); double root = sqr(B2)*C2*D2*sqr(s012) + 2.*A*B2*B*C2*C*D*c01*s012 + 2.*A*B2*B*C*D2*D*c01*s012 + 4.*A2*B2*C2*D2*c012 - A2*B2*C2*D2*s012 - A2*B2*sqr(D2)*s012 - sqr(B2)*sqr(C2)*s012 - sqr(B2)*C2*D2*s012 - 4.*A2*A*B*C*D2*D*c01 - 4.*A*B2*B*C2*C*D*c01 + sqr(A2)*sqr(D2) + 2.*A2*B2*C2*D2 + sqr(B2)*sqr(C2); if(root<0.) return make_pair(phi0,phi0+Constants::pi); root = sqrt(root); double denom = (-2.*A*B*C*D*c01 + A2*D2 + B2*C2); double denom2 = (-B*C*c01 + A*D); double num = B2*C*D*s012; double y1 = B*s01*(-C*(num + root) + D*denom) / denom2; double y2 = B*s01*(-C*(num - root) + D*denom) / denom2; double x1 = -(num + root ); double x2 = -(num - root ); if(denom<0.) { y1*=-1.; y2*=-1.; x1*=-1.; x2*=-1.; } return make_pair(atan2(y1,x1) + phi0,atan2(y2,x2) + phi0); } } double SudakovFormFactor::generatePhiForward(ShowerParticle & particle, const IdList & ids, ShoKinPtr kinematics, const RhoDMatrix & rho) { // no correlations, return flat phi if(! dynamic_ptr_cast(ShowerHandler::currentHandler())->correlations()) return Constants::twopi*UseRandom::rnd(); // get the kinematic variables double z = kinematics->z(); Energy2 t = z*(1.-z)*sqr(kinematics->scale()); Energy pT = kinematics->pT(); // if soft correlations Energy2 pipj,pik; bool canBeSoft[2] = {ids[1]->id()==ParticleID::g || ids[1]->id()==ParticleID::gamma, ids[2]->id()==ParticleID::g || ids[2]->id()==ParticleID::gamma }; array pjk; array Ek; Energy Ei,Ej; Energy2 m12(ZERO),m22(ZERO); InvEnergy2 aziMax(ZERO); bool softAllowed = dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()&& (canBeSoft[0] || canBeSoft[1]); if(softAllowed) { // find the partner for the soft correlations tShowerParticlePtr partner=findCorrelationPartner(particle,true,splittingFn()->interactionType()); // remember we want the softer gluon bool swapOrder = !canBeSoft[1] || (canBeSoft[0] && canBeSoft[1] && z < 0.5); double zFact = !swapOrder ? (1.-z) : z; // compute the transforms to the shower reference frame // first the boost Lorentz5Momentum pVect = particle.showerBasis()->pVector(); Lorentz5Momentum nVect = particle.showerBasis()->nVector(); Boost beta_bb; if(particle.showerBasis()->frame()==ShowerBasis::BackToBack) { beta_bb = -(pVect + nVect).boostVector(); } else if(particle.showerBasis()->frame()==ShowerBasis::Rest) { beta_bb = -pVect.boostVector(); } else assert(false); pVect.boost(beta_bb); nVect.boost(beta_bb); Axis axis; if(particle.showerBasis()->frame()==ShowerBasis::BackToBack) { axis = pVect.vect().unit(); } else if(particle.showerBasis()->frame()==ShowerBasis::Rest) { axis = nVect.vect().unit(); } else assert(false); // and then the rotation LorentzRotation rot; if(axis.perp2()>0.) { double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); } else if(axis.z()<0.) { rot.rotate(Constants::pi,Axis(1.,0.,0.)); } rot.invert(); pVect *= rot; nVect *= rot; // shower parameters Energy2 pn = pVect*nVect, m2 = pVect.m2(); double alpha0 = particle.showerParameters().alpha; double beta0 = 0.5/alpha0/pn* (sqr(particle.dataPtr()->mass())-sqr(alpha0)*m2+sqr(particle.showerParameters().pt)); Lorentz5Momentum qperp0(particle.showerParameters().ptx, particle.showerParameters().pty,ZERO,ZERO); assert(partner); Lorentz5Momentum pj = partner->momentum(); pj.boost(beta_bb); pj *= rot; // compute the two phi independent dot products pik = 0.5*zFact*(sqr(alpha0)*m2 - sqr(particle.showerParameters().pt) + 2.*alpha0*beta0*pn ) +0.5*sqr(pT)/zFact; Energy2 dot1 = pj*pVect; Energy2 dot2 = pj*nVect; Energy2 dot3 = pj*qperp0; pipj = alpha0*dot1+beta0*dot2+dot3; // compute the constants for the phi dependent dot product pjk[0] = zFact*(alpha0*dot1+dot3-0.5*dot2/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) +0.5*sqr(pT)*dot2/pn/zFact/alpha0; pjk[1] = (pj.x() - dot2/alpha0/pn*qperp0.x())*pT; pjk[2] = (pj.y() - dot2/alpha0/pn*qperp0.y())*pT; m12 = sqr(particle.dataPtr()->mass()); m22 = sqr(partner->dataPtr()->mass()); if(swapOrder) { pjk[1] *= -1.; pjk[2] *= -1.; } Ek[0] = zFact*(alpha0*pVect.t()-0.5*nVect.t()/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) +0.5*sqr(pT)*nVect.t()/pn/zFact/alpha0; Ek[1] = -nVect.t()/alpha0/pn*qperp0.x()*pT; Ek[2] = -nVect.t()/alpha0/pn*qperp0.y()*pT; if(swapOrder) { Ek[1] *= -1.; Ek[2] *= -1.; } Energy mag2=sqrt(sqr(Ek[1])+sqr(Ek[2])); Ei = alpha0*pVect.t()+beta0*nVect.t(); Ej = pj.t(); double phi0 = atan2(-pjk[2],-pjk[1]); if(phi0<0.) phi0 += Constants::twopi; double phi1 = atan2(-Ek[2],-Ek[1]); if(phi1<0.) phi1 += Constants::twopi; double xi_min = pik/Ei/(Ek[0]+mag2), xi_max = pik/Ei/(Ek[0]-mag2), xi_ij = pipj/Ei/Ej; if(xi_min>xi_max) swap(xi_min,xi_max); if(xi_min>xi_ij) softAllowed = false; Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2])); if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { double A = (pipj*Ek[0]- Ej*pik)/Ej/sqr(Ej); double B = -sqrt(sqr(pipj)*(sqr(Ek[1])+sqr(Ek[2])))/Ej/sqr(Ej); double C = pjk[0]/sqr(Ej); double D = -sqrt(sqr(pjk[1])+sqr(pjk[2]))/sqr(Ej); pair minima = softPhiMin(phi0,phi1,A,B,C,D); aziMax = 0.5/pik/(Ek[0]-mag2)*(Ei-m12*(Ek[0]-mag2)/pik + max(Ej*(A+B*cos(minima.first -phi1))/(C+D*cos(minima.first -phi0)), Ej*(A+B*cos(minima.second-phi1))/(C+D*cos(minima.second-phi0)))); } else assert(false); } // if spin correlations vector > wgts; if(dynamic_ptr_cast(ShowerHandler::currentHandler())->spinCorrelations()) { // calculate the weights wgts = splittingFn()->generatePhiForward(z,t,ids,rho); } else { wgts = {{ {0, 1.} }}; } // generate the azimuthal angle double phi,wgt; static const Complex ii(0.,1.); unsigned int ntry(0); double phiMax(0.),wgtMax(0.); do { phi = Constants::twopi*UseRandom::rnd(); // first the spin correlations bit (gives 1 if correlations off) Complex spinWgt = 0.; for(unsigned int ix=0;ix1e-10) { generator()->log() << "Forward spin weight problem " << wgt << " " << wgt-1. << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n"; generator()->log() << "Weights \n"; for(unsigned int ix=0;ixlog() << wgts[ix].first << " " << wgts[ix].second << "\n"; } // soft correlations bit double aziWgt = 1.; if(softAllowed) { Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi); Energy Eg = Ek[0]+Ek[1]*cos(phi)+Ek[2]*sin(phi); if(pipj*Eg>pik*Ej) { if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { aziWgt = max(ZERO,0.5/pik/Eg*(Ei-m12*Eg/pik + (pipj*Eg - Ej*pik)/dot)/aziMax); } if(aziWgt-1.>1e-10||aziWgt<-1e-10) { generator()->log() << "Forward soft weight problem " << aziWgt << " " << aziWgt-1. << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n"; } } else { aziWgt = 0.; } } wgt *= aziWgt; if(wgt>wgtMax) { phiMax = phi; wgtMax = wgt; } ++ntry; } while(wgtlog() << "Too many tries to generate phi in forward evolution\n"; phi = phiMax; } // return the azimuthal angle return phi; } double SudakovFormFactor::generatePhiBackward(ShowerParticle & particle, const IdList & ids, ShoKinPtr kinematics, const RhoDMatrix & rho) { // no correlations, return flat phi if(! dynamic_ptr_cast(ShowerHandler::currentHandler())->correlations()) return Constants::twopi*UseRandom::rnd(); // get the kinematic variables double z = kinematics->z(); Energy2 t = (1.-z)*sqr(kinematics->scale())/z; Energy pT = kinematics->pT(); // if soft correlations bool softAllowed = dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations() && (ids[2]->id()==ParticleID::g || ids[2]->id()==ParticleID::gamma); Energy2 pipj,pik,m12(ZERO),m22(ZERO); array pjk; Energy Ei,Ej,Ek; InvEnergy2 aziMax(ZERO); if(softAllowed) { // find the partner for the soft correlations tShowerParticlePtr partner=findCorrelationPartner(particle,false,splittingFn()->interactionType()); double zFact = (1.-z); // compute the transforms to the shower reference frame // first the boost Lorentz5Momentum pVect = particle.showerBasis()->pVector(); Lorentz5Momentum nVect = particle.showerBasis()->nVector(); assert(particle.showerBasis()->frame()==ShowerBasis::BackToBack); Boost beta_bb = -(pVect + nVect).boostVector(); pVect.boost(beta_bb); nVect.boost(beta_bb); Axis axis = pVect.vect().unit(); // and then the rotation LorentzRotation rot; if(axis.perp2()>0.) { double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); } else if(axis.z()<0.) { rot.rotate(Constants::pi,Axis(1.,0.,0.)); } rot.invert(); pVect *= rot; nVect *= rot; // shower parameters Energy2 pn = pVect*nVect; Energy2 m2 = pVect.m2(); double alpha0 = particle.x(); double beta0 = -0.5/alpha0/pn*sqr(alpha0)*m2; Lorentz5Momentum pj = partner->momentum(); pj.boost(beta_bb); pj *= rot; double beta2 = 0.5*(1.-zFact)*(sqr(alpha0*zFact/(1.-zFact))*m2+sqr(pT))/alpha0/zFact/pn; // compute the two phi independent dot products Energy2 dot1 = pj*pVect; Energy2 dot2 = pj*nVect; pipj = alpha0*dot1+beta0*dot2; pik = alpha0*(alpha0*zFact/(1.-zFact)*m2+pn*(beta2+zFact/(1.-zFact)*beta0)); // compute the constants for the phi dependent dot product pjk[0] = alpha0*zFact/(1.-zFact)*dot1+beta2*dot2; pjk[1] = pj.x()*pT; pjk[2] = pj.y()*pT; m12 = ZERO; m22 = sqr(partner->dataPtr()->mass()); Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2])); if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { Ek = alpha0*zFact/(1.-zFact)*pVect.t()+beta2*nVect.t(); Ei = alpha0*pVect.t()+beta0*nVect.t(); Ej = pj.t(); if(pipj*Ek> Ej*pik) { aziMax = 0.5/pik/Ek*(Ei-m12*Ek/pik + (pipj*Ek- Ej*pik)/(pjk[0]-mag)); } else { aziMax = 0.5/pik/Ek*(Ei-m12*Ek/pik); } } else { assert(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==0); } } // if spin correlations vector > wgts; if(dynamic_ptr_cast(ShowerHandler::currentHandler())->spinCorrelations()) { // get the weights wgts = splittingFn()->generatePhiBackward(z,t,ids,rho); } else { wgts = {{ {0, 1.} }}; } // generate the azimuthal angle double phi,wgt; static const Complex ii(0.,1.); unsigned int ntry(0); double phiMax(0.),wgtMax(0.); do { phi = Constants::twopi*UseRandom::rnd(); Complex spinWgt = 0.; for(unsigned int ix=0;ix1e-10) { generator()->log() << "Backward weight problem " << wgt << " " << wgt-1. << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << z << " " << phi << "\n"; generator()->log() << "Weights \n"; for(unsigned int ix=0;ixlog() << wgts[ix].first << " " << wgts[ix].second << "\n"; } // soft correlations bit double aziWgt = 1.; if(softAllowed) { Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi); if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { aziWgt = max(ZERO,0.5/pik/Ek*(Ei-m12*Ek/pik + pipj*Ek/dot - Ej*pik/dot)/aziMax); } if(aziWgt-1.>1e-10||aziWgt<-1e-10) { generator()->log() << "Backward soft weight problem " << aziWgt << " " << aziWgt-1. << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n"; } } wgt *= aziWgt; if(wgt>wgtMax) { phiMax = phi; wgtMax = wgt; } ++ntry; } while(wgtlog() << "Too many tries to generate phi in backward evolution\n"; phi = phiMax; } // return the azimuthal angle return phi; } double SudakovFormFactor::generatePhiDecay(ShowerParticle & particle, const IdList & ids, ShoKinPtr kinematics, const RhoDMatrix &) { // only soft correlations in this case // no correlations, return flat phi if( !(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations() && (ids[2]->id()==ParticleID::g || ids[2]->id()==ParticleID::gamma ))) return Constants::twopi*UseRandom::rnd(); // get the kinematic variables double z = kinematics->z(); Energy pT = kinematics->pT(); // if soft correlations // find the partner for the soft correlations tShowerParticlePtr partner = findCorrelationPartner(particle,true,splittingFn()->interactionType()); double zFact(1.-z); // compute the transforms to the shower reference frame // first the boost Lorentz5Momentum pVect = particle.showerBasis()->pVector(); Lorentz5Momentum nVect = particle.showerBasis()->nVector(); assert(particle.showerBasis()->frame()==ShowerBasis::Rest); Boost beta_bb = -pVect.boostVector(); pVect.boost(beta_bb); nVect.boost(beta_bb); Axis axis = nVect.vect().unit(); // and then the rotation LorentzRotation rot; if(axis.perp2()>0.) { double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); } else if(axis.z()<0.) { rot.rotate(Constants::pi,Axis(1.,0.,0.)); } rot.invert(); pVect *= rot; nVect *= rot; // shower parameters Energy2 pn = pVect*nVect; Energy2 m2 = pVect.m2(); double alpha0 = particle.showerParameters().alpha; double beta0 = 0.5/alpha0/pn* (sqr(particle.dataPtr()->mass())-sqr(alpha0)*m2+sqr(particle.showerParameters().pt)); Lorentz5Momentum qperp0(particle.showerParameters().ptx, particle.showerParameters().pty,ZERO,ZERO); Lorentz5Momentum pj = partner->momentum(); pj.boost(beta_bb); pj *= rot; // compute the two phi independent dot products Energy2 pik = 0.5*zFact*(sqr(alpha0)*m2 - sqr(particle.showerParameters().pt) + 2.*alpha0*beta0*pn ) +0.5*sqr(pT)/zFact; Energy2 dot1 = pj*pVect; Energy2 dot2 = pj*nVect; Energy2 dot3 = pj*qperp0; Energy2 pipj = alpha0*dot1+beta0*dot2+dot3; // compute the constants for the phi dependent dot product array pjk; pjk[0] = zFact*(alpha0*dot1+dot3-0.5*dot2/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) +0.5*sqr(pT)*dot2/pn/zFact/alpha0; pjk[1] = (pj.x() - dot2/alpha0/pn*qperp0.x())*pT; pjk[2] = (pj.y() - dot2/alpha0/pn*qperp0.y())*pT; Energy2 m12 = sqr(particle.dataPtr()->mass()); Energy2 m22 = sqr(partner->dataPtr()->mass()); Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2])); InvEnergy2 aziMax; array Ek; Energy Ei,Ej; if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { Ek[0] = zFact*(alpha0*pVect.t()+-0.5*nVect.t()/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) +0.5*sqr(pT)*nVect.t()/pn/zFact/alpha0; Ek[1] = -nVect.t()/alpha0/pn*qperp0.x()*pT; Ek[2] = -nVect.t()/alpha0/pn*qperp0.y()*pT; Energy mag2=sqrt(sqr(Ek[1])+sqr(Ek[2])); Ei = alpha0*pVect.t()+beta0*nVect.t(); Ej = pj.t(); aziMax = 0.5/pik/(Ek[0]-mag2)*(Ei-m12*(Ek[0]-mag2)/pik + pipj*(Ek[0]+mag2)/(pjk[0]-mag) - Ej*pik/(pjk[0]-mag) ); } else assert(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==0); // generate the azimuthal angle double phi,wgt(0.); unsigned int ntry(0); double phiMax(0.),wgtMax(0.); do { phi = Constants::twopi*UseRandom::rnd(); Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi); if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { wgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { if(qperp0.m2()==ZERO) { wgt = 1.; } else { Energy Eg = Ek[0]+Ek[1]*cos(phi)+Ek[2]*sin(phi); wgt = max(ZERO,0.5/pik/Eg*(Ei-m12*Eg/pik + (pipj*Eg - Ej*pik)/dot)/aziMax); } } if(wgt-1.>1e-10||wgt<-1e-10) { generator()->log() << "Decay soft weight problem " << wgt << " " << wgt-1. << " " << ids[0]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << phi << "\n"; } if(wgt>wgtMax) { phiMax = phi; wgtMax = wgt; } ++ntry; } while(wgtlog() << "Too many tries to generate phi\n"; } // return the azimuthal angle return phi; } Energy SudakovFormFactor::calculateScale(double zin, Energy pt, IdList ids, unsigned int iopt) { Energy2 tmin; initialize(ids,tmin); // final-state branching if(iopt==0) { Energy2 scale=(sqr(pt)+masssquared_[1]*(1.-zin)+masssquared_[2]*zin); if(ids[0]->id()!=ParticleID::g) scale -= zin*(1.-zin)*masssquared_[0]; scale /= sqr(zin*(1-zin)); return scale<=ZERO ? sqrt(tmin) : sqrt(scale); } else if(iopt==1) { Energy2 scale=(sqr(pt)+zin*masssquared_[2])/sqr(1.-zin); return scale<=ZERO ? sqrt(tmin) : sqrt(scale); } else if(iopt==2) { Energy2 scale = (sqr(pt)+zin*masssquared_[2])/sqr(1.-zin)+masssquared_[0]; return scale<=ZERO ? sqrt(tmin) : sqrt(scale); } else { throw Exception() << "Unknown option in SudakovFormFactor::calculateScale() " << "iopt = " << iopt << Exception::runerror; } } diff --git a/Shower/ShowerAlpha.h b/Shower/ShowerAlpha.h --- a/Shower/ShowerAlpha.h +++ b/Shower/ShowerAlpha.h @@ -1,150 +1,179 @@ // -*- C++ -*- // // ShowerAlpha.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_ShowerAlpha_H #define HERWIG_ShowerAlpha_H // // This is the declaration of the ShowerAlpha class. // #include "ThePEG/Interface/Interfaced.h" #include "ShowerAlpha.fh" namespace Herwig { using namespace ThePEG; /** \ingroup Shower * * This class is the abstract class from which all types of running couplings * used in the Showering derive from. * The main purpose of this class, and the ones that derive from it, is * to allow systematic uncertainties for the initial-state radiation and, * independently, the final-state radiation effects, to be evaluated. * * This is achieved by allowing a multiplicative factor, * which is 1.0 for the "central value", * for the scale argument, \f$\mu^2\f$, of the running coupling. This * factor, \f$f\f$ is given by the scaleFactor() member and the coupling * returned by the value() member is such that * \f[\alpha(\mu^2)\to \alpha(f\times\mu^2).\f] * This scale factor is a parameter which is settable by the user, via the * interface. * Although, of course, it is not clear my how much we should scale * in order to get a \f$1\sigma\f$ systematic error (but factors: * 1/2 and 2 are quite common), this method provides a double side error, * and it appears more sensible than the rough and one-sided evaluation * obtained * via turning off the I.S.R. and/or F.S.R. (possibilities which are, * anyway, provided by Herwig). * * @see \ref ShowerAlphaInterfaces "The interfaces" * defined for ShowerAlpha. */ class ShowerAlpha: public Interfaced { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ ShowerAlpha() : _scaleFactor( 1.0 ) {} //@} public: /** * Methods to return the coupling and the scaleFactor */ //@{ /** * Pure virtual method that is supposed to return the * running alpha value evaluated at the input scale. * @param scale The scale * @return The coupling */ virtual double value(const Energy2 scale) const = 0; /** * Virtual method, which * should be overridden in a derived class to provide an * overestimate approximation of the alpha value. */ virtual double overestimateValue() const = 0; /** * Virtual method which returns the ratio of the running alpha * value at the input scale to the overestimated value. * @param scale The scale * @return The ratio */ virtual double ratio(const Energy2 scale,double factor=1.) const = 0; /** * It returns the factor that multiplies the * scale argument, \f$\mu^2\f$, of the running coupling. * This is supposed to be 1.0 in normal conditions (central values) * whereas different values can be useful for systematics evaluation * for Initial State radiation or Final State radiation effects. */ double scaleFactor() const {return _scaleFactor;} /** + * Virtual method that is supposed to return the + * running alpha value evaluated at the input scale. + * @param scale The scale + * @return The coupling + */ + virtual inline double showerValue(const Energy2 scale) const { + return value(scale); + } + + /** + * Virtual method, which + * should be overridden in a derived class to provide an + * overestimate approximation of the alpha value. + */ + virtual inline double showerOverestimateValue() const { + return overestimateValue(); + } + + /** + * Virtual method which returns the ratio of the running alpha + * value at the input scale to the overestimated value. + * @param scale The scale + * @return The ratio + */ + virtual inline double showerRatio(const Energy2 scale,double factor=1.) const { + return ratio(scale,factor); + } + + /** * Initialize this coupling. */ virtual void initialize () {} //@} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ ShowerAlpha & operator=(const ShowerAlpha &) = delete; private: /** * The scale factor */ double _scaleFactor; }; } #endif /* HERWIG_ShowerAlpha_H */ diff --git a/Shower/ShowerHandler.cc b/Shower/ShowerHandler.cc --- a/Shower/ShowerHandler.cc +++ b/Shower/ShowerHandler.cc @@ -1,1132 +1,1147 @@ // -*- C++ -*- // // ShowerHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ShowerHandler class. // #include "ShowerHandler.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Command.h" #include "ThePEG/PDF/PartonExtractor.h" #include "ThePEG/PDF/PartonBinInstance.h" #include "Herwig/PDT/StandardMatchers.h" #include "ThePEG/Cuts/Cuts.h" #include "ThePEG/Handlers/StandardXComb.h" #include "ThePEG/Utilities/Throw.h" #include "ThePEG/Utilities/StringUtils.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Repository/EventGenerator.h" #include "Herwig/Utilities/EnumParticles.h" #include "Herwig/PDF/MPIPDF.h" #include "Herwig/PDF/MinBiasPDF.h" #include "ThePEG/Handlers/EventHandler.h" #include "Herwig/PDF/HwRemDecayer.h" #include #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Decay/DecayIntegrator.h" #include "Herwig/Decay/PhaseSpaceMode.h" using namespace Herwig; DescribeClass describeShowerHandler ("Herwig::ShowerHandler","HwShower.so"); ShowerHandler::~ShowerHandler() {} tShowerHandlerPtr ShowerHandler::currentHandler_ = tShowerHandlerPtr(); void ShowerHandler::doinit() { CascadeHandler::doinit(); // copy particles to decay before showering from input vector to the // set used in the simulation if ( particlesDecayInShower_.empty() ) { for(unsigned int ix=0;ixunrestrictedPhasespace() && restrictPhasespace() ) { generator()->log() << "ShowerApproximation warning: The scale profile chosen requires an unrestricted phase space,\n" << "however, the phase space was set to be restricted. Will switch to unrestricted phase space.\n" << flush; restrictPhasespace_ = false; } } } IBPtr ShowerHandler::clone() const { return new_ptr(*this); } IBPtr ShowerHandler::fullclone() const { return new_ptr(*this); } ShowerHandler::ShowerHandler() : maxtry_(10),maxtryMPI_(10),maxtryDP_(10),maxtryDecay_(100), factorizationScaleFactor_(1.0), renormalizationScaleFactor_(1.0), hardScaleFactor_(1.0), restrictPhasespace_(true), maxPtIsMuF_(false), spinOpt_(1), pdfFreezingScale_(2.5*GeV), doFSR_(true), doISR_(true), splitHardProcess_(true), includeSpaceTime_(false), vMin_(0.1*GeV2), reweight_(1.0) { inputparticlesDecayInShower_.push_back( 6 ); // top inputparticlesDecayInShower_.push_back( 23 ); // Z0 inputparticlesDecayInShower_.push_back( 24 ); // W+/- inputparticlesDecayInShower_.push_back( 25 ); // h0 } void ShowerHandler::doinitrun(){ CascadeHandler::doinitrun(); //can't use isMPIOn here, because the EventHandler is not set at that stage if(MPIHandler_) { MPIHandler_->initialize(); if(MPIHandler_->softInt()) remDec_->initSoftInteractions(MPIHandler_->Ptmin(), MPIHandler_->beta()); } } void ShowerHandler::dofinish() { CascadeHandler::dofinish(); if(MPIHandler_) MPIHandler_->finalize(); } void ShowerHandler::persistentOutput(PersistentOStream & os) const { os << remDec_ << ounit(pdfFreezingScale_,GeV) << maxtry_ << maxtryMPI_ << maxtryDP_ << maxtryDecay_ << inputparticlesDecayInShower_ << particlesDecayInShower_ << MPIHandler_ << PDFA_ << PDFB_ << PDFARemnant_ << PDFBRemnant_ << includeSpaceTime_ << ounit(vMin_,GeV2) << factorizationScaleFactor_ << renormalizationScaleFactor_ << hardScaleFactor_ << restrictPhasespace_ << maxPtIsMuF_ << hardScaleProfile_ << showerVariations_ << doFSR_ << doISR_ << splitHardProcess_ - << spinOpt_ << useConstituentMasses_; + << spinOpt_ << useConstituentMasses_ << tagIntermediates_; } void ShowerHandler::persistentInput(PersistentIStream & is, int) { is >> remDec_ >> iunit(pdfFreezingScale_,GeV) >> maxtry_ >> maxtryMPI_ >> maxtryDP_ >> maxtryDecay_ >> inputparticlesDecayInShower_ >> particlesDecayInShower_ >> MPIHandler_ >> PDFA_ >> PDFB_ >> PDFARemnant_ >> PDFBRemnant_ >> includeSpaceTime_ >> iunit(vMin_,GeV2) >> factorizationScaleFactor_ >> renormalizationScaleFactor_ >> hardScaleFactor_ >> restrictPhasespace_ >> maxPtIsMuF_ >> hardScaleProfile_ >> showerVariations_ >> doFSR_ >> doISR_ >> splitHardProcess_ - >> spinOpt_ >> useConstituentMasses_; + >> spinOpt_ >> useConstituentMasses_ >> tagIntermediates_; } void ShowerHandler::Init() { static ClassDocumentation documentation ("Main driver class for the showering."); static Reference interfaceRemDecayer("RemDecayer", "A reference to the Remnant Decayer object", &Herwig::ShowerHandler::remDec_, false, false, true, false); static Parameter interfacePDFFreezingScale ("PDFFreezingScale", "The PDF freezing scale", &ShowerHandler::pdfFreezingScale_, GeV, 2.5*GeV, 2.0*GeV, 10.0*GeV, false, false, Interface::limited); static Parameter interfaceMaxTry ("MaxTry", "The maximum number of attempts for the main showering loop", &ShowerHandler::maxtry_, 10, 1, 100, false, false, Interface::limited); static Parameter interfaceMaxTryMPI ("MaxTryMPI", "The maximum number of regeneration attempts for an additional scattering", &ShowerHandler::maxtryMPI_, 10, 0, 100, false, false, Interface::limited); static Parameter interfaceMaxTryDP ("MaxTryDP", "The maximum number of regeneration attempts for an additional hard scattering", &ShowerHandler::maxtryDP_, 10, 0, 100, false, false, Interface::limited); static ParVector interfaceDecayInShower ("DecayInShower", "PDG codes of the particles to be decayed in the shower", &ShowerHandler::inputparticlesDecayInShower_, -1, 0l, -10000000l, 10000000l, false, false, Interface::limited); static Reference interfaceMPIHandler ("MPIHandler", "The object that administers all additional scatterings.", &ShowerHandler::MPIHandler_, false, false, true, true); static Reference interfacePDFA ("PDFA", "The PDF for beam particle A. Overrides the particle's own PDF setting." "By default used for both the shower and forced splitting in the remnant", &ShowerHandler::PDFA_, false, false, true, true, false); static Reference interfacePDFB ("PDFB", "The PDF for beam particle B. Overrides the particle's own PDF setting." "By default used for both the shower and forced splitting in the remnant", &ShowerHandler::PDFB_, false, false, true, true, false); static Reference interfacePDFARemnant ("PDFARemnant", "The PDF for beam particle A used to generate forced splittings of the remnant." " This overrides both the particle's own PDF setting and the value set by PDFA if used.", &ShowerHandler::PDFARemnant_, false, false, true, true, false); static Reference interfacePDFBRemnant ("PDFBRemnant", "The PDF for beam particle B used to generate forced splittings of the remnant." " This overrides both the particle's own PDF setting and the value set by PDFB if used.", &ShowerHandler::PDFBRemnant_, false, false, true, true, false); static Switch interfaceIncludeSpaceTime ("IncludeSpaceTime", "Whether to include the model for the calculation of space-time distances", &ShowerHandler::includeSpaceTime_, false, false, false); static SwitchOption interfaceIncludeSpaceTimeYes (interfaceIncludeSpaceTime, "Yes", "Include the model", true); static SwitchOption interfaceIncludeSpaceTimeNo (interfaceIncludeSpaceTime, "No", "Only include the displacement from the particle-s lifetime for decaying particles", false); static Parameter interfaceMinimumVirtuality ("MinimumVirtuality", "The minimum virtuality for the space-time model", &ShowerHandler::vMin_, GeV2, 0.1*GeV2, 0.0*GeV2, 1000.0*GeV2, false, false, Interface::limited); static Parameter interfaceFactorizationScaleFactor ("FactorizationScaleFactor", "The factorization scale factor.", &ShowerHandler::factorizationScaleFactor_, 1.0, 0.0, 0, false, false, Interface::lowerlim); static Parameter interfaceRenormalizationScaleFactor ("RenormalizationScaleFactor", "The renormalization scale factor.", &ShowerHandler::renormalizationScaleFactor_, 1.0, 0.0, 0, false, false, Interface::lowerlim); static Parameter interfaceHardScaleFactor ("HardScaleFactor", "The hard scale factor.", &ShowerHandler::hardScaleFactor_, 1.0, 0.0, 0, false, false, Interface::lowerlim); static Parameter interfaceMaxTryDecay ("MaxTryDecay", "The maximum number of attempts to generate a decay", &ShowerHandler::maxtryDecay_, 200, 10, 0, false, false, Interface::lowerlim); static Reference interfaceHardScaleProfile ("HardScaleProfile", "The hard scale profile to use.", &ShowerHandler::hardScaleProfile_, false, false, true, true, false); static Switch interfaceMaxPtIsMuF ("MaxPtIsMuF", "", &ShowerHandler::maxPtIsMuF_, false, false, false); static SwitchOption interfaceMaxPtIsMuFYes (interfaceMaxPtIsMuF, "Yes", "", true); static SwitchOption interfaceMaxPtIsMuFNo (interfaceMaxPtIsMuF, "No", "", false); static Switch interfaceRestrictPhasespace ("RestrictPhasespace", "Switch on or off phasespace restrictions", &ShowerHandler::restrictPhasespace_, true, false, false); static SwitchOption interfaceRestrictPhasespaceYes (interfaceRestrictPhasespace, "Yes", "Perform phasespace restrictions", true); static SwitchOption interfaceRestrictPhasespaceNo (interfaceRestrictPhasespace, "No", "Do not perform phasespace restrictions", false); static Command interfaceAddVariation ("AddVariation", "Add a shower variation.", &ShowerHandler::doAddVariation, false); static Switch interfaceDoFSR ("DoFSR", "Switch on or off final state radiation.", &ShowerHandler::doFSR_, true, false, false); static SwitchOption interfaceDoFSRYes (interfaceDoFSR, "Yes", "Switch on final state radiation.", true); static SwitchOption interfaceDoFSRNo (interfaceDoFSR, "No", "Switch off final state radiation.", false); static Switch interfaceDoISR ("DoISR", "Switch on or off initial state radiation.", &ShowerHandler::doISR_, true, false, false); static SwitchOption interfaceDoISRYes (interfaceDoISR, "Yes", "Switch on initial state radiation.", true); static SwitchOption interfaceDoISRNo (interfaceDoISR, "No", "Switch off initial state radiation.", false); static Switch interfaceSplitHardProcess ("SplitHardProcess", "Whether or not to try and split the hard process into production and decay processes", &ShowerHandler::splitHardProcess_, true, false, false); static SwitchOption interfaceSplitHardProcessYes (interfaceSplitHardProcess, "Yes", "Split the hard process", true); static SwitchOption interfaceSplitHardProcessNo (interfaceSplitHardProcess, "No", "Don't split the hard process", false); static Switch interfaceSpinCorrelations ("SpinCorrelations", "Treatment of spin correlations in the parton shower", &ShowerHandler::spinOpt_, 1, false, false); static SwitchOption interfaceSpinCorrelationsNo (interfaceSpinCorrelations, "No", "No spin correlations", 0); static SwitchOption interfaceSpinCorrelationsSpin (interfaceSpinCorrelations, "Yes", "Include the azimuthal spin correlations", 1); static Switch interfaceUseConstituentMasses ("UseConstituentMasses", "Whether or not to use constituent masses for the reconstruction of the particle after showering.", &ShowerHandler::useConstituentMasses_, true, false, false); static SwitchOption interfaceUseConstituentMassesYes (interfaceUseConstituentMasses, "Yes", "Use constituent masses.", true); static SwitchOption interfaceUseConstituentMassesNo (interfaceUseConstituentMasses, "No", "Don't use constituent masses.", false); + static Parameter interfaceTagIntermediates + ("TagIntermediates", + "Tag particles after shower with the given status code; if zero, no tagging will be performed.", + &ShowerHandler::tagIntermediates_, 0, 0, 0, + false, false, Interface::lowerlim); + } Energy ShowerHandler::hardScale() const { assert(false); return ZERO; } void ShowerHandler::cascade() { useMe(); // Initialise the weights in the event object // so that any variations are output regardless of // whether showering occurs for the given event initializeWeights(); // get the PDF's from ThePEG (if locally overridden use the local versions) tcPDFPtr first = PDFA_ ? tcPDFPtr(PDFA_) : firstPDF().pdf(); tcPDFPtr second = PDFB_ ? tcPDFPtr(PDFB_) : secondPDF().pdf(); resetPDFs(make_pair(first,second)); // set the PDFs for the remnant if( ! rempdfs_.first) rempdfs_.first = PDFARemnant_ ? PDFPtr(PDFARemnant_) : const_ptr_cast(first); if( ! rempdfs_.second) rempdfs_.second = PDFBRemnant_ ? PDFPtr(PDFBRemnant_) : const_ptr_cast(second); // get the incoming partons tPPair incomingPartons = eventHandler()->currentCollision()->primarySubProcess()->incoming(); // and the parton bins PBIPair incomingBins = make_pair(lastExtractor()->partonBinInstance(incomingPartons.first), lastExtractor()->partonBinInstance(incomingPartons.second)); // and the incoming hadrons tPPair incomingHadrons = eventHandler()->currentCollision()->incoming(); remnantDecayer()->setHadronContent(incomingHadrons); // check if incoming hadron == incoming parton // and get the incoming hadron if exists or parton otherwise incoming_ = make_pair(incomingBins.first ? incomingBins.first ->particle() : incomingPartons.first, incomingBins.second ? incomingBins.second->particle() : incomingPartons.second); // check the collision is of the beam particles // and if not boost collision to the right frame // i.e. the hadron-hadron CMF of the collision bool btotal(false); LorentzRotation rtotal; if(incoming_.first != incomingHadrons.first || incoming_.second != incomingHadrons.second ) { btotal = true; boostCollision(false); } // set the current ShowerHandler setCurrentHandler(); // first shower the hard process try { SubProPtr sub = eventHandler()->currentCollision()->primarySubProcess(); incomingPartons = cascade(sub,lastXCombPtr()); } catch(ShowerTriesVeto &veto){ throw Exception() << "Failed to generate the shower after " << veto.tries << " attempts in ShowerHandler::cascade()" << Exception::eventerror; } if(showerHardProcessVeto()) throw Veto(); + // tag particles as after shower + if ( tagIntermediates_ > 0 ) { + ParticleSet original = newStep()->particles(); + for ( auto p : original ) { + //if ( !p->data().coloured() ) continue; + newStep()->copyParticle(p); + p->status(tagIntermediates_); + } + } // if a non-hadron collision return (both incoming non-hadronic) if( ( !incomingBins.first|| !isResolvedHadron(incomingBins.first ->particle()))&& ( !incomingBins.second|| !isResolvedHadron(incomingBins.second->particle()))) { // boost back to lab if needed if(btotal) boostCollision(true); // perform the reweighting for the hard process shower combineWeights(); // unset the current ShowerHandler unSetCurrentHandler(); return; } // get the remnants for hadronic collision pair remnants(getRemnants(incomingBins)); // set the starting scale of the forced splitting to the PDF freezing scale remnantDecayer()->initialize(remnants, incoming_, *currentStep(), pdfFreezingScale()); // do the first forcedSplitting try { remnantDecayer()->doSplit(incomingPartons, make_pair(rempdfs_.first,rempdfs_.second), true); } catch (ExtraScatterVeto) { throw Exception() << "Remnant extraction failed in " << "ShowerHandler::cascade() from primary interaction. " << "Please check the PDFs you are using and set/unset them if necessary." << Exception::eventerror; } // perform the reweighting for the hard process shower combineWeights(); // if no MPI return if( !isMPIOn() ) { remnantDecayer()->finalize(); // boost back to lab if needed if(btotal) boostCollision(true); // unset the current ShowerHandler unSetCurrentHandler(); return; } // generate the multiple scatters use modified pdf's now: setMPIPDFs(); // additional "hard" processes unsigned int tries(0); // This is the loop over additional hard scatters (most of the time // only one, but who knows...) for(unsigned int i=1; i <= getMPIHandler()->additionalHardProcs(); i++){ //counter for regeneration unsigned int multSecond = 0; // generate the additional scatters while( multSecond < getMPIHandler()->multiplicity(i) ) { // generate the hard scatter tStdXCombPtr lastXC = getMPIHandler()->generate(i); SubProPtr sub = lastXC->construct(); // add to the Step newStep()->addSubProcess(sub); // increment the counters tries++; multSecond++; if(tries == maxtryDP_) throw Exception() << "Failed to establish the requested number " << "of additional hard processes. If this error " << "occurs often, your selection of additional " << "scatter is probably unphysical" << Exception::eventerror; // Generate the shower. If not possible veto the event try { incomingPartons = cascade(sub,lastXC); } catch(ShowerTriesVeto &veto){ throw Exception() << "Failed to generate the shower of " << "a secondary hard process after " << veto.tries << " attempts in Evolver::showerHardProcess()" << Exception::eventerror; } try { // do the forcedSplitting remnantDecayer()->doSplit(incomingPartons, make_pair(remmpipdfs_.first,remmpipdfs_.second), false); } catch(ExtraScatterVeto){ //remove all particles associated with the subprocess newStep()->removeParticle(incomingPartons.first); newStep()->removeParticle(incomingPartons.second); //remove the subprocess from the list newStep()->removeSubProcess(sub); //regenerate the scattering multSecond--; continue; } // connect with the remnants but don't set Remnant colour, // because that causes problems due to the multiple colour lines. if ( !remnants.first ->extract(incomingPartons.first , false) || !remnants.second->extract(incomingPartons.second, false) ) throw Exception() << "Remnant extraction failed in " << "ShowerHandler::cascade() for additional scatter" << Exception::runerror; } // perform the reweighting for the additional hard scatter shower combineWeights(); } // the underlying event processes unsigned int ptveto(1), veto(0); unsigned int max(getMPIHandler()->multiplicity()); for(unsigned int i=0; i maxtryMPI_) break; //generate PSpoint tStdXCombPtr lastXC = getMPIHandler()->generate(); SubProPtr sub = lastXC->construct(); //If Algorithm=1 additional scatters of the signal type // with pt > ptmin have to be vetoed //with probability 1/(m+1), where m is the number of occurances in this event if( getMPIHandler()->Algorithm() == 1 ){ //get the pT Energy pt = sub->outgoing().front()->momentum().perp(); if(pt > getMPIHandler()->PtForVeto() && UseRandom::rnd() < 1./(ptveto+1) ){ ptveto++; i--; continue; } } // add to the SubProcess to the step newStep()->addSubProcess(sub); // Run the Shower. If not possible veto the scattering try { incomingPartons = cascade(sub,lastXC); } // discard this extra scattering, but try the next one catch(ShowerTriesVeto) { newStep()->removeSubProcess(sub); //regenerate the scattering veto++; i--; continue; } try{ //do the forcedSplitting remnantDecayer()->doSplit(incomingPartons, make_pair(remmpipdfs_.first,remmpipdfs_.second), false); } catch (ExtraScatterVeto) { //remove all particles associated with the subprocess newStep()->removeParticle(incomingPartons.first); newStep()->removeParticle(incomingPartons.second); //remove the subprocess from the list newStep()->removeSubProcess(sub); //regenerate the scattering veto++; i--; continue; } //connect with the remnants but don't set Remnant colour, //because that causes problems due to the multiple colour lines. if ( !remnants.first ->extract(incomingPartons.first , false) || !remnants.second->extract(incomingPartons.second, false) ) throw Exception() << "Remnant extraction failed in " << "ShowerHandler::cascade() for MPI hard scattering" << Exception::runerror; //reset veto counter veto = 0; // perform the reweighting for the MPI process shower combineWeights(); } // finalize the remnants remnantDecayer()->finalize(getMPIHandler()->colourDisrupt(), getMPIHandler()->softMultiplicity()); // boost back to lab if needed if(btotal) boostCollision(true); // unset the current ShowerHandler unSetCurrentHandler(); getMPIHandler()->clean(); resetPDFs(make_pair(first,second)); } void ShowerHandler::initializeWeights() { if ( !showerVariations().empty() ) { tEventPtr event = eventHandler()->currentEvent(); for ( map::const_iterator var = showerVariations().begin(); var != showerVariations().end(); ++var ) { // Check that this is behaving as intended //map::iterator wi = event->optionalWeights().find(var->first); //assert(wi == event->optionalWeights().end() ); event->optionalWeights()[var->first] = 1.0; currentWeights_[var->first] = 1.0; } } reweight_ = 1.0; } void ShowerHandler::resetWeights() { for ( map::iterator w = currentWeights_.begin(); w != currentWeights_.end(); ++w ) { w->second = 1.0; } reweight_ = 1.0; } void ShowerHandler::combineWeights() { tEventPtr event = eventHandler()->currentEvent(); for ( map::const_iterator w = currentWeights_.begin(); w != currentWeights_.end(); ++w ) { map::iterator ew = event->optionalWeights().find(w->first); if ( ew != event->optionalWeights().end() ) ew->second *= w->second; else { assert(false && "Weight name unknown."); //event->optionalWeights()[w->first] = w->second; } } if ( reweight_ != 1.0 ) { Ptr::tptr eh = dynamic_ptr_cast::tptr>(eventHandler()); if ( !eh ) { throw Exception() << "ShowerHandler::combineWeights() : Cross section reweighting " << "through the shower is currently only available with standard " << "event generators" << Exception::runerror; } eh->reweight(reweight_); } } string ShowerHandler::doAddVariation(string in) { if ( in.empty() ) return "expecting a name and a variation specification"; string name = StringUtils::car(in); ShowerVariation var; string res = var.fromInFile(StringUtils::cdr(in)); if ( res.empty() ) { if ( !var.firstInteraction && !var.secondaryInteractions ) { // TODO what about decay showers? return "variation does not apply to any shower"; } if ( var.renormalizationScaleFactor == 1.0 && var.factorizationScaleFactor == 1.0 ) { return "variation does not vary anything"; } /* Repository::clog() << "adding a variation with tag '" << name << "' using\nxir = " << var.renormalizationScaleFactor << " xif = " << var.factorizationScaleFactor << "\napplying to:\n" << "first interaction = " << var.firstInteraction << " " << "secondary interactions = " << var.secondaryInteractions << "\n" << flush; */ showerVariations()[name] = var; } return res; } tPPair ShowerHandler::cascade(tSubProPtr, XCPtr) { assert(false); return tPPair(); } ShowerHandler::RemPair ShowerHandler::getRemnants(PBIPair incomingBins) { RemPair remnants; // first beam particle if(incomingBins.first&&!incomingBins.first->remnants().empty()) { remnants.first = dynamic_ptr_cast(incomingBins.first->remnants()[0] ); if(remnants.first) { ParticleVector children=remnants.first->children(); for(unsigned int ix=0;ixdataPtr()==remnants.first->dataPtr()) remnants.first = dynamic_ptr_cast(children[ix]); } //remove existing colour lines from the remnants if(remnants.first->colourLine()) remnants.first->colourLine()->removeColoured(remnants.first); if(remnants.first->antiColourLine()) remnants.first->antiColourLine()->removeAntiColoured(remnants.first); } } // seconnd beam particle if(incomingBins.second&&!incomingBins. second->remnants().empty()) { remnants.second = dynamic_ptr_cast(incomingBins.second->remnants()[0] ); if(remnants.second) { ParticleVector children=remnants.second->children(); for(unsigned int ix=0;ixdataPtr()==remnants.second->dataPtr()) remnants.second = dynamic_ptr_cast(children[ix]); } //remove existing colour lines from the remnants if(remnants.second->colourLine()) remnants.second->colourLine()->removeColoured(remnants.second); if(remnants.second->antiColourLine()) remnants.second->antiColourLine()->removeAntiColoured(remnants.second); } } assert(remnants.first || remnants.second); return remnants; } namespace { void addChildren(tPPtr in,set & particles) { particles.insert(in); for(unsigned int ix=0;ixchildren().size();++ix) addChildren(in->children()[ix],particles); } } void ShowerHandler::boostCollision(bool boost) { // calculate boost from lab to rest if(!boost) { Lorentz5Momentum ptotal=incoming_.first ->momentum()+incoming_.second->momentum(); boost_ = LorentzRotation(-ptotal.boostVector()); Axis axis((boost_*incoming_.first ->momentum()).vect().unit()); if(axis.perp2()>0.) { double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); boost_.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); } } // first call performs the boost and second inverse // get the particles to be boosted set particles; addChildren(incoming_.first,particles); addChildren(incoming_.second,particles); // apply the boost for(set::const_iterator cit=particles.begin(); cit!=particles.end();++cit) { (*cit)->transform(boost_); } if(!boost) boost_.invert(); } void ShowerHandler::setMPIPDFs() { if ( !mpipdfs_.first ) { // first have to check for MinBiasPDF tcMinBiasPDFPtr first = dynamic_ptr_cast(firstPDF().pdf()); if(first) mpipdfs_.first = new_ptr(MPIPDF(first->originalPDF())); else mpipdfs_.first = new_ptr(MPIPDF(firstPDF().pdf())); } if ( !mpipdfs_.second ) { tcMinBiasPDFPtr second = dynamic_ptr_cast(secondPDF().pdf()); if(second) mpipdfs_.second = new_ptr(MPIPDF(second->originalPDF())); else mpipdfs_.second = new_ptr(MPIPDF(secondPDF().pdf())); } if( !remmpipdfs_.first ) { tcMinBiasPDFPtr first = dynamic_ptr_cast(rempdfs_.first); if(first) remmpipdfs_.first = new_ptr(MPIPDF(first->originalPDF())); else remmpipdfs_.first = new_ptr(MPIPDF(rempdfs_.first)); } if( !remmpipdfs_.second ) { tcMinBiasPDFPtr second = dynamic_ptr_cast(rempdfs_.second); if(second) remmpipdfs_.second = new_ptr(MPIPDF(second->originalPDF())); else remmpipdfs_.second = new_ptr(MPIPDF(rempdfs_.second)); } // reset the PDFs stored in the base class resetPDFs(mpipdfs_); } bool ShowerHandler::isResolvedHadron(tPPtr particle) { if(!HadronMatcher::Check(particle->data())) return false; for(unsigned int ix=0;ixchildren().size();++ix) { if(particle->children()[ix]->id()==ParticleID::Remnant) return true; } return false; } namespace { bool decayProduct(tSubProPtr subProcess, tPPtr particle) { // must be time-like and not incoming if(particle->momentum().m2()<=ZERO|| particle == subProcess->incoming().first|| particle == subProcess->incoming().second) return false; // if only 1 outgoing and this is it if(subProcess->outgoing().size()==1 && subProcess->outgoing()[0]==particle) return true; // must not be the s-channel intermediate otherwise if(find(subProcess->incoming().first->children().begin(), subProcess->incoming().first->children().end(),particle)!= subProcess->incoming().first->children().end()&& find(subProcess->incoming().second->children().begin(), subProcess->incoming().second->children().end(),particle)!= subProcess->incoming().second->children().end()&& subProcess->incoming().first ->children().size()==1&& subProcess->incoming().second->children().size()==1) return false; // if non-coloured this is enough if(!particle->dataPtr()->coloured()) return true; // if coloured must be unstable if(particle->dataPtr()->stable()) return false; // must not have same particle type as a child int id = particle->id(); for(unsigned int ix=0;ixchildren().size();++ix) if(particle->children()[ix]->id()==id) return false; // otherwise its a decaying particle return true; } PPtr findParent(PPtr original, bool & isHard, set outgoingset, tSubProPtr subProcess) { PPtr parent=original; isHard |=(outgoingset.find(original) != outgoingset.end()); if(!original->parents().empty()) { PPtr orig=original->parents()[0]; if(decayProduct(subProcess,orig)) parent=findParent(orig,isHard,outgoingset,subProcess); } return parent; } } void ShowerHandler::findDecayProducts(PPtr in,PerturbativeProcessPtr hard, DecayProcessMap & decay) const { ParticleVector children=in->children(); for(ParticleVector::const_iterator it=children.begin(); it!=children.end();++it) { // if decayed or should be decayed in shower make the PerturbaitveProcess bool radiates = false; if(!(**it).children().empty()) { // remove d,u,s,c,b quarks and leptons other than on-shell taus if( StandardQCDPartonMatcher::Check((**it).id()) || ( LeptonMatcher::Check((**it).id()) && !(abs((**it).id())==ParticleID::tauminus && abs((**it).mass()-(**it).dataPtr()->mass())id()==(**it).id()) { foundParticle = true; } else if((**it).children()[iy]->id()==ParticleID::g || (**it).children()[iy]->id()==ParticleID::gamma) { foundGauge = true; } } radiates = foundParticle && foundGauge; } } if(radiates) { findDecayProducts(*it,hard,decay); } else if(!(**it).children().empty()|| (decaysInShower((**it).id())&&!(**it).dataPtr()->stable())) { createDecayProcess(*it,hard,decay); } else { hard->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr())); } } } void ShowerHandler::splitHardProcess(tPVector tagged, PerturbativeProcessPtr & hard, DecayProcessMap & decay) const { // temporary storage of the particles set hardParticles; // tagged particles in a set set outgoingset(tagged.begin(),tagged.end()); bool isHard=false; // loop over the tagged particles for (tParticleVector::const_iterator taggedP = tagged.begin(); taggedP != tagged.end(); ++taggedP) { // skip remnants if (eventHandler()->currentCollision()&& eventHandler()->currentCollision()->isRemnant(*taggedP)) continue; // find the parent and whether its a decaying particle bool isDecayProd=false; // check if hard isHard |=(outgoingset.find(*taggedP) != outgoingset.end()); if(splitHardProcess_) { tPPtr parent = *taggedP; // check if from s channel decaying colourless particle while(parent&&!parent->parents().empty()&&!isDecayProd) { parent = parent->parents()[0]; if(parent == subProcess_->incoming().first || parent == subProcess_->incoming().second ) break; isDecayProd = decayProduct(subProcess_,parent); } if (isDecayProd) hardParticles.insert(findParent(parent,isHard,outgoingset,subProcess_)); } if (!isDecayProd) hardParticles.insert(*taggedP); } // there must be something to shower if(hardParticles.empty()) throw Exception() << "No particles to shower in " << "ShowerHandler::splitHardProcess()" << Exception::eventerror; // must be a hard process if(!isHard) throw Exception() << "Starting on decay not yet implemented in " << "ShowerHandler::splitHardProcess()" << Exception::runerror; // create the hard process hard = new_ptr(PerturbativeProcess()); // incoming particles hard->incoming().push_back(make_pair(subProcess_->incoming().first ,PerturbativeProcessPtr())); hard->incoming().push_back(make_pair(subProcess_->incoming().second,PerturbativeProcessPtr())); // outgoing particles for(set::const_iterator it=hardParticles.begin();it!=hardParticles.end();++it) { // if decayed or should be decayed in shower make the tree PPtr orig = *it; bool radiates = false; if(!orig->children().empty()) { // remove d,u,s,c,b quarks and leptons other than on-shell taus if( StandardQCDPartonMatcher::Check(orig->id()) || ( LeptonMatcher::Check(orig->id()) && !(abs(orig->id())==ParticleID::tauminus && abs(orig->mass()-orig->dataPtr()->mass())children().size();++iy) { if(orig->children()[iy]->id()==orig->id()) { foundParticle = true; } else if(orig->children()[iy]->id()==ParticleID::g || orig->children()[iy]->id()==ParticleID::gamma) { foundGauge = true; } } radiates = foundParticle && foundGauge; } } if(radiates) { findDecayProducts(orig,hard,decay); } else if(!(**it).children().empty()|| (decaysInShower((**it).id())&&!(**it).dataPtr()->stable())) { createDecayProcess(*it,hard,decay); } else { hard->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr())); } } } void ShowerHandler::createDecayProcess(PPtr in,PerturbativeProcessPtr hard, DecayProcessMap & decay) const { // there must be an incoming particle assert(in); // create the new process and connect with the parent PerturbativeProcessPtr newDecay=new_ptr(PerturbativeProcess()); newDecay->incoming().push_back(make_pair(in,hard)); Energy width=in->dataPtr()->generateWidth(in->mass()); decay.insert(make_pair(width,newDecay)); hard->outgoing().push_back(make_pair(in,newDecay)); // we need to deal with the decay products if decayed ParticleVector children = in->children(); if(!children.empty()) { for(ParticleVector::const_iterator it = children.begin(); it!= children.end(); ++it) { // if decayed or should be decayed in shower make the tree in->abandonChild(*it); bool radiates = false; if(!(**it).children().empty()) { if(StandardQCDPartonMatcher::Check((**it).id())|| (LeptonMatcher::Check((**it).id())&& !(abs((**it).id())==ParticleID::tauminus && abs((**it).mass()-(**it).dataPtr()->mass())id()==(**it).id()) { foundParticle = true; } else if((**it).children()[iy]->id()==ParticleID::g || (**it).children()[iy]->id()==ParticleID::gamma) { foundGauge = true; } } radiates = foundParticle && foundGauge; } // finally assume all non-decaying particles are in this class // pr 27/11/15 not sure about this bit // if(!radiates) { // radiates = !decaysInShower((**it).id()); // } } if(radiates) { findDecayProducts(*it,newDecay,decay); } else if(!(**it).children().empty()|| (decaysInShower((**it).id())&&!(**it).dataPtr()->stable())) { createDecayProcess(*it,newDecay,decay); } else { newDecay->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr())); } } } } tDMPtr ShowerHandler::decay(PerturbativeProcessPtr process, DecayProcessMap & decayMap, bool radPhotons ) const { PPtr parent = process->incoming()[0].first; assert(parent); if(parent->spinInfo()) parent->spinInfo()->decay(true); unsigned int ntry = 0; ParticleVector children; tDMPtr dm = DMPtr(); while (true) { // exit if fails if (++ntry>=maxtryDecay_) throw Exception() << "Failed to perform decay in ShowerHandler::decay()" << " after " << maxtryDecay_ << " attempts for " << parent->PDGName() << Exception::eventerror; // select decay mode dm = parent->data().selectMode(*parent); if(!dm) throw Exception() << "Failed to select decay mode in ShowerHandler::decay()" << "for " << parent->PDGName() << Exception::eventerror; if(!dm->decayer()) throw Exception() << "No Decayer for selected decay mode " << " in ShowerHandler::decay()" << Exception::runerror; // start of try block try { children = dm->decayer()->decay(*dm, *parent); // if no children have another go if(children.empty()) continue; if(radPhotons){ // generate radiation in the decay tDecayIntegratorPtr hwdec=dynamic_ptr_cast(dm->decayer()); if (hwdec && hwdec->canGeneratePhotons()) children = hwdec->generatePhotons(*parent,children); } // set up parent parent->decayMode(dm); // add children for (unsigned int i = 0, N = children.size(); i < N; ++i ) { children[i]->setLabVertex(parent->labDecayVertex()); //parent->addChild(children[i]); } // if succeeded break out of loop break; } catch(Veto) { } } assert(!children.empty()); for(ParticleVector::const_iterator it = children.begin(); it!= children.end(); ++it) { if(!(**it).children().empty()|| (decaysInShower((**it).id())&&!(**it).dataPtr()->stable())) { createDecayProcess(*it,process,decayMap); } else { process->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr())); } } return dm; } // Note: The tag must be constructed from an ordered particle container. tDMPtr ShowerHandler::findDecayMode(const string & tag) const { static map cache; map::const_iterator pos = cache.find(tag); if ( pos != cache.end() ) return pos->second; tDMPtr dm = CurrentGenerator::current().findDecayMode(tag); cache[tag] = dm; return dm; } /** * Operator for the particle ordering * @param p1 The first ParticleData object * @param p2 The second ParticleData object */ bool ShowerHandler::ParticleOrdering::operator() (tcPDPtr p1, tcPDPtr p2) const { return abs(p1->id()) > abs(p2->id()) || ( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) || ( p1->id() == p2->id() && p1->fullName() > p2->fullName() ); } diff --git a/Shower/ShowerHandler.h b/Shower/ShowerHandler.h --- a/Shower/ShowerHandler.h +++ b/Shower/ShowerHandler.h @@ -1,898 +1,908 @@ // -*- C++ -*- // // ShowerHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_ShowerHandler_H #define HERWIG_ShowerHandler_H // // This is the declaration of the ShowerHandler class. // #include "ThePEG/Handlers/EventHandler.h" #include "ThePEG/Handlers/CascadeHandler.h" #include "ShowerVariation.h" #include "Herwig/PDF/HwRemDecayer.fh" #include "ThePEG/EventRecord/RemnantParticle.fh" #include "UEBase.h" #include "PerturbativeProcess.h" #include "Herwig/MatrixElement/Matchbox/Matching/HardScaleProfile.h" #include "ShowerHandler.fh" namespace Herwig { using namespace ThePEG; /** \ingroup Shower * * This class is the main driver of the shower: it is responsible for * the proper handling of all other specific collaborating classes * and for the storing of the produced particles in the event record. * * @see \ref ShowerHandlerInterfaces "The interfaces" * * @see ThePEG::CascadeHandler * @see MPIHandler * @see HwRemDecayer */ class ShowerHandler: public CascadeHandler { public: /** * Typedef for a pair of ThePEG::RemnantParticle pointers. */ typedef pair RemPair; public: /** * Default constructor */ ShowerHandler(); /** * Destructor */ virtual ~ShowerHandler(); public: /** * The main method which manages the multiple interactions and starts * the shower by calling cascade(sub, lastXC). */ virtual void cascade(); /** * pointer to "this", the current ShowerHandler. */ static const tShowerHandlerPtr currentHandler() { assert(currentHandler_); return currentHandler_; } /** * pointer to "this", the current ShowerHandler. */ static bool currentHandlerIsSet() { return currentHandler_; } public: /** * Hook to allow vetoing of event after showering hard sub-process * as in e.g. MLM merging. */ virtual bool showerHardProcessVeto() const { return false; } /** * Return true, if this cascade handler will perform reshuffling from hard * process masses. */ virtual bool isReshuffling() const { return true; } /** * Return true, if this cascade handler will put the final state * particles to their constituent mass. If false the nominal mass is used. */ virtual bool retConstituentMasses() const { return useConstituentMasses_; } /** * Return true, if the shower handler can generate a truncated * shower for POWHEG style events generated using Matchbox */ virtual bool canHandleMatchboxTrunc() const { return false; } /** * Get the PDF freezing scale */ Energy pdfFreezingScale() const { return pdfFreezingScale_; } /** * Get the local PDFs. */ PDFPtr getPDFA() const {return PDFA_;} /** * Get the local PDFs. */ PDFPtr getPDFB() const {return PDFB_;} /** * Return true if currently the primary subprocess is showered. */ bool firstInteraction() const { if (!eventHandler()->currentCollision())return true; return ( subProcess_ == eventHandler()->currentCollision()->primarySubProcess() ); } /** * Return the remnant decayer. */ tHwRemDecPtr remnantDecayer() const { return remDec_; } /** * Split the hard process into production and decays * @param tagged The tagged particles from the StepHandler * @param hard The hard perturbative process * @param decay The decay particles */ void splitHardProcess(tPVector tagged, PerturbativeProcessPtr & hard, DecayProcessMap & decay) const; /** * Information if the Showerhandler splits the hard process. */ bool doesSplitHardProcess()const {return splitHardProcess_;} /** * Decay a particle. * radPhotons switches the generation of photon * radiation on/off. * Required for Dipole Shower but not QTilde Shower. */ tDMPtr decay(PerturbativeProcessPtr, DecayProcessMap & decay, bool radPhotons = false) const; /** * Cached lookup of decay modes. * Generator::findDecayMode() is not efficient. */ tDMPtr findDecayMode(const string & tag) const; /** * A struct to order the particles in the same way as in the DecayMode's */ struct ParticleOrdering { bool operator() (tcPDPtr p1, tcPDPtr p2) const; }; /** * A container for ordered particles required * for constructing tags for decay mode lookup. */ typedef multiset OrderedParticles; public: /** * @name Switches for initial- and final-state radiation */ //@{ /** * Switch for any radiation */ bool doRadiation() const {return doFSR_ || doISR_;} /** * Switch on or off final state radiation. */ bool doFSR() const { return doFSR_;} /** * Switch on or off initial state radiation. */ bool doISR() const { return doISR_;} //@} public: /** * @name Switches for scales */ //@{ /** * Return true if maximum pt should be deduced from the factorization scale */ bool hardScaleIsMuF() const { return maxPtIsMuF_; } /** * The factorization scale factor. */ double factorizationScaleFactor() const { return factorizationScaleFactor_; } /** * The renormalization scale factor. */ double renFac() const { return renormalizationScaleFactor_; } /** * The factorization scale factor. */ double facFac() const { return factorizationScaleFactor_; } /** * The renormalization scale factor. */ double renormalizationScaleFactor() const { return renormalizationScaleFactor_; } /** * The scale factor for the hard scale */ double hardScaleFactor() const { return hardScaleFactor_; } /** * Return true, if the phase space restrictions of the dipole shower should * be applied. */ bool restrictPhasespace() const { return restrictPhasespace_; } /** * Return profile scales */ Ptr::tptr profileScales() const { return hardScaleProfile_; } /** * Return the relevant hard scale to be used in the profile scales */ virtual Energy hardScale() const; /** * Return information about shower phase space choices */ virtual int showerPhaseSpaceOption() const { assert(false && "not implemented in general"); return -1; } //@} public: /** * Access the shower variations */ map& showerVariations() { return showerVariations_; } /** * Return the shower variations */ const map& showerVariations() const { return showerVariations_; } /** * Access the current Weights */ map& currentWeights() { return currentWeights_; } /** * Return the current Weights */ const map& currentWeights() const { return currentWeights_; } /** * Change the current reweighting factor */ void reweight(double w) { reweight_ = w; } /** * Return the current reweighting factor */ double reweight() const { return reweight_; } public : /** * Access to switches for spin correlations */ //@{ /** * Spin Correlations */ unsigned int spinCorrelations() const { return spinOpt_; } /** * Any correlations */ virtual bool correlations() const { return spinOpt_!=0; } //@} public: /** * struct that is used to catch exceptions which are thrown * due to energy conservation issues of additional scatters */ struct ExtraScatterVeto {}; /** * struct that is used to catch exceptions which are thrown * due to fact that the Shower has been invoked more than * a defined threshold on a certain configuration */ struct ShowerTriesVeto { /** variable to store the number of attempts */ const int tries; /** constructor */ ShowerTriesVeto(int t) : tries(t) {} }; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Functions to perform the cascade */ //@{ /** * The main method which manages the showering of a subprocess. */ virtual tPPair cascade(tSubProPtr sub, XCPtr xcomb); /** * Set up for the cascade */ void prepareCascade(tSubProPtr sub) { current_ = currentStep(); subProcess_ = sub; } /** * Boost all the particles in the collision so that the collision always occurs * in the rest frame with the incoming particles along the z axis */ void boostCollision(bool boost); //@} protected: /** * Set/unset the current shower handler */ //@{ /** * Set the current handler */ void setCurrentHandler() { currentHandler_ = tShowerHandlerPtr(this); } /** * Unset the current handler */ void unSetCurrentHandler() { currentHandler_ = tShowerHandlerPtr(); } //@} protected: /** * @name Members relating to the underlying event and MPI */ //@{ /** * Return true if multiple parton interactions are switched on * and can be used for this beam setup. */ bool isMPIOn() const { return MPIHandler_ && MPIHandler_->beamOK(); } /** * Access function for the MPIHandler, it should only be called after * checking with isMPIOn. */ tUEBasePtr getMPIHandler() const { assert(MPIHandler_); return MPIHandler_; } /** * Is a beam particle where hadronic structure is resolved */ bool isResolvedHadron(tPPtr); /** * Get the remnants from the ThePEG::PartonBinInstance es and * do some checks. */ RemPair getRemnants(PBIPair incbins); /** * Reset the PDF's after the hard collision has been showered */ void setMPIPDFs(); //@} public: /** * Check if a particle decays in the shower * @param id The PDG code for the particle */ bool decaysInShower(long id) const { return ( particlesDecayInShower_.find( abs(id) ) != particlesDecayInShower_.end() ); } protected: /** * Members to handle splitting up of hard process and decays */ //@{ /** * Find decay products from the hard process and create decay processes * @param parent The parent particle * @param hard The hard process * @param decay The decay processes */ void findDecayProducts(PPtr parent, PerturbativeProcessPtr hard, DecayProcessMap & decay) const; /** * Find decay products from the hard process and create decay processes * @param parent The parent particle * @param hard The parent hard process * @param decay The decay processes */ void createDecayProcess(PPtr parent,PerturbativeProcessPtr hard, DecayProcessMap & decay) const; //@} /** * @name Functions to return information relevant to the process being showered */ //@{ /** * Return the currently used SubProcess. */ tSubProPtr currentSubProcess() const { assert(subProcess_); return subProcess_; } /** * Access to the incoming beam particles */ tPPair incomingBeams() const { return incoming_; } //@} protected: /** * Weight handling for shower variations */ //@ /** * Combine the variation weights which have been encountered */ void combineWeights(); /** * Initialise the weights in currentEvent() */ void initializeWeights(); /** * Reset the current weights */ void resetWeights(); //@} protected: /** * Return the maximum number of attempts for showering * a given subprocess. */ unsigned int maxtry() const { return maxtry_; } protected: /** * Parameters for the space-time model */ //@{ /** * Whether or not to include spa-cetime distances in the shower */ bool includeSpaceTime() const {return includeSpaceTime_;} /** * The minimum virtuality for the space-time model */ Energy2 vMin() const {return vMin_;} //@} protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @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(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); /** * Finalize this object. Called in the run phase just after a * run has ended. Used eg. to write out statistics. */ virtual void dofinish(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ ShowerHandler & operator=(const ShowerHandler &) = delete; private: /** * pointer to "this", the current ShowerHandler. */ static tShowerHandlerPtr currentHandler_; /** * a MPIHandler to administer the creation of several (semihard) * partonic interactions. */ UEBasePtr MPIHandler_; /** * Pointer to the HwRemDecayer */ HwRemDecPtr remDec_; private: /** * Maximum tries for various stages of the showering process */ //@{ /** * Maximum number of attempts for the * main showering loop */ unsigned int maxtry_; /** * Maximum number of attempts for the regeneration of an additional * scattering, before the number of scatters is reduced. */ unsigned int maxtryMPI_; /** * Maximum number of attempts for the regeneration of an additional * hard scattering, before this event is vetoed. */ unsigned int maxtryDP_; /** * Maximum number of attempts to generate a decay */ unsigned int maxtryDecay_; //@} private: /** * Factors for the various scales */ //@{ /** * The factorization scale factor. */ double factorizationScaleFactor_; /** * The renormalization scale factor. */ double renormalizationScaleFactor_; /** * The scale factor for the hard scale */ double hardScaleFactor_; /** * True, if the phase space restrictions of the dipole shower should * be applied. */ bool restrictPhasespace_; /** * True if maximum pt should be deduced from the factorization scale */ bool maxPtIsMuF_; /** * The profile scales */ Ptr::ptr hardScaleProfile_; //@} /** * Option to include spin correlations */ unsigned int spinOpt_; private: /** * Storage of information about the current event */ //@{ /** * The incoming beam particles for the current collision */ tPPair incoming_; /** * Boost to get back to the lab */ LorentzRotation boost_; /** * Const pointer to the currently handeled ThePEG::SubProcess */ tSubProPtr subProcess_; /** * Const pointer to the current step */ tcStepPtr current_; //@} private: /** * PDFs to be used for the various stages and related parameters */ //@{ /** * The PDF freezing scale */ Energy pdfFreezingScale_; /** * PDFs to be used for the various stages and related parameters */ //@{ /** * The PDF for beam particle A. Overrides the particle's own PDF setting. */ PDFPtr PDFA_; /** * The PDF for beam particle B. Overrides the particle's own PDF setting. */ PDFPtr PDFB_; /** * The PDF for beam particle A for remnant splitting. Overrides the particle's own PDF setting. */ PDFPtr PDFARemnant_; /** * The PDF for beam particle B for remnant splitting. Overrides the particle's own PDF setting. */ PDFPtr PDFBRemnant_; /** * The MPI PDF's to be used for secondary scatters. */ pair mpipdfs_; /** * The MPI PDF's to be used for secondary scatters. */ pair rempdfs_; /** * The MPI PDF's to be used for secondary scatters. */ pair remmpipdfs_; //@} private: /** * @name Parameters for initial- and final-state radiation */ //@{ /** * Switch on or off final state radiation. */ bool doFSR_; /** * Switch on or off initial state radiation. */ bool doISR_; //@} private: /** * @name Parameters for particle decays */ //@{ /** * Whether or not to split into hard and decay trees */ bool splitHardProcess_; /** * PDG codes of the particles which decay during showering * this is fast storage for use during running */ set particlesDecayInShower_; /** * PDG codes of the particles which decay during showering * this is a vector that is interfaced so they can be changed */ vector inputparticlesDecayInShower_; //@} private: /** * Parameters for the space-time model */ //@{ /** - * Whether or not to include spa-cetime distances in the shower + * Whether or not to include space-time distances in the shower */ bool includeSpaceTime_; /** * The minimum virtuality for the space-time model */ Energy2 vMin_; //@} private: /** * Parameters for the constituent mass treatment. */ - //@{ - // True if shower should return constituent masses. - bool useConstituentMasses_=true; + //@{ + + /** + * True if shower should return constituent masses. + */ + bool useConstituentMasses_ = true; + + /** + * Status code to tag intermediate state as after the showering; if zero no tagging will be performed + */ + int tagIntermediates_ = 0; + //@} + private: /** * Parameters relevant for reweight and variations */ //@{ /** * The shower variations */ map showerVariations_; /** * Command to add a shower variation */ string doAddVariation(string); /** * A reweighting factor applied by the showering */ double reweight_; /** * The shower variation weights */ map currentWeights_; //@} }; } #endif /* HERWIG_ShowerHandler_H */ diff --git a/Utilities/Makefile.am b/Utilities/Makefile.am --- a/Utilities/Makefile.am +++ b/Utilities/Makefile.am @@ -1,46 +1,47 @@ SUBDIRS = XML Statistics noinst_LTLIBRARIES = libHwUtils.la libHwUtils_la_SOURCES = \ EnumParticles.h \ Interpolator.tcc Interpolator.h \ Kinematics.cc Kinematics.h \ Progress.h Progress.cc \ Maths.h Maths.cc \ StandardSelectors.cc StandardSelectors.h\ Histogram.cc Histogram.fh Histogram.h \ GaussianIntegrator.cc GaussianIntegrator.h \ GaussianIntegrator.tcc \ Statistic.h HerwigStrategy.cc HerwigStrategy.h \ GSLIntegrator.h GSLIntegrator.tcc \ GSLBisection.h GSLBisection.tcc GSLHelper.h \ +Reshuffler.h Reshuffler.cc \ expm-1.h \ HiggsLoopFunctions.h AlphaS.h nodist_libHwUtils_la_SOURCES = hgstamp.inc BUILT_SOURCES = hgstamp.inc CLEANFILES = hgstamp.inc HGVERSION := $(shell hg -R $(top_srcdir) parents --template '"Herwig {node|short} ({branch})"' 2> /dev/null || echo \"$(PACKAGE_STRING)\" || true ) .PHONY: update_hgstamp hgstamp.inc: update_hgstamp @[ -f $@ ] || touch $@ @echo '$(HGVERSION)' | cmp -s $@ - || echo '$(HGVERSION)' > $@ libHwUtils_la_LIBADD = \ XML/libHwXML.la \ Statistics/libHwStatistics.la check_PROGRAMS = utilities_test utilities_test_SOURCES = \ tests/utilitiesTestsMain.cc \ tests/utilitiesTestsGlobalFixture.h \ tests/utilitiesTestsKinematics.h \ tests/utilitiesTestMaths.h \ tests/utilitiesTestsStatistic.h utilities_test_LDADD = $(BOOST_UNIT_TEST_FRAMEWORK_LIBS) $(THEPEGLIB) -ldl libHwUtils.la utilities_test_LDFLAGS = $(AM_LDFLAGS) -export-dynamic $(BOOST_UNIT_TEST_FRAMEWORK_LDFLAGS) $(THEPEGLDFLAGS) utilities_test_CPPFLAGS = $(AM_CPPFLAGS) $(BOOST_CPPFLAGS) TESTS = utilities_test diff --git a/Utilities/Reshuffler.cc b/Utilities/Reshuffler.cc new file mode 100644 --- /dev/null +++ b/Utilities/Reshuffler.cc @@ -0,0 +1,85 @@ +// -*- C++ -*- +// +// Reshuffler.h is a part of Herwig - A multi-purpose Monte Carlo event generator +// Copyright (C) 2002-2019 The Herwig Collaboration +// +// Herwig is licenced under version 3 of the GPL, see COPYING for details. +// Please respect the MCnet academic guidelines, see GUIDELINES for details. +// +// +// This is the implementation of the non-inlined, non-templated member +// functions of the Reshuffler class. +// + +#include +#include "Reshuffler.h" +#include "Herwig/Utilities/GSLBisection.h" +#include "Herwig/Shower/ShowerHandler.h" + +using namespace Herwig; + +void Reshuffler::reshuffle(const PVector& particles, + const vector& masses) const { + + Energy zero (0.*GeV); + Lorentz5Momentum Q (zero,zero,zero,zero); + + for (PVector::const_iterator p = particles.begin(); + p != particles.end(); ++p) { + Q += (**p).momentum(); + } + + Boost beta = Q.findBoostToCM(); + + vector mbackup; + + bool need_boost = (beta.mag2() > Constants::epsilon); + + if (need_boost) { + + for (PVector::const_iterator p = particles.begin(); + p != particles.end(); ++p) { + Lorentz5Momentum mom = (**p).momentum(); + mbackup.push_back(mom); + (**p).set5Momentum(mom.boost(beta)); + } + + } + + double xi; + + ReshuffleEquation::const_iterator> + solve (Q.m(),particles.begin(),particles.end(), + masses.begin(),masses.end()); + + GSLBisection solver(1e-10,1e-8,10000); + + try { + xi = solver.value(solve,0.0,1.1); + } catch (GSLBisection::GSLerror) { + throw Exception("Failed to reshuffle.",Exception::eventerror); + } catch (GSLBisection::IntervalError) { + throw Exception("Failed to reshuffle.",Exception::eventerror); + } + + PVector::const_iterator p = particles.begin(); + vector::const_iterator m = masses.begin(); + for (; p != particles.end(); ++p, ++m) { + + Lorentz5Momentum rm = Lorentz5Momentum (xi*(**p).momentum().x(), + xi*(**p).momentum().y(), + xi*(**p).momentum().z(), + sqrt(sqr(*m) + + xi*xi*(sqr((**p).momentum().t())-sqr((**p).dataPtr()->mass())))); + + rm.rescaleMass(); + + if (need_boost) { + rm.boost(-beta); + } + + (**p).set5Momentum(rm); + + } + +} diff --git a/Utilities/Reshuffler.h b/Utilities/Reshuffler.h new file mode 100644 --- /dev/null +++ b/Utilities/Reshuffler.h @@ -0,0 +1,86 @@ +// -*- C++ -*- +// +// Reshuffler.h is a part of Herwig - A multi-purpose Monte Carlo event generator +// Copyright (C) 2002-2019 The Herwig Collaboration +// +// Herwig is licenced under version 3 of the GPL, see COPYING for details. +// Please respect the MCnet academic guidelines, see GUIDELINES for details. +// +#ifndef HERWIG_Reshuffler_H +#define HERWIG_Reshuffler_H +// +// This is the declaration of the Reshuffler class. +// + +#include "ThePEG/Config/ThePEG.h" + +namespace Herwig { + +using namespace ThePEG; + +/** + * \author Simon Platzer, Stephen Webster + * + * \brief The Reshuffler class implements reshuffling + * of partons on their nominal mass shell to their constituent + * mass shells. + * + */ +class Reshuffler { + +protected: + + /** + * The function object defining the equation + * to be solved. + */ + template + struct ReshuffleEquation { + + ReshuffleEquation (Energy q, + PIterator pp_begin, + PIterator pp_end, + MIterator mm_begin, + MIterator mm_end) + : w(q), p_begin(pp_begin), p_end(pp_end), + m_begin(mm_begin), m_end(mm_end) {} + + typedef double ArgType; + typedef double ValType; + + static double aUnit() { return 1.; } + static double vUnit() { return 1.; } + + double operator() (double xi) const { + double r = - w/GeV; + PIterator p = p_begin; + MIterator m = m_begin; + for (; p != p_end; ++p, ++m) { + r += sqrt(sqr(*m) + + xi*xi*(sqr((**p).momentum().t())-sqr((**p).dataPtr()->mass()))) / GeV; + } + return r; + } + + Energy w; + + PIterator p_begin; + PIterator p_end; + + MIterator m_begin; + MIterator m_end; + + }; + + /** + * Reshuffle to consitutent masses + */ + void reshuffle(const PVector& particles, + const vector& masses) const; + +}; + +} + +#endif /* HERWIG_Reshuffler_H */ + diff --git a/patch_added_newColourReconnection.diff b/patch_added_newColourReconnection.diff deleted file mode 100644 --- a/patch_added_newColourReconnection.diff +++ /dev/null @@ -1,1680 +0,0 @@ -diff -r 1781af1e3113 -r 5d17715d8d26 Hadronization/ColourReconnector.cc ---- a/Hadronization/ColourReconnector.cc Mon Jul 18 09:45:26 2022 +0100 -+++ b/Hadronization/ColourReconnector.cc Wed Jul 20 18:37:09 2022 +0200 -@@ -51,9 +51,18 @@ - - // do the colour reconnection - switch (_algorithm) { -- case 0: _doRecoPlain(clusters); break; -- case 1: _doRecoStatistical(clusters); break; -- case 2: _doRecoBaryonic(clusters); break; -+ case 0: -+ _doRecoPlain(clusters); -+ break; -+ case 1: -+ _doRecoStatistical(clusters); -+ break; -+ case 2: -+ _doRecoBaryonic(clusters); -+ break; -+ case 3: -+ _doRecoBaryonicMesonic(clusters); -+ break; - } - } - -@@ -68,7 +77,40 @@ - return sum; - } - -+double ColourReconnector::_displacement(tcPPtr p, tcPPtr q) const { -+ double deltaRap = (p->rapidity() - q->rapidity()); -+ double deltaPhi = (p->momentum().phi() - q->momentum().phi()); - -+ return sqrt(deltaRap * deltaRap + deltaPhi * deltaPhi); -+} -+ -+ -+double ColourReconnector::_displacementBaryonic(tcPPtr q1, tcPPtr q2, tcPPtr q3) const { -+ if (_junctionMBCR) { -+ /** -+ * Junction-like option i.e. displacement -+ * from "junction centre" (mean rapidity/phi) -+ */ -+ double rap1=q1->rapidity(); -+ double rap2=q2->rapidity(); -+ double rap3=q3->rapidity(); -+ -+ double phi1=q1->momentum().phi(); -+ double phi2=q2->momentum().phi(); -+ double phi3=q3->momentum().phi(); -+ double meanRap=(rap1 + rap2 + rap3)/3.0; -+ double meanPhi=(phi1 + phi2 + phi3)/3.0; -+ double delR; -+ -+ delR = sqrt( (rap1-meanRap)*(rap1-meanRap) + (phi1-meanPhi)*(phi1-meanPhi) ); -+ delR += sqrt( (rap2-meanRap)*(rap2-meanRap) + (phi2-meanPhi)*(phi2-meanPhi) ); -+ delR += sqrt( (rap3-meanRap)*(rap3-meanRap) + (phi3-meanPhi)*(phi3-meanPhi) ); -+ return delR; -+ } else { -+ /* just summing up all possible 2 quark displacements */ -+ return _displacement(q1, q2) + _displacement(q1, q3) + _displacement(q2, q3); -+ } -+} - bool ColourReconnector::_containsColour8(const ClusterVector & cv, - const vector & P) const { - assert (P.size() == cv.size()); -@@ -202,13 +244,13 @@ - } - - namespace { -- inline bool hasDiquark(CluVecIt cit) { -- for(int i = 0; i<(*cit)->numComponents(); i++) { -+inline bool hasDiquark(CluVecIt cit) { -+ for (int i = 0; i<(*cit)->numComponents(); i++) { - if (DiquarkMatcher::Check(*((*cit)->particle(i)->dataPtr()))) - return true; - } - return false; -- } -+} - } - - -@@ -221,7 +263,7 @@ - - // try to avoid systematic errors by randomising the reconnection order - long (*p_irnd)(long) = UseRandom::irnd; -- random_shuffle( newcv.begin(), newcv.end(), p_irnd ); -+ random_shuffle(newcv.begin(), newcv.end(), p_irnd); - - // iterate over all clusters - for (CluVecIt cit = newcv.begin(); cit != newcv.end(); ++cit) { -@@ -245,16 +287,16 @@ - baryonic1, baryonic2); - - // skip this cluster if no possible reconnection partner can be found -- if ( !isBaryonicCandidate && candidate==cit ) -+ if (!isBaryonicCandidate && candidate==cit) - continue; - -- if ( isBaryonicCandidate -- && UseRandom::rnd() < _precoBaryonic ) { -+ if (isBaryonicCandidate -+ && UseRandom::rnd() < _precoBaryonic) { - deleted.push_back(*baryonic2); - - // Function that does the reconnection from 3 -> 2 clusters - ClusterPtr b1, b2; -- _makeBaryonicClusters(*cit,*baryonic1,*baryonic2, b1, b2); -+ _makeBaryonicClusters(*cit, *baryonic1, *baryonic2, b1, b2); - - *cit = b1; - *baryonic1 = b2; -@@ -263,9 +305,9 @@ - } - - // Normal 2->2 Colour reconnection -- if ( !isBaryonicCandidate -- && UseRandom::rnd() < _preco ) { -- auto reconnected = _reconnectBaryonic(*cit, *candidate); -+ if (!isBaryonicCandidate -+ && UseRandom::rnd() < _preco) { -+ auto reconnected = _reconnect(*cit, *candidate); - *cit = reconnected.first; - *candidate = reconnected.second; - } -@@ -274,18 +316,254 @@ - // create a new vector of clusters except for the ones which are "deleted" during - // baryonic reconnection - ClusterVector clustervector; -- for ( const auto & cluster : newcv ) -- if ( find(deleted.begin(), -- deleted.end(), cluster) == deleted.end() ) -+ for (const auto & cluster : newcv) -+ if (find(deleted.begin(), -+ deleted.end(), cluster) == deleted.end()) - clustervector.push_back(cluster); - -- swap(cv,clustervector); -+ swap(cv, clustervector); - - - - } - - -+bool ColourReconnector::_clustersFarApart( const std::vector & clu ) const { -+ int Ncl=clu.size(); -+ assert(Ncl<=3); -+ if (Ncl==1) { -+ return false; -+ } else if (Ncl==2) { -+ // TODO: keep turned off all until things are more clear -+ // BUG: space-time difference compared to _maxDistance -+ // if (((*clu[0])->vertex()-(*clu[1])->vertex()).m() >_maxDistance) return true; -+ // Causal selection if desired -+ // if (((*clu[0])->vertex()-(*clu[1])->vertex()).m() >ZERO) return true; -+ } else if (Ncl==3) { -+ // TODO: keep turned off all until things are more clear -+ // BUG: space-time difference compared to _maxDistance -+ // if (((*clu[0])->vertex()-(*clu[1])->vertex()).m()> _maxDistance) return true; -+ // Causal selection if desired -+ // if (((*clu[0])->vertex()-(*clu[1])->vertex()).m()> ZERO) return true; -+ // if (((*clu[1])->vertex()-(*clu[2])->vertex()).m()> _maxDistance) return true; -+ // if (((*clu[1])->vertex()-(*clu[2])->vertex()).m()> ZERO) return true; -+ // if (((*clu[0])->vertex()-(*clu[2])->vertex()).m()> _maxDistance) return true; -+ // if (((*clu[0])->vertex()-(*clu[2])->vertex()).m()> ZERO) return true; -+ } -+ -+ return false; -+} -+ -+ -+ -+void ColourReconnector::_doReco2BeamClusters(ClusterVector & cv) const { -+ // try other option -+ PPtr p1Di=(cv[0])->colParticle(); -+ PPtr p2Di=(cv[1])->colParticle(); -+ -+ PPtr p1Q=(cv[0])->antiColParticle(); -+ PPtr p2Q=(cv[1])->antiColParticle(); -+ -+ double min_dist=_displacement(p1Di,p1Q)+_displacement(p2Di,p2Q); -+ -+ if ((_displacement(p1Di,p2Q)+_displacement(p1Di,p2Q)) sel; -+ -+ unsigned number_of_tries = _stepFactor*cv.size()*cv.size(); -+ if (number_of_tries<1) number_of_tries=1; -+ -+ long (*p_irnd)(long) = UseRandom::irnd; -+ for (unsigned reconnections_tries = 0; reconnections_tries < number_of_tries; reconnections_tries++) { -+ num_meson = 0; -+ num_baryon = 0; -+ -+ // flag if we are able to find a suitable combinations of clusters -+ bool _found = false; -+ -+ // Shuffle list of clusters to avoid systematic bias in cluster selection -+ random_shuffle(newcv.begin(), newcv.end(), p_irnd); -+ -+ // loop over clustervector to find CR candidates -+ for (CluVecIt cit = newcv.begin(); cit != newcv.end(); ++cit) { -+ -+ // skip the clusters to be deleted later from 3->2 cluster CR -+ if (find(deleted.begin(), deleted.end(), *cit) != deleted.end()) continue; -+ -+ // avoid clusters already containing diuarks -+ if (hasDiquark(cit)) continue; -+ -+ // add to selection -+ sel.push_back(cit); -+ -+ if (_clustersFarApart(sel)) { -+ // reject far appart CR -+ // TODO: after discussion maybe to be omitted -+ sel.pop_back(); -+ continue; -+ } -+ -+ bool isMeson=((*cit)->numComponents() == 2); -+ -+ if ( isMeson && (num_meson ==0|| num_meson==1) && num_baryon ==0) { -+ num_meson++; -+ /** -+ * now we habe either 1 or 2 mesonic clusters and have to continue -+ */ -+ continue; -+ } else if ( isMeson && (num_baryon == 1 || num_meson ==2)) { -+ num_meson++; -+ _found = true; -+ /** -+ * we have either 3 mesonic or 1 mesonic and 1 baryonic cluster -+ * and try to colour reconnect -+ */ -+ break; -+ } else if (num_baryon ==0 && num_meson==0) { -+ num_baryon++; -+ /** -+ * now we have 1 baryonic cluster and have to continue -+ */ -+ continue; -+ } else if (num_meson == 2) { -+ /** -+ * we already have 2 mesonic clusters and dont want a baryonic one -+ * since this is an invalid selection -+ */ -+ // remove previously added cluster -+ sel.pop_back(); -+ continue; -+ } else { -+ num_baryon++; -+ _found = true; -+ /** -+ * now we have either 2 baryonic clusters or 1 mesonic and 1 baryonic cluster -+ * and try to colour reconnect -+ */ -+ break; -+ } -+ } -+ -+ // added for more efficent rejection if some reco probabilities are 0 -+ if ( _found ) { -+ -+ // reject MBtoMB candidates if _precoMB_MB=0 -+ if ( _precoMB_MB == 0 && (num_baryon == 1 && num_meson == 1) ) { -+ _found=false; -+ } -+ -+ // reject BbarBto3M candidates if _precoBbarB_3M=0 -+ if ( _precoBbarB_3M== 0 && num_baryon == 2 ) { -+ bool isBbarBto3Mcandidate=( -+ (*sel[0])->particle(0)->hasColour() && (*sel[1])->particle(0)->hasColour(true) ) -+ || ( (*sel[0])->particle(0)->hasColour(true) && (*sel[1])->particle(0)->hasColour() ); -+ -+ if ( isBbarBto3Mcandidate) _found=false; -+ } -+ -+ // reject 2Bto2B candidates if _preco2B_2B=0 -+ if ( _preco2B_2B == 0 && num_baryon == 2 ) { -+ bool is2Bto2Bcandidate=( -+ (*sel[0])->particle(0)->hasColour() && (*sel[1])->particle(0)->hasColour() ) -+ || ( (*sel[0])->particle(0)->hasColour(true) && (*sel[1])->particle(0)->hasColour(true) ); -+ -+ if ( is2Bto2Bcandidate ) _found=false; -+ } -+ } -+ // were we able to find a combination? -+ if (_found==false) { -+ // clear the selection if we did not find a valid set of clusters -+ sel.erase(sel.begin(), sel.end()); -+ continue; -+ } -+ assert(sel.size()<4); -+ assert(sel.size()>1); -+ -+ string kind_of_reco = ""; -+ int reco_info[3]; -+ -+ // find best CR option for the selection -+ _findbestreconnectionoption(sel, num_baryon, kind_of_reco, reco_info); -+ -+ if (kind_of_reco == "") { -+ // no reconnection was found -+ sel.erase(sel.begin(), sel.end()); -+ continue; -+ } else if (kind_of_reco == "3Mto3M" && UseRandom::rnd() < _preco3M_3M) { -+ // 3Mto3M colour reconnection -+ auto reconnected = _reconnect3Mto3M(*sel[0], *sel[1], *sel[2], -+ reco_info); -+ (*sel[0]) = std::get<0>(reconnected); -+ (*sel[1]) = std::get<1>(reconnected); -+ (*sel[2]) = std::get<2>(reconnected); -+ } else if (kind_of_reco=="2Bto3M" && UseRandom::rnd() < _precoBbarB_3M) { -+ // antibaryonic and baryonic to 3 mesonic reconnecion -+ auto reconnected = _reconnectBbarBto3M(*sel[0], *sel[1], -+ reco_info[0], reco_info[1], reco_info[2]); -+ (*sel[0]) = std::get<0>(reconnected); -+ (*sel[1]) = std::get<1>(reconnected); -+ newcv.push_back(std::get<2>(reconnected)); -+ } else if (kind_of_reco=="3Mto2B" && UseRandom::rnd() < _preco3M_BBbar) { -+ // 3 mesonic to antibaryonic and baryonic reconnection -+ ClusterPtr b1, b2; -+ _makeBaryonicClusters(*sel[0], *sel[1], *sel[2], b1, b2); -+ (*sel[0]) = b1; -+ (*sel[1]) = b2; -+ deleted.push_back(*sel[2]); -+ } else if (kind_of_reco=="2Bto2B" && UseRandom::rnd() < _preco2B_2B) { -+ // 2 (anti)baryonic to 2 (anti)baryonic reconnection -+ auto reconnected = _reconnect2Bto2B(*sel[0], *sel[1], -+ reco_info[0], reco_info[1]); -+ (*sel[0]) = reconnected.first; -+ (*sel[1]) = reconnected.second; -+ } else if (kind_of_reco=="MBtoMB" && UseRandom::rnd() < _precoMB_MB) { -+ // (anti)baryonic and mesonic to (anti)baryonic and mesonic reconnection -+ auto reconnected = _reconnectMBtoMB(*sel[0], *sel[1], -+ reco_info[0]); -+ (*sel[0]) = reconnected.first; -+ (*sel[1]) = reconnected.second; -+ } -+ // erase the sel-vector -+ sel.erase(sel.begin(), sel.end()); -+ } -+ -+ // write to clustervector new CR'd clusters and deleting -+ // all deleted clusters -+ ClusterVector clustervector; -+ for (const auto & cluster : newcv) -+ if (find(deleted.begin(), deleted.end(), cluster) == deleted.end()) -+ clustervector.push_back(cluster); -+ -+ swap(cv, clustervector); -+} - - namespace { - -@@ -314,6 +592,84 @@ - } - - -+void ColourReconnector::_findbestreconnectionoption(std::vector & cls, const unsigned & baryonic, -+ string & kind_of_reco, int (&reco_info)[3]) const { -+ double min_displacement; -+ if (baryonic==0) { -+ // case with 3 mesonic clusters -+ assert(cls.size()==3); -+ -+ // calculate the initial displacement sum -+ min_displacement = _mesonToBaryonFactor * _displacement((*cls[0])->particle(0), (*cls[0])->particle(1)); -+ min_displacement += _mesonToBaryonFactor * _displacement((*cls[1])->particle(0), (*cls[1])->particle(1)); -+ min_displacement += _mesonToBaryonFactor * _displacement((*cls[2])->particle(0), (*cls[2])->particle(1)); -+ -+ // find best CR reco_info and kind_of_reco -+ _3MtoXreconnectionfinder(cls, -+ reco_info[0], reco_info[1], reco_info[2], min_displacement, kind_of_reco); -+ /** -+ * kind_of_reco either "3Mto3M" or "3Mto2B" (or "" if no better configuration is found) -+ * case 3Mto3M: the coloured particle of the i-th cluster forms a new cluster with the -+ * antiparticle of the reco_info[i]-th cluster -+ * case 3MtoBbarB: all 3 (anti)coloured particle form a new (anti)baryonic cluster -+ */ -+ } else if (baryonic == 1) { -+ // case 1 baryonic and 1 mesonic cluster -+ assert(cls.size()==2); -+ -+ // make mesonic cluster always the cls[0] -+ if ((*cls[0])->numComponents() == 3) { -+ ClusterPtr zw = *cls[0]; -+ *cls[0] = *cls[1]; -+ *cls[1] = zw; -+ } -+ -+ // calculate the initial displacement sum -+ min_displacement = _mesonToBaryonFactor *_displacement((*cls[0])->particle(0), (*cls[0])->particle(1)); -+ min_displacement += _displacementBaryonic((*cls[1])->particle(0), (*cls[1])->particle(1), (*cls[1])->particle(2)); -+ -+ // find best CR reco_info and kind_of_reco -+ _BMtoBMreconnectionfinder(*cls[0], *cls[1], -+ reco_info[0], min_displacement, kind_of_reco); -+ /** -+ * reco_info[0] is the index of the (anti)quarks of the baryonic cluster cls[1], which should -+ * be swapped with the (anti)quarks of the mesonic cluster cls[0] -+ */ -+ -+ } else { -+ assert(baryonic==2); -+ assert(cls.size()==2); -+ -+ // calculate the initial displacement sum -+ min_displacement = _displacementBaryonic((*cls[0])->particle(0), (*cls[0])->particle(1), (*cls[0])->particle(2)); -+ min_displacement += _displacementBaryonic((*cls[1])->particle(0), (*cls[1])->particle(1), (*cls[1])->particle(2)); -+ -+ // case 2 (anti)baryonic clusters to 2 other (anti)baryonic clusters -+ if ( ( (*cls[0])->particle(0)->hasColour() && (*cls[1])->particle(0)->hasColour() ) -+ || ( (*cls[0])->particle(0)->hasColour(true) && (*cls[1])->particle(0)->hasColour(true) ) ) { -+ // find best CR reco_info and kind_of_reco -+ _2Bto2BreconnectionFinder(*cls[0], *cls[1], -+ reco_info[0], reco_info[1], min_displacement, kind_of_reco); -+ /** -+ * swap the reco_info[0]-th particle of the first cluster in the vector with the -+ * reco_info[1]-th particle of the second cluster -+ */ -+ } else { -+ // case 1 baryonic and 1 antibaryonic cluster to 3 mesonic clusters -+ -+ // find best CR reco_info and kind_of_reco -+ _BbarBto3MreconnectionFinder(*cls[0], *cls[1], -+ reco_info[0], reco_info[1], reco_info[2], min_displacement, kind_of_reco); -+ /** -+ * the i-th particle of the first cluster form a new mesonic cluster with the -+ * reco_info[i]-th particle of the second cluster -+ */ -+ } -+ } -+ return; -+} -+ -+ - CluVecIt ColourReconnector::_findPartnerBaryonic( - CluVecIt cl, ClusterVector & cv, - bool & baryonicCand, -@@ -349,7 +705,6 @@ - p1col.boost(boostv); - p1anticol.boost(boostv); - -- - for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) { - //avoid looping over clusters containing diquarks - if ( hasDiquark(cit) ) continue; -@@ -364,8 +719,10 @@ - continue; - - // stop it putting far apart clusters together -- if ( ( (**cl).vertex()-(**cit).vertex() ).m() >_maxDistance ) -- continue; -+ // BUG: space-time difference compared to _maxDistance -+ // if (((**cl).vertex()-(**cit).vertex()).m() >_maxDistance) -+ // Causal selection if desired -+ // if (((**cl).vertex()-(**cit).vertex()).m() >ZERO) continue; - - const bool Colour8 = - _isColour8( (*cl)->colParticle(), (*cit)->antiColParticle() ) -@@ -412,7 +769,7 @@ - // first candidate gets here. If second baryonic candidate has higher Ysum than the first - // one, the second candidate becomes the first one and the first the second. - if (sumrap > maxsum) { -- if(maxsum != 0){ -+ if (maxsum != 0) { - baryonic2 = baryonic1; - baryonic1 = cit; - bcand = true; -@@ -431,7 +788,7 @@ - - } - -- if(bcand == true){ -+ if (bcand == true) { - baryonicCand = true; - } - -@@ -447,7 +804,7 @@ - for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) { - - // don't even look at original cluster -- if(cit==cl) continue; -+ if (cit==cl) continue; - - // don't allow colour octet clusters - if ( _isColour8( (*cl)->colParticle(), -@@ -458,10 +815,13 @@ - } - - // stop it putting beam remnants together -- if((*cl)->isBeamCluster() && (*cit)->isBeamCluster()) continue; -+ if ((*cl)->isBeamCluster() && (*cit)->isBeamCluster()) continue; - - // stop it putting far apart clusters together -- if(((**cl).vertex()-(**cit).vertex()).m()>_maxDistance) continue; -+ // BUG: space-time difference compared to _maxDistance -+ // if (((**cl).vertex()-(**cit).vertex()).m()>_maxDistance) continue; -+ // Causal selection if desired -+ // if (((**cl).vertex()-(**cit).vertex()).m()>ZERO) continue; - - // momenta of the old clusters - Lorentz5Momentum p1 = (*cl)->colParticle()->momentum() + -@@ -493,7 +853,7 @@ - ClusterPtr &c1, ClusterPtr &c2, - ClusterPtr &c3, - ClusterPtr &newcluster1, -- ClusterPtr &newcluster2) const{ -+ ClusterPtr &newcluster2) const { - - //make sure they all have 2 components - assert(c1->numComponents()==2); -@@ -522,6 +882,105 @@ - } - - pair -+ColourReconnector::_reconnect2Bto2B(ClusterPtr &c1, ClusterPtr &c2, const int s1, const int s2) const { -+ -+ // form the first new cluster -+ -+ // separate the quarks from their original cluster -+ c1->particleB((s1+1)%3)->abandonChild(c1); -+ c1->particleB((s1+2)%3)->abandonChild(c1); -+ c2->particleB(s2)->abandonChild(c2); -+ -+ // now the new cluster -+ ClusterPtr newCluster1 = new_ptr(Cluster(c1->particleB((s1+1)%3), c1->particleB((s1+2)%3), c2->particleB(s2))); -+ -+ c1->particleB((s1+1)%3)->addChild(newCluster1); -+ c1->particleB((s1+2)%3)->addChild(newCluster1); -+ c2->particleB(s2)->addChild(newCluster1); -+ -+ // set new vertex -+ newCluster1->setVertex(LorentzPoint()); -+ -+ // set beam remnants for new cluster -+ if (c1->isBeamRemnant((s1+1)%3)) newCluster1->setBeamRemnant(0, true); -+ if (c1->isBeamRemnant((s1+2)%3)) newCluster1->setBeamRemnant(1, true); -+ if (c2->isBeamRemnant(s2)) newCluster1->setBeamRemnant(2, true); -+ -+ // for the second cluster same procedure -+ c2->particleB((s2+1)%3)->abandonChild(c2); -+ c2->particleB((s2+2)%3)->abandonChild(c2); -+ c1->particleB(s1)->abandonChild(c1); -+ -+ ClusterPtr newCluster2 = new_ptr(Cluster(c2->particleB((s2+1)%3), c2->particleB((s2+2)%3), c1->particleB(s1))); -+ -+ c2->particleB((s2+1)%3)->addChild(newCluster2); -+ c2->particleB((s2+2)%3)->addChild(newCluster2); -+ c1->particleB(s1)->addChild(newCluster2); -+ -+ newCluster2->setVertex(LorentzPoint()); -+ -+ if (c2->isBeamRemnant((s2+1)%3)) newCluster2->setBeamRemnant(0, true); -+ if (c2->isBeamRemnant((s2+2)%3)) newCluster2->setBeamRemnant(1, true); -+ if (c1->isBeamRemnant(s1)) newCluster2->setBeamRemnant(2, true); -+ -+ return pair (newCluster1, newCluster2); -+} -+ -+ -+std::tuple -+ColourReconnector::_reconnectBbarBto3M(ClusterPtr & c1, ClusterPtr & c2, const int s0, const int s1, const int s2) const { -+ // make sure they all have 3 components -+ assert(c1->numComponents()==3); -+ assert(c2->numComponents()==3); -+ -+ // first Cluster -+ c1->particleB(0)->abandonChild(c1); -+ c2->particleB(s0)->abandonChild(c2); -+ -+ ClusterPtr newCluster1 = new_ptr(Cluster(c1->particleB(0), c2->particleB(s0))); -+ -+ c1->particleB(0)->addChild(newCluster1); -+ c2->particleB(s0)->addChild(newCluster1); -+ -+ // set new vertex -+ newCluster1->setVertex(0.5*(c1->particleB(0)->vertex() + c2->particleB(s0)->vertex())); -+ -+ // set beam remnants for new cluster -+ if (c1->isBeamRemnant(0)) newCluster1->setBeamRemnant(0, true); -+ if (c2->isBeamRemnant(s0)) newCluster1->setBeamRemnant(1, true); -+ -+ // same for second cluster -+ c1->particleB(1)->abandonChild(c1); -+ c2->particleB(s1)->abandonChild(c2); -+ -+ ClusterPtr newCluster2 = new_ptr(Cluster(c1->particleB(1), c2->particleB(s1))); -+ -+ c1->particleB(1)->addChild(newCluster2); -+ c2->particleB(s1)->addChild(newCluster2); -+ -+ newCluster2->setVertex(0.5*(c1->particleB(1)->vertex() + c2->particleB(s1)->vertex())); -+ -+ if (c1->isBeamRemnant(1)) newCluster2->setBeamRemnant(0, true); -+ if (c2->isBeamRemnant(s1)) newCluster2->setBeamRemnant(1, true); -+ -+ // same for third cluster -+ c1->particleB(2)->abandonChild(c1); -+ c2->particleB(s2)->abandonChild(c2); -+ -+ ClusterPtr newCluster3 = new_ptr(Cluster(c1->particleB(2), c2->particleB(s2))); -+ -+ c1->particleB(2)->addChild(newCluster3); -+ c2->particleB(s2)->addChild(newCluster3); -+ -+ newCluster3->setVertex(0.5*(c1->particleB(2)->vertex() + c2->particleB(s2)->vertex())); -+ -+ if (c1->isBeamRemnant(2)) newCluster3->setBeamRemnant(0, true); -+ if (c2->isBeamRemnant(s2)) newCluster3->setBeamRemnant(1, true); -+ -+ return std::tuple (newCluster1, newCluster2, newCluster3); -+} -+ -+pair - ColourReconnector::_reconnect(ClusterPtr &c1, ClusterPtr &c2) const { - - // choose the other possibility to form two clusters from the given -@@ -530,7 +989,7 @@ - assert(c1->numComponents()==2); - assert(c2->numComponents()==2); - int c1_col(-1),c1_anti(-1),c2_col(-1),c2_anti(-1); -- for(unsigned int ix=0;ix<2;++ix) { -+ for(unsigned int ix=0; ix<2; ++ix) { - if (c1->particle(ix)->hasColour(false)) c1_col = ix; - else if(c1->particle(ix)->hasColour(true )) c1_anti = ix; - if (c2->particle(ix)->hasColour(false)) c2_col = ix; -@@ -538,50 +997,8 @@ - } - assert(c1_col>=0&&c2_col>=0&&c1_anti>=0&&c2_anti>=0); - -- ClusterPtr newCluster1 -- = new_ptr( Cluster( c1->colParticle(), c2->antiColParticle() ) ); -- -- newCluster1->setVertex(0.5*( c1->colParticle()->vertex() + -- c2->antiColParticle()->vertex() )); -- -- if(c1->isBeamRemnant(c1_col )) newCluster1->setBeamRemnant(0,true); -- if(c2->isBeamRemnant(c2_anti)) newCluster1->setBeamRemnant(1,true); -- -- ClusterPtr newCluster2 -- = new_ptr( Cluster( c2->colParticle(), c1->antiColParticle() ) ); -- -- newCluster2->setVertex(0.5*( c2->colParticle()->vertex() + -- c1->antiColParticle()->vertex() )); -- -- if(c2->isBeamRemnant(c2_col )) newCluster2->setBeamRemnant(0,true); -- if(c1->isBeamRemnant(c1_anti)) newCluster2->setBeamRemnant(1,true); -- -- return pair (newCluster1, newCluster2); --} -- -- -- -- -- --pair --ColourReconnector::_reconnectBaryonic(ClusterPtr &c1, ClusterPtr &c2) const { -- -- // choose the other possibility to form two clusters from the given -- // constituents -- -- assert(c1->numComponents()==2); -- assert(c2->numComponents()==2); -- int c1_col(-1),c1_anti(-1),c2_col(-1),c2_anti(-1); -- for(unsigned int ix=0;ix<2;++ix) { -- if (c1->particle(ix)->hasColour(false)) c1_col = ix; -- else if(c1->particle(ix)->hasColour(true )) c1_anti = ix; -- if (c2->particle(ix)->hasColour(false)) c2_col = ix; -- else if(c2->particle(ix)->hasColour(true )) c2_anti = ix; -- } -- assert(c1_col>=0&&c2_col>=0&&c1_anti>=0&&c2_anti>=0); -- --c1->colParticle()->abandonChild(c1); --c2->antiColParticle()->abandonChild(c2); -+ c1->colParticle()->abandonChild(c1); -+ c2->antiColParticle()->abandonChild(c2); - - ClusterPtr newCluster1 - = new_ptr( Cluster( c1->colParticle(), c2->antiColParticle() ) ); -@@ -589,8 +1006,11 @@ - c1->colParticle()->addChild(newCluster1); - c2->antiColParticle()->addChild(newCluster1); - -- newCluster1->setVertex(0.5*( c1->colParticle()->vertex() + -- c2->antiColParticle()->vertex() )); -+ /* -+ * TODO: Questionable setting of the vertex -+ * */ -+ newCluster1->setVertex(0.5*(c1->colParticle()->vertex() + -+ c2->antiColParticle()->vertex())); - - if(c1->isBeamRemnant(c1_col )) newCluster1->setBeamRemnant(0,true); - if(c2->isBeamRemnant(c2_anti)) newCluster1->setBeamRemnant(1,true); -@@ -604,8 +1024,11 @@ - c1->antiColParticle()->addChild(newCluster2); - c2->colParticle()->addChild(newCluster2); - -- newCluster2->setVertex(0.5*( c2->colParticle()->vertex() + -- c1->antiColParticle()->vertex() )); -+ /* -+ * TODO: Questionable setting of the vertex -+ * */ -+ newCluster2->setVertex(0.5*(c2->colParticle()->vertex() + -+ c1->antiColParticle()->vertex())); - - if(c2->isBeamRemnant(c2_col )) newCluster2->setBeamRemnant(0,true); - if(c1->isBeamRemnant(c1_anti)) newCluster2->setBeamRemnant(1,true); -@@ -613,9 +1036,289 @@ - return pair (newCluster1, newCluster2); - } - -+std::tuple -+ColourReconnector::_reconnect3Mto3M(ClusterPtr & c1, ClusterPtr & c2, ClusterPtr & c3, const int infos [3]) const { -+ // check if mesonic clusters -+ assert(c1->numComponents()==2); -+ assert(c2->numComponents()==2); -+ assert(c3->numComponents()==2); -+ -+ ClusterVector oldclusters = {c1, c2, c3}; -+ ClusterVector newclusters; -+ -+ for (int i=0; i<3; i++) { -+ -+ int c1_col=-1; -+ int c2_anticol=-1; -+ -+ // get which index is coloured and which anticolour -+ for (unsigned int ix=0; ix<2; ++ix) { -+ if (oldclusters[i]->particle(ix)->hasColour(false)) c1_col = ix; -+ if (oldclusters[infos[i]]->particle(ix)->hasColour(true)) c2_anticol = ix; -+ } -+ -+ assert(c1_col>=0); -+ assert(c2_anticol>=0); -+ -+ oldclusters[i]->colParticle()->abandonChild(oldclusters[i]); -+ oldclusters[infos[i]]->antiColParticle()->abandonChild(oldclusters[infos[i]]); -+ -+ // form new cluster -+ ClusterPtr newCluster = new_ptr(Cluster(oldclusters[i]->colParticle(), oldclusters[infos[i]]->antiColParticle())); -+ -+ oldclusters[i]->colParticle()->addChild(newCluster); -+ oldclusters[infos[i]]->antiColParticle()->addChild(newCluster); -+ -+ // set new vertex -+ newCluster->setVertex(0.5*(oldclusters[i]->colParticle()->vertex() + -+ oldclusters[infos[i]]->antiColParticle()->vertex())); -+ -+ // set beam remnants for new cluster -+ if (oldclusters[i]->isBeamRemnant(c1_col)) newCluster->setBeamRemnant(0, true); -+ if (oldclusters[infos[i]]->isBeamRemnant(c2_anticol)) newCluster->setBeamRemnant(1, true); -+ newclusters.push_back(newCluster); -+ } -+ return std::tuple (newclusters[0], newclusters[1], newclusters[2]); -+} -+ -+ -+pair -+ColourReconnector::_reconnectMBtoMB(ClusterPtr & c1, ClusterPtr & c2, const int s0) const { -+ // make c1 the mesonic cluster -+ if (c1->numComponents()==2) { -+ assert(c2->numComponents()==3); -+ } else { -+ return _reconnectMBtoMB(c2,c1,s0); -+ } -+ -+ int c1_col=-1; -+ int c1_anti=-1; -+ // get which index is coloured and which anticolour -+ for (unsigned int ix=0; ix<2; ++ix) { -+ if (c1->particle(ix)->hasColour(false)) c1_col = ix; -+ else if (c1->particle(ix)->hasColour(true)) c1_anti = ix; -+ -+ } -+ assert(c1_col>=0); -+ assert(c1_anti>=0); -+ -+ // pointers for the new clusters -+ ClusterPtr newCluster1; -+ ClusterPtr newCluster2; -+ if (c2->particle(0)->hasColour()==true) { -+ // first case: we have a baryonic clusters -+ -+ // first make the new mesonic cluster -+ c1->antiColParticle()->abandonChild(c1); -+ c2->particleB(s0)->abandonChild(c2); -+ -+ newCluster1 = new_ptr(Cluster(c1->antiColParticle(), c2->particleB(s0))); -+ -+ c1->antiColParticle()->addChild(newCluster1); -+ c2->particleB(s0)->addChild(newCluster1); -+ -+ // set new vertex -+ newCluster1->setVertex(0.5*(c1->antiColParticle()->vertex() + -+ c2->particleB(s0)->vertex())); -+ -+ // set beam remnants for new cluster -+ if (c1->isBeamRemnant(c1_anti)) newCluster1->setBeamRemnant(0, true); -+ if (c2->isBeamRemnant(s0)) newCluster1->setBeamRemnant(1, true); -+ -+ // then the baryonic one -+ c1->colParticle()->abandonChild(c1); -+ c2->particleB((s0+1)%3)->abandonChild(c2); -+ c2->particleB((s0+2)%3)->abandonChild(c2); -+ -+ newCluster2 = new_ptr(Cluster(c1->colParticle(), c2->particleB((s0+1)%3), c2->particleB((s0+2)%3))); -+ -+ c1->colParticle()->addChild(newCluster2); -+ c2->particleB((s0+1)%3)->addChild(newCluster2); -+ c2->particleB((s0+2)%3)->addChild(newCluster2); -+ -+ // set new vertex -+ newCluster2->setVertex(LorentzPoint()); -+ } else { -+ // second case we have an antibaryonic cluster -+ -+ // first make the new mesonic cluster -+ c1->colParticle()->abandonChild(c1); -+ c2->particleB(s0)->abandonChild(c2); -+ -+ newCluster1 = new_ptr(Cluster(c1->colParticle(), c2->particleB(s0))); -+ -+ c1->colParticle()->addChild(newCluster1); -+ c2->particleB(s0)->addChild(newCluster1); -+ -+ // set new vertex -+ newCluster1->setVertex(0.5*(c1->colParticle()->vertex() + -+ c2->particleB(s0)->vertex())); -+ -+ // set beam remnants for new cluster -+ if (c1->isBeamRemnant(c1_col)) newCluster1->setBeamRemnant(0, true); -+ if (c2->isBeamRemnant(s0)) newCluster1->setBeamRemnant(1, true); -+ -+ // then the baryonic one -+ c1->antiColParticle()->abandonChild(c1); -+ c2->particleB((s0+1)%3)->abandonChild(c2); -+ c2->particleB((s0+2)%3)->abandonChild(c2); -+ -+ newCluster2 = new_ptr(Cluster(c1->antiColParticle(), c2->particleB((s0+1)%3), c2->particleB((s0+2)%3))); -+ -+ c1->antiColParticle()->addChild(newCluster2); -+ c2->particleB((s0+1)%3)->addChild(newCluster2); -+ c2->particleB((s0+2)%3)->addChild(newCluster2); -+ -+ // set new vertex -+ newCluster2->setVertex(LorentzPoint()); -+ } -+ return pair (newCluster1, newCluster2); -+} -+ -+void ColourReconnector::_2Bto2BreconnectionFinder(ClusterPtr & c1, ClusterPtr & c2, -+ int & bswap1, int & bswap2, double min_displ_sum, string & kind_of_reco) const { -+ double tmp_delta; -+ for (int i=0; i<3; i++) { -+ for (int j=0; j<3; j++) { -+ // try swapping particle i of c1 with particle j of c2 -+ tmp_delta = _displacementBaryonic(c2->particle(j), c1->particle((i+1)%3), c1->particle((i+2)%3)); -+ tmp_delta += _displacementBaryonic(c1->particle(i), c2->particle((j+1)%3), c2->particle((j+2)%3)); -+ -+ if (tmp_delta < min_displ_sum) { -+ // if minimal displacement select the 2Bto2B CR option -+ min_displ_sum = tmp_delta; -+ bswap1 = i; -+ bswap2 = j; -+ kind_of_reco = "2Bto2B"; -+ } -+ } -+ } -+ -+} -+ -+void ColourReconnector::_BbarBto3MreconnectionFinder(ClusterPtr & c1, ClusterPtr & c2, int & mswap0, int & mswap1, int & mswap2, -+ double min_displ_sum, string & kind_of_reco) const { -+ double pre_tmp_delta; -+ double tmp_delta; -+ for (int p1=0; p1 <3; p1++) { -+ // make sure not to form a mesonic octet -+ if (_isColour8(c1->particle(0), c2->particle(p1))) continue; -+ -+ pre_tmp_delta = _displacement(c1->particle(0), c2->particle(p1)); -+ for (int p2=1; p2<3; p2++) { -+ -+ // make sure not to form a mesonic octet -+ if (_isColour8(c1->particle(1), c2->particle((p1+p2)%3))) continue; -+ if (_isColour8(c1->particle(2), c2->particle(3-p1-((p1+p2)%3)))) continue; -+ -+ tmp_delta = pre_tmp_delta + _displacement(c1->particle(1), c2->particle((p1+p2)%3)); -+ tmp_delta += _displacement(c1->particle(2), c2->particle(3-p1-((p1+p2)%3))); -+ -+ // factor _mesonToBaryonFactor to compare Baryonic an mesonic cluster -+ tmp_delta *=_mesonToBaryonFactor; -+ -+ if (tmp_delta < min_displ_sum) { -+ // if minimal displacement select the 2Bto3M CR option -+ min_displ_sum = tmp_delta; -+ mswap0 = p1; -+ mswap1 = (p1+p2)%3; -+ mswap2 = 3-p1-((p1+p2)%3); -+ kind_of_reco = "2Bto3M"; -+ -+ } -+ } -+ } -+} -+ -+void ColourReconnector::_BMtoBMreconnectionfinder(ClusterPtr & c1, ClusterPtr & c2, int & swap, double min_displ_sum, -+ string & kind_of_reco) const { -+ assert(c1->numComponents()==2); -+ assert(c2->numComponents()==3); -+ double tmp_displ = 0; -+ for (int i=0; i<3; i++) { -+ // Differ if the second cluster is baryonic or antibaryonic -+ if (c2->particle(0)->hasColour()) { -+ // c2 is baryonic -+ -+ // veto mesonic octets -+ if (_isColour8(c2->particle(i), c1->antiColParticle())) continue; -+ -+ // factor _mesonToBaryonFactor to compare Baryonic an mesonic cluster -+ tmp_displ = _mesonToBaryonFactor * _displacement(c2->particle(i), c1->antiColParticle()); -+ tmp_displ += _displacementBaryonic(c1->colParticle(), c2->particle((i+1)%3), c2->particle((i+2)%3)); -+ } else { -+ // c2 is antibaryonic -+ -+ // veto mesonic octets -+ if (_isColour8(c2->particle(i), c1->colParticle())) continue; -+ -+ // factor _mesonToBaryonFactor to compare Baryonic an mesonic cluster -+ tmp_displ = _mesonToBaryonFactor * _displacement(c2->particle(i), c1->colParticle()); -+ tmp_displ *= _displacementBaryonic(c1->antiColParticle(), c2->particle((i+1)%3), c2->particle((i+2)%3)); -+ } -+ if (tmp_displ < min_displ_sum) { -+ // if minimal displacement select the MBtoMB CR option -+ min_displ_sum = tmp_displ; -+ swap = i; -+ kind_of_reco = "MBtoMB"; -+ } -+ } -+ return; -+} -+ -+void ColourReconnector::_3MtoXreconnectionfinder(std::vector & cv, int & swap0, int & swap1, -+ int & swap2, double min_displ_sum, string & kind_of_reco) const { -+ // case of 3M->BbarB CR -+ double _tmp_displ; -+ _tmp_displ = _displacementBaryonic((*cv[0])->colParticle(), (*cv[1])->colParticle(), (*cv[2])->colParticle()); -+ _tmp_displ += _displacementBaryonic((*cv[0])->antiColParticle(), (*cv[1])->antiColParticle(), (*cv[2])->antiColParticle()); -+ if (_tmp_displ < min_displ_sum) { -+ // if minimal displacement select the 3Mto2B CR option -+ kind_of_reco = "3Mto2B"; -+ min_displ_sum = _tmp_displ; -+ } -+ // case for 3M->3M CR -+ /** -+ * if 3Mto3M reco probability (_preco3M_3M) is 0 we skip this loop -+ * since no 3Mto3M CR shall be performed -+ */ -+ int i,j; -+ int i1,i2,i3; -+ for (i = 0; _preco3M_3M && i<3; i++) { -+ // veto mesonic octets -+ if (_isColour8((*cv[0])->colParticle(), (*cv[i])->antiColParticle())) continue; -+ -+ // factor _mesonToBaryonFactor to compare baryonic an mesonic cluster -+ _tmp_displ = _mesonToBaryonFactor * _displacement((*cv[0])->colParticle(), (*cv[i])->antiColParticle()); -+ for (j=1; j<3; j++) { -+ // i1, i2, i3 are pairwise distinct -+ i1=i; -+ i2=((j+i)%3); -+ if (i1==0 && i2==1) continue; -+ i3=(3-i-((j+i)%3)); -+ -+ // veto mesonic octets -+ if (_isColour8((*cv[1])->colParticle(), (*cv[i2])->antiColParticle())) continue; -+ if (_isColour8((*cv[2])->colParticle(), (*cv[i3])->antiColParticle())) continue; -+ -+ _tmp_displ += _mesonToBaryonFactor * _displacement((*cv[1])->colParticle(), (*cv[i2])->antiColParticle()); -+ _tmp_displ += _mesonToBaryonFactor * _displacement((*cv[2])->colParticle(), (*cv[i3])->antiColParticle()); -+ -+ if (_tmp_displ < min_displ_sum) { -+ // if minimal displacement select the 3Mto3M CR option -+ kind_of_reco = "3Mto3M"; -+ min_displ_sum = _tmp_displ; -+ swap0 = i1; -+ swap1 = i2; -+ swap2 = i3; -+ } -+ } -+ } -+} -+ - - pair ColourReconnector::_shuffle -- (const PVector & q, const PVector & aq, unsigned maxtries) const { -+(const PVector & q, const PVector & aq, unsigned maxtries) const { - - const size_t nclusters = q.size(); - assert (nclusters > 1); -@@ -628,7 +1331,9 @@ - do { - // find two different random integers in the range [0, nclusters) - i = UseRandom::irnd( nclusters ); -- do { j = UseRandom::irnd( nclusters ); } while (i == j); -+ do { -+ j = UseRandom::irnd( nclusters ); -+ } while (i == j); - - // check if one of the two potential clusters would be a colour octet state - octet = _isColour8( q[i], aq[j] ) || _isColour8( q[j], aq[i] ) ; -@@ -666,8 +1371,7 @@ - if(p->hasColour() && q->hasAntiColour()) { - cline = p-> colourLine(); - aline = q->antiColourLine(); -- } -- else { -+ } else { - cline = q-> colourLine(); - aline = p->antiColourLine(); - } -@@ -703,20 +1407,57 @@ - - - void ColourReconnector::persistentOutput(PersistentOStream & os) const { -- os << _clreco << _preco << _precoBaryonic << _algorithm << _initTemp << _annealingFactor -- << _annealingSteps << _triesPerStepFactor << ounit(_maxDistance,femtometer) -- << _octetOption; -+ os -+ << _clreco -+ << _algorithm -+ << _annealingFactor -+ << _annealingSteps -+ << _triesPerStepFactor -+ << _initTemp -+ << _preco -+ << _precoBaryonic -+ << _preco3M_3M -+ << _preco3M_BBbar -+ << _precoBbarB_3M -+ << _preco2B_2B -+ << _precoMB_MB -+ << _stepFactor -+ << _mesonToBaryonFactor -+ << ounit(_maxDistance, femtometer) -+ << _octetOption -+ << _cr2BeamClusters -+ << _debug -+ << _junctionMBCR -+ ; - } - - void ColourReconnector::persistentInput(PersistentIStream & is, int) { -- is >> _clreco >> _preco >> _precoBaryonic >> _algorithm >> _initTemp >> _annealingFactor -- >> _annealingSteps >> _triesPerStepFactor >> iunit(_maxDistance,femtometer) -- >> _octetOption; -+ is -+ >> _clreco -+ >> _algorithm -+ >> _annealingFactor -+ >> _annealingSteps -+ >> _triesPerStepFactor -+ >> _initTemp -+ >> _preco -+ >> _precoBaryonic -+ >> _preco3M_3M -+ >> _preco3M_BBbar -+ >> _precoBbarB_3M -+ >> _preco2B_2B -+ >> _precoMB_MB -+ >> _stepFactor -+ >> _mesonToBaryonFactor -+ >> iunit(_maxDistance, femtometer) -+ >> _octetOption -+ >> _cr2BeamClusters -+ >> _debug -+ >> _junctionMBCR -+ ; - } - - - void ColourReconnector::Init() { -- - static ClassDocumentation documentation - ("This class is responsible of the colour reconnection."); - -@@ -736,8 +1477,36 @@ - "Colour reconnections on", - 1); - -+ // Algorithm interface -+ static Switch interfaceAlgorithm -+ ("Algorithm", -+ "Specifies the colour reconnection algorithm", -+ &ColourReconnector::_algorithm, 0, true, false); -+ static SwitchOption interfaceAlgorithmPlain -+ (interfaceAlgorithm, -+ "Plain", -+ "Plain colour reconnection as in Herwig 2.5.0", -+ 0); -+ static SwitchOption interfaceAlgorithmStatistical -+ (interfaceAlgorithm, -+ "Statistical", -+ "Statistical colour reconnection using simulated annealing", -+ 1); -+ static SwitchOption interfaceAlgorithmBaryonic -+ (interfaceAlgorithm, -+ "Baryonic", -+ "Baryonic cluster reconnection", -+ 2); -+ static SwitchOption interfaceAlgorithmBaryonicMesonic -+ (interfaceAlgorithm, -+ "BaryonicMesonic", -+ "Baryonic cluster reconnection with reconnections to and from Mesonic Clusters", -+ 3); - -- static Parameter interfaceMtrpAnnealingFactor -+ -+ -+ // Statistical CR Parameters: -+ static Parameter interfaceMtrpAnnealingFactor - ("AnnealingFactor", - "The annealing factor is the ratio of the temperatures in two successive " - "temperature steps.", -@@ -766,50 +1535,75 @@ - false, false, Interface::limited); - - -- static Parameter interfaceRecoProb -+ -+ -+ // Plain and Baryonic CR Paramters -+ static Parameter interfaceRecoProb - ("ReconnectionProbability", -- "Probability that a found reconnection possibility is actually accepted", -+ "Probability that a found two meson to two meson reconnection possibility is actually accepted (used in Plain & Baryonic)", - &ColourReconnector::_preco, 0.5, 0.0, 1.0, - false, false, Interface::limited); - - static Parameter interfaceRecoProbBaryonic - ("ReconnectionProbabilityBaryonic", -- "Probability that a found reconnection possibility is actually accepted", -+ "Probability that a found reconnection possibility is actually accepted (used in Baryonic)", - &ColourReconnector::_precoBaryonic, 0.5, 0.0, 1.0, - false, false, Interface::limited); - - -- static Switch interfaceAlgorithm -- ("Algorithm", -- "Specifies the colour reconnection algorithm", -- &ColourReconnector::_algorithm, 0, true, false); -- static SwitchOption interfaceAlgorithmPlain -- (interfaceAlgorithm, -- "Plain", -- "Plain colour reconnection as in Herwig 2.5.0", -- 0); -- static SwitchOption interfaceAlgorithmStatistical -- (interfaceAlgorithm, -- "Statistical", -- "Statistical colour reconnection using simulated annealing", -- 1); -- static SwitchOption interfaceAlgorithmBaryonic -- (interfaceAlgorithm, -- "Baryonic", -- "Baryonic cluster reconnection", -- 2); -+ // BaryonicMesonic CR Paramters -+ static Parameter interfaceReconnectionProbability3Mto3M -+ ("ReconnectionProbability3Mto3M", -+ "Probability that a reconnection candidate is accepted for reconnecting 3M -> 3M\'", -+ &ColourReconnector::_preco3M_3M, 0.5, 0.0, 1.0, -+ false, false, Interface::limited); -+ static Parameter interfaceReconnectionProbability3MtoBBbar -+ ("ReconnectionProbability3MtoBBbar", -+ "Probability that a reconnection candidate is accepted for reconnecting 3M -> B,Bbar", -+ &ColourReconnector::_preco3M_BBbar, 0.5, 0.0, 1.0, -+ false, false, Interface::limited); -+ static Parameter interfaceReconnectionProbabilityBbarBto3M -+ ("ReconnectionProbabilityBbarBto3M", -+ "Probability that a reconnection candidate is accepted for reconnecting B,Bbar -> 3M", -+ &ColourReconnector::_precoBbarB_3M, 0.5, 0.0, 1.0, -+ false, false, Interface::limited); -+ static Parameter interfaceReconnectionProbability2Bto2B -+ ("ReconnectionProbability2Bto2B", -+ "Probability that a reconnection candidate is accepted for reconnecting 2B -> 2B\' or 2Bbar -> 2Bbar\'", -+ &ColourReconnector::_preco2B_2B, 0.5, 0.0, 1.0, -+ false, false, Interface::limited); -+ static Parameter interfaceReconnectionProbabilityMBtoMB -+ ("ReconnectionProbabilityMBtoMB", -+ "Probability that a reconnection candidate is accepted for reconnecting M,B -> M\',B\' or M,Bbar -> M\',Bbar\'", -+ &ColourReconnector::_precoMB_MB, 0.5, 0.0, 1.0, -+ false, false, Interface::limited); - -- static Parameter interfaceMaxDistance -+ static Parameter interfaceFactorforStep -+ ("StepFactor", -+ "Factor for how many reconnection-tries are made in the BaryonicMesonic algorithm", -+ &ColourReconnector::_stepFactor, 1.0, 0.11111, 10., -+ false, false, Interface::limited);// at least 3 Clusters -> _stepFactorMin=1/9 -+ -+ static Parameter interfaceMesonToBaryonFactor -+ ("MesonToBaryonFactor", -+ "Factor for comparing mesonic clusters to baryonic clusters in the displacement if BaryonicMesonic CR model is chosen", -+ &ColourReconnector::_mesonToBaryonFactor, 2.0, 0.5, 3.0, -+ false, false, Interface::limited); -+ -+ -+ -+ // General Parameters and switches -+ static Parameter interfaceMaxDistance - ("MaxDistance", - "Maximum distance between the clusters at which to consider rearrangement" -- " to avoid colour reconneections of displaced vertices", -+ " to avoid colour reconneections of displaced vertices (used in all Algorithms). No unit means femtometer", - &ColourReconnector::_maxDistance, femtometer, 1000.*femtometer, 0.0*femtometer, 1e100*femtometer, - false, false, Interface::limited); - - -- static Switch interfaceOctetTreatment -+ static Switch interfaceOctetTreatment - ("OctetTreatment", -- "Which octets are not allowed to be reconnected", -+ "Which octets are not allowed to be reconnected (used in all Algorithms)", - &ColourReconnector::_octetOption, 0, false, false); - static SwitchOption interfaceOctetTreatmentFinal - (interfaceOctetTreatment, -@@ -821,6 +1615,51 @@ - "All", - "Prevent for all octets", - 1); -+ static Switch interfaceCR2BeamClusters -+ ("CR2BeamClusters", -+ "Option for colour reconnecting 2 beam remnant clusters if the number of clusters is 2.", -+ &ColourReconnector::_cr2BeamClusters, 0, true, false); -+ static SwitchOption interfaceCR2BeamClustersYes -+ (interfaceCR2BeamClusters, -+ "Yes", -+ "If possible CR 2 beam clusters", -+ 1); -+ static SwitchOption interfaceCR2BeamClustersNo -+ (interfaceCR2BeamClusters, -+ "No", -+ "If possible do not CR 2 beam clusters", -+ 0); -+ static Switch interfaceJunction -+ ("Junction", -+ "Option for using Junction-like displacement in rapidity-phi plane to compare baryonic cluster " -+ "instead of pairwise distance (for BaryonicMesonic model)", -+ &ColourReconnector::_junctionMBCR, 1, true, false); -+ static SwitchOption interfaceJunctionYes -+ (interfaceJunction, -+ "Yes", -+ "Using junction-like model instead of pairwise distance model", -+ 1); -+ static SwitchOption interfaceJunctionNo -+ (interfaceJunction, -+ "No", -+ "Using pairwise distance model instead of junction-like model", -+ 0); -+ -+ // Debug -+ static Switch interfaceDebug -+ ("Debug", -+ "Make a file with some Information of the BaryonicMesonic Algorithm", -+ &ColourReconnector::_debug, 0, true, false); -+ static SwitchOption interfaceDebugNo -+ (interfaceDebug, -+ "No", -+ "Debug Information for ColourReconnector Off", -+ 0); -+ static SwitchOption interfaceDebugYes -+ (interfaceDebug, -+ "Yes", -+ "Debug Information for ColourReconnector On", -+ 1); - - } - -diff -r 1781af1e3113 -r 5d17715d8d26 Hadronization/ColourReconnector.h ---- a/Hadronization/ColourReconnector.h Mon Jul 18 09:45:26 2022 +0100 -+++ b/Hadronization/ColourReconnector.h Wed Jul 20 18:37:09 2022 +0200 -@@ -59,6 +59,43 @@ - Energy2 _clusterMassSum(const PVector & q, const PVector & aq) const; - - -+ -+ /** -+ * @brief calculates the "euclidean" distance of two quarks in the -+ * rapdidity-phi plane -+ * @arguments p, q the two quarks -+ * @return the dimensionless distance: -+ * \deltaR_12=sqrt(\deltaY_12^2 + \delta\phi_12^2) -+ */ -+ double _displacement(tcPPtr p, tcPPtr q) const; -+ -+ -+ /** -+ * @brief calculates the "euclidean" distance of a the 3 (anti)quarks -+ * (anti)baryonic cluster in the rapdidity-phi plane -+ * @arguments p, q the two quarks -+ * @return the dimensionless distance: -+ * if Junction is enabled the difference of all 3 quarks -+ * with respect to their mean point is calculated: -+ * = sum_i Y_i/3 -+ * <\phi> = sum_i \phi_i/3 -+ * \deltaR = sum_i sqrt( (Y_i - )^2 + (\phi_i - )^2) -+ * -+ * if Junction is disabled the difference of all distinct -+ * pairing of the 3 quarks is computed: -+ * \deltaR_ij = sqrt(\deltaY_ij^2 + \delta\phi_ij^2) -+ * \deltaR_tot = sum_{i,j &cls, -+ const unsigned &baryonic, -+ string &kind_of_reco, -+ int (&reco_info)[3]) const; - - /** - * @brief Reconnects the constituents of the given clusters to the (only) - * other possible cluster combination. - * @return pair of pointers to the two new clusters -+ * Used for Plain and Baryonic Colour Reconnection models - */ - pair _reconnect(ClusterPtr &c1, ClusterPtr &c2) const; - - /** -- * Reconnection method for baryonic reconenction model -+ * @brief Reconnects (2B->2B) the constituents of the given clusters to -+ another possible cluster combination whose information is given -+ in s1 and s2. -+ * @arguments c1 and c2 are baryonic clusters and s1 and s2 are the respective -+ indices of the constituents of c1 and c2 respectively -+ * @return pair of pointers to the two new clusters -+ * Used only in BaryonicMesonic algorithm and will exchange constituent s1 of -+ * c1 with constituent s2 of c2 -+ */ -+ pair _reconnect2Bto2B(ClusterPtr &c1, ClusterPtr &c2, const int s1, const int s2) const; -+ -+ /** -+ * @brief Reconnects (B,Bbar->3M) the constituents of the given clusters to -+ another possible cluster combination whose information is given in -+ s0, s1 and s2. -+ * @arguments c1 and c2 are baryonic clusters. s0, s1 and s2 are the respective -+ indices which determine the CR -+ * @return tuple of pointers to the 3 new mesonic clusters -+ * Used only in BaryonicMesonic algorithm and will form 3 mesonic clusters according -+ * to the indices s0, s1 and s2. The i-th constituent of c1 is connected to the si-th -+ * constituent of c2 -+ */ -+ std::tuple _reconnectBbarBto3M(ClusterPtr &c1, ClusterPtr &c2, const int s0, const int s1, const int s2 ) const; -+ -+ /** -+ * @brief Reconnects (3M->3M) the constituents of the given clusters to -+ another possible cluster combination whose information is given in -+ sinfos -+ * @arguments c1, c2 and c3 are mesonic clusters. infos[3] are the respective -+ indices which determine the CR -+ * @return tuple of pointers to the 3 CR'd mesonic clusters -+ * Used only in BaryonicMesonic algorithm and will reconnect 3 mesonic clusters according -+ * to the infos, which determine the CR. The coloured quark of the i-th cluster forms -+ * a new cluster with the anticoloured quark of the info[i]-th cluster -+ */ -+ std::tuple _reconnect3Mto3M(ClusterPtr &c1, ClusterPtr &c2, ClusterPtr &c3, const int infos[3] ) const; -+ -+ -+ /** -+ * Reconnection method for a Baryonic and a Mesonic Cluster to a Baryonic and a Mesonic Cluster -+ * s0 is the Number of the (Anti)Patrticle of the Baryonic Cluster , which should be swapped with the Anti(Particle) of the Mesonic Cluster -+ */ -+ /** -+ * @brief Reconnects the constituents of the given clusters to -+ another possible cluster combination whose information -+ is given in s0. -+ * @arguments c1 and c2 are one baryonic and one mesonic cluster respectively -+ and s0 is the respective index of the constituents of the baryonic -+ cluster which is to be exchangeed with the mesonic cluster. -+ * @return pair of pointers to the two new clusters -+ * Used only in BaryonicMesonic algorithm and will exchange constituent s0 of -+ * the baryonic cluster with the (anti)-quark of the mesonic cluster - */ -- pair _reconnectBaryonic(ClusterPtr &c1, ClusterPtr &c2) const; -+ pair _reconnectMBtoMB(ClusterPtr &c1, ClusterPtr &c2, const int s0) const; -+ -+ -+ -+ -+ /** -+ * Methods for the BaryonicMesonic CR algorithm -+ * Find the best reconnection option for the respective cluster-combination -+ * -+ */ -+ -+ /** -+ * @brief veto for far apart clusters -+ * @arguments expects at most 3 CluVecIt in clu vector -+ * @return returns true if clusters are more distant than _maxDistance -+ * in space -+ * TODO: problematic maybe add option to turn off -+ */ -+ bool _clustersFarApart( const std::vector & clu ) const; -+ -+ /** -+ * @brief Does reconnect 2 beam clusters for BaryonicMesonic CR model -+ * if option CR2BeamClusters is enabled -+ * @arguments cv cluster vector -+ */ -+ void _doReco2BeamClusters(ClusterVector & cv) const; -+ - -+ /** -+ * @brief finds the best reconnection option and stores it in bswap1 -+ * and bswap2 (2B->2B colour reconnection) -+ * @arguments c1 and c2 cluster pointers and kind_of_reco will output -+ * the best found reconnection option for c1 and c2 -+ */ -+ void _2Bto2BreconnectionFinder(ClusterPtr &c1, ClusterPtr &c2, int &bswap1, int &bswap2, double mindisplsum, string &kind_of_reco ) const; -+ -+ /** -+ * @brief finds the best reconnection option and stores it in mswap0 -+ * mswap1 and mswap2 (BbarB->3M colour reconnection) -+ * @arguments c1 and c2 cluster pointers and kind_of_reco will output -+ * the best found reconnection option for c1 and c2 -+ */ -+ void _BbarBto3MreconnectionFinder(ClusterPtr &c1, ClusterPtr &c2, int &mswap0, int &mswap1, int &mswap2, -+ double min_displ_sum, string & kind_of_reco) const; -+ -+ /** -+ * @brief finds the best reconnection option and stores it in swap -+ * (BM->BM colour reconnection) -+ * @arguments c1 and c2 cluster pointers and kind_of_reco will output -+ * the best found reconnection option for c1 and c2 -+ */ -+ void _BMtoBMreconnectionfinder(ClusterPtr &c1, ClusterPtr &c2, int &swap, -+ double min_displ_sum, string &kind_of_reco) const; -+ -+ /** -+ * @brief finds the best reconnection option and stores it in swap0, -+ * swap1 and swap2 (3M->{3M,BbarB} colour reconnection) -+ * @arguments c1 and c2 cluster pointers and kind_of_reco will output -+ * the best found reconnection option for c1 and c2 -+ */ -+ void _3MtoXreconnectionfinder(std::vector &cv, int &swap0, int &swap1, -+ int &swap2, double min_displ_sum, string &kind_of_reco) const; - - /** - * @brief At random, swap two antiquarks, if not excluded by the -@@ -156,6 +329,25 @@ - */ - int _clreco = 0; - -+ /** -+ * Do we want debug informations? -+ */ -+ int _debug = 0; -+ -+ /** -+ * @brief Junction-like model for BaryonicMesonic model -+ * Do we want to use the junction-like model for -+ * computing the displacements of BaryonicMesonic model? -+ * otherwise pairwise distances are used. -+ * If _junctionMBCR is activated the displacements are computed in the -+ * rapidity-Phi plane by difference to the average rapidity and phi: -+ * DeltaR_i^2 = (rapidity_i - meanRap)^2 + (phi_i - meanPhi)^2 -+ * DeltaR = sum_i DeltaR_i -+ * if _junctionMBCR=0 the displacements are computed: -+ * DeltaR_ij^2 = (rapidity_i - rapidity_j)^2 + (phi_i - phi_j)^2 -+ * DeltaR = sum_i,j 3M' -+ * used in BaryonicMesonic -+ * NOTE: if 0 this type of reconnection is not even tried -+ */ -+ double _preco3M_3M = 0.5; -+ -+ /** -+ * Probability that a found reconnection possibility is actually accepted. -+ * For reconnecting 3M -> B,Bbar -+ * used in BaryonicMesonic -+ */ -+ double _preco3M_BBbar = 0.5; -+ -+ /** -+ * Probability that a found reconnection possibility is actually accepted. -+ * For reconnecting Bbar,B -> 3M -+ * used in BaryonicMesonic - */ -+ double _precoBbarB_3M = 0.5; -+ -+ /** -+ * Probability that a found reconnection possibility is actually accepted. -+ * For reconnecting 2B -> 2B' or 2Bbar -> 2Bbar' -+ * used in BaryonicMesonic -+ * NOTE: if 0 this type of reconnection is not even tried -+ */ -+ double _preco2B_2B = 0.5; -+ -+ /** -+ * Probability that a found reconnection possibility is actually accepted. -+ * For reconnecting M,B -> M',B' or M,Bbar -> M',Bbar' -+ * used in BaryonicMesonic -+ * NOTE: if 0 this type of reconnection is not even tried -+ */ -+ double _precoMB_MB = 0.5; -+ -+ /** -+ * For the BaryonicMesonic algorithm -+ * How many times do suggest cluster for reconnection? -+ * n(reconnectionstries) = _stepFactor * n(clusters)*n(clusters); -+ */ -+ double _stepFactor = 5.0; -+ -+ /** -+ * Factor for comparing mesonic clusters to baryonic clusters -+ */ -+ double _mesonToBaryonFactor = 2.0; - - /** - * Maximium distance for reconnections -+ * TODO: remove if issues with anticausality are discussed and resolved - */ -- Length _maxDistance = picometer; -+ Length _maxDistance = femtometer; - - /** - * @return true, if the two partons are splitting products of the same -@@ -213,6 +456,11 @@ - */ - unsigned int _octetOption = 0; - -+ /** -+ * Option for colour reconnecting 2 Beam Clusters if no others are present -+ */ -+ int _cr2BeamClusters = 0; -+ - public: - - /** @name Functions used by the persistent I/O system. */ -diff -r 1781af1e3113 -r 5d17715d8d26 src/defaults/Hadronization.in ---- a/src/defaults/Hadronization.in Mon Jul 18 09:45:26 2022 +0100 -+++ b/src/defaults/Hadronization.in Wed Jul 20 18:37:09 2022 +0200 -@@ -41,11 +41,37 @@ - newdef LightClusterDecayer:HadronSelector HadronSelector - newdef ClusterDecayer:HadronSelector HadronSelector - -+# ColourReconnector Default Parameters - newdef ColourReconnector:ColourReconnection Yes --newdef ColourReconnector:ReconnectionProbabilityBaryonic 0.7 -+newdef ColourReconnector:Algorithm Baryonic -+ -+# Statistical CR Parameters: -+newdef ColourReconnector:AnnealingFactor 0.9 -+newdef ColourReconnector:AnnealingSteps 50 -+newdef ColourReconnector:TriesPerStepFactor 5.0 -+newdef ColourReconnector:InitialTemperature 0.1 -+ -+# Plain and Baryonic CR Paramters - newdef ColourReconnector:ReconnectionProbability 0.95 --newdef ColourReconnector:Algorithm Baryonic -+newdef ColourReconnector:ReconnectionProbabilityBaryonic 0.7 -+ -+# BaryonicMesonic and BaryonicMesonic CR Paramters -+newdef ColourReconnector:ReconnectionProbability3Mto3M 0.5 -+newdef ColourReconnector:ReconnectionProbability3MtoBBbar 0.5 -+newdef ColourReconnector:ReconnectionProbabilityBbarBto3M 0.5 -+newdef ColourReconnector:ReconnectionProbability2Bto2B 0.05 -+newdef ColourReconnector:ReconnectionProbabilityMBtoMB 0.5 -+newdef ColourReconnector:StepFactor 1.0 -+newdef ColourReconnector:MesonToBaryonFactor 1.333 -+ -+# General Parameters and switches -+newdef ColourReconnector:MaxDistance 1.0e50 - newdef ColourReconnector:OctetTreatment All -+newdef ColourReconnector:CR2BeamClusters No -+newdef ColourReconnector:Junction Yes -+# Debugging -+newdef ColourReconnector:Debug No -+ - - # Clustering parameters for light quarks - newdef ClusterFissioner:ClMaxLight 3.649