diff --git a/Hadronization/SpinHadronizer.cc b/Hadronization/SpinHadronizer.cc --- a/Hadronization/SpinHadronizer.cc +++ b/Hadronization/SpinHadronizer.cc @@ -1,149 +1,209 @@ // -*- 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 "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/RSSpinorWaveFunction.h" #include "Cluster.h" using namespace Herwig; void SpinHadronizer:: handle(EventHandler &, const tPVector & tagged,const Hint & ) { for(const tPPtr & hadron : tagged) { // mesons if(MesonMatcher::Check(hadron->data())) { continue; } // baryons else if(BaryonMatcher::Check(hadron->data())) { baryonSpin(hadron); } 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); - int prim_quark = (abs(baryon->id())/1000)%10; + 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_quark<3) return; + 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,outgoing); SpinorWaveFunction::constructSpinInfo(waves,baryon,outgoing,true); } else if(baryon->dataPtr()->iSpin()==PDT::Spin3Half) { vector waves; RhoDMatrix rho; RSSpinorWaveFunction::calculateWaveFunctions(waves,rho,baryon,outgoing); RSSpinorWaveFunction::constructSpinInfo(waves,baryon,outgoing,true); } // can't handle spin 5/2 > 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 - vector polB(baryon->dataPtr()->iSpin(),1./double(baryon->dataPtr()->iSpin())); 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) { baryon->spinInfo()->rhoMatrix()(0,0) = 0.375*(1.-pol)*omegaHalf_; baryon->spinInfo()->rhoMatrix()(1,1) = 0.5*(1.-pol)-omegaHalf_/6.*(3.-5.*pol); baryon->spinInfo()->rhoMatrix()(2,2) = 0.5*(1.+pol)-omegaHalf_/6.*(3.+5.*pol); baryon->spinInfo()->rhoMatrix()(3,3) = 0.375*(1.+pol)*omegaHalf_; } else 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"; } 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-Psekin parameter", + "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_ + << "must be lower the than maximum flavour " << maxFlav_ + << " in SpinHadronizer::doinit() " + << Exception::runerror; +} + +void SpinHadronizer::dofinish() { + StepHandler::dofinish(); + for(unsigned int ix=0;ix<3;++ix) { + cerr << "Average polarization of " << getParticleData(long(3+ix))->PDGName() << " antiquarks " + << qPol_[ix].first/qPol_[ix].second << "\n"; + cerr << "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,138 +1,181 @@ // -*- 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.) + SpinHadronizer() : omegaHalf_(2./3.), 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); //@} +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_; + + /** + * 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 */