Provides access to the PYTHIA hadronization models.
Is a concrete implementation of the HadronizationModelI interface. More...

#include <PythiaHadronization.h>

Inheritance diagram for genie::PythiaHadronization:

Public Member Functions
	PythiaHadronization ()

	PythiaHadronization (string config)

virtual	~PythiaHadronization ()

void	Initialize (void) const

TClonesArray *	Hadronize (const Interaction *) const

double	Weight (void) const

PDGCodeList *	SelectParticles (const Interaction *) const

TH1D *	MultiplicityProb (const Interaction , Option_t opt="") const

void	Configure (const Registry &config)

void	Configure (string config)

Public Member Functions inherited from genie::HadronizationModelI
virtual	~HadronizationModelI ()

Public Member Functions inherited from genie::Algorithm
virtual	~Algorithm ()

virtual void	FindConfig (void)

virtual const Registry &	GetConfig (void) const

Registry *	GetOwnedConfig (void)

virtual const AlgId &	Id (void) const
	Get algorithm ID. More...

virtual AlgStatus_t	GetStatus (void) const
	Get algorithm status. More...

virtual bool	AllowReconfig (void) const

virtual AlgCmp_t	Compare (const Algorithm *alg) const
	Compare with input algorithm. More...

virtual void	SetId (const AlgId &id)
	Set algorithm ID. More...

virtual void	SetId (string name, string config)

const Algorithm *	SubAlg (const RgKey &registry_key) const

void	AdoptConfig (void)

void	AdoptSubstructure (void)

virtual void	Print (ostream &stream) const
	Print algorithm info. More...

Private Member Functions
void	LoadConfig (void)

bool	AssertValidity (const Interaction *i) const

Private Attributes
TPythia6 *	fPythia
	PYTHIA6 wrapper class. More...

const DecayModelI *	fDecayer

double	fSSBarSuppression
	ssbar suppression More...

double	fGaussianPt2
	gaussian pt2 distribution width More...

double	fNonGaussianPt2Tail
	non gaussian pt2 tail parameterization More...

double	fRemainingECutoff
	remaining E cutoff for stopping fragmentation More...

double	fDiQuarkSuppression
	di-quark suppression parameter More...

double	fLightVMesonSuppression
	Light vector meon suppression. More...

double	fSVMesonSuppression
	Strange vector meson suppression. More...

double	fLunda
	Lund a parameter. More...

double	fLundb
	Lund b parameter. More...

double	fLundaDiq
	adjustment of Lund a for di-quark More...

Additional Inherited Members
Protected Member Functions inherited from genie::HadronizationModelBase
	HadronizationModelBase ()

	HadronizationModelBase (string name)

	HadronizationModelBase (string name, string config)

virtual	~HadronizationModelBase ()

double	Wmin (void) const
	Various utility methods common to hadronization models. More...

double	MaxMult (const Interaction *i) const

void	ApplyRijk (const Interaction i, bool norm, TH1D mp) const

TH1D *	CreateMultProbHist (double maxmult) const

Protected Member Functions inherited from genie::HadronizationModelI
	HadronizationModelI ()

	HadronizationModelI (string name)

	HadronizationModelI (string name, string config)

Protected Member Functions inherited from genie::Algorithm
	Algorithm ()

	Algorithm (string name)

	Algorithm (string name, string config)

void	Initialize (void)

void	DeleteConfig (void)

void	DeleteSubstructure (void)

Registry *	ExtractLocalConfig (const Registry &in) const

Registry *	ExtractLowerConfig (const Registry &in, const string &alg_key) const
	Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key. More...

template<class T >
bool	GetParam (const RgKey &name, T &p, bool is_top_call=true) const

template<class T >
bool	GetParamDef (const RgKey &name, T &p, const T &def) const

template<class T >
bool	GetParamVect (const std::string &comm_name, std::vector< T > &v, unsigned int max, bool is_top_call=true) const

int	AddTopRegistry (Registry *rp, bool owns=true)
	add registry with top priority, also update ownership More...

int	AddLowRegistry (Registry *rp, bool owns=true)
	add registry with lowest priority, also update ownership More...

int	MergeTopRegistry (const Registry &r)

int	AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
	Add registries with top priority, also udated Ownerships. More...

Protected Attributes inherited from genie::HadronizationModelBase
double	fWcut
	configuration data common to all hadronizers More...

double	fRvpCCm2
	neugen's Rijk: vp, CC, multiplicity = 2 More...

double	fRvpCCm3
	neugen's Rijk: vp, CC, multiplicity = 3 More...

