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

Event generator for SuSAv2 1p1h interactions. More...

#include <QELEventGeneratorSuSA.h>

Inheritance diagram for genie::QELEventGeneratorSuSA:
genie::KineGeneratorWithCache genie::EventRecordVisitorI genie::Algorithm

Public Member Functions

 QELEventGeneratorSuSA ()
 
 QELEventGeneratorSuSA (string config)
 
 ~QELEventGeneratorSuSA ()
 
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 AddTargetNucleusRemnant (GHepRecord *event) const
 
void SelectLeptonKinematics (GHepRecord *event) const
 
void GenerateNucleon (GHepRecord *event) const
 
double ComputeMaxXSec (const Interaction *in) const
 

Private Attributes

const NuclearModelIfNuclModel
 
double fQ3Max
 
bool fForceBound
 
bool fForceEbFromModel
 
bool fForceFixEb
 
double fEbOR
 
double fEbC
 
const EventRecordVisitorIfFreeNucleonEventGenerator
 

Additional Inherited Members

- Static Public Member Functions inherited from genie::Algorithm
static string BuildParamVectKey (const std::string &comm_name, unsigned int i)
 
static string BuildParamVectSizeKey (const std::string &comm_name)
 
- Protected Member Functions inherited from genie::KineGeneratorWithCache
 KineGeneratorWithCache ()
 
 KineGeneratorWithCache (string name)
 
 KineGeneratorWithCache (string name, string config)
 
 ~KineGeneratorWithCache ()
 
virtual double MaxXSec (GHepRecord *evrec) const
 
virtual double FindMaxXSec (const Interaction *in) const
 
virtual void CacheMaxXSec (const Interaction *in, double xsec) const
 
virtual double Energy (const Interaction *in) const
 
virtual CacheBranchFxAccessCacheBranch (const Interaction *in) const
 
virtual void AssertXSecLimits (const Interaction *in, double xsec, double xsec_max) const
 
- Protected Member Functions inherited from genie::EventRecordVisitorI
 EventRecordVisitorI ()
 
 EventRecordVisitorI (string name)
 
 EventRecordVisitorI (string name, string config)
 
- Protected Member Functions inherited from genie::Algorithm
 Algorithm ()
 
 Algorithm (string name)
 
 Algorithm (string name, string config)
 
void Initialize (void)
 
void DeleteConfig (void)
 
void DeleteSubstructure (void)
 
RegistryExtractLocalConfig (const Registry &in) const
 
RegistryExtractLowerConfig (const Registry &in, const string &alg_key) const
 Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key. More...
 
template<class T >
bool GetParam (const RgKey &name, T &p, bool is_top_call=true) const
 
template<class T >
bool GetParamDef (const RgKey &name, T &p, const T &def) const
 
template<class T >
int GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
 Handle to load vectors of parameters. More...
 
int GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
 
int AddTopRegistry (Registry *rp, bool owns=true)
 add registry with top priority, also update ownership More...
 
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership More...
 
int MergeTopRegistry (const Registry &r)
 
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships. More...
 
- Protected Attributes inherited from genie::KineGeneratorWithCache
const XSecAlgorithmIfXSecModel
 
double fSafetyFactor
 maxxsec -> maxxsec * safety_factor More...
 
double fMaxXSecDiffTolerance
 max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec>maxxsec More...
 
double fEMin
 min E for which maxxsec is cached - forcing explicit calc. More...
 
bool fGenerateUniformly
 uniform over allowed phase space + event weight? More...
 
- Protected Attributes inherited from genie::Algorithm
bool fAllowReconfig
 
bool fOwnsSubstruc
 true if it owns its substructure (sub-algs,...) More...
 
AlgId fID
 algorithm name and configuration set More...
 
vector< Registry * > fConfVect
 
vector< boolfOwnerships
 ownership for every registry in fConfVect More...
 
AlgStatus_t fStatus
 algorithm execution status More...
 
AlgMapfOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool) More...
 

Detailed Description

Event generator for SuSAv2 1p1h interactions.

Author
Stephen Dolan <stephen.joseph.dolan cern.ch> European Organization for Nuclear Research (CERN)

Steven Gardiner <gardiner fnal.gov> Fermi National Accelerator Laboratory

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

Definition at line 32 of file QELEventGeneratorSuSA.h.

Constructor & Destructor Documentation

QELEventGeneratorSuSA::QELEventGeneratorSuSA ( )

Definition at line 51 of file QELEventGeneratorSuSA.cxx.

51  :
52 KineGeneratorWithCache("genie::QELEventGeneratorSuSA")
53 {
54 
55 }
QELEventGeneratorSuSA::QELEventGeneratorSuSA ( string  config)

Definition at line 57 of file QELEventGeneratorSuSA.cxx.

57  :
58 KineGeneratorWithCache("genie::QELEventGeneratorSuSA", config)
59 {
60 
61 }
static Config * config
Definition: config.cpp:1054
QELEventGeneratorSuSA::~QELEventGeneratorSuSA ( )

Definition at line 63 of file QELEventGeneratorSuSA.cxx.

