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

Generates values for the kinematic variables describing neutrino-nucleus single kaon production events. Is a concrete implementation of the EventRecordVisitorI interface. More...

#include <SKKinematicsGenerator.h>

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

Public Member Functions

 SKKinematicsGenerator ()
 
 SKKinematicsGenerator (string config)
 
 ~SKKinematicsGenerator ()
 
void ProcessEventRecord (GHepRecord *event_rec) const
 
void Configure (const Registry &config)
 
void Configure (string config)
 
void LoadConfig (void)
 
void CalculateKin_AtharSingleKaon (GHepRecord *event_rec) const
 
double ComputeMaxXSec (const Interaction *in) const
 
double Energy (const Interaction *in) const
 
- 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 Attributes

double fMinLog1MinusCosTheta
 

Additional Inherited Members

- Static Public Member Functions inherited from genie::Algorithm
static string BuildParamVectKey (const std::string &comm_name, unsigned int i)
 
static string BuildParamVectSizeKey (const std::string &comm_name)
 
static string BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)
 
static string BuildParamMatRowSizeKey (const std::string &comm_name)
 
static string BuildParamMatColSizeKey (const std::string &comm_name)
 
- Protected Member Functions inherited from genie::KineGeneratorWithCache
 KineGeneratorWithCache ()
 
 KineGeneratorWithCache (string name)
 
 KineGeneratorWithCache (string name, string config)
 
 ~KineGeneratorWithCache ()
 
virtual double ComputeMaxXSec (const Interaction *in, const int nkey) const
 
virtual double MaxXSec (GHepRecord *evrec, const int nkey=0) const
 
virtual double FindMaxXSec (const Interaction *in, const int nkey=0) const
 
virtual void CacheMaxXSec (const Interaction *in, double xsec, const int nkey=0) const
 
virtual CacheBranchFxAccessCacheBranch (const Interaction *in, const int nkey=0) 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
 
template<class T >
int GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
 Handle to load matrix of parameters. More...
 
template<class T >
int GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const
 
int GetParamMatKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
 
int AddTopRegistry (Registry *rp, bool owns=true)
 add registry with top priority, also update ownership More...
 
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership More...
 
int MergeTopRegistry (const Registry &r)
 
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships. More...
 
- Protected Attributes inherited from genie::KineGeneratorWithCache
const XSecAlgorithmIfXSecModel
 
double fSafetyFactor
 ComputeMaxXSec -> ComputeMaxXSec * fSafetyFactor. More...
 
std::vector< double > vSafetyFactors
 MaxXSec -> MaxXSec * fSafetyFactors[nkey]. More...
 
int fNumOfSafetyFactors
 Number of given safety factors. More...
 
std::vector< string > vInterpolatorTypes
 Type of interpolator for each key in a branch. More...
 
int fNumOfInterpolatorTypes
 Number of given interpolators types. 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< bool > fOwnerships
 ownership for every registry in fConfVect More...
 
AlgStatus_t fStatus
 algorithm execution status More...
 
AlgMapfOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool) More...
 

Detailed Description

Generates values for the kinematic variables describing neutrino-nucleus single kaon production events. Is a concrete implementation of the EventRecordVisitorI interface.

Author
Chris Marshall <marshall pas.rochester.edu> University of Rochester
Created:
October 03, 2014
License:
Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 28 of file SKKinematicsGenerator.h.

Constructor & Destructor Documentation

SKKinematicsGenerator::SKKinematicsGenerator ( )

Definition at line 42 of file SKKinematicsGenerator.cxx.

42  :
43 KineGeneratorWithCache("genie::SKKinematicsGenerator")
44 {
45  //fEnvelope = 0;
46 }
SKKinematicsGenerator::SKKinematicsGenerator ( string  config)

Definition at line 48 of file SKKinematicsGenerator.cxx.

48  :
49 KineGeneratorWithCache("genie::SKKinematicsGenerator", config)
50 {
51  //fEnvelope = 0;
52 }
SKKinematicsGenerator::~SKKinematicsGenerator ( )

Definition at line 54 of file SKKinematicsGenerator.cxx.

55 {
56  //if(fEnvelope) delete fEnvelope;
57 }