double	fRvpNCm2
	neugen's Rijk: vp, NC, multiplicity = 2 More...

double	fRvpNCm3
	neugen's Rijk: vp, NC, multiplicity = 3 More...

double	fRvnCCm2
	neugen's Rijk: vn, CC, multiplicity = 2 More...

double	fRvnCCm3
	neugen's Rijk: vn, CC, multiplicity = 3 More...

double	fRvnNCm2
	neugen's Rijk: vn, NC, multiplicity = 2 More...

double	fRvnNCm3
	neugen's Rijk: vn, NC, multiplicity = 3 More...

double	fRvbpCCm2
	neugen's Rijk: vbp, CC, multiplicity = 2 More...

double	fRvbpCCm3
	neugen's Rijk: vbp, CC, multiplicity = 3 More...

double	fRvbpNCm2
	neugen's Rijk: vbp, NC, multiplicity = 2 More...

double	fRvbpNCm3
	neugen's Rijk: vbp, NC, multiplicity = 3 More...

double	fRvbnCCm2
	neugen's Rijk: vbn, CC, multiplicity = 2 More...

double	fRvbnCCm3
	neugen's Rijk: vbn, CC, multiplicity = 3 More...

double	fRvbnNCm2
	neugen's Rijk: vbn, NC, multiplicity = 2 More...

double	fRvbnNCm3
	neugen's Rijk: vbn, NC, multiplicity = 3 More...

Protected Attributes inherited from genie::Algorithm
bool	fAllowReconfig

bool	fOwnsSubstruc
	true if it owns its substructure (sub-algs,...) More...

AlgId	fID
	algorithm name and configuration set More...

vector< Registry * >	fConfVect

vector< bool >	fOwnerships
	ownership for every registry in fConfVect More...

AlgStatus_t	fStatus
	algorithm execution status More...

AlgMap *	fOwnedSubAlgMp
	local pool for owned sub-algs (taken out of the factory pool) More...

Detailed Description

Provides access to the PYTHIA hadronization models.
Is a concrete implementation of the HadronizationModelI interface.

Author: Costas Andreopoulos <costas.andreopoulos stfc.ac.uk> University of Liverpool & STFC Rutherford Appleton Lab

August 17, 2004

Definition at line 30 of file PythiaHadronization.h.

Constructor & Destructor Documentation

PythiaHadronization::PythiaHadronization ( )

Definition at line 46 of file PythiaHadronization.cxx.

                                          :
 HadronizationModelBase("genie::PythiaHadronization")
 {
   this->Initialize();
 }

PythiaHadronization::PythiaHadronization ( string config )

Definition at line 52 of file PythiaHadronization.cxx.

                                                       :
 HadronizationModelBase("genie::PythiaHadronization", config)
 {
   this->Initialize();
 }

PythiaHadronization::~PythiaHadronization ( )

virtual

Definition at line 58 of file PythiaHadronization.cxx.

 {
 
 }

Member Function Documentation

bool PythiaHadronization::AssertValidity ( const Interaction * i ) const

private