64 {
65 
66 }

Member Function Documentation

void QELEventGeneratorSuSA::AddTargetNucleusRemnant ( GHepRecord event) const
private

Definition at line 287 of file QELEventGeneratorSuSA.cxx.

288 {
289  // add the remnant nuclear target at the GHEP record
290 
291  LOG("QELEvent", pINFO) << "Adding final state nucleus";
292 
293  double Px = 0;
294  double Py = 0;
295  double Pz = 0;
296  double E = 0;
297 
298  GHepParticle * nucleus = event->TargetNucleus();
299  bool have_nucleus = nucleus != 0;
300  if (!have_nucleus) return;
301 
302  int A = nucleus->A();
303  int Z = nucleus->Z();
304 
305  int fd = nucleus->FirstDaughter();
306  int ld = nucleus->LastDaughter();
307 
308  for(int id = fd; id <= ld; id++) {
309 
310  // compute A,Z for final state nucleus & get its PDG code and its mass
311  GHepParticle * particle = event->Particle(id);
312  assert(particle);
313  int pdgc = particle->Pdg();
314  bool is_p = pdg::IsProton (pdgc);
315  bool is_n = pdg::IsNeutron(pdgc);
316 
317  if (is_p) Z--;
318  if (is_p || is_n) A--;
319 
320  Px += particle->Px();
321  Py += particle->Py();
322  Pz += particle->Pz();
323  E += particle->E();
324 
325  }//daughters
326 
327  TParticlePDG * remn = 0;
328  int ipdgc = pdg::IonPdgCode(A, Z);
329  remn = PDGLibrary::Instance()->Find(ipdgc);
330  if(!remn) {
331  LOG("HadronicVtx", pFATAL)
332  << "No particle with [A = " << A << ", Z = " << Z
333  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
334  assert(remn);
335  }
336 
337  double Mi = nucleus->Mass();
338  Px *= -1;
339  Py *= -1;
340  Pz *= -1;
341  E = Mi-E;
342 
343  // Add the nucleus to the event record
344  LOG("QELEvent", pINFO)
345  << "Adding nucleus [A = " << A << ", Z = " << Z
346  << ", pdgc = " << ipdgc << "]";
347 
348  int imom = event->TargetNucleusPosition();
349  event->AddParticle(
350  ipdgc,kIStStableFinalState, imom,-1,-1,-1, Px,Py,Pz,E, 0,0,0,0);
351 
352  LOG("QELEvent", pINFO) << "Done";
353  LOG("QELEvent", pINFO) << *event;
354 }
int Z(void) const
double E(void) const
Get energy.
Definition: GHepParticle.h:91
int FirstDaughter(void) const
Definition: GHepParticle.h:68
#define pFATAL
Definition: Messenger.h:56
double Mass(void) const
Mass that corresponds to the PDG code.
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
int Pdg(void) const
Definition: GHepParticle.h:63
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:338
int LastDaughter(void) const
Definition: GHepParticle.h:69
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:333
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
#define pINFO
Definition: Messenger.h:62
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:57
E
Definition: 018_def.c:13
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:68
#define A
Definition: memgrp.cpp:38
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:75
int A(void) const
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
double Py(void) const
Get Py.
Definition: GHepParticle.h:89
double QELEventGeneratorSuSA::ComputeMaxXSec ( const Interaction in) const
privatevirtual

Implements genie::KineGeneratorWithCache.

Definition at line 643 of file QELEventGeneratorSuSA.cxx.

