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

Generates values for the kinematic variables describing QEL neutrino interaction events for Smith-Moniz model. Is a concrete implementation of the EventRecordVisitorI interface. More...

#include <QELEventGeneratorSM.h>

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

Public Member Functions

 QELEventGeneratorSM ()
 
 QELEventGeneratorSM (string config)
 
 ~QELEventGeneratorSM ()
 
void ProcessEventRecord (GHepRecord *event_rec) 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)
 
double ComputeMaxXSec (const Interaction *in) const
 
void AddTargetNucleusRemnant (GHepRecord *evrec) const
 add a recoiled nucleus remnant More...
 
double ComputeMaxXSec2 (const Interaction *in) const
 
double MaxXSec2 (GHepRecord *evrec) const
 
double FindMaxXSec2 (const Interaction *in) const
 
void CacheMaxXSec2 (const Interaction *in, double xsec) const
 
CacheBranchFxAccessCacheBranch2 (const Interaction *in) const
 
double ComputeMaxDiffv (const Interaction *in) const
 
double MaxDiffv (GHepRecord *evrec) const
 
double FindMaxDiffv (const Interaction *in) const
 
void CacheMaxDiffv (const Interaction *in, double xsec) const
 
CacheBranchFxAccessCacheBranchDiffv (const Interaction *in) const
 

Private Attributes

SmithMonizUtilssm_utils
 
KinePhaseSpace_t fkps
 
bool fGenerateNucleonInNucleus
 generate struck nucleon in nucleus More...
 
double fQ2Min
 Q2-threshold for seeking the second maximum. More...
 
double fSafetyFacor_nu
 

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

Generates values for the kinematic variables describing QEL neutrino interaction events for Smith-Moniz model. Is a concrete implementation of the EventRecordVisitorI interface.

[1] R.A.Smith and E.J.Moniz, Nuclear Physics B43, (1972) 605-622
[2] K.S. Kuzmin, V.V. Lyubushkin, V.A.Naumov Eur. Phys. J. C54, (2008) 517-538

Author
Igor Kakorin kakor.nosp@m.in@j.nosp@m.inr.r.nosp@m.u
Joint Institute for Nuclear Research
adapted from fortran code provided by:

Konstantin Kuzmin kkuzm.nosp@m.in@t.nosp@m.heor..nosp@m.jinr.nosp@m..ru,
Joint Institute for Nuclear Research, Institute for Theoretical and Experimental Physics
Vladimir Lyubushkin,
Joint Institute for Nuclear Research
Vadim Naumov vnaum.nosp@m.ov@t.nosp@m.heor..nosp@m.jinr.nosp@m..ru,
Joint Institute for Nuclear Research
based on code of: Costas Andreopoulos <constantinos.andreopoulos cern.ch> University of Liverpool & STFC Rutherford Appleton Laboratory

May 05, 2017

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

Definition at line 48 of file QELEventGeneratorSM.h.

Constructor & Destructor Documentation

QELEventGeneratorSM::QELEventGeneratorSM ( )

Definition at line 63 of file QELEventGeneratorSM.cxx.

63  :
64 KineGeneratorWithCache("genie::QELEventGeneratorSM")
65 {
66 
67 }
QELEventGeneratorSM::QELEventGeneratorSM ( string  config)

Definition at line 69 of file QELEventGeneratorSM.cxx.

69  :
70 KineGeneratorWithCache("genie::QELEventGeneratorSM", config)
71 {
72 
73 }
static Config * config
Definition: config.cpp:1054
QELEventGeneratorSM::~QELEventGeneratorSM ( )

Definition at line 75 of file QELEventGeneratorSM.cxx.

76 {
77 
78 }

Member Function Documentation

CacheBranchFx * QELEventGeneratorSM::AccessCacheBranch2 ( const Interaction in) const
private

Definition at line 641 of file QELEventGeneratorSM.cxx.

643 {
644 // Returns the cache branch for this algorithm and this interaction. If no
645 // branch is found then one is created.
646 
647  Cache * cache = Cache::Instance();
648 
649  // build the cache branch key as: namespace::algorithm/config/interaction
650  string algkey = this->Id().Key();
651  string intkey = interaction->AsString();
652  string key = cache->CacheBranchKey(algkey, intkey, "2nd");
653 
654  CacheBranchFx * cache_branch =
655  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
656  if(!cache_branch) {
657  //-- create the cache branch at the first pass
658  LOG("Kinematics", pINFO) << "No Max d^nXSec/d{K}^n cache branch found";
659  LOG("Kinematics", pINFO) << "Creating cache branch - key = " << key;
660 
661  cache_branch = new CacheBranchFx("max[d^nXSec/d^n{K}] over phase space");
662  cache->AddCacheBranch(key, cache_branch);
663  }
664  assert(cache_branch);
665 
666  return cache_branch;
667 }
void AddCacheBranch(string key, CacheBranchI *branch)
Definition: Cache.cxx:88
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
def key(type, name=None)
Definition: graph.py:13
string CacheBranchKey(string k0, string k1="", string k2="") const
Definition: Cache.cxx:93
#define pINFO
Definition: Messenger.h:62
CacheBranchI * FindCacheBranch(string key)
finding/adding cache branches
Definition: Cache.cxx:80
GENIE Cache Memory.
Definition: Cache.h:38
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:97
static Cache * Instance(void)
Definition: Cache.cxx:67
A simple cache branch storing the cached data in a TNtuple.
Definition: CacheBranchFx.h:37
string Key(void) const
Definition: AlgId.h:46
CacheBranchFx * QELEventGeneratorSM::AccessCacheBranchDiffv ( const Interaction in) const
private

