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:
genie::EventRecordVisitorI genie::Algorithm

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

Private Member Functions

void LoadConfig (void)
 
void Initialize (void) const
 
TClonesArray * Hadronize (const Interaction *) const
 
double Weight (void) const
 
PDGCodeListSelectParticles (const Interaction *) const
 
TH1D * MultiplicityProb (const Interaction *, Option_t *opt="") const
 
bool AssertValidity (const Interaction *i) const
 
PDGCodeListGenerateHadronCodes (int mult, int maxQ, double W) const
 
int GenerateBaryonPdgCode (int mult, int maxQ, double W) const
 
int HadronShowerCharge (const Interaction *) const
 
double KNO (int nu, int nuc, double z) const
 
double AverageChMult (int nu, int nuc, double W) const
 
void HandleDecays (TClonesArray *particle_list) const
 
double ReWeightPt2 (const PDGCodeList &pdgcv) const
 
double MaxMult (const Interaction *i) const
 
TH1D * CreateMultProbHist (double maxmult) const
 
void ApplyRijk (const Interaction *i, bool norm, TH1D *mp) const
 
double Wmin (void) const
 
TClonesArray * DecayMethod1 (double W, const PDGCodeList &pdgv, bool reweight_decays) const
 
TClonesArray * DecayMethod2 (double W, const PDGCodeList &pdgv, bool reweight_decays) const
 
TClonesArray * DecayBackToBack (double W, const PDGCodeList &pdgv) const
 
bool PhaseSpaceDecay (TClonesArray &pl, TLorentzVector &pd, const PDGCodeList &pdgv, int offset=0, bool reweight=false) const
 

Private Attributes

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

Friends

class KNOTunedQPMDISPXSec
 

Additional Inherited Members

- Static Public Member Functions inherited from genie::Algorithm
static string BuildParamVectKey (const std::string &comm_name, unsigned int i)
 
static string BuildParamVectSizeKey (const std::string &comm_name)
 
- Protected Member Functions inherited from genie::EventRecordVisitorI
 EventRecordVisitorI ()
 
 EventRecordVisitorI (string name)
 
 EventRecordVisitorI (string name, string config)
 
- Protected Member Functions inherited from genie::Algorithm
 Algorithm ()
 
 Algorithm (string name)
 
 Algorithm (string name, string config)
 
void Initialize (void)
 
void DeleteConfig (void)
 
void DeleteSubstructure (void)
 
RegistryExtractLocalConfig (const Registry &in) const
 
RegistryExtractLowerConfig (const Registry &in, const string &alg_key) const
 Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key. More...
 
template<class T >
bool GetParam (const RgKey &name, T &p, bool is_top_call=true) const
 
template<class T >
bool GetParamDef (const RgKey &name, T &p, const T &def) const
 
template<class T >
int GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
 Handle to load vectors of parameters. More...
 
int GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
 
int AddTopRegistry (Registry *rp, bool owns=true)
 add registry with top priority, also update ownership More...
 
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership More...
 
int MergeTopRegistry (const Registry &r)
 
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships. More...
 
- Protected Attributes inherited from genie::Algorithm
bool fAllowReconfig
 
bool fOwnsSubstruc
 true if it owns its substructure (sub-algs,...) More...
 
AlgId fID
 algorithm name and configuration set More...
 
vector< Registry * > fConfVect
 
vector< boolfOwnerships
 ownership for every registry in fConfVect More...
 
AlgStatus_t fStatus
 algorithm execution status More...
 
AlgMapfOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool) More...
 

Detailed Description

A KNO-based hadronization model.

Is a concrete implementation of the EventRecordVisitorI interface.

Author
The main authors of this model are:
      - Costas Andreopoulos <constantinos.andreopoulos \at cern.ch>
        University of Liverpool & STFC Rutherford Appleton Laboratory

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

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

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

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

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

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

August 17, 2004

Copyright (c) 2003-2020, 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.

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

Definition at line 72 of file AGKYLowW2019.cxx.

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

Definition at line 80 of file AGKYLowW2019.cxx.

