GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Private Member Functions | Private Attributes | Friends | List of all members
genie::AGKYLowW2019 Class Reference

A KNO-based hadronization model. More...

#include <AGKYLowW2019.h>

Inheritance diagram for genie::AGKYLowW2019:
Inheritance graph
[legend]
Collaboration diagram for genie::AGKYLowW2019:
Collaboration graph
[legend]

Public Member Functions

 AGKYLowW2019 ()
 
 AGKYLowW2019 (string config)
 
virtual ~AGKYLowW2019 ()
 
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 Registry & GetConfig (void) const
 
Registry * GetOwnedConfig (void)
 
virtual const AlgId & Id (void) const
 Get algorithm ID. More...
 
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status. More...
 
virtual bool AllowReconfig (void) const
 
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm. More...
 
virtual void SetId (const AlgId &id)
 Set algorithm ID. More...
 
virtual void SetId (string name, string config)
 
const Algorithm * SubAlg (const RgKey &registry_key) const
 
void AdoptConfig (void)
 
void AdoptSubstructure (void)
 
virtual void Print (ostream &stream) const
 Print algorithm info. More...
 

Private Member Functions

void LoadConfig (void)
 
void Initialize (void) const
 
TClonesArray * Hadronize (const Interaction *) const
 
double Weight (void) const
 
PDGCodeList * SelectParticles (const Interaction *) const
 
TH1D * MultiplicityProb (const Interaction *, Option_t *opt="") const
 
bool AssertValidity (const Interaction *i) const
 
PDGCodeList * GenerateHadronCodes (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 fRvpEMm2
 Rijk: vp, EM, multiplicity = 2. More...
 
double fRvpEMm3
 Rijk: vp, EM, 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 fRvnEMm2
 Rijk: vn, EM, multiplicity = 2. More...
 
double fRvnEMm3
 Rijk: vn, EM, 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 fRvbpEMm2
 Rijk: vbp, EM, multiplicity = 2. More...
 
double fRvbpEMm3
 Rijk: vbp, EM, 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...
 
double fRvbnEMm2
 Rijk: vbn, EM, multiplicity = 2. More...
 
double fRvbnEMm3
 Rijk: vbn, EM, 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)
 
static string BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)
 
static string BuildParamMatRowSizeKey (const std::string &comm_name)
 
static string BuildParamMatColSizeKey (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)
 
Registry * ExtractLocalConfig (const Registry &in) const
 
Registry * ExtractLowerConfig (const Registry &in, const string &alg_key) const
 Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key. More...
 
template<class T >
bool GetParam (const RgKey &name, T &p, bool is_top_call=true) const
 
template<class T >
bool GetParamDef (const RgKey &name, T &p, const T &def) const
 
template<class T >
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
 
template<class T >
int GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
 Handle to load matrix of parameters. More...
 
template<class T >
int GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const
 
int GetParamMatKeys (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< bool > fOwnerships
 ownership for every registry in fConfVect More...
 
AlgStatus_t fStatus
 algorithm execution status More...
 
AlgMap * fOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool) More...
 

Detailed Description

A KNO-based hadronization model. 

      Is a concrete implementation of the EventRecordVisitorI interface.
Author
The main authors of this model are: 
      - Costas Andreopoulos <c.andreopoulos \at cern.ch>
        University of Liverpool

      - 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)
Created:
August 17, 2004
License:
Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 58 of file AGKYLowW2019.h.

Constructor & Destructor Documentation

AGKYLowW2019::AGKYLowW2019 ( )

Definition at line 64 of file AGKYLowW2019.cxx.

References fBaryonPT2pdf, and fBaryonXFpdf.

64  :
65 EventRecordVisitorI("genie::AGKYLowW2019")
66 {
67  fBaryonXFpdf = 0;
68  fBaryonPT2pdf = 0;
69 //fKNO = 0;
70 }
genie::AGKYLowW2019::fBaryonXFpdf
TF1 * fBaryonXFpdf
baryon xF PDF
Definition: AGKYLowW2019.h:141
genie::AGKYLowW2019::fBaryonPT2pdf
TF1 * fBaryonPT2pdf
baryon pT^2 PDF
Definition: AGKYLowW2019.h:142
AGKYLowW2019::AGKYLowW2019 ( string  config)

Definition at line 72 of file AGKYLowW2019.cxx.

References fBaryonPT2pdf, and fBaryonXFpdf.

72  :
73 EventRecordVisitorI("genie::AGKYLowW2019", config)
74 {
75  fBaryonXFpdf = 0;
76  fBaryonPT2pdf = 0;
77 //fKNO = 0;
78 }
genie::AGKYLowW2019::fBaryonXFpdf
TF1 * fBaryonXFpdf
baryon xF PDF
Definition: AGKYLowW2019.h:141
genie::AGKYLowW2019::fBaryonPT2pdf
TF1 * fBaryonPT2pdf
baryon pT^2 PDF
Definition: AGKYLowW2019.h:142
AGKYLowW2019::~AGKYLowW2019 ( )
virtual

Definition at line 80 of file AGKYLowW2019.cxx.

References fBaryonPT2pdf, and fBaryonXFpdf.

81 {
82  if (fBaryonXFpdf ) delete fBaryonXFpdf;
83  if (fBaryonPT2pdf) delete fBaryonPT2pdf;
84 //if (fKNO ) delete fKNO;
85 }
genie::AGKYLowW2019::fBaryonXFpdf
TF1 * fBaryonXFpdf
baryon xF PDF
Definition: AGKYLowW2019.h:141
genie::AGKYLowW2019::fBaryonPT2pdf
TF1 * fBaryonPT2pdf
baryon pT^2 PDF
Definition: AGKYLowW2019.h:142

Member Function Documentation

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

Definition at line 1578 of file AGKYLowW2019.cxx.