Definition at line 500 of file PythiaHadronization.cxx.

 {
   // check that there is no charm production 
   // (GENIE uses a special model for these cases)
   if(interaction->ExclTag().IsCharmEvent()) {
      LOG("PythiaHad", pWARN) << "Can't hadronize charm events";
      return false;
   }
   // check the available mass
   double W = utils::kinematics::W(interaction);
   if(W < this->Wmin()) {
      LOG("PythiaHad", pWARN) << "Low invariant mass, W = " << W << " GeV!!";
      return false;
   }
 
   const InitialState & init_state = interaction->InitState();
   const ProcessInfo &  proc_info  = interaction->ProcInfo();
   const Target &       target     = init_state.Tgt();
 
   if( ! target.HitQrkIsSet() ) {
      LOG("PythiaHad", pWARN) << "Hit quark was not set!";
      return false;
   }
 
   int  probe       = init_state.ProbePdg();
   int  hit_nucleon = target.HitNucPdg();
   int  hit_quark   = target.HitQrkPdg();
 //bool from_sea    = target.HitSeaQrk();
 
   // check hit-nucleon assignment, input neutrino & weak current
   bool isp  = pdg::IsProton           (hit_nucleon);
   bool isn  = pdg::IsNeutron          (hit_nucleon);
   bool isv  = pdg::IsNeutrino         (probe);
   bool isvb = pdg::IsAntiNeutrino     (probe);
   bool isdm = pdg::IsDarkMatter         (probe);
   bool isl  = pdg::IsNegChargedLepton (probe);
   bool islb = pdg::IsPosChargedLepton (probe);
   bool iscc = proc_info.IsWeakCC      ();
   bool isnc = proc_info.IsWeakNC      ();
   bool isdmi = proc_info.IsDarkMatter  ();
   bool isem = proc_info.IsEM          ();
   if( !(iscc||isnc||isem||isdmi) ) {
     LOG("PythiaHad", pWARN) 
        << "Can only handle electro-weak interactions";
     return false;
   }
   if( !(isp||isn) || !(isv||isvb||isl||islb||isdm) ) {
     LOG("PythiaHad", pWARN) 
       << "Invalid initial state: probe = " 
       << probe << ", hit_nucleon = " << hit_nucleon;
     return false;
   }
 
   // assert that the interaction mode is allowed
   bool isu  = pdg::IsUQuark     (hit_quark);
   bool isd  = pdg::IsDQuark     (hit_quark);
   bool iss  = pdg::IsSQuark     (hit_quark);
   bool isub = pdg::IsAntiUQuark (hit_quark);
   bool isdb = pdg::IsAntiDQuark (hit_quark);
   bool issb = pdg::IsAntiSQuark (hit_quark);
 
   bool allowed = (iscc && isv  && (isd||isub||iss))  ||
                  (iscc && isvb && (isu||isdb||issb)) ||
                  (isnc && (isv||isvb) && (isu||isd||isub||isdb||iss||issb)) ||
                  (isdmi && isdm && (isu||isd||isub||isdb||iss||issb)) ||
                  (isem && (isl||islb) && (isu||isd||isub||isdb||iss||issb));
   if(!allowed) {
     LOG("PythiaHad", pWARN) 
       << "Impossible interaction type / probe / hit quark combination!";
     return false;
   }
 
   return true;
 }

void PythiaHadronization::Configure ( const Registry & config )

virtual

Configure the algorithm with an external registry The registry is merged with the top level registry if it is owned, Otherwise a copy of it is added with the highest priority

Reimplemented from genie::Algorithm.

Definition at line 423 of file PythiaHadronization.cxx.

 {
   Algorithm::Configure(config);
   this->LoadConfig();
 }

void PythiaHadronization::Configure ( string config )

virtual

Configure the algorithm from the AlgoConfigPool based on param_set string given in input An algorithm contains a vector of registries coming from different xml configuration files, which are loaded according a very precise prioriy This methods will load a number registries in order of priority: 1) "Tunable" parameter set from CommonParametes. This is loaded with the highest prioriry and it is designed to be used for tuning procedure Usage not expected from the user. 2) For every string defined in "CommonParame" the corresponding parameter set will be loaded from CommonParameter.xml 3) parameter set specified by the config string and defined in the xml file of the algorithm 4) if config is not "Default" also the Default parameter set from the same xml file will be loaded Effectively this avoids the repetion of a parameter when it is not changed in the requested configuration

Reimplemented from genie::Algorithm.

Definition at line 429 of file PythiaHadronization.cxx.

 {
   Algorithm::Configure(config);
   this->LoadConfig();
 }

TClonesArray * PythiaHadronization::Hadronize ( const Interaction * interaction ) const

virtual

Implements genie::HadronizationModelBase.

