Public Member Functions | Private Member Functions | Private Attributes | Friends | List of all members
genie::KNOHadronization Class Reference

A KNO-based hadronization model. More...

#include <KNOHadronization.h>

Inheritance diagram for genie::KNOHadronization:
genie::EventRecordVisitorI genie::Algorithm

Public Member Functions

 KNOHadronization ()
 
 KNOHadronization (string config)
 
virtual ~KNOHadronization ()
 
void ProcessEventRecord (GHepRecord *event) const
 
virtual void Configure (const Registry &config)
 
virtual void Configure (string config)
 
- Public Member Functions inherited from genie::EventRecordVisitorI
virtual ~EventRecordVisitorI ()
 
- Public Member Functions inherited from genie::Algorithm
virtual ~Algorithm ()
 
virtual void FindConfig (void)
 
virtual const RegistryGetConfig (void) const
 
RegistryGetOwnedConfig (void)
 
virtual const AlgIdId (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 AlgorithmSubAlg (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)
 
void Initialize (void) const
 
TClonesArray * Hadronize (const Interaction *) const
 
double Weight (void) const
 
PDGCodeListSelectParticles (const Interaction *) const
 
TH1D * MultiplicityProb (const Interaction *, Option_t *opt="") const
 
bool AssertValidity (const Interaction *i) const
 
PDGCodeListGenerateHadronCodes (int mult, int maxQ, double W) const
 
int GenerateBaryonPdgCode (int mult, int maxQ, double W) const
 
int HadronShowerCharge (const Interaction *) const
 
double KNO (int nu, int nuc, double z) const
 
double AverageChMult (int nu, int nuc, double W) const
 
void HandleDecays (TClonesArray *particle_list) const
 
double ReWeightPt2 (const PDGCodeList &pdgcv) const
 
double MaxMult (const Interaction *i) const
 
TH1D * CreateMultProbHist (double maxmult) const
 
void ApplyRijk (const Interaction *i, bool norm, TH1D *mp) const
 
double Wmin (void) const
 
TClonesArray * DecayMethod1 (double W, const PDGCodeList &pdgv, bool reweight_decays) const
 
TClonesArray * DecayMethod2 (double W, const PDGCodeList &pdgv, bool reweight_decays) const
 
TClonesArray * DecayBackToBack (double W, const PDGCodeList &pdgv) const
 
bool PhaseSpaceDecay (TClonesArray &pl, TLorentzVector &pd, const PDGCodeList &pdgv, int offset=0, bool reweight=false) const
 

Private Attributes

TGenPhaseSpace fPhaseSpaceGenerator
 a phase space generator More...
 
double fWeight
 weight for generated event More...
 
bool fForceNeuGenLimit
 force upper hadronic multiplicity to NeuGEN limit More...
 
bool fUseIsotropic2BDecays
 force isotropic, non-reweighted 2-body decays for consistency with neugen/daikon More...
 
bool fUseBaryonXfPt2Param
 Generate baryon xF,pT2 from experimental parameterization? More...
 
bool fReWeightDecays
 Reweight phase space decays? More...
 
bool fForceDecays
 force decays of unstable hadrons produced? More...
 
bool fForceMinMult
 force minimum multiplicity if (at low W) generated less? More...
 
bool fGenerateWeighted
 generate weighted events? More...
 
double fPhSpRwA
 parameter for phase space decay reweighting More...
 
double fPpi0
 {pi0 pi0 } production probability More...
 
double fPpic
 {pi+ pi- } production probability More...
 
double fPKc
 {K+ K- } production probability More...
 
double fPK0
 {K0 K0bar} production probability More...
 
double fPpi0eta
 {Pi0 eta} production probability More...
 
double fPeta
 {eta eta} production probability More...
 
double fAvp
 offset in average charged hadron multiplicity = f(W) relation for vp More...
 
double fAvn
 offset in average charged hadron multiplicity = f(W) relation for vn More...
 
double fAvbp
 offset in average charged hadron multiplicity = f(W) relation for vbp More...
 
double fAvbn
 offset in average charged hadron multiplicity = f(W) relation for vbn More...
 
double fBvp
 slope in average charged hadron multiplicity = f(W) relation for vp More...
 
double fBvn
 slope in average charged hadron multiplicity = f(W) relation for vn More...
 
double fBvbp
 slope in average charged hadron multiplicity = f(W) relation for vbp More...
 
double fBvbn
 slope in average charged hadron multiplicity = f(W) relation for vbn More...
 
double fAhyperon
 parameter controlling strange baryon production probability via associated production (P=a+b*lnW^2) More...
 
double fBhyperon
 see above More...
 
double fCvp
 Levy function parameter for vp. More...
 
double fCvn
 Levy function parameter for vn. More...
 
double fCvbp
 Levy function parameter for vbp. More...
 
double fCvbn
 Levy function parameter for vbn. More...
 
TF1 * fBaryonXFpdf
 baryon xF PDF More...
 
TF1 * fBaryonPT2pdf
 baryon pT^2 PDF More...
 
double fWcut
 Rijk applied for W<Wcut (see DIS/RES join scheme) More...
 
double fRvpCCm2
 Rijk: vp, CC, multiplicity = 2. More...
 
double fRvpCCm3
 Rijk: vp, CC, multiplicity = 3. More...
 
double fRvpNCm2
 Rijk: vp, NC, multiplicity = 2. More...
 
double fRvpNCm3
 Rijk: vp, NC, multiplicity = 3. More...
 
double fRvnCCm2
 Rijk: vn, CC, multiplicity = 2. More...
 
double fRvnCCm3
 Rijk: vn, CC, multiplicity = 3. More...
 
double fRvnNCm2
 Rijk: vn, NC, multiplicity = 2. More...
 
double fRvnNCm3
 Rijk: vn, NC, multiplicity = 3. More...
 
double fRvbpCCm2
 Rijk: vbp, CC, multiplicity = 2. More...
 
double fRvbpCCm3
 Rijk: vbp, CC, multiplicity = 3. More...
 
double fRvbpNCm2
 Rijk: vbp, NC, multiplicity = 2. More...
 
double fRvbpNCm3
 Rijk: vbp, NC, multiplicity = 3. More...
 
double fRvbnCCm2
 Rijk: vbn, CC, multiplicity = 2. More...
 
double fRvbnCCm3
 Rijk: vbn, CC, multiplicity = 3. More...
 
double fRvbnNCm2
 Rijk: vbn, NC, multiplicity = 2. More...
 
double fRvbnNCm3
 Rijk: vbn, NC, multiplicity = 3. More...
 

Friends

class KNOTunedQPMDISPXSec
 

Additional Inherited Members

- Static Public Member Functions inherited from genie::Algorithm
static string BuildParamVectKey (const std::string &comm_name, unsigned int i)
 
static string BuildParamVectSizeKey (const std::string &comm_name)
 
- Protected Member Functions inherited from genie::EventRecordVisitorI
 EventRecordVisitorI ()
 
 EventRecordVisitorI (string name)
 
 EventRecordVisitorI (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)
 
RegistryExtractLocalConfig (const Registry &in) const
 
RegistryExtractLowerConfig (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 >
int GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
 Handle to load vectors of parameters. More...
 
int GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, 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::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< boolfOwnerships
 ownership for every registry in fConfVect More...
 
AlgStatus_t fStatus
 algorithm execution status More...
 
AlgMapfOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool) More...
 

Detailed Description

A KNO-based hadronization model.

Is a concrete implementation of the EventRecordVisitorI interface.

Author
The main authors of this model are:
      - Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
        University of Liverpool & STFC Rutherford Appleton Lab

      - Hugh Gallagher <gallag@minos.phy.tufts.edu>
        Tufts University

      - Tinjun Yang <tjyang@stanford.edu>
        Stanford University

      This is an improved version of the legacy neugen3 KNO-based model.
      Giles Barr, Geoff Pearce, and Hugh Gallagher were listed as authors
      of the original neugen3 model.

      Strange baryon production was implemented by Keith Hofmann and
      Hugh Gallagher (Tufts)

      Production of etas was added by Ji Liu (W&M)

      Changes required to implement the GENIE Boosted Dark Matter module
      were installed by Josh Berger (Univ. of Wisconsin)

August 17, 2004

Copyright (c) 2003-2019, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org or see $GENIE/LICENSE

Definition at line 57 of file KNOHadronization.h.

Constructor & Destructor Documentation

KNOHadronization::KNOHadronization ( )

Definition at line 67 of file KNOHadronization.cxx.

67  :
68 EventRecordVisitorI("genie::KNOHadronization")
69 {
70  fBaryonXFpdf = 0;
71  fBaryonPT2pdf = 0;
72 //fKNO = 0;
73 }
TF1 * fBaryonPT2pdf
baryon pT^2 PDF
TF1 * fBaryonXFpdf
baryon xF PDF
KNOHadronization::KNOHadronization ( string  config)

Definition at line 75 of file KNOHadronization.cxx.

75  :
76 EventRecordVisitorI("genie::KNOHadronization", config)
77 {
78  fBaryonXFpdf = 0;
79  fBaryonPT2pdf = 0;
80 //fKNO = 0;
81 }
TF1 * fBaryonPT2pdf
baryon pT^2 PDF
static Config * config
Definition: config.cpp:1054
TF1 * fBaryonXFpdf
baryon xF PDF
KNOHadronization::~KNOHadronization ( )
virtual

Definition at line 83 of file KNOHadronization.cxx.

84 {
85  if (fBaryonXFpdf ) delete fBaryonXFpdf;
86  if (fBaryonPT2pdf) delete fBaryonPT2pdf;
87 //if (fKNO ) delete fKNO;
88 }
TF1 * fBaryonPT2pdf
baryon pT^2 PDF
TF1 * fBaryonXFpdf
baryon xF PDF

Member Function Documentation

void KNOHadronization::ApplyRijk ( const Interaction i,
bool  norm,
TH1D *  mp 
) const
private

Definition at line 1570 of file KNOHadronization.cxx.

1572 {
1573  // Apply the NEUGEN multiplicity probability scaling factors
1574  //
1575  if(!mp) return;
1576 
1577  const InitialState & init_state = interaction->InitState();
1578  int probe_pdg = init_state.ProbePdg();
1579  int nuc_pdg = init_state.Tgt().HitNucPdg();
1580 
1581  const ProcessInfo & proc_info = interaction->ProcInfo();
1582  bool is_CC = proc_info.IsWeakCC();
1583  bool is_NC = proc_info.IsWeakNC();
1584  bool is_EM = proc_info.IsEM();
1585  // EDIT
1586  bool is_dm = proc_info.IsDarkMatter();
1587 
1588  //
1589  // get the R2, R3 factors
1590  //
1591 
1592  double R2=1., R3=1.;
1593 
1594  // weak CC or NC case
1595  // EDIT
1596  if(is_CC || is_NC || is_dm) {
1597  bool is_nu = pdg::IsNeutrino (probe_pdg);
1598  bool is_nubar = pdg::IsAntiNeutrino (probe_pdg);
1599  bool is_p = pdg::IsProton (nuc_pdg);
1600  bool is_n = pdg::IsNeutron (nuc_pdg);
1601  bool is_dmi = pdg::IsDarkMatter (probe_pdg); // EDIT
1602  if((is_nu && is_p) || (is_dmi && is_p)) {
1603  R2 = (is_CC) ? fRvpCCm2 : fRvpNCm2;
1604  R3 = (is_CC) ? fRvpCCm3 : fRvpNCm3;
1605  } else
1606  if((is_nu && is_n) || (is_dmi && is_n)) {
1607  R2 = (is_CC) ? fRvnCCm2 : fRvnNCm2;
1608  R3 = (is_CC) ? fRvnCCm3 : fRvnNCm3;
1609  } else
1610  if(is_nubar && is_p) {
1611  R2 = (is_CC) ? fRvbpCCm2 : fRvbpNCm2;
1612  R3 = (is_CC) ? fRvbpCCm3 : fRvbpNCm3;
1613  } else
1614  if(is_nubar && is_n) {
1615  R2 = (is_CC) ? fRvbnCCm2 : fRvbnNCm2;
1616  R3 = (is_CC) ? fRvbnCCm3 : fRvbnNCm3;
1617  } else {
1618  LOG("Hadronization", pERROR)
1619  << "Invalid initial state: " << init_state;
1620  }
1621  }//cc||nc?
1622 
1623  // EM case (apply the NC tuning factors)
1624 
1625  if(is_EM) {
1626  bool is_l = pdg::IsNegChargedLepton (probe_pdg);
1627  bool is_lbar = pdg::IsPosChargedLepton (probe_pdg);
1628  bool is_p = pdg::IsProton (nuc_pdg);
1629  bool is_n = pdg::IsNeutron (nuc_pdg);
1630  if(is_l && is_p) {
1631  R2 = fRvpNCm2;
1632  R3 = fRvpNCm3;
1633  } else
1634  if(is_l && is_n) {
1635  R2 = fRvnNCm2;
1636  R3 = fRvnNCm3;
1637  } else
1638  if(is_lbar && is_p) {
1639  R2 = fRvbpNCm2;
1640  R3 = fRvbpNCm3;
1641  } else
1642  if(is_lbar && is_n) {
1643  R2 = fRvbnNCm2;
1644  R3 = fRvbnNCm3;
1645  } else {
1646  LOG("Hadronization", pERROR)
1647  << "Invalid initial state: " << init_state;
1648  }
1649  }//em?
1650 
1651  //
1652  // Apply to the multiplicity probability distribution
1653  //
1654 
1655  int nbins = mp->GetNbinsX();
1656  for(int i = 1; i <= nbins; i++) {
1657  int n = TMath::Nint( mp->GetBinCenter(i) );
1658 
1659  double R=1;
1660  if (n==2) R=R2;
1661  else if (n==3) R=R3;
1662 
1663  if(n==2 || n==3) {
1664  double P = mp->GetBinContent(i);
1665  double Psc = R*P;
1666  LOG("Hadronization", pDEBUG)
1667  << "n=" << n << "/ Scaling factor R = "
1668  << R << "/ P " << P << " --> " << Psc;
1669  mp->SetBinContent(i, Psc);
1670  }
1671  if(n>3) break;
1672  }
1673 
1674  // renormalize the histogram?
1675  if(norm) {
1676  double histo_norm = mp->Integral("width");
1677  if(histo_norm>0) mp->Scale(1.0/histo_norm);
1678  }
1679 }
double fRvnCCm2
Rijk: vn, CC, multiplicity = 2.
bool IsWeakCC(void) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
#define pERROR
Definition: Messenger.h:60
double fRvbnNCm2
Rijk: vbn, NC, multiplicity = 2.
double fRvbnCCm2
Rijk: vbn, CC, multiplicity = 2.
int HitNucPdg(void) const
Definition: Target.cxx:321
double fRvnNCm2
Rijk: vn, NC, multiplicity = 2.
bool IsDarkMatter(int pdgc)
Definition: PDGUtils.cxx:125
double fRvbpCCm3
Rijk: vbp, CC, multiplicity = 3.
double fRvbpNCm2
Rijk: vbp, NC, multiplicity = 2.
double fRvpCCm2
Rijk: vp, CC, multiplicity = 2.
double fRvbnCCm3
Rijk: vbn, CC, multiplicity = 3.
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:310
bool IsPosChargedLepton(int pdgc)
Definition: PDGUtils.cxx:146
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:305
bool IsWeakNC(void) const
double fRvbpCCm2
Rijk: vbp, CC, multiplicity = 2.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:43
std::void_t< T > n
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:116
int ProbePdg(void) const
Definition: InitialState.h:65
double fRvbnNCm3
Rijk: vbn, NC, multiplicity = 3.
double fRvbpNCm3
Rijk: vbp, NC, multiplicity = 3.
auto norm(Vector const &v)
Return norm of the specified vector.
bool IsEM(void) const
bool IsDarkMatter(void) const
double fRvpNCm3
Rijk: vp, NC, multiplicity = 3.
double fRvnCCm3
Rijk: vn, CC, multiplicity = 3.
double fRvnNCm3
Rijk: vn, NC, multiplicity = 3.
double fRvpCCm3
Rijk: vp, CC, multiplicity = 3.
const Target & Tgt(void) const
Definition: InitialState.h:67
bool IsNegChargedLepton(int pdgc)
Definition: PDGUtils.cxx:137
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
double fRvpNCm2
Rijk: vp, NC, multiplicity = 2.
bool KNOHadronization::AssertValidity ( const Interaction i) const
private

Definition at line 1536 of file KNOHadronization.cxx.

1537 {
1538  if(interaction->ExclTag().IsCharmEvent()) {
1539  LOG("KNOHad", pWARN) << "Can't hadronize charm events";
1540  return false;
1541  }
1543  if(W < this->Wmin()) {
1544  LOG("KNOHad", pWARN) << "Low invariant mass, W = " << W << " GeV!!";
1545  return false;
1546  }
1547  return true;
1548 }
double W(const Interaction *const i)
Definition: KineUtils.cxx:1031
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
double Wmin(void) const
#define pWARN
Definition: Messenger.h:61
double KNOHadronization::AverageChMult ( int  nu,
int  nuc,
double  W 
) const
private

Definition at line 666 of file KNOHadronization.cxx.

668 {
669 // computes the average charged multiplicity
670 //
671  bool is_p = pdg::IsProton (nuc_pdg);
672  bool is_n = pdg::IsNeutron (nuc_pdg);
673  bool is_nu = pdg::IsNeutrino (probe_pdg);
674  bool is_nubar = pdg::IsAntiNeutrino (probe_pdg);
675  bool is_l = pdg::IsNegChargedLepton (probe_pdg);
676  bool is_lbar = pdg::IsPosChargedLepton (probe_pdg);
677  // EDIT
678  bool is_dm = pdg::IsDarkMatter (probe_pdg);
679 
680  double a=0, b=0; // params controlling average multiplicity
681 
682  if ( is_p && (is_nu || is_l ) ) { a=fAvp; b=fBvp; }
683  else if ( is_n && (is_nu || is_l ) ) { a=fAvn; b=fBvn; }
684  else if ( is_p && (is_nubar || is_lbar) ) { a=fAvbp; b=fBvbp; }
685  else if ( is_n && (is_nubar || is_lbar) ) { a=fAvbn; b=fBvbn; }
686  // EDIT: assume it's neutrino-like for now...
687  else if ( is_p && is_dm ) { a=fAvp; b=fBvp; }
688  else if ( is_n && is_dm ) { a=fAvn; b=fBvn; }
689  else {
690  LOG("KNOHad", pERROR)
691  << "Invalid initial state (probe = " << probe_pdg << ", "
692  << "hit nucleon = " << nuc_pdg << ")";
693  return 0;
694  }
695 
696  double av_nch = a + b * 2*TMath::Log(W);
697  return av_nch;
698 }
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
#define pERROR
Definition: Messenger.h:60
double fAvn
offset in average charged hadron multiplicity = f(W) relation for vn
bool IsDarkMatter(int pdgc)
Definition: PDGUtils.cxx:125
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:310
bool IsPosChargedLepton(int pdgc)
Definition: PDGUtils.cxx:146
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:305
double fBvbp
slope in average charged hadron multiplicity = f(W) relation for vbp
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
double fAvp
offset in average charged hadron multiplicity = f(W) relation for vp
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:116
double fAvbn
offset in average charged hadron multiplicity = f(W) relation for vbn
const GenericPointer< typename T::ValueType > T2 T::AllocatorType & a
Definition: pointer.h:1124
double fBvn
slope in average charged hadron multiplicity = f(W) relation for vn
double fBvp
slope in average charged hadron multiplicity = f(W) relation for vp
double fBvbn
slope in average charged hadron multiplicity = f(W) relation for vbn
static bool * b
Definition: config.cpp:1043
bool IsNegChargedLepton(int pdgc)
Definition: PDGUtils.cxx:137
double fAvbp
offset in average charged hadron multiplicity = f(W) relation for vbp
void KNOHadronization::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 490 of file KNOHadronization.cxx.

491 {
493  this->LoadConfig();
494 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:71
void KNOHadronization::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 496 of file KNOHadronization.cxx.

497 {
499  this->LoadConfig();
500 }
static Config * config
Definition: config.cpp:1054
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:71
TH1D * KNOHadronization::CreateMultProbHist ( double  maxmult) const
private

Definition at line 1558 of file KNOHadronization.cxx.

1559 {
1560  double minmult = 2;
1561  int nbins = TMath::Nint(maxmult-minmult+1);
1562 
1563  TH1D * mult_prob = new TH1D("mult_prob",
1564  "hadronic multiplicity distribution", nbins, minmult-0.5, maxmult+0.5);
1565  mult_prob->SetDirectory(0);
1566 
1567  return mult_prob;
1568 }
TClonesArray * KNOHadronization::DecayBackToBack ( double  W,
const PDGCodeList pdgv 
) const
private

Definition at line 871 of file KNOHadronization.cxx.

873 {
874 // Handles a special case (only two particles) of the 2nd decay method
875 //
876 
877  LOG("KNOHad", pINFO) << "Generating two particles back-to-back";
878 
879  assert(pdgv.size()==2);
880 
881  RandomGen * rnd = RandomGen::Instance();
882 
883  // Create the particle list
884  TClonesArray * plist = new TClonesArray("genie::GHepParticle", pdgv.size());
885 
886  // Get xF,pT2 distribution (y-) maxima for the rejection method
887  double xFo = 1.1 * fBaryonXFpdf ->GetMaximum(-1,1);
888  double pT2o = 1.1 * fBaryonPT2pdf->GetMaximum( 0,1);
889 
890  TLorentzVector p4(0,0,0,W); // 2-body hadronic system 4p
891 
892  // Do the 2-body decay
893  bool accepted = false;
894  while(!accepted) {
895 
896  // Find an allowed (unweighted) phase space decay for the 2 particles
897  // and add them to the list
898  bool ok = this->PhaseSpaceDecay(*plist, p4, pdgv, 0, false);
899 
900  // If the decay isn't allowed clean-up and return NULL
901  if(!ok) {
902  LOG("KNOHad", pERROR) << "*** Decay forbidden by kinematics! ***";
903  plist->Delete();
904  delete plist;
905  return 0;
906  }
907 
908  // If the decay was allowed, then compute the baryon xF,pT2 and accept/
909  // reject the phase space decays so as to reproduce the xF,pT2 PDFs
910 
911  GHepParticle * baryon = (GHepParticle *) (*plist)[0];
912  assert(pdg::IsNeutronOrProton(baryon->Pdg()));
913 
914  double px = baryon->Px();
915  double py = baryon->Py();
916  double pz = baryon->Pz();
917 
918  double pT2 = px*px + py*py;
919  double pL = pz;
920  double xF = pL/(W/2);
921 
922  double pT2rnd = pT2o * rnd->RndHadro().Rndm();
923  double xFrnd = xFo * rnd->RndHadro().Rndm();
924 
925  double pT2pdf = fBaryonPT2pdf->Eval(pT2);
926  double xFpdf = fBaryonXFpdf ->Eval(xF );
927 
928  LOG("KNOHad", pINFO) << "baryon xF = " << xF << ", pT2 = " << pT2;
929 
930  accepted = (xFrnd < xFpdf && pT2rnd < pT2pdf);
931 
932  LOG("KNOHad", pINFO) << ((accepted) ? "Decay accepted":"Decay rejected");
933  }
934  return plist;
935 }
#define pERROR
Definition: Messenger.h:60
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
TF1 * fBaryonPT2pdf
baryon pT^2 PDF
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:91
double Px(void) const
Get Px.
Definition: GHepParticle.h:89
int Pdg(void) const
Definition: GHepParticle.h:64
bool PhaseSpaceDecay(TClonesArray &pl, TLorentzVector &pd, const PDGCodeList &pdgv, int offset=0, bool reweight=false) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
size_t size
Definition: lodepng.cpp:55
#define pINFO
Definition: Messenger.h:63
bool IsNeutronOrProton(int pdgc)
Definition: PDGUtils.cxx:320
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
TF1 * fBaryonXFpdf
baryon xF PDF
double Py(void) const
Get Py.
Definition: GHepParticle.h:90
TClonesArray * KNOHadronization::DecayMethod1 ( double  W,
const PDGCodeList pdgv,
bool  reweight_decays 
) const
private

Definition at line 732 of file KNOHadronization.cxx.

734 {
735 // Simple phase space decay including all generated particles.
736 // The old NeuGEN decay strategy.
737 
738  LOG("KNOHad", pINFO) << "** Using Hadronic System Decay method 1";
739 
740  TLorentzVector p4had(0,0,0,W);
741  TClonesArray * plist = new TClonesArray("genie::GHepParticle", pdgv.size());
742 
743  // do the decay
744  bool ok = this->PhaseSpaceDecay(*plist, p4had, pdgv, 0, reweight_decays);
745 
746  // clean-up and return NULL
747  if(!ok) {
748  plist->Delete();
749  delete plist;
750  return 0;
751  }
752  return plist;
753 }
bool PhaseSpaceDecay(TClonesArray &pl, TLorentzVector &pd, const PDGCodeList &pdgv, int offset=0, bool reweight=false) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
size_t size
Definition: lodepng.cpp:55
#define pINFO
Definition: Messenger.h:63
TClonesArray * KNOHadronization::DecayMethod2 ( double  W,
const PDGCodeList pdgv,
bool  reweight_decays 
) const
private

Definition at line 755 of file KNOHadronization.cxx.

757 {
758 // Generate the baryon based on experimental pT^2 and xF distributions
759 // Then pass the remaining system of N-1 particles to a phase space decayer.
760 // The strategy adopted at the July-2006 hadronization model mini-workshop.
761 
762  LOG("KNOHad", pINFO) << "** Using Hadronic System Decay method 2";
763 
764  // If only 2 particles are input then don't call the phase space decayer
765  if(pdgv.size() == 2) return this->DecayBackToBack(W,pdgv);
766 
767  // Now handle the more general case:
768 
769  // Take the baryon
770  int baryon = pdgv[0];
771  double MN = PDGLibrary::Instance()->Find(baryon)->Mass();
772  double MN2 = TMath::Power(MN, 2);
773 
774  // Check baryon code
775  // ...
776 
777  // Strip the PDG list from the baryon
778  bool allowdup = true;
779  PDGCodeList pdgv_strip(pdgv.size()-1, allowdup);
780  for(unsigned int i=1; i<pdgv.size(); i++) pdgv_strip[i-1] = pdgv[i];
781 
782  // Get the sum of all masses for the particles in the stripped list
783  double mass_sum = 0;
784  vector<int>::const_iterator pdg_iter = pdgv_strip.begin();
785  for( ; pdg_iter != pdgv_strip.end(); ++pdg_iter) {
786  int pdgc = *pdg_iter;
787  mass_sum += PDGLibrary::Instance()->Find(pdgc)->Mass();
788  }
789 
790  // Create the particle list
791  TClonesArray * plist = new TClonesArray("genie::GHepParticle", pdgv.size());
792 
793  RandomGen * rnd = RandomGen::Instance();
794  TLorentzVector p4had(0,0,0,W);
795  TLorentzVector p4N (0,0,0,0);
796  TLorentzVector p4d;
797 
798  // generate the N 4-p independently
799 
800  bool got_baryon_4p = false;
801  bool got_hadsyst_4p = false;
802 
803  while(!got_hadsyst_4p) {
804 
805  LOG("KNOHad", pINFO) << "Generating p4 for baryon with pdg = " << baryon;
806 
807  while(!got_baryon_4p) {
808 
809  //-- generate baryon xF and pT2
810  double xf = fBaryonXFpdf ->GetRandom();
811  double pt2 = fBaryonPT2pdf->GetRandom();
812 
813  //-- generate baryon px,py,pz
814  double pt = TMath::Sqrt(pt2);
815  double phi = (2*kPi) * rnd->RndHadro().Rndm();
816  double px = pt * TMath::Cos(phi);
817  double py = pt * TMath::Sin(phi);
818  double pz = xf*W/2;
819  double p2 = TMath::Power(pz,2) + pt2;
820  double E = TMath::Sqrt(p2+MN2);
821 
822  p4N.SetPxPyPzE(px,py,pz,E);
823 
824  LOG("KNOHad", pDEBUG) << "Trying nucleon xF= "<< xf<< ", pT2= "<< pt2;
825 
826  //-- check whether there is phase space for the remnant N-1 system
827  p4d = p4had-p4N; // 4-momentum vector for phase space decayer
828  double Mav = p4d.Mag();
829 
830  got_baryon_4p = (Mav > mass_sum);
831 
832  } // baryon xf,pt2 seletion
833 
834  LOG("KNOHad", pINFO)
835  << "Generated baryon with P4 = " << utils::print::P4AsString(&p4N);
836 
837  // Insert the baryon at the event record
838  new ((*plist)[0]) GHepParticle(
839  baryon,kIStStableFinalState, -1,-1,-1,-1,
840  p4N.Px(),p4N.Py(),p4N.Pz(),p4N.Energy(), 0,0,0,0
841  );
842 
843  // Do a phase space decay for the N-1 particles and add them to the list
844  LOG("KNOHad", pINFO)
845  << "Generating p4 for the remaining hadronic system";
846  LOG("KNOHad", pINFO)
847  << "Remaining system: Available mass = " << p4d.Mag()
848  << ", Particle masses = " << mass_sum;
849 
850  bool is_ok = this->PhaseSpaceDecay(
851  *plist, p4d, pdgv_strip, 1, reweight_decays);
852 
853  got_hadsyst_4p = is_ok;
854 
855  if(!got_hadsyst_4p) {
856  got_baryon_4p = false;
857  plist->Delete();
858  }
859  }
860 
861  // clean-up and return NULL
862  if(0) {
863  LOG("KNOHad", pERROR) << "*** Decay forbidden by kinematics! ***";
864  plist->Delete();
865  delete plist;
866  return 0;
867  }
868  return plist;
869 }
#define pERROR
Definition: Messenger.h:60
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:34
TF1 * fBaryonPT2pdf
baryon pT^2 PDF
intermediate_table::const_iterator const_iterator
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
A list of PDG codes.
Definition: PDGCodeList.h:33
bool PhaseSpaceDecay(TClonesArray &pl, TLorentzVector &pd, const PDGCodeList &pdgv, int offset=0, bool reweight=false) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
size_t size
Definition: lodepng.cpp:55
#define pINFO
Definition: Messenger.h:63
TClonesArray * DecayBackToBack(double W, const PDGCodeList &pdgv) const
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
static const double kPi
Definition: Constants.h:38
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
TF1 * fBaryonXFpdf
baryon xF PDF
#define pDEBUG
Definition: Messenger.h:64
int KNOHadronization::GenerateBaryonPdgCode ( int  mult,
int  maxQ,
double  W 
) const
private

Definition at line 1383 of file KNOHadronization.cxx.

1385 {
1386 // Selection of main target fragment (identical as in NeuGEN).
1387 // Assign baryon as p or n. Force it for ++ and - I=3/2 at mult. = 2
1388 
1389  RandomGen * rnd = RandomGen::Instance();
1390  double x = rnd->RndHadro().Rndm();
1391  double y = rnd->RndHadro().Rndm();
1392 
1393  // initialize to neutron & then change it to proton if you must
1394  int pdgc = kPdgNeutron;
1395 
1396  // Assign a probability for the given W for the baryon to become strange
1397  // using a function derived from a fit to the data in Jones et al. (1993)
1398  // Don't let the probability be larger than 1.
1399  double Pstr = fAhyperon + fBhyperon * TMath::Log(W*W);
1400  Pstr = TMath::Min(1.,Pstr);
1401  Pstr = TMath::Max(0.,Pstr);
1402 
1403  // Available hadronic system charge = 2
1404  if(maxQ == 2) {
1405  //for multiplicity ==2, force it to p
1406  if(multiplicity ==2 ) pdgc = kPdgProton;
1407  else {
1408  if(x < 0.66667) pdgc = kPdgProton;
1409  }
1410  }
1411  // Available hadronic system charge = 1
1412  if(maxQ == 1) {
1413  if(multiplicity == 2) {
1414  if(x < 0.33333) pdgc = kPdgProton;
1415  } else {
1416  if(x < 0.50000) pdgc = kPdgProton;
1417  }
1418  }
1419 
1420  // Available hadronic system charge = 0
1421  if(maxQ == 0) {
1422  if(multiplicity == 2) {
1423  if(x < 0.66667) pdgc = kPdgProton;
1424  } else {
1425  if(x < 0.50000) pdgc = kPdgProton;
1426  }
1427  }
1428  // Available hadronic system charge = -1
1429  if(maxQ == -1) {
1430  // for multiplicity == 2, force it to n
1431  if(multiplicity != 2) {
1432  if(x < 0.33333) pdgc = kPdgProton;
1433  }
1434  }
1435 
1436  // For neutrino interactions turn protons and neutrons to Sigma+ and
1437  // Lambda respectively (Lambda and Sigma- respectively for anti-neutrino
1438  // interactions).
1439  if(pdgc == kPdgProton && y < Pstr && maxQ > 0) {
1440  pdgc = kPdgSigmaP;
1441  }
1442  else if(pdgc == kPdgProton && y < Pstr && maxQ <= 0) {
1443  pdgc = kPdgLambda;
1444  }
1445  else if(pdgc == kPdgNeutron && y < Pstr && maxQ > 0) {
1446  pdgc = kPdgLambda;
1447  }
1448  else if(pdgc == kPdgNeutron && y < Pstr && maxQ <= 0) {
1449  pdgc = kPdgSigmaM;
1450  }
1451 
1452  if(pdgc == kPdgProton)
1453  LOG("KNOHad", pDEBUG) << " -> Adding a proton";
1454  if(pdgc == kPdgNeutron)
1455  LOG("KNOHad", pDEBUG) << " -> Adding a neutron";
1456  if(pdgc == kPdgSigmaP)
1457  LOG("KNOHad", pDEBUG) << " -> Adding a sigma+";
1458  if(pdgc == kPdgLambda)
1459  LOG("KNOHad", pDEBUG) << " -> Adding a lambda";
1460  if(pdgc == kPdgSigmaM)
1461  LOG("KNOHad", pDEBUG) << " -> Adding a sigma-";
1462 
1463  return pdgc;
1464 }
const int kPdgLambda
Definition: PDGCodes.h:69
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
double y
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
double fBhyperon
see above
const int kPdgSigmaM
Definition: PDGCodes.h:73
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:54
double fAhyperon
parameter controlling strange baryon production probability via associated production (P=a+b*lnW^2) ...
const int kPdgSigmaP
Definition: PDGCodes.h:71
list x
Definition: train.py:276
const int kPdgProton
Definition: PDGCodes.h:65
const int kPdgNeutron
Definition: PDGCodes.h:67
#define pDEBUG
Definition: Messenger.h:64
PDGCodeList * KNOHadronization::GenerateHadronCodes ( int  mult,
int  maxQ,
double  W 
) const
private

Definition at line 1109 of file KNOHadronization.cxx.

1111 {
1112 // Selection of fragments (identical as in NeuGEN).
1113 
1114  // Get PDG library and rnd num generator
1116  RandomGen * rnd = RandomGen::Instance();
1117 
1118  // Create vector to add final state hadron PDG codes
1119  bool allowdup=true;
1120  PDGCodeList * pdgc = new PDGCodeList(allowdup);
1121  //pdgc->reserve(multiplicity);
1122  int hadrons_to_add = multiplicity;
1123 
1124  //
1125  // Assign baryon as p, n, Sigma+, Sigma- or Lambda
1126  //
1127 
1128  int baryon_code = this->GenerateBaryonPdgCode(multiplicity, maxQ, W);
1129  pdgc->push_back(baryon_code);
1130 
1131  bool baryon_is_strange = (baryon_code == kPdgSigmaP ||
1132  baryon_code == kPdgLambda ||
1133  baryon_code == kPdgSigmaM);
1134  bool baryon_chg_is_pos = (baryon_code == kPdgProton ||
1135  baryon_code == kPdgSigmaP);
1136  bool baryon_chg_is_neg = (baryon_code == kPdgSigmaM);
1137 
1138  // Update number of hadrons to add, available shower charge & invariant mass
1139  if(baryon_chg_is_pos) maxQ -= 1;
1140  if(baryon_chg_is_neg) maxQ += 1;
1141  hadrons_to_add--;
1142  W -= pdg->Find( (*pdgc)[0] )->Mass();
1143 
1144  //
1145  // Assign remaining hadrons up to n = multiplicity
1146  //
1147 
1148  // Conserve strangeness
1149  if(baryon_is_strange) {
1150  LOG("KNOHad", pDEBUG)
1151  << " Remnant baryon is strange. Conserving strangeness...";
1152 
1153  //conserve strangeness and handle charge imbalance with one particle
1154  if(multiplicity == 2) {
1155  if(maxQ == 1) {
1156  LOG("KNOHad", pDEBUG) << " -> Adding a K+";
1157  pdgc->push_back( kPdgKP );
1158 
1159  // update n-of-hadrons to add, avail. shower charge & invariant mass
1160  maxQ -= 1;
1161  hadrons_to_add--;
1162  W -= pdg->Find(kPdgKP)->Mass();
1163  }
1164  else if(maxQ == 0) {
1165  LOG("KNOHad", pDEBUG) << " -> Adding a K0";
1166  pdgc->push_back( kPdgK0 );
1167 
1168  // update n-of-hadrons to add, avail. shower charge & invariant mass
1169  hadrons_to_add--;
1170  W -= pdg->Find(kPdgK0)->Mass();
1171  }
1172  }
1173 
1174  //only two particles left to balance charge
1175  else if(multiplicity == 3 && maxQ == 2) {
1176  LOG("KNOHad", pDEBUG) << " -> Adding a K+";
1177  pdgc->push_back( kPdgKP );
1178 
1179  // update n-of-hadrons to add, avail. shower charge & invariant mass
1180  maxQ -= 1;
1181  hadrons_to_add--;
1182  W -= pdg->Find(kPdgKP)->Mass();
1183  }
1184  else if(multiplicity == 3 && maxQ == -1) { //adding K+ makes it impossible to balance charge
1185  LOG("KNOHad", pDEBUG) << " -> Adding a K0";
1186  pdgc->push_back( kPdgK0 );
1187 
1188  // update n-of-hadrons to add, avail. shower charge & invariant mass
1189  hadrons_to_add--;
1190  W -= pdg->Find(kPdgK0)->Mass();
1191  }
1192 
1193  //simply conserve strangeness, without regard to charge
1194  else {
1195  double y = rnd->RndHadro().Rndm();
1196  if(y < 0.5) {
1197  LOG("KNOHad", pDEBUG) <<" -> Adding a K+";
1198  pdgc->push_back( kPdgKP );
1199 
1200  // update n-of-hadrons to add, avail. shower charge & invariant mass
1201  maxQ -= 1;
1202  hadrons_to_add--;
1203  W -= pdg->Find(kPdgKP)->Mass();
1204  }
1205  else {
1206  LOG("KNOHad", pDEBUG) <<" -> Adding a K0";
1207  pdgc->push_back( kPdgK0 );
1208 
1209  // update n-of-hadrons to add, avail. shower charge & invariant mass
1210  hadrons_to_add--;
1211  W -= pdg->Find(kPdgK0)->Mass();
1212  }
1213  }
1214  }//if the baryon is strange
1215 
1216  // Handle charge imbalance
1217  while(maxQ != 0) {
1218 
1219  if (maxQ < 0) {
1220  // Need more negative charge
1221  LOG("KNOHad", pDEBUG) << "Need more negative charge -> Adding a pi-";
1222  pdgc->push_back( kPdgPiM );
1223 
1224  // update n-of-hadrons to add, avail. shower charge & invariant mass
1225  maxQ += 1;
1226  hadrons_to_add--;
1227 
1228  W -= pdg->Find(kPdgPiM)->Mass();
1229 
1230  } else if (maxQ > 0) {
1231  // Need more positive charge
1232  LOG("KNOHad", pDEBUG) << "Need more positive charge -> Adding a pi+";
1233  pdgc->push_back( kPdgPiP );
1234 
1235  // update n-of-hadrons to add, avail. shower charge & invariant mass
1236  maxQ -= 1;
1237  hadrons_to_add--;
1238 
1239  W -= pdg->Find(kPdgPiP)->Mass();
1240  }
1241  }
1242 
1243  // Add remaining neutrals or pairs up to the generated multiplicity
1244  if(maxQ == 0) {
1245 
1246  LOG("KNOHad", pDEBUG)
1247  << "Hadronic charge balanced. Now adding only neutrals or +- pairs";
1248 
1249  // Final state has correct charge.
1250  // Now add pi0 or pairs (pi0 pi0 / pi+ pi- / K+ K- / K0 K0bar) only
1251 
1252  // Masses of particle pairs
1253  double M2pi0 = 2 * pdg -> Find (kPdgPi0) -> Mass();
1254  double M2pic = pdg -> Find (kPdgPiP) -> Mass() +
1255  pdg -> Find (kPdgPiM) -> Mass();
1256  double M2Kc = pdg -> Find (kPdgKP ) -> Mass() +
1257  pdg -> Find (kPdgKM ) -> Mass();
1258  double M2K0 = 2 * pdg -> Find (kPdgK0 ) -> Mass();
1259  double M2Eta = 2 * pdg -> Find (kPdgEta) -> Mass();
1260  double Mpi0eta = pdg -> Find (kPdgPi0) -> Mass() +
1261  pdg -> Find (kPdgEta) -> Mass();
1262 
1263  // Prevent multiplicity overflow.
1264  // Check if we have an odd number of hadrons to add.
1265  // If yes, add a single pi0 and then go on and add pairs
1266 
1267  if( hadrons_to_add > 0 && hadrons_to_add % 2 == 1 ) {
1268 
1269  LOG("KNOHad", pDEBUG)
1270  << "Odd number of hadrons left to add -> Adding a pi0";
1271  pdgc->push_back( kPdgPi0 );
1272 
1273  // update n-of-hadrons to add & available invariant mass
1274  hadrons_to_add--;
1275  W -= pdg->Find(kPdgPi0)->Mass();
1276  }
1277 
1278  // Now add pairs (pi0 pi0 / pi+ pi- / K+ K- / K0 K0bar)
1279  assert( hadrons_to_add % 2 == 0 ); // even number
1280  LOG("KNOHad", pDEBUG)
1281  <<" hadrons_to_add = "<<hadrons_to_add<<" W= "<<W<<" M2pi0 = "<<M2pi0<<" M2pic = "<<M2pic<<" M2Kc = "<<M2Kc<<" M2K0= "<<M2K0<<" M2Eta= "<<M2Eta;
1282 
1283  while(hadrons_to_add > 0 && W >= M2pi0) {
1284 
1285  double x = rnd->RndHadro().Rndm();
1286  LOG("KNOHad", pDEBUG) << "rndm = " << x;
1287  // Add a pi0 pair
1288  if (x >= 0 && x < fPpi0) {
1289  LOG("KNOHad", pDEBUG) << " -> Adding a pi0pi0 pair";
1290  pdgc->push_back( kPdgPi0 );
1291  pdgc->push_back( kPdgPi0 );
1292  hadrons_to_add -= 2; // update the number of hadrons to add
1293  W -= M2pi0; // update the available invariant mass
1294  }
1295 
1296  // Add a pi+ pi- pair
1297  else if (x < fPpi0 + fPpic) {
1298  if(W >= M2pic) {
1299  LOG("KNOHad", pDEBUG) << " -> Adding a pi+pi- pair";
1300  pdgc->push_back( kPdgPiP );
1301  pdgc->push_back( kPdgPiM );
1302  hadrons_to_add -= 2; // update the number of hadrons to add
1303  W -= M2pic; // update the available invariant mass
1304  } else {
1305  LOG("KNOHad", pDEBUG)
1306  << "Not enough mass for a pi+pi-: trying something else";
1307  }
1308  }
1309 
1310  // Add a K+ K- pair
1311  else if (x < fPpi0 + fPpic + fPKc) {
1312  if(W >= M2Kc) {
1313  LOG("KNOHad", pDEBUG) << " -> Adding a K+K- pair";
1314  pdgc->push_back( kPdgKP );
1315  pdgc->push_back( kPdgKM );
1316  hadrons_to_add -= 2; // update the number of hadrons to add
1317  W -= M2Kc; // update the available invariant mass
1318  } else {
1319  LOG("KNOHad", pDEBUG)
1320  << "Not enough mass for a K+K-: trying something else";
1321  }
1322  }
1323 
1324  // Add a K0 - \bar{K0} pair
1325  else if (x <= fPpi0 + fPpic + fPKc + fPK0) {
1326  if( W >= M2K0 ) {
1327  LOG("KNOHad", pDEBUG) << " -> Adding a K0 K0bar pair";
1328  pdgc->push_back( kPdgK0 );
1329  pdgc->push_back( kPdgAntiK0 );
1330  hadrons_to_add -= 2; // update the number of hadrons to add
1331  W -= M2K0; // update the available invariant mass
1332  } else {
1333  LOG("KNOHad", pDEBUG)
1334  << "Not enough mass for a K0 K0bar: trying something else";
1335  }
1336  }
1337 
1338  // Add a Pi0-Eta pair
1339  else if (x <= fPpi0 + fPpic + fPKc + fPK0 + fPpi0eta) {
1340  if( W >= Mpi0eta ) {
1341  LOG("KNOHad", pDEBUG) << " -> Adding a Pi0-Eta pair";
1342  pdgc->push_back( kPdgPi0 );
1343  pdgc->push_back( kPdgEta );
1344  hadrons_to_add -= 2; // update the number of hadrons to add
1345  W -= Mpi0eta; // update the available invariant mass
1346  } else {
1347  LOG("KNOHad", pDEBUG)
1348  << "Not enough mass for a Pi0-Eta pair: trying something else";
1349  }
1350  }
1351 
1352  //Add a Eta pair
1353  else if(x <= fPpi0 + fPpic + fPKc + fPK0 + fPpi0eta + fPeta) {
1354  if( W >= M2Eta ){
1355  LOG("KNOHad", pDEBUG) << " -> Adding a eta-eta pair";
1356  pdgc->push_back( kPdgEta );
1357  pdgc->push_back( kPdgEta );
1358  hadrons_to_add -= 2; // update the number of hadrons to add
1359  W -= M2Eta; // update the available invariant mass
1360  } else {
1361  LOG("KNOHad", pDEBUG)
1362  << "Not enough mass for a Eta-Eta pair: trying something else";
1363  }
1364 
1365  } else {
1366  LOG("KNOHad", pERROR)
1367  << "Hadron Assignment Probabilities do not add up to 1!!";
1368  exit(1);
1369  }
1370 
1371  // make sure it has enough invariant mass to reach the
1372  // given multiplicity, even by adding only the lightest
1373  // hadron pairs (pi0's)
1374  // Otherwise force a lower multiplicity.
1375  if(W < M2pi0) hadrons_to_add = 0;
1376 
1377  } // while there are more hadrons to add
1378  } // if charge is balanced (maxQ == 0)
1379 
1380  return pdgc;
1381 }
double fPpi0eta
{Pi0 eta} production probability
#define pERROR
Definition: Messenger.h:60
const int kPdgLambda
Definition: PDGCodes.h:69
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
double fPeta
{eta eta} production probability
double Mass(Resonance_t res)
resonance mass (GeV)
double fPpi0
{pi0 pi0 } production probability
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
double fPK0
{K0 K0bar} production probability
A list of PDG codes.
Definition: PDGCodeList.h:33
const int kPdgK0
Definition: PDGCodes.h:155
double y
int GenerateBaryonPdgCode(int mult, int maxQ, double W) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const int kPdgKM
Definition: PDGCodes.h:154
const int kPdgKP
Definition: PDGCodes.h:153
const int kPdgEta
Definition: PDGCodes.h:142
const int kPdgPiP
Definition: PDGCodes.h:139
const int kPdgPi0
Definition: PDGCodes.h:141
double fPKc
{K+ K- } production probability
const int kPdgAntiK0
Definition: PDGCodes.h:156
const int kPdgSigmaM
Definition: PDGCodes.h:73
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:54
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:30
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
const int kPdgPiM
Definition: PDGCodes.h:140
const int kPdgSigmaP
Definition: PDGCodes.h:71
list x
Definition: train.py:276
const int kPdgProton
Definition: PDGCodes.h:65
double fPpic
{pi+ pi- } production probability
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
#define pDEBUG
Definition: Messenger.h:64
TClonesArray * KNOHadronization::Hadronize ( const Interaction interaction) const
private

Definition at line 172 of file KNOHadronization.cxx.

174 {
175 // Generate the hadronic system in a neutrino interaction using a KNO-based
176 // model.
177 
178  if(!this->AssertValidity(interaction)) {
179  LOG("KNOHad", pWARN) << "Returning a null particle list!";
180  return 0;
181  }
182  fWeight=1;
183 
184  double W = utils::kinematics::W(interaction);
185  LOG("KNOHad", pINFO) << "W = " << W << " GeV";
186 
187  //-- Select hadronic shower particles
188  PDGCodeList * pdgcv = this->SelectParticles(interaction);
189 
190  if(!pdgcv) {
191  LOG("KNOHad", pNOTICE)
192  << "Failed selecting particles for " << *interaction;
193  return 0;
194  }
195 
196  //-- Decay the hadronic final state
197  // Two strategies are considered (for N particles):
198  // 1- N (>=2) particles get passed to the phase space decayer. This is the
199  // old NeuGEN strategy.
200  // 2- decay strategy adopted at the July-2006 hadronization model mini-workshop
201  // (C.Andreopoulos, H.Gallagher, P.Kehayias, T.Yang)
202  // The generated baryon P4 gets selected from from experimental xF and pT^2
203  // distributions and the remaining N-1 particles are passed to the phase space
204  // decayer, with P4 = P4(Sum_Hadronic) - P4(Baryon).
205  // For N=2, generate a phase space decay and keep the solution according to its
206  // likelihood calculated based on the baryon xF and pT pdfs. Especially for N=2
207  // keep the option of using simple phase space decay with reweighting switched
208  // off (for consistency with the neugen/daikon version).
209  //
210  TClonesArray * particle_list = 0;
211  bool reweight_decays = fReWeightDecays;
213  bool use_isotropic_decay = (pdgcv->size()==2 && fUseIsotropic2BDecays);
214  if(use_isotropic_decay) {
215  particle_list = this->DecayMethod1(W,*pdgcv,false);
216  } else {
217  particle_list = this->DecayMethod2(W,*pdgcv,reweight_decays);
218  }
219  } else {
220  particle_list = this->DecayMethod1(W,*pdgcv,reweight_decays);
221  }
222 
223  if(!particle_list) {
224  LOG("KNOHad", pNOTICE)
225  << "Failed decaying a hadronic system @ W=" << W
226  << "with multiplicity=" << pdgcv->size();
227 
228  // clean-up and exit
229  delete pdgcv;
230  return 0;
231  }
232 
233  //-- Handle unstable particle decays (if requested)
234  this->HandleDecays(particle_list);
235 
236  //-- The container 'owns' its elements
237  particle_list->SetOwner(true);
238 
239  delete pdgcv;
240 
241  return particle_list;
242 }
bool fReWeightDecays
Reweight phase space decays?
TClonesArray * DecayMethod2(double W, const PDGCodeList &pdgv, bool reweight_decays) const
PDGCodeList * SelectParticles(const Interaction *) const
A list of PDG codes.
Definition: PDGCodeList.h:33
void HandleDecays(TClonesArray *particle_list) const
double W(const Interaction *const i)
Definition: KineUtils.cxx:1031
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
size_t size
Definition: lodepng.cpp:55
#define pINFO
Definition: Messenger.h:63
bool AssertValidity(const Interaction *i) const
double fWeight
weight for generated event
#define pWARN
Definition: Messenger.h:61
bool fUseIsotropic2BDecays
force isotropic, non-reweighted 2-body decays for consistency with neugen/daikon
bool fUseBaryonXfPt2Param
Generate baryon xF,pT2 from experimental parameterization?
TClonesArray * DecayMethod1(double W, const PDGCodeList &pdgv, bool reweight_decays) const
#define pNOTICE
Definition: Messenger.h:62
int KNOHadronization::HadronShowerCharge ( const Interaction interaction) const
private

Definition at line 700 of file KNOHadronization.cxx.

701 {
702 // Returns the hadron shower charge in units of +e
703 // HadronShowerCharge = Q{initial} - Q{final state primary lepton}
704 // eg in v p -> l- X the hadron shower charge is +2
705 // in v n -> l- X the hadron shower charge is +1
706 // in v n -> v X the hadron shower charge is 0
707 //
708  int hadronShowerCharge = 0;
709 
710  // find out the charge of the final state lepton
711  double ql = interaction->FSPrimLepton()->Charge() / 3.;
712 
713  // find out the charge of the probe
714  double qp = interaction->InitState().Probe()->Charge() / 3.;
715 
716  // get the initial state, ask for the hit-nucleon and get
717  // its charge ( = initial state charge for vN interactions)
718  const InitialState & init_state = interaction->InitState();
719  int hit_nucleon = init_state.Tgt().HitNucPdg();
720 
721  assert( pdg::IsProton(hit_nucleon) || pdg::IsNeutron(hit_nucleon) );
722 
723  // Ask PDGLibrary for the nucleon charge
724  double qnuc = PDGLibrary::Instance()->Find(hit_nucleon)->Charge() / 3.;
725 
726  // calculate the hadron shower charge
727  hadronShowerCharge = (int) ( qp + qnuc - ql );
728 
729  return hadronShowerCharge;
730 }
int HitNucPdg(void) const
Definition: Target.cxx:321
TParticlePDG * Probe(void) const
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:310
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:305
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
const InitialState & InitState(void) const
Definition: Interaction.h:68
const Target & Tgt(void) const
Definition: InitialState.h:67
Initial State information.
Definition: InitialState.h:49
void KNOHadronization::HandleDecays ( TClonesArray *  particle_list) const
private

Definition at line 1466 of file KNOHadronization.cxx.

1467 {
1468 // Handle decays of unstable particles if requested through the XML config.
1469 // The default is not to decay the particles at this stage (during event
1470 // generation, the UnstableParticleDecayer event record visitor decays what
1471 // is needed to be decayed later on). But, when comparing various models
1472 // (eg PYTHIA vs KNO) independently and not within the full MC simulation
1473 // framework it might be necessary to force the decays at this point.
1474 
1475  // if (fForceDecays) {
1476  // assert(fDecayer);
1477  //
1478  // //-- loop through the fragmentation event record & decay unstables
1479  // int idecaying = -1; // position of decaying particle
1480  // GHepParticle * p = 0; // current particle
1481  //
1482  // TIter piter(plist);
1483  // while ( (p = (GHepParticle *) piter.Next()) ) {
1484  //
1485  // idecaying++;
1486  // int status = p->Status();
1487  //
1488  // // bother for final state particle only
1489  // if(status < 10) {
1490  //
1491  // // until ROOT's T(MC)Particle(PDG) Lifetime() is fixed, decay only
1492  // // pi^0's
1493  // if ( p->Pdg() == kPdgPi0 ) {
1494  //
1495  // DecayerInputs_t dinp;
1496  //
1497  // TLorentzVector p4;
1498  // p4.SetPxPyPzE(p->Px(), p->Py(), p->Pz(), p->Energy());
1499  //
1500  // dinp.PdgCode = p->Pdg();
1501  // dinp.P4 = &p4;
1502  //
1503  // TClonesArray * decay_products = fDecayer->Decay(dinp);
1504  //
1505  // if(decay_products) {
1506  // //-- mark the parent particle as decayed & set daughters
1507  // p->SetStatus(kIStNucleonTarget);
1508  //
1509  // int nfp = plist->GetEntries(); // n. fragm. products
1510  // int ndp = decay_products->GetEntries(); // n. decay products
1511  //
1512  // p->SetFirstDaughter ( nfp ); // decay products added at
1513  // p->SetLastDaughter ( nfp + ndp -1 ); // the end of the fragm.rec.
1514  //
1515  // //-- add decay products to the fragmentation record
1516  // GHepParticle * dp = 0;
1517  // TIter dpiter(decay_products);
1518  //
1519  // while ( (dp = (GHepParticle *) dpiter.Next()) ) {
1520  //
1521  // dp->SetFirstMother(idecaying);
1522  // new ( (*plist)[plist->GetEntries()] ) GHepParticle(*dp);
1523  // }
1524  //
1525  // //-- clean up decay products
1526  // decay_products->Delete();
1527  // delete decay_products;
1528  // }
1529  //
1530  // } // particle is to be decayed
1531  // } // KS < 10 : final state particle (as in PYTHIA LUJETS record)
1532  // } // particles in fragmentation record
1533  // } // force decay
1534 }
void KNOHadronization::Initialize ( void  ) const
private

Definition at line 92 of file KNOHadronization.cxx.

93 {
94 
95 }
double KNOHadronization::KNO ( int  nu,
int  nuc,
double  z 
) const
private

Definition at line 631 of file KNOHadronization.cxx.

632 {
633 // Computes <n>P(n) for the input reduced multiplicity z=n/<n>
634 
635  bool is_p = pdg::IsProton (nuc_pdg);
636  bool is_n = pdg::IsNeutron (nuc_pdg);
637  bool is_nu = pdg::IsNeutrino (probe_pdg);
638  bool is_nubar = pdg::IsAntiNeutrino (probe_pdg);
639  bool is_l = pdg::IsNegChargedLepton (probe_pdg);
640  bool is_lbar = pdg::IsPosChargedLepton (probe_pdg);
641  // EDIT
642  bool is_dm = pdg::IsDarkMatter (probe_pdg);
643 
644  double c=0; // Levy function parameter
645 
646  if ( is_p && (is_nu || is_l ) ) c=fCvp;
647  else if ( is_n && (is_nu || is_l ) ) c=fCvn;
648  else if ( is_p && (is_nubar || is_lbar) ) c=fCvbp;
649  else if ( is_n && (is_nubar || is_lbar) ) c=fCvbn;
650  // EDIT: assume it's neutrino-like for now...
651  else if ( is_p && is_dm ) c=fCvp;
652  else if ( is_n && is_dm ) c=fCvn;
653  else {
654  LOG("KNOHad", pERROR)
655  << "Invalid initial state (probe = " << probe_pdg << ", "
656  << "hit nucleon = " << nuc_pdg << ")";
657  return 0;
658  }
659 
660  double x = c*z+1;
661  double kno = 2*TMath::Exp(-c)*TMath::Power(c,x)/TMath::Gamma(x);
662 
663  return kno;
664 }
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
#define pERROR
Definition: Messenger.h:60
double fCvp
Levy function parameter for vp.
double fCvbp
Levy function parameter for vbp.
bool IsDarkMatter(int pdgc)
Definition: PDGUtils.cxx:125
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:310
bool IsPosChargedLepton(int pdgc)
Definition: PDGUtils.cxx:146
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:305
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
double fCvbn
Levy function parameter for vbn.
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:116
double z
double fCvn
Levy function parameter for vn.
list x
Definition: train.py:276
bool IsNegChargedLepton(int pdgc)
Definition: PDGUtils.cxx:137
void KNOHadronization::LoadConfig ( void  )
private

Definition at line 504 of file KNOHadronization.cxx.

505 {
506  // Force decays of unstable hadronization products?
507  //GetParamDef( "ForceDecays", fForceDecays, false ) ;
508 
509  // Force minimum multiplicity (if generated less than that) or abort?
510  GetParamDef( "ForceMinMultiplicity", fForceMinMult, true ) ;
511 
512  // Generate the baryon xF and pT^2 using experimental data as PDFs?
513  // In this case, only the N-1 other particles would be fed into the phase
514  // space decayer. This seems to improve hadronic system features such as
515  // bkw/fwd xF hemisphere average multiplicities.
516  // Note: not in the legacy KNO model (NeuGEN). Switch this feature off for
517  // comparisons or for reproducing old simulations.
518  GetParam( "KNO-UseBaryonPdfs-xFpT2", fUseBaryonXfPt2Param ) ;
519 
520  // Reweight the phase space decayer events to reproduce the experimentally
521  // measured pT^2 distributions.
522  // Note: not in the legacy KNO model (NeuGEN). Switch this feature off for
523  // comparisons or for reproducing old simulations.
524  GetParam( "KNO-PhaseSpDec-Reweight", fReWeightDecays ) ;
525 
526  // Parameter for phase space re-weighting. See ReWeightPt2()
527  GetParam( "KNO-PhaseSpDec-ReweightParm", fPhSpRwA ) ;
528 
529  // use isotropic non-reweighted 2-body phase space decays for consistency
530  // with neugen/daikon
531  GetParam( "KNO-UseIsotropic2BodyDec", fUseIsotropic2BDecays ) ;
532 
533  // Generated weighted or un-weighted hadronic systems?
534  GetParamDef( "GenerateWeighted", fGenerateWeighted, false ) ;
535 
536 
537  // Probabilities for producing hadron pairs
538 
539  //-- pi0 pi0
540  GetParam( "KNO-ProbPi0Pi0", fPpi0 ) ;
541  //-- pi+ pi-
542  GetParam( "KNO-ProbPiplusPiminus", fPpic ) ;
543  //-- K+ K-
544  GetParam( "KNO-ProbKplusKminus", fPKc ) ;
545  //-- K0 K0bar
546  GetParam( "KNO-ProbK0K0bar", fPK0 ) ;
547  //-- pi0 eta
548  GetParam( "KNO-ProbPi0Eta", fPpi0eta ) ;
549  //-- eta eta
550  GetParam( "KNO-ProbEtaEta", fPeta ) ;
551 
552  double fsum = fPeta + fPpi0eta + fPK0 + fPKc + fPpic + fPpi0;
553  double diff = TMath::Abs(1.-fsum);
554  if(diff>0.001) {
555  LOG("KNOHad", pWARN) << "KNO Probabilities do not sum to unity! Renormalizing..." ;
556  fPpi0 = fPpi0/fsum;
557  fPpic = fPpic/fsum;
558  fPKc = fPKc/fsum;
559  fPK0 = fPK0/fsum;
560  fPpi0eta = fPpi0eta/fsum;
561  fPeta = fPeta/fsum;
562  }
563 
564 
565  // Baryon pT^2 and xF parameterizations used as PDFs
566 
567  if (fBaryonXFpdf ) delete fBaryonXFpdf;
568  if (fBaryonPT2pdf) delete fBaryonPT2pdf;
569 
570  fBaryonXFpdf = new TF1("fBaryonXFpdf",
571  "0.083*exp(-0.5*pow(x+0.385,2.)/0.131)",-1,0.5);
572  fBaryonPT2pdf = new TF1("fBaryonPT2pdf",
573  "exp(-0.214-6.625*x)",0,0.6);
574  // stop ROOT from deleting these object of its own volition
575  gROOT->GetListOfFunctions()->Remove(fBaryonXFpdf);
576  gROOT->GetListOfFunctions()->Remove(fBaryonPT2pdf);
577 
578 
579  // Load parameters determining the average charged hadron multiplicity
580  GetParam( "KNO-Alpha-vp", fAvp ) ;
581  GetParam( "KNO-Alpha-vn", fAvn ) ;
582  GetParam( "KNO-Alpha-vbp", fAvbp ) ;
583  GetParam( "KNO-Alpha-vbn", fAvbn ) ;
584  GetParam( "KNO-Beta-vp", fBvp ) ;
585  GetParam( "KNO-Beta-vn", fBvn ) ;
586  GetParam( "KNO-Beta-vbp", fBvbp ) ;
587  GetParam( "KNO-Beta-vbn", fBvbn ) ;
588 
589  // Load parameters determining the prob of producing a strange baryon
590  // via associated production
591  GetParam( "KNO-Alpha-Hyperon", fAhyperon ) ;
592  GetParam( "KNO-Beta-Hyperon", fBhyperon ) ;
593 
594  // Load the Levy function parameter
595  GetParam( "KNO-LevyC-vp", fCvp ) ;
596  GetParam( "KNO-LevyC-vn", fCvn ) ;
597  GetParam( "KNO-LevyC-vbp", fCvbp ) ;
598  GetParam( "KNO-LevyC-vbn", fCvbn ) ;
599 
600  // Check whether to generate weighted or unweighted particle decays
601  fGenerateWeighted = false ;
602  //this->GetParam("GenerateWeighted", fGenerateWeighted, false);{
603 
604  // Load Wcut determining the phase space area where the multiplicity prob.
605  // scaling factors would be applied -if requested-
606  this->GetParam( "Wcut", fWcut ) ;
607 
608  // Load NEUGEN multiplicity probability scaling parameters Rijk
609  // neutrinos
610  this->GetParam( "DIS-HMultWgt-vp-CC-m2", fRvpCCm2 ) ;
611  this->GetParam( "DIS-HMultWgt-vp-CC-m3", fRvpCCm3 ) ;
612  this->GetParam( "DIS-HMultWgt-vp-NC-m2", fRvpNCm2 ) ;
613  this->GetParam( "DIS-HMultWgt-vp-NC-m3", fRvpNCm3 ) ;
614  this->GetParam( "DIS-HMultWgt-vn-CC-m2", fRvnCCm2 ) ;
615  this->GetParam( "DIS-HMultWgt-vn-CC-m3", fRvnCCm3 ) ;
616  this->GetParam( "DIS-HMultWgt-vn-NC-m2", fRvnNCm2 ) ;
617  this->GetParam( "DIS-HMultWgt-vn-NC-m3", fRvnNCm3 ) ;
618  //Anti-neutrinos
619  this->GetParam( "DIS-HMultWgt-vbp-CC-m2", fRvbpCCm2 ) ;
620  this->GetParam( "DIS-HMultWgt-vbp-CC-m3", fRvbpCCm3 ) ;
621  this->GetParam( "DIS-HMultWgt-vbp-NC-m2", fRvbpNCm2 ) ;
622  this->GetParam( "DIS-HMultWgt-vbp-NC-m3", fRvbpNCm3 ) ;
623  this->GetParam( "DIS-HMultWgt-vbn-CC-m2", fRvbnCCm2 ) ;
624  this->GetParam( "DIS-HMultWgt-vbn-CC-m3", fRvbnCCm3 ) ;
625  this->GetParam( "DIS-HMultWgt-vbn-NC-m2", fRvbnNCm2 ) ;
626  this->GetParam( "DIS-HMultWgt-vbn-NC-m3", fRvbnNCm3 ) ;
627 
628 
629 }
double fRvnCCm2
Rijk: vn, CC, multiplicity = 2.
double fPpi0eta
{Pi0 eta} production probability
double fRvbnNCm2
Rijk: vbn, NC, multiplicity = 2.
double fRvbnCCm2
Rijk: vbn, CC, multiplicity = 2.
double fWcut
Rijk applied for W<Wcut (see DIS/RES join scheme)
double fAvn
offset in average charged hadron multiplicity = f(W) relation for vn
double fPeta
{eta eta} production probability
double fCvp
Levy function parameter for vp.
double fCvbp
Levy function parameter for vbp.
double fRvnNCm2
Rijk: vn, NC, multiplicity = 2.
bool fReWeightDecays
Reweight phase space decays?
TF1 * fBaryonPT2pdf
baryon pT^2 PDF
double fRvbpCCm3
Rijk: vbp, CC, multiplicity = 3.
double fPpi0
{pi0 pi0 } production probability
double fPK0
{K0 K0bar} production probability
double fRvbpNCm2
Rijk: vbp, NC, multiplicity = 2.
double fRvpCCm2
Rijk: vp, CC, multiplicity = 2.
double fRvbnCCm3
Rijk: vbn, CC, multiplicity = 3.
double fRvbpCCm2
Rijk: vbp, CC, multiplicity = 2.
double fBvbp
slope in average charged hadron multiplicity = f(W) relation for vbp
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
double fAvp
offset in average charged hadron multiplicity = f(W) relation for vp
double fCvbn
Levy function parameter for vbn.
double fBhyperon
see above
bool fGenerateWeighted
generate weighted events?
double fAvbn
offset in average charged hadron multiplicity = f(W) relation for vbn
double fPKc
{K+ K- } production probability
double fRvbnNCm3
Rijk: vbn, NC, multiplicity = 3.
#define pWARN
Definition: Messenger.h:61
double fRvbpNCm3
Rijk: vbp, NC, multiplicity = 3.
double fCvn
Levy function parameter for vn.
double fPhSpRwA
parameter for phase space decay reweighting
bool fUseIsotropic2BDecays
force isotropic, non-reweighted 2-body decays for consistency with neugen/daikon
double fBvn
slope in average charged hadron multiplicity = f(W) relation for vn
double fBvp
slope in average charged hadron multiplicity = f(W) relation for vp
double fAhyperon
parameter controlling strange baryon production probability via associated production (P=a+b*lnW^2) ...
double fRvpNCm3
Rijk: vp, NC, multiplicity = 3.
double fBvbn
slope in average charged hadron multiplicity = f(W) relation for vbn
bool fUseBaryonXfPt2Param
Generate baryon xF,pT2 from experimental parameterization?
double fRvnCCm3
Rijk: vn, CC, multiplicity = 3.
double fRvnNCm3
Rijk: vn, NC, multiplicity = 3.
bool fForceMinMult
force minimum multiplicity if (at low W) generated less?
double fRvpCCm3
Rijk: vp, CC, multiplicity = 3.
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
double fPpic
{pi+ pi- } production probability
TF1 * fBaryonXFpdf
baryon xF PDF
double fRvpNCm2
Rijk: vp, NC, multiplicity = 2.
double fAvbp
offset in average charged hadron multiplicity = f(W) relation for vbp
double KNOHadronization::MaxMult ( const Interaction i) const
private

Definition at line 1550 of file KNOHadronization.cxx.

1551 {
1552  double W = interaction->Kine().W();
1553 
1554  double maxmult = TMath::Floor(1 + (W-kNeutronMass)/kPionMass);
1555  return maxmult;
1556 }
static const double kNeutronMass
Definition: Constants.h:77
static const double kPionMass
Definition: Constants.h:74
TH1D * KNOHadronization::MultiplicityProb ( const Interaction interaction,
Option_t *  opt = "" 
) const
private

Definition at line 372 of file KNOHadronization.cxx.

374 {
375 // Returns a multiplicity probability distribution for the input interaction.
376 // The input option (Default: "") can contain (combinations) of these strings:
377 // - "+LowMultSuppr": applies NeuGEN Rijk factors suppresing the low multipl.
378 // (1-pion and 2-pion) states as part of the DIS/RES joining scheme.
379 // - "+Renormalize": renormalizes the probability distribution after applying
380 // the NeuGEN scaling factors: Eg, when used as a hadronic multiplicity pdf
381 // the output hadronic multiplicity probability histogram needs to be re-
382 // normalized. But, when this method is called from a DIS cross section
383 // algorithm using the integrated probability reduction as a cross section
384 // section reduction factor then the output histogram should not be re-
385 // normalized after applying the scaling factors.
386 
387  if(!this->AssertValidity(interaction)) {
388  LOG("KNOHad", pWARN)
389  << "Returning a null multiplicity probability distribution!";
390  return 0;
391  }
392 
393  const InitialState & init_state = interaction->InitState();
394  int nu_pdg = init_state.ProbePdg();
395  int nuc_pdg = init_state.Tgt().HitNucPdg();
396 
397  // Compute the average charged hadron multiplicity as: <n> = a + b*ln(W^2)
398  // Calculate avergage hadron multiplicity (= 1.5 x charged hadron mult.)
399 
400  double W = utils::kinematics::W(interaction);
401  double avnch = this->AverageChMult(nu_pdg, nuc_pdg, W);
402  double avn = 1.5*avnch;
403 
404  SLOG("KNOHad", pINFO)
405  << "Average hadronic multiplicity (W=" << W << ") = " << avn;
406 
407  // Find the max possible multiplicity as W = Mneutron + (maxmult-1)*Mpion
408  double maxmult = this->MaxMult(interaction);
409 
410  // If required force the NeuGEN maximum multiplicity limit (10)
411  // Note: use for NEUGEN/GENIE comparisons, not physics MC production
412  if(fForceNeuGenLimit && maxmult>10) maxmult=10;
413 
414  // Set maximum multiplicity so that it does not exceed the max number of
415  // particles accepted by the ROOT phase space decayer (18)
416  // Change this if ROOT authors remove the TGenPhaseSpace limitation.
417  if(maxmult>18) maxmult=18;
418 
419  SLOG("KNOHad", pDEBUG) << "Computed maximum multiplicity = " << maxmult;
420 
421  if(maxmult<2) {
422  LOG("KNOHad", pWARN) << "Low maximum multiplicity! Quiting.";
423  return 0;
424  }
425 
426  // Create multiplicity probability histogram
427  TH1D * mult_prob = this->CreateMultProbHist(maxmult);
428 
429  // Compute the multiplicity probabilities values up to the bin corresponding
430  // to the computed maximum multiplicity
431 
432  if(maxmult>2) {
433  int nbins = mult_prob->FindBin(maxmult);
434 
435  for(int i = 1; i <= nbins; i++) {
436  // KNO distribution is <n>*P(n) vs n/<n>
437  double n = mult_prob->GetBinCenter(i); // bin centre
438  double z = n/avn; // z=n/<n>
439  double avnP = this->KNO(nu_pdg,nuc_pdg,z); // <n>*P(n)
440  double P = avnP / avn; // P(n)
441 
442  SLOG("KNOHad", pDEBUG)
443  << "n = " << n << " (n/<n> = " << z
444  << ", <n>*P = " << avnP << ") => P = " << P;
445 
446  mult_prob->Fill(n,P);
447  }
448  } else {
449  SLOG("KNOHad", pDEBUG) << "Fixing multiplicity to 2";
450  mult_prob->Fill(2,1.);
451  }
452 
453  double integral = mult_prob->Integral("width");
454  if(integral>0) {
455  // Normalize the probability distribution
456  mult_prob->Scale(1.0/integral);
457  } else {
458  SLOG("KNOHad", pWARN) << "probability distribution integral = 0";
459  return mult_prob;
460  }
461 
462  string option(opt);
463 
464  bool apply_neugen_Rijk = option.find("+LowMultSuppr") != string::npos;
465  bool renormalize = option.find("+Renormalize") != string::npos;
466 
467  // Apply the NeuGEN probability scaling factors -if requested-
468  if(apply_neugen_Rijk) {
469  SLOG("KNOHad", pINFO) << "Applying NeuGEN scaling factors";
470  // Only do so for W<Wcut
471  if(W<fWcut) {
472  this->ApplyRijk(interaction, renormalize, mult_prob);
473  } else {
474  SLOG("KNOHad", pDEBUG)
475  << "W = " << W << " < Wcut = " << fWcut
476  << " - Will not apply scaling factors";
477  }//<wcut?
478  }//apply?
479 
480  return mult_prob;
481 }
double fWcut
Rijk applied for W<Wcut (see DIS/RES join scheme)
int HitNucPdg(void) const
Definition: Target.cxx:321
opt
Definition: train.py:196
TH1D * CreateMultProbHist(double maxmult) const
bool fForceNeuGenLimit
force upper hadronic multiplicity to NeuGEN limit
double W(const Interaction *const i)
Definition: KineUtils.cxx:1031
double MaxMult(const Interaction *i) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
std::void_t< T > n
int ProbePdg(void) const
Definition: InitialState.h:65
double z
double KNO(int nu, int nuc, double z) const
#define pINFO
Definition: Messenger.h:63
bool AssertValidity(const Interaction *i) const
#define pWARN
Definition: Messenger.h:61
double AverageChMult(int nu, int nuc, double W) const
const InitialState & InitState(void) const
Definition: Interaction.h:68
const Target & Tgt(void) const
Definition: InitialState.h:67
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
void ApplyRijk(const Interaction *i, bool norm, TH1D *mp) const
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
bool KNOHadronization::PhaseSpaceDecay ( TClonesArray &  pl,
TLorentzVector &  pd,
const PDGCodeList pdgv,
int  offset = 0,
bool  reweight = false 
) const
private

Definition at line 937 of file KNOHadronization.cxx.

940 {
941 // General method decaying the input particle system 'pdgv' with available 4-p
942 // given by 'pd'. The decayed system is used to populate the input GHepParticle
943 // array starting from the slot 'offset'.
944 //
945  LOG("KNOHad", pINFO) << "*** Performing a Phase Space Decay";
946  LOG("KNOHad", pINFO) << "pT reweighting is " << (reweight ? "on" : "off");
947 
948  assert ( offset >= 0);
949  assert ( pdgv.size() > 1);
950 
951  // Get the decay product masses
952 
954  int i = 0;
955  double * mass = new double[pdgv.size()];
956  double sum = 0;
957  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
958  int pdgc = *pdg_iter;
959  double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
960  mass[i++] = m;
961  sum += m;
962  }
963 
964  LOG("KNOHad", pINFO)
965  << "Decaying N = " << pdgv.size() << " particles / total mass = " << sum;
966  LOG("KNOHad", pINFO)
967  << "Decaying system p4 = " << utils::print::P4AsString(&pd);
968 
969  // Set the decay
970  bool permitted = fPhaseSpaceGenerator.SetDecay(pd, pdgv.size(), mass);
971  if(!permitted) {
972  LOG("KNOHad", pERROR)
973  << " *** Phase space decay is not permitted \n"
974  << " Total particle mass = " << sum << "\n"
975  << " Decaying system p4 = " << utils::print::P4AsString(&pd);
976 
977  // clean-up and return
978  delete [] mass;
979  return false;
980  }
981 
982  // Get the maximum weight
983  //double wmax = fPhaseSpaceGenerator.GetWtMax();
984  double wmax = -1;
985  for(int idec=0; idec<200; idec++) {
986  double w = fPhaseSpaceGenerator.Generate();
987  if(reweight) { w *= this->ReWeightPt2(pdgv); }
988  wmax = TMath::Max(wmax,w);
989  }
990  assert(wmax>0);
991 
992  LOG("KNOHad", pNOTICE)
993  << "Max phase space gen. weight @ current hadronic system: " << wmax;
994 
995  // Generate a weighted or unweighted decay
996 
997  RandomGen * rnd = RandomGen::Instance();
998 
1000  {
1001  // *** generating weighted decays ***
1002  double w = fPhaseSpaceGenerator.Generate();
1003  if(reweight) { w *= this->ReWeightPt2(pdgv); }
1004  fWeight *= TMath::Max(w/wmax, 1.);
1005  }
1006  else
1007  {
1008  // *** generating un-weighted decays ***
1009  wmax *= 2.3;
1010  bool accept_decay=false;
1011  unsigned int itry=0;
1012 
1013  while(!accept_decay)
1014  {
1015  itry++;
1016 
1017  if(itry>kMaxUnweightDecayIterations) {
1018  // report, clean-up and return
1019  LOG("KNOHad", pWARN)
1020  << "Couldn't generate an unweighted phase space decay after "
1021  << itry << " attempts";
1022  delete [] mass;
1023  return false;
1024  }
1025 
1026  double w = fPhaseSpaceGenerator.Generate();
1027  if(reweight) { w *= this->ReWeightPt2(pdgv); }
1028  if(w > wmax) {
1029  LOG("KNOHad", pWARN)
1030  << "Decay weight = " << w << " > max decay weight = " << wmax;
1031  }
1032  double gw = wmax * rnd->RndHadro().Rndm();
1033  accept_decay = (gw<=w);
1034 
1035  LOG("KNOHad", pINFO)
1036  << "Decay weight = " << w << " / R = " << gw
1037  << " - accepted: " << accept_decay;
1038 
1039  bool return_after_not_accepted_decay = false;
1040  if(return_after_not_accepted_decay && !accept_decay) {
1041  LOG("KNOHad", pWARN)
1042  << "Was instructed to return after a not-accepted decay";
1043  delete [] mass;
1044  return false;
1045  }
1046  }
1047  }
1048 
1049  // Insert final state products into a TClonesArray of GHepParticle's
1050 
1051  i=0;
1052  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
1053 
1054  //-- current PDG code
1055  int pdgc = *pdg_iter;
1056 
1057  //-- get the 4-momentum of the i-th final state particle
1058  TLorentzVector * p4fin = fPhaseSpaceGenerator.GetDecay(i);
1059 
1060  new ( plist[offset+i] ) GHepParticle(
1061  pdgc, /* PDG Code */
1062  kIStStableFinalState, /* GHepStatus_t */
1063  -1, /* first mother particle */
1064  -1, /* second mother particle */
1065  -1, /* first daughter particle */
1066  -1, /* last daughter particle */
1067  p4fin->Px(), /* 4-momentum: px component */
1068  p4fin->Py(), /* 4-momentum: py component */
1069  p4fin->Pz(), /* 4-momentum: pz component */
1070  p4fin->Energy(), /* 4-momentum: E component */
1071  0, /* production vertex 4-vector: vx */
1072  0, /* production vertex 4-vector: vy */
1073  0, /* production vertex 4-vector: vz */
1074  0 /* production vertex 4-vector: time */
1075  );
1076  i++;
1077  }
1078 
1079  // Clean-up
1080  delete [] mass;
1081 
1082  return true;
1083 }
static const double m
Definition: Units.h:79
double ReWeightPt2(const PDGCodeList &pdgcv) const
#define pERROR
Definition: Messenger.h:60
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:34
static const unsigned int kMaxUnweightDecayIterations
Definition: Controls.h:62
intermediate_table::const_iterator const_iterator
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
size_t size
Definition: lodepng.cpp:55
bool fGenerateWeighted
generate weighted events?
#define pINFO
Definition: Messenger.h:63
double fWeight
weight for generated event
#define pWARN
Definition: Messenger.h:61
TGenPhaseSpace fPhaseSpaceGenerator
a phase space generator
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:54
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
#define pNOTICE
Definition: Messenger.h:62
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
void KNOHadronization::ProcessEventRecord ( GHepRecord event) const
virtual