References fRvbnCCm2, fRvbnCCm3, fRvbnEMm2, fRvbnEMm3, fRvbnNCm2, fRvbnNCm3, fRvbpCCm2, fRvbpCCm3, fRvbpEMm2, fRvbpEMm3, fRvbpNCm2, fRvbpNCm3, fRvnCCm2, fRvnCCm3, fRvnEMm2, fRvnEMm3, fRvnNCm2, fRvnNCm3, fRvpCCm2, fRvpCCm3, fRvpEMm2, fRvpEMm3, fRvpNCm2, fRvpNCm3, genie::Target::HitNucPdg(), genie::Interaction::InitState(), genie::pdg::IsAntiNeutrino(), genie::pdg::IsDarkMatter(), genie::ProcessInfo::IsDarkMatter(), genie::ProcessInfo::IsEM(), genie::pdg::IsNegChargedLepton(), genie::pdg::IsNeutrino(), genie::pdg::IsNeutron(), genie::pdg::IsPosChargedLepton(), genie::pdg::IsProton(), genie::ProcessInfo::IsWeakCC(), genie::ProcessInfo::IsWeakNC(), LOG, pDEBUG, pERROR, genie::InitialState::ProbePdg(), genie::Interaction::ProcInfo(), and genie::InitialState::Tgt().

Referenced by MultiplicityProb().

