Andreopoulos - Gallagher (AG) GENIE Charm Hadronization model. More...

#include <AGCharm2019.h>

Inheritance diagram for genie::AGCharm2019:

Collaboration diagram for genie::AGCharm2019:

Public Member Functions
	AGCharm2019 ()

	AGCharm2019 (string config)

virtual	~AGCharm2019 ()

void	ProcessEventRecord (GHepRecord *event) const

void	Configure (const Registry &config)

void	Configure (string config)

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

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

virtual void	FindConfig (void)

virtual const Registry &	GetConfig (void) const

Registry *	GetOwnedConfig (void)

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

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

virtual bool	AllowReconfig (void) const

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

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

virtual void	SetId (string name, string config)

const Algorithm *	SubAlg (const RgKey &registry_key) const

void	AdoptConfig (void)

void	AdoptSubstructure (void)

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

Private Member Functions
void	LoadConfig (void)

void	Initialize (void) const

TClonesArray *	Hadronize (const Interaction *) const

int	GenerateCharmHadron (int nupdg, double EvLab) const

double	Weight (void) const

Private Attributes
TGenPhaseSpace	fPhaseSpaceGenerator
	a phase space generator More...

bool	fCharmOnly
	don't hadronize non-charm blob More...

TF1 *	fCharmPT2pdf
	charm hadron pT^2 pdf More...

const FragmentationFunctionI *	fFragmFunc
	charm hadron fragmentation func More...

double	fFracMaxEnergy
	Maximum energy available for the Meson fractions. More...

Spline *	fD0FracSpl
	nu charm fraction vs Ev: D0 More...

Spline *	fDpFracSpl
	nu charm fraction vs Ev: D+ More...

Spline *	fDsFracSpl
	nu charm fraction vs Ev: Ds+ More...

double	fD0BarFrac
	nubar {D0} charm fraction More...

double	fDmFrac
	nubar D- charm fraction More...

TPythia6 *	fPythia
	remnant (non-charm) hadronizer More...

Additional Inherited Members
Static Public Member Functions inherited from genie::Algorithm
static string	BuildParamVectKey (const std::string &comm_name, unsigned int i)

static string	BuildParamVectSizeKey (const std::string &comm_name)

static string	BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)

static string	BuildParamMatRowSizeKey (const std::string &comm_name)

static string	BuildParamMatColSizeKey (const std::string &comm_name)

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

	EventRecordVisitorI (string name)

	EventRecordVisitorI (string name, string config)

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

	Algorithm (string name)

	Algorithm (string name, string config)

void	Initialize (void)

void	DeleteConfig (void)

void	DeleteSubstructure (void)

Registry *	ExtractLocalConfig (const Registry &in) const

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

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

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

template<class T >
int	GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
	Handle to load vectors of parameters. More...

int	GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const

template<class T >
int	GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
	Handle to load matrix of parameters. More...

template<class T >
int	GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const

int	GetParamMatKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const

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

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

int	MergeTopRegistry (const Registry &r)

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

Protected Attributes inherited from genie::Algorithm
bool	fAllowReconfig

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

AlgId	fID
	algorithm name and configuration set More...

vector< Registry * >	fConfVect

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

AlgStatus_t	fStatus
	algorithm execution status More...

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

Detailed Description

Andreopoulos - Gallagher (AG) GENIE Charm Hadronization model.

      The model relies on empirical charm fragmentation and pT functions,
      as well as on experimentally-determined charm fractions, to produce
      the ID and 4-momentum of charmed hadron in charm production events.

      The remnant (non-charm) system is hadronised by a call to PYTHIA.

      Is a concrete implementation of the EventRecordVisitorI interface.

Author: Costas Andreopoulos <c.andreopoulos cern.ch> University of Liverpool

Hugh Gallagher galla.nosp@m.g@mi.nosp@m.nos.p.nosp@m.hy.t.nosp@m.ufts..nosp@m.edu Tufts University

Created:: August 17, 2004

License:: Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 45 of file AGCharm2019.h.

Constructor & Destructor Documentation

AGCharm2019::AGCharm2019 ( )

Definition at line 55 of file AGCharm2019.cxx.

References Initialize().

                          :
 EventRecordVisitorI("genie::AGCharm2019")
 {
   this->Initialize();
 }

AGCharm2019::AGCharm2019 ( string config )

Definition at line 61 of file AGCharm2019.cxx.

References Initialize().

                                       :
 EventRecordVisitorI("genie::AGCharm2019", config)
 {
   this->Initialize();
 }

AGCharm2019::~AGCharm2019 ( )

virtual

Definition at line 67 of file AGCharm2019.cxx.

References fCharmPT2pdf, fD0FracSpl, fDpFracSpl, and fDsFracSpl.

 {
   delete fCharmPT2pdf;
   fCharmPT2pdf = 0;
 
   delete fD0FracSpl;
   fD0FracSpl = 0;
 
   delete fDpFracSpl;
   fDpFracSpl = 0;
 
   delete fDsFracSpl;
   fDsFracSpl = 0;
 }

Member Function Documentation

void AGCharm2019::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 785 of file AGCharm2019.cxx.

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

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

void AGCharm2019::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 791 of file AGCharm2019.cxx.

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

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

int AGCharm2019::GenerateCharmHadron	(	int	nupdg,
		double	EvLab
	)		const

private

Definition at line 744 of file AGCharm2019.cxx.

References genie::Spline::Evaluate(), fD0BarFrac, fD0FracSpl, fDmFrac, fDpFracSpl, fDsFracSpl, fFracMaxEnergy, genie::RandomGen::Instance(), genie::pdg::IsAntiNeutrino(), genie::pdg::IsNeutrino(), genie::kPdgAntiD0, genie::kPdgD0, genie::kPdgDM, genie::kPdgDMs, genie::kPdgDP, genie::kPdgDPs, genie::kPdgLambdaPc, LOG, pERROR, and genie::RandomGen::RndHadro().