645 {
646  // Computes the maximum differential cross section in the requested phase
647  // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
648  // method and the value is cached at a circular cache branch for retrieval
649  // during subsequent event generation.
650  // The computed max differential cross section does not need to be the exact
651  // maximum. The number used in the rejection method will be scaled up by a
652  // safety factor. But it needs to be fast - do not use very small steps.
653 
654  double max_xsec = utils::mec::GetMaxXSecTlctl( *fXSecModel, *interaction );
655 
656  // Apply safety factor, since value retrieved from the cache might
657  // correspond to a slightly different energy
658  max_xsec *= fSafetyFactor;
659 
660  #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
661  SLOG("QELKinematics", pDEBUG) << interaction->AsString();
662  SLOG("QELKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec;
663  SLOG("QELKinematics", pDEBUG) << "Computed using alg = " << *fXSecModel;
664  #endif
665 
666  return max_xsec;
667 }
double fSafetyFactor
maxxsec -> maxxsec * safety_factor
double GetMaxXSecTlctl(const XSecAlgorithmI &xsec_model, const Interaction &inter, const double tolerance=0.01, const double safety_factor=1.2, const int max_n_layers=100)
Definition: MECUtils.cxx:334
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:84
#define pDEBUG
Definition: Messenger.h:63
void QELEventGeneratorSuSA::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 582 of file QELEventGeneratorSuSA.cxx.

583 {
584  Algorithm::Configure(config);
585  this->LoadConfig();
586 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void QELEventGeneratorSuSA::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 588 of file QELEventGeneratorSuSA.cxx.

589 {
591  this->LoadConfig();
592 }
static Config * config
Definition: config.cpp:1054
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void QELEventGeneratorSuSA::GenerateNucleon ( GHepRecord event) const
private

Definition at line 356 of file QELEventGeneratorSuSA.cxx.

357 {
358  // We need a kinematic limits accept/reject loop here, so generating the
359  // initial hadrons is combined with generating the recoil hadrons...
360 
361  LOG("QELEvent",pDEBUG) << "Generate Nucleon - Start";
362 
363  Interaction * interaction = event->Summary();
364  GHepParticle * neutrino = event->Probe();
365  assert(neutrino);
366  TLorentzVector p4nu(*neutrino->P4());
367 
368  // get final state primary lepton & its 4-momentum
369  GHepParticle * fsl = event->FinalStatePrimaryLepton();
370  assert(fsl);
371  TLorentzVector p4l(*fsl->P4());
372 
373  // calculate 4-momentum transfer from these
374  TLorentzVector Q4 = p4nu - p4l;
375  //Q4.Print();
376 
377  // get the target nucleon and nucleus.
378  // the remnant nucleus is apparently set, except for its momentum.
379  GHepParticle * target_nucleus = event->TargetNucleus();
380  bool have_nucleus = (target_nucleus != 0);
381  //assert(target_nucleus);
382  GHepParticle * initial_nucleon = event->HitNucleon();
383  assert(initial_nucleon);
384  GHepParticle * remnant_nucleus = event->RemnantNucleus();
385  //assert(remnant_nucleus);
386 
387  // instantiate an empty local target nucleus, so I can use existing methods
388  // to get a momentum from the prevailing Fermi-motion distribution
389  int tgtpdg;
390  if(have_nucleus) tgtpdg = target_nucleus->Pdg();
391  else tgtpdg = kPdgTgtFreeP;
392  Target tgt(tgtpdg);
393 
394  // These things need to be saved through to the end of the accept loop.
395  bool accept = false;
396  TVector3 p3i;
397  unsigned int iter = 0;
398 
399  int initial_nucleon_pdg = initial_nucleon->Pdg();
400  int final_nucleon_pdg = interaction->RecoilNucleonPdg();
401 
402  TLorentzVector p4initial_nucleon;
403  TLorentzVector p4final_nucleon;
404  double removalenergy;
405 
406  //remnant kinematic alterations
407  double pxb = 0;
408  double pyb = 0;
409  double pzb = 0;
410 
411  //===========================================================================
412  // Choose nucleons from the prevailing fermi-motion distribution.
413  // Possible to produce kinematically unallowed (Pauli blocked) nucleon.
414  // Find out, and if so choose them again with this accept/reject loop.
415  // Pauli blocking was included in the SuSAv2 tensor tables, so it should not
416  // be allowed to affect the inclusive xsec.
417  while(!accept){
418  iter++;
419  if(iter > kRjMaxIterations) {
420  // error if try too many times
421  LOG("QELEvent", pWARN)
422  << "Couldn't select a valid nucleon after "
423  << iter << " iterations";
424  event->EventFlags()->SetBitNumber(kKineGenErr, true);
426  exception.SetReason("Couldn't select initial hadron kinematics");
427  exception.SwitchOnFastForward();
428  throw exception;
429  }
430 
431  // generate nucleons
432  // tgt is a Target object for local use, just waiting to be filled.
433  // this sets the pdg of each nucleon and its momentum from user chosen nuclear model
434 
435  double hitNucPos = initial_nucleon->X4()->Vect().Mag();
436  tgt.SetHitNucPdg(initial_nucleon_pdg);
437  fNuclModel->GenerateNucleon(tgt,hitNucPos);
438  p3i = fNuclModel->Momentum3();
439 
440  // Default: Calculate the removal energy as in Guille's thesis - this is a simplicification of
441  // a fairly complex aproach employed in SuSAv2, but we expect it to work pretty well.
442  // We should write something about this in the implementation technical paper ...
443  // IMPORTANT CAVEAT: By default we choose to allow the binding energy to depend on the interaction
444  // (as it should), but this means we don't corrolate the chosen Eb with the intial nucleon
445  // momentum. Therefore we can sometimes have initial state nucleons with KE>Eb. This isn't
446  // great, but fot the moment it's what we have (working on improvments).
447 
448  double q3 = Q4.Vect().Mag();
449 
450  if(!have_nucleus){
451  // For elastic case no Fermi momentum and no Eb
452  removalenergy = 0.0;
453  p3i.SetXYZ(0.0,0.0,0.0);
454  }
455  else if(fForceEbFromModel) removalenergy = fNuclModel->RemovalEnergy();
456  else if(fForceFixEb) removalenergy = fEbOR;
457  else{
458  if(q3<0.827){
459  removalenergy = -0.017687 + 0.0564*q3;
460  }
461  else{
462  removalenergy = 0.0289558;
463  }
464  // The above is for Carbon, should shift for different targets
465  removalenergy += (fNuclModel->RemovalEnergy()) - fEbC;
466  if(removalenergy<0.005) removalenergy=0.005;
467  }
468 
469  //removalenergy=0.0;
470  //initial_nucleon->SetRemovalEnergy(removalenergy);
471 
472  LOG("QELEvent",pDEBUG) << "q3 for this event:" << q3;
473  LOG("QELEvent",pDEBUG) << "Removal energy:" << removalenergy;
474 
475  // Now write down the initial nucleon four-vector for this choice
476  double mass = PDGLibrary::Instance()->Find( initial_nucleon_pdg )->Mass();
477  double mass2 = mass*mass;
478  double energy = TMath::Sqrt(p3i.Mag2() + mass2);
479  p4initial_nucleon.SetPxPyPzE(p3i.Px(),p3i.Py(),p3i.Pz(),energy);
480 
481  //One rather unsubtle option for making sure nucleons remain bound.
482  //This will give us bound nucleons with a sensible missing eneergy
483  //but the distribution of Fermi motion will look crazy.
484  // Anything we do is wrong (without semi-inclusive inputs) you just
485  // have to decide what is less wrong!
486  if(fForceBound && (energy-mass>removalenergy)) continue;
487 
488  // cast the removal energy as the energy component of a 4-vector for later.
489  TLorentzVector tLVebind(0., 0., 0., -1.0 * (removalenergy));
490 
491  // Taken from MEC event generator:
492  // calculate recoil nucleon cluster 4-momentum (tLVebind is negative)
493  p4final_nucleon = p4initial_nucleon + Q4 + tLVebind;
494 
495  // Put on shell as in the Aggregator
496  // This is a bit of a horrible approximation but it is hard to think
497  // of anything simple and better without semi-inclusive model predictions.
498  // However, we are working on an improvments.
499  if(have_nucleus){
500  double En = p4final_nucleon.E();
501  double M = PDGLibrary::Instance()->Find( final_nucleon_pdg )->Mass();
502  double pmag_old = p4final_nucleon.Vect().Mag();
503 
504  double pmag_new = TMath::Sqrt(utils::math::NonNegative(En*En-M*M));
505  double scale = pmag_new / pmag_old;
506 
507  double pxn = scale * p4final_nucleon.Px();
508  double pyn = scale * p4final_nucleon.Py();
509  double pzn = scale * p4final_nucleon.Pz();
510 
511  p4final_nucleon.SetPxPyPzE(pxn,pyn,pzn,En);
512  }
513 
514  // Extra momentum xfer to keep nucleon on shell - I'll give this to the remnant
515 
516  pxb = p4nu.Px()-p4l.Px()-p4final_nucleon.Px();
517  pyb = p4nu.Py()-p4l.Py()-p4final_nucleon.Py();
518  pzb = p4nu.Pz()-p4l.Pz()-p4final_nucleon.Pz();
519 
520  LOG("QELEvent",pDEBUG) << "Remnant momentum is: (" << pxb << ", " << pyb << ", " << pzb << ")";
521 
522 
523  // Pauli blocking check:
524  // Test if the resulting four-vector corresponds to a high-enough energy.
525  // Fail the accept if we couldn't put this thing on-shell.
526  // Basically: is energy of the nucleon positive after we subtracting Eb
527  if (p4final_nucleon.E() < PDGLibrary::Instance()->Find(final_nucleon_pdg)->Mass()) {
528  accept = false;
529  LOG("QELEvent",pDEBUG) << "Rejected nucleon, can't be put on-shell";
530  LOG("QELEvent",pDEBUG) << "Nucleon invariant mass:" << p4final_nucleon.M();
531  LOG("QELEvent",pDEBUG) << "Nucleon real mass:" << PDGLibrary::Instance()->Find(final_nucleon_pdg)->Mass();
532  LOG("QELEvent",pDEBUG) << "Nucleon 4 momenutum:";
533  //p4final_nucleon.Print();
534  LOG("QELEvent",pDEBUG) << "Removal energy:" << removalenergy;
535  LOG("QELEvent",pDEBUG) << "Q4 transfer:";
536  //Q4.Print();
537  }
538  else {
539  accept = true;
540  LOG("QELEvent",pDEBUG) << "Nucleon accepted, Q4 is";
541  //Q4.Print();
542  LOG("QELEvent",pDEBUG) << "Initial nucleon mass is" << sqrt((p4initial_nucleon.E()*p4initial_nucleon.E())-(p4initial_nucleon.Vect().Mag()*p4initial_nucleon.Vect().Mag()));
543  LOG("QELEvent",pDEBUG) << "Final nucleon mass is" << sqrt((p4final_nucleon.E()*p4final_nucleon.E())-(p4final_nucleon.Vect().Mag()*p4final_nucleon.Vect().Mag()));
544  LOG("QELEvent",pDEBUG) << "Nucleon invariant mass:" << p4final_nucleon.M();
545  LOG("QELEvent",pDEBUG) << "Nucleon real mass:" << PDGLibrary::Instance()->Find(final_nucleon_pdg)->Mass();
546  LOG("QELEvent",pDEBUG) << "Nucleon 4 momenutum:";
547  //p4final_nucleon.Print();
548  LOG("QELEvent",pDEBUG) << "Removal energy:" << removalenergy;
549  LOG("QELEvent",pDEBUG) << "Q4 transfer:";
550  }
551 
552  } // end accept loop
553 
554  // we got here if we accepted the final state hadron system
555  // so now we need to write everything to ghep
556 
557  // First the initial state nucleon.
558  initial_nucleon->SetMomentum(p4initial_nucleon);
559 
560  // and the remnant nucleus
561  if(have_nucleus){
562  double Mi = PDGLibrary::Instance()->Find(target_nucleus->Pdg() )-> Mass();
563  remnant_nucleus->SetMomentum(pxb,pyb,pzb, Mi - p4initial_nucleon.E() + removalenergy);
564  }
565 
566  // Getting this v4 is a position, i.e. a position within the nucleus (?)
567  // possibly it takes on a non-trivial value only for local Fermi gas
568  // or for sophisticated treatments of intranuclear rescattering.
569  TLorentzVector v4(*neutrino->X4());
570 
571  // Now add the final nucleon
572 
573  interaction->KinePtr()->SetHadSystP4(p4final_nucleon);
574 
575  GHepStatus_t ist = (tgt.IsNucleus()) ? kIStHadronInTheNucleus : kIStStableFinalState;
576  event->AddParticle(interaction->RecoilNucleonPdg(), ist, event->HitNucleonPosition(),-1,-1,-1, interaction->KinePtr()->HadSystP4(), v4);
577 
578  LOG("QELEvent",pDEBUG) << "Generate Nucleon - End";
579 
580 }
double RemovalEnergy(void) const
Definition: NuclearModelI.h:65
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
int RecoilNucleonPdg(void) const
recoil nucleon pdg
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:66
double Mass(Resonance_t res)
resonance mass (GeV)
void SetMomentum(const TLorentzVector &p4)
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:75
int Pdg(void) const
Definition: GHepParticle.h:63
Summary information for an interaction.
Definition: Interaction.h:56
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
const int kPdgTgtFreeP
Definition: PDGCodes.h:199
#define pWARN
Definition: Messenger.h:60
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:57
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
void SetHadSystP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:307
virtual bool GenerateNucleon(const Target &) const =0
double NonNegative(double x)
Definition: MathUtils.cxx:273
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:75
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
const NuclearModelI * fNuclModel
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
enum genie::EGHepStatus GHepStatus_t
#define pDEBUG
Definition: Messenger.h:63
void QELEventGeneratorSuSA::LoadConfig ( void  )
private

Definition at line 594 of file QELEventGeneratorSuSA.cxx.

595 {
596  fNuclModel = 0;
597  RgKey nuclkey = "NuclearModel";
598  fNuclModel = dynamic_cast<const NuclearModelI *> ( this->SubAlg(nuclkey) );
599  assert( fNuclModel );
600 
601  // Sub-algorithm for generating events on free nucleon targets
602  // (not handled by the SuSAv2 calculation)
603  fFreeNucleonEventGenerator = dynamic_cast<const EventRecordVisitorI* >(
604  this->SubAlg("FreeNucleonEventGenerator") );
605  assert( fFreeNucleonEventGenerator );
606 
607  //-- Maximum q3 in input hadron tensors
608  GetParam( "QEL-Q3Max", fQ3Max ) ;
609 
610  //-- Whether to force nucleons to be bound
611  GetParam( "QEL-ForceBound", fForceBound) ;
612 
613  //-- Whether to force Eb to come from the nuclear model
614  GetParam( "QEL-ForceEbFromModel", fForceEbFromModel) ;
615 
616  //-- Whether to force some fixed Eb
617  GetParam( "QEL-ForceFixedEb", fForceFixEb) ;
618  GetParam( "QEL-EbOR", fEbOR) ;
619 
620  //-- Carbon Eb - needed for scaling
621  this->GetParam( "RFG-NucRemovalE@Pdg=1000060120", fEbC);
622 
623  //-- Safety factor for the maximum differential cross section
624  GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor , 5.00 ) ;
625 
626  //-- Minimum energy for which max xsec would be cached, forcing explicit
627  // calculation for lower energies
628  //-- I've set this to an extremely high value to avoid problems with the
629  // SuSAv2 model seen during testing. Everything seems to work OK when
630  // the cache is disabled. - S. Gardiner, 1 July 2020
631  GetParamDef( "Cache-MinEnergy", fEMin, 1000.00 ) ;
632 
633  //-- Maximum allowed fractional cross section deviation from maxim cross
634  // section used in rejection method
635  GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ;
636  assert( fMaxXSecDiffTolerance >= 0. );
637 
638  //-- Generate kinematics uniformly over allowed phase space and compute
639  // an event weight? NOT IMPLEMENTED FOR SUSA YET!
640  GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false ) ;
641 }
bool fGenerateUniformly
uniform over allowed phase space + event weight?
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
double fSafetyFactor
maxxsec -> maxxsec * safety_factor
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec>maxxsec
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:46
const EventRecordVisitorI * fFreeNucleonEventGenerator
string RgKey
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const NuclearModelI * fNuclModel
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
void QELEventGeneratorSuSA::ProcessEventRecord ( GHepRecord event) const
virtual