1580 {
1581  // Apply the NEUGEN multiplicity probability scaling factors
1582  //
1583  if(!mp) return;
1584 
1585  const InitialState & init_state = interaction->InitState();
1586  int probe_pdg = init_state.ProbePdg();
1587  int nuc_pdg = init_state.Tgt().HitNucPdg();
1588 
1589  const ProcessInfo & proc_info = interaction->ProcInfo();
1590  bool is_CC = proc_info.IsWeakCC();
1591  bool is_NC = proc_info.IsWeakNC();
1592  bool is_EM = proc_info.IsEM();
1593  // EDIT
1594  bool is_dm = proc_info.IsDarkMatter();
1595 
1596  //
1597  // get the R2, R3 factors
1598  //
1599 
1600  double R2=1., R3=1.;
1601 
1602  // weak CC or NC case
1603  // EDIT
1604  if(is_CC || is_NC || is_dm) {
1605  bool is_nu = pdg::IsNeutrino (probe_pdg);
1606  bool is_nubar = pdg::IsAntiNeutrino (probe_pdg);
1607  bool is_p = pdg::IsProton (nuc_pdg);
1608  bool is_n = pdg::IsNeutron (nuc_pdg);
1609  bool is_dmi = pdg::IsDarkMatter (probe_pdg); // EDIT
1610  if((is_nu && is_p) || (is_dmi && is_p)) {
1611  R2 = (is_CC) ? fRvpCCm2 : fRvpNCm2;
1612  R3 = (is_CC) ? fRvpCCm3 : fRvpNCm3;
1613  } else
1614  if((is_nu && is_n) || (is_dmi && is_n)) {
1615  R2 = (is_CC) ? fRvnCCm2 : fRvnNCm2;
1616  R3 = (is_CC) ? fRvnCCm3 : fRvnNCm3;
1617  } else
1618  if(is_nubar && is_p) {
1619  R2 = (is_CC) ? fRvbpCCm2 : fRvbpNCm2;
1620  R3 = (is_CC) ? fRvbpCCm3 : fRvbpNCm3;
1621  } else
1622  if(is_nubar && is_n) {
1623  R2 = (is_CC) ? fRvbnCCm2 : fRvbnNCm2;
1624  R3 = (is_CC) ? fRvbnCCm3 : fRvbnNCm3;
1625  } else {
1626  LOG("Hadronization", pERROR)
1627  << "Invalid initial state: " << init_state;
1628  }
1629  }//cc||nc?
1630 
1631  // EM case (apply the NC tuning factors)
1632 
1633  if(is_EM) {
1634  bool is_l = pdg::IsNegChargedLepton (probe_pdg);
1635  bool is_lbar = pdg::IsPosChargedLepton (probe_pdg);
1636  bool is_p = pdg::IsProton (nuc_pdg);
1637  bool is_n = pdg::IsNeutron (nuc_pdg);
1638  if(is_l && is_p) {
1639  R2 = fRvpEMm2;
1640  R3 = fRvpEMm3;
1641  } else
1642  if(is_l && is_n) {
1643  R2 = fRvnEMm2;
1644  R3 = fRvnEMm3;
1645  } else
1646  if(is_lbar && is_p) {
1647  R2 = fRvbpEMm2;
1648  R3 = fRvbpEMm3;
1649  } else
1650  if(is_lbar && is_n) {
1651  R2 = fRvbnEMm2;
1652  R3 = fRvbnEMm3;
1653  } else {
1654  LOG("Hadronization", pERROR)
1655  << "Invalid initial state: " << init_state;
1656  }
1657  }//em?
1658 
1659  //
1660  // Apply to the multiplicity probability distribution
1661  //
1662 
1663  int nbins = mp->GetNbinsX();
1664  for(int i = 1; i <= nbins; i++) {
1665  int n = TMath::Nint( mp->GetBinCenter(i) );
1666 
1667  double R=1;
1668  if (n==2) R=R2;
1669  else if (n==3) R=R3;
1670 
1671  if(n==2 || n==3) {
1672  double P = mp->GetBinContent(i);
1673  double Psc = R*P;
1674  LOG("Hadronization", pDEBUG)
1675  << "n=" << n << "/ Scaling factor R = "
1676  << R << "/ P " << P << " --> " << Psc;
1677  mp->SetBinContent(i, Psc);
1678  }
1679  if(n>3) break;
1680  }
1681 
1682  // renormalize the histogram?
1683  if(norm) {
1684  double histo_norm = mp->Integral("width");
1685  if(histo_norm>0) mp->Scale(1.0/histo_norm);
1686  }
1687 }
genie::AGKYLowW2019::fRvnEMm3
double fRvnEMm3
Rijk: vn, EM, multiplicity = 3.
Definition: AGKYLowW2019.h:157
genie::ProcessInfo::IsWeakCC
bool IsWeakCC(void) const
Definition: ProcessInfo.cxx:203
genie::pdg::IsNeutrino
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
pERROR
#define pERROR
Definition: Messenger.h:59
genie::Target::HitNucPdg
int HitNucPdg(void) const
Definition: Target.cxx:304
genie::AGKYLowW2019::fRvbpCCm3
double fRvbpCCm3
Rijk: vbp, CC, multiplicity = 3.
Definition: AGKYLowW2019.h:159
genie::AGKYLowW2019::fRvbnEMm2
double fRvbnEMm2
Rijk: vbn, EM, multiplicity = 2.
Definition: AGKYLowW2019.h:168
genie::AGKYLowW2019::fRvnCCm2
double fRvnCCm2
Rijk: vn, CC, multiplicity = 2.
Definition: AGKYLowW2019.h:152
genie::AGKYLowW2019::fRvpCCm2
double fRvpCCm2
Rijk: vp, CC, multiplicity = 2.
Definition: AGKYLowW2019.h:146
genie::AGKYLowW2019::fRvbpNCm2
double fRvbpNCm2
Rijk: vbp, NC, multiplicity = 2.
Definition: AGKYLowW2019.h:160
genie::pdg::IsDarkMatter
bool IsDarkMatter(int pdgc)
Definition: PDGUtils.cxx:127
genie::AGKYLowW2019::fRvpNCm2
double fRvpNCm2
Rijk: vp, NC, multiplicity = 2.
Definition: AGKYLowW2019.h:148
genie::AGKYLowW2019::fRvbnCCm3
double fRvbnCCm3
Rijk: vbn, CC, multiplicity = 3.
Definition: AGKYLowW2019.h:165
genie::AGKYLowW2019::fRvnNCm3
double fRvnNCm3
Rijk: vn, NC, multiplicity = 3.
Definition: AGKYLowW2019.h:155
genie::AGKYLowW2019::fRvbpNCm3
double fRvbpNCm3
Rijk: vbp, NC, multiplicity = 3.
Definition: AGKYLowW2019.h:161
genie::pdg::IsNeutron
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
genie::pdg::IsPosChargedLepton
bool IsPosChargedLepton(int pdgc)
Definition: PDGUtils.cxx:148
genie::AGKYLowW2019::fRvbpEMm2
double fRvbpEMm2
Rijk: vbp, EM, multiplicity = 2.
Definition: AGKYLowW2019.h:162
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
genie::ProcessInfo::IsWeakNC
bool IsWeakNC(void) const
Definition: ProcessInfo.cxx:208
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::AGKYLowW2019::fRvbpCCm2
double fRvbpCCm2
Rijk: vbp, CC, multiplicity = 2.
Definition: AGKYLowW2019.h:158
genie::AGKYLowW2019::fRvbpEMm3
double fRvbpEMm3
Rijk: vbp, EM, multiplicity = 3.
Definition: AGKYLowW2019.h:163
genie::ProcessInfo
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
genie::AGKYLowW2019::fRvpCCm3
double fRvpCCm3
Rijk: vp, CC, multiplicity = 3.
Definition: AGKYLowW2019.h:147
genie::pdg::IsAntiNeutrino
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
genie::AGKYLowW2019::fRvbnCCm2
double fRvbnCCm2
Rijk: vbn, CC, multiplicity = 2.
Definition: AGKYLowW2019.h:164
genie::InitialState::ProbePdg
int ProbePdg(void) const
Definition: InitialState.h:64
genie::AGKYLowW2019::fRvpNCm3
double fRvpNCm3
Rijk: vp, NC, multiplicity = 3.
Definition: AGKYLowW2019.h:149
genie::AGKYLowW2019::fRvnNCm2
double fRvnNCm2
Rijk: vn, NC, multiplicity = 2.
Definition: AGKYLowW2019.h:154
genie::AGKYLowW2019::fRvpEMm3
double fRvpEMm3
Rijk: vp, EM, multiplicity = 3.
Definition: AGKYLowW2019.h:151
genie::AGKYLowW2019::fRvbnEMm3
double fRvbnEMm3
Rijk: vbn, EM, multiplicity = 3.
Definition: AGKYLowW2019.h:169
genie::AGKYLowW2019::fRvbnNCm2
double fRvbnNCm2
Rijk: vbn, NC, multiplicity = 2.
Definition: AGKYLowW2019.h:166
genie::AGKYLowW2019::fRvnEMm2
double fRvnEMm2
Rijk: vn, EM, multiplicity = 2.
Definition: AGKYLowW2019.h:156
genie::ProcessInfo::IsEM
bool IsEM(void) const
Definition: ProcessInfo.cxx:193
genie::AGKYLowW2019::fRvpEMm2
double fRvpEMm2
Rijk: vp, EM, multiplicity = 2.
Definition: AGKYLowW2019.h:150
genie::AGKYLowW2019::fRvbnNCm3
double fRvbnNCm3
Rijk: vbn, NC, multiplicity = 3.
Definition: AGKYLowW2019.h:167
genie::ProcessInfo::IsDarkMatter
bool IsDarkMatter(void) const
Definition: ProcessInfo.cxx:213
genie::InitialState::Tgt
const Target & Tgt(void) const
Definition: InitialState.h:66
genie::AGKYLowW2019::fRvnCCm3
double fRvnCCm3
Rijk: vn, CC, multiplicity = 3.
Definition: AGKYLowW2019.h:153
genie::pdg::IsNegChargedLepton
bool IsNegChargedLepton(int pdgc)
Definition: PDGUtils.cxx:139
genie::InitialState
Initial State information.
Definition: InitialState.h:48
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
bool AGKYLowW2019::AssertValidity ( const Interaction *  i) const
private

Definition at line 1544 of file AGKYLowW2019.cxx.

References genie::Interaction::ExclTag(), genie::XclsTag::IsCharmEvent(), LOG, pWARN, genie::utils::kinematics::W(), and Wmin().

Referenced by Hadronize(), MultiplicityProb(), and SelectParticles().

1545 {
1546  if(interaction->ExclTag().IsCharmEvent()) {
1547  LOG("KNOHad", pWARN) << "Can't hadronize charm events";
1548  return false;
1549  }
1550  double W = utils::kinematics::W(interaction);
1551  if(W < this->Wmin()) {
1552  LOG("KNOHad", pWARN) << "Low invariant mass, W = " << W << " GeV!!";
1553  return false;
1554  }
1555  return true;
1556 }
genie::AGKYLowW2019::Wmin
double Wmin(void) const
Definition: AGKYLowW2019.cxx:1689
genie::utils::kinematics::W
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
pWARN
#define pWARN
Definition: Messenger.h:60
double AGKYLowW2019::AverageChMult ( int  nu,
int  nuc,
double  W 
) const
private