Referenced by Hadronize().

 {
   // generate a charmed hadron pdg code using a charm fraction table
 
   RandomGen * rnd = RandomGen::Instance();
   double r  = rnd->RndHadro().Rndm();
 
   // the ratios are giving up to a certain value.
   // The rations saturates at high energies, so the values used above that enery
   // are the evaluated at
 
   if(pdg::IsNeutrino(nu_pdg)) {
 
     // the ratios are giving up to a certain value.
     // The rations saturates at high energies, so the values used above that enery
     // are the evaluated at the maximum energies avaiable for the ratios
 
     EvLab = TMath::Min( EvLab, fFracMaxEnergy ) ;
 
     double tf = 0;
     if      (r < (tf+=fD0FracSpl->Evaluate(EvLab)))  return kPdgD0;       // D^0
     else if (r < (tf+=fDpFracSpl->Evaluate(EvLab)))  return kPdgDP;       // D^+
     else if (r < (tf+=fDsFracSpl->Evaluate(EvLab)))  return kPdgDPs;      // Ds^+
     else                                             return kPdgLambdaPc; // Lamda_c^+
 
   } else if(pdg::IsAntiNeutrino(nu_pdg)) {
     if      (r < fD0BarFrac)          return kPdgAntiD0;
     else if (r < fD0BarFrac+fDmFrac)  return kPdgDM;
     else                              return kPdgDMs;
   }
 
   LOG("CharmHad", pERROR) << "Could not generate a charm hadron!";
   return 0;
 }

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

private

Definition at line 172 of file AGCharm2019.cxx.

References fCharmOnly, fCharmPT2pdf, fFragmFunc, genie::PDGLibrary::Find(), genie::GHepParticle::FirstDaughter(), genie::GHepParticle::FirstMother(), fPhaseSpaceGenerator, fPythia, genie::Interaction::FSPrimLepton(), GenerateCharmHadron(), genie::FragmentationFunctionI::GenerateZ(), genie::Kinematics::HadSystP4(), genie::Target::HitNucPdg(), genie::RandomGen::Instance(), genie::PDGLibrary::Instance(), genie::pdg::IsAntiNeutrino(), genie::pdg::IsDarkMatter(), genie::pdg::IsNeutrino(), genie::pdg::IsNeutron(), genie::pdg::IsProton(), genie::controls::kASmallNum, genie::kIStNucleonTarget, genie::kIStStableFinalState, genie::controls::kMaxUnweightDecayIterations, genie::constants::kNucleonMass, genie::kPdgAntiD0, genie::kPdgAntiDQuark, genie::kPdgD0, genie::kPdgDDDiquarkS1, genie::kPdgDM, genie::kPdgDMs, genie::kPdgDP, genie::kPdgDPs, genie::kPdgDQuark, genie::kPdgHadronicBlob, genie::kPdgLambdaPc, genie::kPdgNeutron, genie::kPdgP33m1232_Delta0, genie::kPdgP33m1232_DeltaM, genie::kPdgP33m1232_DeltaP, genie::kPdgP33m1232_DeltaPP, genie::kPdgPi0, genie::kPdgPiM, genie::kPdgPiP, genie::kPdgProton, genie::kPdgSQuark, genie::kPdgUDDiquarkS0, genie::kPdgUDDiquarkS1, genie::kPdgUQuark, genie::kPdgUUDiquarkS1, genie::constants::kPi, genie::constants::kPionMass, genie::kRfLab, genie::controls::kRjMaxIterations, genie::GHepParticle::LastDaughter(), LOG, genie::utils::print::P4AsString(), pDEBUG, pERROR, pINFO, pNOTICE, genie::InitialState::ProbeE(), genie::InitialState::ProbePdg(), genie::PDGCodeList::push_back(), pWARN, py2ent_(), genie::RandomGen::RndHadro(), genie::InitialState::Tgt(), genie::Kinematics::W(), and genie::utils::kinematics::W().