Definition at line 797 of file QELEventGeneratorSM.cxx.

798 {
799 // Returns the cache branch for this algorithm and this interaction. If no
800 // branch is found then one is created.
801 
802  Cache * cache = Cache::Instance();
803 
804  // build the cache branch key as: namespace::algorithm/config/interaction
805  string algkey = this->Id().Key();
806  string intkey = interaction->AsString();
807  string key = cache->CacheBranchKey(algkey, intkey, "diffv");
808 
809  CacheBranchFx * cache_branch =
810  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
811  if(!cache_branch)
812  {
813  //-- create the cache branch at the first pass
814  LOG("Kinematics", pINFO) << "No Max vmax(Q2)-vmin(Q2) cache branch found";
815  LOG("Kinematics", pINFO) << "Creating cache branch - key = " << key;
816 
817  cache_branch = new CacheBranchFx("max[vmax(Q2)-vmin(Q2)] over phase space");
818  cache->AddCacheBranch(key, cache_branch);
819  }
820  assert(cache_branch);
821 
822  return cache_branch;
823 }
void AddCacheBranch(string key, CacheBranchI *branch)
Definition: Cache.cxx:88
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
def key(type, name=None)
Definition: graph.py:13
string CacheBranchKey(string k0, string k1="", string k2="") const
Definition: Cache.cxx:93
#define pINFO
Definition: Messenger.h:62
CacheBranchI * FindCacheBranch(string key)
finding/adding cache branches
Definition: Cache.cxx:80
GENIE Cache Memory.
Definition: Cache.h:38
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:97
static Cache * Instance(void)
Definition: Cache.cxx:67
A simple cache branch storing the cached data in a TNtuple.
Definition: CacheBranchFx.h:37
string Key(void) const
Definition: AlgId.h:46
void QELEventGeneratorSM::AddTargetNucleusRemnant ( GHepRecord evrec) const
private

add a recoiled nucleus remnant

Definition at line 337 of file QELEventGeneratorSM.cxx.

338 {
339 // add the remnant nuclear target at the GHEP record
340 
341  LOG("QELEvent", pINFO) << "Adding final state nucleus";
342 
343  double Px = 0;
344  double Py = 0;
345  double Pz = 0;
346  double E = 0;
347 
348  GHepParticle * nucleus = evrec->TargetNucleus();
349  int A = nucleus->A();
350  int Z = nucleus->Z();
351 
352  int fd = nucleus->FirstDaughter();
353  int ld = nucleus->LastDaughter();
354 
355  for(int id = fd; id <= ld; id++)
356  {
357 
358  // compute A,Z for final state nucleus & get its PDG code and its mass
359  GHepParticle * particle = evrec->Particle(id);
360  assert(particle);
361  int pdgc = particle->Pdg();
362  bool is_p = pdg::IsProton (pdgc);
363  bool is_n = pdg::IsNeutron(pdgc);
364 
365  if (is_p) Z--;
366  if (is_p || is_n) A--;
367 
368  Px += particle->Px();
369  Py += particle->Py();
370  Pz += particle->Pz();
371  E += particle->E();
372 
373  }//daughters
374 
375  TParticlePDG * remn = 0;
376  int ipdgc = pdg::IonPdgCode(A, Z);
377  remn = PDGLibrary::Instance()->Find(ipdgc);
378  if(!remn)
379  {
380  LOG("HadronicVtx", pFATAL)
381  << "No particle with [A = " << A << ", Z = " << Z
382  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
383  assert(remn);
384  }
385 
386  double Mi = nucleus->Mass();
387  Px *= -1;
388  Py *= -1;
389  Pz *= -1;
390  E = Mi-E;
391 
392  // Add the nucleus to the event record
393  LOG("QELEvent", pINFO)
394  << "Adding nucleus [A = " << A << ", Z = " << Z
395  << ", pdgc = " << ipdgc << "]";
396 
397  int imom = evrec->TargetNucleusPosition();
398  evrec->AddParticle(
399  ipdgc,kIStStableFinalState, imom,-1,-1,-1, Px,Py,Pz,E, 0,0,0,0);
400 
401  LOG("QELEvent", pINFO) << "Done";
402  LOG("QELEvent", pINFO) << *evrec;
403 }
int Z(void) const
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:104
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
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:286
#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
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:491
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
virtual int TargetNucleusPosition(void) const
Definition: GHepRecord.cxx:362
double Py(void) const
Get Py.
Definition: GHepParticle.h:89
void QELEventGeneratorSM::CacheMaxDiffv ( const Interaction in,
double  xsec 
) const
private

Definition at line 774 of file QELEventGeneratorSM.cxx.

775 {
776  LOG("Kinematics", pINFO)
777  << "Adding the computed max{vmax(Q2)-vmin(Q2)} value to cache";
779 
780  double E = this->Energy(interaction);
781  if(max_diffv>0) cb->AddValues(E,max_diffv);
782 
783  if(! cb->Spl() )
784  {
785  if( cb->Map().size() > 40 ) cb->CreateSpline();
786  }
787 
788  if( cb->Spl() )
789  {
790  if( E < cb->Spl()->XMin() || E > cb->Spl()->XMax() )
791  {
792  cb->CreateSpline();
793  }
794  }
795 }
CacheBranchFx * AccessCacheBranchDiffv(const Interaction *in) const
void AddValues(double x, double y)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
Spline * Spl(void) const
Definition: CacheBranchFx.h:47
#define pINFO
Definition: Messenger.h:62
double XMax(void) const
Definition: Spline.h:77
const map< double, double > & Map(void) const
Definition: CacheBranchFx.h:46
virtual double Energy(const Interaction *in) const
A simple cache branch storing the cached data in a TNtuple.
Definition: CacheBranchFx.h:37
void QELEventGeneratorSM::CacheMaxXSec2 ( const Interaction in,
double  xsec 
) const
private