Implements genie::EventRecordVisitorI.

Definition at line 68 of file QELEventGeneratorSuSA.cxx.

69 {
70  // If we're working with a free nucleon target, then the SuSAv2 calculation
71  // isn't set up to handle it. In these cases, we delegate the work to another
72  // EventRecordVisitorI object (likely QELEventGenerator) configured to run as
73  // a sub-algorithm.
74  if ( !event->Summary()->InitState().Tgt().IsNucleus() ) {
76  }
77 
78  //-- Access cross section algorithm for running thread
80  const EventGeneratorI * evg = rtinfo->RunningThread();
81  fXSecModel = evg->CrossSectionAlg();
82  //event->Print();
83  this -> SelectLeptonKinematics(event);
84  //event->Print();
85  this -> AddTargetNucleusRemnant(event);
86  //event->Print();
87  this -> GenerateNucleon(event);
88  //event->Print();
89 }
virtual void ProcessEventRecord(GHepRecord *event_rec) const =0
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:91
bool IsNucleus(void) const
Definition: Target.cxx:272
Defines the EventGeneratorI interface.
void SelectLeptonKinematics(GHepRecord *event) const
void AddTargetNucleusRemnant(GHepRecord *event) const
void GenerateNucleon(GHepRecord *event) const
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
static RunningThreadInfo * Instance(void)
const EventRecordVisitorI * fFreeNucleonEventGenerator
const InitialState & InitState(void) const
Definition: Interaction.h:69
const Target & Tgt(void) const
Definition: InitialState.h:66
const EventGeneratorI * RunningThread(void)
Keep info on the event generation thread currently on charge. This is used so that event generation m...
void QELEventGeneratorSuSA::SelectLeptonKinematics ( GHepRecord event) const
private