Definition at line 72 of file PythiaHadronization.cxx.

 {
   LOG("PythiaHad", pNOTICE) << "Running PYTHIA hadronizer";
 
   if(!this->AssertValidity(interaction)) {
      LOG("PythiaHad", pERROR) << "Returning a null particle list!";
      return 0;
   }
 
   // get kinematics / init-state / process-info
 
   const Kinematics &   kinematics = interaction->Kine();
   const InitialState & init_state = interaction->InitState();
   const ProcessInfo &  proc_info  = interaction->ProcInfo();
   const Target &       target     = init_state.Tgt();
 
   assert(target.HitQrkIsSet()); 
 
   double W = kinematics.W();
 
   int  probe       = init_state.ProbePdg();
   int  hit_nucleon = target.HitNucPdg();
   int  hit_quark   = target.HitQrkPdg();
   bool from_sea    = target.HitSeaQrk();
 
   LOG("PythiaHad", pNOTICE)
           << "Hit nucleon pdgc = " << hit_nucleon << ", W = " << W;
   LOG("PythiaHad", pNOTICE)
             << "Selected hit quark pdgc = " << hit_quark
                            << ((from_sea) ? "[sea]" : "[valence]");
 
   // check hit-nucleon assignment, input neutrino & interaction type
   bool isp  = pdg::IsProton           (hit_nucleon);
   bool isn  = pdg::IsNeutron          (hit_nucleon);
   bool isv  = pdg::IsNeutrino         (probe);
   bool isvb = pdg::IsAntiNeutrino     (probe);
 //bool isl  = pdg::IsNegChargedLepton (probe);
 //bool islb = pdg::IsPosChargedLepton (probe);
   bool iscc = proc_info.IsWeakCC      ();
   bool isnc = proc_info.IsWeakNC      ();
   bool isdm = proc_info.IsDarkMatter  ();
   bool isem = proc_info.IsEM          ();
   bool isu  = pdg::IsUQuark           (hit_quark);
   bool isd  = pdg::IsDQuark           (hit_quark);
   bool iss  = pdg::IsSQuark           (hit_quark);
   bool isub = pdg::IsAntiUQuark       (hit_quark);
   bool isdb = pdg::IsAntiDQuark       (hit_quark);
   bool issb = pdg::IsAntiSQuark       (hit_quark);
 
   //
   // Generate the quark system (q + qq) initiating the hadronization
   //
 
   int  final_quark = 0; // leading quark (hit quark after the interaction)
   int  diquark     = 0; // remnant diquark (xF<0 at hadronic CMS)
 
   // Figure out the what happens to the hit quark after the interaction
   if (isnc || isem || isdm) {
     // NC, EM
     final_quark = hit_quark;
   } else {
     // CC
     if      (isv  && isd ) final_quark = kPdgUQuark;
     else if (isv  && iss ) final_quark = kPdgUQuark;
     else if (isv  && isub) final_quark = kPdgAntiDQuark;
     else if (isvb && isu ) final_quark = kPdgDQuark;
     else if (isvb && isdb) final_quark = kPdgAntiUQuark;
     else if (isvb && issb) final_quark = kPdgAntiUQuark;
     else {
       LOG("PythiaHad", pERROR)
         << "Not allowed mode. Refused to make a final quark assignment!";
       return 0;
     }
   }//CC
 
   // Figure out what the remnant diquark is.
   // Note from Hugh, following a conversation with his local HEP theorist 
   // (Gary Goldstein): "I am told that the probability of finding the diquark 
   // in the singlet vs. triplet states is 50-50."  
 
   // hit quark = valence quark
   if(!from_sea) {
     if (isp && isu) diquark = kPdgUDDiquarkS1; /* u(->q) + ud */
     if (isp && isd) diquark = kPdgUUDiquarkS1; /* d(->q) + uu */
     if (isn && isu) diquark = kPdgDDDiquarkS1; /* u(->q) + dd */
     if (isn && isd) diquark = kPdgUDDiquarkS1; /* d(->q) + ud */
   }
   // hit quark = sea quark
   else {
     if(isp && isu) diquark = kPdgUDDiquarkS1; /* u(->q) + bar{u} uud (=ud) */
     if(isp && isd) diquark = kPdgUUDiquarkS1; /* d(->q) + bar{d} uud (=uu) */
     if(isn && isu) diquark = kPdgDDDiquarkS1; /* u(->q) + bar{u} udd (=dd) */
     if(isn && isd) diquark = kPdgUDDiquarkS1; /* d(->q) + bar{d} udd (=ud) */
 
     // The following section needs revisiting.
 
     // The lepton is scattered off a sea antiquark, materializing its quark
     // partner and leaving me with a 5q system ( <qbar + q> + qqq(valence) )
     // I will force few qbar+q annhilations below to get my quark/diquark system
     // Probably it is best to leave the qqq system in the final state and then
     // just do the fragmentation of the qbar q system? But how do I figure out
     // how to split the available energy?
 
     /* bar{u} (-> bar{d}) + u uud => u + uu */
     if(isp && isub && iscc)         {final_quark = kPdgUQuark; diquark = kPdgUUDiquarkS1;}
     /* bar{u} (-> bar{u}) + u uud => u + ud */
     if(isp && isub && (isnc||isem||isdm)) {final_quark = kPdgUQuark; diquark = kPdgUDDiquarkS1;}
     /* bar{d} (-> bar{u}) + d uud => d + ud */
     if(isp && isdb && iscc)         {final_quark = kPdgDQuark; diquark = kPdgUDDiquarkS1;}
     /* bar{d} (-> bar{d}) + d uud => d + uu */
     if(isp && isdb && (isnc||isem||isdm)) {final_quark = kPdgDQuark; diquark = kPdgUUDiquarkS1;}
     /* bar{u} (-> bar{d}) + u udd => u + ud */
     if(isn && isub && iscc)         {final_quark = kPdgUQuark; diquark = kPdgUDDiquarkS1;}
     /* bar{u} (-> bar{u}) + u udd => u + dd */
     if(isn && isub && (isnc||isem||isdm)) {final_quark = kPdgUQuark; diquark = kPdgDDDiquarkS1;}
     /* bar{d} (-> bar{u}) + d udd => d + dd */
     if(isn && isdb && iscc)         {final_quark = kPdgDQuark; diquark = kPdgDDDiquarkS1;}
     /* bar{d} (-> bar{d}) + d udd => d + ud */
     if(isn && isdb && (isnc||isem||isdm)) {final_quark = kPdgDQuark; diquark = kPdgUDDiquarkS1;}
 
     // The neutrino is scatterred off s or sbar sea quarks 
     // For the time being I will handle s like d and sbar like dbar (copy & paste
     // from above) so that I conserve charge. 
 
     if(iss || issb) {
        LOG("PythiaHad", pNOTICE) 
                  << "Can not really handle a hit s or sbar quark / Faking it";
 
        if(isp && iss) { diquark = kPdgUUDiquarkS1; }
        if(isn && iss) { diquark = kPdgUDDiquarkS1; }
 
        if(isp && issb && iscc)         {final_quark = kPdgDQuark; diquark = kPdgUDDiquarkS1;}
        if(isp && issb && (isnc||isem||isdm)) {final_quark = kPdgDQuark; diquark = kPdgUUDiquarkS1;}
        if(isn && issb && iscc)         {final_quark = kPdgDQuark; diquark = kPdgDDDiquarkS1;}
        if(isn && issb && (isnc||isem||isdm)) {final_quark = kPdgDQuark; diquark = kPdgUDDiquarkS1;}
     }
  
     // if the diquark is a ud, switch it to the singlet state with 50% probability
     if(diquark == kPdgUDDiquarkS1) {
       RandomGen * rnd = RandomGen::Instance();
       double Rqq = rnd->RndHadro().Rndm();
       if(Rqq<0.5) diquark = kPdgUDDiquarkS0;
     }
   }
   assert(diquark!=0);
 
   //
   // PYTHIA -> HADRONIZATION
   //
 
   LOG("PythiaHad", pNOTICE)
         << "Fragmentation / Init System: "
         << "q = " << final_quark << ", qq = " << diquark;
   int ip = 0;
 
   // Determine how jetset treats un-stable particles appearing in hadronization
 
   int pi0_decflag = fPythia->GetMDCY(fPythia->Pycomp(kPdgPi0),              1); 
   int K0_decflag  = fPythia->GetMDCY(fPythia->Pycomp(kPdgK0),               1); 
   int K0b_decflag = fPythia->GetMDCY(fPythia->Pycomp(kPdgAntiK0),           1); 
   int L0_decflag  = fPythia->GetMDCY(fPythia->Pycomp(kPdgLambda),           1); 
   int L0b_decflag = fPythia->GetMDCY(fPythia->Pycomp(kPdgAntiLambda),       1);
   int Dm_decflag  = fPythia->GetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaM),  1); 
   int D0_decflag  = fPythia->GetMDCY(fPythia->Pycomp(kPdgP33m1232_Delta0),  1); 
   int Dp_decflag  = fPythia->GetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaP),  1); 
   int Dpp_decflag = fPythia->GetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaPP), 1); 
 
 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
   LOG("PythiaHad", pDEBUG) << "Original decay flag for pi0           =  " << pi0_decflag;
   LOG("PythiaHad", pDEBUG) << "Original decay flag for K0            =  " << K0_decflag;
   LOG("PythiaHad", pDEBUG) << "Original decay flag for \bar{K0}      =  " << K0b_decflag;
   LOG("PythiaHad", pDEBUG) << "Original decay flag for Lambda        =  " << L0_decflag;
   LOG("PythiaHad", pDEBUG) << "Original decay flag for \bar{Lambda0} =  " << L0b_decflag;
   LOG("PythiaHad", pDEBUG) << "Original decay flag for D-            =  " << Dm_decflag;
   LOG("PythiaHad", pDEBUG) << "Original decay flag for D0            =  " << D0_decflag;
   LOG("PythiaHad", pDEBUG) << "Original decay flag for D+            =  " << Dp_decflag;
   LOG("PythiaHad", pDEBUG) << "Original decay flag for D++           =  " << Dpp_decflag;
 #endif
 
   fPythia->SetMDCY(fPythia->Pycomp(kPdgPi0),               1,0); // don't decay pi0
   fPythia->SetMDCY(fPythia->Pycomp(kPdgK0),                1,0); // don't decay K0
   fPythia->SetMDCY(fPythia->Pycomp(kPdgAntiK0),            1,0); // don't decay \bar{K0}
   fPythia->SetMDCY(fPythia->Pycomp(kPdgLambda),            1,0); // don't decay Lambda0
   fPythia->SetMDCY(fPythia->Pycomp(kPdgAntiLambda),        1,0); // don't decay \bar{Lambda0}
   fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaM),   1,1); // decay Delta-
   fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_Delta0),   1,1); // decay Delta0
   fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaP),   1,1); // decay Delta+
   fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaPP),  1,1); // decay Delta++
 
   // -- hadronize --
   py2ent_(&ip, &final_quark, &diquark, &W); // hadronizer
 
   // restore pythia decay settings so as not to interfere with decayer 
   fPythia->SetMDCY(fPythia->Pycomp(kPdgPi0),              1, pi0_decflag);
   fPythia->SetMDCY(fPythia->Pycomp(kPdgK0),               1, K0_decflag);
   fPythia->SetMDCY(fPythia->Pycomp(kPdgAntiK0),           1, K0b_decflag);
   fPythia->SetMDCY(fPythia->Pycomp(kPdgLambda),           1, L0_decflag); 
   fPythia->SetMDCY(fPythia->Pycomp(kPdgAntiLambda),       1, L0b_decflag); 
   fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaM),  1, Dm_decflag); 
   fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_Delta0),  1, D0_decflag); 
   fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaP),  1, Dp_decflag); 
   fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaPP), 1, Dpp_decflag); 
 
   // get LUJETS record
   fPythia->GetPrimaries();
   TClonesArray * pythia_particles =
        (TClonesArray *) fPythia->ImportParticles("All");
 
   // copy PYTHIA container to a new TClonesArray so as to transfer ownership
   // of the container and of its elements to the calling method
 
   int np = pythia_particles->GetEntries();
   assert(np>0);
   TClonesArray * particle_list = new TClonesArray("TMCParticle", np);
   particle_list->SetOwner(true);
 
   unsigned int i = 0;
   TMCParticle * particle = 0;
   TIter particle_iter(pythia_particles);
 
   while( (particle = (TMCParticle *) particle_iter.Next()) ) {
      LOG("PythiaHad", pDEBUG)
           << "Adding final state particle pdgc = " << particle->GetKF() 
           << " with status = " << particle->GetKS();
 
      if(particle->GetKS() == 1) {
         if( pdg::IsQuark  (particle->GetKF()) || 
             pdg::IsDiQuark(particle->GetKF()) ) {
                 LOG("PythiaHad", pERROR)
                   << "Hadronization failed! Bare quark/di-quarks appear in final state!";
             particle_list->Delete();
             delete particle_list;
             return 0;            
         }
      }
 
      // fix numbering scheme used for mother/daughter assignments
      particle->SetParent     (particle->GetParent()     - 1);
      particle->SetFirstChild (particle->GetFirstChild() - 1);
      particle->SetLastChild  (particle->GetLastChild()  - 1);
 
      // insert the particle in the list
      new ( (*particle_list)[i++] ) TMCParticle(*particle);
   }
 
   utils::fragmrec::Print(particle_list);
   return particle_list;
 }