Definition at line 620 of file QELEventGeneratorSM.cxx.

622 {
623  LOG("Kinematics", pINFO)
624  << "Adding the computed max{dxsec/dK} value to cache";
626 
627  double E = this->Energy(interaction);
628  if(max_xsec>0) cb->AddValues(E,max_xsec);
629 
630  if(! cb->Spl() ) {
631  if( cb->Map().size() > 40 ) cb->CreateSpline();
632  }
633 
634  if( cb->Spl() ) {
635  if( E < cb->Spl()->XMin() || E > cb->Spl()->XMax() ) {
636  cb->CreateSpline();
637  }
638  }
639 }
void AddValues(double x, double y)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
Spline * Spl(void) const
Definition: CacheBranchFx.h:47
#define pINFO
Definition: Messenger.h:62
double XMax(void) const
Definition: Spline.h:77
const map< double, double > & Map(void) const
Definition: CacheBranchFx.h:46
virtual double Energy(const Interaction *in) const
A simple cache branch storing the cached data in a TNtuple.
Definition: CacheBranchFx.h:37
CacheBranchFx * AccessCacheBranch2(const Interaction *in) const
double QELEventGeneratorSM::ComputeMaxDiffv ( const Interaction in) const
private

Definition at line 669 of file QELEventGeneratorSM.cxx.

670 {
671  double max_diffv = -1;
672  const int N_Q2 = 10;
673 
675  for (int Q2_n = 0; Q2_n<N_Q2; Q2_n++) // Scan around Q2
676  {
677  double Q2 = rQ2.min + 1.*Q2_n * (rQ2.max-rQ2.min)/N_Q2;
678  Range1D_t rv = sm_utils->vQES_SM_lim(Q2);
679  if (rv.max-rv.min > max_diffv)
680  max_diffv = rv.max-rv.min;
681  } // Done with Q2 scan
682  max_diffv *= fSafetyFactor;
683  return max_diffv;
684 
685 }
Range1D_t Q2QES_SM_lim(void) const
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
double fSafetyFactor
maxxsec -> maxxsec * safety_factor
A simple [min,max] interval for doubles.
Definition: Range1.h:42
Range1D_t vQES_SM_lim(double Q2) const
double max
Definition: Range1.h:53
double min
Definition: Range1.h:52
double QELEventGeneratorSM::ComputeMaxXSec ( const Interaction in) const
privatevirtual

Implements genie::KineGeneratorWithCache.

Definition at line 451 of file QELEventGeneratorSM.cxx.

452 {
453  double xsec_max = -1;
454  const int N_Q2 = 8;
455  const int N_v = 8;
456 
458  const double logQ2min = TMath::Log(TMath::Max(rQ2.min, eps));
459  const double logQ2max = TMath::Log(TMath::Min(rQ2.max, fQ2Min));
460 
461  double tmp_xsec_max = -1;
462  // Now scan through kinematical variables Q2,v
463  for (int Q2_n=0; Q2_n < N_Q2; Q2_n++)
464  {
465  // Scan around Q2
466  double Q2 = TMath::Exp(Q2_n * (logQ2max-logQ2min)/N_Q2 + logQ2min);
467  Range1D_t rv = sm_utils->vQES_SM_lim(Q2);
468  const double logvmin = TMath::Log(TMath::Max(rv.min, eps));
469  const double logvmax = TMath::Log(TMath::Max(rv.max, TMath::Max(rv.min, eps)));
470  for (int v_n=0; v_n < N_v; v_n++)
471  {
472  // Scan around v
473  double v = TMath::Exp(v_n * (logvmax-logvmin)/N_v + logvmin);
474  Kinematics * kinematics = interaction->KinePtr();
475  kinematics->SetKV(kKVQ2, Q2);
476  kinematics->SetKV(kKVv, v);
477  // Compute the QE cross section for the current kinematics
478  double xs = fXSecModel->XSec(interaction, fkps);
479  if (xs > tmp_xsec_max)
480  tmp_xsec_max = xs;
481  } // Done with v scan
482  }// Done with Q2 scan
483 
484  xsec_max = tmp_xsec_max;
485  // Apply safety factor, since value retrieved from the cache might
486  // correspond to a slightly different value
487  xsec_max *= fSafetyFactor;
488  return xsec_max;
489 
490 }
Range1D_t Q2QES_SM_lim(void) const
double fQ2Min
Q2-threshold for seeking the second maximum.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
double fSafetyFactor
maxxsec -> maxxsec * safety_factor
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
A simple [min,max] interval for doubles.
Definition: Range1.h:42
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
Range1D_t vQES_SM_lim(double Q2) const
void SetKV(KineVar_t kv, double value)
Definition: Kinematics.cxx:335
double max
Definition: Range1.h:53
double min
Definition: Range1.h:52
double QELEventGeneratorSM::ComputeMaxXSec2 ( const Interaction in) const
private

Definition at line 492 of file QELEventGeneratorSM.cxx.