Referenced by ProcessEventRecord().

 {
   LOG("CharmHad", pNOTICE) << "** Running CHARM hadronizer";
 
   PDGLibrary * pdglib = PDGLibrary::Instance();
   RandomGen *  rnd    = RandomGen::Instance();
 
   // ....................................................................
   // Get information on the input event
   //
   const InitialState & init_state = interaction -> InitState();
   const Kinematics &   kinematics = interaction -> Kine();
   const Target &       target     = init_state.Tgt();
 
   const TLorentzVector & p4Had = kinematics.HadSystP4();
 
   double Ev = init_state.ProbeE(kRfLab);
   double W  = kinematics.W(true);
 
   TVector3 beta = -1 * p4Had.BoostVector(); // boost vector for LAB' -> HCM'
   TLorentzVector p4H(0,0,0,W);              // hadronic system 4p @ HCM'
 
   double Eh = p4Had.Energy();
 
   LOG("CharmHad", pNOTICE) << "Ehad (LAB) = " << Eh << ", W = " << W;
 
   int  nu_pdg  = init_state.ProbePdg();
   int  nuc_pdg = target.HitNucPdg();
 //int  qpdg    = target.HitQrkPdg();
 //bool sea     = target.HitSeaQrk();
   bool isp     = pdg::IsProton (nuc_pdg);
   bool isn     = pdg::IsNeutron(nuc_pdg);
   bool isnu    = pdg::IsNeutrino(nu_pdg);
   bool isnub   = pdg::IsAntiNeutrino(nu_pdg);
   bool isdm    = pdg::IsDarkMatter(nu_pdg);
 
   // ....................................................................
   // Attempt to generate a charmed hadron & its 4-momentum
   //
 
   TLorentzVector p4C(0,0,0,0);
   int ch_pdg = -1;
 
   bool got_charmed_hadron = false;
   unsigned int itry=0;
 
   while(itry++ < kRjMaxIterations && !got_charmed_hadron) {
 
      // Generate a charmed hadron PDG code
      int    pdg = this->GenerateCharmHadron(nu_pdg,Ev); // generate hadron
      double mc  = pdglib->Find(pdg)->Mass();           // lookup mass
 
      LOG("CharmHad", pNOTICE)
          << "Trying charm hadron = " << pdg << "(m = " << mc << ")";
 
      if(mc>=W) continue; // dont' accept
 
      // Generate the charmed hadron energy based on the input
      // fragmentation function
      double z   = fFragmFunc->GenerateZ();  // generate z(=Eh/Ev)
      double Ec  = z*Eh;                     // @ LAB'
      double mc2 = TMath::Power(mc,2);
      double Ec2 = TMath::Power(Ec,2);
      double pc2 = Ec2-mc2;
 
      LOG("CharmHad", pINFO)
          << "Trying charm hadron z = " << z << ", E = " << Ec;
 
      if(pc2<=0) continue;
 
      // Generate the charm hadron pT^2 and pL^2 (with respect to the
      // hadronic system direction @ the LAB)
      double ptc2 = fCharmPT2pdf->GetRandom();
      double plc2 = Ec2 - ptc2 - mc2;
      LOG("CharmHad", pINFO)
            << "Trying charm hadron pT^2 (tranv to pHad) = " << ptc2;
      if(plc2<0) continue;
 
      // Generate the charm hadron momentum components (@ LAB', z:\vec{pHad})
      double ptc = TMath::Sqrt(ptc2);
      double plc = TMath::Sqrt(plc2);
      double phi = (2*kPi) * rnd->RndHadro().Rndm();
      double pxc = ptc * TMath::Cos(phi);
      double pyc = ptc * TMath::Sin(phi);
      double pzc = plc;
 
      p4C.SetPxPyPzE(pxc,pyc,pzc,Ec); // @ LAB'
 
      // Boost charm hadron 4-momentum from the LAB' to the HCM' frame
      //
      LOG("CharmHad", pDEBUG)
        << "Charm hadron p4 (@LAB') = " << utils::print::P4AsString(&p4C);
 
      p4C.Boost(beta);
 
      LOG("CharmHad", pDEBUG)
         << "Charm hadron p4 (@HCM') = " << utils::print::P4AsString(&p4C);
 
      // Hadronic non-charm remnant 4p at HCM'
      TLorentzVector p4 = p4H - p4C;
      double wr = p4.M();
      LOG("CharmHad", pINFO)
              << "Invariant mass of remnant hadronic system= " << wr;
      if(wr < kNucleonMass + kPionMass + kASmallNum) {
         LOG("CharmHad", pINFO) << "Too small hadronic remnant mass!";
         continue;
      }
 
      ch_pdg = pdg;
      got_charmed_hadron = true;
 
      LOG("CharmHad", pNOTICE)
            << "Generated charm hadron = " << pdg << "(m = " <<  mc << ")";
      LOG("CharmHad", pNOTICE)
            << "Generated charm hadron z = " << z << ", E = " << Ec;
   }
 
   // ....................................................................
   // Check whether the code above had difficulty generating the charmed
   // hadron near the W - threshold.
   // If yes, attempt a phase space decay of a low mass charm hadron + nucleon
   // pair that maintains the charge.
   // That's a desperate solution but don't want to quit too early as that
   // would distort the generated dsigma/dW distribution near threshold.
   //
   bool used_lowW_strategy = false;
   int  fs_nucleon_pdg = -1;
   if(ch_pdg==-1 && W < 3.){
      LOG("CharmHad", pNOTICE)
          << "Had difficulty generating charm hadronic system near W threshold";
      LOG("CharmHad", pNOTICE)
          << "Trying an alternative strategy";
 
      double qfsl  = interaction->FSPrimLepton()->Charge() / 3.;
      double qinit = pdglib->Find(nuc_pdg)->Charge() / 3.;
      int qhad  = (int) (qinit - qfsl);
 
      int remn_pdg = -1;
      int chrm_pdg = -1;
 
      //cc-only: qhad(nu) = +1,+2, qhad(nubar)= -1,0
      //
      if(qhad == 2) {
          chrm_pdg = kPdgDP; remn_pdg = kPdgProton;
      } else if(qhad ==  1) {
          if(rnd->RndHadro().Rndm() > 0.5) {
             chrm_pdg = kPdgD0; remn_pdg = kPdgProton;
          } else {
             chrm_pdg = kPdgDP; remn_pdg = kPdgNeutron;
          }
      } else if(qhad ==  0) {
          chrm_pdg = kPdgAntiD0; remn_pdg = kPdgNeutron;
      } else if(qhad == -1) {
          chrm_pdg = kPdgDM; remn_pdg = kPdgNeutron;
      }
 
      double mc  = pdglib->Find(chrm_pdg)->Mass();
      double mn  = pdglib->Find(remn_pdg)->Mass();
 
      if(mc+mn < W) {
         // Set decay
         double mass[2] = {mc, mn};
         bool permitted = fPhaseSpaceGenerator.SetDecay(p4H, 2, mass);
         assert(permitted);
 
         // Get the maximum weight
         double wmax = -1;
         for(int i=0; i<200; i++) {
            double w = fPhaseSpaceGenerator.Generate();
            wmax = TMath::Max(wmax,w);
         }
 
         if(wmax>0) {
            wmax *= 2;
 
            // Generate unweighted decay
            bool accept_decay=false;
            unsigned int idecay_try=0;
            while(!accept_decay)
            {
              idecay_try++;
 
              if(idecay_try>kMaxUnweightDecayIterations) {
                 LOG("CharmHad", pWARN)
                   << "Couldn't generate an unweighted phase space decay after "
                   << idecay_try << " attempts";
              }
              double w  = fPhaseSpaceGenerator.Generate();
              if(w > wmax) {
                 LOG("CharmHad", pWARN)
                  << "Decay weight = " << w << " > max decay weight = " << wmax;
              }
              double gw = wmax * rnd->RndHadro().Rndm();
              accept_decay = (gw<=w);
 
              if(accept_decay) {
                  used_lowW_strategy = true;
                  TLorentzVector * p4 = fPhaseSpaceGenerator.GetDecay(0);
                  p4C            = *p4;
                  ch_pdg         = chrm_pdg;
                  fs_nucleon_pdg = remn_pdg;
              }
            } // decay loop
         }//wmax>0
 
      }// allowed decay
   } // alt low-W strategy
 
   // ....................................................................
   // Check success in generating the charm hadron & compute 4p for
   // remnant system
   //
   if(ch_pdg==-1){
      LOG("CharmHad", pWARN)
          << "Couldn't generate charm hadron for: " << *interaction;
      return 0;
   }
 
   TLorentzVector p4R = p4H - p4C;
   double WR = p4R.M();
   //double MC = pdglib->Find(ch_pdg)->Mass();
 
   LOG("CharmHad", pNOTICE) << "Remnant hadronic system mass = " << WR;
 
   // ....................................................................
   // Handle case where the user doesn't want to remnant system to be
   // hadronized (add as 'hadronic blob')
   //
   if(fCharmOnly) {
     // Create particle list (fragmentation record)
     TClonesArray * particle_list = new TClonesArray("genie::GHepParticle", 2);
     particle_list->SetOwner(true);
 
     // insert the generated particles
     new ((*particle_list)[0]) GHepParticle (ch_pdg,kIStStableFinalState,
              -1,-1,-1,-1, p4C.Px(),p4C.Py(),p4C.Pz(),p4C.E(), 0,0,0,0);
     new ((*particle_list)[1]) GHepParticle (kPdgHadronicBlob,kIStStableFinalState,
              -1,-1,-1,-1, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(), 0,0,0,0);
 
     return particle_list;
   }
 
   // ....................................................................
   // Handle case where the remnant system is already known and doesn't
   // have to be hadronized. That happens when (close to the W threshold)
   // the hadronic system was generated by a simple 2-body decay
   //
   if(used_lowW_strategy) {
     // Create particle list (fragmentation record)
     TClonesArray * particle_list = new TClonesArray("genie::GHepParticle", 3);
     particle_list->SetOwner(true);
 
     // insert the generated particles
     new ((*particle_list)[0]) GHepParticle (ch_pdg,kIStStableFinalState,
             -1,-1,-1,-1, p4C.Px(),p4C.Py(),p4C.Pz(),p4C.E(), 0,0,0,0);
     new ((*particle_list)[1]) GHepParticle (kPdgHadronicBlob,kIStNucleonTarget,
             -1,-1,2,2, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(), 0,0,0,0);
     new ((*particle_list)[2]) GHepParticle (fs_nucleon_pdg,kIStStableFinalState,
             1,1,-1,-1, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(), 0,0,0,0);
 
     return particle_list;
   }
 
   // ....................................................................
   // --------------------------------------------------------------------
   // Hadronize non-charm hadronic blob using PYTHIA/JETSET
   // --------------------------------------------------------------------
   // ....................................................................
 
   // Create output event record
   // Insert the generated charm hadron & the hadronic (non-charm) blob.
   // In this case the hadronic blob is entered as a pre-fragm. state.
 
   TClonesArray * particle_list = new TClonesArray("genie::GHepParticle");
   particle_list->SetOwner(true);
 
   new ((*particle_list)[0]) GHepParticle (ch_pdg,kIStStableFinalState,
           -1,-1,-1,-1,  p4C.Px(),p4C.Py(),p4C.Pz(),p4C.E(), 0,0,0,0);
   new ((*particle_list)[1]) GHepParticle (kPdgHadronicBlob,kIStNucleonTarget,
           -1,-1,2,3, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(), 0,0,0,0);
 
   unsigned int rpos =2; // offset in event record
 
   bool use_pythia = (WR>1.5);
 
 /*
   // Determining quark systems to input to PYTHIA based on simple quark model
   // arguments
   //
   // Neutrinos
   // ------------------------------------------------------------------
   // Scattering off valence q
   // ..................................................................
   // p: [uu]+d
   //         |--> c --> D0       <c+\bar(u)>  : [u]
   //                --> D+       <c+\bar(d)>  : [d]
   //                --> Ds+      <c+\bar(s)>  : [s]
   //                --> Lamda_c+ <c+ud     >  : [\bar(ud)]
   //
   // (for n: [uu] -> 50%[ud]_{0} + 50%[ud]_{1})
   //
   // Scattering off sea q
   // ..................................................................
   // p: [uud] + [\bar(d)]d (or)
   //            [\bar(s)]s
   //                     |--> c --> D0       <c+\bar(u)>  : [u]
   //                            --> D+       <c+\bar(d)>  : [d]
   //                            --> Ds+      <c+\bar(s)>  : [s]
   //                            --> Lamda_c+ <c+ud     >  : [\bar(ud)]
   // Anti-Neutrinos
   // ------------------------------------------------------------------
   // Scattering off sea q
   // ..................................................................
   // p: [uud] + [d] \bar(d) (or)
   //            [s] \bar(s)
   //                   |----> \bar(c) --> \bar(D0) <\bar(c)+u>  : [\bar(u)]
   //                                  --> D-       <\bar(c)+d>  : [\bar(d)]
   //                                  --> Ds-      <\bar(c)+s>  : [\bar(s)]
   // [Summary]
   //                                                                                    Qq
   // | v        + p  [val/d]        -->    D0       +  {           u          uu       }(+2) / u,uu
   // | v        + p  [val/d]        -->    D+       +  {           d          uu       }(+1) / d,uu
   // | v        + p  [val/d]        -->    Ds+      +  {           s          uu       }(+1) / s,uu
   // | v        + p  [val/d]        -->    Lc+      +  {           \bar(ud)   uu       }(+1) / \bar(d),u
   // | v        + n  [val/d]        -->    D0       +  {           u          ud       }(+1) / u,ud
   // | v        + n  [val/d]        -->    D+       +  {           d          ud       }( 0) / d,ud
   // | v        + n  [val/d]        -->    Ds+      +  {           s          ud       }( 0) / s,ud
   // | v        + n  [val/d]        -->    Lc+      +  {           \bar(ud)   ud       }( 0) / \bar(d),d
   // | v        + p  [sea/d]        -->    D0       +  {  uud      \bar(d)    u        }(+2) / u,uu
   // | v        + p  [sea/d]        -->    D+       +  {  uud      \bar(d)    d        }(+1) / d,uu
   // | v        + p  [sea/d]        -->    Ds+      +  {  uud      \bar(d)    s        }(+1) / s,uu
   // | v        + p  [sea/d]        -->    Lc+      +  {  uud      \bar(d)    \bar(ud) }(+1) / \bar(d),u
   // | v        + n  [sea/d]        -->    D0       +  {  udd      \bar(d)    u        }(+1) / u,ud
   // | v        + n  [sea/d]        -->    D+       +  {  udd      \bar(d)    d        }( 0) / d,ud
   // | v        + n  [sea/d]        -->    Ds+      +  {  udd      \bar(d)    s        }( 0) / s,ud
   // | v        + n  [sea/d]        -->    Lc+      +  {  udd      \bar(d)    \bar(ud) }( 0) / \bar(d),d
   // | v        + p  [sea/s]        -->    D0       +  {  uud      \bar(s)    u        }(+2) / u,uu
   // | v        + p  [sea/s]        -->    D+       +  {  uud      \bar(s)    d        }(+1) / d,uu
   // | v        + p  [sea/s]        -->    Ds+      +  {  uud      \bar(s)    s        }(+1) / s,uu
   // | v        + p  [sea/s]        -->    Lc+      +  {  uud      \bar(s)    \bar(ud) }(+1) / \bar(d),u
   // | v        + n  [sea/s]        -->    D0       +  {  udd      \bar(s)    u        }(+1) / u,ud
   // | v        + n  [sea/s]        -->    D+       +  {  udd      \bar(s)    d        }( 0) / d,ud
   // | v        + n  [sea/s]        -->    Ds+      +  {  udd      \bar(s)    s        }( 0) / s,ud
   // | v        + n  [sea/s]        -->    Lc+      +  {  udd      \bar(s)    \bar(ud) }( 0) / \bar(d),d
 
   // | \bar(v)  + p  [sea/\bar(d)]  -->    \bar(D0) +  {  uud      d          \bar(u)  }( 0) / d,ud
   // | \bar(v)  + p  [sea/\bar(d)]  -->    D-       +  {  uud      d          \bar(d)  }(+1) / d,uu
   // | \bar(v)  + p  [sea/\bar(d)]  -->    Ds-      +  {  uud      d          \bar(s)  }(+1) / d,uu
   // | \bar(v)  + n  [sea/\bar(d)]  -->    \bar(D0) +  {  udd      d          \bar(u)  }(-1) / d,dd
   // | \bar(v)  + n  [sea/\bar(d)]  -->    D-       +  {  udd      d          \bar(d)  }( 0) / d,ud
   // | \bar(v)  + n  [sea/\bar(d)]  -->    Ds-      +  {  udd      d          \bar(s)  }( 0) / d,ud
   // | \bar(v)  + p  [sea/\bar(s)]  -->    \bar(D0) +  {  uud      s          \bar(u)  }( 0) / d,ud
   // | \bar(v)  + p  [sea/\bar(s)]  -->    D-       +  {  uud      s          \bar(d)  }(+1) / d,uu
   // | \bar(v)  + p  [sea/\bar(s)]  -->    Ds-      +  {  uud      s          \bar(s)  }(+1) / d,uu
   // | \bar(v)  + n  [sea/\bar(s)]  -->    \bar(D0) +  {  udd      s          \bar(u)  }(-1) / d,dd
   // | \bar(v)  + n  [sea/\bar(s)]  -->    D-       +  {  udd      s          \bar(d)  }( 0) / d,ud
   // | \bar(v)  + n  [sea/\bar(s)]  -->    Ds-      +  {  udd      s          \bar(s)  }( 0) / d,ud
   //
   //
   // Taking some short-cuts below :
   // In the process of obtaining 2 q systems (while conserving the charge) I might tread
   // d,s or \bar(d),\bar(s) as the same
   // In the future I should perform the first steps of the multi-q hadronization manualy
   // (to reduce the number of q's input to PYTHIA) or use py3ent_(), py4ent_() ...
   //
 */
 
   if(use_pythia) {
     int  qrkSyst1 = 0;
     int  qrkSyst2 = 0;
     if(isnu||isdm) { // neutrinos
        if(ch_pdg==kPdgLambdaPc) {
           if(isp) { qrkSyst1 = kPdgAntiDQuark; qrkSyst2 = kPdgUQuark; }
           if(isn) { qrkSyst1 = kPdgAntiDQuark; qrkSyst2 = kPdgDQuark; }
        } else {
           if(isp) { qrkSyst2 = kPdgUUDiquarkS1; }
           if(isn) { qrkSyst2 = (rnd->RndHadro().Rndm()<0.5) ?  kPdgUDDiquarkS0 :  kPdgUDDiquarkS1; }
           if (ch_pdg==kPdgD0      ) { qrkSyst1 = kPdgUQuark; }
           if (ch_pdg==kPdgDP      ) { qrkSyst1 = kPdgDQuark; }
           if (ch_pdg==kPdgDPs     ) { qrkSyst1 = kPdgSQuark; }
        }
      }
      if(isnub) { // antineutrinos
        qrkSyst1 = kPdgDQuark;
        if (isp && ch_pdg==kPdgAntiD0) { qrkSyst2 = (rnd->RndHadro().Rndm()<0.5) ? kPdgUDDiquarkS0 : kPdgUDDiquarkS1; }
        if (isp && ch_pdg==kPdgDM    ) { qrkSyst2 = kPdgUUDiquarkS1; }
        if (isp && ch_pdg==kPdgDMs   ) { qrkSyst2 = kPdgUUDiquarkS1; }
        if (isn && ch_pdg==kPdgAntiD0) { qrkSyst2 = kPdgDDDiquarkS1; }
        if (isn && ch_pdg==kPdgDM    ) { qrkSyst2 = (rnd->RndHadro().Rndm()<0.5) ? kPdgUDDiquarkS0 : kPdgUDDiquarkS1; }
        if (isn && ch_pdg==kPdgDMs   ) { qrkSyst2 = (rnd->RndHadro().Rndm()<0.5) ? kPdgUDDiquarkS0 : kPdgUDDiquarkS1; }
      }
      if(qrkSyst1 == 0 && qrkSyst2 == 0) {
          LOG("CharmHad", pWARN)
              << "Couldn't generate quark systems for PYTHIA in: " << *interaction;
          return 0;
      }
 
      //
      // Run PYTHIA for the hadronization of remnant system
      //
      fPythia->SetMDCY(fPythia->Pycomp(kPdgPi0),              1,0); // don't decay pi0
      fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaM),  1,1); // decay Delta+
      fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_Delta0),  1,1); // decay Delta++
      fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaP),  1,1); // decay Delta++
      fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaPP), 1,1); // decay Delta++
 //   fPythia->SetMDCY(fPythia->Pycomp(kPdgDeltaP),  1,1); // decay Delta+
 //   fPythia->SetMDCY(fPythia->Pycomp(kPdgDeltaPP), 1,1); // decay Delta++
 
      int ip = 0;
      py2ent_(&ip, &qrkSyst1, &qrkSyst2, &WR); // hadronize
 
      fPythia->SetMDCY(fPythia->Pycomp(kPdgPi0),1,1); // restore
 
      //-- Get PYTHIA's LUJETS event record
      TClonesArray * pythia_remnants = 0;
      fPythia->GetPrimaries();
      pythia_remnants = dynamic_cast<TClonesArray *>(fPythia->ImportParticles("All"));
      if(!pythia_remnants) {
          LOG("CharmHad", pWARN) << "Couldn't hadronize (non-charm) remnants!";
          return 0;
       }
 
      int np = pythia_remnants->GetEntries();
      assert(np>0);
 
       // PYTHIA performs the hadronization at the *remnant hadrons* centre of mass
       // frame  (not the hadronic centre of mass frame).
       // Boost all hadronic blob fragments to the HCM', fix their mother/daughter
       // assignments and add them to the fragmentation record.
 
       TVector3 rmnbeta = +1 * p4R.BoostVector(); // boost velocity
 
       TMCParticle * pythia_remn  = 0; // remnant
       GHepParticle * bremn = 0; // boosted remnant
       TIter remn_iter(pythia_remnants);
       while( (pythia_remn = (TMCParticle *) remn_iter.Next()) ) {
 
          // insert and get a pointer to inserted object for mods
         bremn = new ((*particle_list)[rpos++]) GHepParticle ( pythia_remn->GetKF(),                // pdg
                                                               GHepStatus_t(pythia_remn->GetKS()),  // status
                                                               pythia_remn->GetParent(),            // first parent
                                                               -1,                                  // second parent
                                                               pythia_remn->GetFirstChild(),        // first daughter
                                                               pythia_remn->GetLastChild(),         // second daughter
                                                               pythia_remn -> GetPx(),              // px
                                                               pythia_remn -> GetPy(),              // py
                                                               pythia_remn -> GetPz(),              // pz
                                                               pythia_remn -> GetEnergy(),          // e
                                                               pythia_remn->GetVx(),                // x
                                                               pythia_remn->GetVy(),                // y
                                                               pythia_remn->GetVz(),                // z
                                                               pythia_remn->GetTime()               // t
                                                               );
 
         // boost
         bremn -> P4() -> Boost( rmnbeta ) ;
 
          // handle insertion of charmed hadron
          int jp  = bremn->FirstMother();
          int ifc = bremn->FirstDaughter();
          int ilc = bremn->LastDaughter();
 
          bremn -> SetFirstMother( (jp  == 0 ?  1 : jp +1) );
          bremn -> SetFirstDaughter ( (ifc == 0 ? -1 : ifc+1) );
          bremn -> SetLastDaughter  ( (ilc == 0 ? -1 : ilc+1) );
       }
   } // use_pythia
 
   // ....................................................................
   // Hadronizing low-W non-charm hadronic blob using a phase space decay
   // ....................................................................
 
   else {
      // Just a small fraction of events (low-W remnant syste) causing trouble
      // to PYTHIA/JETSET
      // Set a 2-body N+pi system that matches the remnant system charge and
      // do a simple phase space decay
      //
      // q(remn)  remn/syst
      // +2    :  (p pi+)
      // +1    :  50%(p pi0) + 50% n pi+
      //  0    :  50%(p pi-) + 50% n pi0
      // -1    :  (n pi-)
      //
      double qfsl  = interaction->FSPrimLepton()->Charge() / 3.;
      double qinit = pdglib->Find(nuc_pdg)->Charge() / 3.;
      double qch   = pdglib->Find(ch_pdg)->Charge() / 3.;
      int Q = (int) (qinit - qfsl - qch); // remnant hadronic system charge
 
      bool allowdup=true;
      PDGCodeList pd(allowdup);
      if(Q==2) {
             pd.push_back(kPdgProton);  pd.push_back(kPdgPiP);  }
      else if (Q==1) {
        if(rnd->RndHadro().Rndm()<0.5) {
             pd.push_back(kPdgProton);  pd.push_back(kPdgPi0);  }
        else {
             pd.push_back(kPdgNeutron); pd.push_back(kPdgPiP);  }
      }
      else if (Q==0) {
         if(rnd->RndHadro().Rndm()<0.5) {
             pd.push_back(kPdgProton);  pd.push_back(kPdgPiM);  }
         else {
            pd.push_back(kPdgNeutron);  pd.push_back(kPdgPi0);  }
      }
      else if (Q==-1) {
            pd.push_back(kPdgNeutron);  pd.push_back(kPdgPiM);  }
 
      double mass[2] = {
        pdglib->Find(pd[0])->Mass(), pdglib->Find(pd[1])->Mass()
      };
 
      // Set the decay
      bool permitted = fPhaseSpaceGenerator.SetDecay(p4R, 2, mass);
      if(!permitted) {
        LOG("CharmHad", pERROR) << " *** Phase space decay is not permitted";
        return 0;
      }
      // Get the maximum weight
      double wmax = -1;
      for(int i=0; i<200; i++) {
        double w = fPhaseSpaceGenerator.Generate();
        wmax = TMath::Max(wmax,w);
      }
      if(wmax<=0) {
        LOG("CharmHad", pERROR) << " *** Non-positive maximum weight";
        LOG("CharmHad", pERROR) << " *** Can not generate an unweighted phase space decay";
        return 0;
      }
 
      LOG("CharmHad", pINFO)
         << "Max phase space gen. weight @ current hadronic system: " << wmax;
 
       // *** generating an un-weighted decay ***
       wmax *= 1.3;
       bool accept_decay=false;
       unsigned int idectry=0;
       while(!accept_decay)
       {
        idectry++;
        if(idectry>kMaxUnweightDecayIterations) {
          // report, clean-up and return
          LOG("Char,Had", pWARN)
              << "Couldn't generate an unweighted phase space decay after "
              << itry << " attempts";
          return 0;
       }
       double w = fPhaseSpaceGenerator.Generate();
       if(w > wmax) {
           LOG("CharmHad", pWARN)
              << "Decay weight = " << w << " > max decay weight = " << wmax;
        }
        double gw = wmax * rnd->RndHadro().Rndm();
        accept_decay = (gw<=w);
        LOG("CharmHad", pDEBUG)
           << "Decay weight = " << w << " / R = " << gw << " - accepted: " << accept_decay;
      }
      for(unsigned int i=0; i<2; i++) {
         int pdgc = pd[i];
         TLorentzVector * p4d = fPhaseSpaceGenerator.GetDecay(i);
         new ( (*particle_list)[rpos+i] ) GHepParticle(
            pdgc,kIStStableFinalState,1,1,-1,-1,p4d->Px(),p4d->Py(),p4d->Pz(),p4d->Energy(),
            0,0,0,0);
      }
   }
 
   //-- Print & return the fragmentation record
   //utils::fragmrec::Print(particle_list);
   return particle_list;
 }

