diff --git a/Hadronization/CheckId.cc b/Hadronization/CheckId.cc --- a/Hadronization/CheckId.cc +++ b/Hadronization/CheckId.cc @@ -1,160 +1,164 @@ // -*- C++ -*- // // CheckId.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the CheckId class. // #include "CheckId.h" #include using namespace Herwig; namespace { - + /** - * Return true if the particle pointer corresponds to a diquark + * Return true if the particle pointer corresponds to a diquark * or anti-diquark carrying b flavour; false otherwise. */ inline bool isDiquarkWithB(tcPDPtr par1) { if (!par1) return false; long id1 = par1->id(); return DiquarkMatcher::Check(id1) && (abs(id1)/1000)%10 == ParticleID::b; } - + /** * Return true if the particle pointer corresponds to a diquark * or anti-diquark carrying c flavour; false otherwise. */ inline bool isDiquarkWithC(tcPDPtr par1) { if (!par1) return false; long id1 = par1->id(); - return ( DiquarkMatcher::Check(id1) && - ( (abs(id1)/1000)%10 == ParticleID::c + return ( DiquarkMatcher::Check(id1) && + ( (abs(id1)/1000)%10 == ParticleID::c || (abs(id1)/100)%10 == ParticleID::c ) ); } } -long CheckId::makeDiquarkID(long id1, long id2) { +long CheckId::makeDiquarkID(long id1, long id2, long pspin) { - assert( id1 * id2 > 0 - && QuarkMatcher::Check(id1) + assert( id1 * id2 > 0 + && QuarkMatcher::Check(id1) && QuarkMatcher::Check(id2)) ; long ida = abs(id1); long idb = abs(id2); if (ida < idb) swap(ida,idb); - - long idnew = ida*1000 + idb*100 + 1; - // Diquarks made of quarks of the same type: uu, dd, ss, cc, bb, + if (pspin != 1 && pspin != 3) assert(false); + long idnew = ida*1000 + idb*100 + pspin; + // Diquarks made of quarks of the same type: uu, dd, ss, cc, bb, // have spin 1, and therefore the less significant digit (which // corresponds to 2*J+1) is 3 rather than 1 as all other Diquarks. - if (id1 == id2) idnew += 2; + if (id1 == id2 && pspin == 1) { + cerr<<"WARNING: spin-0 diquiark of the same type cannot exist." + <<" Switching to spin-1 diquark.\n"; + idnew = ida*1000 + idb*100 + 3; + } return id1 > 0 ? idnew : -idnew; } -PDPtr CheckId::makeDiquark(tcPDPtr par1, tcPDPtr par2) { +PDPtr CheckId::makeDiquark(tcPDPtr par1, tcPDPtr par2, long pspin) { long id1 = par1->id(); long id2 = par2->id(); - long idnew = makeDiquarkID(id1,id2); + long idnew = makeDiquarkID(id1,id2, pspin); assert(!CurrentGenerator::isVoid()); return CurrentGenerator::current().getParticleData(idnew); } bool CheckId::canBeMeson(tcPDPtr par1,tcPDPtr par2) { assert(par1 && par2); long id1 = par1->id(); long id2 = par2->id(); // a Meson must not have any diquarks if(DiquarkMatcher::Check(id1) || DiquarkMatcher::Check(id2)) return false; - return ( abs(int(par1->iColour()))== 3 && - abs(int(par2->iColour())) == 3 && + return ( abs(int(par1->iColour()))== 3 && + abs(int(par2->iColour())) == 3 && id1*id2 < 0); } bool CheckId::canBeBaryon(tcPDPtr par1, tcPDPtr par2 , tcPDPtr par3) { assert(par1 && par2); long id1 = par1->id(), id2 = par2->id(); if (!par3) { if( id1*id2 < 0) return false; if(DiquarkMatcher::Check(id1)) -return abs(int(par2->iColour())) == 3 && !DiquarkMatcher::Check(id2); +return abs(int(par2->iColour())) == 3 && !DiquarkMatcher::Check(id2); if(DiquarkMatcher::Check(id2)) return abs(int(par1->iColour())) == 3; return false; - } + } else { // In this case, to be a baryon, all three components must be (anti-)quarks // and with the same sign. return (par1->iColour() == 3 && par2->iColour() == 3 && par3->iColour() == 3) || (par1->iColour() == -3 && par2->iColour() == -3 && par3->iColour() == -3); } } - + bool CheckId::hasBottom(tcPDPtr par1, tcPDPtr par2, tcPDPtr par3) { long id1 = par1 ? par1->id() : 0; if ( !par2 && !par3 ) { - return + return abs(id1) == ThePEG::ParticleID::b || isDiquarkWithB(par1) || - ( MesonMatcher::Check(id1) + ( MesonMatcher::Check(id1) && (abs(id1)/100)%10 == ThePEG::ParticleID::b ) || - ( BaryonMatcher::Check(id1) + ( BaryonMatcher::Check(id1) && (abs(id1)/1000)%10 == ThePEG::ParticleID::b ); - } + } else { long id2 = par2 ? par2->id() : 0; long id3 = par3 ? par3->id() : 0; - return - abs(id1) == ThePEG::ParticleID::b || isDiquarkWithB(par1) || - abs(id2) == ThePEG::ParticleID::b || isDiquarkWithB(par2) || - abs(id3) == ThePEG::ParticleID::b || isDiquarkWithB(par3); + return + abs(id1) == ThePEG::ParticleID::b || isDiquarkWithB(par1) || + abs(id2) == ThePEG::ParticleID::b || isDiquarkWithB(par2) || + abs(id3) == ThePEG::ParticleID::b || isDiquarkWithB(par3); } } bool CheckId::hasCharm(tcPDPtr par1, tcPDPtr par2, tcPDPtr par3) { long id1 = par1 ? par1->id(): 0; if (!par2 && !par3) { return abs(id1) == ThePEG::ParticleID::c || isDiquarkWithC(par1) || - ( MesonMatcher::Check(id1) && + ( MesonMatcher::Check(id1) && ((abs(id1)/100)%10 == ThePEG::ParticleID::c || (abs(id1)/10)%10 == ThePEG::ParticleID::c) ) || - ( BaryonMatcher::Check(id1) && + ( BaryonMatcher::Check(id1) && ((abs(id1)/1000)%10 == ThePEG::ParticleID::c || (abs(id1)/100)%10 == ThePEG::ParticleID::c || (abs(id1)/10)%10 == ThePEG::ParticleID::c) ); - } + } else { long id2 = par2 ? par1->id(): 0; long id3 = par3 ? par1->id(): 0; - return - abs(id1) == ThePEG::ParticleID::c || isDiquarkWithC(par1) || - abs(id2) == ThePEG::ParticleID::c || isDiquarkWithC(par2) || - abs(id3) == ThePEG::ParticleID::c || isDiquarkWithC(par3); + return + abs(id1) == ThePEG::ParticleID::c || isDiquarkWithC(par1) || + abs(id2) == ThePEG::ParticleID::c || isDiquarkWithC(par2) || + abs(id3) == ThePEG::ParticleID::c || isDiquarkWithC(par3); } -} +} bool CheckId::isExotic(tcPDPtr par1, tcPDPtr par2, tcPDPtr par3) { /// \todo make this more general long id1 = par1 ? par1->id(): 0; long id2 = par2 ? par2->id(): 0; long id3 = par3 ? par3->id(): 0; -return +return ( (id1/1000000)% 10 != 0 && (id1/1000000)% 10 != 9 ) || ( (id2/1000000)% 10 != 0 && (id2/1000000)% 10 != 9 ) || ( (id3/1000000)% 10 != 0 && (id3/1000000)% 10 != 9 ) || abs(id1)==6||abs(id2)==6; } diff --git a/Hadronization/CheckId.h b/Hadronization/CheckId.h --- a/Hadronization/CheckId.h +++ b/Hadronization/CheckId.h @@ -1,108 +1,108 @@ // -*- C++ -*- // // CheckId.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_CheckId_H #define HERWIG_CheckId_H // // This is the declaration of the CheckId class. #include "ThePEG/PDT/StandardMatchers.h" #include "ThePEG/PDT/ParticleData.h" #include #include "ThePEG/Repository/CurrentGenerator.h" namespace Herwig { using namespace ThePEG; /** \ingroup Hadronization * * This is a namespace which provides some useful methods * for checking the PDG id of particles. - * Notice that the name "quark" in the methods below means any of + * Notice that the name "quark" in the methods below means any of * the following: * d , u , s , c , b anti-d , anti-u , anti-s , anti-c , anti-b - * that is we do not distinguish between particle or antiparticle - * ( because can be done directly by the user: + * that is we do not distinguish between particle or antiparticle + * ( because can be done directly by the user: * id > 0 for particles ; id < 0 for anti-particles ) * and we do not include t (and anti-t ) because we are interested * in components of mesons and baryons. * Similarly for the name "diquark" which include all diquarks (id > 0) - * and anti-diquarks (id < 0) not made with t ( anti-t ) component. + * and anti-diquarks (id < 0) not made with t ( anti-t ) component. * * NB) For Other useful methods (even some implemented in CheckId) * @see StandardMatchers */ namespace CheckId { - + /** - * Return the id of the diquark (anti-diquark) made by the two + * Return the id of the diquark (anti-diquark) made by the two * quarks (antiquarks) of id specified in input (id1, id2). * Caller must ensure that id1 and id2 are quarks. */ - long makeDiquarkID(long id1, long id2); - + long makeDiquarkID(long id1, long id2, long pspin); + /** - * Return the particle data of the diquark (anti-diquark) made by the two + * Return the particle data of the diquark (anti-diquark) made by the two * quarks (antiquarks) par1, par2. * @param par1 (anti-)quark data pointer * @param par2 (anti-)quark data pointer */ - PDPtr makeDiquark(tcPDPtr par1, tcPDPtr par2); + PDPtr makeDiquark(tcPDPtr par1, tcPDPtr par2, long pspin); /** * Return true if the two particles in input can be the components of a meson; *false otherwise. */ bool canBeMeson(tcPDPtr par1,tcPDPtr par2); /** - * Return true if the two or three particles in input can be the components + * Return true if the two or three particles in input can be the components * of a baryon; false otherwise. */ bool canBeBaryon(tcPDPtr par1, tcPDPtr par2 , tcPDPtr par3 = PDPtr()); - + /** - * Return true if the two or three particles in input can be the components + * Return true if the two or three particles in input can be the components * of a hadron; false otherwise. */ inline bool canBeHadron(tcPDPtr par1, tcPDPtr par2 , tcPDPtr par3 = PDPtr()) { return (!par3 && canBeMeson(par1,par2)) || canBeBaryon(par1,par2,par3); } - + /** * Return true if any of the possible three input particles has - * b-flavour; + * b-flavour; * false otherwise. In the case that only the first particle is specified, * it can be: an (anti-)quark, an (anti-)diquark * an (anti-)meson, an (anti-)baryon; in the other cases, each pointer * is assumed to be either (anti-)quark or (anti-)diquark. */ bool hasBottom(tcPDPtr par1, tcPDPtr par2 = PDPtr(), tcPDPtr par3 = PDPtr()); /** - * Return true if any of the possible three input particles has - * c-flavour; + * Return true if any of the possible three input particles has + * c-flavour; * false otherwise.In the case that only the first pointer is specified, * it can be: a (anti-)quark, a (anti-)diquark * a (anti-)meson, a (anti-)baryon; in the other cases, each pointer * is assumed to be either (anti-)quark or (anti-)diquark. */ bool hasCharm(tcPDPtr par1, tcPDPtr par2 = PDPtr(), tcPDPtr par3 = PDPtr()); /** * Return true, if any of the possible input particle pointer is an exotic quark, e.g. Susy quark; - * false otherwise. + * false otherwise. */ bool isExotic(tcPDPtr par1, tcPDPtr par2 = PDPtr(), tcPDPtr par3 = PDPtr()); } - + } #endif /* HERWIG_CheckId_H */ diff --git a/Hadronization/ClusterFinder.cc b/Hadronization/ClusterFinder.cc --- a/Hadronization/ClusterFinder.cc +++ b/Hadronization/ClusterFinder.cc @@ -1,485 +1,486 @@ // -*- C++ -*- // // ClusterFinder.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ClusterFinder class. // #include "ClusterFinder.h" #include #include #include #include #include #include #include #include #include "CheckId.h" #include "Herwig/Utilities/EnumParticles.h" #include "Herwig/Utilities/Kinematics.h" #include "Cluster.h" #include using namespace Herwig; DescribeClass describeClusterFinder("Herwig::ClusterFinder","Herwig.so"); IBPtr ClusterFinder::clone() const { return new_ptr(*this); } IBPtr ClusterFinder::fullclone() const { return new_ptr(*this); } void ClusterFinder::persistentOutput(PersistentOStream & os) const { os << heavyDiquarks_ << diQuarkSelection_ << diQuarkOnShell_; } void ClusterFinder::persistentInput(PersistentIStream & is, int) { is >> heavyDiquarks_ >> diQuarkSelection_ >> diQuarkOnShell_; } void ClusterFinder::Init() { static ClassDocumentation documentation ("This class is responsible of finding clusters."); static Switch interfaceHeavyDiquarks ("HeavyDiquarks", "How to treat heavy quarks in baryon number violating clusters", &ClusterFinder::heavyDiquarks_, 2, false, false); static SwitchOption interfaceHeavyDiquarksDefault (interfaceHeavyDiquarks, "Allow", "No special treatment, allow both heavy and doubly heavy diquarks", 0); static SwitchOption interfaceHeavyDiquarksNoDoublyHeavy (interfaceHeavyDiquarks, "NoDoublyHeavy", "Avoid having diquarks with twoo heavy quarks", 1); static SwitchOption interfaceHeavyDiquarksNoHeavy (interfaceHeavyDiquarks, "NoHeavy", "Try and avoid diquarks contain c and b altogether", 2); static Switch interfaceDiQuarkSelection ("DiQuarkSelection", "Option controlling the selection of quarks to merge into a diquark in baryon-number violating clusters", &ClusterFinder::diQuarkSelection_, 1, false, false); static SwitchOption interfaceDiQuarkSelectionRandom (interfaceDiQuarkSelection, "Random", "Randomly pick a pair to combine", 0); static SwitchOption interfaceDiQuarkSelectionLowestMass (interfaceDiQuarkSelection, "LowestMass", "Combine the lowest mass pair", 1); static Switch interfaceDiQuarkOnShell ("DiQuarkOnShell", "Force the diquark produced in baryon-number violating clusters to be on-shell", &ClusterFinder::diQuarkOnShell_, false, false, false); static SwitchOption interfaceDiQuarkOnShellYes (interfaceDiQuarkOnShell, "Yes", "Force to be on-shell", true); static SwitchOption interfaceDiQuarkOnShellNo (interfaceDiQuarkOnShell, "No", "Leave off-shell", false); } ClusterVector ClusterFinder::formClusters(const PVector & partons) { set examinedSet; // colour particles already included in a cluster - map > quarkQuark; // quark quark + map > quarkQuark; // quark quark map > aQuarkQuark; // anti quark anti quark ParticleSet inputParticles(partons.begin(),partons.end()); ClusterVector clusters; - + // Loop over all current particles. for(PVector::const_iterator pit=partons.begin();pit!=partons.end();++pit){ // Skip to the next particle if it is not coloured or already examined. assert(*pit); assert((*pit)->dataPtr()); - if(!(**pit).data().coloured() + if(!(**pit).data().coloured() || examinedSet.find(*pit) != examinedSet.end()) { - continue; + continue; } // We assume that a cluster is made of, at most, 3 constituents although // in most cases the number will be 2; however, for baryon violating decays // (for example in Susy model without R parity conservation) we can have 3 // constituents. In the latter case, a quark (antiquark) do not have an - // anticolour (colour) partner as usual, but its colour line either stems + // anticolour (colour) partner as usual, but its colour line either stems // from a colour source, or ends in a colour sink. In the case of double - // baryon violating decays, but with overall baryon conservation + // baryon violating decays, but with overall baryon conservation // ( for instance: - // tilde_u_R -> dbar_1 + dbar_2 - // tilde_u_R_star -> d1 + d2 + // tilde_u_R -> dbar_1 + dbar_2 + // tilde_u_R_star -> d1 + d2 // where tilde_u_R and tilde_u_R_star are colour connected ) // a special treatment is needed, because first we have to process all // partons in the current step, and then for each left pair of quarks which - // stem from a colour source we have to find the corresponding pair of + // stem from a colour source we have to find the corresponding pair of // anti-quarks which ends in a colour sink and is connected with the // above colour source. These special pairs are kept into the maps: // spec/CluHadConfig.hialQuarkQuarkMap and specialAntiQuarkAntiQuarkMap. tParticleVector connected(3); - int iElement = 0; + int iElement = 0; connected[iElement++] = *pit; bool specialCase = false; if((*pit)->hasColour()) { - tPPtr partner = + tPPtr partner = (*pit)->colourLine()->getColouredParticle(partons.begin(), partons.end(), true); if(partner) { connected[iElement++]= partner; - } + } // colour source : baryon-violating process else { if((*pit)->colourLine()->sourceNeighbours() != tColinePair()) { tColinePair sourcePair = (*pit)->colourLine()->sourceNeighbours(); tColinePtr intCL = tColinePtr(); for(int i = 0; i < 2; ++i) { tColinePtr pLine = i==0 ? sourcePair.first : sourcePair.second; int saveNumElements = iElement; - for(tPVector::const_iterator cit = pLine->coloured().begin(); + for(tPVector::const_iterator cit = pLine->coloured().begin(); cit != pLine->coloured().end(); ++cit ) { ParticleSet::const_iterator cjt = inputParticles.find(*cit); if(cjt!=inputParticles.end()) connected[iElement++]= (*cit); - } + } if(iElement == saveNumElements) intCL = pLine; } if(intCL && iElement == 2) { specialCase = true; pair qp=pair(connected[0],connected[1]); - quarkQuark.insert(pair >(intCL,qp)); + quarkQuark.insert(pair >(intCL,qp)); } else if(iElement != 3) { - throw Exception() << "Colour connections fail in the hadronization for " + throw Exception() << "Colour connections fail in the hadronization for " << **pit << "in ClusterFinder::formClusters" << " for a coloured particle." << " Failed to find particles from a source" << Exception::runerror; } } else { - throw Exception() << "Colour connections fail in the hadronization for " + throw Exception() << "Colour connections fail in the hadronization for " << **pit << "in ClusterFinder::formClusters for" << " a coloured particle" << Exception::runerror; } } } if((*pit)->hasAntiColour()) { - tPPtr partner = + tPPtr partner = (*pit)->antiColourLine()->getColouredParticle(partons.begin(), partons.end(), false); if(partner) { - connected[iElement++]=partner; + connected[iElement++]=partner; } // colour sink : baryon-violating process else { if((*pit)->antiColourLine()->sinkNeighbours() != tColinePair()) { tColinePair sinkPair = (*pit)->antiColourLine()->sinkNeighbours(); tColinePtr intCL = tColinePtr(); for(int i = 0; i < 2; ++i) { tColinePtr pLine = i==0 ? sinkPair.first : sinkPair.second; int saveNumElements = iElement; - for(tPVector::const_iterator cit = pLine->antiColoured().begin(); + for(tPVector::const_iterator cit = pLine->antiColoured().begin(); cit != pLine->antiColoured().end(); ++cit ) { ParticleSet::const_iterator cjt = inputParticles.find(*cit); if(cjt!=inputParticles.end()) connected[iElement++]= (*cit); } if(iElement == saveNumElements) intCL = pLine; } if(intCL && iElement == 2) { specialCase = true; pair aqp=pair(connected[0],connected[1]); aQuarkQuark.insert(pair >(intCL,aqp)); } else if( iElement !=3) { throw Exception() << "Colour connections fail in the hadronization for " << **pit << "in ClusterFinder::formClusters for" << " an anti-coloured particle." << " Failed to find particles from a sink" << Exception::runerror; } } else { - throw Exception() << "Colour connections fail in the hadronization for " + throw Exception() << "Colour connections fail in the hadronization for " << **pit << "in ClusterFinder::formClusters for" << " an anti-coloured particle" << Exception::runerror; } } } if(!specialCase) { // Tag the components of the found cluster as already examined. for (int i=0; iaddChild(cluPtr); connected[1]->addChild(cluPtr); if(connected[2]) connected[2]->addChild(cluPtr); clusters.push_back(cluPtr); // Check if any of the components is a beam remnant, and if this - // is the case then inform the cluster. - // this will only work for baryon collisions + // is the case then inform the cluster. + // this will only work for baryon collisions for (int i=0; iid()==ParticleID::Remnant) { fromRemnant = true; break; } parent = parent->parents().empty() ? tPPtr() : parent->parents()[0]; } if(fromRemnant&&DiquarkMatcher::Check(connected[i]->id())) cluPtr->isBeamCluster(connected[i]); } } } // Treat now the special cases, if any. The idea is to find for each pair // of quarks coming from a common colour source the corresponding pair of // antiquarks coming from a common colour sink, connected to the above // colour source via the same colour line. Then, randomly couple one of // the two quarks with one of the two antiquarks, and do the same with the // quark and antiquark left. - for(map >::const_iterator + for(map >::const_iterator cit = quarkQuark.begin(); cit != quarkQuark.end(); ++cit ) { tColinePtr coline = cit->first; pair quarkPair = cit->second; if(aQuarkQuark.find( coline ) != aQuarkQuark.end()) { pair antiQuarkPair = aQuarkQuark.find(coline)->second; ClusterPtr cluPtr1, cluPtr2; - if ( UseRandom::rndbool() ) { + if ( UseRandom::rndbool() ) { cluPtr1 = new_ptr(Cluster(quarkPair.first , antiQuarkPair.first)); cluPtr2 = new_ptr(Cluster(quarkPair.second , antiQuarkPair.second)); quarkPair.first->addChild(cluPtr1); antiQuarkPair.first->addChild(cluPtr1); quarkPair.second->addChild(cluPtr2); antiQuarkPair.second->addChild(cluPtr2); } else { cluPtr1 = new_ptr(Cluster(quarkPair.first , antiQuarkPair.second)); cluPtr2 = new_ptr(Cluster(quarkPair.second , antiQuarkPair.first)); quarkPair.second->addChild(cluPtr2); antiQuarkPair.first->addChild(cluPtr2); quarkPair.first->addChild(cluPtr1); antiQuarkPair.second->addChild(cluPtr1); } clusters.push_back(cluPtr1); clusters.push_back(cluPtr2); - } + } else { - throw Exception() << "ClusterFinder::formClusters : " + throw Exception() << "ClusterFinder::formClusters : " << "***Skip event: unable to match pairs in " << "Baryon-violating processes***" << Exception::eventerror; } } return clusters; } namespace { bool PartOrdering(tPPtr p1,tPPtr p2) { return abs(p1->id())id()); } } void ClusterFinder::reduceToTwoComponents(ClusterVector & clusters) { - // In order to preserve all of the information, we do not modify the + // In order to preserve all of the information, we do not modify the // directly the 3-component clusters, but instead we define new clusters, // which are related to the original ones by a child-parent relationship, // by considering two randomly chosen components as a diquark (or anti-diquark). // These new clusters are first added to the vector vecNewRedefinedCluPtr, - // and at the end, when all input clusters have been examined, the elements of - // this vector will be copied in collecCluPtr (the reason is that it is not + // and at the end, when all input clusters have been examined, the elements of + // this vector will be copied in collecCluPtr (the reason is that it is not // allowed to modify a STL container while iterating over it). - vector redefinedClusters; - for(ClusterVector::iterator cluIter = clusters.begin() ; + vector redefinedClusters; + for(ClusterVector::iterator cluIter = clusters.begin() ; cluIter != clusters.end() ; ++cluIter) { tParticleVector vec; if ( (*cluIter)->numComponents() != 3 || ! (*cluIter)->isAvailable() ) continue; - + tPPtr other; for(unsigned int i = 0; i<(*cluIter)->numComponents(); i++) { tPPtr part = (*cluIter)->particle(i); if(!DiquarkMatcher::Check(*(part->dataPtr()))) vec.push_back(part); else other = part; } if(vec.size()<2) { throw Exception() << "Could not make a diquark for a baryonic cluster decay from " << (*cluIter)->particle(0)->PDGName() << " " << (*cluIter)->particle(1)->PDGName() << " " << (*cluIter)->particle(2)->PDGName() << " " << " in ClusterFinder::reduceToTwoComponents()." << Exception::eventerror; } // order the vector so heaviest at the end std::sort(vec.begin(),vec.end(),PartOrdering); // Special treatment of heavy quarks // avoid doubly heavy diquarks if(heavyDiquarks_>=1 && vec.size()>2 && abs(vec[1]->id())>3 && abs(vec[0]->id())<=3) { if(UseRandom::rndbool()) swap(vec[1],vec[2]); other = vec[2]; vec.pop_back(); } // avoid singly heavy diquarks if(heavyDiquarks_==2 && vec.size()>2 && abs(vec[2]->id())>3 && abs(vec[1]->id())<=3) { other = vec[2]; vec.pop_back(); } // if there's a choice pick the pair to make a diquark from if(vec.size()>2) { unsigned int ichoice(0); // random choice if(diQuarkSelection_==0) { ichoice = UseRandom::rnd3(1.0, 1.0, 1.0); } // pick the lightest quark pair else if(diQuarkSelection_==1) { Energy m12 = (vec[0]->momentum()+vec[1]->momentum()).m(); Energy m13 = (vec[0]->momentum()+vec[2]->momentum()).m(); Energy m23 = (vec[1]->momentum()+vec[2]->momentum()).m(); if (m13<=m12&&m13<=m23) ichoice = 2; else if(m23<=m12&&m23<=m13) ichoice = 1; } else assert(false); // make the swaps so select pair first switch (ichoice) { case 0: break; case 1: swap(vec[2],vec[0]); break; case 2: swap(vec[2],vec[1]); break; } } // set up tcPDPtr temp1 = vec[0]->dataPtr(); tcPDPtr temp2 = vec[1]->dataPtr(); if(!other) other = vec[2]; - tcPDPtr dataDiquark = CheckId::makeDiquark(temp1,temp2); - - if(!dataDiquark) + long rndSpin = UseRandom::rnd() > 0.5 ? 1 : 3; + tcPDPtr dataDiquark = CheckId::makeDiquark(temp1,temp2,rndSpin); + + if(!dataDiquark) throw Exception() << "Could not make a diquark from" << temp1->PDGName() << " and " - << temp2->PDGName() + << temp2->PDGName() << " in ClusterFinder::reduceToTwoComponents()" << Exception::eventerror; - + // Create the new cluster (with two components) and assign to it the same // momentum and position of the original (with three components) one. // Furthermore, assign to the diquark component a momentum given by the // sum of the two original components from which has been formed; for the // position, we are assuming, very simply, that the diquark position is // the average positions of the two original components. // Notice that the mass (5-th component of the 5-momentum) of the diquark // is set by hand to the constituent mass of the diquark (which is equal // to the sum of the constituent masses of the two quarks which form the // diquark) because the sum of 5-component vectors do add only the "normal" // 4-components, not the 5-th one. After that, the 5-momentum of the diquark // is in an inconsistent state, because the mass (5-th component) is not // equal to the invariant mass obtained from the 4-momemtum. This is not // unique to this kind of component (all perturbative components are in // a similar situation), but it is not harmful. - + // construct the diquark PPtr diquark = dataDiquark->produceParticle(); vec[0]->addChild(diquark); vec[1]->addChild(diquark); diquark->set5Momentum(Lorentz5Momentum(vec[0]->momentum() + vec[1]->momentum(), dataDiquark->constituentMass())); // use the same method as for cluster to determine the diquark position diquark->setVertex(Cluster::calculateX(vec[0],vec[1])); // put on-shell if required if(diQuarkOnShell_) { Lorentz5Momentum psum = diquark->momentum()+other->momentum(); psum.rescaleMass(); Boost boost = psum.boostVector(); Lorentz5Momentum pother = other->momentum(); Lorentz5Momentum pdiquark = diquark->momentum(); pother.boost(-boost); pdiquark.boost(-boost); Energy pcm = Kinematics::pstarTwoBodyDecay(psum.mass(), other->dataPtr()->constituentMass(), diquark->dataPtr()->constituentMass()); if(pcm>ZERO) { double fact = pcm/pother.vect().mag(); pother *= fact; - pdiquark *= fact; + pdiquark *= fact; pother .setMass(other->dataPtr()->constituentMass()); pdiquark.setMass(dataDiquark ->constituentMass()); pother .rescaleEnergy(); pdiquark.rescaleEnergy(); pother .boost(boost); pdiquark.boost(boost); other->set5Momentum(pother); diquark->set5Momentum(pdiquark); } } // make the new cluster ClusterPtr nclus = new_ptr(Cluster(other,diquark)); //vec[0]->addChild(nclus); //diquark->addChild(nclus); - // Set the parent/children relationship between the original cluster + // Set the parent/children relationship between the original cluster // (the one with three components) with the new one (the one with two components) // and add the latter to the vector of new redefined clusters. (*cluIter)->addChild(nclus); redefinedClusters.push_back(nclus); } - // Add to collecCluPtr all of the redefined new clusters (indeed the + // Add to collecCluPtr all of the redefined new clusters (indeed the // pointers to them are added) contained in vecNewRedefinedCluPtr. /// \todo why do we keep the original of the redefined clusters? for (tClusterVector::const_iterator it = redefinedClusters.begin(); it != redefinedClusters.end(); ++it) { clusters.push_back(*it); } } diff --git a/Hadronization/ClusterFissioner.cc b/Hadronization/ClusterFissioner.cc --- a/Hadronization/ClusterFissioner.cc +++ b/Hadronization/ClusterFissioner.cc @@ -1,1121 +1,1122 @@ // -*- C++ -*- // // ClusterFissioner.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // Thisk is the implementation of the non-inlined, non-templated member // functions of the ClusterFissioner class. // #include "ClusterFissioner.h" #include #include #include #include #include #include #include #include "Herwig/Utilities/Kinematics.h" #include "CheckId.h" #include "Cluster.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include using namespace Herwig; DescribeClass describeClusterFissioner("Herwig::ClusterFissioner","Herwig.so"); ClusterFissioner::ClusterFissioner() : _clMaxLight(3.35*GeV), _clMaxBottom(3.35*GeV), _clMaxCharm(3.35*GeV), _clMaxExotic(3.35*GeV), _clPowLight(2.0), _clPowBottom(2.0), _clPowCharm(2.0), _clPowExotic(2.0), _pSplitLight(1.0), _pSplitBottom(1.0), _pSplitCharm(1.0), _pSplitExotic(1.0), _fissionPwtUquark(1), _fissionPwtDquark(1), _fissionPwtSquark(0.5), _fissionCluster(0), _btClM(1.0*GeV), _iopRem(1), _kappa(1.0e15*GeV/meter), _enhanceSProb(0), _m0Fission(2.*GeV), _massMeasure(0) {} IBPtr ClusterFissioner::clone() const { return new_ptr(*this); } IBPtr ClusterFissioner::fullclone() const { return new_ptr(*this); } void ClusterFissioner::persistentOutput(PersistentOStream & os) const { os << _hadronsSelector << ounit(_clMaxLight,GeV) << ounit(_clMaxBottom,GeV) << ounit(_clMaxCharm,GeV) << ounit(_clMaxExotic,GeV) << _clPowLight << _clPowBottom << _clPowCharm << _clPowExotic << _pSplitLight << _pSplitBottom << _pSplitCharm << _pSplitExotic << _fissionCluster << _fissionPwtUquark << _fissionPwtDquark << _fissionPwtSquark << ounit(_btClM,GeV) << _iopRem << ounit(_kappa, GeV/meter) << _enhanceSProb << ounit(_m0Fission,GeV) << _massMeasure; } void ClusterFissioner::persistentInput(PersistentIStream & is, int) { is >> _hadronsSelector >> iunit(_clMaxLight,GeV) >> iunit(_clMaxBottom,GeV) >> iunit(_clMaxCharm,GeV) >> iunit(_clMaxExotic,GeV) >> _clPowLight >> _clPowBottom >> _clPowCharm >> _clPowExotic >> _pSplitLight >> _pSplitBottom >> _pSplitCharm >> _pSplitExotic >> _fissionCluster >> _fissionPwtUquark >> _fissionPwtDquark >> _fissionPwtSquark >> iunit(_btClM,GeV) >> _iopRem >> iunit(_kappa, GeV/meter) >> _enhanceSProb >> iunit(_m0Fission,GeV) >> _massMeasure; } void ClusterFissioner::Init() { static ClassDocumentation documentation ("Class responsibles for chopping up the clusters"); static Reference interfaceHadronSelector("HadronSelector", "A reference to the HadronSelector object", &Herwig::ClusterFissioner::_hadronsSelector, false, false, true, false); // ClMax for light, Bottom, Charm and exotic (e.g. Susy) quarks static Parameter interfaceClMaxLight ("ClMaxLight","cluster max mass for light quarks (unit [GeV])", &ClusterFissioner::_clMaxLight, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); static Parameter interfaceClMaxBottom ("ClMaxBottom","cluster max mass for b quarks (unit [GeV])", &ClusterFissioner::_clMaxBottom, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); static Parameter interfaceClMaxCharm ("ClMaxCharm","cluster max mass for c quarks (unit [GeV])", &ClusterFissioner::_clMaxCharm, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); static Parameter interfaceClMaxExotic ("ClMaxExotic","cluster max mass for exotic quarks (unit [GeV])", &ClusterFissioner::_clMaxExotic, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); // ClPow for light, Bottom, Charm and exotic (e.g. Susy) quarks static Parameter interfaceClPowLight ("ClPowLight","cluster mass exponent for light quarks", &ClusterFissioner::_clPowLight, 0, 2.0, 0.0, 10.0,false,false,false); static Parameter interfaceClPowBottom ("ClPowBottom","cluster mass exponent for b quarks", &ClusterFissioner::_clPowBottom, 0, 2.0, 0.0, 10.0,false,false,false); static Parameter interfaceClPowCharm ("ClPowCharm","cluster mass exponent for c quarks", &ClusterFissioner::_clPowCharm, 0, 2.0, 0.0, 10.0,false,false,false); static Parameter interfaceClPowExotic ("ClPowExotic","cluster mass exponent for exotic quarks", &ClusterFissioner::_clPowExotic, 0, 2.0, 0.0, 10.0,false,false,false); // PSplit for light, Bottom, Charm and exotic (e.g. Susy) quarks static Parameter interfacePSplitLight ("PSplitLight","cluster mass splitting param for light quarks", &ClusterFissioner::_pSplitLight, 0, 1.0, 0.0, 10.0,false,false,false); static Parameter interfacePSplitBottom ("PSplitBottom","cluster mass splitting param for b quarks", &ClusterFissioner::_pSplitBottom, 0, 1.0, 0.0, 10.0,false,false,false); static Parameter interfacePSplitCharm ("PSplitCharm","cluster mass splitting param for c quarks", &ClusterFissioner::_pSplitCharm, 0, 1.0, 0.0, 10.0,false,false,false); static Parameter interfacePSplitExotic ("PSplitExotic","cluster mass splitting param for exotic quarks", &ClusterFissioner::_pSplitExotic, 0, 1.0, 0.0, 10.0,false,false,false); static Switch interfaceFission ("Fission", "Option for different Fission options", &ClusterFissioner::_fissionCluster, 1, false, false); static SwitchOption interfaceFissionDefault (interfaceFission, "default", "Normal cluster fission which depends on the hadron selector class.", 0); static SwitchOption interfaceFissionNew (interfaceFission, "new", "Alternative cluster fission which does not depend on the hadron selector class", 1); static Parameter interfaceFissionPwtUquark ("FissionPwtUquark", "Weight for fission in U quarks", &ClusterFissioner::_fissionPwtUquark, 1, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceFissionPwtDquark ("FissionPwtDquark", "Weight for fission in D quarks", &ClusterFissioner::_fissionPwtDquark, 1, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceFissionPwtSquark ("FissionPwtSquark", "Weight for fission in S quarks", &ClusterFissioner::_fissionPwtSquark, 0.5, 0.0, 1.0, false, false, Interface::limited); static Switch interfaceRemnantOption ("RemnantOption", "Option for the treatment of remnant clusters", &ClusterFissioner::_iopRem, 1, false, false); static SwitchOption interfaceRemnantOptionSoft (interfaceRemnantOption, "Soft", "Both clusters produced in the fission of the beam cluster" " are treated as soft clusters.", 0); static SwitchOption interfaceRemnantOptionHard (interfaceRemnantOption, "Hard", "Only the cluster containing the remnant is treated as a soft cluster.", 1); static SwitchOption interfaceRemnantOptionVeryHard (interfaceRemnantOption, "VeryHard", "Even remnant clusters are treated as hard, i.e. all clusters the same", 2); static Parameter interfaceBTCLM ("SoftClusterFactor", "Parameter for the mass spectrum of remnant clusters", &ClusterFissioner::_btClM, GeV, 1.*GeV, 0.1*GeV, 10.0*GeV, false, false, Interface::limited); static Parameter interfaceStringTension ("StringTension", "String tension used in vertex displacement calculation", &ClusterFissioner::_kappa, GeV/meter, 1.0e15*GeV/meter, ZERO, ZERO, false, false, Interface::lowerlim); static Switch interfaceEnhanceSProb ("EnhanceSProb", "Option for enhancing strangeness", &ClusterFissioner::_enhanceSProb, 0, false, false); static SwitchOption interfaceEnhanceSProbNo (interfaceEnhanceSProb, "No", "No strangeness enhancement.", 0); static SwitchOption interfaceEnhanceSProbScaled (interfaceEnhanceSProb, "Scaled", "Scaled strangeness enhancement", 1); static SwitchOption interfaceEnhanceSProbExponential (interfaceEnhanceSProb, "Exponential", "Exponential strangeness enhancement", 2); static Switch interfaceMassMeasure ("MassMeasure", "Option to use different mass measures", &ClusterFissioner::_massMeasure,0,false,false); static SwitchOption interfaceMassMeasureMass (interfaceMassMeasure, "Mass", "Mass Measure", 0); static SwitchOption interfaceMassMeasureLambda (interfaceMassMeasure, "Lambda", "Lambda Measure", 1); static Parameter interfaceFissionMassScale ("FissionMassScale", "Cluster fission mass scale", &ClusterFissioner::_m0Fission, GeV, 2.0*GeV, 0.1*GeV, 50.*GeV, false, false, Interface::limited); } tPVector ClusterFissioner::fission(ClusterVector & clusters, bool softUEisOn) { // return if no clusters if (clusters.empty()) return tPVector(); /***************** * Loop over the (input) collection of cluster pointers, and store in * the vector splitClusters all the clusters that need to be split * (these are beam clusters, if soft underlying event is off, and * heavy non-beam clusters). ********************/ stack splitClusters; for(ClusterVector::iterator it = clusters.begin() ; it != clusters.end() ; ++it) { /************** * Skip 3-component clusters that have been redefined (as 2-component * clusters) or not available clusters. The latter check is indeed * redundant now, but it is used for possible future extensions in which, * for some reasons, some of the clusters found by ClusterFinder are tagged * straight away as not available. **************/ if((*it)->isRedefined() || !(*it)->isAvailable()) continue; // if the cluster is a beam cluster add it to the vector of clusters // to be split or if it is heavy if((*it)->isBeamCluster() || isHeavy(*it)) splitClusters.push(*it); } tPVector finalhadrons; cut(splitClusters, clusters, finalhadrons, softUEisOn); return finalhadrons; } void ClusterFissioner::cut(stack & clusterStack, ClusterVector &clusters, tPVector & finalhadrons, bool softUEisOn) { /************************************************** * This method does the splitting of the cluster pointed by cluPtr * and "recursively" by all of its cluster children, if heavy. All of these * new children clusters are added (indeed the pointers to them) to the * collection of cluster pointers collecCluPtr. The method works as follows. * Initially the vector vecCluPtr contains just the input pointer to the * cluster to be split. Then it will be filled "recursively" by all * of the cluster's children that are heavy enough to require, in their turn, * to be split. In each loop, the last element of the vector vecCluPtr is * considered (only once because it is then removed from the vector). * This approach is conceptually recursive, but avoid the overhead of * a concrete recursive function. Furthermore it requires minimal changes * in the case that the fission of an heavy cluster could produce more * than two cluster children as assumed now. * * Draw the masses: for normal, non-beam clusters a power-like mass dist * is used, whereas for beam clusters a fast-decreasing exponential mass * dist is used instead (to avoid many iterative splitting which could * produce an unphysical large transverse energy from a supposed soft beam * remnant process). ****************************************/ // Here we recursively loop over clusters in the stack and cut them while (!clusterStack.empty()) { // take the last element of the vector ClusterPtr iCluster = clusterStack.top(); clusterStack.pop(); // split it cutType ct = iCluster->numComponents() == 2 ? cutTwo(iCluster, finalhadrons, softUEisOn) : cutThree(iCluster, finalhadrons, softUEisOn); // There are cases when we don't want to split, even if it fails mass test if(!ct.first.first || !ct.second.first) { // if an unsplit beam cluster leave if for the underlying event if(iCluster->isBeamCluster() && softUEisOn) iCluster->isAvailable(false); continue; } // check if clusters ClusterPtr one = dynamic_ptr_cast(ct.first.first); ClusterPtr two = dynamic_ptr_cast(ct.second.first); // is a beam cluster must be split into two clusters if(iCluster->isBeamCluster() && (!one||!two) && softUEisOn) { iCluster->isAvailable(false); continue; } // There should always be a intermediate quark(s) from the splitting assert(ct.first.second && ct.second.second); /// \todo sort out motherless quark pairs here. Watch out for 'quark in final state' errors iCluster->addChild(ct.first.first); // iCluster->addChild(ct.first.second); // ct.first.second->addChild(ct.first.first); iCluster->addChild(ct.second.first); // iCluster->addChild(ct.second.second); // ct.second.second->addChild(ct.second.first); // Sometimes the clusters decay C -> H + C' rather then C -> C' + C'' if(one) { clusters.push_back(one); if(one->isBeamCluster() && softUEisOn) one->isAvailable(false); if(isHeavy(one) && one->isAvailable()) clusterStack.push(one); } if(two) { clusters.push_back(two); if(two->isBeamCluster() && softUEisOn) two->isAvailable(false); if(isHeavy(two) && two->isAvailable()) clusterStack.push(two); } } } namespace { /** * Check if can't make a hadron from the partons */ bool cantMakeHadron(tcPPtr p1, tcPPtr p2) { return ! CheckId::canBeHadron(p1->dataPtr(), p2->dataPtr()); } /** * Check if can't make a diquark from the partons */ bool cantMakeDiQuark(tcPPtr p1, tcPPtr p2) { long id1 = p1->id(), id2 = p2->id(); return ! (QuarkMatcher::Check(id1) && QuarkMatcher::Check(id2) && id1*id2>0); } } ClusterFissioner::cutType ClusterFissioner::cutTwo(ClusterPtr & cluster, tPVector & finalhadrons, bool softUEisOn) { // need to make sure only 2-cpt clusters get here assert(cluster->numComponents() == 2); tPPtr ptrQ1 = cluster->particle(0); tPPtr ptrQ2 = cluster->particle(1); Energy Mc = cluster->mass(); assert(ptrQ1); assert(ptrQ2); // And check if those particles are from a beam remnant bool rem1 = cluster->isBeamRemnant(0); bool rem2 = cluster->isBeamRemnant(1); // workout which distribution to use bool soft1(false),soft2(false); switch (_iopRem) { case 0: soft1 = rem1 || rem2; soft2 = rem2 || rem1; break; case 1: soft1 = rem1; soft2 = rem2; break; } // Initialization for the exponential ("soft") mass distribution. static const int max_loop = 1000; int counter = 0; Energy Mc1 = ZERO, Mc2 = ZERO,m1=ZERO,m2=ZERO,m=ZERO; tcPDPtr toHadron1, toHadron2; PPtr newPtr1 = PPtr (); PPtr newPtr2 = PPtr (); bool succeeded = false; 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()); + long rndSpin = UseRandom::rnd() > 0.5 ? 1 : 3; + diquark = CheckId::makeDiquark(ptrQ[iq2]->dataPtr(),newPtr2->dataPtr(),rndSpin); 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 = (_maxScale < scale) ? 0. : pow(_hadronsSelector->pwtSquark(),scale); else if (_enhanceSProb == 2) prob_s = (_maxScale < scale) ? 0. : exp(-scale); break; /* 1: ClusterFissioner uses its own unique set of weights, i.e. decoupled from ClusterDecayer */ case 1: prob_d = _fissionPwtDquark; prob_u = _fissionPwtUquark; if (_enhanceSProb == 1) prob_s = (_maxScale < scale) ? 0. : pow(_fissionPwtSquark,scale); else if (_enhanceSProb == 2) prob_s = (_maxScale < scale) ? 0. : exp(-scale); break; } int choice = UseRandom::rnd3(prob_u, prob_d, prob_s); long idNew = 0; switch (choice) { case 0: idNew = ThePEG::ParticleID::u; break; case 1: idNew = ThePEG::ParticleID::d; break; case 2: idNew = ThePEG::ParticleID::s; break; } newPtrPos = getParticle(idNew); newPtrNeg = getParticle(-idNew); assert (newPtrPos); assert(newPtrNeg); assert (newPtrPos->dataPtr()); assert(newPtrNeg->dataPtr()); } Energy2 ClusterFissioner::clustermass(const ClusterPtr & cluster){ 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/HadronSelector.cc b/Hadronization/HadronSelector.cc --- a/Hadronization/HadronSelector.cc +++ b/Hadronization/HadronSelector.cc @@ -1,1018 +1,1022 @@ // -*- C++ -*- // // HadronSelector.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the HadronSelector class. // #include "HadronSelector.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include #include #include #include #include "CheckId.h" #include using namespace Herwig; DescribeAbstractClass describeHadronSelector("Herwig::HadronSelector","Herwig.so"); namespace { // // debug helper // void dumpTable(const HadronSelector::HadronTable & tbl) { // typedef HadronSelector::HadronTable::const_iterator TableIter; // for (TableIter it = tbl.begin(); it != tbl.end(); ++it) { // cerr << it->first.first << ' ' // << it->first.second << '\n'; // for (HadronSelector::KupcoData::const_iterator jt = it->second.begin(); // jt != it->second.end(); ++jt) { // cerr << '\t' << *jt << '\n'; // } // } // } bool weightIsLess (pair a, pair b) { return a.second < b.second; } - - } ostream & operator<< (ostream & os, const HadronSelector::HadronInfo & hi ) { os << std::scientific << std::showpoint << std::setprecision(4) << setw(2) << hi.id << '\t' // << hi.ptrData << ' ' << hi.swtef << '\t' << hi.wt << '\t' << hi.overallWeight << '\t' << ounit(hi.mass,GeV); return os; } HadronSelector::HadronSelector(unsigned int opt) : _pwtDquark( 1.0 ),_pwtUquark( 1.0 ),_pwtSquark( 1.0 ),_pwtCquark( 0.0 ), _pwtBquark( 0.0 ),_pwtDIquarkS0( 1.0 ),_pwtDIquarkS1( 1.0 ), _weight1S0(Nmax,1.),_weight3S1(Nmax,1.),_weight1P1(Nmax,1.),_weight3P0(Nmax,1.), _weight3P1(Nmax,1.),_weight3P2(Nmax,1.),_weight1D2(Nmax,1.),_weight3D1(Nmax,1.), _weight3D2(Nmax,1.),_weight3D3(Nmax,1.), _repwt(Lmax,vector >(Jmax,vector(Nmax))), _sngWt( 1.0 ),_decWt( 1.0 ), _topt(opt),_trial(0), _limBottom(), _limCharm(), _limExotic(), belowThreshold_(0) { // The mixing angles // the ideal mixing angle const double idealAngleMix = atan( sqrt(0.5) ) * 180.0 / Constants::pi; // \eta-\eta' mixing angle _etamix = -23.0; // phi-omega mixing angle _phimix = +36.0; // h_1'-h_1 mixing angle _h1mix = idealAngleMix; // f_0(1710)-f_0(1370) mixing angle _f0mix = idealAngleMix; // f_1(1420)-f_1(1285)\f$ mixing angle _f1mix = idealAngleMix; // f'_2-f_2\f$ mixing angle _f2mix = +26.0; // eta_2(1870)-eta_2(1645) mixing angle _eta2mix = idealAngleMix; // phi(???)-omega(1650) mixing angle _omhmix = idealAngleMix; // phi_3-omega_3 mixing angle _ph3mix = +28.0; // eta(1475)-eta(1295) mixing angle _eta2Smix = idealAngleMix; // phi(1680)-omega(1420) mixing angle _phi2Smix = idealAngleMix; } void HadronSelector::persistentOutput(PersistentOStream & os) const { os << _partons << _pwtDquark << _pwtUquark << _pwtSquark << _pwtCquark << _pwtBquark << _pwtDIquarkS0 << _pwtDIquarkS1 << _etamix << _phimix << _h1mix << _f0mix << _f1mix << _f2mix << _eta2mix << _omhmix << _ph3mix << _eta2Smix << _phi2Smix << _weight1S0 << _weight3S1 << _weight1P1 << _weight3P0 << _weight3P1 << _weight3P2 << _weight1D2 << _weight3D1 << _weight3D2 << _weight3D3 << _forbidden << _sngWt << _decWt << _repwt << _pwt << _limBottom << _limCharm << _limExotic << belowThreshold_ << _table; } void HadronSelector::persistentInput(PersistentIStream & is, int) { is >> _partons >> _pwtDquark >> _pwtUquark >> _pwtSquark >> _pwtCquark >> _pwtBquark >> _pwtDIquarkS0 >> _pwtDIquarkS1 >> _etamix >> _phimix >> _h1mix >> _f0mix >> _f1mix >> _f2mix >> _eta2mix >> _omhmix >> _ph3mix >> _eta2Smix >> _phi2Smix >> _weight1S0 >> _weight3S1 >> _weight1P1 >> _weight3P0 >> _weight3P1 >> _weight3P2 >> _weight1D2 >> _weight3D1 >> _weight3D2 >> _weight3D3 >> _forbidden >> _sngWt >> _decWt >> _repwt >> _pwt >> _limBottom >> _limCharm >> _limExotic >> belowThreshold_ >> _table; } void HadronSelector::Init() { static ClassDocumentation documentation ("There is no documentation for the HadronSelector class"); static Parameter interfacePwtDquark("PwtDquark","Weight for choosing a quark D", &HadronSelector::_pwtDquark, 0, 1.0, 0.0, 10.0, false,false,false); static Parameter interfacePwtUquark("PwtUquark","Weight for choosing a quark U", &HadronSelector::_pwtUquark, 0, 1.0, 0.0, 10.0, false,false,false); static Parameter interfacePwtSquark("PwtSquark","Weight for choosing a quark S", &HadronSelector::_pwtSquark, 0, 1.0, 0.0, 10.0, false,false,false); static Parameter interfacePwtCquark("PwtCquark","Weight for choosing a quark C", &HadronSelector::_pwtCquark, 0, 0.0, 0.0, 10.0, false,false,false); static Parameter interfacePwtBquark("PwtBquark","Weight for choosing a quark B", &HadronSelector::_pwtBquark, 0, 0.0, 0.0, 10.0, false,false,false); static Parameter interfacePwtDIquarkS0("PwtDIquarkS0","Weight for choosing a spin-0 DIquark", &HadronSelector::_pwtDIquarkS0, 0, 1.0, 0.0, 100.0, false,false,false); static Parameter interfacePwtDIquarkS1("PwtDIquarkS1","Weight for choosing a spin-1 DIquark", &HadronSelector::_pwtDIquarkS1, 0, 1.0, 0.0, 100.0, false,false,false); static Parameter interfaceSngWt("SngWt","Weight for singlet baryons", &HadronSelector::_sngWt, 0, 1.0, 0.0, 10.0, false,false,false); static Parameter interfaceDecWt("DecWt","Weight for decuplet baryons", &HadronSelector::_decWt, 0, 1.0, 0.0, 10.0, false,false,false); static RefVector interfacePartons ("Partons", "The partons which are to be considered as the consistuents of the hadrons.", &HadronSelector::_partons, -1, false, false, true, false, false); static RefVector interfaceForbidden ("Forbidden", "The PDG codes of the particles which cannot be produced in the hadronization.", &HadronSelector::_forbidden, -1, false, false, true, false, false); // // mixing angles // // the ideal mixing angle const double idealAngleMix = atan( sqrt(0.5) ) * 180.0 / Constants::pi; static Parameter interface11S0Mixing ("11S0Mixing", "The mixing angle for the I=0 mesons from the 1 1S0 multiplet," " i.e. eta and etaprime.", &HadronSelector::_etamix, -23., -180., 180., false, false, Interface::limited); static Parameter interface13S1Mixing ("13S1Mixing", "The mixing angle for the I=0 mesons from the 1 3S1 multiplet," " i.e. phi and omega.", &HadronSelector::_phimix, +36., -180., 180., false, false, Interface::limited); static Parameter interface11P1Mixing ("11P1Mixing", "The mixing angle for the I=0 mesons from the 1 1P1 multiplet," " i.e. h_1' and h_1.", &HadronSelector::_h1mix, idealAngleMix, -180., 180., false, false, Interface::limited); static Parameter interface13P0Mixing ("13P0Mixing", "The mixing angle for the I=0 mesons from the 1 3P0 multiplet," " i.e. f_0(1710) and f_0(1370).", &HadronSelector::_f0mix, idealAngleMix, -180., 180., false, false, Interface::limited); static Parameter interface13P1Mixing ("13P1Mixing", "The mixing angle for the I=0 mesons from the 1 3P1 multiplet," " i.e. f_1(1420) and f_1(1285).", &HadronSelector::_f1mix, idealAngleMix, -180., 180., false, false, Interface::limited); static Parameter interface13P2Mixing ("13P2Mixing", "The mixing angle for the I=0 mesons from the 1 3P2 multiplet," " i.e. f'_2 and f_2.", &HadronSelector::_f2mix, 26.0, -180., 180., false, false, Interface::limited); static Parameter interface11D2Mixing ("11D2Mixing", "The mixing angle for the I=0 mesons from the 1 1D2 multiplet," " i.e. eta_2(1870) and eta_2(1645).", &HadronSelector::_eta2mix, idealAngleMix, -180., 180., false, false, Interface::limited); static Parameter interface13D0Mixing ("13D0Mixing", "The mixing angle for the I=0 mesons from the 1 3D0 multiplet," " i.e. eta_2(1870) phi(?) and omega(1650).", &HadronSelector::_omhmix, idealAngleMix, -180., 180., false, false, Interface::limited); static Parameter interface13D1Mixing ("13D1Mixing", "The mixing angle for the I=0 mesons from the 1 3D1 multiplet," " i.e. phi_3 and omega_3.", &HadronSelector::_ph3mix, 28.0, -180., 180., false, false, Interface::limited); static Parameter interface21S0Mixing ("21S0Mixing", "The mixing angle for the I=0 mesons from the 2 1S0 multiplet," " i.e. eta(1475) and eta(1295).", &HadronSelector::_eta2Smix, idealAngleMix, -180., 180., false, false, Interface::limited); static Parameter interface23S1Mixing ("23S1Mixing", "The mixing angle for the I=0 mesons from the 1 3S1 multiplet," " i.e. phi(1680) and omega(1420).", &HadronSelector::_phi2Smix, idealAngleMix, -180., 180., false, false, Interface::limited); // // the meson weights // static ParVector interface1S0Weights ("1S0Weights", "The weights for the 1S0 multiplets start with n=1.", &HadronSelector::_weight1S0, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface3S1Weights ("3S1Weights", "The weights for the 3S1 multiplets start with n=1.", &HadronSelector::_weight3S1, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface1P1Weights ("1P1Weights", "The weights for the 1P1 multiplets start with n=1.", &HadronSelector::_weight1P1, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface3P0Weights ("3P0Weights", "The weights for the 3P0 multiplets start with n=1.", &HadronSelector::_weight3P0, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface3P1Weights ("3P1Weights", "The weights for the 3P1 multiplets start with n=1.", &HadronSelector::_weight3P1, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface3P2Weights ("3P2Weights", "The weights for the 3P2 multiplets start with n=1.", &HadronSelector::_weight3P2, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface1D2Weights ("1D2Weights", "The weights for the 1D2 multiplets start with n=1.", &HadronSelector::_weight1D2, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface3D1Weights ("3D1Weights", "The weights for the 3D1 multiplets start with n=1.", &HadronSelector::_weight3D1, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface3D2Weights ("3D2Weights", "The weights for the 3D2 multiplets start with n=1.", &HadronSelector::_weight3D2, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface3D3Weights ("3D3Weights", "The weights for the 3D3 multiplets start with n=1.", &HadronSelector::_weight3D3, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static Switch interfaceTrial ("Trial", "A Debugging option to only produce certain types of hadrons", &HadronSelector::_trial, 0, false, false); static SwitchOption interfaceTrialAll (interfaceTrial, "All", "Produce all the hadrons", 0); static SwitchOption interfaceTrialPions (interfaceTrial, "Pions", "Only produce pions", 1); static SwitchOption interfaceTrialSpin2 (interfaceTrial, "Spin2", "Only mesons with spin less than or equal to two are produced", 2); static SwitchOption interfaceTrialSpin3 (interfaceTrial, "Spin3", "Only hadrons with spin less than or equal to three are produced", 3); static Parameter interfaceSingleHadronLimitBottom ("SingleHadronLimitBottom", "Threshold for one-hadron decay of b-cluster", &HadronSelector::_limBottom, 0, 0.0, 0.0, 100.0,false,false,false); static Parameter interfaceSingleHadronLimitCharm ("SingleHadronLimitCharm", "threshold for one-hadron decay of c-cluster", &HadronSelector::_limCharm, 0, 0.0, 0.0, 100.0,false,false,false); static Parameter interfaceSingleHadronLimitExotic ("SingleHadronLimitExotic", "threshold for one-hadron decay of exotic cluster", &HadronSelector::_limExotic, 0, 0.0, 0.0, 100.0,false,false,false); static Switch interfaceBelowThreshold ("BelowThreshold", "Option fo the selection of the hadrons if the cluster is below the pair threshold", &HadronSelector::belowThreshold_, 0, false, false); static SwitchOption interfaceBelowThresholdLightest (interfaceBelowThreshold, "Lightest", "Force cluster to decay to the lightest hadron with the appropriate flavours", 0); static SwitchOption interfaceBelowThresholdAll (interfaceBelowThreshold, "All", "Select from all the hadrons below the two hadron threshold according to their spin weights", 1); } double HadronSelector::mixingStateWeight(long id) const { switch(id) { case ParticleID::eta: return 0.5*probabilityMixing(_etamix ,1); case ParticleID::etaprime: return 0.5*probabilityMixing(_etamix ,2); case ParticleID::phi: return 0.5*probabilityMixing(_phimix ,1); case ParticleID::omega: return 0.5*probabilityMixing(_phimix ,2); case ParticleID::hprime_1: return 0.5*probabilityMixing(_h1mix ,1); case ParticleID::h_1: return 0.5*probabilityMixing(_h1mix ,2); case 10331: return 0.5*probabilityMixing(_f0mix ,1); case 10221: return 0.5*probabilityMixing(_f0mix ,2); case ParticleID::fprime_1: return 0.5*probabilityMixing(_f1mix ,1); case ParticleID::f_1: return 0.5*probabilityMixing(_f1mix ,2); case ParticleID::fprime_2: return 0.5*probabilityMixing(_f2mix ,1); case ParticleID::f_2: return 0.5*probabilityMixing(_f2mix ,2); case 10335: return 0.5*probabilityMixing(_eta2mix ,1); case 10225: return 0.5*probabilityMixing(_eta2mix ,2); // missing phi member of 13D1 should be here case 30223: return 0.5*probabilityMixing(_omhmix ,2); case 337: return 0.5*probabilityMixing(_ph3mix ,1); case 227: return 0.5*probabilityMixing(_ph3mix ,2); case 100331: return 0.5*probabilityMixing(_eta2mix ,1); case 100221: return 0.5*probabilityMixing(_eta2mix ,2); case 100333: return 0.5*probabilityMixing(_phi2Smix,1); case 100223: return 0.5*probabilityMixing(_phi2Smix,2); default: return 1./3.; } } void HadronSelector::doinit() { Interfaced::doinit(); // the default partons allowed // the quarks for ( int ix=1; ix<=5; ++ix ) { _partons.push_back(getParticleData(ix)); } // the diquarks for(unsigned int ix=1;ix<=5;++ix) { for(unsigned int iy=1; iy<=ix;++iy) { - _partons.push_back(getParticleData(CheckId::makeDiquarkID(ix,iy))); + _partons.push_back(getParticleData(CheckId::makeDiquarkID(ix,iy,long(3)))); + if(ix!=iy) + _partons.push_back(getParticleData(CheckId::makeDiquarkID(ix,iy,long(1)))); } } // set the weights for the various excited mesons // set all to one to start with for (int l = 0; l < Lmax; ++l ) { for (int j = 0; j < Jmax; ++j) { for (int n = 0; n < Nmax; ++n) { _repwt[l][j][n] = 1.0; } } } // set the others from the relevant vectors for( int ix=0;ixid()]=1.; } _pwt[1] = _pwtDquark; _pwt[2] = _pwtUquark; _pwt[3] = _pwtSquark; _pwt[4] = _pwtCquark; _pwt[5] = _pwtBquark; _pwt[1103] = _pwtDIquarkS1 * _pwtDquark * _pwtDquark; _pwt[2101] = 0.5 * _pwtDIquarkS0 * _pwtUquark * _pwtDquark; + _pwt[2103] = 0.5 * _pwtDIquarkS1 * _pwtUquark * _pwtDquark; _pwt[2203] = _pwtDIquarkS1 * _pwtUquark * _pwtUquark; _pwt[3101] = 0.5 * _pwtDIquarkS0 * _pwtSquark * _pwtDquark; + _pwt[3103] = 0.5 * _pwtDIquarkS1 * _pwtSquark * _pwtDquark; _pwt[3201] = 0.5 * _pwtDIquarkS0 * _pwtSquark * _pwtUquark; + _pwt[3203] = 0.5 * _pwtDIquarkS1 * _pwtSquark * _pwtUquark; _pwt[3303] = _pwtDIquarkS1 * _pwtSquark * _pwtSquark; + // Commenting out heavy di-quark weights _pwt[4101] = 0.0; _pwt[4201] = 0.0; _pwt[4301] = 0.0; _pwt[4403] = 0.0; _pwt[5101] = 0.0; _pwt[5201] = 0.0; _pwt[5301] = 0.0; _pwt[5401] = 0.0; _pwt[5503] = 0.0; // find the maximum map::iterator pit = max_element(_pwt.begin(),_pwt.end(),weightIsLess); const double pmax = pit->second; for(pit=_pwt.begin(); pit!=_pwt.end(); ++pit) { pit->second/=pmax; } // construct the hadron tables constructHadronTable(); - // for debugging // dumpTable(table()); } void HadronSelector::constructHadronTable() { // initialise the table _table.clear(); for(unsigned int ix=0; ix<_partons.size(); ++ix) { for(unsigned int iy=0; iy<_partons.size(); ++iy) { if (!(DiquarkMatcher::Check(_partons[ix]->id()) && DiquarkMatcher::Check(_partons[iy]->id()))) _table[make_pair(_partons[ix]->id(),_partons[iy]->id())] = KupcoData(); } } // get the particles from the event generator ParticleMap particles = generator()->particles(); // loop over the particles double maxdd(0.),maxss(0.),maxrest(0.); for(ParticleMap::iterator it=particles.begin(); it!=particles.end(); ++it) { long pid = it->first; tPDPtr particle = it->second; int pspin = particle->iSpin(); // Don't include hadrons which are explicitly forbidden if(find(_forbidden.begin(),_forbidden.end(),particle)!=_forbidden.end()) continue; // Don't include non-hadrons or antiparticles if(pid < 100) continue; // remove diffractive particles if(pspin == 0) continue; // K_0S and K_0L not made make K0 and Kbar0 if(pid==ParticleID::K_S0||pid==ParticleID::K_L0) continue; // Debugging options // Only include those with 2J+1 less than...5 if(_trial==2 && pspin >= 5) continue; // Only include those with 2J+1 less than...7 if(_trial==3 && pspin >= 7) continue; // Only include pions if(_trial==1 && pid!=111 && pid!=211) continue; // shouldn't be coloured if(particle->coloured()) continue; // Get the flavours const int x4 = (pid/1000)%10; const int x3 = (pid/100 )%10; const int x2 = (pid/10 )%10; const int x7 = (pid/1000000)%10; const bool wantSusy = x7 == 1 || x7 == 2; int flav1; int flav2; // Skip non-hadrons (susy particles, etc...) if(x3 == 0 || x2 == 0) continue; else if(x4 == 0) { // meson flav1 = x2; flav2 = x3; } else { // baryon - flav1 = CheckId::makeDiquarkID(x2,x3); + long rndSpin = UseRandom::rnd() > 0.5 ? 1 : 3; + flav1 = CheckId::makeDiquarkID(x2,x3,rndSpin); flav2 = x4; } if (wantSusy) flav2 += 1000000 * x7; HadronInfo a(pid, particle, specialWeight(pid), particle->mass()); // set the weight to the number of spin states a.overallWeight = pspin; // identical light flavours if(flav1 == flav2 && flav1<=3) { // ddbar> uubar> admixture states if(flav1==1) { if(_topt != 0) a.overallWeight *= 0.5*a.swtef; _table[make_pair(1,1)].insert(a); _table[make_pair(2,2)].insert(a); if(_topt == 0 && a.overallWeight > maxdd) maxdd = a.overallWeight; } // load up ssbar> uubar> ddbar> admixture states else { a.wt = mixingStateWeight(pid); a.overallWeight *= a.wt; if(_topt != 0) a.overallWeight *= a.swtef; _table[make_pair(1,1)].insert(a); _table[make_pair(2,2)].insert(a); if(_topt == 0 && a.overallWeight > maxdd) maxdd = a.overallWeight; a.wt = (_topt != 0) ? 1.- 2.*a.wt : 1 - a.wt; if(a.wt > 0) { a.overallWeight = a.wt * a.swtef * pspin; _table[make_pair(3,3)].insert(a); if(_topt == 0 && a.overallWeight > maxss) maxss = a.overallWeight; } } } // light baryons with all quarks identical else if((flav1 == 1 && flav2 == 1103) || (flav1 == 1103 && flav2 == 1) || (flav1 == 2 && flav2 == 2203) || (flav1 == 2203 && flav2 == 2) || (flav1 == 3 && flav2 == 3303) || (flav1 == 3303 && flav2 == 3)) { if(_topt != 0) a.overallWeight *= 1.5*a.swtef; _table[make_pair(flav1,flav2)].insert(a); _table[make_pair(flav2,flav1)].insert(a); if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight; } // all other cases else { if(_topt != 0) a.overallWeight *=a.swtef; _table[make_pair(flav1,flav2)].insert(a); if(flav1 != flav2) _table[make_pair(flav2,flav1)].insert(a); if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight; } } // Account for identical combos of diquark/quarks and symmetrical elements // e.g. U UD = D UU HadronTable::iterator tit; for(tit=_table.begin();tit!=_table.end();++tit) { if(tit->first.first>ParticleID::c) continue; if(!DiquarkMatcher::Check(tit->first.second)) continue; long k, l, sub; if(tit->first.second>=ParticleID::bd_0) { k = ParticleID::b; sub = ParticleID::bd_0/100; } else if(tit->first.second>=ParticleID::cd_0) { k = ParticleID::c; sub = ParticleID::cd_0/100; } else if(tit->first.second>=ParticleID::sd_0) { k = ParticleID::s; sub = ParticleID::sd_0/100; } else if(tit->first.second>=ParticleID::ud_0) { k = ParticleID::u; sub = ParticleID::ud_0/100; } else if(tit->first.second==ParticleID::dd_1) { k = ParticleID::d; sub = ParticleID::dd_1/100; } else continue; sub=tit->first.second/100-sub+1; if(sub > tit->first.first) { l = 1000*sub+100*tit->first.first+1; } else if(sub==tit->first.first) { l = 1000*sub+ 100*tit->first.first+3; } else { l = 100*sub +1000*tit->first.first+1; } if(tit->second.empty()) { pair newpair(k,l); tit->second=_table[newpair]; newpair=make_pair(tit->first.second,tit->first.first); _table[newpair]=tit->second; }; } // normalise weights to one for first option if(_topt == 0) { HadronTable::const_iterator tit; KupcoData::iterator it; for(tit=_table.begin();tit!=_table.end();++tit) { double weight; if(tit->first.first==tit->first.second) { if(tit->first.first==1||tit->first.first==2) weight=1./maxdd; else if (tit->first.first==3) weight=1./maxss; else weight=1./maxrest; } else weight=1./maxrest; for(it = tit->second.begin(); it!=tit->second.end(); ++it) { it->rescale(weight); } } } } double HadronSelector::specialWeight(long id) const { const int pspin = id % 10; // Only K0L and K0S have pspin == 0, should // not get them until Decay step assert( pspin != 0 ); // Baryon : J = 1/2 or 3/2 if(pspin == 2) { // Singlet (Lambda-like) baryon if( (id/100)%10 < (id/10 )%10 ) return sqr(_sngWt); // octet else return 1.; } // Decuplet baryon else if (pspin == 4) { return sqr(_decWt); } // Meson else if(pspin % 2 == 1) { // Total angular momentum int j = (pspin - 1) / 2; // related to Orbital angular momentum l int nl = (id/10000 )%10; int l = -999; int n = (id/100000)%10; // Radial excitation if(j == 0) l = nl; else if(nl == 0) l = j - 1; else if(nl == 1 || nl == 2) l = j; else if(nl == 3) l = j + 1; // Angular or Radial excited meson if((l||j||n) && l>=0 && l= 5/2 (ispin >= 6), haven't got those return 1.0; } int HadronSelector::signHadron(tcPDPtr idQ1, tcPDPtr idQ2, tcPDPtr hadron) const { // This method receives in input three PDG ids, whose the // first two have proper signs (corresponding to particles, id > 0, // or antiparticles, id < 0 ), whereas the third one must // be always positive (particle not antiparticle), // corresponding to: // --- quark-antiquark, or antiquark-quark, or // quark-diquark, or diquark-quark, or // antiquark-antidiquark, or antidiquark-antiquark // for the first two input (idQ1, idQ2); // --- meson or baryon for the third input (idHad): // The method returns: // --- + 1 if the two partons (idQ1, idQ2) are exactly // the constituents for the hadron idHad; // --- - 1 if the two partons (idQ1, idQ2) are exactly // the constituents for the anti-hadron -idHad; // --- + 0 otherwise. // The method it is therefore useful to decide the // sign of the id of the produced hadron as appeared // in the vector _vecHad (where only hadron idHad > 0 are present) // given the two constituent partons. int sign = 0; long idHad = hadron->id(); assert(idHad > 0); int chargeIn = idQ1->iCharge() + idQ2->iCharge(); int chargeOut = hadron->iCharge(); // same charge if( chargeIn == chargeOut && chargeIn !=0 ) sign = +1; else if(chargeIn == -chargeOut && chargeIn !=0 ) sign = -1; else if(chargeIn == 0 && chargeOut == 0 ) { // In the case of same null charge, there are four cases: // i) K0-like mesons, B0-like mesons, Bs-like mesons // the PDG convention is to consider them "antiparticle" (idHad < 0) // if the "dominant" (heavier) flavour (respectively, s, b) // is a quark (idQ > 0): for instance, B0s = (b, sbar) has id < 0 // Remember that there is an important exception for K0L (id=130) and // K0S (id=310): they don't have antiparticles, therefore idHad > 0 // always. We use below the fact that K0L and K0S are the unique // hadrons having 0 the first (less significant) digit of their id. // 2) D0-like mesons: the PDG convention is to consider them "particle" // (idHad > 0) if the charm flavour is carried by a c: (c,ubar) has id>0 // 3) the remaining mesons should not have antiparticle, therefore their // sign is always positive. // 4) for baryons, that is when one of idQ1 and idQ2 is a (anti-) quark and // the other one is a (anti-) diquark the sign is negative when both // constituents are "anti", that is both with id < 0; positive otherwise. // meson if(abs(int(idQ1->iColour()))== 3 && abs(int(idQ2->iColour())) == 3 && !DiquarkMatcher::Check(idQ1->id()) && !DiquarkMatcher::Check(idQ2->id())) { int idQa = abs(idQ1->id()); int idQb = abs(idQ2->id()); int dominant = idQ2->id(); if(idQa > idQb) { swap(idQa,idQb); dominant = idQ1->id(); } if((idQa==ParticleID::d && idQb==ParticleID::s) || (idQa==ParticleID::d && idQb==ParticleID::b) || (idQa==ParticleID::s && idQb==ParticleID::b)) { // idHad%10 is zero for K0L,K0S if (dominant < 0 || idHad%10 == 0) sign = +1; else if(dominant > 0) sign = -1; } else if((idQa==ParticleID::u && idQb==ParticleID::c) || (idQa==ParticleID::u && idQb==ParticleID::t) || (idQa==ParticleID::c && idQb==ParticleID::t)) { if (dominant > 0) sign = +1; else if(dominant < 0) sign = -1; } else if(idQa==idQb) sign = +1; // sets sign for Susy particles else sign = (dominant > 0) ? +1 : -1; } // baryon else if(DiquarkMatcher::Check(idQ1->id()) || DiquarkMatcher::Check(idQ2->id())) { if (idQ1->id() > 0 && idQ2->id() > 0) sign = +1; else if(idQ1->id() < 0 && idQ2->id() < 0) sign = -1; } } if (sign == 0) { cerr << "Could not work out sign for " << idQ1->PDGName() << ' ' << idQ2->PDGName() << " => " << hadron->PDGName() << '\n'; assert(false); } return sign; } pair HadronSelector::lightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr ptr3) const { // throw exception of id3!=0 as doesn't work if ( ptr3 ) throw Exception() << "ptr3!=0 not yet implemented in HadronSelector::lightestHadronPair" << Exception::abortnow; // charge int totalcharge = ptr1->iCharge() + ptr2->iCharge(); if ( ptr3 ) totalcharge += ptr3->iCharge(); tcPDPtr vIdHad1[2]={tcPDPtr(),tcPDPtr()},vIdHad2[2]={tcPDPtr(),tcPDPtr()}; bool vOk[2] = {false, false}; Energy vMassPair[2] = { ZERO, ZERO }; for (int i = 0; i < 2; i++) { tcPDPtr idPartner = i==0 ? getParticleData(ParticleID::d) : getParticleData(ParticleID::u); // Change sign to idPartner (transform it into a anti-quark) if it is not // possible to form a meson or a baryon. assert (ptr1 && idPartner); if (!CheckId::canBeHadron(ptr1, idPartner)) idPartner = idPartner->CC(); vIdHad1[i] = lightestHadron(ptr1, idPartner); vIdHad2[i] = lightestHadron(ptr2, idPartner->CC()); if ( vIdHad1[i] && vIdHad2[i] && vIdHad1[i]->iCharge() + vIdHad2[i]->iCharge() == totalcharge ) { vOk[i] = true; vMassPair[i] = vIdHad1[i]->mass() + vIdHad2[i]->mass(); } } // Take the lightest pair compatible with charge conservation. if ( vOk[0] && ( ! vOk[1] || vMassPair[0] <= vMassPair[1] ) ) { return make_pair(vIdHad1[0],vIdHad2[0]); } else if ( vOk[1] && ( ! vOk[0] || vMassPair[1] < vMassPair[0] ) ) { return make_pair(vIdHad1[1],vIdHad2[1]); } else { return make_pair(tcPDPtr(),tcPDPtr()); } } Energy HadronSelector::massLightestBaryonPair(tcPDPtr ptr1, tcPDPtr ptr2) const { // Make sure that we don't have any diquarks as input, return arbitrarily // large value if we do Energy currentSum = Constants::MaxEnergy; for(unsigned int ix=0; ix<_partons.size(); ++ix) { if(!DiquarkMatcher::Check(_partons[ix]->id())) continue; HadronTable::const_iterator tit1=_table.find(make_pair(abs(ptr1->id()),_partons[ix]->id())), tit2=_table.find(make_pair(_partons[ix]->id(),abs(ptr2->id()))); if( tit1==_table.end() || tit2==_table.end()) continue; if(tit1->second.empty()||tit2->second.empty()) continue; Energy s = tit1->second.begin()->mass + tit2->second.begin()->mass; if(currentSum > s) currentSum = s; } return currentSum; } tcPDPtr HadronSelector::lightestHadron(tcPDPtr ptr1, tcPDPtr ptr2, #ifndef NDEBUG tcPDPtr ptr3) const { #else tcPDPtr ) const { #endif // The method assumes ptr3 == 0 rest not implemented assert(ptr1 && ptr2 && !ptr3); // find entry in the table pair ids = make_pair(abs(ptr1->id()),abs(ptr2->id())); HadronTable::const_iterator tit=_table.find(ids); // throw exception if flavours wrong if (tit==_table.end()) throw Exception() << "Could not find " << ids.first << ' ' << ids.second << " in _table. " << "In HadronSelector::lightestHadron()" << Exception::eventerror; if(tit->second.empty()) throw Exception() << "HadronSelector::lightestHadron " << "could not find any hadrons containing " << ptr1->id() << ' ' << ptr2->id() << '\n' << tit->first.first << ' ' << tit->first.second << Exception::eventerror; // find the lightest hadron int sign = signHadron(ptr1,ptr2,tit->second.begin()->ptrData); tcPDPtr candidate = sign > 0 ? tit->second.begin()->ptrData : tit->second.begin()->ptrData->CC(); // \todo 20 GeV limit is temporary fudge to let SM particles go through. // \todo Use isExotic instead? if (candidate->mass() > 20*GeV && candidate->mass() < ptr1->constituentMass() + ptr2->constituentMass()) { generator()->log() << "HadronSelector::lightestHadron: " << "chosen candidate " << candidate->PDGName() << " is lighter than its constituents " << ptr1->PDGName() << ", " << ptr2->PDGName() << '\n' << candidate->mass()/GeV << " < " << ptr1->constituentMass()/GeV << " + " << ptr2->constituentMass()/GeV << '\n' << "Check your particle data tables.\n"; assert(false); } return candidate; } vector > HadronSelector::hadronsBelowThreshold(Energy threshold, tcPDPtr ptr1, tcPDPtr ptr2, #ifndef NDEBUG tcPDPtr ptr3) const { #else tcPDPtr ) const { #endif // The method assumes ptr3 == 0 rest not implemented assert(ptr1 && ptr2 && !ptr3); // find entry in the table pair ids = make_pair(abs(ptr1->id()),abs(ptr2->id())); HadronTable::const_iterator tit=_table.find(ids); // throw exception if flavours wrong if (tit==_table.end()) throw Exception() << "Could not find " << ids.first << ' ' << ids.second << " in _table. " << "In HadronSelector::hadronsBelowThreshold()" << Exception::eventerror; if(tit->second.empty()) throw Exception() << "HadronSelector::hadronsBelowThreshold() " << "could not find any hadrons containing " << ptr1->id() << ' ' << ptr2->id() << '\n' << tit->first.first << ' ' << tit->first.second << Exception::eventerror; vector > candidates; KupcoData::const_iterator hit = tit->second.begin(); // find the hadrons while(hit!=tit->second.end()&&hit->massptrData); tcPDPtr candidate = sign > 0 ? hit->ptrData : hit->ptrData->CC(); // \todo 20 GeV limit is temporary fudge to let SM particles go through. // \todo Use isExotic instead? if (candidate->mass() > 20*GeV && candidate->mass() < ptr1->constituentMass() + ptr2->constituentMass()) { generator()->log() << "HadronSelector::hadronsBelowTheshold: " << "chosen candidate " << candidate->PDGName() << " is lighter than its constituents " << ptr1->PDGName() << ", " << ptr2->PDGName() << '\n' << candidate->mass()/GeV << " < " << ptr1->constituentMass()/GeV << " + " << ptr2->constituentMass()/GeV << '\n' << "Check your particle data tables.\n"; assert(false); } candidates.push_back(make_pair(candidate,hit->overallWeight)); ++hit; } return candidates; } tcPDPtr HadronSelector::chooseSingleHadron(tcPDPtr par1, tcPDPtr par2, Energy mass) const { // Determine the sum of the nominal masses of the two lightest hadrons // with the right flavour numbers as the cluster under consideration. // Notice that we don't need real masses (drawn by a Breit-Wigner // distribution) because the lightest pair of hadrons does not involve // any broad resonance. Energy threshold = massLightestHadronPair(par1,par2); // Special: it allows one-hadron decays also above threshold. if (CheckId::isExotic(par1,par2)) threshold *= (1.0 + UseRandom::rnd()*_limExotic); else if (CheckId::hasBottom(par1,par2)) threshold *= (1.0 + UseRandom::rnd()*_limBottom); else if (CheckId::hasCharm(par1,par2)) threshold *= (1.0 + UseRandom::rnd()*_limCharm); // only do one hadron decay is mass less than the threshold if(mass>=threshold) return tcPDPtr(); // select the hadron tcPDPtr hadron; // old option pick the lightest hadron if(belowThreshold_ == 0) { hadron= lightestHadron(par1,par2); } // new option select from those available else if(belowThreshold_ == 1) { vector > hadrons = hadronsBelowThreshold(threshold,par1,par2); if(hadrons.size()==1) { hadron = hadrons[0].first; } else if(hadrons.empty()) { hadron= lightestHadron(par1,par2); } else { double totalWeight=0.; for(unsigned int ix=0;ix #include using namespace Herwig; DescribeClass describeHwppSelector("Herwig::HwppSelector","Herwig.so"); 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 << _scHadronWtFactor << _sbHadronWtFactor; } void HwppSelector::persistentInput(PersistentIStream & is, int) { is >> _mode >> _enhanceSProb >> iunit(_m0Decay,GeV) >> _massMeasure >> _scHadronWtFactor >> _sbHadronWtFactor; } 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 interfacescHadronWtFactor ("scHadronWtFactor", "Wight factor for strenge-charm heavy hadrns", &HwppSelector::_scHadronWtFactor, 1., 0., 10., false, false, Interface::limited); static Parameter interfacesbHadronWtFactor ("sbHadronWtFactor", "Wight factor for strenge-bottom heavy hadrns", &HwppSelector::_sbHadronWtFactor, 1., 0., 10., false, false, Interface::limited); 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.+pwtDIquarkS0()+pwtDIquarkS1()) && 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; // quark weight double quarkWeight = pwt(quarktopick->id()); if(abs(quarktopick->id()) == 3) { // Decoupling the weight of heavy strenge hadrons if(_enhanceSProb == 0 && abs(par1->id()) == 4) { quarkWeight = pwt(quarktopick->id())*_scHadronWtFactor; } else if(_enhanceSProb == 0 && abs(par1->id()) == 5) { quarkWeight = pwt(quarktopick->id())*_sbHadronWtFactor; } // Scaling strangeness enhancement else if(_enhanceSProb == 1) { double scale = double(sqr(_m0Decay/cluMass)); quarkWeight = (_maxScale < scale) ? 0. : pow(quarkWeight,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); quarkWeight = (_maxScale < scale) ? 0. : exp(-scale); } } // 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; weight = quarkWeight * H1->overallWeight * H2->overallWeight * Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass); - + //cerr<id<<" "<id<<" "<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())); }