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

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

#include <AGCharm2019.h>

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

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 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
 
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 FragmentationFunctionIfFragmFunc
 charm hadron fragmentation func More...
 
double fFracMaxEnergy
 Maximum energy available for the Meson fractions. More...
 
SplinefD0FracSpl
 nu charm fraction vs Ev: D0 More...
 
SplinefDpFracSpl
 nu charm fraction vs Ev: D+ More...
 
SplinefDsFracSpl
 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)
 
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
 
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...
 
AlgMapfOwnedSubAlgMp
 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().

55  :
56 EventRecordVisitorI("genie::AGCharm2019")
57 {
58  this->Initialize();
59 }
void Initialize(void) const
Definition: AGCharm2019.cxx:82
AGCharm2019::AGCharm2019 ( string  config)

Definition at line 61 of file AGCharm2019.cxx.

References Initialize().

61  :
62 EventRecordVisitorI("genie::AGCharm2019", config)
63 {
64  this->Initialize();
65 }
void Initialize(void) const
Definition: AGCharm2019.cxx:82
AGCharm2019::~AGCharm2019 ( )
virtual

Definition at line 67 of file AGCharm2019.cxx.

References fCharmPT2pdf, fD0FracSpl, fDpFracSpl, and fDsFracSpl.

68 {
69  delete fCharmPT2pdf;
70  fCharmPT2pdf = 0;
71 
72  delete fD0FracSpl;
73  fD0FracSpl = 0;
74 
75  delete fDpFracSpl;
76  fDpFracSpl = 0;
77 
78  delete fDsFracSpl;
79  fDsFracSpl = 0;
80 }
Spline * fDpFracSpl
nu charm fraction vs Ev: D+
Definition: AGCharm2019.h:82
TF1 * fCharmPT2pdf
charm hadron pT^2 pdf
Definition: AGCharm2019.h:76
Spline * fD0FracSpl
nu charm fraction vs Ev: D0
Definition: AGCharm2019.h:81
Spline * fDsFracSpl
nu charm fraction vs Ev: Ds+
Definition: AGCharm2019.h:83

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().