Definition at line 91 of file QELEventGeneratorSuSA.cxx.

92 {
93 
94  // -- Event Properties -----------------------------//
95  Interaction * interaction = event->Summary();
96  Kinematics * kinematics = interaction->KinePtr();
97 
98  // Choose the appropriate minimum Q^2 value based on the interaction
99  // mode (this is important for EM interactions since the differential
100  // cross section blows up as Q^2 --> 0)
101  double Q2min = genie::controls::kMinQ2Limit; // CC/NC limit
102  if ( interaction->ProcInfo().IsEM() ) Q2min = genie::utils::kinematics
103  ::electromagnetic::kMinQ2Limit; // EM limit
104 
105  // The SuSA 1p1h model kinematics works in a system where
106  // the whole nuclear target system has no momentum.
107  double Enu = interaction->InitState().ProbeE(kRfLab);
108 
109  int NuPDG = interaction->InitState().ProbePdg();
110  int TgtPDG = interaction->InitState().TgtPdg();
111  // interacton vtx
112  TLorentzVector v4(*event->Probe()->X4());
113  TLorentzVector tempp4(0.,0.,0.,0.);
114 
115  GHepParticle * nucleus = event->TargetNucleus();
116  bool have_nucleus = nucleus != 0;
117 
118  // -- Lepton Kinematic Limits ----------------------------------------- //
119  double Costh = 0.0; // lepton angle
120  double CosthMax = 1.0;
121  double CosthMin = -1.0;
122 
123  double T = 0.0; // lepton kinetic energy
124  double TMax = std::numeric_limits<double>::max();
125  double TMin = 0.0;
126 
127  double Plep = 0.0; // lepton 3 momentum
128  double Elep = 0.0; // lepton energy
129  double LepMass = interaction->FSPrimLepton()->Mass();
130 
131  double Q0 = 0.0; // energy component of q four vector
132  double Q3 = 0.0; // magnitude of transfered 3 momentum
133  double Q2 = 0.0; // properly Q^2 (Q squared) - transfered 4 momentum.
134 
135  // Set lepton KE TMax for for throwing rndm in the accept/reject loop.
136  // We can accidentally set it too high, because the xsec will return zero.
137  // This way if someone reuses this code, they are not tripped up by it.
138  TMax = Enu - LepMass;
139 
140  // Set Tmin for throwing rndm in the accept/reject loop
141  // the hadron tensors we expect will be limited in q3
142  // therefore also the outgoing lepton KE can't be too low or costheta too backward
143  // make the accept/reject loop more efficient by using Min values.
144  if(Enu < fQ3Max){
145  TMin = 0 ;
146  CosthMin = -1 ;
147  } else {
148  TMin = TMath::Sqrt(TMath::Power(LepMass, 2) + TMath::Power((Enu - fQ3Max), 2)) - LepMass;
149  CosthMin = TMath::Sqrt(1 - TMath::Power((fQ3Max / Enu ), 2));
150  }
151 
152 
153  // -- Generate and Test the Kinematics----------------------------------//
154 
155  RandomGen * rnd = RandomGen::Instance();
156  bool accept = false;
157  unsigned int iter = 0;
158  unsigned int maxIter = kRjMaxIterations * 1000;
159 
160  //e-scat xsecs blow up close to theta=0, MC methods won't work so well...
161  // NOTE: SuSAv2 1p1h e-scatting has not been validated yet, use with caution
162  if ( NuPDG == 11 ) maxIter *= 1000;
163 
164  // Get Max XSec:
165  double XSecMax = this->MaxXSec( event );
166 
167  LOG("Kinematics", pDEBUG) << "Max XSec = " << XSecMax;
168 
169  // loop over different (randomly) selected T and Costh
170  while (!accept) {
171  iter++;
172  if(iter > maxIter) {
173  // error if try too many times
174  LOG("QELEvent", pERROR)
175  << "Couldn't select a valid Tmu, CosTheta pair after "
176  << iter << " iterations";
177  event->EventFlags()->SetBitNumber(kKineGenErr, true);
179  exception.SetReason("Couldn't select lepton kinematics");
180  exception.SwitchOnFastForward();
181  throw exception;
182  }
183 
184  // generate random kinetic energy T and Costh
185  T = TMin + (TMax-TMin)*rnd->RndKine().Rndm();
186  Costh = CosthMin + (CosthMax-CosthMin)*rnd->RndKine().Rndm();
187 
188  // Calculate useful values for judging this choice
189  Plep = TMath::Sqrt( T * (T + (2.0 * LepMass))); // ok is sqrt(E2 - m2)
190  Q3 = TMath::Sqrt(Plep*Plep + Enu*Enu - 2.0 * Plep * Enu * Costh);
191 
192  // Calculate what we need to enforce the minimum Q^2 cut
193  Q0 = Enu - (T + LepMass);
194  Q2 = Q3*Q3 - Q0*Q0;
195 
196  // Anti neutrino elastic scattering case - include delta in xsec
197  if (!have_nucleus){
198  genie::utils::mec::Getq0q3FromTlCostl(T, Costh, Enu, LepMass, Q0, Q3);
199  Q3 = sqrt(Q0*Q0+2*kNucleonMass*Q0);
200  genie::utils::mec::GetTlCostlFromq0q3(Q0, Q3, Enu, LepMass, T, Costh);
201  LOG("QELEvent", pDEBUG) << " Anu elastic case. T, Costh: " << T << ", " << Costh ;
202  LOG("QELEvent", pDEBUG) << " Anu elastic case. Q0, Q3: " << Q0 << ", " << Q3 ;
203  }
204 
205  // Don't bother doing hard work if the selected Q3 is greater than Q3Max
206  // or if Q^2 is less than the minimum allowed value
207  if ( Q3 < fQ3Max && Q2 >= Q2min ){
208 
209  kinematics->SetKV(kKVTl, T);
210  kinematics->SetKV(kKVctl, Costh);
211  LOG("QELEvent", pDEBUG) << " T, Costh, Q2: " << T << ", " << Costh << ", " << Q2;
212 
213  double XSec = fXSecModel->XSec(interaction, kPSTlctl);
214 
215  // Some debugging if things go wrong here...
216  if (XSec > XSecMax) {
217  LOG("QELEvent", pDEBUG) << "XSec in cm2 is " << XSec/(units::cm2);
218  LOG("QELEvent", pDEBUG) << "XSec in cm2 /neutron is " << XSec/(units::cm2*pdg::IonPdgCodeToZ(TgtPDG));
219  LOG("QELEvent", pDEBUG) << "XSecMax in cm2 /neutron is " << XSecMax/(units::cm2*pdg::IonPdgCodeToZ(TgtPDG));
220  LOG("QELEvent", pERROR) << " T, Costh, Q2: " << T << ", " << Costh << ", " << Q2;
221  LOG("QELEvent", pERROR) << "XSec is > XSecMax for nucleus " << TgtPDG << " "
222  << XSec << " > " << XSecMax
223  << " don't let this happen.";
224  }
225  // decide whether to accept or reject these kinematics
226  this->AssertXSecLimits( interaction, XSec, XSecMax );
227  accept = XSec > XSecMax*rnd->RndKine().Rndm();
228  LOG("QELEvent", pINFO) << "Xsec, Max, Accept: " << XSec << ", "
229  << XSecMax << ", " << accept;
230  LOG("QELEvent", pDEBUG) << "XSec in cm2 /neutron is " << XSec/(units::cm2*pdg::IonPdgCodeToZ(TgtPDG));
231  LOG("QELEvent", pDEBUG) << "XSecMax in cm2 /neutron is " << XSecMax/(units::cm2*pdg::IonPdgCodeToZ(TgtPDG));
232 
233 
234  }// end if passes q3 test
235  } // end main accept-reject loop
236 
237  // -- finish lepton kinematics
238  // If the code got here, then we accepted some kinematics
239  // and we can proceed to generate the final state.
240 
241  // define coordinate system wrt neutrino: z along neutrino, xy perp
242 
243  // Cos theta gives us z, the rest in xy:
244  double PlepZ = Plep * Costh;
245  double PlepXY = Plep * TMath::Sqrt(1. - TMath::Power(Costh,2));
246 
247  // random rotation about unit vector for phi direction
248  double phi= 2 * kPi * rnd->RndLep().Rndm();
249  // now fill x and y from PlepXY
250  double PlepX = PlepXY * TMath::Cos(phi);
251  double PlepY = PlepXY * TMath::Sin(phi);
252 
253  // Rotate lepton momentum vector from the reference frame (x'y'z') where
254  // {z':(neutrino direction), z'x':(theta plane)} to the LAB
255  TVector3 unit_nudir = event->Probe()->P4()->Vect().Unit();
256  TVector3 p3l(PlepX, PlepY, PlepZ);
257  p3l.RotateUz(unit_nudir);
258 
259  // Lepton 4-momentum in LAB
260  Elep = TMath::Sqrt(LepMass*LepMass + PlepX*PlepX + PlepY*PlepY + PlepZ*PlepZ);
261  TLorentzVector p4l(p3l,Elep);
262 
263  // Figure out the final-state primary lepton PDG code
264  int pdgc = interaction->FSPrimLepton()->PdgCode();
265  int momidx = event->ProbePosition();
266 
267  // -- Store Values ------------------------------------------//
268  // -- Interaction: Q2
269  Q0 = Enu - Elep;
270  Q2 = Q3*Q3 - Q0*Q0;
271  double gy = Q0 / Enu;
272  double gx = kinematics::Q2YtoX(Enu, 2 * kNucleonMass, Q2, gy);
273  double gW = kinematics::XYtoW(Enu, 2 * kNucleonMass, gx, gy);
274 
275  interaction->KinePtr()->SetQ2(Q2, true);
276  interaction->KinePtr()->Sety(gy, true);
277  interaction->KinePtr()->Setx(gx, true);
278  interaction->KinePtr()->SetW(gW, true);
279  interaction->KinePtr()->SetFSLeptonP4(p4l);
280 
281  // -- Lepton
282  event->AddParticle( pdgc, kIStStableFinalState, momidx, -1, -1, -1, p4l, v4);
283 
284  LOG("QELEvent",pDEBUG) << "~~~ LEPTON DONE ~~~";
285 }
virtual double MaxXSec(GHepRecord *evrec) const
bool GetTlCostlFromq0q3(double q0, double q3, double Enu, double ml, double &Tl, double &costl)
Definition: MECUtils.cxx:82
#define pERROR
Definition: Messenger.h:59
TRandom3 & RndLep(void) const
rnd number generator used by final state primary lepton generators
Definition: RandomGen.h:62
static const double kNucleonMass
Definition: Constants.h:77
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
virtual void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
double XYtoW(double Ev, double M, double x, double y)
Definition: KineUtils.cxx:1179
static const double kMinQ2Limit
Definition: Controls.h:41
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:277
Summary information for an interaction.
Definition: Interaction.h:56
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void SetFSLeptonP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:297
static constexpr double cm2
Definition: Units.h:69
int ProbePdg(void) const
Definition: InitialState.h:64
bool Getq0q3FromTlCostl(double Tl, double costl, double Enu, double ml, double &q0, double &q3)
Definition: MECUtils.cxx:121
#define pINFO
Definition: Messenger.h:62
static int max(int a, int b)
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:50
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool IsEM(void) const
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
void SetKV(KineVar_t kv, double value)
Definition: Kinematics.cxx:335
Kinematical utilities.
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
double Q2YtoX(double Ev, double M, double Q2, double y)
Definition: KineUtils.cxx:1209
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
int TgtPdg(void) const
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:52
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
double ProbeE(RefFrame_t rf) const
static const double kPi
Definition: Constants.h:37
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
#define pDEBUG
Definition: Messenger.h:63

