diff --git a/Decay/WeakCurrents/OmegaPiPiCurrent.cc b/Decay/WeakCurrents/OmegaPiPiCurrent.cc --- a/Decay/WeakCurrents/OmegaPiPiCurrent.cc +++ b/Decay/WeakCurrents/OmegaPiPiCurrent.cc @@ -1,328 +1,314 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the OmegaPiPiCurrent class. // #include "OmegaPiPiCurrent.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; -OmegaPiPiCurrent::OmegaPiPiCurrent() :mRes_(1.62*GeV) { +OmegaPiPiCurrent::OmegaPiPiCurrent() { + mRes_ = 1.67*GeV; wRes_ = 0.288*GeV; - gRes_ = 2.83; + gRes_ = 0.469*GeV; + mSigma_ = 0.6*GeV; wSigma_ = 1.0*GeV; - gSigma_ = 1.0; mf0_ = 0.98*GeV; - gPiPi_ = 0.331; - gKK_ = 0.144; - gf0_ = 0.85; + gSigma_ = 1.0 *GeV2; + gf0_ = 0.85*GeV2; + gPiPi_ = 0.165*GeV2; + gKK_ = 0.695*GeV2; addDecayMode(1,-1); addDecayMode(1,-1); setInitialModes(2); } IBPtr OmegaPiPiCurrent::clone() const { return new_ptr(*this); } IBPtr OmegaPiPiCurrent::fullclone() const { return new_ptr(*this); } void OmegaPiPiCurrent::persistentOutput(PersistentOStream & os) const { - os << ounit(mRes_,GeV) << ounit(wRes_,GeV) << gRes_ + os << ounit(mRes_,GeV) << ounit(wRes_,GeV) << ounit(gRes_,GeV) << ounit(mSigma_,GeV) << ounit(wSigma_,GeV) << ounit(mf0_,GeV) - << gPiPi_ << gKK_ << gSigma_ << gf0_; + << ounit(gPiPi_,GeV2) << ounit(gKK_,GeV2) << ounit(gSigma_,GeV2) << ounit(gf0_,GeV2); } void OmegaPiPiCurrent::persistentInput(PersistentIStream & is, int) { - is >> iunit(mRes_,GeV) >> iunit(wRes_,GeV) >> gRes_ + is >> iunit(mRes_,GeV) >> iunit(wRes_,GeV) >> iunit(gRes_,GeV) >> iunit(mSigma_,GeV) >> iunit(wSigma_,GeV) >> iunit(mf0_,GeV) - >> gPiPi_ >> gKK_ >> gSigma_ >> gf0_; + >> iunit(gPiPi_,GeV2) >> iunit(gKK_,GeV2) >> iunit(gSigma_,GeV2) >> iunit(gf0_,GeV2); } void OmegaPiPiCurrent::doinit() { WeakCurrent::doinit(); } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigOmegaPiPiCurrent("Herwig::OmegaPiPiCurrent", "HwWeakCurrents.so"); void OmegaPiPiCurrent::Init() { static ClassDocumentation documentation ("The OmegaPiPiCurrent class provides the current for I=0 omega pi pi"); static Parameter interfacemRes ("mRes", "The mass of the s-channel resonance", &OmegaPiPiCurrent::mRes_, GeV, 1.62*GeV, 0.*GeV, 10.*GeV, false, false, Interface::limited); static Parameter interfacewRes ("wRes", "The width of the s-channel resonance", &OmegaPiPiCurrent::wRes_, GeV, 0.288*GeV, 0.0*GeV, 10.0*GeV, false, false, Interface::limited); - static Parameter interfacegRes + static Parameter interfacegRes ("gRes", "The coupling of the s-channel resonance", - &OmegaPiPiCurrent::gRes_, 2.83, 0.0, 10.0, + &OmegaPiPiCurrent::gRes_, GeV, 2.83*GeV, 0.0*GeV, 10.0*GeV, false, false, Interface::limited); static Parameter interfacemSigma ("mSigma", "The mass of the Sigma", &OmegaPiPiCurrent::mSigma_, GeV, 0.6*GeV, 0.0*GeV, 10.0*GeV, false, false, Interface::limited); static Parameter interfacewSigma ("wSigma", "The width of the Sigma", &OmegaPiPiCurrent::wSigma_, GeV, 1.0*GeV, 0.0*GeV, 10.0*GeV, false, false, Interface::limited); - - static Parameter interfacegSigma + + static Parameter interfacegSigma ("gSigma", "The coupling of the Sigma resonance", - &OmegaPiPiCurrent::gSigma_, 1.0, 0.0, 10.0, + &OmegaPiPiCurrent::gSigma_, GeV2, 1.0*GeV2, 0.0*GeV2, 10.0*GeV2, false, false, Interface::limited); static Parameter interfacemf0 ("mf0", "The mass of the f_0(980)", &OmegaPiPiCurrent::mf0_, GeV, 0.98*GeV, 0.0*GeV, 10.0*GeV, false, false, Interface::limited); - - static Parameter interfacegf0 + + static Parameter interfacegf0 ("gf0", - "The coupling of the f_0(980) meson", - &OmegaPiPiCurrent::gf0_, 0.85, 0.0, 10.0, + "The coupling of the f0(980) resonance", + &OmegaPiPiCurrent::gf0_, GeV2, 1.0*GeV2, 0.0*GeV2, 10.0*GeV2, false, false, Interface::limited); - static Parameter interfacegPiPi + static Parameter interfacegPiPi ("gPiPi", "The coupling of the f_0(980) to pipi", - &OmegaPiPiCurrent::gPiPi_, .331, 0.0, 10.0, + &OmegaPiPiCurrent::gPiPi_, GeV2, 0.165*GeV2, 0.0*GeV2, 10.0*GeV2, false, false, Interface::limited); - static Parameter interfacegKK + static Parameter interfacegKK ("gKK", "The coupling of the f_0(980) to KK", - &OmegaPiPiCurrent::gKK_, .144, 0.0, 10.0, + &OmegaPiPiCurrent::gKK_, GeV2, 0.695*GeV2, 0.0*GeV2, 10.0*GeV2, false, false, Interface::limited); - } // complete the construction of the decay mode for integration bool OmegaPiPiCurrent::createMode(int icharge, tcPDPtr resonance, FlavourInfo flavour, unsigned int, PhaseSpaceModePtr mode, unsigned int iloc,int ires, PhaseSpaceChannel phase, Energy upp ) { // check the charge if(icharge!=0) return false; // check the total isospin if(flavour.I!=IsoSpin::IUnknown && flavour.I!=IsoSpin::IZero) return false; // check I_3 if(flavour.I3!=IsoSpin::I3Unknown && flavour.I3!=IsoSpin::I3Zero) return false; // and other flavour if(flavour.strange != Strangeness::Unknown and flavour.strange != Strangeness::Zero) return false; if(flavour.charm != Charm::Unknown and flavour.charm != Charm::Zero ) return false; if(flavour.bottom != Beauty::Unknown and flavour.bottom !=Beauty::Zero ) return false; // check that the mode is are kinematical allowed Energy min = getParticleData(ParticleID::omega)->massMin()+ 2.*getParticleData(ParticleID::pi0)->mass(); if(min>upp) return false; // resonances for the intermediate channels tPDVector res = {getParticleData(30223)}; tPDVector res2 = {getParticleData(9000221),getParticleData(9010221)}; // set up the integration channels; for(unsigned int ix=0;ixaddChannel((PhaseSpaceChannel(phase),ires,res[ix], ires+1,iloc+1,ires+1,res2[iy],ires+2,iloc+2,ires+2,iloc+3)); } } return true; } // the particles produced by the current tPDVector OmegaPiPiCurrent::particles(int icharge, unsigned int imode,int,int) { assert(icharge==0 && imode<=1); if(imode==0) return {getParticleData(ParticleID::omega), getParticleData(ParticleID::piplus), getParticleData(ParticleID::piminus)}; else if(imode==1) return {getParticleData(ParticleID::omega), getParticleData(ParticleID::pi0), getParticleData(ParticleID::pi0)}; else assert(false); } void OmegaPiPiCurrent::constructSpinInfo(ParticleVector decay) const { vector temp(3); for(unsigned int ix=0;ix<3;++ix) { temp[ix] = HelicityFunctions::polarizationVector(-decay[0]->momentum(), ix,Helicity::outgoing); } VectorWaveFunction::constructSpinInfo(temp,decay[0], outgoing,true,true); for(unsigned int ix=1;ix<3;++ix) ScalarWaveFunction::constructSpinInfo(decay[ix],outgoing,true); } // the hadronic currents vector OmegaPiPiCurrent::current(tcPDPtr resonance, FlavourInfo flavour, const int, const int ichan, Energy & scale, const tPDVector & , const vector & momenta, DecayIntegrator::MEOption) const { // no isospin/flavour here if(flavour.I!=IsoSpin::IUnknown && flavour.I!=IsoSpin::IZero) return vector(); if(flavour.I3!=IsoSpin::I3Unknown && flavour.I3!=IsoSpin::I3Zero) return vector(); if(flavour.strange != Strangeness::Unknown and flavour.strange != Strangeness::Zero) return vector(); if(flavour.charm != Charm::Unknown and flavour.charm != Charm::Zero ) return vector(); if(flavour.bottom != Beauty::Unknown and flavour.bottom !=Beauty::Zero ) return vector(); if(resonance and resonance->id()!=30223) return vector(); useMe(); // polarization vectors of the omega vector temp(3); for(unsigned int ix=0;ix<3;++ix) { temp[ix] = HelicityFunctions::polarizationVector(-momenta[0],ix,Helicity::outgoing); } // total momentum of the system Lorentz5Momentum q(momenta[0]+momenta[1]+momenta[2]); // overall hadronic mass q.rescaleMass(); scale=q.mass(); + Complex ii(0.,1.); Energy2 q2(q.m2()); // resonance factor for s channel resonance Energy2 mR2=sqr(mRes_); - Complex pre= mR2*gRes_/(q2-mR2 + Complex(0.,1.)*scale*wRes_); - // compute the form factor - complex formFactor(ZERO); + complex pre= mR2*gRes_/(q2-mR2 - ii*scale*wRes_); + + //cerr << "pre factor " << scale/GeV << " " << q2/GeV2 << " " << pre/GeV << "\n"; + // for(auto p : momenta) cerr << p/GeV << " " << p.m()/GeV << "\n"; - //formFactor = pre*scale*gSigma_; - - // needs to be multiplied by thing inside || in 1.9 - //cm energy for intermediate f0 channel + + // virtual mass energy for intermediate f0 channel Energy2 s1 = (momenta[1]+momenta[2]).m2(); Energy sqrs1 = sqrt(s1); //sigma meson Energy2 mSigma2 = sqr(mSigma_); - Complex Sigma_form = mSigma2/(s1-mSigma2 + Complex(0.,1.)*sqrs1*wSigma_); + Complex Sigma_form = gSigma_/(mSigma2 -s1 - ii*sqrs1*wSigma_); + + // compute the form factor + Energy2 mf02 = sqr(mf0_); + Energy mPi = getParticleData(211)->mass(); + Energy2 mPi2 = sqr(mPi); + + complex mGamma = gPiPi_*sqrt(max(0.,1.-4.*mPi2/s1)); + // cerr << "testing pi " << mGamma/GeV2 << " " << gPiPi_/GeV2 << " " << sqrt(max(0.,1.-4.*mPi2/s1))<< "\n"; + Energy2 mKp2 = sqr(getParticleData(321)->mass()); + double val = 1.-4.*mKp2/s1; + if(val>=0.) mGamma += 0.5*gKK_* sqrt( val); + else mGamma += 0.5*gKK_*ii*sqrt(-val); + Energy2 mK02 = sqr(getParticleData(311)->mass()); + val = 1.-4.*mK02/s1; + if(val>=0.) mGamma += 0.5*gKK_* sqrt( val); + else mGamma += 0.5*gKK_*ii*sqrt(-val); - //f0(980) following Phys. Lett. B 63, 224 (1976) - Energy2 mf02 = sqr(mf0_); - Energy mPi = getParticleData(211)->mass(); - Energy2 mPi2 = sqr(mPi); - Energy pcm_pipi = 0.5*sqrt(s1-4*mPi2); - Energy pcm_f0 = 0.5*sqrt(mf0_*mf0_-4*mPi2); - Energy GPiPi = gPiPi_*pcm_pipi; - Energy G0 = gPiPi_*pcm_f0; - Energy m_kaon0 = getParticleData(311)->mass(); - Energy2 m_kaon02 = sqr(m_kaon0); - Energy m_kaonP = getParticleData(321)->mass(); - Energy2 m_kaonP2 = sqr(m_kaonP); - Energy GKK; - Complex f0_form; - if(0.25*s1>m_kaon02){ - GKK = 0.5*gKK_*(sqrt(0.25*s1-m_kaon02)+sqrt(0.25*s1-m_kaonP2)); - f0_form = gf0_*mf0_*sqrt(G0*GPiPi)/(s1-mf02+Complex(0.,1.)*mf0_*(GPiPi+GKK)); - } - else{ - if(0.25*s1>m_kaonP2){ - f0_form = gf0_*mf0_*sqrt(G0*GPiPi)/(s1-mf02+Complex(0.,1.)*mf0_*(GPiPi+0.5*gKK_*(Complex(0.,1.)*sqrt(m_kaon02-0.25*s1)+sqrt(0.25*s1-m_kaonP2)))); - } - else{ - GKK = 0.5*gKK_*(sqrt(m_kaon02-0.25*s1)+sqrt(m_kaonP2-0.25*s1)); - f0_form = gf0_*mf0_*sqrt(G0*GPiPi)/(s1-mf02+Complex(0.,1.)*mf0_*(GPiPi+Complex(0.,1.)*GKK)); - } - } + Complex f0_form = gf0_/(mf02-s1-Complex(0.,1.)*mGamma); + // cerr << "f0 pieces " << mGamma/GeV2 << "\n"; + // cerr << "testing form factor " << s1/GeV2 << " " << f0_form << " " << Sigma_form << "\n"; + complex formFactor(ZERO); if(ichan<0) - formFactor = (Sigma_form+f0_form)*pre*scale*gSigma_; + formFactor = (Sigma_form+f0_form)*pre; else if(ichan==0) - formFactor = Sigma_form*pre*scale*gSigma_; + formFactor = Sigma_form*pre; else if(ichan==1) - formFactor = f0_form*pre*scale*gSigma_; -// // loop over the resonances -// for(unsigned int ix=imin;ix ret(3); for(unsigned int ix=0;ix<3;++ix) { ret[ix] = formFactor*temp[ix]; } return ret; } bool OmegaPiPiCurrent::accept(vector id) { if(id.size()!=3) return false; unsigned int nomega(0),npip(0),npim(0),npi0(0); for(unsigned int ix=0;ix id) { unsigned int npi0(0); for(unsigned int ix=0;ix current(tcPDPtr resonance, FlavourInfo flavour, const int imode, const int ichan,Energy & scale, const tPDVector & outgoing, const vector & momenta, DecayIntegrator::MEOption meopt) const; /** * Construct the SpinInfo for the decay products */ virtual void constructSpinInfo(ParticleVector decay) const; /** * Accept the decay. Checks the meson against the list * @param id The id's of the particles in the current. * @return Can this current have the external particles specified. */ virtual bool accept(vector id); /** * Return the decay mode number for a given set of particles in the current. * Checks the meson against the list * @param id The id's of the particles in the current. * @return The number of the mode */ virtual unsigned int decayMode(vector id); /** * Output the setup information for the particle database * @param os The stream to output the information to * @param header Whether or not to output the information for MySQL * @param create Whether or not to add a statement creating the object */ virtual void dataBaseOutput(ofstream & os,bool header,bool create) const; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ OmegaPiPiCurrent & operator=(const OmegaPiPiCurrent &) = delete; private: /** * Parameters for the \f$\omega(1650)\f$ */ //@{ /** * Mass of the resonance */ Energy mRes_; /** * Width of the resonance */ Energy wRes_; /** * Coupling of the resonance */ - double gRes_; + Energy gRes_; //@} /** * Parameters for the \f$f_0\f$ resonances */ //@{ /** * Mass of the \f$\sigma\f$ */ Energy mSigma_; /** * Width of the \f$\sigma\f$ */ Energy wSigma_; /** * Mass of the \f$f_0(980)\f$ */ Energy mf0_; /** * \f$f_0\f$ coupling to \f$\pi\pi\f$ */ - double gPiPi_; + Energy2 gPiPi_; /** * \f$f_0\f$ coupling to KK */ - double gKK_; + Energy2 gKK_; /** * Sigma coupling */ - double gSigma_; + Energy2 gSigma_; /** * f_0 coupling */ - double gf0_; + Energy2 gf0_; //@} }; } #endif /* Herwig_OmegaPiPiCurrent_H */ diff --git a/Hadronization/HwppSelector.cc b/Hadronization/HwppSelector.cc --- a/Hadronization/HwppSelector.cc +++ b/Hadronization/HwppSelector.cc @@ -1,251 +1,245 @@ // -*- C++ -*- // // HwppSelector.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the HwppSelector class. // #include "HwppSelector.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Utilities/Kinematics.h" #include "ThePEG/Utilities/Selector.h" #include "ThePEG/Repository/UseRandom.h" #include "CheckId.h" #include #include using namespace Herwig; DescribeClass describeHwppSelector("Herwig::HwppSelector","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; } void HwppSelector::persistentInput(PersistentIStream & is, int) { is >> _mode >> _enhanceSProb >> iunit(_m0Decay,GeV) >> _massMeasure; } void HwppSelector::Init() { static ClassDocumentation documentation ("The HwppSelector class implements the Herwig algorithm for selecting" " the hadrons", "The hadronization used the selection algorithm described in \\cite{Kupco:1998fx}.", "%\\cite{Kupco:1998fx}\n" "\\bibitem{Kupco:1998fx}\n" " A.~Kupco,\n" " ``Cluster hadronization in HERWIG 5.9,''\n" " arXiv:hep-ph/9906412.\n" " %%CITATION = HEP-PH/9906412;%%\n" ); // put useMe() only in correct place! static Switch interfaceMode ("Mode", "Which algorithm to use", &HwppSelector::_mode, 1, false, false); static SwitchOption interfaceModeKupco (interfaceMode, "Kupco", "Use the Kupco approach", 0); static SwitchOption interfaceModeHwpp (interfaceMode, "Hwpp", "Use the Herwig approach", 1); static Switch interfaceEnhanceSProb ("EnhanceSProb", "Option for enhancing strangeness", &HwppSelector::_enhanceSProb, 0, false, false); static SwitchOption interfaceEnhanceSProbNo (interfaceEnhanceSProb, "No", "No strangeness enhancement.", 0); static SwitchOption interfaceEnhanceSProbScaled (interfaceEnhanceSProb, "Scaled", "Scaled strangeness enhancement", 1); static SwitchOption interfaceEnhanceSProbExponential (interfaceEnhanceSProb, "Exponential", "Exponential strangeness enhancement", 2); static Switch interfaceMassMeasure ("MassMeasure", "Option to use different mass measures", &HwppSelector::_massMeasure,0,false,false); static SwitchOption interfaceMassMeasureMass (interfaceMassMeasure, "Mass", "Mass Measure", 0); static SwitchOption interfaceMassMeasureLambda (interfaceMassMeasure, "Lambda", "Lambda Measure", 1); static Parameter interfaceDecayMassScale ("DecayMassScale", "Cluster decay mass scale", &HwppSelector::_m0Decay, GeV, 1.0*GeV, 0.1*GeV, 50.*GeV, false, false, Interface::limited); } pair HwppSelector::chooseHadronPair(const Energy cluMass,tcPDPtr par1, tcPDPtr par2,tcPDPtr ) const { // if either of the input partons is a diquark don't allow diquarks to be // produced bool diquark = !(DiquarkMatcher::Check(par1->id()) || DiquarkMatcher::Check(par2->id())); bool quark = true; // if the Herwig algorithm if(_mode ==1) { if(UseRandom::rnd() > 1./(1.+pwtDIquark()) &&cluMass > massLightestBaryonPair(par1,par2)) { diquark = true; quark = false; } else { useMe(); diquark = false; quark = true; } } // weights for the different possibilities Energy weight, wgtsum(ZERO); // loop over all hadron pairs with the allowed flavours static vector hadrons; hadrons.clear(); for(unsigned int ix=0;ixiColour())) == 3 && !DiquarkMatcher::Check(quarktopick->id())) continue; if(!diquark && abs(int(quarktopick->iColour())) == 3 && DiquarkMatcher::Check(quarktopick->id())) continue; HadronTable::const_iterator tit1 = table().find(make_pair(abs(par1->id()),quarktopick->id())); HadronTable::const_iterator tit2 = table().find(make_pair(quarktopick->id(),abs(par2->id()))); // If not in table skip if(tit1 == table().end()||tit2==table().end()) continue; // tables empty skip const KupcoData & T1 = tit1->second; const KupcoData & T2 = tit2->second; if(T1.empty()||T2.empty()) continue; // if too massive skip if(cluMass <= T1.begin()->mass + T2.begin()->mass) continue; + // quark weight + double quarkWeight = pwt(quarktopick->id()); + if (quarktopick->id() == 3) { + // Scaling strangeness enhancement + 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; - // calculate the weight - double pwtstrange; - if (quarktopick->id() == 3) { - // Strangeness weight takes the automatic flat weight - pwtstrange = pwt(3); - // Scaling strangeness enhancement - if (_enhanceSProb == 1){ - double scale = double(sqr(_m0Decay/cluMass)); - pwtstrange = (_maxScale < scale) ? 0. : pow(pwtstrange,scale); - } - // Exponential strangeness enhancement - else if (_enhanceSProb == 2){ - Energy2 mass2; - Energy endpointmass = par1->mass() + par2->mass(); - // Choose to use either the cluster mass - // or to use the lambda measure - mass2 = (_massMeasure == 0) ? sqr(cluMass) : - sqr(cluMass) - sqr(endpointmass); - double scale = double(sqr(_m0Decay)/mass2); - pwtstrange = (_maxScale < scale) ? 0. : exp(-scale); - } - weight = pwtstrange * H1->overallWeight * H2->overallWeight * - Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass ); - } - else { - weight = pwt(quarktopick->id()) * H1->overallWeight * H2->overallWeight * - Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass ); - } + weight = quarkWeight * H1->overallWeight * H2->overallWeight * + Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass ); int signQ = 0; assert (par1 && quarktopick); assert (par2); assert(quarktopick->CC()); if(CheckId::canBeHadron(par1, quarktopick->CC()) && CheckId::canBeHadron(quarktopick, par2)) signQ = +1; else if(CheckId::canBeHadron(par1, quarktopick) && CheckId::canBeHadron(quarktopick->CC(), par2)) signQ = -1; else { cerr << "Could not make sign for" << par1->id()<< " " << quarktopick->id() << " " << par2->id() << "\n"; assert(false); } if (signQ == -1) quarktopick = quarktopick->CC(); // construct the object with the info Kupco a(quarktopick, H1->ptrData, H2->ptrData, weight); hadrons.push_back(a); wgtsum += weight; } } } if (hadrons.empty()) return make_pair(tcPDPtr(),tcPDPtr()); // select the hadron wgtsum *= UseRandom::rnd(); unsigned int ix=0; do { wgtsum-= hadrons[ix].weight; ++ix; } while(wgtsum > ZERO && ix < hadrons.size()); if(ix == hadrons.size() && wgtsum > ZERO) return make_pair(tcPDPtr(),tcPDPtr()); --ix; assert(hadrons[ix].idQ); int signHad1 = signHadron(par1, hadrons[ix].idQ->CC(), hadrons[ix].hadron1); int signHad2 = signHadron(par2, hadrons[ix].idQ, hadrons[ix].hadron2); assert( signHad1 != 0 && signHad2 != 0 ); return make_pair ( signHad1 > 0 ? hadrons[ix].hadron1 : tcPDPtr(hadrons[ix].hadron1->CC()), signHad2 > 0 ? hadrons[ix].hadron2 : tcPDPtr(hadrons[ix].hadron2->CC())); } diff --git a/MatrixElement/MEDM2Mesons.cc b/MatrixElement/MEDM2Mesons.cc --- a/MatrixElement/MEDM2Mesons.cc +++ b/MatrixElement/MEDM2Mesons.cc @@ -1,246 +1,231 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the MEDM2Mesons class. // #include "MEDM2Mesons.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h" using namespace Herwig; typedef LorentzVector > LorentzPolarizationVectorInvE; MEDM2Mesons::MEDM2Mesons() { cSMmed_ = {1.0,1.0,1.0}; } Energy2 MEDM2Mesons::scale() const { return sHat(); } unsigned int MEDM2Mesons::orderInAlphaS() const { return 0; } unsigned int MEDM2Mesons::orderInAlphaEW() const { return 0; } IBPtr MEDM2Mesons::clone() const { return new_ptr(*this); } IBPtr MEDM2Mesons::fullclone() const { return new_ptr(*this); } void MEDM2Mesons::doinit() { // make sure the current got initialised current_->init(); // max energy Energy Emax = generator()->maximumCMEnergy(); // loop over the modes int nmode=0; for(unsigned int imode=0;imodenumberOfModes();++imode) { // get the external particles for this mode int iq(0),ia(0); tPDVector out = current_->particles(0,imode,iq,ia); current_->decayModeInfo(imode,iq,ia); if(iq==2&&ia==-2) continue; PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(incomingA_,out,1., incomingB_,Emax)); PhaseSpaceChannel channel(mode); if(!current_->createMode(0,tcPDPtr(), FlavourInfo(), imode,mode,0,-1,channel,Emax)) continue; modeMap_[imode] = nmode; addMode(mode); ++nmode; } MEMultiChannel::doinit(); } void MEDM2Mesons::persistentOutput(PersistentOStream & os) const { os << current_ << modeMap_ << incomingA_ << incomingB_ << Mediator_ << cDMmed_ << cSMmed_ << ounit(MMed_,GeV); } void MEDM2Mesons::persistentInput(PersistentIStream & is, int) { is >> current_ >> modeMap_ >> incomingA_ >> incomingB_ >> Mediator_ >> cDMmed_ >> cSMmed_ >> iunit(MMed_,GeV); } //The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigMEDM2Mesons("Herwig::MEDM2Mesons", "Herwig.so"); void MEDM2Mesons::Init() { static ClassDocumentation documentation ("The MEDM2Mesons class simulates the annhilation of" " DM particles to mesons at low energy"); static Reference interfaceWeakCurrent ("WeakCurrent", "The reference for the decay current to be used.", &MEDM2Mesons::current_, false, false, true, false, false); static Reference interfaceIncomingA ("IncomingA", "First incoming particle", &MEDM2Mesons::incomingA_, false, false, true, false, false); static Reference interfaceIncomingB ("IncomingB", "Second incoming particle", &MEDM2Mesons::incomingB_, false, false, true, false, false); static Reference interfaceMediator_ ("Mediator", "DM mediator", &MEDM2Mesons::Mediator_, false, false, true, false, false); static Parameter interfacecDMmed ("cDMmed", "coupling of DM to dark mediator", &MEDM2Mesons::cDMmed_, 1.0, 0., 10., false, false, Interface::limited); static ParVector interfacecSMmed ("cSMmed", "coupling of SM to dark mediator", &MEDM2Mesons::cSMmed_, -1 , 1.0 , 0. , 10. , false, false, Interface::limited); } double MEDM2Mesons::me2(const int ichan) const { // compute the incoming current LorentzPolarizationVectorInvE lepton[2][2]; if(incomingA_->iSpin()==PDT::Spin1Half && incomingB_->iSpin()==PDT::Spin1Half) { SpinorWaveFunction em_in( meMomenta()[0],mePartonData()[0],incoming); SpinorBarWaveFunction ep_in( meMomenta()[1],mePartonData()[1],incoming); vector f1; vector a1; for(unsigned int ix=0;ix<2;++ix) { em_in.reset(ix); f1.push_back(em_in); ep_in.reset(ix); a1.push_back(ep_in); } // this should be coupling of DM to mediator/ mediator propagator complex mmed = Mediator_->mass(); complex mmed2 = sqr(mmed); complex mwid = Mediator_->width(); complex prop = sHat()-mmed2+Complex(0.,1.)*mmed*mwid; complex pre = cDMmed_/prop; for(unsigned ix=0;ix<2;++ix) { for(unsigned iy=0;iy<2;++iy) { lepton[ix][iy]= pre*f1[ix].dimensionedWave().vectorCurrent(a1[iy].dimensionedWave()); } } } // TODO think about other spins for the DM else assert(false); // work out the mapping for the hadron vector int nOut = int(meMomenta().size())-2; vector constants(nOut+1); vector iSpin(nOut); vector hadrons(nOut); int itemp(1); int ix(nOut); do { --ix; iSpin[ix] = mePartonData()[ix+2]->iSpin(); itemp *= iSpin[ix]; constants[ix] = itemp; hadrons[ix] = mePartonData()[ix+2]->id(); } while(ix>0); constants[nOut] = 1; // calculate the matrix element me_.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,iSpin)); // calculate the hadron current unsigned int imode = current_->decayMode(hadrons); Energy q = sqrt(sHat()); vector momenta(meMomenta() .begin()+2, meMomenta().end()); tPDVector out = mode(modeMap_.at(imode))->outgoing(); if(ichan<0) iMode(modeMap_.at(imode)); // get the hadronic currents for the I=1 and I=0 components vector hadronI0(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::Zero), imode,ichan,q,out,momenta,DecayIntegrator::Calculate)); vector hadronI1(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IOne, IsoSpin::I3Zero,Strangeness::Zero), imode,ichan,q,out,momenta,DecayIntegrator::Calculate)); vector - hadronssbar(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IOne, IsoSpin::I3Zero,Strangeness::ssbar), + hadronssbar(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::ssbar), imode,ichan,q,out,momenta,DecayIntegrator::Calculate)); + + + // compute the matrix element vector ihel(meMomenta().size()); double output(0.); int hI0_size = hadronI0.size(); int hI1_size = hadronI1.size(); int hss_size = hadronssbar.size(); int maxsize = max(max(hadronI0.size(),hadronI1.size()),hss_size); for(unsigned int hhel=0;hhel0;--ix) { ihel[ix+1]=(hhel%constants[ix-1])/constants[ix]; } // loop over the helicities of the incoming particles for(ihel[1]=0;ihel[1]<2;++ihel[1]){ for(ihel[0]=0;ihel[0]<2;++ihel[0]) { - Complex amp; + Complex amp = 0.; // work on coefficients for the I1 and I0 bits - if(hI0_size != 0 && hI1_size !=0){ - if(hss_size !=0){ - amp = lepton[ihel[0]][ihel[1]].dot((cSMmed_[0]-cSMmed_[1])*hadronI0[hhel]/sqrt(2)+(cSMmed_[0]+cSMmed_[1])*hadronI1[hhel]/sqrt(2)+cSMmed_[2]*hadronssbar[hhel]); - } - else { - amp = lepton[ihel[0]][ihel[1]].dot((cSMmed_[0]-cSMmed_[1])*hadronI0[hhel]/sqrt(2)+(cSMmed_[0]+cSMmed_[1])*hadronI1[hhel]/sqrt(2)); - } - } - else if(hI0_size != 0 && hI1_size == 0){ - if(hss_size !=0){ - amp = lepton[ihel[0]][ihel[1]].dot((cSMmed_[0]-cSMmed_[1])*hadronI0[hhel]/sqrt(2)+cSMmed_[2]*hadronssbar[hhel]); - } - else { - amp = lepton[ihel[0]][ihel[1]].dot((cSMmed_[0]-cSMmed_[1])*hadronI0[hhel]/sqrt(2)); - } - } - else { - if(hss_size !=0){ - amp = lepton[ihel[0]][ihel[1]].dot((cSMmed_[0]+cSMmed_[1])*hadronI1[hhel]/sqrt(2)+cSMmed_[2]*hadronssbar[hhel]); - } - else{ - amp = lepton[ihel[0]][ihel[1]].dot(cSMmed_[1]*hadronI1[hhel]); - } - } + if(hI0_size != 0 ) + amp += (cSMmed_[0]-cSMmed_[1])/sqrt(2.)*(lepton[ihel[0]][ihel[1]].dot(hadronI0[hhel])); + if(hI1_size !=0) + amp += (cSMmed_[0]+cSMmed_[1])/sqrt(2.)*(lepton[ihel[0]][ihel[1]].dot(hadronI1[hhel])); + if(hss_size !=0) + amp += cSMmed_[2]* (lepton[ihel[0]][ihel[1]].dot(hadronssbar[hhel])); me_(ihel)= amp; output += std::norm(amp); } } } // symmetry factors map ncount; double symmetry(1.); for(tPDPtr o : out) ncount[o->id()]+=1; for(map::const_iterator it=ncount.begin();it!=ncount.end();++it) { symmetry *= it->second; } // prefactors output *= 0.25*sqr(pow(sqrt(sHat())/q,int(momenta.size()-2))); return output/symmetry; } void MEDM2Mesons::constructVertex(tSubProPtr) { } diff --git a/Tests/Rivet/EE/EE-91.in b/Tests/Rivet/EE/EE-91.in --- a/Tests/Rivet/EE/EE-91.in +++ b/Tests/Rivet/EE/EE-91.in @@ -1,137 +1,138 @@ # -*- ThePEG-repository -*- ################################################## # LEP physics parameters (override defaults) ################################################## set /Herwig/Generators/EventGenerator:EventHandler:LuminosityFunction:Energy 91.2 ################################################## # select the analyses ################################################## # Validated ################################################## # ALEPH charged particle multiplicity insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_1991_S2435284 # OPAL charged particle multiplicity insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1992_I321190 # DELPHI charged particle multiplicity insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1991_I301657 # ALEPH main LEP I QCD summary paper insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_1996_S3486095 # ALEPH D* insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_1999_S4193598 # OPAL D* insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1995_I382219 # OPAL charged hadron analysis insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1994_S2927284 # OPAL Delta++ analysis insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1995_S3198391 # OPAL J/Psi analysis analysis insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1996_S3257789 # ALEPH eta/omega analysis insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_2002_S4823664 # OPAL K*0 analysis insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1997_S3608263 # OPAL flavour specific charged multiplicities etc insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1998_S3780481 # OPAL f_0,f_2 and phi production insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1998_S3702294 # OPAL gamma,pi0,eta,eta',rho+/-,a0+/- insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1998_S3749908 # OPAL K0 insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2000_S4418603 # OPAL K* +/- insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1993_I342766 # SLD flavour specific charged multiplicities etc insert /Herwig/Analysis/RivetAnalysis:Analyses 0 SLD_1996_S3398250 # SLD flavour specific charged multiplicities etc insert /Herwig/Analysis/RivetAnalysis:Analyses 0 SLD_1999_S3743934 # SLD flavour specific charged multiplicities etc insert /Herwig/Analysis/RivetAnalysis:Analyses 0 SLD_2004_S5693039 # OPAL event shapes and multiplicities at different energies insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2004_S6132243 # ALEPH jet and event shapes at many energies insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_2004_S5765862 # OPAL/JADE jet rates at many energies insert /Herwig/Analysis/RivetAnalysis:Analyses 0 JADE_OPAL_2000_S4300807 # DELPHI strange baryon production insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1995_S3137023 insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_2006_I719387 insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_2000_I524694 # DELPHI f_0, rho_0 and f_2 production insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1999_S3960137 # DELPHI K*0 and phi production insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1996_I420528 # DELPHI f'_2 production insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1996_I416741 # DELPHI pi0 insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1996_I401100 # DELPHI K- K*+/- insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1995_I377487 # DELPHI K+/-, p pbar insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1995_I394052 # DELPHI Delta++ insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1995_I399737 # DELPHI pi, K p, flavour seperated insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1998_I473409 # OPAL strange baryon production insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1997_S3396100 insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1997_I421977 # DELPHI tuning paper insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1996_S3430090 # DELPHI b quark insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_2011_I890503 # ALEPH b quark insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_2001_S4656318 # SLD b quark insert /Herwig/Analysis/RivetAnalysis:Analyses 0 SLD_2002_S4869273 # OPAL b quark insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2003_I599181 # PDG hadron multiplicities and ratios insert /Herwig/Analysis/RivetAnalysis:Analyses 0 PDG_HADRON_MULTIPLICITIES insert /Herwig/Analysis/RivetAnalysis:Analyses 0 PDG_HADRON_MULTIPLICITIES_RATIOS # OPAL from gluon paper insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2004_I648738 insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2004_I631361:PROCESS=QQ # L3 eta insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_1992_I336180 # L3 jet rates insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_2004_I652683 # ALEPH pi+-, K+- and (p, anti-p) insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_1995_I382179 # L3 eta' and omega specta insert /Herwig/Analysis/RivetAnalysis:Analyses 0 L3_1997_I427107 # charged mult, different # of jets insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1992_I334948 # charged mult, different rapidity regions insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1991_I324035 # lambda, lambdabar correlations insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1993_I360638 # delphi flavour sep charged # lambda, lambdabar correlations insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1997_I428178 # ALEPH pi0 insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_1997_I427131 # ALEPH Lambda polarization insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_1996_I415745 # DELPHI Lambda/K asymmetry insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1995_I382285 # OPAL b baryon polarization insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1998_I474012 # DELPHI b baryon polarization insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_2000_I513614 # ALEPH b baryon polarization insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ALEPH_1996_I402895 # OPAL rho and comega polarization insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2000_I502750 # DELPHI charged particle distributions insert /Herwig/Analysis/RivetAnalysis:Analyses 0 DELPHI_1999_I448370 # OPAL K0 spectrum insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_1995_I393503 ################################################## # unvalidated ################################################## # OPAL 4 jet angles insert /Herwig/Analysis/RivetAnalysis:Analyses 0 OPAL_2001_S4553896 ################################################## # MC ################################################## insert /Herwig/Analysis/RivetAnalysis:Analyses 0 MC_Bottom_Hadrons +insert /Herwig/Analysis/RivetAnalysis:Analyses 0 MC_Particle_Multiplicities diff --git a/Tests/Rivet/EE/EE-JPsi.in b/Tests/Rivet/EE/EE-JPsi.in --- a/Tests/Rivet/EE/EE-JPsi.in +++ b/Tests/Rivet/EE/EE-JPsi.in @@ -1,21 +1,23 @@ # -*- ThePEG-repository -*- # e+ e- -> J/Psi create Herwig::MEee2VectorMeson /Herwig/MatrixElements/MEJpsi HwMELepton.so set /Herwig/MatrixElements/MEJpsi:VectorMeson /Herwig/Particles/Jpsi set /Herwig/MatrixElements/MEJpsi:Coupling 11.43148 set /Herwig/MatrixElements/SubProcess:MatrixElements 0 /Herwig/MatrixElements/MEJpsi set EventGenerator:EventHandler:LuminosityFunction:Energy 3.096916 set /Herwig/Generators/EventGenerator:EventHandler:Cuts:MHatMin 0.2 set /Herwig/Particles/e-:PDF /Herwig/Partons/NoPDF set /Herwig/Particles/e+:PDF /Herwig/Partons/NoPDF do /Herwig/Particles/Jpsi:SelectDecayModes /Herwig/Particles/Jpsi/Jpsi->n0,nbar0; /Herwig/Particles/Jpsi/Jpsi->p+,pbar-; /Herwig/Particles/Jpsi/Jpsi->Sigma0,Sigmabar0; /Herwig/Particles/Jpsi/Jpsi->Lambda0,Lambdabar0; /Herwig/Particles/Jpsi/Jpsi->Sigma*-,Sigma*bar+; /Herwig/Particles/Jpsi/Jpsi->Sigma*0,Sigma*bar0; /Herwig/Particles/Jpsi/Jpsi->Sigma*+,Sigma*bar-; /Herwig/Particles/Jpsi/Jpsi->Xi-,Xibar+; /Herwig/Particles/Jpsi/Jpsi->Sigma*0,Sigma*bar0; /Herwig/Particles/Jpsi/Jpsi->Xi0,Xibar0; # J/psi-> lambda anti-lambda and sigma anti-sigma insert /Herwig/Analysis/RivetAnalysis:Analyses 0 BESIII_2017_I1510563 # J/psi -> p pbar and n nbar insert /Herwig/Analysis/RivetAnalysis:Analyses 0 BESIII_2012_I1113599 # J/Psi -> xi- and Sigma*+/- insert /Herwig/Analysis/RivetAnalysis:Analyses 0 BESIII_2016_I1422780 # J/Psi -> xi0 and Sigma*0 insert /Herwig/Analysis/RivetAnalysis:Analyses 0 BESIII_2017_I1506414 # J/psi-> lambda anti-lambda -insert /Herwig/Analysis/RivetAnalysis:Analyses 0 BESIII_2019_I1691850 \ No newline at end of file +insert /Herwig/Analysis/RivetAnalysis:Analyses 0 BESIII_2019_I1691850 +# J/psi-> lambda anti-Sigma0 +cc +insert /Herwig/Analysis/RivetAnalysis:Analyses 0 BESIII_2012_I1121378 diff --git a/Tests/python/LowEnergy-EE-NonPerturbative.in b/Tests/python/LowEnergy-EE-NonPerturbative.in --- a/Tests/python/LowEnergy-EE-NonPerturbative.in +++ b/Tests/python/LowEnergy-EE-NonPerturbative.in @@ -1,99 +1,103 @@ # -*- ThePEG-repository -*- read snippets/EECollider.in ################################################## # Selected the hard process ################################################## # leading-order processes ################################################## cd /Herwig/MatrixElements # e+ e- -> pi+pi- create Herwig::MEee2Mesons MEee2Pions HwMELeptonLowEnergy.so create Herwig::TwoPionCzyzCurrent /Herwig/Decays/TwoPionCzyzCurrent set MEee2Pions:WeakCurrent /Herwig/Decays/TwoPionCzyzCurrent # e+ e- -> K+K-/ K0K0 create Herwig::MEee2Mesons MEee2Kaons HwMELeptonLowEnergy.so create Herwig::TwoKaonCzyzCurrent /Herwig/Decays/TwoKaonCzyzCurrent set MEee2Kaons:WeakCurrent /Herwig/Decays/TwoKaonCzyzCurrent # e+ e- -> pi+ pi- pi0 create Herwig::MEee2Mesons MEee3Pions HwMELeptonLowEnergy.so create Herwig::ThreePionCzyzCurrent /Herwig/Decays/ThreePionCzyzCurrent set MEee3Pions:WeakCurrent /Herwig/Decays/ThreePionCzyzCurrent # e+ e- -> 2pi+ 2pi-, 2pi0, pi+ pi- create Herwig::MEee2Mesons MEee4Pions HwMELeptonLowEnergy.so create Herwig::FourPionCzyzCurrent /Herwig/Decays/FourPionCzyzCurrent set MEee4Pions:WeakCurrent /Herwig/Decays/FourPionCzyzCurrent # e+ e- -> eta pi+ pi- create Herwig::MEee2Mesons MEee2EtaPiPi HwMELeptonLowEnergy.so create Herwig::EtaPiPiCurrent /Herwig/Decays/EtaPiPiCurrent set MEee2EtaPiPi:WeakCurrent /Herwig/Decays/EtaPiPiCurrent # e+ e- -> eta' pi+ pi- create Herwig::MEee2Mesons MEee2EtaPrimePiPi HwMELeptonLowEnergy.so create Herwig::EtaPrimePiPiCurrent /Herwig/Decays/EtaPrimePiPiCurrent set MEee2EtaPrimePiPi:WeakCurrent /Herwig/Decays/EtaPrimePiPiCurrent # e+ e- -> omega pi (omega -> pi0 gamma) create Herwig::MEee2Mesons MEee2OmegaPi HwMELeptonLowEnergy.so create Herwig::TwoPionPhotonSNDCurrent /Herwig/Decays/OmegaPiCurrent #create Herwig::TwoPionPhotonCurrent /Herwig/Decays/OmegaPiCurrent set MEee2OmegaPi:WeakCurrent /Herwig/Decays/OmegaPiCurrent # e+ e- > pi0 gamma create Herwig::MEee2Mesons MEee2PiGamma HwMELeptonLowEnergy.so create Herwig::PionPhotonCurrent /Herwig/Decays/PiGammaCurrent set MEee2PiGamma:WeakCurrent /Herwig/Decays/PiGammaCurrent # e+e- -> eta gamma create Herwig::MEee2Mesons MEee2EtaGamma HwMELeptonLowEnergy.so create Herwig::EtaPhotonCurrent /Herwig/Decays/EtaGammaCurrent set MEee2EtaGamma:WeakCurrent /Herwig/Decays/EtaGammaCurrent # e+e- -> eta phi create Herwig::MEee2Mesons MEee2EtaPhi HwMELeptonLowEnergy.so create Herwig::EtaPhiCurrent /Herwig/Decays/EtaPhiCurrent set MEee2EtaPhi:WeakCurrent /Herwig/Decays/EtaPhiCurrent # e+e- -> eta omega create Herwig::MEee2Mesons MEee2EtaOmega HwMELeptonLowEnergy.so create Herwig::EtaOmegaCurrent /Herwig/Decays/EtaOmegaCurrent set MEee2EtaOmega:WeakCurrent /Herwig/Decays/EtaOmegaCurrent # e+e- > p pbar create Herwig::MEee2Mesons MEee2ppbar HwMELeptonLowEnergy.so create Herwig::WeakBaryonCurrent /Herwig/Decays/CzyzCurrent create Herwig::CzyzNucleonFormFactor /Herwig/Decays/CzyzFormFactor set /Herwig/Decays/CzyzCurrent:FormFactor /Herwig/Decays/CzyzFormFactor set MEee2ppbar:WeakCurrent /Herwig/Decays/CzyzCurrent # e+e- -> KKpi create Herwig::MEee2Mesons MEee2KKPi HwMELeptonLowEnergy.so create Herwig::KKPiCurrent /Herwig/Decays/KKPiCurrent set MEee2KKPi:WeakCurrent /Herwig/Decays/KKPiCurrent # e+e- -> phi pi create Herwig::MEee2Mesons MEee2PhiPi HwMELeptonLowEnergy.so create Herwig::PhiPiCurrent /Herwig/Decays/PhiPiCurrent set MEee2PhiPi:WeakCurrent /Herwig/Decays/PhiPiCurrent +# e+ e- -> omega pi pi +create Herwig::MEee2Mesons MEee2OmegaPiPi HwMELeptonLowEnergy.so +create Herwig::OmegaPiPiCurrent /Herwig/Decays/OmegaPiPiCurrent +set MEee2OmegaPiPi:WeakCurrent /Herwig/Decays/OmegaPiPiCurrent # default e+e- > q qbar (5 flavours d,u,s,c,b) ${processes} ################################################## # LEP physics parameters (override defaults) ################################################## cd /Herwig/Generators set EventGenerator:EventHandler:LuminosityFunction:Energy ${ECMS} set /Herwig/Generators/EventGenerator:EventHandler:Cuts:MHatMin 0.2 set /Herwig/Particles/e-:PDF /Herwig/Partons/NoPDF set /Herwig/Particles/e+:PDF /Herwig/Partons/NoPDF set /Herwig/Particles/pi0:Stable Stable set /Herwig/Particles/K_S0:Stable Stable cd /Herwig/Generators ################################################## ## prepare for Rivet analysis or HepMC output ## when running with parton shower ################################################## create ThePEG::RivetAnalysis /Herwig/Analysis/Rivet RivetAnalysis.so insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/Rivet ${ANALYSES} ################################################### # Save run for later usage with 'Herwig run' ################################################## set EventGenerator:MaxErrors 10000 set EventGenerator:EventHandler:StatLevel Full set EventGenerator:EventHandler:CascadeHandler NULL saverun Rivet-LowEnergy-EE-NonPerturbative-${ECMS} EventGenerator diff --git a/Tests/python/LowEnergy.py b/Tests/python/LowEnergy.py --- a/Tests/python/LowEnergy.py +++ b/Tests/python/LowEnergy.py @@ -1,388 +1,390 @@ #! /usr/bin/env python import yoda,os,math,subprocess,optparse from string import Template # get the path for the rivet data p = subprocess.Popen(["rivet-config", "--datadir"],stdout=subprocess.PIPE) path=p.communicate()[0].strip() #Define the arguments op = optparse.OptionParser(usage=__doc__) op.add_option("--process" , dest="processes" , default=[], action="append") op.add_option("--path" , dest="path" , default=path) op.add_option("--non-perturbative", dest="nonPerturbative" , default=False, action="store_true") op.add_option("--perturbative" , dest="perturbative" , default=False, action="store_true") op.add_option("--dest" , dest="dest" , default="Rivet") op.add_option("--list" , dest="list" , default=False, action="store_true") op.add_option("--flavour" , dest="flavour" , default="All" ) op.add_option("--plots" , dest="plot" , default=False, action="store_true") opts, args = op.parse_args() path=opts.path thresholds = [0.7,2.*.5,2.*1.87,2.*5.28] # the list of analyses and processes analyses={ 'KK' : {} , 'pipi' : {}, 'ppbar' : {}, "3pi" : {}, "etapipi" : {}, "etaprimepipi" : {} , "4pi" : {}, "etaPhi" : {}, "etaOmega" : {}, "2K1pi" : {} ,"2K2pi" : {} , "4K" : {}, "6m" : {}, "omegapi" : {} , "pigamma" : {}, "etagamma" : {}, "phipi" : {}, "omegapipi" : {} , "DD" : {} , "BB" : {}, "5pi" :{} , "LL" : {} } # pi+pi- analyses["pipi"]["KLOE_2009_I797438"] = ["d02-x01-y01"] analyses["pipi"]["KLOE_2005_I655225"] = ["d02-x01-y01"] analyses["pipi"]["KLOE2_2017_I1634981"] = ["d01-x01-y01"] analyses["pipi"]["BABAR_2009_I829441"] = ["d01-x01-y01"] analyses["pipi"]["DM1_1978_I134061"] = ["d01-x01-y01"] analyses["pipi"]["DM2_1989_I267118"] = ["d01-x01-y01"] analyses["pipi"]["CMD2_2007_I728302"] = ["d02-x01-y01"] analyses["pipi"]["CMD2_2006_I728191"] = ["d03-x01-y01"] analyses["pipi"]["BESIII_2016_I1385603"] = ["d01-x01-y01"] analyses["pipi"]["SND_2005_I686349"] = ["d01-x01-y01"] analyses["pipi"]["CMD_1985_I221309"] = ["d01-x01-y01","d02-x01-y01"] analyses["pipi"]["CMD2_2002_I568807"] = ["d01-x01-y02"] analyses["pipi"]["CMD2_1999_I498859"] = ["d01-x01-y01"] analyses['pipi']["CLEOC_2005_I693873"] = ["d01-x01-y01"] analyses['pipi']["ND_1991_I321108"] = ["d11-x01-y01"] analyses['pipi']["OLYA_1984_I208231"] = ["d01-x01-y01"] # K+K- and K_S^0 K_L^0 analyses['KK']["BESIII_2018_I1704558"] = ["d01-x01-y01"] analyses['KK']["BABAR_2013_I1238807"] = ["d01-x01-y01"] analyses['KK']["DM1_1981_I156053"] = ["d01-x01-y01"] analyses['KK']["DM1_1981_I156054"] = ["d01-x01-y01"] analyses['KK']["CLEOC_2005_I693873"] = ["d01-x01-y02"] analyses['KK']["BABAR_2015_I1383130"] = ["d01-x01-y04"] analyses['KK']["DM2_1988_I262690" ] = ["d01-x01-y01"] analyses['KK']["SND_2007_I755881"] = ["d01-x01-y01"] analyses['KK']["SND_2001_I533574"] = ["d01-x01-y01","d01-x01-y02","d01-x01-y03", "d02-x01-y01","d02-x01-y02","d02-x01-y03"] analyses['KK']["SND_2006_I720035"] = ["d01-x01-y01"] analyses['KK']["BABAR_2014_I1287920"] = ["d09-x01-y01"] analyses['KK']["CMD2_2003_I601222"] = ["d01-x01-y01"] analyses['KK']["CMD3_2016_I1444990"] = ["d01-x01-y06"] analyses['KK']["CMD2_1995_I406880"] = ["d01-x01-y01","d01-x01-y02"] analyses['KK']["CMD2_1999_I502164"] = ["d01-x01-y01","d02-x01-y01", "d03-x01-y01","d04-x01-y01"] analyses['KK']["CMD2_2008_I782516"] = ["d01-x01-y01","d02-x01-y01"] analyses['KK']["ND_1991_I321108"] = ["d12-x01-y01","d13-x01-y01"] analyses['KK']["OLYA_1981_I173076"] = ["d01-x01-y01"] analyses['KK']["SND_2016_I1484677"] = ["d01-x01-y01","d02-x01-y01"] # proton-antiproton analyses['ppbar']["BESIII_2019_I1736235"] = ["d01-x01-y01"] analyses['ppbar']["BESIII_2019_I1718337"] = ["d01-x01-y01"] analyses['ppbar']["BESIII_2015_I1358937"] = ["d01-x01-y05"] analyses['ppbar']["BABAR_2013_I1217421"] = ["d01-x01-y01"] analyses['ppbar']["SND_2014_I1321689"] = ["d01-x01-y01","d02-x01-y01"] analyses['ppbar']["CMD3_2016_I1385598"] = ["d01-x01-y06"] analyses['ppbar']["CLEOC_2005_I693873"] = ["d01-x01-y03"] analyses['ppbar']["BABAR_2006_I700020"] = ["d01-x01-y01","d02-x01-y01"] analyses['ppbar']["DM2_1983_I190558"] = ["d01-x01-y01"] analyses["ppbar"]["DM2_1990_I297706"] = ["d01-x01-y01"] analyses["ppbar"]["DM1_1979_I141565"] = ["d01-x01-y01"] analyses["ppbar"]["FENICE_1998_I471263"] = ["d01-x01-y01"] analyses["ppbar"]["FENICE_1994_I377833"] = ["d01-x01-y01"] analyses['ppbar']["BESII_2005_I685906"] = ["d01-x01-y01"] analyses['ppbar']["BESIII_2014_I1286898"] = ["d01-x01-y06"] # pi0 gamma analyses["pigamma"]["SND_2018_I1694988"] = ["d01-x01-y01"] analyses["pigamma"]["SND_2016_I1418483"] = ["d01-x01-y05"] analyses["pigamma"]["SND_2003_I612867"] = ["d01-x01-y01"] analyses["pigamma"]["CMD2_2005_I658856"] = ["d02-x01-y01"] analyses["pigamma"]["SND_2000_I524221"] = ["d01-x01-y02"] # eta gamma analyses["etagamma"]["CMD2_2005_I658856"] = ["d01-x01-y01"] analyses["etagamma"]["SND_2006_I717778" ] = ["d01-x01-y01","d02-x01-y01"] analyses["etagamma"]["SND_2014_I1275333"] = ["d01-x01-y01"] analyses["etagamma"]["SND_2000_I524221"] = ["d01-x01-y01"] analyses["etagamma"]["CMD2_1999_I503154"] = ["d01-x01-y01"] analyses["etagamma"]["CMD2_2001_I554522"] = ["d01-x01-y01"] analyses['etagamma']["CMD2_1995_I406880"] = ["d01-x01-y04"] analyses['etagamma']["BABAR_2006_I716277"] = ["d01-x01-y01"] # 3 meson analyses["3pi"]["BABAR_2004_I656680"] = ["d01-x01-y01"] analyses["3pi"]["SND_2002_I582183"] = ["d01-x01-y01"] analyses["3pi"]["SND_2003_I619011"] = ["d01-x01-y01"] analyses["3pi"]["SND_1999_I508003"] = ["d01-x01-y01"] analyses["3pi"]["SND_2001_I533574"] = ["d01-x01-y04","d02-x01-y04"] analyses["3pi"]["CMD2_2000_I523691"] = ["d01-x01-y01"] analyses["3pi"]["CMD2_1998_I480170"] = ["d01-x01-y01"] analyses['3pi']["CMD2_1995_I406880"] = ["d01-x01-y03"] analyses['3pi']["DM2_1992_I339265" ] = ["d01-x01-y01"] analyses['3pi']["DM1_1980_I140174" ] = ["d01-x01-y01"] analyses['3pi']["ND_1991_I321108"] = ["d05-x01-y01","d10-x01-y04"] analyses['3pi']["GAMMAGAMMA_1981_I158474"] = ["d01-x01-y01"] analyses["3pi"]["CLEO_2006_I691720"] = ["d01-x01-y01"] analyses["3pi"]["SND_2015_I1389908"] = ["d01-x01-y01"] # etapipi analyses["etapipi"]["BABAR_2018_I1700745"] = ["d02-x01-y01"] analyses["etapipi"]["BABAR_2018_I1647139"] = ["d01-x01-y01"] analyses["etapipi"]["SND_2015_I1332929"] = ["d01-x01-y01"] analyses["etapipi"]["SND_2018_I1638368"] = ["d01-x01-y01"] analyses["etapipi"]["BABAR_2007_I758568"] = ["d01-x01-y01","d02-x01-y01"] analyses["etapipi"]["CMD2_2000_I532970"] = ["d02-x01-y01"] analyses["etapipi"]["DM2_1988_I264144"] = ["d01-x01-y01"] analyses['etapipi']["ND_1991_I321108"] = ["d06-x01-y01","d14-x01-y01"] # eta Phi analyses["etaPhi"]["BABAR_2008_I765258"] = ["d04-x01-y01","d05-x01-y01"] analyses["etaPhi"]["SND_2018_I1693737"] = ["d01-x01-y01"] analyses["etaPhi"]["BABAR_2017_I1511276"] = ["d03-x01-y01"] analyses["etaPhi"]["SND_2019_I1726419"] = ["d01-x01-y03"] # eta Omega analyses["etaOmega"]["SND_2016_I1473343"] = ["d01-x01-y01"] analyses["etaOmega"]["BABAR_2006_I709730"] = ["d02-x01-y01"] analyses["etaOmega"]["SND_2019_I1726419"] = ["d01-x01-y01","d01-x01-y02"] analyses["etaOmega"]["CMD3_2017_I1606078"] = ["d01-x01-y01","d01-x01-y02"] # 4 pions analyses["4pi"]["BABAR_2017_I1621593"] = ["d01-x01-y01","d02-x01-y01"] analyses["4pi"]["BABAR_2012_I1086164"] = ["d01-x01-y01"] analyses["4pi"]["CMD2_2000_I531667"] = ["d01-x01-y01"] analyses["4pi"]["CMD2_2004_I648023"] = ["d01-x01-y01"] analyses["4pi"]["BABAR_2005_I676691"] = ["d01-x01-y01"] analyses["4pi"]["CMD2_2000_I511375"] = ["d01-x01-y01"] analyses["4pi"]["CMD2_1999_I483994"] = ["d01-x01-y01","d02-x01-y01"] analyses["4pi"]["BESII_2008_I801210"] = ["d01-x01-y01"] analyses["4pi"]["KLOE_2008_I791841"] = ["d01-x01-y01"] analyses['4pi']["ND_1991_I321108"] = ["d07-x01-y01","d08-x01-y01","d10-x01-y01","d10-x01-y02"] analyses['4pi']["BESII_2007_I750713"] = ["d01-x01-y03"] analyses['4pi']["SND_2001_I579319"] = ["d01-x01-y01","d02-x01-y01"] analyses['4pi']["DM1_1982_I168552"] = ["d01-x01-y01"] analyses['4pi']["DM1_1979_I132828"] = ["d01-x01-y01"] analyses['4pi']["GAMMAGAMMA_1980_I153382"] = ["d01-x01-y01"] analyses['4pi']["GAMMAGAMMA_1981_I158474"] = ["d01-x01-y02"] # (these are omega(-> pi0 gamma) pi0) analyses["omegapi"]["SND_2016_I1489182" ] = ["d01-x01-y01"] analyses["omegapi"]["SND_2000_I527752" ] = ["d01-x01-y01"] analyses["omegapi"]["SND_2000_I503946" ] = ["d01-x01-y01"] analyses["omegapi"]["CMD2_2003_I616446" ] = ["d01-x01-y01"] # non omega analyses["omegapi"]["SND_2002_I587084" ] = ["d01-x01-y01"] analyses["omegapi"]["CMD2_2004_I630009" ] = ["d01-x01-y01"] analyses["omegapi"]["KLOE_2008_I791841"] = ["d02-x01-y01"] # from 4pion analyses["omegapi"]["CMD2_1999_I483994" ] = ["d03-x01-y01"] analyses['omegapi']["ND_1991_I321108"] = ["d01-x01-y01","d02-x01-y01","d03-x01-y01", "d04-x01-y01","d10-x01-y03"] # 5 pion and related -analyses["omegapipi"]["DM1_1981_I166964"] = ["d01-x01-y01"] -analyses["omegapipi"]["DM2_1992_I339265"]= ["d02-x01-y01"] -analyses["omegapipi"]["CMD2_2000_I532970"] = ["d01-x01-y01"] +analyses["omegapipi"]["DM1_1981_I166964" ] = ["d01-x01-y01"] +analyses["omegapipi"]["DM2_1992_I339265" ] = ["d02-x01-y01"] +analyses["omegapipi"]["CMD2_2000_I532970" ] = ["d01-x01-y01"] analyses["omegapipi"]["BABAR_2018_I1700745"] = ["d03-x01-y01"] +analyses["omegapipi"]["BABAR_2007_I758568" ] = ["d03-x01-y01","d04-x01-y01"] analyses["5pi"]["CMD2_2000_I532970"] = ["d03-x01-y01"] analyses["5pi"]["BABAR_2007_I758568"] = ["d01-x01-y01"] analyses['5pi']["ND_1991_I321108"] = ["d14-x01-y01"] analyses['5pi']["GAMMAGAMMA_1981_I158474"] = ["d01-x01-y03"] analyses["5pi"]["BABAR_2018_I1700745"] = ["d01-x01-y01"] # 2K 1 pi analyses["2K1pi"]["BABAR_2008_I765258"] = ["d01-x01-y01","d02-x01-y01","d03-x01-y01"] analyses["2K1pi"]["DM1_1982_I176801"] = ["d01-x01-y01"] analyses["2K1pi"]["DM2_1991_I318558"] = ["d01-x01-y01","d02-x01-y01"] analyses["2K1pi"]["BESII_2008_I801208"] = ["d01-x01-y01"] analyses["2K1pi"]["SND_2018_I1637194"] = ["d01-x01-y01"] analyses["2K1pi"]["BESIII_2018_I1691798"] = ["d01-x01-y01"] analyses["2K1pi"]["BABAR_2017_I1511276"] = ["d01-x01-y01"] analyses["phipi"]["BABAR_2017_I1511276"] = ["d01-x01-y01","d02-x01-y01"] analyses["phipi"]["BABAR_2008_I765258"] = ["d02-x01-y01","d03-x01-y01"] # 2K 2 pi analyses["2K2pi"]["DM1_1982_I169382"] = ["d01-x01-y01"] analyses["2K2pi"]["BABAR_2005_I676691"] = ["d02-x01-y01"] analyses["2K2pi"]["BABAR_2014_I1287920"] = ["d09-x01-y01","d10-x01-y01","d11-x01-y01"] analyses["2K2pi"]["BABAR_2012_I892684"] = ["d01-x01-y01","d02-x01-y01","d03-x01-y01", "d04-x01-y01","d05-x01-y01", "d06-x01-y01","d07-x01-y01"] analyses["2K2pi"]["BABAR_2007_I747875"] = ["d01-x01-y01","d02-x01-y01","d03-x01-y01", "d04-x01-y01","d05-x01-y01","d07-x01-y01"] analyses["2K2pi"]["BESII_2008_I801210"] = ["d01-x01-y02"] analyses["2K2pi"]["BESII_2008_I801208"] = ["d01-x01-y02"] analyses["2K2pi"]["BELLE_2009_I809630"] = ["d01-x01-y01"] analyses["2K2pi"]["CMD3_2016_I1395968"] = ["d01-x01-y01"] analyses['2K2pi']["BESII_2007_I750713"] = ["d01-x01-y04"] analyses["2K2pi"]["BABAR_2017_I1511276"] = ["d03-x01-y01","d04-x01-y01"] analyses["2K2pi"]["BABAR_2017_I1591716"] = ["d01-x01-y01","d02-x01-y01"] analyses['2K2pi']["BESIII_2018_I1699641"] = ["d01-x01-y01","d02-x01-y01"] # 4K analyses["4K"]["BESIII_2019_I1743841"] = ["d01-x01-y01","d02-x01-y01"] analyses["4K"]["BABAR_2005_I676691"] = ["d03-x01-y01"] analyses["4K"]["BABAR_2014_I1287920"] = ["d12-x01-y01"] analyses["4K"]["BABAR_2012_I892684"] = ["d08-x01-y01"] analyses["4K"]["BABAR_2007_I747875"] = ["d07-x01-y01"] analyses['4K']["BESII_2007_I750713"] = ["d01-x01-y06","d01-x01-y07"] # 6 mesons analyses["6m"]["CMD3_2013_I1217420"] = ["d01-x01-y01"] analyses["6m"]["SND_2019_I1726419"] = ["d01-x01-y01","d01-x01-y04"] analyses["6m"]["CMD3_2017_I1606078"] = ["d01-x01-y03","d01-x01-y04"] analyses["6m"]["CMD3_2019_I1720610"] = ["d01-x01-y01","d01-x01-y02","d01-x01-y03"] analyses["6m"]["BABAR_2018_I1700745"] = ["d04-x01-y01","d05-x01-y01"] analyses["6m"]["SND_2016_I1471515"] = ["d01-x01-y06"] analyses["6m"]["DM1_1981_I166353"] = ["d01-x01-y01"] analyses["6m"]["BABAR_2006_I709730"] = ["d01-x01-y01","d02-x01-y01","d03-x01-y01"] -analyses["6m"]["BABAR_2007_I758568"] = ["d03-x01-y01","d04-x01-y01","d05-x01-y01","d07-x01-y01", +analyses["6m"]["BABAR_2007_I758568"] = ["d05-x01-y01","d07-x01-y01", "d08-x01-y01","d09-x01-y01","d10-x01-y01","d11-x01-y01"] analyses["etaprimepipi"]["BABAR_2007_I758568"] = ["d05-x01-y01","d06-x01-y01"] analyses["6m"]["BESII_2007_I763880"] = ["d01-x01-y01","d01-x01-y02","d01-x01-y03","d01-x01-y04", "d01-x01-y05","d01-x01-y06","d01-x01-y07"] analyses["6m"]["BESII_2007_I762901"] = ["d01-x01-y01","d01-x01-y02","d01-x01-y03","d01-x01-y04", "d01-x01-y06","d01-x01-y07","d01-x01-y08","d01-x01-y09","d01-x01-y10"] analyses["6m"]["CLEO_2006_I691720"] = ["d01-x01-y02","d01-x01-y03","d01-x01-y04","d01-x01-y05", "d01-x01-y07","d01-x01-y08","d01-x01-y09","d01-x01-y10","d01-x01-y11", "d01-x01-y12","d01-x01-y13","d01-x01-y14","d01-x01-y15","d01-x01-y17"] analyses["6m"]["BESII_2008_I801210"] = ["d01-x01-y03","d01-x01-y04","d01-x01-y05"] analyses["6m"]["BESII_2008_I801208"] = ["d01-x01-y03","d01-x01-y04","d01-x01-y05","d01-x01-y06"] analyses["6m"]["JADE_1979_I142874"] = ["d02-x01-y01"] analyses["6m"]["MARKI_1982_I169326"] = ["d06-x01-y01"] analyses["6m"]["MARKI_1975_I100592"] = ["d01-x01-y01","d02-x01-y01"] analyses['6m']["BESII_2007_I750713"] = ["d01-x01-y08","d01-x01-y09","d01-x01-y11", "d01-x01-y12","d01-x01-y13","d01-x01-y14", "d01-x01-y15","d01-x01-y16","d01-x01-y17","d01-x01-y18"] analyses['6m']["SND_2016_I1473343"] = ["d01-x01-y01"] # DD analyses["DD"]["BELLE_2007_I723333" ] = ["d01-x01-y01","d02-x01-y01"] analyses["DD"]["BELLE_2007_I756012" ] = ["d01-x01-y01"] analyses["DD"]["BELLE_2007_I756643" ] = ["d01-x01-y01"] analyses["DD"]["BELLE_2008_I757220" ] = ["d01-x01-y01","d02-x01-y01"] analyses["DD"]["BELLE_2008_I759073" ] = ["d01-x01-y01"] analyses["DD"]["BABAR_2008_I776519" ] = ["d01-x01-y01","d01-x01-y02"] analyses["DD"]["BELLE_2008_I791660" ] = ["d01-x01-y01"] analyses["DD"]["BELLE_2013_I1225975"] = ["d01-x01-y01"] analyses["DD"]["BELLE_2014_I1282602"] = ["d01-x01-y01"] analyses["DD"]["BELLE_2015_I1324785"] = ["d01-x01-y01"] analyses["DD"]["BESIII_2016_I1457597"] = ["d01-x01-y07"] analyses["DD"]["BESIII_2015_I1355215"] = ["d01-x01-y10"] analyses["DD"]["BESIII_2015_I1377204"] = ["d01-x01-y10"] analyses["DD"]["BESIII_2016_I1495838"] = ["d01-x01-y01","d02-x01-y01"] analyses["DD"]["CRYSTAL_BALL_1986_I238081"] = ["d02-x01-y01"] analyses["DD"]["CLEOC_2008_I777917"] = ["d01-x01-y01","d01-x01-y02","d01-x01-y03", "d02-x01-y01","d02-x01-y02","d02-x01-y03", "d03-x01-y01","d03-x01-y02","d03-x01-y03", "d04-x01-y01","d04-x01-y02", "d05-x01-y01","d05-x01-y02"] analyses["DD"]["BELLE_2017_I1613517"] = ["d01-x01-y01","d01-x01-y02"] analyses["DD"]["BESIII_2014_I1323621"] = ["d01-x01-y01"] analyses["DD"]["BESIII_2015_I1406939"] = ["d02-x01-y06","d03-x01-y06"] analyses["DD"]["BESIII_2017_I1628093"] = ["d01-x01-y01"] analyses["DD"]["BESIII_2019_I1723934"] = ["d01-x01-y01"] - +analyses["DD"]["BESIII_2019_I1756876"] = ["d01-x01-y09","d01-x01-y10"] # BB analyses["BB"]["BELLE_2016_I1389855"] = ["d01-x01-y02","d01-x01-y03"] analyses["BB"]["BELLE_2008_I764099"] = ["d01-x01-y01","d02-x01-y01", "d03-x01-y01","d04-x01-y01"] analyses["BB"]["CLEO_1999_I474676"] = ["d01-x01-y01","d01-x01-y02"] analyses["BB"]["BABAR_2001_I558091"] = ["d01-x01-y01"] analyses["BB"]["CLEO_1991_I29927"] = ["d01-x01-y01"] analyses["LL"]["BESIII_2018_I1627871"] = ["d01-x01-y01"] analyses["LL"]["DM2_1990_I297706"] = ["d02-x01-y01"] analyses["LL"]["BABAR_2007_I760730"] = ["d01-x01-y01","d02-x01-y01","d03-x01-y01"] # list the analysis if required and quit() if "all" in opts.processes : processes = sorted(list(analyses.keys())) else : processes = sorted(list(set(opts.processes))) if(opts.list) : for process in processes : print " ".join(analyses[process]) quit() if(opts.plot) : output="" for process in processes: for analysis in analyses[process] : for plot in analyses[process][analysis]: output+= " -m/%s/%s" % (analysis,plot) print output quit() # mapping of process to me to use -me = { "pipi" : "MEee2Pions", - "KK" : "MEee2Kaons", - "3pi" : "MEee3Pions", - "4pi" : "MEee4Pions", - "etapipi" : "MEee2EtaPiPi", - "etaprimepipi" : "MEee2EtaPrimePiPi", - "etaPhi" : "MEee2EtaPhi", - "etaOmega" : "MEee2EtaOmega", - "omegapi" : "MEee2OmegaPi", - "phipi" : "MEee2PhiPi", - "pigamma" : "MEee2PiGamma", - "etagamma" : "MEee2EtaGamma", - "ppbar" : "MEee2ppbar", - "2K1pi" : "MEee2KKPi"} +me = { "pipi" : "MEee2Pions", + "KK" : "MEee2Kaons", + "3pi" : "MEee3Pions", + "4pi" : "MEee4Pions", + "etapipi" : "MEee2EtaPiPi", + "etaprimepipi" : "MEee2EtaPrimePiPi", + "etaPhi" : "MEee2EtaPhi", + "etaOmega" : "MEee2EtaOmega", + "omegapi" : "MEee2OmegaPi", + "omegapipi" : "MEee2OmegaPiPi", + "phipi" : "MEee2PhiPi", + "pigamma" : "MEee2PiGamma", + "etagamma" : "MEee2EtaGamma", + "ppbar" : "MEee2ppbar", + "2K1pi" : "MEee2KKPi"} # energies we need energies={} def nearestEnergy(en) : Emin=0 delta=1e30 anals=[] for val in energies : if(abs(val-en)200) : energy *= 0.001 emin,delta,anals = nearestEnergy(energy) if(energy in energies) : if(analysis not in energies[energy][1]) : energies[energy][1].append(analysis) if(matrix!="" and matrix not in energies[energy][0]) : energies[energy][0].append(matrix) elif(delta<1e-7) : if(analysis not in anals[1]) : anals[1].append(analysis) if(matrix!="" and matrix not in anals[0]) : anals[0].append(matrix) else : if(matrix=="") : energies[energy]=[[],[analysis]] else : energies[energy]=[[matrix],[analysis]] with open("python/LowEnergy-EE-Perturbative.in", 'r') as f: templateText = f.read() perturbative=Template( templateText ) with open("python/LowEnergy-EE-NonPerturbative.in", 'r') as f: templateText = f.read() nonPerturbative=Template( templateText ) targets="" for energy in sorted(energies) : anal="" for analysis in energies[energy][1]: anal+= "insert /Herwig/Analysis/Rivet:Analyses 0 %s\n" %analysis proc="" for matrix in energies[energy][0] : proc+="insert SubProcess:MatrixElements 0 %s\n" % matrix proc+="set %s:Flavour %s\n" % (matrix,opts.flavour) maxflavour =5 if(energy thresholds[0]) : inputPerturbative = perturbative.substitute({"ECMS" : "%8.6f" % energy, "ANALYSES" : anal, "lepton" : "", "maxflavour" : maxflavour}) with open(opts.dest+"/Rivet-LowEnergy-EE-Perturbative-%8.6f.in" % energy ,'w') as f: f.write(inputPerturbative) targets += "Rivet-LowEnergy-EE-Perturbative-%8.6f.yoda " % energy # input file for currents if(opts.nonPerturbative and proc!="") : inputNonPerturbative = nonPerturbative.substitute({"ECMS" : "%8.6f" % energy, "ANALYSES" : anal, "processes" : proc}) with open(opts.dest+"/Rivet-LowEnergy-EE-NonPerturbative-%8.6f.in" % energy ,'w') as f: f.write(inputNonPerturbative) targets += "Rivet-LowEnergy-EE-NonPerturbative-%8.6f.yoda " % energy print targets diff --git a/Tests/python/mergeLowEnergy.py b/Tests/python/mergeLowEnergy.py --- a/Tests/python/mergeLowEnergy.py +++ b/Tests/python/mergeLowEnergy.py @@ -1,67 +1,72 @@ #! /usr/bin/env python import yoda,glob,math,optparse op = optparse.OptionParser(usage=__doc__) opts, args = op.parse_args() if(len(args)!=1) : print 'Must be one and only 1 name' quit() for runType in ["NonPerturbative","Perturbative"]: outhistos={} for fileName in glob.glob("Rivet-LowEnergy-EE-%s-*.yoda" % runType): energy = float(fileName.split("-")[-1].strip(".yoda")) energyMeV = energy*1000. aos = yoda.read(fileName) for hpath,histo in aos.iteritems(): if("/_" in hpath or "TMP" in hpath) : continue if(type(histo)==yoda.core.Histo1D) : outhistos[hpath] = histo continue # create histo if it doesn't exist elif(hpath not in outhistos) : - outhistos[hpath] = yoda.core.Scatter2D(histo.path, - histo.title) + title="" + path="" + if hasattr(histo, 'title'): + title=histo.title() + if hasattr(histo, 'path'): + path=histo.path() + outhistos[hpath] = yoda.core.Scatter2D(path,title) matched = False - for i in range(0,aos[hpath].numPoints) : - x = aos[hpath].points[i].x + for i in range(0,aos[hpath].numPoints()) : + x = aos[hpath].points()[i].x() delta=1e-5 if("KLOE_2009_I797438" in hpath or "KLOE_2005_I655225" in hpath or "KLOE2_2017_I1634981" in hpath or "FENICE_1994_I377833" in hpath or "FENICE_1996_I426675" in hpath): x=math.sqrt(x) delta=1e-3 if(abs(x-energy)<1e-3*delta or abs(x-energyMeV) xmin and energy < xmax) or (energyMeV > xmin and energyMeV < xmax) ) : duplicate = False - for j in range(0,outhistos[hpath].numPoints) : - if(outhistos[hpath].points[j].x==aos[hpath].points[i].x) : + for j in range(0,outhistos[hpath].numPoints()) : + if(outhistos[hpath].points()[j].x()==aos[hpath].points()[i].x()) : duplicate = True break if(not duplicate) : - outhistos[hpath].addPoint(aos[hpath].points[i]) + outhistos[hpath].addPoint(aos[hpath].points()[i]) break if len(outhistos) == 0: continue yoda.writeYODA(outhistos,"LowEnergy-EE-%s-%s.yoda" % (runType,args[0]))