diff --git a/Hadronization/ClusterFissioner.cc b/Hadronization/ClusterFissioner.cc --- a/Hadronization/ClusterFissioner.cc +++ b/Hadronization/ClusterFissioner.cc @@ -1,1112 +1,1121 @@ // -*- C++ -*- // // ClusterFissioner.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // 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",""); 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), _fissionCluster(0), _fissionPwtUquark(1), _fissionPwtDquark(1), _fissionPwtSquark(0.5), _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; do { succeeded = false; ++counter; if (_enhanceSProb == 0){ drawNewFlavour(newPtr1,newPtr2); } else { Energy2 mass2 = clustermass(cluster); drawNewFlavourEnhanced(newPtr1,newPtr2,mass2); } // 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; if (CheckId::isExotic(ptrQ2->dataPtr())) exp2 = _pSplitExotic; else if(CheckId::hasBottom(ptrQ2->dataPtr())) exp2 = _pSplitBottom; else if(CheckId::hasCharm(ptrQ2->dataPtr())) exp2 = _pSplitCharm; // 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); 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. ****************************/ toHadron1 = _hadronsSelector->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc1); if(toHadron1) Mc1 = toHadron1->mass(); toHadron2 = _hadronsSelector->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc2); if(toHadron2) Mc2 = toHadron2->mass(); // 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(); } else if(!toHadron2) { toHadron2 = _hadronsSelector->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc-Mc1); if(toHadron2) Mc2 = toHadron2->mass(); } } 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 /****************** * 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; do { succeeded = false; ++counter; if (_enhanceSProb == 0) { drawNewFlavour(newPtr1,newPtr2); } else { Energy2 mass2 = clustermass(cluster); drawNewFlavourEnhanced(newPtr1,newPtr2, mass2); } // 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; // 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); 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(); // 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(); 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(); } else if(!toDiQuark) { Mc2 = _hadronsSelector->massLightestHadron(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr()); 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 = (20. < scale || scale < 0.) ? 0. : pow(_hadronsSelector->pwtSquark(),scale); + prob_s = (_maxScale < scale) ? 0. : pow(_hadronsSelector->pwtSquark(),scale); else if (_enhanceSProb == 2) - prob_s = (20. < scale || scale < 0.) ? 0. : exp(-scale); + 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 = (20. < scale || scale < 0.) ? 0. : pow(_fissionPwtSquark,scale); + prob_s = (_maxScale < scale) ? 0. : pow(_fissionPwtSquark,scale); else if (_enhanceSProb == 2) - prob_s = (20. < scale || scale < 0.) ? 0. : exp(-scale); + 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){ 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,412 +1,439 @@ // -*- C++ -*- // // ClusterFissioner.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_ClusterFissioner_H #define HERWIG_ClusterFissioner_H #include #include "CluHadConfig.h" #include "HadronSelector.h" #include "ClusterFissioner.fh" 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: /** * 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; + /** + * 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; - Energy2 clustermass(const ClusterPtr & cluster); - + /** + * 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; -protected: - void drawNewFlavourEnhanced(PPtr& newPtrPos,PPtr& newPtrNeg, Energy2 mass2) const; + /** + * 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/HwppSelector.cc b/Hadronization/HwppSelector.cc --- a/Hadronization/HwppSelector.cc +++ b/Hadronization/HwppSelector.cc @@ -1,248 +1,253 @@ // -*- C++ -*- // // HwppSelector.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the HwppSelector class. // #include "HwppSelector.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Utilities/Kinematics.h" #include "ThePEG/Utilities/Selector.h" #include "ThePEG/Repository/UseRandom.h" #include "CheckId.h" #include #include using namespace Herwig; DescribeClass describeHwppSelector("Herwig::HwppSelector",""); IBPtr HwppSelector::clone() const { return new_ptr(*this); } IBPtr HwppSelector::fullclone() const { return new_ptr(*this); } void HwppSelector::doinit() { HadronSelector::doinit(); } void HwppSelector::persistentOutput(PersistentOStream & os) const { os << _mode << _enhanceSProb << ounit(_m0Decay,GeV) << _massMeasure; } void HwppSelector::persistentInput(PersistentIStream & is, int) { is >> _mode >> _enhanceSProb >> iunit(_m0Decay,GeV) >> _massMeasure; } void HwppSelector::Init() { static ClassDocumentation documentation ("The HwppSelector class implements the Herwig algorithm for selecting" " the hadrons", "The hadronization used the selection algorithm described in \\cite{Kupco:1998fx}.", "%\\cite{Kupco:1998fx}\n" "\\bibitem{Kupco:1998fx}\n" " A.~Kupco,\n" " ``Cluster hadronization in HERWIG 5.9,''\n" " arXiv:hep-ph/9906412.\n" " %%CITATION = HEP-PH/9906412;%%\n" ); // put useMe() only in correct place! static Switch interfaceMode ("Mode", "Which algorithm to use", &HwppSelector::_mode, 1, false, false); static SwitchOption interfaceModeKupco (interfaceMode, "Kupco", "Use the Kupco approach", 0); static SwitchOption interfaceModeHwpp (interfaceMode, "Hwpp", "Use the Herwig approach", 1); static Switch interfaceEnhanceSProb ("EnhanceSProb", "Option for enhancing strangeness", &HwppSelector::_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", &HwppSelector::_massMeasure,0,false,false); static SwitchOption interfaceMassMeasureMass (interfaceMassMeasure, "Mass", "Mass Measure", 0); static SwitchOption interfaceMassMeasureLambda (interfaceMassMeasure, "Lambda", "Lambda Measure", 1); static Parameter interfaceDecayMassScale ("DecayMassScale", "Cluster decay mass scale", &HwppSelector::_m0Decay, GeV, 1.0*GeV, 0.1*GeV, 50.*GeV, false, false, Interface::limited); } pair HwppSelector::chooseHadronPair(const Energy cluMass,tcPDPtr par1, tcPDPtr par2,tcPDPtr ) const { // if either of the input partons is a diquark don't allow diquarks to be // produced bool diquark = !(DiquarkMatcher::Check(par1->id()) || DiquarkMatcher::Check(par2->id())); bool quark = true; // if the Herwig algorithm if(_mode ==1) { if(UseRandom::rnd() > 1./(1.+pwtDIquark()) &&cluMass > massLightestBaryonPair(par1,par2)) { diquark = true; quark = false; } else { useMe(); diquark = false; quark = true; } } // weights for the different possibilities Energy weight, wgtsum(ZERO); // loop over all hadron pairs with the allowed flavours static vector hadrons; hadrons.clear(); for(unsigned int ix=0;ixiColour())) == 3 && !DiquarkMatcher::Check(quarktopick->id())) continue; if(!diquark && abs(int(quarktopick->iColour())) == 3 && DiquarkMatcher::Check(quarktopick->id())) continue; HadronTable::const_iterator tit1 = table().find(make_pair(abs(par1->id()),quarktopick->id())); HadronTable::const_iterator tit2 = table().find(make_pair(quarktopick->id(),abs(par2->id()))); // If not in table skip if(tit1 == table().end()||tit2==table().end()) continue; // tables empty skip const KupcoData & T1 = tit1->second; const KupcoData & T2 = tit2->second; if(T1.empty()||T2.empty()) continue; // if too massive skip if(cluMass <= T1.begin()->mass + T2.begin()->mass) continue; // loop over the hadrons KupcoData::const_iterator H1,H2; for(H1 = T1.begin();H1 != T1.end(); ++H1) { for(H2 = T2.begin();H2 != T2.end(); ++H2) { // break if cluster too light if(cluMass < H1->mass + H2->mass) break; // calculate the weight double pwtstrange; if (quarktopick->id() == 3){ + // Strangeness weight takes the automatic flat weight pwtstrange = pwt(3); + // Scaling strangeness enhancement if (_enhanceSProb == 1){ double scale = double(sqr(_m0Decay/cluMass)); - pwtstrange = (20. < scale || scale < 0.) ? 0. : pow(pwtstrange,scale); + pwtstrange = (_maxScale < scale) ? 0. : pow(pwtstrange,scale); } + // Exponential strangeness enhancement else if (_enhanceSProb == 2){ Energy2 mass2; Energy endpointmass = par1->mass() + par2->mass(); + // Choose to use either the cluster mass + // or to use the lambda measure mass2 = (_massMeasure == 0) ? sqr(cluMass) : sqr(cluMass) - sqr(endpointmass); double scale = double(sqr(_m0Decay)/mass2); - pwtstrange = (20. < scale || scale < 0.) ? 0. : exp(-scale); + pwtstrange = (_maxScale < scale) ? 0. : exp(-scale); } weight = pwtstrange * H1->overallWeight * H2->overallWeight * Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass ); } else { weight = pwt(quarktopick->id()) * H1->overallWeight * H2->overallWeight * Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass ); } int signQ = 0; assert (par1 && quarktopick); assert (par2); assert(quarktopick->CC()); if(CheckId::canBeHadron(par1, quarktopick->CC()) && CheckId::canBeHadron(quarktopick, par2)) signQ = +1; else if(CheckId::canBeHadron(par1, quarktopick) && CheckId::canBeHadron(quarktopick->CC(), par2)) signQ = -1; else { cerr << "Could not make sign for" << par1->id()<< " " << quarktopick->id() << " " << par2->id() << "\n"; assert(false); } if (signQ == -1) quarktopick = quarktopick->CC(); // construct the object with the info Kupco a(quarktopick, H1->ptrData, H2->ptrData, weight); hadrons.push_back(a); wgtsum += weight; } } } if (hadrons.empty()) return make_pair(tcPDPtr(),tcPDPtr()); // select the hadron wgtsum *= UseRandom::rnd(); unsigned int ix=0; do { wgtsum-= hadrons[ix].weight; ++ix; } while(wgtsum > ZERO && ix < hadrons.size()); if(ix == hadrons.size() && wgtsum > ZERO) return make_pair(tcPDPtr(),tcPDPtr()); --ix; assert(hadrons[ix].idQ); int signHad1 = signHadron(par1, hadrons[ix].idQ->CC(), hadrons[ix].hadron1); int signHad2 = signHadron(par2, hadrons[ix].idQ, hadrons[ix].hadron2); assert( signHad1 != 0 && signHad2 != 0 ); return make_pair ( signHad1 > 0 ? hadrons[ix].hadron1 : tcPDPtr(hadrons[ix].hadron1->CC()), signHad2 > 0 ? hadrons[ix].hadron2 : tcPDPtr(hadrons[ix].hadron2->CC())); } diff --git a/Hadronization/HwppSelector.h b/Hadronization/HwppSelector.h --- a/Hadronization/HwppSelector.h +++ b/Hadronization/HwppSelector.h @@ -1,151 +1,168 @@ // -*- C++ -*- // // HwppSelector.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_HwppSelector_H #define HERWIG_HwppSelector_H // // This is the declaration of the HwppSelector class. // #include "HadronSelector.h" #include "HwppSelector.fh" namespace Herwig { using namespace ThePEG; /** \ingroup hadronization * The HwppSelector class selects the hadrons produced in cluster decay using * the Herwig variant of the cluster model. * * @see \ref HwppSelectorInterfaces "The interfaces" * defined for HwppSelector. */ class HwppSelector: public HadronSelector { public: /** * The default constructor. */ HwppSelector() : HadronSelector(1), _mode(1), _enhanceSProb(0), _m0Decay(1.*GeV) {} /** * * This method is used to choose a pair of hadrons. * * Given the mass of a cluster and the particle pointers of its * two (or three) constituents, this returns the pair of particle pointers of * the two hadrons with proper flavour numbers. * Furthermore, the first of the two hadron must have the * constituent with par1, and the second must have the constituent with par2. * At the moment it does *nothing* in the case that also par3 is present. * * Kupco's method is used, rather than one used in FORTRAN HERWIG * The idea is to build on the fly a table of all possible pairs * of hadrons (Had1,Had2) (that we can call "cluster decay channels") * which are kinematically above threshold and have flavour * Had1=(par1,quarktopick->CC()), Had2=(quarktopick,par2), where quarktopick * is the poniter of: * --- d, u, s, c, b * if either par1 or par2 is a diquark; * --- d, u, s, c, b, dd, ud, uu, sd, su, ss, * cd, cu, cs, cc, bd, bu, bs, bc, bb * if both par1 and par2 are quarks. * The weight associated with each channel is given by the product * of: the phase space available including the spin factor 2*J+1, * the constant weight factor for chosen idQ, * the octet-singlet isoscalar mixing factor, and finally * the singlet-decuplet weight factor. */ pair chooseHadronPair(const Energy cluMass,tcPDPtr par1, tcPDPtr par2,tcPDPtr par3 = PDPtr()) const ; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @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(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ HwppSelector & operator=(const HwppSelector &) = delete; private: /** * Which algorithm to use */ unsigned int _mode; + /** + * 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 _m0Decay; + /** + * 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 that stops the scale in strangeness enhancement from + * becoming too large + */ + const double _maxScale = 20.; + }; } #endif /* HERWIG_HwppSelector_H */ diff --git a/Hadronization/PartonSplitter.cc b/Hadronization/PartonSplitter.cc --- a/Hadronization/PartonSplitter.cc +++ b/Hadronization/PartonSplitter.cc @@ -1,555 +1,555 @@ // -*- C++ -*- // // PartonSplitter.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the PartonSplitter class. // #include "PartonSplitter.h" #include #include #include #include #include #include #include #include "ThePEG/Interface/Parameter.h" #include #include #include "ThePEG/Repository/UseRandom.h" #include "Herwig/Utilities/Kinematics.h" #include #include "ClusterHadronizationHandler.h" #include #include #include "CheckId.h" using namespace Herwig; IBPtr PartonSplitter::clone() const { return new_ptr(*this); } IBPtr PartonSplitter::fullclone() const { return new_ptr(*this); } void PartonSplitter::persistentOutput(PersistentOStream & os) const { os << _quarkSelector << ounit(_gluonDistance,femtometer) << _splitGluon << _splitPwtUquark << _splitPwtDquark << _splitPwtSquark << _enhanceSProb << ounit(_m0,GeV) << _massMeasure; } void PartonSplitter::persistentInput(PersistentIStream & is, int) { is >> _quarkSelector >> iunit(_gluonDistance,femtometer) >> _splitGluon >> _splitPwtUquark >> _splitPwtDquark >> _splitPwtSquark >>_enhanceSProb >> iunit(_m0,GeV) >> _massMeasure; } DescribeClass describePartonSplitter("Herwig::PartonSplitter",""); void PartonSplitter::Init() { static ClassDocumentation documentation ("This class is reponsible of the nonperturbative splitting of partons"); static Switch interfaceSplit ("Split", "Option for different splitting options", &PartonSplitter::_splitGluon, 0, false, false); static SwitchOption interfaceSplitDefault (interfaceSplit, "ud", "Normal cluster splitting where only u and d quarks are drawn is used.", 0); static SwitchOption interfaceSplitAll (interfaceSplit, "uds", "Alternative cluster splitting where all light quark pairs (u, d, s) can be drawn.", 1); static Parameter interfaceSplitPwtUquark ("SplitPwtUquark", "Weight for splitting in U quarks", &PartonSplitter::_splitPwtUquark, 1, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceSplitPwtDquark ("SplitPwtDquark", "Weight for splitting in D quarks", &PartonSplitter::_splitPwtDquark, 1, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceSplitPwtSquark ("SplitPwtSquark", "Weight for splitting in S quarks", &PartonSplitter::_splitPwtSquark, 0.5, 0.0, 1.0, false, false, Interface::limited); static Switch interfaceEnhanceSProb ("EnhanceSProb", "Option to enhance the strangeness weight (MassSplit Switch needs to be on)", &PartonSplitter::_enhanceSProb,0,false,false); static SwitchOption interfaceEnhanceSProbNo (interfaceEnhanceSProb, "No", "No scaling for strangeness", 0); static SwitchOption interfaceEnhanceSProbScale (interfaceEnhanceSProb, "Scaled", "Power-law scaling for strangeness", 1); static SwitchOption interfaceEnhanceSProbExp (interfaceEnhanceSProb, "Exponential", "Exponential suppression for strangeness", 2); static Switch interfaceMassMeasure ("MassMeasure", "Option to use different mass measures", &PartonSplitter::_massMeasure,0,false,false); static SwitchOption interfaceMassMeasureMass (interfaceMassMeasure, "Mass", "Mass Measure", 0); static SwitchOption interfaceMassMeasureLambda (interfaceMassMeasure, "Lambda", "Lambda Measure", 1); static Parameter interfaceMassScale ("MassScale", "Mass scale for g->qqb strangeness enhancement", &PartonSplitter::_m0, GeV, 20.*GeV, 1.*GeV, 1000.*GeV, false, false, Interface::limited); } void PartonSplitter::split(PVector & tagged) { // set the gluon c tau once and for all static bool first = true; if(first) { _gluonDistance = hbarc*getParticleData(ParticleID::g)->constituentMass()/ ClusterHadronizationHandler::currentHandler()->minVirtuality2(); first = false; } // Copy incoming for the (possible sorting and) splitting PVector particles = tagged; // Switch to enhance strangeness if (_enhanceSProb >= 1){ colourSorted(particles); } PVector newtag; Energy2 Q02 = 0.99*sqr(getParticleData(ParticleID::g)->constituentMass()); // Loop over all of the particles in the event. // Loop counter for colourSorted for(PVector::iterator pit = particles.begin(); pit!=particles.end(); ++pit) { // only considering gluons so add other particles to list of particles if( (**pit).data().id() != ParticleID::g ) { newtag.push_back(*pit); continue; } // should not have been called for massless or space-like gluons if((**pit).momentum().m2() <= 0.0*sqr(MeV) ) { throw Exception() << "Spacelike or massless gluon m2= " << (**pit).momentum().m2()/GeV2 << "GeV2 in PartonSplitter::split()" << Exception::eventerror; } // time like gluon gets split PPtr ptrQ = PPtr(); PPtr ptrQbar = PPtr(); if (_enhanceSProb == 0){ splitTimeLikeGluon(*pit,ptrQ,ptrQbar); } else { size_t i = pit - particles.begin(); massSplitTimeLikeGluon(*pit, ptrQ, ptrQbar, i); } ptrQ->scale(Q02); ptrQbar->scale(Q02); (*pit)->colourLine()->addColoured(ptrQ); (*pit)->addChild(ptrQ); newtag.push_back(ptrQ); (*pit)->antiColourLine()->addAntiColoured(ptrQbar); (*pit)->addChild(ptrQbar); newtag.push_back(ptrQbar); // set the life length of gluon Length distance = UseRandom::rndExp(_gluonDistance); (**pit).setLifeLength((distance/(**pit).mass())*(**pit).momentum()); // assume quarks same position as gluon ptrQ ->setVertex((**pit).decayVertex()); ptrQ ->setLifeLength(Lorentz5Distance()); ptrQbar->setVertex((**pit).decayVertex()); ptrQbar->setLifeLength(Lorentz5Distance()); } swap(tagged,newtag); } void PartonSplitter::splitTimeLikeGluon(tcPPtr ptrGluon, PPtr & ptrQ, PPtr & ptrQbar){ // select the quark flavour tPDPtr quark; long idNew=0; switch(_splitGluon){ case 0: quark = _quarkSelector.select(UseRandom::rnd()); break; case 1: if ( ptrGluon->momentum().m() < 2.0 *getParticle(ThePEG::ParticleID::s)->data().constituentMass() ) { throw Exception() << "Impossible Kinematics in PartonSplitter::splitTimeLikeGluon()" << Exception::runerror; } // Only allow light quarks u,d,s with the probabilities double prob_d = _splitPwtDquark; double prob_u = _splitPwtUquark; double prob_s = _splitPwtSquark; int choice = UseRandom::rnd3(prob_u, prob_d, prob_s); switch(choice) { case 0: idNew = ThePEG::ParticleID::u; break; case 1: idNew = ThePEG::ParticleID::d; break; case 2: idNew = ThePEG::ParticleID::s; break; } ptrQ = getParticle(idNew); ptrQbar = getParticle(-idNew); break; } // Solve the kinematics of the two body decay G --> Q + Qbar Lorentz5Momentum momentumQ; Lorentz5Momentum momentumQbar; double cosThetaStar = UseRandom::rnd( -1.0 , 1.0 ); using Constants::pi; double phiStar = UseRandom::rnd( -pi , pi ); Energy constituentQmass; if(_splitGluon==0) { constituentQmass = quark->constituentMass(); }else{ constituentQmass = ptrQ->data().constituentMass(); } if (ptrGluon->momentum().m() < 2.0*constituentQmass) { throw Exception() << "Impossible Kinematics in PartonSplitter::splitTimeLikeGluon()" << Exception::eventerror; } Kinematics::twoBodyDecay(ptrGluon->momentum(), constituentQmass, constituentQmass, cosThetaStar, phiStar, momentumQ, momentumQbar ); // Create quark and anti-quark particles of the chosen flavour // and set they 5-momentum (the mass is the constituent one). if(_splitGluon==0) { ptrQ = new_ptr(Particle(quark )); ptrQbar = new_ptr(Particle(quark->CC())); } ptrQ ->set5Momentum( momentumQ ); ptrQbar ->set5Momentum( momentumQbar ); } void PartonSplitter::doinit() { Interfaced::doinit(); // calculate the probabilties for the gluon to branch into each quark type // based on the available phase-space, as in fortran. Energy mg=getParticleData(ParticleID::g)->constituentMass(); for( int ix=1; ix<6; ++ix ) { PDPtr quark = getParticleData(ix); Energy pcm = Kinematics::pstarTwoBodyDecay(mg,quark->constituentMass(), quark->constituentMass()); if(pcm>ZERO) _quarkSelector.insert(pcm/GeV,quark); } if(_quarkSelector.empty()) throw InitException() << "At least one quark must have constituent mass less " << "then the constituent mass of the gluon in " << "PartonSplitter::doinit()" << Exception::runerror; } // Method to colour sort the event and calculate the masses of the // pre-clusters // Convention is to have the triplet of the colour singlet first, // then all gluons, then the antitriplet (and repeat for all the // colour-singlets in the event) void PartonSplitter::colourSorted(PVector& tagged) { // Set up the output PVector sorted; // Reset the storage variables for doing the mass-based strangeness // enhancement _colSingletSize.resize(0); _colSingletm2.resize(0); // Variable to exit out of gluon loops bool gluonLoop = false; // Partons left to consider PVector left = tagged; // Loop over all the tagged particles while (int(left.size()) > 0){ // Pick the first particle available PPtr p = left[0]; // Temporary holding pen for momenta Lorentz5Momentum tempMom(ZERO, ZERO, ZERO, ZERO); // Don't do anything if the particle has no colour // Simply add it back into the list of particles if ( !p->coloured() ) { sorted.push_back(p); // nparts is the index of the particle after adding it to the sorted list int nparts = 1; // Add on the last entry of colSingletSize if the vector is not empty // This is essentially the index the particle will have once it // Get placed into the new colour sorted event nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back(); tempMom += p->momentum(); Energy2 singletm2 = tempMom.m2(); // Store the number of particles and mass of the colour-singlet _colSingletSize.push_back(nparts); _colSingletm2.push_back(singletm2); // Remove the particle from the list of particles left left.erase(remove(left.begin(),left.end(), p), left.end()); } // Temporary holding pen for partons PVector temp; // Variable to sum end-point masses i.e. triplets and anti-triplets Energy endPointMass = ZERO; // If the particle in question is a gluon, search up it's antiColourLine // until we get to the triplet. // Note there are situations where we have a gluon loop if ( p->hasColour() && p->hasAntiColour() ){ // Search up its anticolour line until we get to the start particle tPPtr partner = p->antiColourLine()->endParticle(); // Dummy index used to loop over the anticolour line trace tPPtr dummy = partner; // Store the partner particle temp.push_back(partner); tempMom += partner->momentum(); left.erase(remove(left.begin(),left.end(), partner), left.end()); // While loop continues while we still reach a particle with with // anti-colour, i.e. a gluon while ( dummy->hasAntiColour() ){ dummy = dummy->antiColourLine()->endParticle(); // Check that we haven't already included it via colour indices // If we have, it is a gluon loop. if ( find(left.begin(), left.end(), dummy) == left.end() ) { gluonLoop = true; break; } // Store the dummy partons in reverse temp.push_back(dummy); tempMom += dummy->momentum(); // Remove counted ones from the remaining list left.erase(remove(left.begin(),left.end(), dummy), left.end()); } // Number of particles in this colour singlets so far int nparts = int(temp.size()); // Insert the new particles in the reverse order sorted.insert(sorted.end(), temp.rbegin(), temp.rend()); endPointMass += ((temp.back())->mass()); // If it is a gluon loop we've already looped over the colour-singlet // in its entirety, so we need to end early if (gluonLoop){ // Store the index of the final entry nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back(); // Insert the new particles in the correct order // i.e. triplet first, then the gluons we have seen so far Energy2 singletm2 = tempMom.m2(); _colSingletSize.push_back(nparts); _colSingletm2.push_back(singletm2); continue; } // If it is not a gluon loop, we now need to trace the colourLine // down, until we reach the triplet. // Works similarly to the first half // Reset the temp PVector temp.resize(0); // Push back the particle in question temp.push_back(p); // tempMom hasn't been reset, add the particle we started with tempMom += p->momentum(); left.erase(remove(left.begin(),left.end(), p), left.end()); // Search down its colour line until we get to the end particle dummy = p->colourLine()->startParticle(); temp.push_back(dummy); tempMom += dummy->momentum(); left.erase(remove(left.begin(),left.end(), dummy), left.end()); while ( dummy->hasColour() ){ dummy = dummy->colourLine()->startParticle(); temp.push_back(dummy); tempMom += dummy->momentum(); left.erase(remove(left.begin(),left.end(), dummy), left.end()); } endPointMass += ((temp.back())->mass()); // Update size of colour singlet nparts += int(temp.size()); nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back(); // Insert the new particles in the correct order Energy2 singletm2 = tempMom.m2(); sorted.insert(sorted.end(), temp.begin(), temp.end()); endPointMass += ((sorted.back())->mass()); _colSingletSize.push_back(nparts); // Chooses whether to use invariant mass of singlet // or to use the lambda measure i.e. m^2 - (\sum m_i)^2 Energy2 m2 = (_massMeasure == 0) ? singletm2 : singletm2 - sqr(endPointMass); _colSingletm2.push_back(m2); } // Else if it's a quark else if ( p->hasColour() ) { // Search up its colour line until we get to the start particle // Works the same way as the second half of the gluon handling tPPtr partner = p->colourLine()->startParticle(); tPPtr dummy = partner; temp.push_back(p); endPointMass += ((temp.back())->mass()); temp.push_back(partner); tempMom += p->momentum(); tempMom += partner->momentum(); left.erase(remove(left.begin(),left.end(), p), left.end()); left.erase(remove(left.begin(),left.end(), partner), left.end()); while ( dummy->hasColour() ){ dummy = dummy->colourLine()->startParticle(); temp.push_back(dummy); tempMom += dummy->momentum(); left.erase(remove(left.begin(),left.end(), dummy), left.end()); } // Number of particles in this colour singlets int nparts = int(temp.size()); nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back(); // Insert the new particles in the correct order Energy2 singletm2 = tempMom.m2(); sorted.insert(sorted.end(), temp.begin(), temp.end()); endPointMass += ((sorted.back())->mass()); _colSingletSize.push_back(nparts); Energy2 m2 = (_massMeasure == 0) ? singletm2 : singletm2 - sqr(endPointMass); _colSingletm2.push_back(m2); } // Else it's an antiquark else if ( p->hasAntiColour() ) { // Search along anti-colour line, storing particles, and reversing the order // at the end // Works in the same way as the first half of the gluon handling tPPtr partner = p->antiColourLine()->endParticle(); tPPtr dummy = partner; temp.push_back(p); endPointMass += ((temp.back())->mass()); temp.push_back(partner); tempMom += p->momentum(); tempMom += partner->momentum(); left.erase(remove(left.begin(),left.end(), p), left.end()); left.erase(remove(left.begin(),left.end(), partner), left.end()); while ( dummy->hasAntiColour() ){ dummy = dummy->antiColourLine()->endParticle(); temp.push_back(dummy); tempMom += dummy->momentum(); left.erase(remove(left.begin(),left.end(), dummy), left.end()); } // Number of particles in this colour singlets int nparts = int(temp.size()); nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back(); // Insert the particles in the reverse order Energy2 singletm2 = tempMom.m2(); sorted.insert(sorted.end(), temp.rbegin(), temp.rend()); endPointMass += ((temp.back())->mass()); _colSingletSize.push_back(nparts); Energy2 m2 = (_massMeasure == 0) ? singletm2 : singletm2 - sqr(endPointMass); _colSingletm2.push_back(m2); } } // Check that the algorithm hasn't missed any particles. assert( sorted.size() == tagged.size() ); swap(sorted,tagged); } void PartonSplitter::massSplitTimeLikeGluon(tcPPtr ptrGluon, PPtr & ptrQ, PPtr & ptrQbar, size_t i){ // select the quark flavour tPDPtr quark; long idNew=0; if ( ptrGluon->momentum().m() < 2.0 *getParticle(ThePEG::ParticleID::s)->data().constituentMass() ) { throw Exception() << "Impossible Kinematics in PartonSplitter::massSplitTimeLikeGluon()" << Exception::runerror; } // Only allow light quarks u,d,s with the probabilities double prob_d = _splitPwtDquark; double prob_u = _splitPwtUquark; double prob_s = enhanceStrange(i); int choice = UseRandom::rnd3(prob_u, prob_d, prob_s); switch(choice) { case 0: idNew = ThePEG::ParticleID::u; break; case 1: idNew = ThePEG::ParticleID::d; break; case 2: idNew = ThePEG::ParticleID::s; break; } ptrQ = getParticle(idNew); ptrQbar = getParticle(-idNew); // Solve the kinematics of the two body decay G --> Q + Qbar Lorentz5Momentum momentumQ; Lorentz5Momentum momentumQbar; double cosThetaStar = UseRandom::rnd( -1.0 , 1.0 ); using Constants::pi; double phiStar = UseRandom::rnd( -pi , pi ); Energy constituentQmass; constituentQmass = ptrQ->data().constituentMass(); if (ptrGluon->momentum().m() < 2.0*constituentQmass) { throw Exception() << "Impossible Kinematics in PartonSplitter::massSplitTimeLikeGluon()" << Exception::eventerror; } Kinematics::twoBodyDecay(ptrGluon->momentum(), constituentQmass, constituentQmass, cosThetaStar, phiStar, momentumQ, momentumQbar ); // Create quark and anti-quark particles of the chosen flavour // and set they 5-momentum (the mass is the constituent one). ptrQ ->set5Momentum( momentumQ ); ptrQbar ->set5Momentum( momentumQbar ); } double PartonSplitter::enhanceStrange(size_t i){ // Get the m2 of the relevant colour singlet // First we need to get the index of the colour-singlet the indexed i // parton has auto const it = lower_bound(_colSingletSize.begin(), _colSingletSize.end(), i); // Get the index of the colourSinglet mass int indx = distance(_colSingletSize.begin(), it); Energy2 mass2 = _colSingletm2[indx]; Energy2 m2 = _m0*_m0; // Scaling strangeness enhancement if (_enhanceSProb == 1){ // If the mass is significantly smaller than the characteristic mass, // just set the prob to 0 double scale = double(m2/mass2); - return (20. < scale || scale < 0.) ? 0. : pow(_splitPwtSquark,scale); + return (_maxScale < scale || scale < 0.) ? 0. : pow(_splitPwtSquark,scale); } // Exponential strangeness enhancement else if (_enhanceSProb == 2){ double scale = double(m2/mass2); - return (20. < scale || scale < 0.) ? 0. : exp(-scale); + return (_maxScale < scale || scale < 0.) ? 0. : exp(-scale); } else return _splitPwtSquark; } diff --git a/Hadronization/PartonSplitter.h b/Hadronization/PartonSplitter.h --- a/Hadronization/PartonSplitter.h +++ b/Hadronization/PartonSplitter.h @@ -1,222 +1,232 @@ // -*- C++ -*- // // PartonSplitter.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_PartonSplitter_H #define HERWIG_PartonSplitter_H #include "CluHadConfig.h" #include #include #include "PartonSplitter.fh" namespace Herwig { using namespace ThePEG; /** \ingroup Hadronization * \class PartonSplitter * \brief This class splits the gluons from the end of the shower. * \author Philip Stephens * \author Alberto Ribon * * This class does all of the nonperturbative parton splittings needed * immediately after the end of the showering (both initial and final), * as very first step of the cluster hadronization. * * the quarks are attributed with different weights for the splitting * by default only the splitting in u and d quarks is allowed * the option "set /Herwig/Hadronization/PartonSplitter:Split 1" * allows for additional splitting into s quarks based on some weight * in order for that to work the mass of the strange quark has to be changed * from the default value s.t. m_g > 2m_s * * * * @see \ref PartonSplitterInterfaces "The interfaces" * defined for PartonSplitter. */ class PartonSplitter: public Interfaced { public: /** * Default constructor */ PartonSplitter() : _splitPwtUquark(1), _splitPwtDquark(1), _splitPwtSquark(0.5), _gluonDistance(ZERO), _splitGluon(0), _enhanceSProb(0), _m0(10.*GeV), _massMeasure(0) {} /** * This method does the nonperturbative splitting of: * time-like gluons. At the end of the shower the gluons should be * on a "physical" mass shell and should therefore be time-like. * @param tagged The tagged particles to be split * @return The particles which were not split and the products of splitting. */ void split(PVector & tagged); 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; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} private: /** * Private and non-existent assignment operator. */ PartonSplitter & operator=(const PartonSplitter &) = delete; /** * Non-perturbatively split a time-like gluon, * if something goes wrong null pointers are returned. * @param gluon The gluon to be split * @param quark The quark produced in the splitting * @param anti The antiquark produced in the splitting */ void splitTimeLikeGluon(tcPPtr gluon, PPtr & quark, PPtr & anti); /** * Non-perturbatively split a time-like gluon using a mass-based * strangeness enhancement, * if something goes wrong null pointers are returned. * @param gluon The gluon to be split * @param quark The quark produced in the splitting * @param anti The antiquark produced in the splitting * @param gluon's index in the colour-sorted Particle vector */ void massSplitTimeLikeGluon(tcPPtr gluon, PPtr & quark, PPtr & anti, size_t i); /** * Colour-sort the event into grouped colour-singlets. * Convention: triplet - gluons - antitriplet, repeat for * all other colour-singlets or colourless particles. */ void colourSorted(PVector& tagged); /** * Method to calculate the probability for producing strange quarks */ double enhanceStrange(size_t i); private: - // probabilities for the different quark types + /** + * Different variables for the weight for splitting into various + * light quark species. + */ double _splitPwtUquark; double _splitPwtDquark; double _splitPwtSquark; /** * The selector to pick the type of quark */ Selector _quarkSelector; /** * c tau for gluon decays */ Length _gluonDistance; /** * Flag used to determine between normal gluon splitting and alternative gluon splitting */ int _splitGluon; /** * Vector that stores the index in the event record of the anti-triplet * of the colour-singlets, or of a colourless particle */ vector _colSingletSize; /** * Vector that stores the masses of the colour-singlets or of a colourless * particle */ vector _colSingletm2; /** * Option to choose which functional form to use for strangeness * production */ int _enhanceSProb; /** - * Characteristic mass scale for mass-based strangeness enhancement + * Characteristic mass scale for strangeness enhancement */ Energy _m0; /** - * Option to choose which mass measure to use - */ + * 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 that stops the scale in strangeness enhancement from + * becoming too large + */ + const double _maxScale = 20.; + }; } #endif /* HERWIG_PartonSplitter_H */