Definition at line 674 of file AGKYLowW2019.cxx.

References a, genie::units::b, fAvbn, fAvbp, fAvn, fAvp, fBvbn, fBvbp, fBvn, fBvp, genie::pdg::IsAntiNeutrino(), genie::pdg::IsDarkMatter(), genie::pdg::IsNegChargedLepton(), genie::pdg::IsNeutrino(), genie::pdg::IsNeutron(), genie::pdg::IsPosChargedLepton(), genie::pdg::IsProton(), LOG, and pERROR.

Referenced by MultiplicityProb().

676 {
677 // computes the average charged multiplicity
678 //
679  bool is_p = pdg::IsProton (nuc_pdg);
680  bool is_n = pdg::IsNeutron (nuc_pdg);
681  bool is_nu = pdg::IsNeutrino (probe_pdg);
682  bool is_nubar = pdg::IsAntiNeutrino (probe_pdg);
683  bool is_l = pdg::IsNegChargedLepton (probe_pdg);
684  bool is_lbar = pdg::IsPosChargedLepton (probe_pdg);
685  // EDIT
686  bool is_dm = pdg::IsDarkMatter (probe_pdg);
687 
688  double a=0, b=0; // params controlling average multiplicity
689 
690  if ( is_p && (is_nu || is_l ) ) { a=fAvp; b=fBvp; }
691  else if ( is_n && (is_nu || is_l ) ) { a=fAvn; b=fBvn; }
692  else if ( is_p && (is_nubar || is_lbar) ) { a=fAvbp; b=fBvbp; }
693  else if ( is_n && (is_nubar || is_lbar) ) { a=fAvbn; b=fBvbn; }
694  // EDIT: assume it's neutrino-like for now...
695  else if ( is_p && is_dm ) { a=fAvp; b=fBvp; }
696  else if ( is_n && is_dm ) { a=fAvn; b=fBvn; }
697  else {
698  LOG("KNOHad", pERROR)
699  << "Invalid initial state (probe = " << probe_pdg << ", "
700  << "hit nucleon = " << nuc_pdg << ")";
701  return 0;
702  }
703 
704  double av_nch = a + b * 2*TMath::Log(W);
705  return av_nch;
706 }
genie::pdg::IsNeutrino
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
pERROR
#define pERROR
Definition: Messenger.h:59
genie::pdg::IsDarkMatter
bool IsDarkMatter(int pdgc)
Definition: PDGUtils.cxx:127
genie::units::b
static constexpr double b
Definition: Units.h:78
genie::utils::kinematics::W
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
genie::pdg::IsNeutron
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
genie::pdg::IsPosChargedLepton
bool IsPosChargedLepton(int pdgc)
Definition: PDGUtils.cxx:148
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
a
const double a
Definition: gUpMuFluxGen.cxx:164
genie::pdg::IsAntiNeutrino
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
genie::AGKYLowW2019::fAvn
double fAvn
offset in average charged hadron multiplicity = f(W) relation for vn
Definition: AGKYLowW2019.h:128
genie::AGKYLowW2019::fAvbp
double fAvbp
offset in average charged hadron multiplicity = f(W) relation for vbp
Definition: AGKYLowW2019.h:129
genie::AGKYLowW2019::fBvp
double fBvp
slope in average charged hadron multiplicity = f(W) relation for vp
Definition: AGKYLowW2019.h:131
genie::AGKYLowW2019::fBvbp
double fBvbp
slope in average charged hadron multiplicity = f(W) relation for vbp
Definition: AGKYLowW2019.h:133
genie::AGKYLowW2019::fBvbn
double fBvbn
slope in average charged hadron multiplicity = f(W) relation for vbn
Definition: AGKYLowW2019.h:134
genie::AGKYLowW2019::fAvbn
double fAvbn
offset in average charged hadron multiplicity = f(W) relation for vbn
Definition: AGKYLowW2019.h:130
genie::AGKYLowW2019::fBvn
double fBvn
slope in average charged hadron multiplicity = f(W) relation for vn
Definition: AGKYLowW2019.h:132
genie::pdg::IsNegChargedLepton
bool IsNegChargedLepton(int pdgc)
Definition: PDGUtils.cxx:139
genie::AGKYLowW2019::fAvp
double fAvp
offset in average charged hadron multiplicity = f(W) relation for vp
Definition: AGKYLowW2019.h:127
void AGKYLowW2019::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 488 of file AGKYLowW2019.cxx.

References genie::Algorithm::Configure(), and LoadConfig().

489 {
490  EventRecordVisitorI::Configure(config);
491  this->LoadConfig();
492 }
genie::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void AGKYLowW2019::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 494 of file AGKYLowW2019.cxx.

References genie::Algorithm::Configure(), and LoadConfig().

495 {
496  EventRecordVisitorI::Configure(config);
497  this->LoadConfig();
498 }
genie::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
TH1D * AGKYLowW2019::CreateMultProbHist ( double  maxmult) const
private

Definition at line 1566 of file AGKYLowW2019.cxx.

Referenced by MultiplicityProb().

1567 {
1568  double minmult = 2;
1569  int nbins = TMath::Nint(maxmult-minmult+1);
1570 
1571  TH1D * mult_prob = new TH1D("mult_prob",
1572  "hadronic multiplicity distribution", nbins, minmult-0.5, maxmult+0.5);
1573  mult_prob->SetDirectory(0);
1574 
1575  return mult_prob;
1576 }
TClonesArray * AGKYLowW2019::DecayBackToBack ( double  W,
const PDGCodeList &  pdgv 
) const
private

Definition at line 879 of file AGKYLowW2019.cxx.

References fBaryonPT2pdf, fBaryonXFpdf, genie::RandomGen::Instance(), genie::pdg::IsNeutronOrProton(), LOG, genie::GHepParticle::Pdg(), pERROR, PhaseSpaceDecay(), pINFO, genie::GHepParticle::Px(), genie::GHepParticle::Py(), and genie::GHepParticle::Pz().