void AGCharm2019::Initialize ( void ) const

private

Definition at line 82 of file AGCharm2019.cxx.

References fPythia.

Referenced by AGCharm2019().

 {
   fPythia = TPythia6::Instance();
 }

void AGCharm2019::LoadConfig ( void )

private

Definition at line 797 of file AGCharm2019.cxx.

References fCharmOnly, fCharmPT2pdf, fD0BarFrac, fD0FracSpl, fDmFrac, fDpFracSpl, fDsFracSpl, fFracMaxEnergy, fFragmFunc, genie::Algorithm::GetParam(), genie::Algorithm::GetParamDef(), genie::Algorithm::GetParamVect(), genie::controls::kASmallNum, LOG, pFATAL, and genie::Algorithm::SubAlg().

Referenced by Configure().

 {
 
   bool hadronize_remnants ;
   GetParamDef( "HadronizeRemnants", hadronize_remnants, true ) ;
 
   fCharmOnly = ! hadronize_remnants ;
 
   //-- Get a fragmentation function
   fFragmFunc = dynamic_cast<const FragmentationFunctionI *> (
     this->SubAlg("FragmentationFunc"));
   assert(fFragmFunc);
 
   string pt_function ;
   this -> GetParam( "PTFunction", pt_function ) ;
 
   fCharmPT2pdf = new TF1("fCharmPT2pdf", pt_function.c_str(),0,0.6);
 
   // stop ROOT from deleting this object of its own volition
   gROOT->GetListOfFunctions()->Remove(fCharmPT2pdf);
 
   // neutrino charm fractions: D^0, D^+, Ds^+ (remainder: Lamda_c^+)
   std::vector<double> ec, d0frac, dpfrac, dsfrac ;
 
   std::string raw ;
   std::vector<std::string> bits ;
 
   bool invalid_configuration = false ;
 
   // load energy points
   this -> GetParamVect( "CharmFrac-E", ec ) ;
   fFracMaxEnergy = ec.back() - controls::kASmallNum ;
 
   // load D0 fractions
   this -> GetParamVect( "CharmFrac-D0", d0frac ) ;
 
   // check the size
   if ( d0frac.size() != ec.size() ) {
     LOG("AGCharm2019", pFATAL) << "E entries don't match D0 fraction entries";
     LOG("AGCharm2019", pFATAL) << "E:  " << ec.size() ;
     LOG("AGCharm2019", pFATAL) << "D0: " << d0frac.size() ;
     invalid_configuration = true ;
   }
 
   // load D+ fractions
   this -> GetParamVect( "CharmFrac-D+", dpfrac ) ;
 
   // check the size
   if ( dpfrac.size() != ec.size() ) {
     LOG("AGCharm2019", pFATAL) << "E entries don't match D+ fraction entries";
     LOG("AGCharm2019", pFATAL) << "E:  " << ec.size() ;
     LOG("AGCharm2019", pFATAL) << "D+: " << dpfrac.size() ;
     invalid_configuration = true ;
   }
 
     // load D_s fractions
   this -> GetParamVect( "CharmFrac-Ds", dsfrac ) ;
 
   // check the size
   if ( dsfrac.size() != ec.size() ) {
     LOG("AGCharm2019", pFATAL) << "E entries don't match Ds fraction entries";
     LOG("AGCharm2019", pFATAL) << "E:  " << ec.size() ;
     LOG("AGCharm2019", pFATAL) << "Ds: " << dsfrac.size() ;
     invalid_configuration = true ;
   }
 
   fD0FracSpl = new Spline( ec.size(), & ec[0], & d0frac[0] );
   fDpFracSpl = new Spline( ec.size(), & ec[0], & dpfrac[0] );
   fDsFracSpl = new Spline( ec.size(), & ec[0], & dsfrac[0] );
 
   // anti-neutrino charm fractions: bar(D^0), D^-, (remainder: Ds^-)
 
   this -> GetParam( "CharmFrac-D0bar", fD0BarFrac ) ;
   this -> GetParam( "CharmFrac-D-",    fDmFrac ) ;
 
   if ( invalid_configuration ) {
 
     LOG("AGCharm2019", pFATAL)
       << "Invalid configuration: Exiting" ;
 
     // From the FreeBSD Library Functions Manual
     //
     // EX_CONFIG (78)   Something was found in an unconfigured or miscon-
     //                  figured state.
 
     exit( 78 ) ;
   }
 }