81 {
82  if (fBaryonXFpdf ) delete fBaryonXFpdf;
83  if (fBaryonPT2pdf) delete fBaryonPT2pdf;
84 //if (fKNO ) delete fKNO;
85 }
TF1 * fBaryonXFpdf
baryon xF PDF
Definition: AGKYLowW2019.h:141
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 1567 of file AGKYLowW2019.cxx.

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

Definition at line 1533 of file AGKYLowW2019.cxx.

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

Definition at line 663 of file AGKYLowW2019.cxx.

665 {
666 // computes the average charged multiplicity
667 //
668  bool is_p = pdg::IsProton (nuc_pdg);
669  bool is_n = pdg::IsNeutron (nuc_pdg);
670  bool is_nu = pdg::IsNeutrino (probe_pdg);
671  bool is_nubar = pdg::IsAntiNeutrino (probe_pdg);
672  bool is_l = pdg::IsNegChargedLepton (probe_pdg);
673  bool is_lbar = pdg::IsPosChargedLepton (probe_pdg);
674  // EDIT
675  bool is_dm = pdg::IsDarkMatter (probe_pdg);
676 
677  double a=0, b=0; // params controlling average multiplicity
678 
679  if ( is_p && (is_nu || is_l ) ) { a=fAvp; b=fBvp; }
680  else if ( is_n && (is_nu || is_l ) ) { a=fAvn; b=fBvn; }
681  else if ( is_p && (is_nubar || is_lbar) ) { a=fAvbp; b=fBvbp; }
682  else if ( is_n && (is_nubar || is_lbar) ) { a=fAvbn; b=fBvbn; }
683  // EDIT: assume it's neutrino-like for now...
684  else if ( is_p && is_dm ) { a=fAvp; b=fBvp; }
685  else if ( is_n && is_dm ) { a=fAvn; b=fBvn; }
686  else {
687  LOG("KNOHad", pERROR)
688  << "Invalid initial state (probe = " << probe_pdg << ", "
689  << "hit nucleon = " << nuc_pdg << ")";
690  return 0;
691  }
692 
693  double av_nch = a + b * 2*TMath::Log(W);
694  return av_nch;
695 }
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:107
#define pERROR
Definition: Messenger.h:59
bool IsDarkMatter(int pdgc)
Definition: PDGUtils.cxx:124
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:338
bool IsPosChargedLepton(int pdgc)
Definition: PDGUtils.cxx:145
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:333
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const double a
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:115
double fAvn
offset in average charged hadron multiplicity = f(W) relation for vn
Definition: AGKYLowW2019.h:128
double fAvbp
offset in average charged hadron multiplicity = f(W) relation for vbp
Definition: AGKYLowW2019.h:129
static bool * b
Definition: config.cpp:1043
double fBvp
slope in average charged hadron multiplicity = f(W) relation for vp
Definition: AGKYLowW2019.h:131
double fBvbp
slope in average charged hadron multiplicity = f(W) relation for vbp
Definition: AGKYLowW2019.h:133
double fBvbn
slope in average charged hadron multiplicity = f(W) relation for vbn
Definition: AGKYLowW2019.h:134
double fAvbn
offset in average charged hadron multiplicity = f(W) relation for vbn
Definition: AGKYLowW2019.h:130
double fBvn
slope in average charged hadron multiplicity = f(W) relation for vn
Definition: AGKYLowW2019.h:132
bool IsNegChargedLepton(int pdgc)
Definition: PDGUtils.cxx:136
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 487 of file AGKYLowW2019.cxx.

488 {
490  this->LoadConfig();
491 }
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 493 of file AGKYLowW2019.cxx.

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

Definition at line 1555 of file AGKYLowW2019.cxx.

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

Definition at line 868 of file AGKYLowW2019.cxx.

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

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

Definition at line 752 of file AGKYLowW2019.cxx.

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

Definition at line 1380 of file AGKYLowW2019.cxx.

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

Definition at line 1106 of file AGKYLowW2019.cxx.

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

Definition at line 169 of file AGKYLowW2019.cxx.

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

Definition at line 697 of file AGKYLowW2019.cxx.

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