786 {
787  Algorithm::Configure(config);
788  this->LoadConfig();
789 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void LoadConfig(void)
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().

792 {
793  Algorithm::Configure(config);
794  this->LoadConfig();
795 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void LoadConfig(void)
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().

745 {
746  // generate a charmed hadron pdg code using a charm fraction table
747 
748  RandomGen * rnd = RandomGen::Instance();
749  double r = rnd->RndHadro().Rndm();
750 
751  // the ratios are giving up to a certain value.
752  // The rations saturates at high energies, so the values used above that enery
753  // are the evaluated at
754 
755  if(pdg::IsNeutrino(nu_pdg)) {
756 
757  // the ratios are giving up to a certain value.
758  // The rations saturates at high energies, so the values used above that enery
759  // are the evaluated at the maximum energies avaiable for the ratios
760 
761  EvLab = TMath::Min( EvLab, fFracMaxEnergy ) ;
762 
763  double tf = 0;
764  if (r < (tf+=fD0FracSpl->Evaluate(EvLab))) return kPdgD0; // D^0
765  else if (r < (tf+=fDpFracSpl->Evaluate(EvLab))) return kPdgDP; // D^+
766  else if (r < (tf+=fDsFracSpl->Evaluate(EvLab))) return kPdgDPs; // Ds^+
767  else return kPdgLambdaPc; // Lamda_c^+
768 
769  } else if(pdg::IsAntiNeutrino(nu_pdg)) {
770  if (r < fD0BarFrac) return kPdgAntiD0;
771  else if (r < fD0BarFrac+fDmFrac) return kPdgDM;
772  else return kPdgDMs;
773  }
774 
775  LOG("CharmHad", pERROR) << "Could not generate a charm hadron!";
776  return 0;
777 }
const int kPdgAntiD0
Definition: PDGCodes.h:184
const int kPdgDPs
Definition: PDGCodes.h:185
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
#define pERROR
Definition: Messenger.h:59
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
Spline * fDpFracSpl
nu charm fraction vs Ev: D+
Definition: AGCharm2019.h:82
double Evaluate(double x) const
Definition: Spline.cxx:363
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
Spline * fD0FracSpl
nu charm fraction vs Ev: D0
Definition: AGCharm2019.h:81
double fFracMaxEnergy
Maximum energy available for the Meson fractions.
Definition: AGCharm2019.h:79
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
double fD0BarFrac
nubar {D0} charm fraction
Definition: AGCharm2019.h:84
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
const int kPdgLambdaPc
Definition: PDGCodes.h:99
const int kPdgDP
Definition: PDGCodes.h:181
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:53
const int kPdgDMs
Definition: PDGCodes.h:186
const int kPdgDM
Definition: PDGCodes.h:182
const int kPdgD0
Definition: PDGCodes.h:183
double fDmFrac
nubar D- charm fraction
Definition: AGCharm2019.h:85
Spline * fDsFracSpl
nu charm fraction vs Ev: Ds+
Definition: AGCharm2019.h:83
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().

174 {
175  LOG("CharmHad", pNOTICE) << "** Running CHARM hadronizer";
176 
177  PDGLibrary * pdglib = PDGLibrary::Instance();
178  RandomGen * rnd = RandomGen::Instance();
179 
180  // ....................................................................
181  // Get information on the input event
182  //
183  const InitialState & init_state = interaction -> InitState();
184  const Kinematics & kinematics = interaction -> Kine();
185  const Target & target = init_state.Tgt();
186 
187  const TLorentzVector & p4Had = kinematics.HadSystP4();
188 
189  double Ev = init_state.ProbeE(kRfLab);
190  double W = kinematics.W(true);
191 
192  TVector3 beta = -1 * p4Had.BoostVector(); // boost vector for LAB' -> HCM'
193  TLorentzVector p4H(0,0,0,W); // hadronic system 4p @ HCM'
194 
195  double Eh = p4Had.Energy();
196 
197  LOG("CharmHad", pNOTICE) << "Ehad (LAB) = " << Eh << ", W = " << W;
198 
199  int nu_pdg = init_state.ProbePdg();
200  int nuc_pdg = target.HitNucPdg();
201 //int qpdg = target.HitQrkPdg();
202 //bool sea = target.HitSeaQrk();
203  bool isp = pdg::IsProton (nuc_pdg);
204  bool isn = pdg::IsNeutron(nuc_pdg);
205  bool isnu = pdg::IsNeutrino(nu_pdg);
206  bool isnub = pdg::IsAntiNeutrino(nu_pdg);
207  bool isdm = pdg::IsDarkMatter(nu_pdg);
208 
209  // ....................................................................
210  // Attempt to generate a charmed hadron & its 4-momentum
211  //
212 
213  TLorentzVector p4C(0,0,0,0);
214  int ch_pdg = -1;
215 
216  bool got_charmed_hadron = false;
217  unsigned int itry=0;
218 
219  while(itry++ < kRjMaxIterations && !got_charmed_hadron) {
220 
221  // Generate a charmed hadron PDG code
222  int pdg = this->GenerateCharmHadron(nu_pdg,Ev); // generate hadron
223  double mc = pdglib->Find(pdg)->Mass(); // lookup mass
224 
225  LOG("CharmHad", pNOTICE)
226  << "Trying charm hadron = " << pdg << "(m = " << mc << ")";
227 
228  if(mc>=W) continue; // dont' accept
229 
230  // Generate the charmed hadron energy based on the input
231  // fragmentation function
232  double z = fFragmFunc->GenerateZ(); // generate z(=Eh/Ev)
233  double Ec = z*Eh; // @ LAB'
234  double mc2 = TMath::Power(mc,2);
235  double Ec2 = TMath::Power(Ec,2);
236  double pc2 = Ec2-mc2;
237 
238  LOG("CharmHad", pINFO)
239  << "Trying charm hadron z = " << z << ", E = " << Ec;
240 
241  if(pc2<=0) continue;
242 
243  // Generate the charm hadron pT^2 and pL^2 (with respect to the
244  // hadronic system direction @ the LAB)
245  double ptc2 = fCharmPT2pdf->GetRandom();
246  double plc2 = Ec2 - ptc2 - mc2;
247  LOG("CharmHad", pINFO)
248  << "Trying charm hadron pT^2 (tranv to pHad) = " << ptc2;
249  if(plc2<0) continue;
250 
251  // Generate the charm hadron momentum components (@ LAB', z:\vec{pHad})
252  double ptc = TMath::Sqrt(ptc2);
253  double plc = TMath::Sqrt(plc2);
254  double phi = (2*kPi) * rnd->RndHadro().Rndm();
255  double pxc = ptc * TMath::Cos(phi);
256  double pyc = ptc * TMath::Sin(phi);
257  double pzc = plc;
258 
259  p4C.SetPxPyPzE(pxc,pyc,pzc,Ec); // @ LAB'
260 
261  // Boost charm hadron 4-momentum from the LAB' to the HCM' frame
262  //
263  LOG("CharmHad", pDEBUG)
264  << "Charm hadron p4 (@LAB') = " << utils::print::P4AsString(&p4C);
265 
266  p4C.Boost(beta);
267 
268  LOG("CharmHad", pDEBUG)
269  << "Charm hadron p4 (@HCM') = " << utils::print::P4AsString(&p4C);
270 
271  // Hadronic non-charm remnant 4p at HCM'
272  TLorentzVector p4 = p4H - p4C;
273  double wr = p4.M();
274  LOG("CharmHad", pINFO)
275  << "Invariant mass of remnant hadronic system= " << wr;
276  if(wr < kNucleonMass + kPionMass + kASmallNum) {
277  LOG("CharmHad", pINFO) << "Too small hadronic remnant mass!";
278  continue;
279  }
280 
281  ch_pdg = pdg;
282  got_charmed_hadron = true;
283 
284  LOG("CharmHad", pNOTICE)
285  << "Generated charm hadron = " << pdg << "(m = " << mc << ")";
286  LOG("CharmHad", pNOTICE)
287  << "Generated charm hadron z = " << z << ", E = " << Ec;
288  }
289 
290  // ....................................................................
291  // Check whether the code above had difficulty generating the charmed
292  // hadron near the W - threshold.
293  // If yes, attempt a phase space decay of a low mass charm hadron + nucleon
294  // pair that maintains the charge.
295  // That's a desperate solution but don't want to quit too early as that
296  // would distort the generated dsigma/dW distribution near threshold.
297  //
298  bool used_lowW_strategy = false;
299  int fs_nucleon_pdg = -1;
300  if(ch_pdg==-1 && W < 3.){
301  LOG("CharmHad", pNOTICE)
302  << "Had difficulty generating charm hadronic system near W threshold";
303  LOG("CharmHad", pNOTICE)
304  << "Trying an alternative strategy";
305 
306  double qfsl = interaction->FSPrimLepton()->Charge() / 3.;
307  double qinit = pdglib->Find(nuc_pdg)->Charge() / 3.;
308  int qhad = (int) (qinit - qfsl);
309 
310  int remn_pdg = -1;
311  int chrm_pdg = -1;
312 
313  //cc-only: qhad(nu) = +1,+2, qhad(nubar)= -1,0
314  //
315  if(qhad == 2) {
316  chrm_pdg = kPdgDP; remn_pdg = kPdgProton;
317  } else if(qhad == 1) {
318  if(rnd->RndHadro().Rndm() > 0.5) {
319  chrm_pdg = kPdgD0; remn_pdg = kPdgProton;
320  } else {
321  chrm_pdg = kPdgDP; remn_pdg = kPdgNeutron;
322  }
323  } else if(qhad == 0) {
324  chrm_pdg = kPdgAntiD0; remn_pdg = kPdgNeutron;
325  } else if(qhad == -1) {
326  chrm_pdg = kPdgDM; remn_pdg = kPdgNeutron;
327  }
328 
329  double mc = pdglib->Find(chrm_pdg)->Mass();
330  double mn = pdglib->Find(remn_pdg)->Mass();
331 
332  if(mc+mn < W) {
333  // Set decay
334  double mass[2] = {mc, mn};
335  bool permitted = fPhaseSpaceGenerator.SetDecay(p4H, 2, mass);
336  assert(permitted);
337 
338  // Get the maximum weight
339  double wmax = -1;
340  for(int i=0; i<200; i++) {
341  double w = fPhaseSpaceGenerator.Generate();
342  wmax = TMath::Max(wmax,w);
343  }
344 
345  if(wmax>0) {
346  wmax *= 2;
347 
348  // Generate unweighted decay
349  bool accept_decay=false;
350  unsigned int idecay_try=0;
351  while(!accept_decay)
352  {
353  idecay_try++;
354 
355  if(idecay_try>kMaxUnweightDecayIterations) {
356  LOG("CharmHad", pWARN)
357  << "Couldn't generate an unweighted phase space decay after "
358  << idecay_try << " attempts";
359  }
360  double w = fPhaseSpaceGenerator.Generate();
361  if(w > wmax) {
362  LOG("CharmHad", pWARN)
363  << "Decay weight = " << w << " > max decay weight = " << wmax;
364  }
365  double gw = wmax * rnd->RndHadro().Rndm();
366  accept_decay = (gw<=w);
367 
368  if(accept_decay) {
369  used_lowW_strategy = true;
370  TLorentzVector * p4 = fPhaseSpaceGenerator.GetDecay(0);
371  p4C = *p4;
372  ch_pdg = chrm_pdg;
373  fs_nucleon_pdg = remn_pdg;
374  }
375  } // decay loop
376  }//wmax>0
377 
378  }// allowed decay
379  } // alt low-W strategy
380 
381  // ....................................................................
382  // Check success in generating the charm hadron & compute 4p for
383  // remnant system
384  //
385  if(ch_pdg==-1){
386  LOG("CharmHad", pWARN)
387  << "Couldn't generate charm hadron for: " << *interaction;
388  return 0;
389  }
390 
391  TLorentzVector p4R = p4H - p4C;
392  double WR = p4R.M();
393  //double MC = pdglib->Find(ch_pdg)->Mass();
394 
395  LOG("CharmHad", pNOTICE) << "Remnant hadronic system mass = " << WR;
396 
397  // ....................................................................
398  // Handle case where the user doesn't want to remnant system to be
399  // hadronized (add as 'hadronic blob')
400  //
401  if(fCharmOnly) {
402  // Create particle list (fragmentation record)
403  TClonesArray * particle_list = new TClonesArray("genie::GHepParticle", 2);
404  particle_list->SetOwner(true);
405 
406  // insert the generated particles
407  new ((*particle_list)[0]) GHepParticle (ch_pdg,kIStStableFinalState,
408  -1,-1,-1,-1, p4C.Px(),p4C.Py(),p4C.Pz(),p4C.E(), 0,0,0,0);
409  new ((*particle_list)[1]) GHepParticle (kPdgHadronicBlob,kIStStableFinalState,
410  -1,-1,-1,-1, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(), 0,0,0,0);
411 
412  return particle_list;
413  }
414 
415  // ....................................................................
416  // Handle case where the remnant system is already known and doesn't
417  // have to be hadronized. That happens when (close to the W threshold)
418  // the hadronic system was generated by a simple 2-body decay
419  //
420  if(used_lowW_strategy) {
421  // Create particle list (fragmentation record)
422  TClonesArray * particle_list = new TClonesArray("genie::GHepParticle", 3);
423  particle_list->SetOwner(true);
424 
425  // insert the generated particles
426  new ((*particle_list)[0]) GHepParticle (ch_pdg,kIStStableFinalState,
427  -1,-1,-1,-1, p4C.Px(),p4C.Py(),p4C.Pz(),p4C.E(), 0,0,0,0);
428  new ((*particle_list)[1]) GHepParticle (kPdgHadronicBlob,kIStNucleonTarget,
429  -1,-1,2,2, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(), 0,0,0,0);
430  new ((*particle_list)[2]) GHepParticle (fs_nucleon_pdg,kIStStableFinalState,
431  1,1,-1,-1, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(), 0,0,0,0);
432 
433  return particle_list;
434  }
435 
436  // ....................................................................
437  // --------------------------------------------------------------------
438  // Hadronize non-charm hadronic blob using PYTHIA/JETSET
439  // --------------------------------------------------------------------
440  // ....................................................................
441 
442  // Create output event record
443  // Insert the generated charm hadron & the hadronic (non-charm) blob.
444  // In this case the hadronic blob is entered as a pre-fragm. state.
445 
446  TClonesArray * particle_list = new TClonesArray("genie::GHepParticle");
447  particle_list->SetOwner(true);
448 
449  new ((*particle_list)[0]) GHepParticle (ch_pdg,kIStStableFinalState,
450  -1,-1,-1,-1, p4C.Px(),p4C.Py(),p4C.Pz(),p4C.E(), 0,0,0,0);
451  new ((*particle_list)[1]) GHepParticle (kPdgHadronicBlob,kIStNucleonTarget,
452  -1,-1,2,3, p4R.Px(),p4R.Py(),p4R.Pz(),p4R.E(), 0,0,0,0);
453 
454  unsigned int rpos =2; // offset in event record
455 
456  bool use_pythia = (WR>1.5);
457 
458 /*
459  // Determining quark systems to input to PYTHIA based on simple quark model
460  // arguments
461  //
462  // Neutrinos
463  // ------------------------------------------------------------------
464  // Scattering off valence q
465  // ..................................................................
466  // p: [uu]+d
467  // |--> c --> D0 <c+\bar(u)> : [u]
468  // --> D+ <c+\bar(d)> : [d]
469  // --> Ds+ <c+\bar(s)> : [s]
470  // --> Lamda_c+ <c+ud > : [\bar(ud)]
471  //
472  // (for n: [uu] -> 50%[ud]_{0} + 50%[ud]_{1})
473  //
474  // Scattering off sea q
475  // ..................................................................
476  // p: [uud] + [\bar(d)]d (or)
477  // [\bar(s)]s
478  // |--> c --> D0 <c+\bar(u)> : [u]
479  // --> D+ <c+\bar(d)> : [d]
480  // --> Ds+ <c+\bar(s)> : [s]
481  // --> Lamda_c+ <c+ud > : [\bar(ud)]
482  // Anti-Neutrinos
483  // ------------------------------------------------------------------
484  // Scattering off sea q
485  // ..................................................................
486  // p: [uud] + [d] \bar(d) (or)
487  // [s] \bar(s)
488  // |----> \bar(c) --> \bar(D0) <\bar(c)+u> : [\bar(u)]
489  // --> D- <\bar(c)+d> : [\bar(d)]
490  // --> Ds- <\bar(c)+s> : [\bar(s)]
491  // [Summary]
492  // Qq
493  // | v + p [val/d] --> D0 + { u uu }(+2) / u,uu
494  // | v + p [val/d] --> D+ + { d uu }(+1) / d,uu
495  // | v + p [val/d] --> Ds+ + { s uu }(+1) / s,uu
496  // | v + p [val/d] --> Lc+ + { \bar(ud) uu }(+1) / \bar(d),u
497  // | v + n [val/d] --> D0 + { u ud }(+1) / u,ud
498  // | v + n [val/d] --> D+ + { d ud }( 0) / d,ud
499  // | v + n [val/d] --> Ds+ + { s ud }( 0) / s,ud
500  // | v + n [val/d] --> Lc+ + { \bar(ud) ud }( 0) / \bar(d),d
501  // | v + p [sea/d] --> D0 + { uud \bar(d) u }(+2) / u,uu
502  // | v + p [sea/d] --> D+ + { uud \bar(d) d }(+1) / d,uu
503  // | v + p [sea/d] --> Ds+ + { uud \bar(d) s }(+1) / s,uu
504  // | v + p [sea/d] --> Lc+ + { uud \bar(d) \bar(ud) }(+1) / \bar(d),u
505  // | v + n [sea/d] --> D0 + { udd \bar(d) u }(+1) / u,ud
506  // | v + n [sea/d] --> D+ + { udd \bar(d) d }( 0) / d,ud
507  // | v + n [sea/d] --> Ds+ + { udd \bar(d) s }( 0) / s,ud
508  // | v + n [sea/d] --> Lc+ + { udd \bar(d) \bar(ud) }( 0) / \bar(d),d
509  // | v + p [sea/s] --> D0 + { uud \bar(s) u }(+2) / u,uu
510  // | v + p [sea/s] --> D+ + { uud \bar(s) d }(+1) / d,uu
511  // | v + p [sea/s] --> Ds+ + { uud \bar(s) s }(+1) / s,uu
512  // | v + p [sea/s] --> Lc+ + { uud \bar(s) \bar(ud) }(+1) / \bar(d),u
513  // | v + n [sea/s] --> D0 + { udd \bar(s) u }(+1) / u,ud
514  // | v + n [sea/s] --> D+ + { udd \bar(s) d }( 0) / d,ud
515  // | v + n [sea/s] --> Ds+ + { udd \bar(s) s }( 0) / s,ud
516  // | v + n [sea/s] --> Lc+ + { udd \bar(s) \bar(ud) }( 0) / \bar(d),d
517 
518  // | \bar(v) + p [sea/\bar(d)] --> \bar(D0) + { uud d \bar(u) }( 0) / d,ud
519  // | \bar(v) + p [sea/\bar(d)] --> D- + { uud d \bar(d) }(+1) / d,uu
520  // | \bar(v) + p [sea/\bar(d)] --> Ds- + { uud d \bar(s) }(+1) / d,uu
521  // | \bar(v) + n [sea/\bar(d)] --> \bar(D0) + { udd d \bar(u) }(-1) / d,dd
522  // | \bar(v) + n [sea/\bar(d)] --> D- + { udd d \bar(d) }( 0) / d,ud
523  // | \bar(v) + n [sea/\bar(d)] --> Ds- + { udd d \bar(s) }( 0) / d,ud
524  // | \bar(v) + p [sea/\bar(s)] --> \bar(D0) + { uud s \bar(u) }( 0) / d,ud
525  // | \bar(v) + p [sea/\bar(s)] --> D- + { uud s \bar(d) }(+1) / d,uu
526  // | \bar(v) + p [sea/\bar(s)] --> Ds- + { uud s \bar(s) }(+1) / d,uu
527  // | \bar(v) + n [sea/\bar(s)] --> \bar(D0) + { udd s \bar(u) }(-1) / d,dd
528  // | \bar(v) + n [sea/\bar(s)] --> D- + { udd s \bar(d) }( 0) / d,ud
529  // | \bar(v) + n [sea/\bar(s)] --> Ds- + { udd s \bar(s) }( 0) / d,ud
530  //
531  //
532  // Taking some short-cuts below :
533  // In the process of obtaining 2 q systems (while conserving the charge) I might tread
534  // d,s or \bar(d),\bar(s) as the same
535  // In the future I should perform the first steps of the multi-q hadronization manualy
536  // (to reduce the number of q's input to PYTHIA) or use py3ent_(), py4ent_() ...
537  //
538 */
539 
540  if(use_pythia) {
541  int qrkSyst1 = 0;
542  int qrkSyst2 = 0;
543  if(isnu||isdm) { // neutrinos
544  if(ch_pdg==kPdgLambdaPc) {
545  if(isp) { qrkSyst1 = kPdgAntiDQuark; qrkSyst2 = kPdgUQuark; }
546  if(isn) { qrkSyst1 = kPdgAntiDQuark; qrkSyst2 = kPdgDQuark; }
547  } else {
548  if(isp) { qrkSyst2 = kPdgUUDiquarkS1; }
549  if(isn) { qrkSyst2 = (rnd->RndHadro().Rndm()<0.5) ? kPdgUDDiquarkS0 : kPdgUDDiquarkS1; }
550  if (ch_pdg==kPdgD0 ) { qrkSyst1 = kPdgUQuark; }
551  if (ch_pdg==kPdgDP ) { qrkSyst1 = kPdgDQuark; }
552  if (ch_pdg==kPdgDPs ) { qrkSyst1 = kPdgSQuark; }
553  }
554  }
555  if(isnub) { // antineutrinos
556  qrkSyst1 = kPdgDQuark;
557  if (isp && ch_pdg==kPdgAntiD0) { qrkSyst2 = (rnd->RndHadro().Rndm()<0.5) ? kPdgUDDiquarkS0 : kPdgUDDiquarkS1; }
558  if (isp && ch_pdg==kPdgDM ) { qrkSyst2 = kPdgUUDiquarkS1; }
559  if (isp && ch_pdg==kPdgDMs ) { qrkSyst2 = kPdgUUDiquarkS1; }
560  if (isn && ch_pdg==kPdgAntiD0) { qrkSyst2 = kPdgDDDiquarkS1; }
561  if (isn && ch_pdg==kPdgDM ) { qrkSyst2 = (rnd->RndHadro().Rndm()<0.5) ? kPdgUDDiquarkS0 : kPdgUDDiquarkS1; }
562  if (isn && ch_pdg==kPdgDMs ) { qrkSyst2 = (rnd->RndHadro().Rndm()<0.5) ? kPdgUDDiquarkS0 : kPdgUDDiquarkS1; }
563  }
564  if(qrkSyst1 == 0 && qrkSyst2 == 0) {
565  LOG("CharmHad", pWARN)
566  << "Couldn't generate quark systems for PYTHIA in: " << *interaction;
567  return 0;
568  }
569 
570  //
571  // Run PYTHIA for the hadronization of remnant system
572  //
573  fPythia->SetMDCY(fPythia->Pycomp(kPdgPi0), 1,0); // don't decay pi0
574  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaM), 1,1); // decay Delta+
575  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_Delta0), 1,1); // decay Delta++
576  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaP), 1,1); // decay Delta++
577  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaPP), 1,1); // decay Delta++
578 // fPythia->SetMDCY(fPythia->Pycomp(kPdgDeltaP), 1,1); // decay Delta+
579 // fPythia->SetMDCY(fPythia->Pycomp(kPdgDeltaPP), 1,1); // decay Delta++
580 
581  int ip = 0;
582  py2ent_(&ip, &qrkSyst1, &qrkSyst2, &WR); // hadronize
583 
584  fPythia->SetMDCY(fPythia->Pycomp(kPdgPi0),1,1); // restore
585 
586  //-- Get PYTHIA's LUJETS event record
587  TClonesArray * pythia_remnants = 0;
588  fPythia->GetPrimaries();
589  pythia_remnants = dynamic_cast<TClonesArray *>(fPythia->ImportParticles("All"));
590  if(!pythia_remnants) {
591  LOG("CharmHad", pWARN) << "Couldn't hadronize (non-charm) remnants!";
592  return 0;
593  }
594 
595  int np = pythia_remnants->GetEntries();
596  assert(np>0);
597 
598  // PYTHIA performs the hadronization at the *remnant hadrons* centre of mass
599  // frame (not the hadronic centre of mass frame).
600  // Boost all hadronic blob fragments to the HCM', fix their mother/daughter
601  // assignments and add them to the fragmentation record.
602 
603  TVector3 rmnbeta = +1 * p4R.BoostVector(); // boost velocity
604 
605  TMCParticle * pythia_remn = 0; // remnant
606  GHepParticle * bremn = 0; // boosted remnant
607  TIter remn_iter(pythia_remnants);
608  while( (pythia_remn = (TMCParticle *) remn_iter.Next()) ) {
609 
610  // insert and get a pointer to inserted object for mods
611  bremn = new ((*particle_list)[rpos++]) GHepParticle ( pythia_remn->GetKF(), // pdg
612  GHepStatus_t(pythia_remn->GetKS()), // status
613  pythia_remn->GetParent(), // first parent
614  -1, // second parent
615  pythia_remn->GetFirstChild(), // first daughter
616  pythia_remn->GetLastChild(), // second daughter
617  pythia_remn -> GetPx(), // px
618  pythia_remn -> GetPy(), // py
619  pythia_remn -> GetPz(), // pz
620  pythia_remn -> GetEnergy(), // e
621  pythia_remn->GetVx(), // x
622  pythia_remn->GetVy(), // y
623  pythia_remn->GetVz(), // z
624  pythia_remn->GetTime() // t
625  );
626 
627  // boost
628  bremn -> P4() -> Boost( rmnbeta ) ;
629 
630  // handle insertion of charmed hadron
631  int jp = bremn->FirstMother();
632  int ifc = bremn->FirstDaughter();
633  int ilc = bremn->LastDaughter();
634 
635  bremn -> SetFirstMother( (jp == 0 ? 1 : jp +1) );
636  bremn -> SetFirstDaughter ( (ifc == 0 ? -1 : ifc+1) );
637  bremn -> SetLastDaughter ( (ilc == 0 ? -1 : ilc+1) );
638  }
639  } // use_pythia
640 
641  // ....................................................................
642  // Hadronizing low-W non-charm hadronic blob using a phase space decay
643  // ....................................................................
644 
645  else {
646  // Just a small fraction of events (low-W remnant syste) causing trouble
647  // to PYTHIA/JETSET
648  // Set a 2-body N+pi system that matches the remnant system charge and
649  // do a simple phase space decay
650  //
651  // q(remn) remn/syst
652  // +2 : (p pi+)
653  // +1 : 50%(p pi0) + 50% n pi+
654  // 0 : 50%(p pi-) + 50% n pi0
655  // -1 : (n pi-)
656  //
657  double qfsl = interaction->FSPrimLepton()->Charge() / 3.;
658  double qinit = pdglib->Find(nuc_pdg)->Charge() / 3.;
659  double qch = pdglib->Find(ch_pdg)->Charge() / 3.;
660  int Q = (int) (qinit - qfsl - qch); // remnant hadronic system charge
661 
662  bool allowdup=true;
663  PDGCodeList pd(allowdup);
664  if(Q==2) {
665  pd.push_back(kPdgProton); pd.push_back(kPdgPiP); }
666  else if (Q==1) {
667  if(rnd->RndHadro().Rndm()<0.5) {
668  pd.push_back(kPdgProton); pd.push_back(kPdgPi0); }
669  else {
670  pd.push_back(kPdgNeutron); pd.push_back(kPdgPiP); }
671  }
672  else if (Q==0) {
673  if(rnd->RndHadro().Rndm()<0.5) {
674  pd.push_back(kPdgProton); pd.push_back(kPdgPiM); }
675  else {
676  pd.push_back(kPdgNeutron); pd.push_back(kPdgPi0); }
677  }
678  else if (Q==-1) {
679  pd.push_back(kPdgNeutron); pd.push_back(kPdgPiM); }
680 
681  double mass[2] = {
682  pdglib->Find(pd[0])->Mass(), pdglib->Find(pd[1])->Mass()
683  };
684 
685  // Set the decay
686  bool permitted = fPhaseSpaceGenerator.SetDecay(p4R, 2, mass);
687  if(!permitted) {
688  LOG("CharmHad", pERROR) << " *** Phase space decay is not permitted";
689  return 0;
690  }
691  // Get the maximum weight
692  double wmax = -1;
693  for(int i=0; i<200; i++) {
694  double w = fPhaseSpaceGenerator.Generate();
695  wmax = TMath::Max(wmax,w);
696  }
697  if(wmax<=0) {
698  LOG("CharmHad", pERROR) << " *** Non-positive maximum weight";
699  LOG("CharmHad", pERROR) << " *** Can not generate an unweighted phase space decay";
700  return 0;
701  }
702 
703  LOG("CharmHad", pINFO)
704  << "Max phase space gen. weight @ current hadronic system: " << wmax;
705 
706  // *** generating an un-weighted decay ***
707  wmax *= 1.3;
708  bool accept_decay=false;
709  unsigned int idectry=0;
710  while(!accept_decay)
711  {
712  idectry++;
713  if(idectry>kMaxUnweightDecayIterations) {
714  // report, clean-up and return
715  LOG("Char,Had", pWARN)
716  << "Couldn't generate an unweighted phase space decay after "
717  << itry << " attempts";
718  return 0;
719  }
720  double w = fPhaseSpaceGenerator.Generate();
721  if(w > wmax) {
722  LOG("CharmHad", pWARN)
723  << "Decay weight = " << w << " > max decay weight = " << wmax;
724  }
725  double gw = wmax * rnd->RndHadro().Rndm();
726  accept_decay = (gw<=w);
727  LOG("CharmHad", pDEBUG)
728  << "Decay weight = " << w << " / R = " << gw << " - accepted: " << accept_decay;
729  }
730  for(unsigned int i=0; i<2; i++) {
731  int pdgc = pd[i];
732  TLorentzVector * p4d = fPhaseSpaceGenerator.GetDecay(i);
733  new ( (*particle_list)[rpos+i] ) GHepParticle(
734  pdgc,kIStStableFinalState,1,1,-1,-1,p4d->Px(),p4d->Py(),p4d->Pz(),p4d->Energy(),
735  0,0,0,0);
736  }
737  }
738 
739  //-- Print & return the fragmentation record
740  //utils::fragmrec::Print(particle_list);
741  return particle_list;
742 }
const int kPdgAntiD0
Definition: PDGCodes.h:184
const int kPdgP33m1232_DeltaPP
Definition: PDGCodes.h:107
const int kPdgDPs
Definition: PDGCodes.h:185
const int kPdgUUDiquarkS1
Definition: PDGCodes.h:58
double W(bool selected=false) const
Definition: Kinematics.cxx:157
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
#define pERROR
Definition: Messenger.h:59
static const double kNucleonMass
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
int HitNucPdg(void) const
Definition: Target.cxx:304
int FirstDaughter(void) const
Definition: GHepParticle.h:68
const int kPdgHadronicBlob
Definition: PDGCodes.h:211
int GenerateCharmHadron(int nupdg, double EvLab) const
const int kPdgUQuark
Definition: PDGCodes.h:42
TF1 * fCharmPT2pdf
charm hadron pT^2 pdf
Definition: AGCharm2019.h:76
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
TGenPhaseSpace fPhaseSpaceGenerator
a phase space generator
Definition: AGCharm2019.h:71
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:27
static const unsigned int kMaxUnweightDecayIterations
Definition: Controls.h:61
bool IsDarkMatter(int pdgc)
Definition: PDGUtils.cxx:127
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:66
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
const int kPdgSQuark
Definition: PDGCodes.h:46
A list of PDG codes.
Definition: PDGCodeList.h:32
const int kPdgP33m1232_DeltaP
Definition: PDGCodes.h:106
bool fCharmOnly
don&#39;t hadronize non-charm blob
Definition: AGCharm2019.h:75
const int kPdgP33m1232_DeltaM
Definition: PDGCodes.h:104
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
int FirstMother(void) const
Definition: GHepParticle.h:66
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
int LastDaughter(void) const
Definition: GHepParticle.h:69
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const int kPdgUDDiquarkS1
Definition: PDGCodes.h:57
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
const int kPdgAntiDQuark
Definition: PDGCodes.h:45
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
int ProbePdg(void) const
Definition: InitialState.h:64
const int kPdgLambdaPc
Definition: PDGCodes.h:99
void py2ent_(int *, int *, int *, double *)
const int kPdgPiP
Definition: PDGCodes.h:158
const int kPdgPi0
Definition: PDGCodes.h:160
const int kPdgDQuark
Definition: PDGCodes.h:44
static const double kASmallNum
Definition: Controls.h:40
#define pINFO
Definition: Messenger.h:62
const int kPdgUDDiquarkS0
Definition: PDGCodes.h:56
#define pWARN
Definition: Messenger.h:60
const int kPdgP33m1232_Delta0
Definition: PDGCodes.h:105
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
const int kPdgDP
Definition: PDGCodes.h:181
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:53
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
const int kPdgDMs
Definition: PDGCodes.h:186
const FragmentationFunctionI * fFragmFunc
charm hadron fragmentation func
Definition: AGCharm2019.h:77
const int kPdgDM
Definition: PDGCodes.h:182
virtual double GenerateZ(void) const =0
const int kPdgPiM
Definition: PDGCodes.h:159
const int kPdgProton
Definition: PDGCodes.h:81
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
#define pNOTICE
Definition: Messenger.h:61
const Target & Tgt(void) const
Definition: InitialState.h:66
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
double ProbeE(RefFrame_t rf) const
const int kPdgNeutron
Definition: PDGCodes.h:83
const int kPdgD0
Definition: PDGCodes.h:183
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
TPythia6 * fPythia
remnant (non-charm) hadronizer
Definition: AGCharm2019.h:86
enum genie::EGHepStatus GHepStatus_t
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
const int kPdgDDDiquarkS1
Definition: PDGCodes.h:55
void AGCharm2019::Initialize ( void  ) const
private

