diff --git a/Hadronization/HadronSelector.cc b/Hadronization/HadronSelector.cc --- a/Hadronization/HadronSelector.cc +++ b/Hadronization/HadronSelector.cc @@ -1,1013 +1,998 @@ // -*- 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",""); 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'; - // } - // } - // } + // 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 ),_pwtDIquark( 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 << _pwtDIquark << _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 >> _pwtDIquark>> _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 interfacePwtDIquark("PwtDIquark","Weight for choosing a DIquark", &HadronSelector::_pwtDIquark, 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))); } } // 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] = _pwtDIquark * _pwtDquark * _pwtDquark; _pwt[2101] = 0.5 * _pwtDIquark * _pwtUquark * _pwtDquark; _pwt[2203] = _pwtDIquark * _pwtUquark * _pwtUquark; _pwt[3101] = 0.5 * _pwtDIquark * _pwtSquark * _pwtDquark; _pwt[3201] = 0.5 * _pwtDIquark * _pwtSquark * _pwtUquark; _pwt[3303] = _pwtDIquark * _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); 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 #include #include #include #include "HadronSelector.fh" namespace Herwig { using namespace ThePEG; /**\ingroup Hadronization * \class HadronSelector * \brief This class selects the hadron flavours of a cluster decay. * \author Philip Stephens * \author Alberto Ribon * \author Peter Richardson * * This is the base class for the selection of either a pair of hadrons, or * in some cases a single hadron. The different approaches which were * previously implemented in this class are now implemented in the * HwppSelector and Hw64Selector which inherit from this class. * * This class implements a number of methods which are needed by all models * and in addition contains the weights for the different meson multiplets and * mixing of the light \f$I=0\f$ mesons. * * @see \ref HadronSelectorInterfaces "The interfaces" * defined for HadronSelector. * @see HwppSelector * @see Hw64Selector */ class HadronSelector: public Interfaced { public: /** \ingroup Hadronization * \class HadronInfo * \brief Class used to store all the hadron information for easy access. * \author Philip Stephens * * Note that: * - the hadrons in _table can be filled in any ordered * w.r.t. the mass value, and flavours for different * groups (for instance, (u,s) hadrons don't need to * be placed after (d,s) or any other flavour), but * all hadrons with the same flavours must be consecutive * ( for instance you cannot alternate hadrons of type * (d,s) with those of flavour (u,s) ). * Furthermore, it is assumed that particle and antiparticle * have the same weights, and therefore only one of them * must be entered in the table: we have chosen to refer * to the particle, defined as PDG id > 0, although if * an anti-particle is provided in input it is automatically * transform to its particle, simply by taking the modulus * of its id. */ class HadronInfo { public: /** * Constructor * @param idin The PDG code of the hadron * @param datain The pointer to the ParticleData object * @param swtin The singlet/decuplet/orbital factor * @param massin The mass of the hadron */ HadronInfo(long idin=0, tPDPtr datain=tPDPtr(), double swtin=1., Energy massin=ZERO) : id(idin), ptrData(datain), swtef(swtin), wt(1.0), overallWeight(0.0), mass(massin) {} /** * Comparision operator on mass */ bool operator<(const HadronInfo &x) const { if(mass!=x.mass) return mass < x.mass; else return id < x.id; } /** * The hadrons id. */ long id; /** * pointer to ParticleData, to get the spin, etc... */ tPDPtr ptrData; /** * singlet/decuplet/orbital factor */ double swtef; /** * mixing factor */ double wt; /** * (2*J+1)*wt*swtef */ double overallWeight; /** * The hadrons mass */ Energy mass; /** * Rescale the weight for a given hadron */ void rescale(double x) const { const_cast(this)->overallWeight *= x; } /** * Friend method used to print the value of a table element. */ friend PersistentOStream & operator<< (PersistentOStream & os, const HadronInfo & hi ) { os << hi.id << hi.ptrData << hi.swtef << hi.wt << hi.overallWeight << ounit(hi.mass,GeV); return os; } /** * debug output */ - friend ostream & operator<< (ostream & os, const HadronInfo & hi ); + friend ostream & operator<< (ostream & os, const HadronInfo & hi ) { + os << std::scientific << std::showpoint + << std::setprecision(4) + << setw(2) + << hi.id << '\t' + << hi.swtef << '\t' + << hi.wt << '\t' + << hi.overallWeight << '\t' + << ounit(hi.mass,GeV); + return os; + } /** * Friend method used to read in the value of a table element. */ friend PersistentIStream & operator>> (PersistentIStream & is, HadronInfo & hi ) { is >> hi.id >> hi.ptrData >> hi.swtef >> hi.wt >> hi.overallWeight >> iunit(hi.mass,GeV); return is; } }; /** \ingroup Hadronization * \class Kupco * \brief Class designed to make STL routines easy to use. * \author Philip Stephens * * This class is used to generate a list of the hadron pairs which can * be produced that allows easy traversal and quick access. */ class Kupco { public: /** * Constructor * @param inidQ PDG code of the quark drawn from the vacuum. * @param inhad1 ParticleData for the first hadron produced. * @param inhad2 ParticleData for the second hadron produced. * @param inwgt The weight for the hadron pair */ Kupco(tcPDPtr inidQ,tcPDPtr inhad1,tcPDPtr inhad2, Energy inwgt) : idQ(inidQ),hadron1(inhad1),hadron2(inhad2),weight(inwgt) {} /** * id of the quark drawn from the vacuum. */ tcPDPtr idQ; /** * The ParticleData object for the first hadron produced. */ tcPDPtr hadron1; /** * The ParticleData object for the second hadron produced. */ tcPDPtr hadron2; /** * Weight factor of this componation. */ Energy weight; }; public: /** * The helper classes */ //@{ /** * The type is used to contain all the hadrons info of a given flavour. */ typedef set KupcoData; //@} /** * The hadron table type. */ typedef map,KupcoData> HadronTable; public: /** * The default constructor. */ HadronSelector(unsigned int); /** * Method to return a pair of hadrons given the PDG codes of * two or three constituents * @param cluMass The mass of the cluster * @param par1 The first constituent * @param par2 The second constituent * @param par3 The third constituent */ virtual pair chooseHadronPair(const Energy cluMass, tcPDPtr par1, tcPDPtr par2,tcPDPtr par3 = PDPtr()) const = 0; /** * Select the single hadron for a cluster decay * return null pointer if not a single hadron decay * @param par1 1st constituent * @param par2 2nd constituent * @param mass Mass of the cluster */ tcPDPtr chooseSingleHadron(tcPDPtr par1, tcPDPtr par2, Energy mass) const; /** * This returns the lightest pair of hadrons given by the flavours. * * Given the two (or three) constituents of a cluster, it returns * the two lightest hadrons with proper flavour numbers. * Furthermore, the first of the two hadrons must have the constituent with * par1, and the second must have the constituent with par2. * \todo At the moment it does *nothing* in the case that also par3 is present. * * The method is implemented by calling twice lightestHadron, * once with (par1,quarktopick->CC()) ,and once with (par2,quarktopick) * where quarktopick is either the pointer to * d or u quarks . In fact, the idea is that whatever the flavour of par1 * and par2, no matter if (anti-)quark or (anti-)diquark, the lightest * pair of hadrons containing flavour par1 and par2 will have either * flavour d or u, being the lightest quarks. * The method returns the pair (PDPtr(),PDPtr()) if anything goes wrong. * * \todo The method assumes par3 == PDPtr() (otherwise we don't know how to proceed: a * possible, trivial way would be to randomly select two of the three * (anti-)quarks and treat them as a (anti-)diquark, reducing the problem * to two components as treated below. * In the normal (two components) situation, the strategy is the following: * treat in the same way the two possibilities: (d dbar) (i=0) and * (u ubar) (i=1) as the pair quark-antiquark necessary to form a * pair of hadrons containing the input flavour par1 and par2; finally, * select the one that produces the lightest pair of hadrons, compatible * with the charge conservation constraint. */ pair lightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr ptr3 = PDPtr ()) const; /** * Returns the mass of the lightest pair of hadrons with the given particles * @param ptr1 is the first constituent * @param ptr2 is the second constituent * @param ptr3 is the third constituent */ Energy massLightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr ptr3 = PDPtr ()) const { pair pairData = lightestHadronPair(ptr1, ptr2, ptr3); if ( ! pairData.first || ! pairData.second ) return ZERO; return ( pairData.first->mass() + pairData.second->mass() ); } /** * Returns the lightest hadron formed by the given particles. * * Given the id of two (or three) constituents of a cluster, it returns * the lightest hadron with proper flavour numbers. * At the moment it does *nothing* in the case that also 'ptr3' present. * @param ptr1 is the first constituent * @param ptr2 is the second constituent * @param ptr3 is the third constituent */ tcPDPtr lightestHadron(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr ptr3 = PDPtr ()) const; /** * Returns the hadrons below the constituent mass threshold formed by the given particles, * together with their total weight * * Given the id of two (or three) constituents of a cluster, it returns * the lightest hadron with proper flavour numbers. * At the moment it does *nothing* in the case that also 'ptr3' present. * @param threshold The theshold * @param ptr1 is the first constituent * @param ptr2 is the second constituent * @param ptr3 is the third constituent */ vector > hadronsBelowThreshold(Energy threshold, tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr ptr3 = PDPtr ()) const; /** * Return the nominal mass of the hadron returned by lightestHadron() * @param ptr1 is the first constituent * @param ptr2 is the second constituent * @param ptr3 is the third constituent */ Energy massLightestHadron(tcPDPtr ptr1, tcPDPtr ptr2, #ifndef NDEBUG tcPDPtr ptr3 = PDPtr ()) const { #else tcPDPtr = PDPtr ()) const { #endif // The method assumes ptr3 == empty assert(!(ptr3)); // find entry in the table pair ids(abs(ptr1->id()),abs(ptr2->id())); HadronTable::const_iterator tit=_table.find(ids); // throw exception if flavours wrong if(tit==_table.end()||tit->second.empty()) throw Exception() << "HadronSelector::massLightestHadron " << "failed for particle" << ptr1->id() << " " << ptr2->id() << Exception::eventerror; // return the mass return tit->second.begin()->mass; } /** * Returns the mass of the lightest pair of baryons. * @param ptr1 is the first constituent * @param ptr2 is the second constituent */ Energy massLightestBaryonPair(tcPDPtr ptr1, tcPDPtr ptr2) const; /** * Return the weights for the different quarks and diquarks */ //@{ /** * The down quark weight. */ double pwtDquark() const { return _pwtDquark; } /** * The up quark weight. */ double pwtUquark() const { return _pwtUquark; } /** * The strange quark weight. */ double pwtSquark() const { return _pwtSquark; } /** * The charm quark weight. */ double pwtCquark() const { return _pwtCquark; } /** * The bottom quark weight. */ double pwtBquark() const { return _pwtBquark; } /** * The diquark weight. */ double pwtDIquark() const { return _pwtDIquark; } //@} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * * The array _repwt is initialized using the interfaces to set different * weights for different meson multiplets and the constructHadronTable() * method called to complete the construction of the hadron tables. * * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} protected: /** * Construct the table of hadron data * This is the main method to initialize the hadron data (mainly the * weights associated to each hadron, taking into account its spin, * eventual isoscalar-octect mixing, singlet-decuplet factor). This is * the method that one should update when new or updated hadron data is * available. * * This class implements the construction of the basic table but can be * overridden if needed in inheriting classes. * * The rationale for factors used for diquarks involving different quarks can * be can be explained by taking a prototype example that in the exact SU(2) limit, * in which: * \f[m_u=m_d\f] * \f[M_p=M_n=M_\Delta\f] * and we will have equal numbers of u and d quarks produced. * Suppose that we weight 1 the diquarks made of the same * quark and 1/2 those made of different quarks, the fractions * of u and d baryons (p, n, Delta) we get are the following: * - \f$\Delta^{++}\f$: 1 possibility only u uu with weight 1 * - \f$\Delta^- \f$: 1 possibility only d dd with weight 1 * - \f$p,\Delta^+ \f$: 2 possibilities u ud with weight 1/2 * d uu with weight 1 * - \f$n,\Delta^0 \f$: 2 possibilities d ud with weight 1/2 * u dd with weight 1 * In the latter two cases, we have to take into account the * fact that p and n have spin 1/2 whereas Delta+ and Delta0 * have spin 3/2 therefore from phase space we get a double weight * for Delta+ and Delta0 relative to p and n respectively. * Therefore the relative amount of these baryons that is * produced is the following: * # p = # n = ( 1/2 + 1 ) * 1/3 = 1/2 * # Delta++ = # Delta- = 1 = ( 1/2 + 1) * 2/3 # Delta+ = # Delta0 * which is correct, and therefore the weight 1/2 for the * diquarks of different types of quarks is justified (at least * in this limit of exact SU(2) ). */ virtual void constructHadronTable(); /** * Access to the table of hadrons */ const HadronTable & table() const { return _table; } /** * Access to the list of partons */ const vector & partons() const { return _partons; } /** * Access the parton weights */ double pwt(long pid) const { map::const_iterator it = _pwt.find(abs(pid)); assert( it != _pwt.end() ); return it->second; } /** * Methods for the mixing of \f$I=0\f$ mesons */ //@{ /** * Return the probability of mixing for Octet-Singlet isoscalar mixing, * the probability of the * \f$\frac1{\sqrt{2}}(|u\bar{u}\rangle + |d\bar{d}\rangle)\f$ component * is returned. * @param angleMix The mixing angle in degrees (not radians) * @param order is 0 for no mixing, 1 for the first resonance of a pair, * 2 for the second one. * The mixing is defined so that for example with \f$\eta-\eta'\f$ mixing where * the mixing angle is \f$\theta=-23^0$ with $\eta\f$ as the first particle * and \f$\eta'\f$ the second one. * The convention used is * \f[\eta = \cos\theta|\eta_{\rm octet }\rangle * -\sin\theta|\eta_{\rm singlet}\rangle\f] * \f[\eta' = \sin\theta|\eta_{\rm octet }\rangle * -\cos\theta|\eta_{\rm singlet}\rangle\f] * with * \f[|\eta_{\rm singlet}\rangle = \frac1{\sqrt{3}} * \left[|u\bar{u}\rangle + |d\bar{d}\rangle + |s\bar{s}\rangle\right]\f] * \f[|\eta_{\rm octet }\rangle = \frac1{\sqrt{6}} * \left[|u\bar{u}\rangle + |d\bar{d}\rangle - 2|s\bar{s}\rangle\right]\f] */ double probabilityMixing(const double angleMix, const int order) const { static double convert=Constants::pi/180.0; if (order == 1) return sqr( cos( angleMix*convert + atan( sqrt(2.0) ) ) ); else if (order == 2) return sqr( sin( angleMix*convert + atan( sqrt(2.0) ) ) ); else return 1.; } /** * Returns the weight of given mixing state. * @param id The PDG code of the meson */ virtual double mixingStateWeight(long id) const; //@} /** * Calculates a special weight specific to a given hadron. * @param id The PDG code of the hadron */ double specialWeight(long id) const; /** * This method returns the proper sign ( > 0 hadron; < 0 anti-hadron ) * for the input PDG id idHad > 0, suppose to be made by the * two constituent particle pointers: par1 and par2 (both with proper sign). */ int signHadron(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr hadron) const; private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ HadronSelector & operator=(const HadronSelector &) = delete; private: /** * The PDG codes of the constituent particles allowed */ vector _partons; /** * The PDG codes of the hadrons which cannot be produced in the hadronization */ vector _forbidden; /** * The weights for the different quarks and diquarks */ //@{ /** * The probability of producting a down quark. */ double _pwtDquark; /** * The probability of producting an up quark. */ double _pwtUquark; /** * The probability of producting a strange quark. */ double _pwtSquark; /** * The probability of producting a charm quark. */ double _pwtCquark; /** * The probability of producting a bottom quark. */ double _pwtBquark; /** * The probability of producting a diquark. */ double _pwtDIquark; /** * Weights for quarks and diquarks. */ map _pwt; //@} /** * The mixing angles for the \f$I=0\f$ mesons containing light quarks */ //@{ /** * The \f$\eta-\eta'\f$ mixing angle */ double _etamix; /** * The \f$\phi-\omega\f$ mixing angle */ double _phimix; /** * The \f$h_1'-h_1\f$ mixing angle */ double _h1mix; /** * The \f$f_0(1710)-f_0(1370)\f$ mixing angle */ double _f0mix; /** * The \f$f_1(1420)-f_1(1285)\f$ mixing angle */ double _f1mix; /** * The \f$f'_2-f_2\f$ mixing angle */ double _f2mix; /** * The \f$\eta_2(1870)-\eta_2(1645)\f$ mixing angle */ double _eta2mix; /** * The \f$\phi(???)-\omega(1650)\f$ mixing angle */ double _omhmix; /** * The \f$\phi_3-\omega_3\f$ mixing angle */ double _ph3mix; /** * The \f$\eta(1475)-\eta(1295)\f$ mixing angle */ double _eta2Smix; /** * The \f$\phi(1680)-\omega(1420)\f$ mixing angle */ double _phi2Smix; //@} /** * The weights for the various meson multiplets to be used to supress the * production of particular states */ //@{ /** * The weights for the \f$\phantom{1}^1S_0\f$ multiplets */ vector _weight1S0; /** * The weights for the \f$\phantom{1}^3S_1\f$ multiplets */ vector _weight3S1; /** * The weights for the \f$\phantom{1}^1P_1\f$ multiplets */ vector _weight1P1; /** * The weights for the \f$\phantom{1}^3P_0\f$ multiplets */ vector _weight3P0; /** * The weights for the \f$\phantom{1}^3P_1\f$ multiplets */ vector _weight3P1; /** * The weights for the \f$\phantom{1}^3P_2\f$ multiplets */ vector _weight3P2; /** * The weights for the \f$\phantom{1}^1D_2\f$ multiplets */ vector _weight1D2; /** * The weights for the \f$\phantom{1}^3D_1\f$ multiplets */ vector _weight3D1; /** * The weights for the \f$\phantom{1}^3D_2\f$ multiplets */ vector _weight3D2; /** * The weights for the \f$\phantom{1}^3D_3\f$ multiplets */ vector _weight3D3; //@} /** * The weights for the excited meson multiplets */ vector > > _repwt; /** * Singlet and Decuplet weights */ //@{ /** * The singlet weight */ double _sngWt; /** * The decuplet weight */ double _decWt; //@} /** * The table of hadron data */ HadronTable _table; /** * Enums so arrays can be statically allocated */ //@{ /** * Defines values for array sizes. L,J,N max values for excited mesons. */ enum MesonMultiplets { Lmax = 3, Jmax = 4, Nmax = 4}; //@} /** * Option for the construction of the tables */ unsigned int _topt; /** * Which particles to produce for debugging purposes */ unsigned int _trial; /** * @name A parameter used for determining when clusters are too light. * * This parameter is used for setting the lower threshold, \f$ t \f$ as * \f[ t' = t(1 + r B^1_{\rm lim}) \f] * where \f$ r \f$ is a random number [0,1]. */ //@{ double _limBottom; double _limCharm; double _limExotic; //@} /** * Option for the selection of hadrons below the pair threshold */ unsigned int belowThreshold_; }; } #endif /* HERWIG_HadronSelector_H */