void PythiaHadronization::Initialize ( void ) const

virtual

Don't implement the HadronizationModelI interface Leave it for the concrete implementations (KNO, Pythia,...)

Implements genie::HadronizationModelBase.

Definition at line 63 of file PythiaHadronization.cxx.

 {
   fPythia = TPythia6::Instance();
 
   // sync GENIE/PYTHIA6 seed number
   RandomGen::Instance();
 }

void PythiaHadronization::LoadConfig ( void )

private

Definition at line 435 of file PythiaHadronization.cxx.

 {
   // the configurable PYTHIA parameters used here are the ones used by NUX 
   // (see A.Rubbia's talk @ NuINT-01)
   // The defaults are the values used by PYTHIA
   // Use the NUX config set to set the tuned values as used in NUX.
 
   GetParam( "PYTHIA-SSBarSuppression", fSSBarSuppression ) ;
   GetParam( "PYTHIA-GaussianPt2",      fGaussianPt2      ) ;
   GetParam( "PYTHIA-NonGaussianPt2Tail", fNonGaussianPt2Tail  ) ;
   GetParam( "PYTHIA-RemainingEnergyCutoff", fRemainingECutoff ) ;
 
   GetParam( "PYTHIA-DiQuarkSuppression", fDiQuarkSuppression ) ;
   GetParam( "PYTHIA-LightVMesonSuppression", fLightVMesonSuppression ) ;
   GetParam( "PYTHIA-SVMesonSuppression", fSVMesonSuppression ) ;
   GetParam( "PYTHIA-Lunda", fLunda ) ;
   GetParam( "PYTHIA-Lundb", fLundb ) ;
   GetParam( "PYTHIA-LundaDiq", fLundaDiq ) ;
 
   fPythia->SetPARJ(1,  fDiQuarkSuppression ) ;
   fPythia->SetPARJ(11, fLightVMesonSuppression ) ;
   fPythia->SetPARJ(12, fSVMesonSuppression ) ;
   fPythia->SetPARJ(41, fLunda ) ;
   fPythia->SetPARJ(42, fLundb ) ;
   fPythia->SetPARJ(45, fLundaDiq ) ;
   
   fPythia->SetPARJ(2,  fSSBarSuppression);
   fPythia->SetPARJ(21, fGaussianPt2);
   fPythia->SetPARJ(23, fNonGaussianPt2Tail);
   fPythia->SetPARJ(33, fRemainingECutoff);
 
   // Load Wcut determining the phase space area where the multiplicity prob.
   // scaling factors would be applied -if requested-
   GetParam( "Wcut", fWcut ) ;
 
   // decayer
   fDecayer = 0;
   if( GetConfig().Exists("Decayer") ) {
      fDecayer = dynamic_cast<const DecayModelI *> (this->SubAlg("Decayer"));
      assert(fDecayer);
   }
 
   // Load NEUGEN multiplicity probability scaling parameters Rijk
    //neutrinos
    GetParam( "DIS-HMultWgt-vp-CC-m2",  fRvpCCm2  ) ;
    GetParam( "DIS-HMultWgt-vp-CC-m3",  fRvpCCm3  ) ;
    GetParam( "DIS-HMultWgt-vp-NC-m2",  fRvpNCm2  ) ;
    GetParam( "DIS-HMultWgt-vp-NC-m3",  fRvpNCm3  ) ;
    GetParam( "DIS-HMultWgt-vn-CC-m2",  fRvnCCm2  ) ;
    GetParam( "DIS-HMultWgt-vn-CC-m3",  fRvnCCm3  ) ;
    GetParam( "DIS-HMultWgt-vn-NC-m2",  fRvnNCm2  ) ;
    GetParam( "DIS-HMultWgt-vn-NC-m3",  fRvnNCm3  ) ;
    //Anti-neutrinos
    GetParam( "DIS-HMultWgt-vbp-CC-m2", fRvbpCCm2 ) ;
    GetParam( "DIS-HMultWgt-vbp-CC-m3", fRvbpCCm3 ) ;
    GetParam( "DIS-HMultWgt-vbp-NC-m2", fRvbpNCm2 ) ;
    GetParam( "DIS-HMultWgt-vbp-NC-m3", fRvbpNCm3 ) ;
    GetParam( "DIS-HMultWgt-vbn-CC-m2", fRvbnCCm2 ) ;
    GetParam( "DIS-HMultWgt-vbn-CC-m3", fRvbnCCm3 ) ;
    GetParam( "DIS-HMultWgt-vbn-NC-m2", fRvbnNCm2 ) ;
    GetParam( "DIS-HMultWgt-vbn-NC-m3", fRvbnNCm3 ) ;
 
   LOG("PythiaHad", pDEBUG) << GetConfig() ;
 }