Implements genie::EventRecordVisitorI.

Definition at line 97 of file KNOHadronization.cxx.

97  {
98 
99  Interaction * interaction = event->Summary();
100  TClonesArray * particle_list = this->Hadronize(interaction);
101 
102  if(! particle_list ) {
103  LOG("KNOHadronization", pWARN) << "Got an empty particle list. Hadronizer failed!";
104  LOG("KNOHadronization", pWARN) << "Quitting the current event generation thread";
105 
106  event->EventFlags()->SetBitNumber(kHadroSysGenErr, true);
107 
109  exception.SetReason("Could not simulate the hadronic system");
110  exception.SwitchOnFastForward();
111  throw exception;
112 
113  return;
114  }
115 
116  int mom = event->FinalStateHadronicSystemPosition();
117  assert(mom!=-1);
118 
119  // find the proper status for the particles we are going to put in event record
120  bool is_nucleus = interaction->InitState().Tgt().IsNucleus();
121  GHepStatus_t istfin = (is_nucleus) ?
123 
124  // retrieve the hadronic blob lorentz boost
125  // Because Hadronize() returned particles not in the LAB reference frame
126  const TLorentzVector * had_syst = event -> Particle(mom) -> P4() ;
127  TVector3 boost = had_syst -> BoostVector() ;
128 
129  GHepParticle * neutrino = event->Probe();
130  const TLorentzVector & vtx = *(neutrino->X4());
131 
132  GHepParticle * particle = 0;
133  TIter particle_iter(particle_list);
134  while ((particle = (GHepParticle *) particle_iter.Next())) {
135 
136  int pdgc = particle -> Pdg() ;
137 
138  // bring the particle in the LAB reference frame
139  particle -> P4() -> Boost( boost ) ;
140 
141  // set the proper status according to a number of things:
142  // interaction on a nucleaus or nucleon, particle type
143  GHepStatus_t ist = ( particle -> Status() ==1 ) ? istfin : kIStDISPreFragmHadronicState;
144 
145  // handle gammas, and leptons that might come from internal pythia decays
146  // mark them as final state particles
147  bool not_hadron = ( pdgc == kPdgGamma ||
148  pdg::IsNeutralLepton(pdgc) ||
149  pdg::IsChargedLepton(pdgc) ) ;
150 
151  if(not_hadron) { ist = kIStStableFinalState; }
152  particle -> SetStatus( ist ) ;
153 
154  int im = mom + 1 + particle -> FirstMother() ;
155  //int ifc = ( particle -> FirstDaughter() == -1) ? -1 : mom + 1 + particle -> FirstDaughter();
156  //int ilc = ( particle -> LastDaughter() == -1) ? -1 : mom + 1 + particle -> LastDaughter();
157 
158  particle -> SetFirstMother( im ) ;
159 
160  particle -> SetPosition( vtx ) ;
161 
162  event->AddParticle(*particle);
163  }
164 
165  delete particle_list ;
166 
167  // update the weight of the event
168  event -> SetWeight ( Weight() * event->Weight() );
169 
170 }
enum genie::EGHepStatus GHepStatus_t
bool IsNucleus(void) const
Definition: Target.cxx:289
TClonesArray * Hadronize(const Interaction *) const
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:99
virtual double Weight(void) const
Definition: GHepRecord.h:125
Summary information for an interaction.
Definition: Interaction.h:55
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const int kPdgGamma
Definition: PDGCodes.h:170
#define pWARN
Definition: Messenger.h:61
bool IsNeutralLepton(int pdgc)
Definition: PDGUtils.cxx:93
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
double Weight(void) const
const InitialState & InitState(void) const
Definition: Interaction.h:68
const Target & Tgt(void) const
Definition: InitialState.h:67
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
double KNOHadronization::ReWeightPt2 ( const PDGCodeList pdgcv) const
private