Referenced by DecayMethod2().

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

Definition at line 740 of file AGKYLowW2019.cxx.

References LOG, PhaseSpaceDecay(), and pINFO.

Referenced by Hadronize().

742 {
743 // Simple phase space decay including all generated particles.
744 // The old NeuGEN decay strategy.
745 
746  LOG("KNOHad", pINFO) << "** Using Hadronic System Decay method 1";
747 
748  TLorentzVector p4had(0,0,0,W);
749  TClonesArray * plist = new TClonesArray("genie::GHepParticle", pdgv.size());
750 
751  // do the decay
752  bool ok = this->PhaseSpaceDecay(*plist, p4had, pdgv, 0, reweight_decays);
753 
754  // clean-up and return NULL
755  if(!ok) {
756  plist->Delete();
757  delete plist;
758  return 0;
759  }
760  return plist;
761 }
genie::utils::kinematics::W
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
pINFO
#define pINFO
Definition: Messenger.h:62
genie::AGKYLowW2019::PhaseSpaceDecay
bool PhaseSpaceDecay(TClonesArray &pl, TLorentzVector &pd, const PDGCodeList &pdgv, int offset=0, bool reweight=false) const
Definition: AGKYLowW2019.cxx:945
TClonesArray * AGKYLowW2019::DecayMethod2 ( double  W,
const PDGCodeList &  pdgv,
bool  reweight_decays 
) const
private

Definition at line 763 of file AGKYLowW2019.cxx.

References DecayBackToBack(), fBaryonPT2pdf, fBaryonXFpdf, genie::PDGLibrary::Find(), genie::RandomGen::Instance(), genie::PDGLibrary::Instance(), genie::kIStStableFinalState, genie::constants::kPi, LOG, genie::utils::print::P4AsString(), pDEBUG, pERROR, PhaseSpaceDecay(), and pINFO.

Referenced by Hadronize().

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

Definition at line 1391 of file AGKYLowW2019.cxx.

References fAhyperon, fBhyperon, genie::RandomGen::Instance(), genie::kPdgLambda, genie::kPdgNeutron, genie::kPdgProton, genie::kPdgSigmaM, genie::kPdgSigmaP, LOG, pDEBUG, and genie::RandomGen::RndHadro().

Referenced by GenerateHadronCodes().

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

Definition at line 1117 of file AGKYLowW2019.cxx.

References genie::PDGLibrary::Find(), fPeta, fPK0, fPKc, fPpi0, fPpi0eta, fPpic, GenerateBaryonPdgCode(), genie::RandomGen::Instance(), genie::PDGLibrary::Instance(), genie::kPdgAntiK0, genie::kPdgEta, genie::kPdgK0, genie::kPdgKM, genie::kPdgKP, genie::kPdgLambda, genie::kPdgPi0, genie::kPdgPiM, genie::kPdgPiP, genie::kPdgProton, genie::kPdgSigmaM, genie::kPdgSigmaP, LOG, genie::utils::res::Mass(), pDEBUG, pERROR, genie::PDGCodeList::push_back(), and genie::RandomGen::RndHadro().

Referenced by SelectParticles().

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

Definition at line 170 of file AGKYLowW2019.cxx.

References AssertValidity(), DecayMethod1(), DecayMethod2(), fReWeightDecays, fUseBaryonXfPt2Param, fUseIsotropic2BDecays, fWeight, HandleDecays(), LOG, pINFO, pNOTICE, pWARN, SelectParticles(), and genie::utils::kinematics::W().

Referenced by ProcessEventRecord().

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

Definition at line 708 of file AGKYLowW2019.cxx.

References genie::PDGLibrary::Find(), genie::Interaction::FSPrimLepton(), genie::Target::HitNucPdg(), genie::Interaction::InitState(), genie::PDGLibrary::Instance(), genie::pdg::IsNeutron(), genie::pdg::IsProton(), genie::InitialState::Probe(), and genie::InitialState::Tgt().

Referenced by SelectParticles().

709 {
710 // Returns the hadron shower charge in units of +e
711 // HadronShowerCharge = Q{initial} - Q{final state primary lepton}
712 // eg in v p -> l- X the hadron shower charge is +2
713 // in v n -> l- X the hadron shower charge is +1
714 // in v n -> v X the hadron shower charge is 0
715 //
716  int hadronShowerCharge = 0;
717 
718  // find out the charge of the final state lepton
719  double ql = interaction->FSPrimLepton()->Charge() / 3.;
720 
721  // find out the charge of the probe
722  double qp = interaction->InitState().Probe()->Charge() / 3.;
723 
724  // get the initial state, ask for the hit-nucleon and get
725  // its charge ( = initial state charge for vN interactions)
726  const InitialState & init_state = interaction->InitState();
727  int hit_nucleon = init_state.Tgt().HitNucPdg();
728 
729  assert( pdg::IsProton(hit_nucleon) || pdg::IsNeutron(hit_nucleon) );
730 
731  // Ask PDGLibrary for the nucleon charge
732  double qnuc = PDGLibrary::Instance()->Find(hit_nucleon)->Charge() / 3.;
733 
734  // calculate the hadron shower charge
735  hadronShowerCharge = (int) ( qp + qnuc - ql );
736 
737  return hadronShowerCharge;
738 }
genie::Target::HitNucPdg
int HitNucPdg(void) const
Definition: Target.cxx:304
genie::InitialState::Probe
TParticlePDG * Probe(void) const
Definition: InitialState.cxx:147
genie::pdg::IsNeutron
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
genie::Interaction::FSPrimLepton
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
Definition: Interaction.cxx:126
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
genie::Interaction::InitState
const InitialState & InitState(void) const
Definition: Interaction.h:69
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
genie::InitialState::Tgt
const Target & Tgt(void) const
Definition: InitialState.h:66
genie::InitialState
Initial State information.
Definition: InitialState.h:48
void AGKYLowW2019::HandleDecays ( TClonesArray *  particle_list) const
private

Definition at line 1474 of file AGKYLowW2019.cxx.