Definition at line 1463 of file AGKYLowW2019.cxx.

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

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

Definition at line 501 of file AGKYLowW2019.cxx.

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

1548 {
1549  double W = interaction->Kine().W();
1550 
1551  double maxmult = TMath::Floor(1 + (W-kNeutronMass)/kPionMass);
1552  return maxmult;
1553 }
static const double kNeutronMass
Definition: Constants.h:76
static const double kPionMass
Definition: Constants.h:73
TH1D * AGKYLowW2019::MultiplicityProb ( const Interaction interaction,
Option_t *  opt = "" 
) const
private

Definition at line 369 of file AGKYLowW2019.cxx.

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

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

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 
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) ?
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 boost = had_syst -> BoostVector() ;
125 
126  GHepParticle * neutrino = event->Probe();
127  const TLorentzVector & vtx = *(neutrino->X4());
128 
129  GHepParticle * particle = 0;
130  TIter particle_iter(particle_list);
131  while ((particle = (GHepParticle *) particle_iter.Next())) {
132 
133  int pdgc = particle -> Pdg() ;
134 
135  // bring the particle in the LAB reference frame
136  particle -> P4() -> Boost( boost ) ;
137 
138  // set the proper status according to a number of things:
139  // interaction on a nucleaus or nucleon, particle type
140  GHepStatus_t ist = ( particle -> Status() ==1 ) ? istfin : kIStDISPreFragmHadronicState;
141 
142  // handle gammas, and leptons that might come from internal pythia decays
143  // mark them as final state particles
144  bool not_hadron = ( pdgc == kPdgGamma ||
145  pdg::IsNeutralLepton(pdgc) ||
146  pdg::IsChargedLepton(pdgc) ) ;
147 
148  if(not_hadron) { ist = kIStStableFinalState; }
149  particle -> SetStatus( ist ) ;
150 
151  int im = mom + 1 + particle -> FirstMother() ;
152  //int ifc = ( particle -> FirstDaughter() == -1) ? -1 : mom + 1 + particle -> FirstDaughter();
153  //int ilc = ( particle -> LastDaughter() == -1) ? -1 : mom + 1 + particle -> LastDaughter();
154 
155  particle -> SetFirstMother( im ) ;
156 
157  particle -> SetPosition( vtx ) ;
158 
159  event->AddParticle(*particle);
160  }
161 
162  delete particle_list ;
163 
164  // update the weight of the event
165  event -> SetWeight ( Weight() * event->Weight() );
166 
167 }
double Weight(void) const
bool IsNucleus(void) const
Definition: Target.cxx:272
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:98
virtual double Weight(void) const
Definition: GHepRecord.h:124
Summary information for an interaction.
Definition: Interaction.h:56
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const int kPdgGamma
Definition: PDGCodes.h:189
#define pWARN
Definition: Messenger.h:60
bool IsNeutralLepton(int pdgc)
Definition: PDGUtils.cxx:92
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
const InitialState & InitState(void) const
Definition: Interaction.h:69
const Target & Tgt(void) const
Definition: InitialState.h:66
TClonesArray * Hadronize(const Interaction *) const
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
enum genie::EGHepStatus GHepStatus_t
double AGKYLowW2019::ReWeightPt2 ( const PDGCodeList pdgcv) const
private

Definition at line 1082 of file AGKYLowW2019.cxx.

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

Definition at line 241 of file AGKYLowW2019.cxx.

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

Definition at line 480 of file AGKYLowW2019.cxx.

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

Definition at line 1678 of file AGKYLowW2019.cxx.

1679 {
1680  return (kNucleonMass+kPionMass);
1681 }
static const double kNucleonMass
Definition: Constants.h:77
static const double kPionMass
Definition: Constants.h:73

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.

double genie::AGKYLowW2019::fAvbn
private

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

Definition at line 130 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fAvbp
private

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

Definition at line 129 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fAvn
private

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

Definition at line 128 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fAvp
private

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

Definition at line 127 of file AGKYLowW2019.h.

TF1* genie::AGKYLowW2019::fBaryonPT2pdf
private

baryon pT^2 PDF

