diff --git a/Hadronization/SpinHadronizer.cc b/Hadronization/SpinHadronizer.cc --- a/Hadronization/SpinHadronizer.cc +++ b/Hadronization/SpinHadronizer.cc @@ -1,210 +1,297 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the SpinHadronizer class. // #include "SpinHadronizer.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/PDT/StandardMatchers.h" -# include "Herwig/Utilities/EnumParticles.h" +#include "Herwig/Utilities/EnumParticles.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/RSSpinorWaveFunction.h" +#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h" +#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h" +#include "ThePEG/Helicity/WaveFunction/TensorWaveFunction.h" #include "Cluster.h" using namespace Herwig; +// choose only baryons void SpinHadronizer:: handle(EventHandler &, const tPVector & tagged,const Hint & ) { for(const tPPtr & hadron : tagged) { // mesons if(MesonMatcher::Check(hadron->data())) { - continue; + mesonSpin(hadron); } // baryons else if(BaryonMatcher::Check(hadron->data())) { baryonSpin(hadron); } - else + else continue; } } void SpinHadronizer::baryonSpin(tPPtr baryon) { // check only one parent if(baryon->parents().size()!=1) return; tPPtr parent = baryon->parents()[0]; // and its a cluster if(parent->id()!=ParticleID::Cluster) return; tClusterPtr cluster = dynamic_ptr_cast(parent); unsigned int prim_quark = (abs(baryon->id())/1000)%10; int sign_quark = baryon->id()>0 ? prim_quark : -prim_quark; // only strange, charm and bottom for the moment if(prim_quarkmaxFlav_ ) return; tPPtr quark; for(unsigned int ix=0;ixnumComponents();++ix) { if(cluster->particle(ix)->id()==sign_quark) { quark = cluster->particle(ix); } } if(!quark) return; if(!quark->spinInfo()) return; tcFermionSpinPtr sp(dynamic_ptr_cast(quark->spinInfo())); // decay it sp->decay(); // create the spin info if(baryon->dataPtr()->iSpin()==PDT::Spin1Half) { vector waves; RhoDMatrix rho; SpinorWaveFunction::calculateWaveFunctions(waves,rho,baryon,baryon->id() >0 ? incoming : outgoing); SpinorWaveFunction::constructSpinInfo(waves,baryon,outgoing,true); } else if(baryon->dataPtr()->iSpin()==PDT::Spin3Half) { vector waves; RhoDMatrix rho; RSSpinorWaveFunction::calculateWaveFunctions(waves,rho,baryon,baryon->id() >0 ? incoming : outgoing); RSSpinorWaveFunction::constructSpinInfo(waves,baryon,outgoing,true); } // can't handle spin > 3/2 else { return; } // extract the polarization of the quark double pol = 2.*sp->rhoMatrix()(1,1).real()-1.; if(sign_quark<0) { qPol_[prim_quark-3].first += pol; qPol_[prim_quark-3].second += 1.; } else { qPol_[prim_quark].first += pol; qPol_[prim_quark].second += 1.; } - - // the different options for different spin types const int mult = prim_quark*1000; int bid = abs(baryon->id()); // lambda and Xi spin 1/2 (spin0 diquark) if(bid== mult+122|| bid== mult+132|| bid== mult+232) { baryon->spinInfo()->rhoMatrix()(0,0) = 0.5*(1.-pol); baryon->spinInfo()->rhoMatrix()(1,1) = 0.5*(1.+pol); } // sigma_b, xi' and omega_b spin 1/2 (spin1 diquark) else if(bid== mult+112|| bid== mult+212|| bid== mult+222|| bid== mult+312|| bid== mult+322|| bid== mult+332) { baryon->spinInfo()->rhoMatrix()(0,0) = 0.5*(1.-pol) +pol*omegaHalf_; baryon->spinInfo()->rhoMatrix()(1,1) = 0.5*(1.+pol) -pol*omegaHalf_; } // sigma*, xi* and omegab* spin 3/2 (spin1 diquark) - else if(bid== mult+114|| bid== mult+214|| bid== mult+224|| bid== mult+334) { + else if(bid== mult+114|| bid== mult+214|| bid== mult+224|| bid== mult+334) { baryon->spinInfo()->rhoMatrix()(0,0) = 0.375*(1.-pol)*omegaHalf_; baryon->spinInfo()->rhoMatrix()(1,1) = 0.5*(1.-pol)-omegaHalf_/8.*(3.-5.*pol); baryon->spinInfo()->rhoMatrix()(2,2) = 0.5*(1.+pol)-omegaHalf_/8.*(3.+5.*pol); baryon->spinInfo()->rhoMatrix()(3,3) = 0.375*(1.+pol)*omegaHalf_; } else - return; + return; // generator()->log() << "Baryon: " << *baryon << "\n"; // generator()->log() << "Parent: " << *cluster << "\n"; // generator()->log() << "Quark: " << *quark << "\n"; // generator()->log() << "Rho\n" << sp->rhoMatrix() << "\n"; // generator()->log() << "testing is decayed " << sp->decayed() <<" \n"; // generator()->log() << baryon->spinInfo()->rhoMatrix() << "\n"; } +void SpinHadronizer::mesonSpin(tPPtr meson) { + // check only one parent + if(meson->parents().size()!=1) return; + tPPtr parent = meson->parents()[0]; + // and its a cluster + if(parent->id()!=ParticleID::Cluster) return; + tClusterPtr cluster = dynamic_ptr_cast(parent); + unsigned int prim_quark = (abs(meson->id())/100)%10; + int sign_quark = meson->id()>0 ? prim_quark : -prim_quark; + // only strange, charm and bottom for the moment + if(prim_quarkmaxFlav_ ) return; + tPPtr quark; + for(unsigned int ix=0;ixnumComponents();++ix) { + if(cluster->particle(ix)->id()==sign_quark) { + quark = cluster->particle(ix); + } + } + if(!quark) return; + if(!quark->spinInfo()) return; + tcFermionSpinPtr sp(dynamic_ptr_cast(quark->spinInfo())); + // decay it + sp->decay(); + // create the spin info + if(meson->dataPtr()->iSpin()==PDT::Spin0) { + RhoDMatrix rho; + ScalarWaveFunction::calculateWaveFunctions(rho,meson,meson->id() > 0 ? incoming : outgoing); + ScalarWaveFunction::constructSpinInfo(meson,outgoing,true); + } + else if(meson->dataPtr()->iSpin()==PDT::Spin1) { + vector waves; + RhoDMatrix rho; + VectorWaveFunction::calculateWaveFunctions(waves,rho,meson,meson->id() > 0 ? incoming : outgoing,true); + VectorWaveFunction::constructSpinInfo(waves,meson,outgoing,true,true); + } + else if(meson->dataPtr()->iSpin()==PDT::Spin2) { + vector waves; + RhoDMatrix rho; + TensorWaveFunction::calculateWaveFunctions(waves,rho,meson,meson->id() > 0 ? incoming : outgoing,true); + TensorWaveFunction::constructSpinInfo(waves,meson,outgoing,true,true); + } + else { + return; + } + // extract the polarization of the quark + double pol = 2.*sp->rhoMatrix()(1,1).real()-1.; + if(sign_quark<0) { + qPol_[prim_quark-3].first += pol; + qPol_[prim_quark-3].second += 1.; + } + else { + qPol_[prim_quark].first += pol; + qPol_[prim_quark].second += 1.; + } + // the different options for different spin types + int bid = abs(meson->id())%1000; + // light-quark spin 1/2+ -> exited spin 0 heavy meson + if(bid==311 || bid==321 || bid==331 || bid==411 || bid==421 || bid==431 + || bid==511 || bid==521 || bid==531) { + return; + } + // light-quark spin 1/2+ + else if(bid==313 || bid==323 || bid==333 || bid==413 || bid==423 || bid==433 + || bid==513 || bid==523 || bid==533) { + // Falk-Peskin "no-win" theorem for non-excited heavy mesons: + // no polarization information would be find in the non-excited meson + // for the excted mesons + meson->spinInfo()->rhoMatrix()(0,0) = (1.-pol)/16. + (omega3Half_/16.)*(3.-5.*pol); + meson->spinInfo()->rhoMatrix()(1,1) = 0.25*(1.-omega3Half_); + meson->spinInfo()->rhoMatrix()(2,2) = (1.+pol)/16. + (omega3Half_/16.)*(3.+5.*pol); + } + // light-quark spin 3/2+ -> exited spin 2 meson + else if(bid==315 || bid==325 || bid==335 || bid==415 || bid==425 || bid==435 + || bid==515 || bid==525 || bid==535) { + meson->spinInfo()->rhoMatrix()(0,0) = 0.25*(1.-pol)*omega3Half_; + meson->spinInfo()->rhoMatrix()(1,1) = 0.1875*(1.-pol)-0.125*(1.-pol)*omega3Half_; + meson->spinInfo()->rhoMatrix()(2,2) = 0.25*(1.-omega3Half_); + meson->spinInfo()->rhoMatrix()(3,3) = 0.1875*(1.+pol)-0.125*(1.+pol)*omega3Half_; + meson->spinInfo()->rhoMatrix()(4,4) = 0.25*(1.+pol)*omega3Half_; + } + else { + return; + } +} + + IBPtr SpinHadronizer::clone() const { return new_ptr(*this); } IBPtr SpinHadronizer::fullclone() const { return new_ptr(*this); } void SpinHadronizer::persistentOutput(PersistentOStream & os) const { os << omegaHalf_; } void SpinHadronizer::persistentInput(PersistentIStream & is, int) { is >> omegaHalf_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigSpinHadronizer("Herwig::SpinHadronizer", "Herwig.so"); void SpinHadronizer::Init() { static ClassDocumentation documentation ("The SpinHadronizer class implements a simple mode for" " the transfer of spin from quarks to hadrons"); static Parameter interfaceOmegaHalf ("OmegaHalf", "The omega_1/2 Falk-Peskin parameter", &SpinHadronizer::omegaHalf_, 2./3., 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceMinimumFlavour ("MinimumFlavour", "The minimum flavour of quark for which to transfer the polarization", &SpinHadronizer::minFlav_, 3, 3, 5, false, false, Interface::limited); static Parameter interfaceMaximumFlavour ("MaximumFlavour", "The maximum flavour of quark for which to transfer the polarization", &SpinHadronizer::maxFlav_, 5, 3, 5, false, false, Interface::limited); static Switch interfaceDebug ("Debug", "Output info on polarizations each for debugging", &SpinHadronizer::debug_, false, false, false); static SwitchOption interfaceDebugYes (interfaceDebug, "Yes", "Debug", true); static SwitchOption interfaceDebugNo (interfaceDebug, "No", "No info", false); } void SpinHadronizer::doinit() { StepHandler::doinit(); if(minFlav_>maxFlav_) - throw InitException() << "The minimum flavour " << minFlav_ + throw InitException() << "The minimum flavour " << minFlav_ << "must be lower the than maximum flavour " << maxFlav_ - << " in SpinHadronizer::doinit() " + << " in SpinHadronizer::doinit() " << Exception::runerror; } void SpinHadronizer::dofinish() { StepHandler::dofinish(); if(debug_) { for(unsigned int ix=0;ix<3;++ix) { if(qPol_[ix].second!=0) generator()->log() << "Average polarization of " << getParticleData(long(3+ix))->PDGName() << " antiquarks " << qPol_[ix].first/qPol_[ix].second << "\n"; if(qPol_[ix+3].second!=0) generator()->log() << "Average polarization of " << getParticleData(long(3+ix))->PDGName() << " quarks " << qPol_[ix+3].first/qPol_[ix+3].second << "\n"; } } } void SpinHadronizer::doinitrun() { StepHandler::doinitrun(); } diff --git a/Hadronization/SpinHadronizer.h b/Hadronization/SpinHadronizer.h --- a/Hadronization/SpinHadronizer.h +++ b/Hadronization/SpinHadronizer.h @@ -1,181 +1,191 @@ // -*- C++ -*- #ifndef Herwig_SpinHadronizer_H #define Herwig_SpinHadronizer_H // // This is the declaration of the SpinHadronizer class. // #include "ThePEG/Handlers/StepHandler.h" namespace Herwig { using namespace ThePEG; /** * The SpinHadronizer class is designed to be used as a post-hadronization handler to * give a simple model of spin transfer between the perturbative and non-perturbative * stages. * * @see \ref SpinHadronizerInterfaces "The interfaces" * defined for SpinHadronizer. */ class SpinHadronizer: public StepHandler { public: - + /** * The default constructor. */ - SpinHadronizer() : omegaHalf_(2./3.), minFlav_(3), maxFlav_(5), debug_(false), - qPol_(6,make_pair(0.,0.)) + SpinHadronizer() : omegaHalf_(2./3.), omega3Half_(0.2), + minFlav_(3), maxFlav_(5), debug_(false), qPol_(6,make_pair(0.,0.)) {} public: /** @name Virtual functions required by the StepHandler class. */ //@{ /** * The main function called by the EventHandler class to * perform a step. Given the current state of an Event, this function * performs the event generation step and includes the result in a new * Step object int the Event record. * @param eh the EventHandler in charge of the Event generation. * @param tagged if not empty these are the only particles which should * be considered by the StepHandler. * @param hint a Hint object with possible information from previously * performed steps. * @throws Veto if the StepHandler requires the current step to be discarded. * @throws Stop if the generation of the current Event should be stopped * after this call. * @throws Exception if something goes wrong. */ virtual void handle(EventHandler & eh, const tPVector & tagged, const Hint & hint); //@} 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: /** * Functions to calculate the spins */ //@{ /** * Calculate the spin of a baryon */ void baryonSpin(tPPtr baryon); - + + /** + * Calculate the spin of a meson + */ + void mesonSpin(tPPtr meson); + //@} 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(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); /** * Finalize this object. Called in the run phase just after a * run has ended. Used eg. to write out statistics. */ virtual void dofinish(); //@} protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ SpinHadronizer & operator=(const SpinHadronizer &); private: /** * Parameters */ //@{ /** * Falk-Peskin \f$\omega_\frac12\f$ parameter */ double omegaHalf_; /** + * Falk-Peskin \f$\omega_\frac32\f$ parameter + */ + double omega3Half_; + + /** * Minimum quark flavour */ unsigned int minFlav_; - + /** * Maximum quark flavour */ unsigned int maxFlav_; /** * Print out debugging info */ bool debug_; /** * Polarization of the quarks */ vector > qPol_; //@} }; } #endif /* Herwig_SpinHadronizer_H */