Definition at line 1085 of file KNOHadronization.cxx.

1086 {
1087 // Phase Space Decay re-weighting to reproduce exp(-pT2/<pT2>) pion pT2
1088 // distributions.
1089 // See: A.B.Clegg, A.Donnachie, A Description of Jet Structure by pT-limited
1090 // Phase Space.
1091 
1092  double w = 1;
1093 
1094  for(unsigned int i = 0; i < pdgcv.size(); i++) {
1095 
1096  //int pdgc = pdgcv[i];
1097  //if(pdgc!=kPdgPiP&&pdgc!=kPdgPiM) continue;
1098 
1099  TLorentzVector * p4 = fPhaseSpaceGenerator.GetDecay(i);
1100  double pt2 = TMath::Power(p4->Px(),2) + TMath::Power(p4->Py(),2);
1101  double wi = TMath::Exp(-fPhSpRwA*TMath::Sqrt(pt2));
1102  //double wi = (9.41 * TMath::Landau(pt2,0.24,0.12));
1103 
1104  w *= wi;
1105  }
1106  return w;
1107 }
size_t size
Definition: lodepng.cpp:55
TGenPhaseSpace fPhaseSpaceGenerator
a phase space generator
double fPhSpRwA
parameter for phase space decay reweighting
PDGCodeList * KNOHadronization::SelectParticles ( const Interaction interaction) const
private