Definition at line 142 of file AGKYLowW2019.h.

TF1* genie::AGKYLowW2019::fBaryonXFpdf
private

baryon xF PDF

Definition at line 141 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fBhyperon
private

see above

Definition at line 136 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fBvbn
private

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

Definition at line 134 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fBvbp
private

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

Definition at line 133 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fBvn
private

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

Definition at line 132 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fBvp
private

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

Definition at line 131 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fCvbn
private

Levy function parameter for vbn.

Definition at line 140 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fCvbp
private

Levy function parameter for vbp.

Definition at line 139 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fCvn
private

Levy function parameter for vn.

Definition at line 138 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fCvp
private

Levy function parameter for vp.

Definition at line 137 of file AGKYLowW2019.h.

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.

bool genie::AGKYLowW2019::fForceNeuGenLimit
private

force upper hadronic multiplicity to NeuGEN limit

Definition at line 112 of file AGKYLowW2019.h.

bool genie::AGKYLowW2019::fGenerateWeighted
private

generate weighted events?

Definition at line 119 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fPeta
private

{eta eta} production probability

Definition at line 126 of file AGKYLowW2019.h.

TGenPhaseSpace genie::AGKYLowW2019::fPhaseSpaceGenerator
mutableprivate

a phase space generator

Definition at line 105 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fPhSpRwA
private

parameter for phase space decay reweighting

Definition at line 120 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fPK0
private

{K0 K0bar} production probability

Definition at line 124 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fPKc
private

{K+ K- } production probability

Definition at line 123 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fPpi0
private

{pi0 pi0 } production probability

Definition at line 121 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fPpi0eta
private

{Pi0 eta} production probability

Definition at line 125 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fPpic
private

{pi+ pi- } production probability

Definition at line 122 of file AGKYLowW2019.h.

bool genie::AGKYLowW2019::fReWeightDecays
private

Reweight phase space decays?

Definition at line 116 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fRvbnCCm2
private

Rijk: vbn, CC, multiplicity = 2.

Definition at line 158 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fRvbnCCm3
private

Rijk: vbn, CC, multiplicity = 3.

Definition at line 159 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fRvbnNCm2
private

Rijk: vbn, NC, multiplicity = 2.

Definition at line 160 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fRvbnNCm3
private

Rijk: vbn, NC, multiplicity = 3.

Definition at line 161 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fRvbpCCm2
private

Rijk: vbp, CC, multiplicity = 2.

Definition at line 154 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fRvbpCCm3
private

Rijk: vbp, CC, multiplicity = 3.

Definition at line 155 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fRvbpNCm2
private

Rijk: vbp, NC, multiplicity = 2.

Definition at line 156 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fRvbpNCm3
private

Rijk: vbp, NC, multiplicity = 3.

Definition at line 157 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fRvnCCm2
private

Rijk: vn, CC, multiplicity = 2.

Definition at line 150 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fRvnCCm3
private

Rijk: vn, CC, multiplicity = 3.

Definition at line 151 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fRvnNCm2
private

Rijk: vn, NC, multiplicity = 2.

Definition at line 152 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fRvnNCm3
private

Rijk: vn, NC, multiplicity = 3.

Definition at line 153 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fRvpCCm2
private

Rijk: vp, CC, multiplicity = 2.

Definition at line 146 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fRvpCCm3
private

Rijk: vp, CC, multiplicity = 3.

Definition at line 147 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fRvpNCm2
private

Rijk: vp, NC, multiplicity = 2.

Definition at line 148 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fRvpNCm3
private

Rijk: vp, NC, multiplicity = 3.

Definition at line 149 of file AGKYLowW2019.h.

bool genie::AGKYLowW2019::fUseBaryonXfPt2Param
private

Generate baryon xF,pT2 from experimental parameterization?

Definition at line 115 of file AGKYLowW2019.h.

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.

double genie::AGKYLowW2019::fWcut
private

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

Definition at line 145 of file AGKYLowW2019.h.

double genie::AGKYLowW2019::fWeight
mutableprivate

weight for generated event

Definition at line 106 of file AGKYLowW2019.h.


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