493 {
494  double xsec_max = -1;
495  const int N_Q2 = 8;
496  const int N_v = 8;
497 
499  if (rQ2.max<fQ2Min) return xsec_max;
500  const double logQ2min = TMath::Log(fQ2Min);
501  const double logQ2max = TMath::Log(rQ2.max);
502 
503  double tmp_xsec_max = -1;
504  // Now scan through kinematical variables Q2,v
505  for (int Q2_n=0; Q2_n < N_Q2; Q2_n++)
506  {
507  // Scan around Q2
508  double Q2 = TMath::Exp(Q2_n * (logQ2max-logQ2min)/N_Q2 + logQ2min);
509  Range1D_t rv = sm_utils->vQES_SM_lim(Q2);
510  const double logvmin = TMath::Log(TMath::Max(rv.min, eps));
511  const double logvmax = TMath::Log(TMath::Max(rv.max, TMath::Max(rv.min, eps)));
512  for (int v_n=0; v_n < N_v; v_n++)
513  {
514  // Scan around v
515  double v = TMath::Exp(v_n * (logvmax-logvmin)/N_v + logvmin);
516  Kinematics * kinematics = interaction->KinePtr();
517  kinematics->SetKV(kKVQ2, Q2);
518  kinematics->SetKV(kKVv, v);
519  // Compute the QE cross section for the current kinematics
520  double xs = fXSecModel->XSec(interaction, fkps);
521  if (xs > tmp_xsec_max)
522  tmp_xsec_max = xs;
523  } // Done with v scan
524  }// Done with Q2 scan
525 
526  xsec_max = tmp_xsec_max;
527  // Apply safety factor, since value retrieved from the cache might
528  // correspond to a slightly different value
529  xsec_max *= fSafetyFactor;
530  return xsec_max;
531 
532 }
Range1D_t Q2QES_SM_lim(void) const
double fQ2Min
Q2-threshold for seeking the second maximum.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
double fSafetyFactor
maxxsec -> maxxsec * safety_factor
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
A simple [min,max] interval for doubles.
Definition: Range1.h:42
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
Range1D_t vQES_SM_lim(double Q2) const
void SetKV(KineVar_t kv, double value)
Definition: Kinematics.cxx:335
double max
Definition: Range1.h:53
double min
Definition: Range1.h:52
void QELEventGeneratorSM::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 405 of file QELEventGeneratorSM.cxx.

406 {
407  Algorithm::Configure(config);
408  this->LoadConfig();
409 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void QELEventGeneratorSM::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 411 of file QELEventGeneratorSM.cxx.

412 {
414 
415  Registry r( "QELEventGeneratorSM_specific", false ) ;
416  r.Set( "sm_utils_algo", RgAlg("genie::SmithMonizUtils","Default") ) ;
417 
419 
420  this->LoadConfig();
421 }
static Config * config
Definition: config.cpp:1054
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
double QELEventGeneratorSM::FindMaxDiffv ( const Interaction in) const
private

Definition at line 727 of file QELEventGeneratorSM.cxx.

728 {
729 // Find a cached maximum of vmax(Q2)-vmin(Q2) for xsec algorithm & interaction and
730 // close to the specified energy
731 
732  // get neutrino energy
733  double E = this->Energy(interaction);
734  LOG("Kinematics", pINFO) << "E = " << E;
735 
736  if(E < fEMin) {
737  LOG("Kinematics", pINFO)
738  << "Below minimum energy - Forcing explicit calculation";
739  return -1.;
740  }
741 
742  // access the the cache branch
744 
745  // if there are enough points stored in the cache buffer to build a
746  // spline, then intepolate
747  if( cb->Spl() ) {
748  if( E >= cb->Spl()->XMin() && E <= cb->Spl()->XMax())
749  {
750  double spl_maxdiffv = cb->Spl()->Evaluate(E);
751  LOG("Kinematics", pINFO)
752  << "\nInterpolated: max vmax(Q2)-vmin(Q2) (E=" << E << ") = " << spl_maxdiffv;
753  return spl_maxdiffv;
754  }
755  LOG("Kinematics", pINFO)
756  << "Outside spline boundaries - Forcing explicit calculation";
757  return -1.;
758  }
759 
760  // if there are not enough points at the cache buffer to have a spline,
761  // look whether there is another point that is sufficiently close
762  double dE = TMath::Min(0.25, 0.05*E);
763  const map<double,double> & fmap = cb->Map();
764  map<double,double>::const_iterator iter = fmap.lower_bound(E);
765  if(iter != fmap.end())
766  {
767  if(TMath::Abs(E - iter->first) < dE) return iter->second;
768  }
769 
770  return -1;
771 
772 }
intermediate_table::const_iterator const_iterator
double Evaluate(double x) const
Definition: Spline.cxx:361
CacheBranchFx * AccessCacheBranchDiffv(const Interaction *in) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
Spline * Spl(void) const
Definition: CacheBranchFx.h:47
#define pINFO
Definition: Messenger.h:62
double XMax(void) const
Definition: Spline.h:77
const map< double, double > & Map(void) const
Definition: CacheBranchFx.h:46
E
Definition: 018_def.c:13
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
virtual double Energy(const Interaction *in) const
A simple cache branch storing the cached data in a TNtuple.
Definition: CacheBranchFx.h:37
double XMin(void) const
Definition: Spline.h:76
double QELEventGeneratorSM::FindMaxXSec2 ( const Interaction in) const
private

Definition at line 574 of file QELEventGeneratorSM.cxx.

576 {
577 // Find a cached max xsec for the specified xsec algorithm & interaction and
578 // close to the specified energy
579 
580  // get neutrino energy
581  double E = this->Energy(interaction);
582  LOG("Kinematics", pINFO) << "E = " << E;
583 
584  if(E < fEMin) {
585  LOG("Kinematics", pINFO)
586  << "Below minimum energy - Forcing explicit calculation";
587  return -1.;
588  }
589 
590  // access the the cache branch
592 
593  // if there are enough points stored in the cache buffer to build a
594  // spline, then intepolate
595  if( cb->Spl() ) {
596  if( E >= cb->Spl()->XMin() && E <= cb->Spl()->XMax()) {
597  double spl_max_xsec = cb->Spl()->Evaluate(E);
598  LOG("Kinematics", pINFO)
599  << "\nInterpolated: max xsec (E=" << E << ") = " << spl_max_xsec;
600  return spl_max_xsec;
601  }
602  LOG("Kinematics", pINFO)
603  << "Outside spline boundaries - Forcing explicit calculation";
604  return -1.;
605  }
606 
607  // if there are not enough points at the cache buffer to have a spline,
608  // look whether there is another point that is sufficiently close
609  double dE = TMath::Min(0.25, 0.05*E);
610  const map<double,double> & fmap = cb->Map();
611  map<double,double>::const_iterator iter = fmap.lower_bound(E);
612  if(iter != fmap.end()) {
613  if(TMath::Abs(E - iter->first) < dE) return iter->second;
614  }
615 
616  return -1;
617 
618 }
intermediate_table::const_iterator const_iterator
double Evaluate(double x) const
Definition: Spline.cxx:361
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
Spline * Spl(void) const
Definition: CacheBranchFx.h:47
#define pINFO
Definition: Messenger.h:62
double XMax(void) const
Definition: Spline.h:77
const map< double, double > & Map(void) const
Definition: CacheBranchFx.h:46
E
Definition: 018_def.c:13
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
virtual double Energy(const Interaction *in) const
A simple cache branch storing the cached data in a TNtuple.
Definition: CacheBranchFx.h:37
double XMin(void) const
Definition: Spline.h:76
CacheBranchFx * AccessCacheBranch2(const Interaction *in) const
void QELEventGeneratorSM::LoadConfig ( void  )
private

Definition at line 423 of file QELEventGeneratorSM.cxx.

424 {
425 
426  // Safety factor for the maximum differential cross section
427  GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 1.2) ;
428  GetParamDef( "MaxDiffv-SafetyFactor",fSafetyFacor_nu, 1.2);
429 
430  // Minimum energy for which max xsec would be cached, forcing explicit
431  // calculation for lower eneries
432  GetParamDef( "Cache-MinEnergy", fEMin, 1.00) ;
433  GetParamDef( "Threshold-Q2", fQ2Min, 2.00);
434 
435  // Maximum allowed fractional cross section deviation from maxim cross
436  // section used in rejection method
437  GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999.);
438  assert(fMaxXSecDiffTolerance>=0);
439 
440  //-- Generate kinematics uniformly over allowed phase space and compute
441  // an event weight?
442  GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false);
443 
444  // Generate nucleon in nucleus?
445  GetParamDef( "IsNucleonInNucleus", fGenerateNucleonInNucleus, true);
446 
447 
448  sm_utils = const_cast<genie::SmithMonizUtils *>(dynamic_cast<const genie::SmithMonizUtils *>( this -> SubAlg("sm_utils_algo") ) ) ;
449 }
double fQ2Min
Q2-threshold for seeking the second maximum.
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double fSafetyFactor
maxxsec -> maxxsec * safety_factor
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec>maxxsec
bool fGenerateNucleonInNucleus
generate struck nucleon in nucleus
Contains auxiliary functions for Smith-Moniz model. .
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
double QELEventGeneratorSM::MaxDiffv ( GHepRecord evrec) const
private