Definition at line 244 of file KNOHadronization.cxx.

246 {
247  if(!this->AssertValidity(interaction)) {
248  LOG("KNOHad", pWARN) << "Returning a null particle list!";
249  return 0;
250  }
251 
252  unsigned int min_mult = 2;
253  unsigned int mult = 0;
254  PDGCodeList * pdgcv = 0;
255 
256  double W = utils::kinematics::W(interaction);
257 
258  //-- Get the charge that the hadron shower needs to have so as to
259  // conserve charge in the interaction
260  int maxQ = this->HadronShowerCharge(interaction);
261  LOG("KNOHad", pINFO) << "Hadron Shower Charge = " << maxQ;
262 
263  //-- Build the multiplicity probabilities for the input interaction
264  LOG("KNOHad", pDEBUG) << "Building Multiplicity Probability distribution";
265  LOG("KNOHad", pDEBUG) << *interaction;
266  Option_t * opt = "+LowMultSuppr+Renormalize";
267  TH1D * mprob = this->MultiplicityProb(interaction,opt);
268 
269  if(!mprob) {
270  LOG("KNOHad", pWARN) << "Null multiplicity probability distribution!";
271  return 0;
272  }
273  if(mprob->Integral("width")<=0) {
274  LOG("KNOHad", pWARN) << "Empty multiplicity probability distribution!";
275  delete mprob;
276  return 0;
277  }
278 
279  //----- FIND AN ALLOWED SOLUTION FOR THE HADRONIC FINAL STATE
280 
281  bool allowed_state=false;
282  unsigned int itry = 0;
283 
284  while(!allowed_state)
285  {
286  itry++;
287 
288  //-- Go in error if a solution has not been found after many attempts
289  if(itry>kMaxKNOHadSystIterations) {
290  LOG("KNOHad", pERROR)
291  << "Couldn't select hadronic shower particles after: "
292  << itry << " attempts!";
293  delete mprob;
294  return 0;
295  }
296 
297  //-- Generate a hadronic multiplicity
298  mult = TMath::Nint( mprob->GetRandom() );
299 
300  LOG("KNOHad", pINFO) << "Hadron multiplicity = " << mult;
301 
302  //-- Check that the generated multiplicity is consistent with the charge
303  // that the hadronic shower is required to have - else retry
304  if(mult < (unsigned int) TMath::Abs(maxQ)) {
305  LOG("KNOHad", pWARN)
306  << "Multiplicity not enough to generate hadronic charge! Retrying.";
307  allowed_state = false;
308  continue;
309  }
310 
311  //-- Force a min multiplicity
312  // This should never happen if the multiplicity probability distribution
313  // was properly built
314  if(mult < min_mult) {
315  if(fForceMinMult) {
316  LOG("KNOHad", pWARN)
317  << "Low generated multiplicity: " << mult
318  << ". Forcing to minimum accepted multiplicity: " << min_mult;
319  mult = min_mult;
320  } else {
321  LOG("KNOHad", pWARN)
322  << "Generated multiplicity: " << mult << " is too low! Quitting";
323  delete mprob;
324  return 0;
325  }
326  }
327 
328  //-- Determine what kind of particles we have in the final state
329  pdgcv = this->GenerateHadronCodes(mult, maxQ, W);
330 
331  LOG("KNOHad", pNOTICE)
332  << "Generated multiplicity (@ W = " << W << "): " << pdgcv->size();
333 
334  // muliplicity might have been forced to smaller value if the invariant
335  // mass of the hadronic system was not sufficient
336  mult = pdgcv->size(); // update for potential change
337 
338  // is it an allowed decay?
339  double msum=0;
341  for(pdg_iter = pdgcv->begin(); pdg_iter != pdgcv->end(); ++pdg_iter) {
342  int pdgc = *pdg_iter;
343  double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
344 
345  msum += m;
346  LOG("KNOHad", pDEBUG) << "- PDGC=" << pdgc << ", m=" << m << " GeV";
347  }
348  bool permitted = (W > msum);
349 
350  if(!permitted) {
351  LOG("KNOHad", pWARN) << "*** Decay forbidden by kinematics! ***";
352  LOG("KNOHad", pWARN) << "sum{mass} = " << msum << ", W = " << W;
353  LOG("KNOHad", pWARN) << "Discarding hadronic system & re-trying!";
354  delete pdgcv;
355  allowed_state = false;
356  continue;
357  }
358 
359  allowed_state = true;
360 
361  LOG("KNOHad", pNOTICE)
362  << "Found an allowed hadronic state @ W=" << W
363  << " multiplicity=" << mult;
364 
365  } // attempts
366 
367  delete mprob;
368 
369  return pdgcv;
370 }
static const double m
Definition: Units.h:79
#define pERROR
Definition: Messenger.h:60
int HadronShowerCharge(const Interaction *) const
opt
Definition: train.py:196
intermediate_table::const_iterator const_iterator
A list of PDG codes.
Definition: PDGCodeList.h:33
double W(const Interaction *const i)
Definition: KineUtils.cxx:1031
PDGCodeList * GenerateHadronCodes(int mult, int maxQ, double W) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
size_t size
Definition: lodepng.cpp:55
TH1D * MultiplicityProb(const Interaction *, Option_t *opt="") const
#define pINFO
Definition: Messenger.h:63
bool AssertValidity(const Interaction *i) const
#define pWARN
Definition: Messenger.h:61
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
bool fForceMinMult
force minimum multiplicity if (at low W) generated less?
static const unsigned int kMaxKNOHadSystIterations
Definition: Controls.h:58
#define pNOTICE
Definition: Messenger.h:62
#define pDEBUG
Definition: Messenger.h:64
double KNOHadronization::Weight ( void  ) const
private