Definition at line 82 of file AGCharm2019.cxx.

References fPythia.

Referenced by AGCharm2019().

83 {
84  fPythia = TPythia6::Instance();
85 }
TPythia6 * fPythia
remnant (non-charm) hadronizer
Definition: AGCharm2019.h:86
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().

798 {
799 
800  bool hadronize_remnants ;
801  GetParamDef( "HadronizeRemnants", hadronize_remnants, true ) ;
802 
803  fCharmOnly = ! hadronize_remnants ;
804 
805  //-- Get a fragmentation function
806  fFragmFunc = dynamic_cast<const FragmentationFunctionI *> (
807  this->SubAlg("FragmentationFunc"));
808  assert(fFragmFunc);
809 
810  string pt_function ;
811  this -> GetParam( "PTFunction", pt_function ) ;
812 
813  fCharmPT2pdf = new TF1("fCharmPT2pdf", pt_function.c_str(),0,0.6);
814 
815  // stop ROOT from deleting this object of its own volition
816  gROOT->GetListOfFunctions()->Remove(fCharmPT2pdf);
817 
818  // neutrino charm fractions: D^0, D^+, Ds^+ (remainder: Lamda_c^+)
819  std::vector<double> ec, d0frac, dpfrac, dsfrac ;
820 
821  std::string raw ;
822  std::vector<std::string> bits ;
823 
824  bool invalid_configuration = false ;
825 
826  // load energy points
827  this -> GetParamVect( "CharmFrac-E", ec ) ;
828  fFracMaxEnergy = ec.back() - controls::kASmallNum ;
829 
830  // load D0 fractions
831  this -> GetParamVect( "CharmFrac-D0", d0frac ) ;
832 
833  // check the size
834  if ( d0frac.size() != ec.size() ) {
835  LOG("AGCharm2019", pFATAL) << "E entries don't match D0 fraction entries";
836  LOG("AGCharm2019", pFATAL) << "E: " << ec.size() ;
837  LOG("AGCharm2019", pFATAL) << "D0: " << d0frac.size() ;
838  invalid_configuration = true ;
839  }
840 
841  // load D+ fractions
842  this -> GetParamVect( "CharmFrac-D+", dpfrac ) ;
843 
844  // check the size
845  if ( dpfrac.size() != ec.size() ) {
846  LOG("AGCharm2019", pFATAL) << "E entries don't match D+ fraction entries";
847  LOG("AGCharm2019", pFATAL) << "E: " << ec.size() ;
848  LOG("AGCharm2019", pFATAL) << "D+: " << dpfrac.size() ;
849  invalid_configuration = true ;
850  }
851 
852  // load D_s fractions
853  this -> GetParamVect( "CharmFrac-Ds", dsfrac ) ;
854 
855  // check the size
856  if ( dsfrac.size() != ec.size() ) {
857  LOG("AGCharm2019", pFATAL) << "E entries don't match Ds fraction entries";
858  LOG("AGCharm2019", pFATAL) << "E: " << ec.size() ;
859  LOG("AGCharm2019", pFATAL) << "Ds: " << dsfrac.size() ;
860  invalid_configuration = true ;
861  }
862 
863  fD0FracSpl = new Spline( ec.size(), & ec[0], & d0frac[0] );
864  fDpFracSpl = new Spline( ec.size(), & ec[0], & dpfrac[0] );
865  fDsFracSpl = new Spline( ec.size(), & ec[0], & dsfrac[0] );
866 
867  // anti-neutrino charm fractions: bar(D^0), D^-, (remainder: Ds^-)
868 
869  this -> GetParam( "CharmFrac-D0bar", fD0BarFrac ) ;
870  this -> GetParam( "CharmFrac-D-", fDmFrac ) ;
871 
872  if ( invalid_configuration ) {
873 
874  LOG("AGCharm2019", pFATAL)
875  << "Invalid configuration: Exiting" ;
876 
877  // From the FreeBSD Library Functions Manual
878  //
879  // EX_CONFIG (78) Something was found in an unconfigured or miscon-
880  // figured state.
881 
882  exit( 78 ) ;
883  }
884 }
Spline * fDpFracSpl
nu charm fraction vs Ev: D+
Definition: AGCharm2019.h:82
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:58
#define pFATAL
Definition: Messenger.h:56
int GetParamVect(const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
Handle to load vectors of parameters.
TF1 * fCharmPT2pdf
charm hadron pT^2 pdf
Definition: AGCharm2019.h:76
Spline * fD0FracSpl
nu charm fraction vs Ev: D0
Definition: AGCharm2019.h:81
bool fCharmOnly
don&#39;t hadronize non-charm blob
Definition: AGCharm2019.h:75
double fFracMaxEnergy
Maximum energy available for the Meson fractions.
Definition: AGCharm2019.h:79
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
double fD0BarFrac
nubar {D0} charm fraction
Definition: AGCharm2019.h:84
static const double kASmallNum
Definition: Controls.h:40
Pure abstract base class. Defines the FragmentationFunctionI interface to be implemented by any algor...
const FragmentationFunctionI * fFragmFunc
charm hadron fragmentation func
Definition: AGCharm2019.h:77
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 fDmFrac
nubar D- charm fraction
Definition: AGCharm2019.h:85
Spline * fDsFracSpl
nu charm fraction vs Ev: Ds+
Definition: AGCharm2019.h:83
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
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().

88 {
89  Interaction * interaction = event->Summary();
90  TClonesArray * particle_list = this->Hadronize(interaction);
91 
92  if(! particle_list ) {
93  LOG("AGCharm2019", pWARN) << "Got an empty particle list. Hadronizer failed!";
94  LOG("AGCharm2019", pWARN) << "Quitting the current event generation thread";
95 
96  event->EventFlags()->SetBitNumber(kHadroSysGenErr, true);
97 
99  exception.SetReason("Could not simulate the hadronic system");
100  exception.SwitchOnFastForward();
101  throw exception;
102 
103  return;
104  }
105 
106  int mom = event->FinalStateHadronicSystemPosition();
107  assert(mom!=-1);
108 
109  // find the proper status for the particles we are going to put in event record
110  bool is_nucleus = interaction->InitState().Tgt().IsNucleus();
111  GHepStatus_t istfin = (is_nucleus) ?
113 
114  // retrieve the hadronic blob lorentz boost
115  // Because Hadronize() returned particles not in the LAB reference frame
116  const TLorentzVector * had_syst = event -> Particle(mom) -> P4() ;
117  TVector3 beta( 0., 0., had_syst -> P()/ had_syst -> Energy() ) ;
118 
119  // Vector defining rotation from LAB to LAB' (z:= \vec{phad})
120  TVector3 unitvq = had_syst -> Vect().Unit();
121 
122  GHepParticle * neutrino = event->Probe();
123  const TLorentzVector & vtx = *(neutrino->X4());
124 
125  GHepParticle * particle = 0;
126  TIter particle_iter(particle_list);
127  while ((particle = (GHepParticle *) particle_iter.Next())) {
128 
129  int pdgc = particle -> Pdg() ;
130 
131  // bring the particle in the LAB reference frame
132  particle -> P4() -> Boost( beta ) ;
133  particle -> P4() -> RotateUz( unitvq ) ;
134 
135  // set the proper status according to a number of things:
136  // interaction on a nucleaus or nucleon, particle type
137  GHepStatus_t ist = ( particle -> Status() ==1 ) ? istfin : kIStDISPreFragmHadronicState;
138 
139  // handle gammas, and leptons that might come from internal pythia decays
140  // mark them as final state particles
141  bool not_hadron = ( pdgc == kPdgGamma ||
142  pdg::IsNeutralLepton(pdgc) ||
143  pdg::IsChargedLepton(pdgc) ) ;
144 
145  if(not_hadron) { ist = kIStStableFinalState; }
146  particle -> SetStatus( ist ) ;
147 
148  int im = mom + 1 + particle -> FirstMother() ;
149  int ifc = ( particle -> FirstDaughter() == -1) ? -1 : mom + 1 + particle -> FirstDaughter();
150  int ilc = ( particle -> LastDaughter() == -1) ? -1 : mom + 1 + particle -> LastDaughter();
151 
152  particle -> SetFirstMother( im ) ;
153  if ( ifc > -1 ) {
154  particle -> SetFirstDaughter( ifc ) ;
155  particle -> SetLastDaughter( ilc ) ;
156  }
157 
158  // the Pythia particle position is overridden
159  particle -> SetPosition( vtx ) ;
160 
161  event->AddParticle(*particle);
162  }
163 
164  particle_list -> Delete() ;
165  delete particle_list ;
166 
167  // update the weight of the event
168  event -> SetWeight ( Weight() * event->Weight() );
169 
170 }
bool IsNucleus(void) const
Definition: Target.cxx:272
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:101
virtual double Weight(void) const
Definition: GHepRecord.h:124
Summary information for an interaction.
Definition: Interaction.h:56
TClonesArray * Hadronize(const Interaction *) const
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
double Weight(void) const
#define pWARN
Definition: Messenger.h:60
bool IsNeutralLepton(int pdgc)
Definition: PDGUtils.cxx:95
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
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
enum genie::EGHepStatus GHepStatus_t
double AGCharm2019::Weight ( void  ) const
private

Definition at line 779 of file AGCharm2019.cxx.

Referenced by ProcessEventRecord().

780 {
781  return 1. ;
782 }

Member Data Documentation

bool genie::AGCharm2019::fCharmOnly
private

don't hadronize non-charm blob

Definition at line 75 of file AGCharm2019.h.

Referenced by Hadronize(), and LoadConfig().

TF1* genie::AGCharm2019::fCharmPT2pdf
private

charm hadron pT^2 pdf

Definition at line 76 of file AGCharm2019.h.

Referenced by Hadronize(), LoadConfig(), and ~AGCharm2019().

double genie::AGCharm2019::fD0BarFrac
private

nubar {D0} charm fraction

Definition at line 84 of file AGCharm2019.h.

Referenced by GenerateCharmHadron(), and LoadConfig().

Spline* genie::AGCharm2019::fD0FracSpl
private

nu charm fraction vs Ev: D0

Definition at line 81 of file AGCharm2019.h.

Referenced by GenerateCharmHadron(), LoadConfig(), and ~AGCharm2019().

double genie::AGCharm2019::fDmFrac
private

nubar D- charm fraction

Definition at line 85 of file AGCharm2019.h.

Referenced by GenerateCharmHadron(), and LoadConfig().

Spline* genie::AGCharm2019::fDpFracSpl
private

nu charm fraction vs Ev: D+

Definition at line 82 of file AGCharm2019.h.

Referenced by GenerateCharmHadron(), LoadConfig(), and ~AGCharm2019().

Spline* genie::AGCharm2019::fDsFracSpl
private

nu charm fraction vs Ev: Ds+

Definition at line 83 of file AGCharm2019.h.

Referenced by GenerateCharmHadron(), LoadConfig(), and ~AGCharm2019().

double genie::AGCharm2019::fFracMaxEnergy
private

Maximum energy available for the Meson fractions.

Definition at line 79 of file AGCharm2019.h.

Referenced by GenerateCharmHadron(), and LoadConfig().

const FragmentationFunctionI* genie::AGCharm2019::fFragmFunc
private

charm hadron fragmentation func

Definition at line 77 of file AGCharm2019.h.

Referenced by Hadronize(), and LoadConfig().

TGenPhaseSpace genie::AGCharm2019::fPhaseSpaceGenerator
mutableprivate

a phase space generator

Definition at line 71 of file AGCharm2019.h.

Referenced by Hadronize().

TPythia6* genie::AGCharm2019::fPythia
mutableprivate

remnant (non-charm) hadronizer

Definition at line 86 of file AGCharm2019.h.

Referenced by Hadronize(), and Initialize().


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