Referenced by Hadronize().

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

Definition at line 89 of file AGKYLowW2019.cxx.

90 {
91 
92 }
double AGKYLowW2019::KNO ( int  nu,
int  nuc,
double  z 
) const
private

Definition at line 639 of file AGKYLowW2019.cxx.

References fCvbn, fCvbp, fCvn, fCvp, genie::pdg::IsAntiNeutrino(), genie::pdg::IsDarkMatter(), genie::pdg::IsNegChargedLepton(), genie::pdg::IsNeutrino(), genie::pdg::IsNeutron(), genie::pdg::IsPosChargedLepton(), genie::pdg::IsProton(), LOG, and pERROR.

Referenced by MultiplicityProb().

640 {
641 // Computes <n>P(n) for the input reduced multiplicity z=n/<n>
642 
643  bool is_p = pdg::IsProton (nuc_pdg);
644  bool is_n = pdg::IsNeutron (nuc_pdg);
645  bool is_nu = pdg::IsNeutrino (probe_pdg);
646  bool is_nubar = pdg::IsAntiNeutrino (probe_pdg);
647  bool is_l = pdg::IsNegChargedLepton (probe_pdg);
648  bool is_lbar = pdg::IsPosChargedLepton (probe_pdg);
649  // EDIT
650  bool is_dm = pdg::IsDarkMatter (probe_pdg);
651 
652  double c=0; // Levy function parameter
653 
654  if ( is_p && (is_nu || is_l ) ) c=fCvp;
655  else if ( is_n && (is_nu || is_l ) ) c=fCvn;
656  else if ( is_p && (is_nubar || is_lbar) ) c=fCvbp;
657  else if ( is_n && (is_nubar || is_lbar) ) c=fCvbn;
658  // EDIT: assume it's neutrino-like for now...
659  else if ( is_p && is_dm ) c=fCvp;
660  else if ( is_n && is_dm ) c=fCvn;
661  else {
662  LOG("KNOHad", pERROR)
663  << "Invalid initial state (probe = " << probe_pdg << ", "
664  << "hit nucleon = " << nuc_pdg << ")";
665  return 0;
666  }
667 
668  double x = c*z+1;
669  double kno = 2*TMath::Exp(-c)*TMath::Power(c,x)/TMath::Gamma(x);
670 
671  return kno;
672 }
genie::pdg::IsNeutrino
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
pERROR
#define pERROR
Definition: Messenger.h:59
genie::pdg::IsDarkMatter
bool IsDarkMatter(int pdgc)
Definition: PDGUtils.cxx:127
genie::AGKYLowW2019::fCvbp
double fCvbp
Levy function parameter for vbp.
Definition: AGKYLowW2019.h:139
genie::pdg::IsNeutron
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
genie::pdg::IsPosChargedLepton
bool IsPosChargedLepton(int pdgc)
Definition: PDGUtils.cxx:148
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::AGKYLowW2019::fCvbn
double fCvbn
Levy function parameter for vbn.
Definition: AGKYLowW2019.h:140
genie::pdg::IsAntiNeutrino
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
genie::AGKYLowW2019::fCvp
double fCvp
Levy function parameter for vp.
Definition: AGKYLowW2019.h:137
genie::pdg::IsNegChargedLepton
bool IsNegChargedLepton(int pdgc)
Definition: PDGUtils.cxx:139
genie::AGKYLowW2019::fCvn
double fCvn
Levy function parameter for vn.
Definition: AGKYLowW2019.h:138
void AGKYLowW2019::LoadConfig ( void  )
private

Definition at line 502 of file AGKYLowW2019.cxx.

References fAhyperon, fAvbn, fAvbp, fAvn, fAvp, fBaryonPT2pdf, fBaryonXFpdf, fBhyperon, fBvbn, fBvbp, fBvn, fBvp, fCvbn, fCvbp, fCvn, fCvp, fForceMinMult, fGenerateWeighted, fPeta, fPhSpRwA, fPK0, fPKc, fPpi0, fPpi0eta, fPpic, fReWeightDecays, fRvbnCCm2, fRvbnCCm3, fRvbnEMm2, fRvbnEMm3, fRvbnNCm2, fRvbnNCm3, fRvbpCCm2, fRvbpCCm3, fRvbpEMm2, fRvbpEMm3, fRvbpNCm2, fRvbpNCm3, fRvnCCm2, fRvnCCm3, fRvnEMm2, fRvnEMm3, fRvnNCm2, fRvnNCm3, fRvpCCm2, fRvpCCm3, fRvpEMm2, fRvpEMm3, fRvpNCm2, fRvpNCm3, fUseBaryonXfPt2Param, fUseIsotropic2BDecays, fWcut, genie::Algorithm::GetParam(), genie::Algorithm::GetParamDef(), LOG, and pWARN.

Referenced by Configure().

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

Definition at line 1558 of file AGKYLowW2019.cxx.

References genie::Interaction::Kine(), genie::constants::kNeutronMass, genie::constants::kPionMass, genie::Kinematics::W(), and genie::utils::kinematics::W().

Referenced by MultiplicityProb().

1559 {
1560  double W = interaction->Kine().W();
1561 
1562  double maxmult = TMath::Floor(1 + (W-kNeutronMass)/kPionMass);
1563  return maxmult;
1564 }
genie::utils::kinematics::W
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
genie::constants::kNeutronMass
static const double kNeutronMass
Definition: Framework/Conventions/Constants.h:76
TH1D * AGKYLowW2019::MultiplicityProb ( const Interaction *  interaction,
Option_t *  opt = "" 
) const
private

Definition at line 370 of file AGKYLowW2019.cxx.

References ApplyRijk(), AssertValidity(), AverageChMult(), CreateMultProbHist(), fForceNeuGenLimit, fWcut, genie::Target::HitNucPdg(), genie::Interaction::InitState(), KNO(), LOG, MaxMult(), pDEBUG, pINFO, genie::InitialState::ProbePdg(), pWARN, SLOG, genie::InitialState::Tgt(), and genie::utils::kinematics::W().

Referenced by genie::KNOTunedQPMDISPXSec::DISRESJoinSuppressionFactor(), and SelectParticles().

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