Definition at line 483 of file KNOHadronization.cxx.

484 {
485  return fWeight;
486 }
double fWeight
weight for generated event
double KNOHadronization::Wmin ( void  ) const
private

Definition at line 1681 of file KNOHadronization.cxx.

1682 {
1683  return (kNucleonMass+kPionMass);
1684 }
static const double kNucleonMass
Definition: Constants.h:78
static const double kPionMass
Definition: Constants.h:74

Friends And Related Function Documentation

friend class KNOTunedQPMDISPXSec
friend

Definition at line 73 of file KNOHadronization.h.

Member Data Documentation

double genie::KNOHadronization::fAhyperon
private

parameter controlling strange baryon production probability via associated production (P=a+b*lnW^2)

Definition at line 134 of file KNOHadronization.h.

double genie::KNOHadronization::fAvbn
private

offset in average charged hadron multiplicity = f(W) relation for vbn

Definition at line 129 of file KNOHadronization.h.

double genie::KNOHadronization::fAvbp
private

offset in average charged hadron multiplicity = f(W) relation for vbp

Definition at line 128 of file KNOHadronization.h.

double genie::KNOHadronization::fAvn
private

offset in average charged hadron multiplicity = f(W) relation for vn

Definition at line 127 of file KNOHadronization.h.