TH1D * PythiaHadronization::MultiplicityProb	(	const Interaction *	interaction,
		Option_t *	opt = `""`
	)		const

virtual

Implements genie::HadronizationModelBase.

Definition at line 352 of file PythiaHadronization.cxx.

 {
 // Similar comments apply as in SelectParticles()
 
   if(!this->AssertValidity(interaction)) {
      LOG("PythiaHad", pWARN) 
                 << "Returning a null multipicity probability distribution!";
      return 0;
   }
   double maxmult   = this->MaxMult(interaction);
   TH1D * mult_prob = this->CreateMultProbHist(maxmult);
 
   const int nev=500;
   TMCParticle * particle = 0;
 
   for(int iev=0; iev<nev; iev++) {
 
      TClonesArray * particle_list = this->Hadronize(interaction);
      double         weight        = this->Weight();
 
      if(!particle_list) { iev--; continue; }
 
      int n = 0;
      TIter particle_iter(particle_list);
      while ((particle = (TMCParticle *) particle_iter.Next())) 
      {
        if (particle->GetKS()==1) n++;
      }   
      particle_list->Delete();
      delete particle_list;
      mult_prob->Fill( (double)n, weight);
   }
 
   double integral = mult_prob->Integral("width");
   if(integral>0) {
     // Normalize the probability distribution
     mult_prob->Scale(1.0/integral);
   } else {
     SLOG("PythiaHad", pWARN) << "probability distribution integral = 0";
     return mult_prob;
   }
 
   string option(opt);
 
   bool apply_neugen_Rijk = option.find("+LowMultSuppr") != string::npos;
   bool renormalize       = option.find("+Renormalize")  != string::npos;
 
   // Apply the NeuGEN probability scaling factors -if requested-
   if(apply_neugen_Rijk) {
     SLOG("KNOHad", pINFO) << "Applying NeuGEN scaling factors";
      // Only do so for W<Wcut
      const Kinematics & kinematics = interaction->Kine();
      double W = kinematics.W();
      if(W<fWcut) {
        this->ApplyRijk(interaction, renormalize, mult_prob);
      } else {
         SLOG("PythiaHad", pDEBUG)
               << "W = " << W << " < Wcut = " << fWcut
                                 << " - Will not apply scaling factors";
      }//<wcut?
   }//apply?
 
   return mult_prob;
 }