Definition at line 945 of file AGKYLowW2019.cxx.

References fGenerateWeighted, genie::PDGLibrary::Find(), fPhaseSpaceGenerator, fWeight, genie::RandomGen::Instance(), genie::PDGLibrary::Instance(), genie::kIStStableFinalState, genie::controls::kMaxUnweightDecayIterations, LOG, genie::units::m, genie::utils::print::P4AsString(), pERROR, pINFO, pNOTICE, pWARN, ReWeightPt2(), and genie::RandomGen::RndHadro().

Referenced by DecayBackToBack(), DecayMethod1(), and DecayMethod2().

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

Implements genie::EventRecordVisitorI.

Definition at line 94 of file AGKYLowW2019.cxx.

References Hadronize(), genie::Interaction::InitState(), genie::pdg::IsChargedLepton(), genie::pdg::IsNeutralLepton(), genie::Target::IsNucleus(), genie::kHadroSysGenErr, genie::kIStDISPreFragmHadronicState, genie::kIStHadronInTheNucleus, genie::kIStStableFinalState, genie::kPdgGamma, LOG, pWARN, genie::exceptions::EVGThreadException::SetReason(), genie::exceptions::EVGThreadException::SwitchOnFastForward(), genie::InitialState::Tgt(), Weight(), genie::GHepRecord::Weight(), and genie::GHepParticle::X4().

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

Definition at line 1093 of file AGKYLowW2019.cxx.

References fPhaseSpaceGenerator, and fPhSpRwA.

Referenced by PhaseSpaceDecay().

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

Definition at line 242 of file AGKYLowW2019.cxx.

References AssertValidity(), fForceMinMult, genie::PDGLibrary::Find(), GenerateHadronCodes(), HadronShowerCharge(), genie::PDGLibrary::Instance(), genie::controls::kMaxKNOHadSystIterations, LOG, genie::units::m, MultiplicityProb(), pDEBUG, pERROR, pINFO, pNOTICE, pWARN, and genie::utils::kinematics::W().

Referenced by Hadronize().

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

Definition at line 481 of file AGKYLowW2019.cxx.

References fWeight.

Referenced by ProcessEventRecord().

482 {
483  return fWeight;
484 }
genie::AGKYLowW2019::fWeight
double fWeight
weight for generated event
Definition: AGKYLowW2019.h:106
double AGKYLowW2019::Wmin ( void  ) const
private

Definition at line 1689 of file AGKYLowW2019.cxx.

References genie::constants::kNucleonMass, and genie::constants::kPionMass.

Referenced by AssertValidity().

1690 {
1691  return (kNucleonMass+kPionMass);
1692 }
genie::constants::kNucleonMass
static const double kNucleonMass
Definition: Framework/Conventions/Constants.h:77

Friends And Related Function Documentation

friend class KNOTunedQPMDISPXSec
friend

Definition at line 74 of file AGKYLowW2019.h.

Member Data Documentation

double genie::AGKYLowW2019::fAhyperon
private

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

Definition at line 135 of file AGKYLowW2019.h.

Referenced by GenerateBaryonPdgCode(), and LoadConfig().

double genie::AGKYLowW2019::fAvbn
private

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

Definition at line 130 of file AGKYLowW2019.h.

Referenced by AverageChMult(), and LoadConfig().

double genie::AGKYLowW2019::fAvbp
private

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

Definition at line 129 of file AGKYLowW2019.h.

Referenced by AverageChMult(), and LoadConfig().

double genie::AGKYLowW2019::fAvn
private

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

Definition at line 128 of file AGKYLowW2019.h.

Referenced by AverageChMult(), and LoadConfig().

double genie::AGKYLowW2019::fAvp
private

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

Definition at line 127 of file AGKYLowW2019.h.

Referenced by AverageChMult(), and LoadConfig().

TF1* genie::AGKYLowW2019::fBaryonPT2pdf
private

baryon pT^2 PDF

Definition at line 142 of file AGKYLowW2019.h.

Referenced by AGKYLowW2019(), DecayBackToBack(), DecayMethod2(), LoadConfig(), and ~AGKYLowW2019().

TF1* genie::AGKYLowW2019::fBaryonXFpdf
private

baryon xF PDF

Definition at line 141 of file AGKYLowW2019.h.

Referenced by AGKYLowW2019(), DecayBackToBack(), DecayMethod2(), LoadConfig(), and ~AGKYLowW2019().

double genie::AGKYLowW2019::fBhyperon
private

see above

Definition at line 136 of file AGKYLowW2019.h.

Referenced by GenerateBaryonPdgCode(), and LoadConfig().

double genie::AGKYLowW2019::fBvbn
private

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

Definition at line 134 of file AGKYLowW2019.h.

Referenced by AverageChMult(), and LoadConfig().

double genie::AGKYLowW2019::fBvbp
private

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

Definition at line 133 of file AGKYLowW2019.h.

Referenced by AverageChMult(), and LoadConfig().

double genie::AGKYLowW2019::fBvn
private

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

Definition at line 132 of file AGKYLowW2019.h.

Referenced by AverageChMult(), and LoadConfig().

double genie::AGKYLowW2019::fBvp
private

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

Definition at line 131 of file AGKYLowW2019.h.

Referenced by AverageChMult(), and LoadConfig().

double genie::AGKYLowW2019::fCvbn
private

Levy function parameter for vbn.

Definition at line 140 of file AGKYLowW2019.h.

Referenced by KNO(), and LoadConfig().

double genie::AGKYLowW2019::fCvbp
private

Levy function parameter for vbp.

Definition at line 139 of file AGKYLowW2019.h.

Referenced by KNO(), and LoadConfig().

double genie::AGKYLowW2019::fCvn
private

Levy function parameter for vn.

Definition at line 138 of file AGKYLowW2019.h.

Referenced by KNO(), and LoadConfig().

double genie::AGKYLowW2019::fCvp
private

Levy function parameter for vp.

Definition at line 137 of file AGKYLowW2019.h.

Referenced by KNO(), and LoadConfig().

bool genie::AGKYLowW2019::fForceDecays
private

force decays of unstable hadrons produced?

