diff --git a/Hadronization/ClusterFissioner.cc b/Hadronization/ClusterFissioner.cc --- a/Hadronization/ClusterFissioner.cc +++ b/Hadronization/ClusterFissioner.cc @@ -1,1121 +1,1103 @@ // -*- C++ -*- // // ClusterFissioner.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. // // // Thisk is the implementation of the non-inlined, non-templated member // functions of the ClusterFissioner class. // #include "ClusterFissioner.h" #include #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 */ +