PDGCodeList * PythiaHadronization::SelectParticles ( const Interaction * interaction ) const

virtual

Implements genie::HadronizationModelBase.

Definition at line 323 of file PythiaHadronization.cxx.

 {
 // Works the opposite way (compared with the KNO hadronization model)
 // Rather than having this method as one of the hadronization model components,
 // we extract the list of particles from the fragmentation record after the
 // hadronization has been completed.
 
   TClonesArray * particle_list = this->Hadronize(interaction);
 
   if(!particle_list) return 0;
 
   bool allowdup=true;
   PDGCodeList * pdgcv = new PDGCodeList(allowdup);
   pdgcv->reserve(particle_list->GetEntries());
 
   TMCParticle * particle = 0;
   TIter particle_iter(particle_list);
 
   while ((particle = (TMCParticle *) particle_iter.Next())) 
   {
     if (particle->GetKS()==1) pdgcv->push_back(particle->GetKF());
   }
   particle_list->Delete();
   delete particle_list;
 
   return pdgcv;
 }

double PythiaHadronization::Weight ( void ) const

virtual

Implements genie::HadronizationModelBase.

Definition at line 418 of file PythiaHadronization.cxx.

 {
   return 1.; // does not generate weighted events
 }

Member Data Documentation

const DecayModelI* genie::PythiaHadronization::fDecayer