Member Function Documentation

void SKKinematicsGenerator::CalculateKin_AtharSingleKaon ( GHepRecord event_rec) const

Definition at line 73 of file SKKinematicsGenerator.cxx.

References genie::KineGeneratorWithCache::AssertXSecLimits(), genie::Kinematics::ClearRunningValues(), genie::GHepRecord::EventFlags(), genie::Interaction::ExclTag(), genie::KineGeneratorWithCache::fGenerateUniformly, genie::PDGLibrary::Find(), fMinLog1MinusCosTheta, genie::Interaction::FSPrimLeptonPdg(), genie::KineGeneratorWithCache::fXSecModel, genie::InitialState::GetProbeP4(), genie::Target::HitNucP4(), genie::Interaction::InitState(), genie::Interaction::InitStatePtr(), genie::RandomGen::Instance(), genie::PDGLibrary::Instance(), genie::utils::mec::J(), genie::Interaction::KinePtr(), genie::kISkipKinematicChk, genie::kISkipProcessChk, genie::kKineGenErr, genie::kKVctl, genie::kKVphikq, genie::kKVSelctl, genie::kKVSelphikq, genie::kKVSelTk, genie::kKVSelTl, genie::kKVTk, genie::kKVTl, genie::constants::kPi, genie::kPSTkTlctl, genie::kRfHitNucRest, genie::controls::kRjMaxIterations, LOG, genie::Target::Mass(), genie::KineGeneratorWithCache::MaxXSec(), pDEBUG, pNOTICE, pWARN, genie::utils::kinematics::Q2(), genie::RandomGen::RndKine(), genie::GHepRecord::SetDiffXSec(), genie::Kinematics::SetKV(), genie::Kinematics::SetQ2(), genie::exceptions::EVGThreadException::SetReason(), genie::Kinematics::SetW(), genie::GHepRecord::SetWeight(), genie::Kinematics::Setx(), genie::Kinematics::Sety(), genie::XclsTag::StrangeHadronPdg(), genie::GHepRecord::Summary(), genie::exceptions::EVGThreadException::SwitchOnFastForward(), genie::InitialState::Tgt(), genie::GHepRecord::Weight(), and genie::XSecAlgorithmI::XSec().

Referenced by ProcessEventRecord().