double genie::KNOHadronization::fAvp
private

offset in average charged hadron multiplicity = f(W) relation for vp

Definition at line 126 of file KNOHadronization.h.

TF1* genie::KNOHadronization::fBaryonPT2pdf
private

baryon pT^2 PDF

Definition at line 141 of file KNOHadronization.h.

TF1* genie::KNOHadronization::fBaryonXFpdf
private

baryon xF PDF

Definition at line 140 of file KNOHadronization.h.

double genie::KNOHadronization::fBhyperon
private

see above

Definition at line 135 of file KNOHadronization.h.

double genie::KNOHadronization::fBvbn
private

slope in average charged hadron multiplicity = f(W) relation for vbn

Definition at line 133 of file KNOHadronization.h.

double genie::KNOHadronization::fBvbp
private

slope in average charged hadron multiplicity = f(W) relation for vbp

Definition at line 132 of file KNOHadronization.h.

double genie::KNOHadronization::fBvn
private

slope in average charged hadron multiplicity = f(W) relation for vn

Definition at line 131 of file KNOHadronization.h.

double genie::KNOHadronization::fBvp
private

slope in average charged hadron multiplicity = f(W) relation for vp

Definition at line 130 of file KNOHadronization.h.

double genie::KNOHadronization::fCvbn
private

