diff --git a/Hadronization/ClusterFissioner.cc b/Hadronization/ClusterFissioner.cc --- a/Hadronization/ClusterFissioner.cc +++ b/Hadronization/ClusterFissioner.cc @@ -1,1164 +1,1188 @@ // -*- 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 "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), _probPowFactor(4.0), _probShift(0.0), - _kinThresholdShift(1.0*sqr(GeV)) + _kinThresholdShift(1.0*sqr(GeV)), + _kinematicThresholdChoice(0) {} IBPtr ClusterFissioner::clone() const { return new_ptr(*this); } IBPtr ClusterFissioner::fullclone() const { return new_ptr(*this); } void ClusterFissioner::persistentOutput(PersistentOStream & os) const { os << _hadronSelector << 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) + << ounit(_btClM,GeV) << _kinematicThresholdChoice << _iopRem << ounit(_kappa, GeV/meter) << _enhanceSProb << ounit(_m0Fission,GeV) << _massMeasure << _probPowFactor << _probShift << ounit(_kinThresholdShift,sqr(GeV)); } void ClusterFissioner::persistentInput(PersistentIStream & is, int) { is >> _hadronSelector >> 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(_btClM,GeV) >> _iopRem >> _kinematicThresholdChoice >> iunit(_kappa, GeV/meter) >> _enhanceSProb >> iunit(_m0Fission,GeV) >> _massMeasure >> _probPowFactor >> _probShift >> iunit(_kinThresholdShift,sqr(GeV)); } 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::_hadronSelector, 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); static Parameter interfaceProbPowFactor ("ProbablityPowerFactor", "Power factor in ClausterFissioner bell probablity function", &ClusterFissioner::_probPowFactor, 2.0, 1.0, 20.0, false, false, Interface::limited); static Parameter interfaceProbShift ("ProbablityShift", "Shifts from the center in ClausterFissioner bell probablity function", &ClusterFissioner::_probShift, 0.0, -10.0, 10.0, false, false, Interface::limited); static Parameter interfaceKineticThresholdShift ("KineticThresholdShift", "Shifts from the kinetic threshold in ClausterFissioner", &ClusterFissioner::_kinThresholdShift, sqr(GeV), 0.*sqr(GeV), -10.0*sqr(GeV), 10.0*sqr(GeV), false, false, Interface::limited); + static Switch interfaceKinematicThreshold + ("KinematicThreshold", + "Option for using static or dynamic kinematic thresholds in cluster splittings", + &ClusterFissioner::_kinematicThresholdChoice, 0, false, false); + static SwitchOption interfaceKinematicThresholdStatic + (interfaceKinematicThreshold, + "Static", + "Set static kinematic thresholds for cluster splittings.", + 0); + static SwitchOption interfaceKinematicThresholdDynamic + (interfaceKinematicThreshold, + "Dynamic", + "Set dynamic kinematic thresholds for cluster splittings.", + 1); + } 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; // 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; pQ1.setMass(m1); pQone.setMass(m); pQtwo.setMass(m); pQ2.setMass(m2); pair res = drawNewMasses(Mc, soft1, soft2, pClu1, pClu2, ptrQ1, pQ1, newPtr1, pQone, newPtr2, pQtwo, ptrQ2, pQ2); // derive the masses of the children Mc1 = res.first; Mc2 = res.second; + // static kinematic threshold + if(_kinematicThresholdChoice == 0) + if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > Mc) continue; + // dynamic kinematic threshold + else if(_kinematicThresholdChoice == 1) { + bool C1 = ( sqr(Mc1) )/( sqr(m1) + sqr(m) + _kinThresholdShift ) < 1.0 ? true : false; + bool C2 = ( sqr(Mc2) )/( sqr(m2) + sqr(m) + _kinThresholdShift ) < 1.0 ? true : false; + bool C3 = ( sqr(Mc1) + sqr(Mc2) )/( sqr(Mc) ) > 1.0 ? true : false; - // old kinematic threshold - //if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > Mc) continue; - - // new kinematic threshold - bool C1 = ( sqr(Mc1) )/( sqr(m1) + sqr(m) + _kinThresholdShift ) < 1.0 ? true : false; - bool C2 = ( sqr(Mc2) )/( sqr(m2) + sqr(m) + _kinThresholdShift ) < 1.0 ? true : false; - bool C3 = ( sqr(Mc1) + sqr(Mc2) )/( sqr(Mc) ) > 1.0 ? true : false; - - if( C1 || C2 || C3 ) continue; + if( C1 || C2 || C3 ) 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 = _hadronSelector->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc1); if(toHadron1) { Mc1 = toHadron1->mass(); pClu1.setMass(Mc1); } toHadron2 = _hadronSelector->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc2); 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 = _hadronSelector->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc-Mc2); if(toHadron1) { Mc1 = toHadron1->mass(); pClu1.setMass(Mc1); } } else if(!toHadron2) { toHadron2 = _hadronSelector->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc-Mc1); 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) calculateKinematics(pClu,p0Q1,toHadron1,toHadron2, 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; // 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; pQ1.setMass(m1); pQone.setMass(m); pQtwo.setMass(m); pQ2.setMass(m2); 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 = _hadronSelector->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),Mc1); if(toHadron) { Mc1 = toHadron->mass(); pClu1.setMass(Mc1); } // check if need to force diquark cluster to be on-shell toDiQuark = false; diquark = _hadronSelector->makeDiquark(ptrQ[iq2]->dataPtr(),newPtr2->dataPtr()); if(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 = _hadronSelector->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),mmax-Mc2); if(toHadron) { Mc1 = toHadron->mass(); pClu1.setMass(Mc1); } } else if(!toDiQuark) { Mc2 = _hadronSelector->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(); 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 = (_hadronSelector->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 = _hadronSelector->pwt(1); prob_u = _hadronSelector->pwt(2); prob_s = _hadronSelector->pwt(3); 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 = _hadronSelector->pwt(1); prob_u = _hadronSelector->pwt(2); /* Strangeness enhancement: Case 1: probability scaling Case 2: Exponential scaling */ if (_enhanceSProb == 1) prob_s = (_maxScale < scale) ? 0. : pow(_hadronSelector->pwt(3),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; 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()); } 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. Energy2 mag2 = u.vect().mag2(); InvEnergy fact = mag2>ZERO ? 1./sqrt(mag2) : 1./GeV; Length x1 = ( 0.25*Mclu + 0.5*( pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa; Length t1 = Mclu/_kappa - x1; LorentzDistance distanceClu1( x1 * fact * u.vect(), 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 * fact * u.vect(), 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::ProbablityFunction(double scale, double threshold) { double cut = UseRandom::rnd(0.0,1.0); return 1./(1.+pow(abs((threshold-_probShift)/scale),_probPowFactor)) > cut ? true : false; } bool ClusterFissioner::isHeavy(tcClusterPtr clu) { // default double clpow = _clPowLight; Energy clmax = _clMaxLight; // particle data for constituents tcPDPtr cptr[3]={tcPDPtr(),tcPDPtr(),tcPDPtr()}; for(unsigned 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; } - //some smooth probablity function to create a dynamic thershold - double scale = pow(clu->mass()/GeV , clpow); - double threshold = pow(clmax/GeV, clpow) - + pow(clu->sumConstituentMasses()/GeV, clpow); - bool aboveCutoff = ProbablityFunction(scale,threshold); - //regular checks - // 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 aboveCutoff, canSplitMinimally; + // static kinematic threshold + if(_kinematicThresholdChoice == 0) { + aboveCutoff = ( + pow(clu->mass()*UnitRemoval::InvE , clpow) + > + pow(clmax*UnitRemoval::InvE, clpow) + + pow(clu->sumConstituentMasses()*UnitRemoval::InvE, clpow) + ); - scale = clu->mass()/GeV; - threshold = clu->sumConstituentMasses()/GeV + 2.0 * minmass/GeV; + canSplitMinimally = clu->mass() > clu->sumConstituentMasses() + 2.0 * minmass; + } + // dynamic kinematic threshold + else if(_kinematicThresholdChoice == 1) { + //some smooth probablity function to create a dynamic thershold + double scale = pow(clu->mass()/GeV , clpow); + double threshold = pow(clmax/GeV, clpow) + + pow(clu->sumConstituentMasses()/GeV, clpow); + bool aboveCutoff = ProbablityFunction(scale,threshold); - bool canSplitMinimally = ProbablityFunction(scale,threshold); - //bool canSplitMinimally = clu->mass() > clu->sumConstituentMasses() + 2.0 * minmass; + scale = clu->mass()/GeV; + threshold = clu->sumConstituentMasses()/GeV + 2.0 * minmass/GeV; + + canSplitMinimally = ProbablityFunction(scale,threshold); + } return aboveCutoff && canSplitMinimally; } diff --git a/Hadronization/ClusterFissioner.h b/Hadronization/ClusterFissioner.h --- a/Hadronization/ClusterFissioner.h +++ b/Hadronization/ClusterFissioner.h @@ -1,511 +1,516 @@ // -*- 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 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; /** * 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; /** * 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; protected: /** @name Access members for child classes. */ //@{ /** * Access to the hadron selector */ HadronSelectorPtr hadronSelector() const {return _hadronSelector;} /** * 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: /** * Smooth probability for dynamic threshold cuts: * @scale the current scale, e.g. the mass of the cluster, * @threshold the physical threshold, */ bool ProbablityFunction(double scale, double threshold); /** * Check if a cluster is heavy enough to split again */ bool isHeavy(tcClusterPtr ); /** * A pointer to a Herwig::HadronSelector object for generating hadrons. */ HadronSelectorPtr _hadronSelector; /** * @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; + /** + * Flag to choose static or dynamic kinematic thresholds in cluster splittings + */ + int _kinematicThresholdChoice; + //@} /** * 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.; /** * Power factor in ClausterFissioner bell probablity function */ double _probPowFactor; /** * Shifts from the center in ClausterFissioner bell probablity function */ double _probShift; /** * Shifts from the kinetic threshold in ClausterFissioner */ Energy2 _kinThresholdShift; }; } #endif /* HERWIG_ClusterFissioner_H */