Definition at line 687 of file QELEventGeneratorSM.cxx.

688 {
689  LOG("Kinematics", pINFO)
690  << "Getting max. vmax(Q2)-vmin(Q2) for the rejection method";
691 
692  double max_diffv = -1;
693  Interaction * interaction = event_rec->Summary();
694 
695  LOG("Kinematics", pINFO)
696  << "Attempting to find a cached max{vmax(Q2)-vmin(Q2)} value";
697  max_diffv = this->FindMaxDiffv(interaction);
698  if(max_diffv>0) return max_diffv;
699 
700  LOG("Kinematics", pINFO)
701  << "Attempting to compute the max{vmax(Q2)-vmin(Q2)} value";
702  max_diffv = this->ComputeMaxDiffv(interaction);
703  if(max_diffv>0) {
704  LOG("Kinematics", pINFO) << "max{vmax(Q2)-vmin(Q2)} = " << max_diffv;
705  this->CacheMaxDiffv(interaction, max_diffv);
706  return max_diffv;
707  }
708 
709  LOG("Kinematics", pNOTICE)
710  << "Can not generate event kinematics (max{vmax(Q2)-vmin(Q2);E}<=0)";
711  // xsec for selected kinematics = 0
712  event_rec->SetDiffXSec(0,kPSNull);
713  // switch on error flag
714  event_rec->EventFlags()->SetBitNumber(kKineGenErr, true);
715  // reset 'trust' bits
716  interaction->ResetBit(kISkipProcessChk);
717  interaction->ResetBit(kISkipKinematicChk);
718  // throw exception
720  exception.SetReason("kinematics generation: max{vmax(Q2)-vmin(Q2);E}<=0");
721  exception.SwitchOnFastForward();
722  throw exception;
723 
724  return 0;
725 }
double ComputeMaxDiffv(const Interaction *in) const
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
#define pINFO
Definition: Messenger.h:62
void CacheMaxDiffv(const Interaction *in, double xsec) const
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
double FindMaxDiffv(const Interaction *in) const
#define pNOTICE
Definition: Messenger.h:61
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
double QELEventGeneratorSM::MaxXSec2 ( GHepRecord evrec) const
private

Definition at line 534 of file QELEventGeneratorSM.cxx.