Levy function parameter for vbn.

Definition at line 139 of file KNOHadronization.h.

double genie::KNOHadronization::fCvbp
private

Levy function parameter for vbp.

Definition at line 138 of file KNOHadronization.h.

double genie::KNOHadronization::fCvn
private

Levy function parameter for vn.

Definition at line 137 of file KNOHadronization.h.

double genie::KNOHadronization::fCvp
private

Levy function parameter for vp.

Definition at line 136 of file KNOHadronization.h.

bool genie::KNOHadronization::fForceDecays
private

force decays of unstable hadrons produced?

Definition at line 116 of file KNOHadronization.h.

bool genie::KNOHadronization::fForceMinMult
private

force minimum multiplicity if (at low W) generated less?

Definition at line 117 of file KNOHadronization.h.

bool genie::KNOHadronization::fForceNeuGenLimit
private

force upper hadronic multiplicity to NeuGEN limit

Definition at line 111 of file KNOHadronization.h.

bool genie::KNOHadronization::fGenerateWeighted
private

generate weighted events?

Definition at line 118 of file KNOHadronization.h.

double genie::KNOHadronization::fPeta
private

{eta eta} production probability

Definition at line 125 of file KNOHadronization.h.

TGenPhaseSpace genie::KNOHadronization::fPhaseSpaceGenerator
mutableprivate