private

Definition at line 59 of file PythiaHadronization.h.

double genie::PythiaHadronization::fDiQuarkSuppression

private

di-quark suppression parameter

Definition at line 68 of file PythiaHadronization.h.

double genie::PythiaHadronization::fGaussianPt2

private

gaussian pt2 distribution width

Definition at line 65 of file PythiaHadronization.h.

double genie::PythiaHadronization::fLightVMesonSuppression

private

Light vector meon suppression.

Definition at line 69 of file PythiaHadronization.h.

double genie::PythiaHadronization::fLunda

private

Lund a parameter.

Definition at line 71 of file PythiaHadronization.h.

double genie::PythiaHadronization::fLundaDiq

private

adjustment of Lund a for di-quark

Definition at line 73 of file PythiaHadronization.h.

double genie::PythiaHadronization::fLundb

private

Lund b parameter.

Definition at line 72 of file PythiaHadronization.h.

double genie::PythiaHadronization::fNonGaussianPt2Tail

private

non gaussian pt2 tail parameterization

Definition at line 66 of file PythiaHadronization.h.

TPythia6* genie::PythiaHadronization::fPythia

mutableprivate

PYTHIA6 wrapper class.

Definition at line 57 of file PythiaHadronization.h.

double genie::PythiaHadronization::fRemainingECutoff

private

remaining E cutoff for stopping fragmentation

Definition at line 67 of file PythiaHadronization.h.

double genie::PythiaHadronization::fSSBarSuppression

private

ssbar suppression

Definition at line 64 of file PythiaHadronization.h.

double genie::PythiaHadronization::fSVMesonSuppression

private

Strange vector meson suppression.

Definition at line 70 of file PythiaHadronization.h.

The documentation for this class was generated from the following files:

Generator/src/Physics/Hadronization/PythiaHadronization.h
Generator/src/Physics/Hadronization/PythiaHadronization.cxx

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation