GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SKKinematicsGenerator.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Chris Marshall <marshall \at pas.rochester.edu>
7  University of Rochester
8 
9  Martti Nirkko
10  University of Berne
11 */
12 //____________________________________________________________________________
13 
14 #include <cstdlib>
15 
16 #include <TMath.h>
17 
19 #include "Framework/Conventions/GBuild.h"
35 
36 using namespace genie;
37 using namespace genie::constants;
38 using namespace genie::controls;
39 using namespace genie::utils;
40 
41 //___________________________________________________________________________
43 KineGeneratorWithCache("genie::SKKinematicsGenerator")
44 {
45  //fEnvelope = 0;
46 }
47 //___________________________________________________________________________
49 KineGeneratorWithCache("genie::SKKinematicsGenerator", config)
50 {
51  //fEnvelope = 0;
52 }
53 //___________________________________________________________________________
55 {
56  //if(fEnvelope) delete fEnvelope;
57 }
58 //___________________________________________________________________________
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 }
72 //___________________________________________________________________________
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 }
246 //___________________________________________________________________________
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 }
339 
340 //___________________________________________________________________________
341 double SKKinematicsGenerator::Energy(const Interaction * interaction) const
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 }
350 //___________________________________________________________________________
352 {
353  Algorithm::Configure(config);
354  this->LoadConfig();
355 }
356 //____________________________________________________________________________
358 {
359  Algorithm::Configure(config);
360  this->LoadConfig();
361 }
362 //____________________________________________________________________________
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 }
381 //____________________________________________________________________________
virtual void SetWeight(double wght)
Definition: GHepRecord.h:130
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
void Configure(const Registry &config)
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:91
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
double fSafetyFactor
ComputeMaxXSec -&gt; ComputeMaxXSec * fSafetyFactor.
virtual void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec&gt;maxxsec
Defines the EventGeneratorI interface.
virtual double MaxXSec(GHepRecord *evrec, const int nkey=0) const
virtual double Weight(void) const
Definition: GHepRecord.h:124
double Mass(void) const
Definition: Target.cxx:224
Abstract class. Provides a data caching mechanism for for concrete implementations of the EventRecord...
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
string AsString(void) const
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
double ComputeMaxXSec(const Interaction *in) const
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
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void ProcessEventRecord(GHepRecord *event_rec) 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 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
static RunningThreadInfo * Instance(void)
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
virtual TBits * EventFlags(void) const
Definition: GHepRecord.h:117
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
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
double Energy(const Interaction *in) const
#define pNOTICE
Definition: Messenger.h:61
void ClearRunningValues(void)
Definition: Kinematics.cxx:347
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Target & Tgt(void) const
Definition: InitialState.h:66
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
void CalculateKin_AtharSingleKaon(GHepRecord *event_rec) const
const EventGeneratorI * RunningThread(void)
#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
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
Keep info on the event generation thread currently on charge. This is used so that event generation m...
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
virtual void SetDiffXSec(double xsec, KinePhaseSpace_t ps)
Definition: GHepRecord.h:133
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63