void AGCharm2019::ProcessEventRecord ( GHepRecord * event ) const

virtual

Implements genie::EventRecordVisitorI.

Definition at line 87 of file AGCharm2019.cxx.

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

 {
   Interaction * interaction = event->Summary();
   TClonesArray * particle_list = this->Hadronize(interaction);
 
   if(! particle_list ) {
     LOG("AGCharm2019", pWARN) << "Got an empty particle list. Hadronizer failed!";
     LOG("AGCharm2019", pWARN) << "Quitting the current event generation thread";
 
     event->EventFlags()->SetBitNumber(kHadroSysGenErr, true);
 
     genie::exceptions::EVGThreadException exception;
     exception.SetReason("Could not simulate the hadronic system");
     exception.SwitchOnFastForward();
     throw exception;
 
     return;
   }
 
   int mom = event->FinalStateHadronicSystemPosition();
   assert(mom!=-1);
 
   // find the proper status for the particles we are going to put in event record
   bool is_nucleus = interaction->InitState().Tgt().IsNucleus();
   GHepStatus_t istfin = (is_nucleus) ?
     kIStHadronInTheNucleus : kIStStableFinalState ;
 
   // retrieve the hadronic blob lorentz boost
   // Because Hadronize() returned particles not in the LAB reference frame
   const TLorentzVector * had_syst = event -> Particle(mom) -> P4() ;
   TVector3 beta( 0., 0., had_syst -> P()/ had_syst -> Energy() ) ;
 
   // Vector defining rotation from LAB to LAB' (z:= \vec{phad})
   TVector3 unitvq = had_syst -> Vect().Unit();
 
   GHepParticle * neutrino  = event->Probe();
   const TLorentzVector & vtx = *(neutrino->X4());
 
   GHepParticle * particle = 0;
   TIter particle_iter(particle_list);
   while ((particle = (GHepParticle *) particle_iter.Next()))  {
 
     int pdgc = particle -> Pdg() ;
 
     //  bring the particle in the LAB reference frame
     particle -> P4() -> Boost( beta ) ;
     particle -> P4() -> RotateUz( unitvq ) ;
 
     // set the proper status according to a number of things:
     // interaction on a nucleaus or nucleon, particle type
     GHepStatus_t ist = ( particle -> Status() ==1 ) ? istfin : kIStDISPreFragmHadronicState;
 
     // handle gammas, and leptons that might come from internal pythia decays
     // mark them as final state particles
     bool not_hadron = ( pdgc == kPdgGamma ||
                         pdg::IsNeutralLepton(pdgc) ||
                         pdg::IsChargedLepton(pdgc) ) ;
 
     if(not_hadron)  { ist = kIStStableFinalState; }
     particle -> SetStatus( ist ) ;
 
     int im  = mom + 1 + particle -> FirstMother() ;
     int ifc = ( particle -> FirstDaughter() == -1) ? -1 : mom + 1 + particle -> FirstDaughter();
     int ilc = ( particle -> LastDaughter()  == -1) ? -1 : mom + 1 + particle -> LastDaughter();
 
     particle -> SetFirstMother( im ) ;
     if ( ifc > -1 ) {
       particle -> SetFirstDaughter( ifc ) ;
       particle -> SetLastDaughter( ilc ) ;
     }
 
     // the Pythia particle position is overridden
     particle -> SetPosition( vtx ) ;
 
     event->AddParticle(*particle);
   }
 
   particle_list -> Delete() ;
   delete particle_list ;
 
   // update the weight of the event
   event -> SetWeight ( Weight() * event->Weight() );
 
 }