a phase space generator

Definition at line 104 of file KNOHadronization.h.

double genie::KNOHadronization::fPhSpRwA
private

parameter for phase space decay reweighting

Definition at line 119 of file KNOHadronization.h.

double genie::KNOHadronization::fPK0
private

{K0 K0bar} production probability

Definition at line 123 of file KNOHadronization.h.

double genie::KNOHadronization::fPKc
private

{K+ K- } production probability

Definition at line 122 of file KNOHadronization.h.

double genie::KNOHadronization::fPpi0
private

{pi0 pi0 } production probability

Definition at line 120 of file KNOHadronization.h.

double genie::KNOHadronization::fPpi0eta
private

{Pi0 eta} production probability

Definition at line 124 of file KNOHadronization.h.

double genie::KNOHadronization::fPpic
private

{pi+ pi- } production probability

Definition at line 121 of file KNOHadronization.h.

bool genie::KNOHadronization::fReWeightDecays
private

Reweight phase space decays?

Definition at line 115 of file KNOHadronization.h.

double genie::KNOHadronization::fRvbnCCm2
private

Rijk: vbn, CC, multiplicity = 2.

Definition at line 157 of file KNOHadronization.h.

double genie::KNOHadronization::fRvbnCCm3
private

Rijk: vbn, CC, multiplicity = 3.

Definition at line 158 of file KNOHadronization.h.

double genie::KNOHadronization::fRvbnNCm2
private

Rijk: vbn, NC, multiplicity = 2.

Definition at line 159 of file KNOHadronization.h.

double genie::KNOHadronization::fRvbnNCm3
private

Rijk: vbn, NC, multiplicity = 3.

Definition at line 160 of file KNOHadronization.h.

double genie::KNOHadronization::fRvbpCCm2
private

Rijk: vbp, CC, multiplicity = 2.

Definition at line 153 of file KNOHadronization.h.

double genie::KNOHadronization::fRvbpCCm3
private

Rijk: vbp, CC, multiplicity = 3.

Definition at line 154 of file KNOHadronization.h.

double genie::KNOHadronization::fRvbpNCm2
private

Rijk: vbp, NC, multiplicity = 2.

Definition at line 155 of file KNOHadronization.h.

double genie::KNOHadronization::fRvbpNCm3
private

Rijk: vbp, NC, multiplicity = 3.

Definition at line 156 of file KNOHadronization.h.

double genie::KNOHadronization::fRvnCCm2
private

Rijk: vn, CC, multiplicity = 2.

Definition at line 149 of file KNOHadronization.h.

double genie::KNOHadronization::fRvnCCm3
private

Rijk: vn, CC, multiplicity = 3.

Definition at line 150 of file KNOHadronization.h.

double genie::KNOHadronization::fRvnNCm2
private

Rijk: vn, NC, multiplicity = 2.

Definition at line 151 of file KNOHadronization.h.

double genie::KNOHadronization::fRvnNCm3
private

Rijk: vn, NC, multiplicity = 3.

Definition at line 152 of file KNOHadronization.h.

double genie::KNOHadronization::fRvpCCm2
private

Rijk: vp, CC, multiplicity = 2.

Definition at line 145 of file KNOHadronization.h.

double genie::KNOHadronization::fRvpCCm3
private

Rijk: vp, CC, multiplicity = 3.

Definition at line 146 of file KNOHadronization.h.

double genie::KNOHadronization::fRvpNCm2
private

Rijk: vp, NC, multiplicity = 2.

Definition at line 147 of file KNOHadronization.h.

double genie::KNOHadronization::fRvpNCm3
private

Rijk: vp, NC, multiplicity = 3.

Definition at line 148 of file KNOHadronization.h.

bool genie::KNOHadronization::fUseBaryonXfPt2Param
private

Generate baryon xF,pT2 from experimental parameterization?

Definition at line 114 of file KNOHadronization.h.

bool genie::KNOHadronization::fUseIsotropic2BDecays
private

force isotropic, non-reweighted 2-body decays for consistency with neugen/daikon

Definition at line 113 of file KNOHadronization.h.

double genie::KNOHadronization::fWcut
private

Rijk applied for W<Wcut (see DIS/RES join scheme)

Definition at line 144 of file KNOHadronization.h.

double genie::KNOHadronization::fWeight
mutableprivate

weight for generated event

Definition at line 105 of file KNOHadronization.h.


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