74 {
75  // Get the Primary Interacton object
76  Interaction * interaction = evrec->Summary();
77  interaction->SetBit(kISkipProcessChk);
78  interaction->SetBit(kISkipKinematicChk);
79 
80  // Initialise a random number generator
82 
83  //-- For the subsequent kinematic selection with the rejection method:
84  // Calculate the max differential cross section or retrieve it from the
85  // cache. Throw an exception and quit the evg thread if a non-positive
86  // value is found.
87  // If the kinematics are generated uniformly over the allowed phase
88  // space the max xsec is irrelevant
89  double xsec_max = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
90 
91  // Determine lepton and kaon masses
92  int leppdg = interaction->FSPrimLeptonPdg();
93  const TLorentzVector pnuc4 = interaction->InitState().Tgt().HitNucP4(); // 4-momentum of struck nucleon in lab frame
94  TVector3 beta = pnuc4.BoostVector();
95  TLorentzVector P4_nu = *(interaction->InitStatePtr()->GetProbeP4(kRfHitNucRest)); // struck nucleon rest frame
96 
97  double enu = P4_nu.E(); // in nucleon rest frame
98  int kaon_pdgc = interaction->ExclTag().StrangeHadronPdg();
99  double mk = PDGLibrary::Instance()->Find(kaon_pdgc)->Mass();
100  double ml = PDGLibrary::Instance()->Find(leppdg)->Mass();
101 
102  // Maximum possible kinetic energy
103  const double Tkmax = enu - mk - ml;
104  const double Tlmax = enu - mk - ml;
105 
106  // Tkmax <= 0 means we are below threshold for sure
107  if( Tkmax <= 0.0 ) {
108  LOG("SKKinematics", pWARN) << "No available phase space";
109  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
111  exception.SetReason("No available phase space");
112  exception.SwitchOnFastForward();
113  throw exception;
114  }
115 
116  const double Tkmin = 0.0;
117  const double Tlmin = 0.0;
118  // for performance, use log( 1 - cos(theta_lepton) ) instead of cos(theta_lepton) because it is sharply peaked near 1.0
119  const double xmin = fMinLog1MinusCosTheta; // set in LoadConfig
120  const double xmax = 0.69314718056; // log(2) is physical boundary
121  const double phikqmin = 0.0;
122  const double phikqmax = 2.0 * kPi;
123  const double dtk = Tkmax - Tkmin;
124  const double dtl = Tlmax - Tlmin;
125  const double dx = xmax - xmin;
126  const double dphikq = phikqmax - phikqmin;
127 
128  //------ Try to select a valid tk, tl, costhetal, phikq quadruplet
129 
130  unsigned int iter = 0;
131  bool accept = false;
132  double xsec = -1; // diff XS
133  double tk = -1; // kaon kinetic energy
134  double tl = -1; // lepton kinetic energy
135  double costhetal = -1; // cosine of lepton angle
136  double phikq = -1; // azimuthal angle between kaon and q vector
137 
138  while(1) {
139  iter++;
140  if(iter > kRjMaxIterations) {
141  LOG("SKKinematics", pWARN)
142  << "*** Could not select a valid (tk, tl, costhetal) triplet after "
143  << iter << " iterations";
144  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
146  exception.SetReason("Couldn't select kinematics");
147  exception.SwitchOnFastForward();
148  throw exception;
149  }
150 
151  if(fGenerateUniformly) {
152  tk = Tkmin + dtk * rnd->RndKine().Rndm();
153  tl = Tlmin + dtl * rnd->RndKine().Rndm();
154  double x = xmin + dx * rnd->RndKine().Rndm(); // log(1-costheta)
155  costhetal = 1.0 - TMath::Exp(x);
156  phikq = phikqmin + dphikq * rnd->RndKine().Rndm();
157  } else {
158  tk = Tkmin + dtk * rnd->RndKine().Rndm();
159  tl = Tlmin + dtl * rnd->RndKine().Rndm();
160  double x = xmin + dx * rnd->RndKine().Rndm(); // log(1-costheta)
161  costhetal = 1.0 - TMath::Exp(x);
162  phikq = phikqmin + dphikq * rnd->RndKine().Rndm();
163  }
164 
165  LOG("SKKinematics", pDEBUG) << "Trying: Tk = " << tk << ", Tl = " << tl << ", cosThetal = " << costhetal << ", phikq = " << phikq;
166 
167  // nucleon rest frame! these need to be boosted to the lab frame before they become actual particles
168  interaction->KinePtr()->SetKV(kKVTk, tk);
169  interaction->KinePtr()->SetKV(kKVTl, tl);
170  interaction->KinePtr()->SetKV(kKVctl, costhetal);
171  interaction->KinePtr()->SetKV(kKVphikq, phikq);
172 
173  // lorentz invariant stuff, but do all the calculations in the nucleon rest frame
174  double el = tl + ml;
175  double pl = TMath::Sqrt(el*el - ml*ml);
176  double M = interaction->InitState().Tgt().Mass();
177  TVector3 lepton_3vector = TVector3(0,0,0);
178  lepton_3vector.SetMagThetaPhi(pl,TMath::ACos(costhetal),0.0);
179  TLorentzVector P4_lep( lepton_3vector, tl+ml );
180  TLorentzVector q = P4_nu - P4_lep;
181  double Q2 = -q.Mag2();
182  double xbj = Q2/(2*M*q.E());
183  double y = q.E()/P4_nu.E();
184  double W2 = (pnuc4+q).Mag2();
185 
186 
187  // computing cross section for the current kinematics
188  xsec = fXSecModel->XSec(interaction, kPSTkTlctl);
189 
190  //-- decide whether to accept the current kinematics
191  if(!fGenerateUniformly) {
192  // Jacobian is 1-costheta for x = log(1-costheta)
193  double max = xsec_max;
194  double t = max * rnd->RndKine().Rndm();
195  double J = TMath::Abs(1. - costhetal);
196 
197  this->AssertXSecLimits(interaction, xsec*J, max);
198 
199  if( xsec*J > xsec_max ) { // freak out if this happens
200  LOG("SKKinematics", pWARN)
201  << "!!!!!!XSEC ABOVE MAX!!!!! xsec= " << xsec << ", J= " << J << ", xsec*J = " << xsec*J << " max= " << xsec_max;
202  }
203 
204  accept = (t< J*xsec);
205  }
206  else {
207  accept = (xsec>0);
208  }
209 
210  //-- If the generated kinematics are accepted, finish-up module's job
211  if(accept) {
212 
213  // calculate the stuff
214 
215  // for uniform kinematics, compute an event weight as
216  // wght = (phase space volume)*(differential xsec)/(event total xsec)
217  if(fGenerateUniformly) {
218  double wght = 1.0; // change this
219  wght *= evrec->Weight();
220  LOG("SKKinematics", pNOTICE) << "Current event wght = " << wght;
221  evrec->SetWeight(wght);
222  }
223  LOG("SKKinematics", pWARN) << "\nLepton energy (rest frame) = " << el << " kaon = " << tl + mk;
224 
225  // reset bits
226  interaction->ResetBit(kISkipProcessChk);
227  interaction->ResetBit(kISkipKinematicChk);
228 
229  interaction->KinePtr()->SetKV(kKVSelTk, tk); // nucleon rest frame
230  interaction->KinePtr()->SetKV(kKVSelTl, tl); // nucleon rest frame
231  interaction->KinePtr()->SetKV(kKVSelctl, costhetal); // nucleon rest frame
232  interaction->KinePtr()->SetKV(kKVSelphikq, phikq); // nucleon rest frame
233  interaction->KinePtr()->SetQ2(Q2, true);
234  interaction->KinePtr()->SetW(TMath::Sqrt(W2), true);
235  interaction->KinePtr()->Setx( xbj, true );
236  interaction->KinePtr()->Sety( y, true );
237  interaction->KinePtr()->ClearRunningValues();
238 
239  // set the cross section for the selected kinematics
240  evrec->SetDiffXSec(xsec*(1.0-costhetal),kPSTkTlctl); // phase space is really log(1-costheta)
241 
242  return;
243  }
244  }// iterations
245 }
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
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 void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
virtual double MaxXSec(GHepRecord *evrec, const int nkey=0) const
double Mass(void) const
Definition: Target.cxx:224
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
Summary information for an interaction.
Definition: Interaction.h:56
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:91
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
int StrangeHadronPdg(void) const
Definition: XclsTag.h:55
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
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 double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
void SetKV(KineVar_t kv, double value)
Definition: Kinematics.cxx:335
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
const InitialState & InitState(void) const
Definition: Interaction.h:69
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
#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
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
#define pDEBUG
Definition: Messenger.h:63
double SKKinematicsGenerator::ComputeMaxXSec ( const Interaction in) const
virtual