535 {
536  LOG("Kinematics", pINFO)
537  << "Getting max. differential xsec for the rejection method";
538 
539  double xsec_max = -1;
540  Interaction * interaction = event_rec->Summary();
541 
542  LOG("Kinematics", pINFO)
543  << "Attempting to find a cached max{dxsec/dK} value";
544  xsec_max = this->FindMaxXSec2(interaction);
545  if(xsec_max>0) return xsec_max;
546 
547  LOG("Kinematics", pINFO)
548  << "Attempting to compute the max{dxsec/dK} value";
549  xsec_max = this->ComputeMaxXSec2(interaction);
550  if(xsec_max>0) {
551  LOG("Kinematics", pINFO) << "max{dxsec/dK} = " << xsec_max;
552  this->CacheMaxXSec2(interaction, xsec_max);
553  return xsec_max;
554  }
555 
556  LOG("Kinematics", pNOTICE)
557  << "Can not generate event kinematics {K} (max_xsec({K};E)<=0)";
558  // xsec for selected kinematics = 0
559  event_rec->SetDiffXSec(0,kPSNull);
560  // switch on error flag
561  event_rec->EventFlags()->SetBitNumber(kKineGenErr, true);
562  // reset 'trust' bits
563  interaction->ResetBit(kISkipProcessChk);
564  interaction->ResetBit(kISkipKinematicChk);
565  // throw exception
567  exception.SetReason("kinematics generation: max_xsec({K};E)<=0");
568  exception.SwitchOnFastForward();
569  throw exception;
570 
571  return 0;
572 }
Summary information for an interaction.
Definition: Interaction.h:56
void CacheMaxXSec2(const Interaction *in, double xsec) 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
double ComputeMaxXSec2(const Interaction *in) const
#define pINFO
Definition: Messenger.h:62
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
double FindMaxXSec2(const Interaction *in) const
#define pNOTICE
Definition: Messenger.h:61
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
void QELEventGeneratorSM::ProcessEventRecord ( GHepRecord event_rec) const
virtual

Implements genie::EventRecordVisitorI.

Definition at line 80 of file QELEventGeneratorSM.cxx.