Definition at line 117 of file AGKYLowW2019.h.

bool genie::AGKYLowW2019::fForceMinMult
private

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

Definition at line 118 of file AGKYLowW2019.h.

Referenced by LoadConfig(), and SelectParticles().

bool genie::AGKYLowW2019::fForceNeuGenLimit
private

force upper hadronic multiplicity to NeuGEN limit

Definition at line 112 of file AGKYLowW2019.h.

Referenced by MultiplicityProb().

bool genie::AGKYLowW2019::fGenerateWeighted
private

generate weighted events?

Definition at line 119 of file AGKYLowW2019.h.

Referenced by LoadConfig(), and PhaseSpaceDecay().

double genie::AGKYLowW2019::fPeta
private

{eta eta} production probability

Definition at line 126 of file AGKYLowW2019.h.

Referenced by GenerateHadronCodes(), and LoadConfig().

TGenPhaseSpace genie::AGKYLowW2019::fPhaseSpaceGenerator
mutableprivate

a phase space generator

Definition at line 105 of file AGKYLowW2019.h.

Referenced by PhaseSpaceDecay(), and ReWeightPt2().

double genie::AGKYLowW2019::fPhSpRwA
private

parameter for phase space decay reweighting

Definition at line 120 of file AGKYLowW2019.h.

Referenced by LoadConfig(), and ReWeightPt2().

double genie::AGKYLowW2019::fPK0
private

{K0 K0bar} production probability

Definition at line 124 of file AGKYLowW2019.h.

Referenced by GenerateHadronCodes(), and LoadConfig().

double genie::AGKYLowW2019::fPKc
private

{K+ K- } production probability

Definition at line 123 of file AGKYLowW2019.h.

Referenced by GenerateHadronCodes(), and LoadConfig().

double genie::AGKYLowW2019::fPpi0
private

{pi0 pi0 } production probability

Definition at line 121 of file AGKYLowW2019.h.

Referenced by GenerateHadronCodes(), and LoadConfig().

double genie::AGKYLowW2019::fPpi0eta
private

{Pi0 eta} production probability

Definition at line 125 of file AGKYLowW2019.h.

Referenced by GenerateHadronCodes(), and LoadConfig().

double genie::AGKYLowW2019::fPpic
private

{pi+ pi- } production probability

Definition at line 122 of file AGKYLowW2019.h.

Referenced by GenerateHadronCodes(), and LoadConfig().

bool genie::AGKYLowW2019::fReWeightDecays
private

Reweight phase space decays?

Definition at line 116 of file AGKYLowW2019.h.

Referenced by Hadronize(), and LoadConfig().

double genie::AGKYLowW2019::fRvbnCCm2
private

Rijk: vbn, CC, multiplicity = 2.

Definition at line 164 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvbnCCm3
private

Rijk: vbn, CC, multiplicity = 3.

Definition at line 165 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvbnEMm2
private

Rijk: vbn, EM, multiplicity = 2.

Definition at line 168 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvbnEMm3
private

Rijk: vbn, EM, multiplicity = 3.

Definition at line 169 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvbnNCm2
private

Rijk: vbn, NC, multiplicity = 2.

Definition at line 166 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvbnNCm3
private

Rijk: vbn, NC, multiplicity = 3.

Definition at line 167 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvbpCCm2
private

Rijk: vbp, CC, multiplicity = 2.

Definition at line 158 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvbpCCm3
private

Rijk: vbp, CC, multiplicity = 3.

Definition at line 159 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvbpEMm2
private

Rijk: vbp, EM, multiplicity = 2.

Definition at line 162 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvbpEMm3
private

Rijk: vbp, EM, multiplicity = 3.

Definition at line 163 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvbpNCm2
private

Rijk: vbp, NC, multiplicity = 2.

Definition at line 160 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvbpNCm3
private

Rijk: vbp, NC, multiplicity = 3.

Definition at line 161 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvnCCm2
private

Rijk: vn, CC, multiplicity = 2.

Definition at line 152 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvnCCm3
private

Rijk: vn, CC, multiplicity = 3.

Definition at line 153 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvnEMm2
private

Rijk: vn, EM, multiplicity = 2.

Definition at line 156 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvnEMm3
private

Rijk: vn, EM, multiplicity = 3.

Definition at line 157 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvnNCm2
private

Rijk: vn, NC, multiplicity = 2.

Definition at line 154 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvnNCm3
private

Rijk: vn, NC, multiplicity = 3.

Definition at line 155 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvpCCm2
private

Rijk: vp, CC, multiplicity = 2.

Definition at line 146 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvpCCm3
private

Rijk: vp, CC, multiplicity = 3.

Definition at line 147 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvpEMm2
private

Rijk: vp, EM, multiplicity = 2.

Definition at line 150 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvpEMm3
private

Rijk: vp, EM, multiplicity = 3.

Definition at line 151 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvpNCm2
private

Rijk: vp, NC, multiplicity = 2.

Definition at line 148 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

double genie::AGKYLowW2019::fRvpNCm3
private

Rijk: vp, NC, multiplicity = 3.

Definition at line 149 of file AGKYLowW2019.h.

Referenced by ApplyRijk(), and LoadConfig().

bool genie::AGKYLowW2019::fUseBaryonXfPt2Param
private

Generate baryon xF,pT2 from experimental parameterization?

Definition at line 115 of file AGKYLowW2019.h.

Referenced by Hadronize(), and LoadConfig().

bool genie::AGKYLowW2019::fUseIsotropic2BDecays
private

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

Definition at line 114 of file AGKYLowW2019.h.

Referenced by Hadronize(), and LoadConfig().

double genie::AGKYLowW2019::fWcut
private

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

Definition at line 145 of file AGKYLowW2019.h.

Referenced by LoadConfig(), and MultiplicityProb().

double genie::AGKYLowW2019::fWeight
mutableprivate

weight for generated event

Definition at line 106 of file AGKYLowW2019.h.

Referenced by Hadronize(), PhaseSpaceDecay(), and Weight().

The documentation for this class was generated from the following files:
Generated on Wed Jul 17 2024 14:46:57 for GENIEGenerator by   doxygen 1.8.5