Implements genie::KineGeneratorWithCache.

Definition at line 247 of file SKKinematicsGenerator.cxx.

References genie::Interaction::AsString(), genie::Interaction::ExclTag(), genie::PDGLibrary::Find(), fMinLog1MinusCosTheta, genie::KineGeneratorWithCache::fSafetyFactor, genie::Interaction::FSPrimLeptonPdg(), genie::KineGeneratorWithCache::fXSecModel, genie::Algorithm::Id(), genie::Interaction::InitState(), genie::PDGLibrary::Instance(), genie::Interaction::KinePtr(), genie::kKVctl, genie::kKVphikq, genie::kKVTk, genie::kKVTl, genie::constants::kPi, genie::kPSTkTlctl, genie::kRfHitNucRest, LOG, pDEBUG, pINFO, genie::InitialState::ProbeE(), genie::Kinematics::SetKV(), SLOG, genie::XclsTag::StrangeHadronPdg(), and genie::XSecAlgorithmI::XSec().

248 {
249 // Computes the maximum differential cross section in the requested phase
250 // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
251 // method and the value is cached at a circular cache branch for retrieval
252 // during subsequent event generation.
253 
254 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
255  SLOG("SKKinematics", pDEBUG)
256  << "Scanning the allowed phase space {K} for the max(dxsec/d{K})";
257 #endif
258 
259  double max_xsec = 0;
260  double max_tk = -1;
261  double max_tl = -1;
262  double max_ctl = -1;
263 
264  const int Ntk = 100;
265  const int Ntl = 100;
266  const int Nctl = 100;
267  // don't do phi_kq -- the maximum will always occur at phi_kq = pi
268 
269  int leppdg = in->FSPrimLeptonPdg();
270  double enu = in->InitState().ProbeE(kRfHitNucRest); // Enu in nucleon rest frame
271  int kaon_pdgc = in->ExclTag().StrangeHadronPdg();
272  double mk = PDGLibrary::Instance()->Find(kaon_pdgc)->Mass();
273  double ml = PDGLibrary::Instance()->Find(leppdg)->Mass();
274 
275  const double Tkmax = enu - mk - ml;
276  const double Tlmax = enu - mk - ml;
277  const double Tkmin = 0.0;
278  const double Tlmin = 0.0;
279  // cross section is sharply peaked at forward lepton
280  // faster to scan over log(1 - cos(theta)) = x
281  const double xmin = fMinLog1MinusCosTheta; // set in LoadConfig
282  const double xmax = 0.69314718056; // physical limit -- this is log(2)
283  const double dtk = (Tkmax - Tkmin)/Ntk;
284  const double dtl = (Tlmax - Tlmin)/Ntl;
285  const double dx = (xmax - xmin)/Nctl;
286 
287  // XS is 0 when the kinetic energies are == 0, so start at 1
288  for(int i = 1; i < Ntk; i++) {
289  double tk = Tkmin + dtk*i;
290  for(int j = 1; j < Ntl; j++) {
291  double tl = Tlmin + dtl*j;
292  for(int k = 0; k < Nctl; k++) {
293  double log_1_minus_cosine_theta_lepton = xmin + dx*k;
294  // XS takes cos(theta_lepton), we are scanning over log(1-that) to save time, convert to cos(theta_lepton) here
295  double ctl = 1.0 - TMath::Exp(log_1_minus_cosine_theta_lepton);
296 
297  // The cross section is 4D, but the maximum differential cross section always occurs at phi_kq = pi (or -pi)
298  // this is because there is more phase space for the kaon when it is opposite the lepton
299  // to save time in this loop, just set phi_kq to pi
300  double phikq = kPi;
301 
302  in->KinePtr()->SetKV(kKVTk, tk);
303  in->KinePtr()->SetKV(kKVTl, tl);
304  in->KinePtr()->SetKV(kKVctl, ctl);
305  in->KinePtr()->SetKV(kKVphikq, phikq);
306 
307  double xsec = fXSecModel->XSec(in, kPSTkTlctl);
308 
309  // xsec returned by model is d4sigma/(dtk dtl dcosthetal dphikq)
310  // convert lepton theta to log(1-costheta) by multiplying by jacobian 1 - costheta
311  xsec *= (1.0 - ctl);
312 
313  if( xsec > max_xsec ) {
314  max_xsec = xsec;
315  max_tk = tk;
316  max_tl = tl;
317  max_ctl = ctl;
318  }
319  }//ctl
320  }//tl
321  }//tk
322 
323  LOG("SKKinematics", pINFO) << "Max XSec is " << max_xsec << " for enu = " << enu << " tk = " << max_tk << " tl = " << max_tl << " cosine theta = " << max_ctl;
324 
325  // Apply safety factor, since value retrieved from the cache might
326  // correspond to a slightly different energy.
327  max_xsec *= fSafetyFactor;
328 
329 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
330  SLOG("SKKinematics", pDEBUG) << in->AsString();
331  SLOG("SKKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec;
332  SLOG("SKKinematics", pDEBUG) << "Computed using alg = " << fXSecModel->Id();
333 #endif
334 
335 
336 
337  return max_xsec;
338 }
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
double fSafetyFactor
ComputeMaxXSec -&gt; ComputeMaxXSec * fSafetyFactor.
string AsString(void) const
int FSPrimLeptonPdg(void) const
final state primary lepton pdg
int StrangeHadronPdg(void) const
Definition: XclsTag.h:55
#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
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
void SetKV(KineVar_t kv, double value)
Definition: Kinematics.cxx:335
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
const InitialState & InitState(void) const
Definition: Interaction.h:69
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:84
double ProbeE(RefFrame_t rf) const
#define pDEBUG
Definition: Messenger.h:63
void SKKinematicsGenerator::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 351 of file SKKinematicsGenerator.cxx.

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

352 {
353  Algorithm::Configure(config);
354  this->LoadConfig();
355 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void SKKinematicsGenerator::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 357 of file SKKinematicsGenerator.cxx.

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

358 {
359  Algorithm::Configure(config);
360  this->LoadConfig();
361 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
double SKKinematicsGenerator::Energy ( const Interaction in) const
virtual

Reimplemented from genie::KineGeneratorWithCache.

Definition at line 341 of file SKKinematicsGenerator.cxx.

References genie::Interaction::InitState(), genie::kRfLab, and genie::InitialState::ProbeE().

342 {
343 // Override the base class Energy() method to cache the max xsec for the
344 // neutrino energy in the LAB rather than in the hit nucleon rest frame.
345 
346  const InitialState & init_state = interaction->InitState();
347  double E = init_state.ProbeE(kRfLab);
348  return E;
349 }
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:48
void SKKinematicsGenerator::LoadConfig ( void  )

Definition at line 363 of file SKKinematicsGenerator.cxx.

References genie::KineGeneratorWithCache::fEMin, genie::KineGeneratorWithCache::fGenerateUniformly, genie::KineGeneratorWithCache::fMaxXSecDiffTolerance, fMinLog1MinusCosTheta, genie::KineGeneratorWithCache::fSafetyFactor, and genie::Algorithm::GetParamDef().

Referenced by Configure().

364 {
365  // max xsec safety factor (for rejection method) and min cached energy
366  this->GetParamDef("MaxXSec-SafetyFactor", fSafetyFactor, 1.5);
367  this->GetParamDef("Cache-MinEnergy", fEMin, 0.6);
368 
369  // Generate kinematics uniformly over allowed phase space and compute
370  // an event weight?
371  this->GetParamDef("UniformOverPhaseSpace", fGenerateUniformly, false);
372 
373  // Maximum allowed fractional cross section deviation from maxim cross
374  // section used in rejection method
375  this->GetParamDef("MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999.);
376  assert(fMaxXSecDiffTolerance>=0);
377 
378  //
379  this->GetParamDef("MaxXSec-MinLog1MinusCosTheta", fMinLog1MinusCosTheta, -20.0);
380 }
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double fSafetyFactor
ComputeMaxXSec -&gt; ComputeMaxXSec * fSafetyFactor.
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec&gt;maxxsec
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
bool GetParamDef(const RgKey &name, T &p, const T &def) const
void SKKinematicsGenerator::ProcessEventRecord ( GHepRecord event_rec) const
virtual

Implements genie::EventRecordVisitorI.

Definition at line 59 of file SKKinematicsGenerator.cxx.

References CalculateKin_AtharSingleKaon(), genie::EventGeneratorI::CrossSectionAlg(), genie::KineGeneratorWithCache::fGenerateUniformly, genie::KineGeneratorWithCache::fXSecModel, genie::RunningThreadInfo::Instance(), LOG, pNOTICE, and genie::RunningThreadInfo::RunningThread().

60 {
61  if(fGenerateUniformly) {
62  LOG("SKKinematics", pNOTICE)
63  << "Generating kinematics uniformly over the allowed phase space";
64  }
65 
66  //-- Access cross section algorithm for running thread
68  const EventGeneratorI * evg = rtinfo->RunningThread();
69  fXSecModel = evg->CrossSectionAlg();
71 }
bool fGenerateUniformly
uniform over allowed phase space + event weight?
Defines the EventGeneratorI interface.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static RunningThreadInfo * Instance(void)
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
#define pNOTICE
Definition: Messenger.h:61
void CalculateKin_AtharSingleKaon(GHepRecord *event_rec) const
const EventGeneratorI * RunningThread(void)
Keep info on the event generation thread currently on charge. This is used so that event generation m...

Member Data Documentation

double genie::SKKinematicsGenerator::fMinLog1MinusCosTheta
private

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