81 {
82  LOG("QELEvent", pINFO) << "Generating QE event kinematics...";
83 
84  if(fGenerateUniformly) {
85  LOG("QELEvent", pNOTICE)
86  << "Generating kinematics uniformly over the allowed phase space";
87  }
88 
89  // Get the interaction and set the 'trust' bits
90  Interaction * interaction = evrec->Summary();
91  const InitialState & init_state = interaction -> InitState();
92  interaction->SetBit(kISkipProcessChk);
93  interaction->SetBit(kISkipKinematicChk);
94 
95  // Skip if no hit nucleon is set
96  if(! evrec->HitNucleon())
97  {
98  LOG("QELEvent", pFATAL) << "No hit nucleon was set";
99  gAbortingInErr = true;
100  exit(1);
101  }
102 
103  // Access the target from the interaction summary
104  Target * tgt = init_state.TgtPtr();
105  GHepParticle * nucleon = evrec->HitNucleon();
106  // Store position of nucleon
107  double hitNucPos = nucleon->X4()->Vect().Mag();
108  tgt->SetHitNucPosition( hitNucPos );
109 
110  // Get the random number generators
111  RandomGen * rnd = RandomGen::Instance();
112 
113  // Access cross section algorithm for running thread
115  const EventGeneratorI * evg = rtinfo->RunningThread();
116  fXSecModel = evg->CrossSectionAlg();
117 
118  // heavy nucleus is nucleus that heavier than hydrogen and deuterium
119  bool isHeavyNucleus = tgt->A()>=3;
120 
121  sm_utils->SetInteraction(interaction);
122  // phase space for heavy nucleus is different from light one
123  fkps = isHeavyNucleus?kPSQ2vfE:kPSQ2fE;
125  // Try to calculate the maximum cross-section in kinematical limits
126  // if not pre-computed already
127  double xsec_max1 = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
128  double xsec_max2 = (fGenerateUniformly) ? -1 : (rQ2.max<fQ2Min)? 0: this->MaxXSec2(evrec);// this make correct calculation of probability
129  double vmax= isHeavyNucleus?this->MaxDiffv(evrec) : 0.;
130 
131 
132  // generate Q2, v
133  double gQ2, v, xsec;
134  unsigned int iter = 0;
135  bool accept = false;
136  TLorentzVector q;
137  while(1)
138  {
139  LOG("QELEvent", pINFO) << "Attempt #: " << iter;
140  if(iter > kRjMaxIterations)
141  {
142  LOG("QELEvent", pWARN)
143  << "Couldn't select a valid kinematics after " << iter << " iterations";
144  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
146  exception.SetReason("Couldn't select kinematics");
147  exception.SwitchOnFastForward();
148  throw exception;
149  }
150 
151  // Pick Q2 and v
152  double xsec_max = 0.;
153  double pth = rnd->RndKine().Rndm();
154  //pth < prob1/(prob1+prob2), where prob1,prob2 - probabilities to generate event in area1 (Q2<fQ2Min) and area2 (Q2>fQ2Min) which are not normalized
155  if (pth <= xsec_max1*(TMath::Min(rQ2.max, fQ2Min)-rQ2.min)/(xsec_max1*(TMath::Min(rQ2.max, fQ2Min)-rQ2.min)+xsec_max2*(rQ2.max-fQ2Min)))
156  {
157  xsec_max = xsec_max1;
158  gQ2 = (rnd->RndKine().Rndm() * (TMath::Min(rQ2.max, fQ2Min)-rQ2.min)) + rQ2.min;
159  }
160  else
161  {
162  xsec_max = xsec_max2;
163  gQ2 = (rnd->RndKine().Rndm() * (rQ2.max-fQ2Min)) + fQ2Min;
164  }
165  Range1D_t rv = sm_utils->vQES_SM_lim(gQ2);
166  // for nuclei heavier than deuterium generate energy transfer in defined energy interval
167  v = 0.;
168  if (isHeavyNucleus)
169  {
170  v = vmax * rnd->RndKine().Rndm();
171  if (v > (rv.max-rv.min))
172  {
173  continue;
174  }
175  }
176  v += rv.min;
177 
178  Kinematics * kinematics = interaction->KinePtr();
179  kinematics->SetKV(kKVQ2, gQ2);
180  kinematics->SetKV(kKVv, v);
181  xsec = fXSecModel->XSec(interaction, fkps);
182 
183  //-- Decide whether to accept the current kinematics
184  if(!fGenerateUniformly)
185  {
186  this->AssertXSecLimits(interaction, xsec, xsec_max);
187  double t = xsec_max * rnd->RndKine().Rndm();
188 
189 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
190  LOG("QELEvent", pDEBUG)<< "xsec= " << xsec << ", J= " << J << ", Rnd= " << t;
191 #endif
192  accept = (t < xsec);
193  }
194  else
195  {
196  accept = (xsec>0);
197  }
198 
199  // If the generated kinematics are accepted, finish-up module's job
200  if(accept)
201  {
202  interaction->ResetBit(kISkipProcessChk);
203  interaction->ResetBit(kISkipKinematicChk);
204  break;
205  }
206  iter++;
207  }
208 
209  // Z-frame is frame where momentum transfer is (v,0,0,qv)
210  double qv = TMath::Sqrt(v*v+gQ2);
211  TLorentzVector transferMom(0, 0, qv, v);
212 
213  Range1D_t rkF = sm_utils->kFQES_SM_lim(gQ2, v);
214  double kF = (rnd->RndKine().Rndm() * (rkF.max-rkF.min)) + rkF.min;
215 
216  // Momentum of initial neutrino in LAB frame
217  TLorentzVector * tempTLV = evrec->Probe()->GetP4();
218  TLorentzVector neutrinoMom = *tempTLV;
219  delete tempTLV;
220 
221  // define all angles in Z frame
222  double theta = neutrinoMom.Vect().Theta();
223  double phi = neutrinoMom.Vect().Phi();
224  double theta_k = sm_utils-> GetTheta_k(v, qv);
225  double costheta_k = TMath::Cos(theta_k);
226  double sintheta_k = TMath::Sin(theta_k);
227  double E_p; //energy of initial nucleon
228  double theta_p = sm_utils-> GetTheta_p(kF, v, qv, E_p);
229  double costheta_p = TMath::Cos(theta_p);
230  double sintheta_p = TMath::Sin(theta_p);
231  double fi_p = 2 * TMath::Pi() * rnd->RndGen().Rndm();
232  double cosfi_p = TMath::Cos(fi_p);
233  double sinfi_p = TMath::Sin(fi_p);
234  double psi = 2 * TMath::Pi() * rnd->RndGen().Rndm();
235 
236  // 4-momentum of initial neutrino in Z-frame
237  TLorentzVector neutrinoMomZ(neutrinoMom.P()*sintheta_k, 0, neutrinoMom.P()*costheta_k, neutrinoMom.E());
238 
239  // 4-momentum of final lepton in Z-frame
240  TLorentzVector outLeptonMom = neutrinoMomZ-transferMom;
241 
242  // 4-momentum of initial nucleon in Z-frame
243  TLorentzVector inNucleonMom(kF*sintheta_p*cosfi_p, kF*sintheta_p*sinfi_p, kF*costheta_p, E_p);
244 
245  // 4-momentum of final nucleon in Z-frame
246  TLorentzVector outNucleonMom = inNucleonMom+transferMom;
247 
248 
249  // Rotate all vectors from Z frame to LAB frame
250  TVector3 yvec (0.0, 1.0, 0.0);
251  TVector3 zvec (0.0, 0.0, 1.0);
252  TVector3 unit_nudir = neutrinoMom.Vect().Unit();
253 
254  outLeptonMom.Rotate(theta-theta_k, yvec);
255  outLeptonMom.Rotate(phi, zvec);
256 
257  inNucleonMom.Rotate(theta-theta_k, yvec);
258  inNucleonMom.Rotate(phi, zvec);
259 
260  outNucleonMom.Rotate(theta-theta_k, yvec);
261  outNucleonMom.Rotate(phi, zvec);
262 
263  outLeptonMom.Rotate(psi, unit_nudir);
264  inNucleonMom.Rotate(psi, unit_nudir);
265  outNucleonMom.Rotate(psi, unit_nudir);
266 
267  // set 4-momentum of struck nucleon
268  TLorentzVector * p4 = tgt->HitNucP4Ptr();
269  p4->SetPx( inNucleonMom.Px() );
270  p4->SetPy( inNucleonMom.Py() );
271  p4->SetPz( inNucleonMom.Pz() );
272  p4->SetE ( inNucleonMom.E() );
273 
274  int rpdgc = interaction->RecoilNucleonPdg();
275  assert(rpdgc);
276  double gW = PDGLibrary::Instance()->Find(rpdgc)->Mass();
277  LOG("QELEvent", pNOTICE) << "Selected: W = "<< gW;
278  double M = init_state.Tgt().HitNucP4().M();
279  double E = init_state.ProbeE(kRfHitNucRest);
280 
281  // (W,Q2) -> (x,y)
282  double gx=0, gy=0;
283  kinematics::WQ2toXY(E,M,gW,gQ2,gx,gy);
284 
285  // lock selected kinematics & clear running values
286  interaction->KinePtr()->SetQ2(gQ2, true);
287  interaction->KinePtr()->SetW (gW, true);
288  interaction->KinePtr()->Setx (gx, true);
289  interaction->KinePtr()->Sety (gy, true);
290  interaction->KinePtr()->ClearRunningValues();
291 
292  // set the cross section for the selected kinematics
293  evrec->SetDiffXSec(xsec,fkps);
295  {
296  double vol = sm_utils->PhaseSpaceVolume(fkps);
297  double totxsec = evrec->XSec();
298  double wght = (vol/totxsec)*xsec;
299  LOG("QELEvent", pNOTICE) << "Kinematics wght = "<< wght;
300 
301  // apply computed weight to the current event weight
302  wght *= evrec->Weight();
303  LOG("QELEvent", pNOTICE) << "Current event wght = " << wght;
304  evrec->SetWeight(wght);
305  }
306  TLorentzVector x4l(*(evrec->Probe())->X4());
307 
308  // Add the final-state lepton to the event record
309  evrec->AddParticle(interaction->FSPrimLeptonPdg(), kIStStableFinalState, evrec->ProbePosition(),-1,-1,-1, outLeptonMom, x4l);
310 
311  // Set its polarization
313 
314  GHepStatus_t ist;
316  ist = kIStStableFinalState;
317  else
319 
320  GHepParticle outNucleon(interaction->RecoilNucleonPdg(), ist, evrec->HitNucleonPosition(),-1,-1,-1, outNucleonMom , x4l);
321  evrec->AddParticle(outNucleon);
322 
323  // Store struck nucleon momentum
324  LOG("QELEvent",pNOTICE) << "pn: " << inNucleonMom.X() << ", " <<inNucleonMom.Y() << ", " <<inNucleonMom.Z() << ", " <<inNucleonMom.E();
325  nucleon->SetMomentum(inNucleonMom);
327 
328  // skip if not a nuclear target
329  if(evrec->Summary()->InitState().Tgt().IsNucleus())
330  // add a recoiled nucleus remnant
331  this->AddTargetNucleusRemnant(evrec);
332 
333 
334  LOG("QELEvent", pINFO) << "Done generating QE event kinematics!";
335 }
void SetInteraction(const Interaction *i)
Range1D_t Q2QES_SM_lim(void) const
virtual double MaxXSec(GHepRecord *evrec) const
double fQ2Min
Q2-threshold for seeking the second maximum.
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
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.
int RecoilNucleonPdg(void) const
recoil nucleon pdg
int A(void) const
Definition: Target.h:70
A simple [min,max] interval for doubles.
Definition: Range1.h:42
virtual void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
#define pFATAL
Definition: Messenger.h:56
bool IsNucleus(void) const
Definition: Target.cxx:272
double MaxDiffv(GHepRecord *evrec) const
Defines the EventGeneratorI interface.
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
void SetHitNucPosition(double r)
Definition: Target.cxx:210
void SetMomentum(const TLorentzVector &p4)
Range1D_t vQES_SM_lim(double Q2) const
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
int FSPrimLeptonPdg(void) const
final state primary lepton pdg
double GetBindingEnergy(void) const
Summary information for an interaction.
Definition: Interaction.h:56
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:91
Range1D_t kFQES_SM_lim(double nu, double Q2) 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
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
Definition: KineUtils.cxx:1119
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
bool fGenerateNucleonInNucleus
generate struck nucleon in nucleus
double MaxXSec2(GHepRecord *evrec) const
#define pINFO
Definition: Messenger.h:62
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:50
#define pWARN
Definition: Messenger.h:60
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
void SetKV(KineVar_t kv, double value)
Definition: Kinematics.cxx:335
void SetRemovalEnergy(double Erm)
static RunningThreadInfo * Instance(void)
double max
Definition: Range1.h:53
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:57
void AddTargetNucleusRemnant(GHepRecord *evrec) const
add a recoiled nucleus remnant
double PhaseSpaceVolume(KinePhaseSpace_t ps) const
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
Target * TgtPtr(void) const
Definition: InitialState.h:67
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:80
double gQ2
Definition: gtestDISSF.cxx:55
double min
Definition: Range1.h:52
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:75
#define pNOTICE
Definition: Messenger.h:61
void ClearRunningValues(void)
Definition: Kinematics.cxx:347
const Target & Tgt(void) const
Definition: InitialState.h:66
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
const EventGeneratorI * RunningThread(void)
bool gAbortingInErr
Definition: Messenger.cxx:34
void SetPrimaryLeptonPolarization(GHepRecord *ev)
double ProbeE(RefFrame_t rf) const
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
Keep info on the event generation thread currently on charge. This is used so that event generation m...
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
enum genie::EGHepStatus GHepStatus_t
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63

Member Data Documentation

bool genie::QELEventGeneratorSM::fGenerateNucleonInNucleus
private

generate struck nucleon in nucleus

Definition at line 86 of file QELEventGeneratorSM.h.

KinePhaseSpace_t genie::QELEventGeneratorSM::fkps
mutableprivate

Definition at line 83 of file QELEventGeneratorSM.h.

double genie::QELEventGeneratorSM::fQ2Min
private

Q2-threshold for seeking the second maximum.

Definition at line 87 of file QELEventGeneratorSM.h.

double genie::QELEventGeneratorSM::fSafetyFacor_nu
private

Definition at line 89 of file QELEventGeneratorSM.h.

SmithMonizUtils* genie::QELEventGeneratorSM::sm_utils
mutableprivate

Definition at line 65 of file QELEventGeneratorSM.h.


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