Member Data Documentation

double genie::QELEventGeneratorSuSA::fEbC
private

Definition at line 64 of file QELEventGeneratorSuSA.h.

double genie::QELEventGeneratorSuSA::fEbOR
private

Definition at line 61 of file QELEventGeneratorSuSA.h.

bool genie::QELEventGeneratorSuSA::fForceBound
private

Definition at line 58 of file QELEventGeneratorSuSA.h.

bool genie::QELEventGeneratorSuSA::fForceEbFromModel
private

Definition at line 59 of file QELEventGeneratorSuSA.h.

bool genie::QELEventGeneratorSuSA::fForceFixEb
private

Definition at line 60 of file QELEventGeneratorSuSA.h.

const EventRecordVisitorI* genie::QELEventGeneratorSuSA::fFreeNucleonEventGenerator
private

Delegate event generation for free nucleon targets (which are not handled by the SuSAv2 calculation) to a different event generator

Definition at line 69 of file QELEventGeneratorSuSA.h.

const NuclearModelI* genie::QELEventGeneratorSuSA::fNuclModel
private

Definition at line 55 of file QELEventGeneratorSuSA.h.

double genie::QELEventGeneratorSuSA::fQ3Max
private

Definition at line 57 of file QELEventGeneratorSuSA.h.


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