GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
COHKinematicsGenerator.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  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <cstdlib>
12 
13 #include <TROOT.h>
14 #include <TMath.h>
15 #include <TF2.h>
16 #include "Math/Minimizer.h"
17 #include "Math/Factory.h"
18 
20 #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::COHKinematicsGenerator")
44 {
45  fEnvelope = 0;
46 }
47 //___________________________________________________________________________
49  KineGeneratorWithCache("genie::COHKinematicsGenerator", config)
50 {
51  fEnvelope = 0;
52 }
53 //___________________________________________________________________________
55 {
56  if(fEnvelope) delete fEnvelope;
57 }
58 //___________________________________________________________________________
60 {
61  if(fGenerateUniformly) {
62  LOG("COHKinematics", 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();
70  if (fXSecModel->Id().Name() == "genie::ReinSehgalCOHPiPXSec") {
72  } else if (fXSecModel->Id().Name() == "genie::BergerSehgalCOHPiPXSec2015") {
74  } else if (fXSecModel->Id().Name() == "genie::BergerSehgalFMCOHPiPXSec2015") {
76  } else if ((fXSecModel->Id().Name() == "genie::AlvarezRusoCOHPiPXSec")) {
78  }
79  else {
80  LOG("COHKinematicsGenerator",pFATAL) <<
81  "ProcessEventRecord >> Cannot calculate kinematics for " <<
82  fXSecModel->Id().Name();
83  }
84 }
85 //___________________________________________________________________________
87 {
88  // Get the Primary Interacton object
89  Interaction * interaction = evrec->Summary();
90  interaction->SetBit(kISkipProcessChk);
91  interaction->SetBit(kISkipKinematicChk); // TODO: Why turn this off?
92 
93  // Initialise a random number generator
95 
96  //-- For the subsequent kinematic selection with the rejection method:
97  // Calculate the max differential cross section or retrieve it from the
98  // cache. Throw an exception and quit the evg thread if a non-positive
99  // value is found.
100  //
101  // TODO: We are not offering the "fGenerateUniformly" option here.
102  double xsec_max = this->MaxXSec(evrec);
103 
104  //-- Get the kinematical limits for the generated x,y
105  const KPhaseSpace & kps = interaction->PhaseSpace();
106  Range1D_t y = kps.YLim();
108  assert(y.min>0. && y.max>0. && y.min<1. && y.max<1. && y.min<y.max);
109 
110  const double ymin = y.min + kASmallNum;
111  const double ymax = y.max - kASmallNum;
112  const double dy = ymax - ymin;
113  const double Q2min = Q2.min + kASmallNum;
114  const double Q2max = Q2.max - kASmallNum;
115  const double dQ2 = Q2max - Q2min;
116 
117  unsigned int iter = 0;
118  bool accept=false;
119  double xsec=-1, gy=-1, gQ2=-1;
120 
121  while(1) {
122  iter++;
123  if(iter > kRjMaxIterations) this->throwOnTooManyIterations(iter,evrec);
124 
125  //-- Select unweighted kinematics using importance sampling method.
126  // TODO: The importance sampling envelope is not used. Currently,
127  // we just employ a standard rejection-method approach.
128 
129  gy = ymin + dy * rnd->RndKine().Rndm();
130  gQ2 = Q2min + dQ2 * rnd->RndKine().Rndm();
131 
132  LOG("COHKinematics", pINFO) <<
133  "Trying: Q^2 = " << gQ2 << ", y = " << gy; /* << ", t = " << gt; */
134 
135  interaction->KinePtr()->Sety(gy);
136  interaction->KinePtr()->SetQ2(gQ2);
137  kinematics::UpdateXFromQ2Y(interaction);
138 
139  // computing cross section for the current kinematics
140  xsec = fXSecModel->XSec(interaction, kPSQ2yfE);
141 
142  //-- decide whether to accept the current kinematics
143  accept = (xsec_max * rnd->RndKine().Rndm() < xsec);
144 
145  //-- If the generated kinematics are accepted, finish-up module's job
146  if(accept) {
147  LOG("COHKinematics", pNOTICE)
148  << "Selected: Q^2 = " << gQ2 << ", y = " << gy; /* << ", t = " << gt; */
149 
150  // the Berger-Sehgal COH cross section should be a triple differential cross section
151  // d^2xsec/dQ2dydt where t is the the square of the 4p transfer to the
152  // nucleus. The cross section used for kinematical selection should have
153  // the t-dependence integrated out. The t-dependence is of the form
154  // ~exp(-bt). Now that the x,y kinematical variables have been selected
155  // we can generate a t using the t-dependence as a PDF.
156  const InitialState & init_state = interaction->InitState();
157  double Ev = init_state.ProbeE(kRfLab);
158  double Epi = gy*Ev; // pion energy
159  double Epi2 = TMath::Power(Epi,2);
160  double pme2 = kPionMass2/Epi2;
161  double gx = interaction->KinePtr()->x(); // utils::kinematics::Q2YtoX(Ev,kNucleonMass,gQ2,gy); // gQ2 / ( 2. * gy * kNucleonMass * Ev);
162  double xME = kNucleonMass*gx/Epi;
163  double tA = 1. + xME - 0.5*pme2;
164  /* Range1D_t tl = kps.TLim(); // this becomes the bounds for t */
165  double tB = TMath::Sqrt(1.+ 2*xME) * TMath::Sqrt(1.-pme2);
166  double tmin = 2*Epi2 * (tA-tB);
167  double tmax = 2*Epi2 * (tA+tB);
168  double A = (double) init_state.Tgt().A();
169  double A13 = TMath::Power(A,1./3.);
170  double R = fRo * A13 * units::fermi; // nuclear radius
171  double R2 = TMath::Power(R,2.);
172  double b = 0.33333 * R2;
173  double tsum = (TMath::Exp(-b*tmin) - TMath::Exp(-b*tmax))/b;
174  double rt = tsum * rnd->RndKine().Rndm();
175  double gt = -1.*TMath::Log(-1.*b*rt + TMath::Exp(-1.*b*tmin))/b;
176 
177  // TODO: If we re-install the fGenerateUniformly option, we
178  // would compute the event weight here.
179 
180  // reset bits
181  interaction->ResetBit(kISkipProcessChk);
182  interaction->ResetBit(kISkipKinematicChk);
183 
184  // lock selected kinematics & clear running values
185  interaction->KinePtr()->Setx(gx, true);
186  interaction->KinePtr()->Sety(gy, true);
187  interaction->KinePtr()->Sett(gt, true);
188  interaction->KinePtr()->SetW(this->pionMass(interaction), true);
189  interaction->KinePtr()->SetQ2(gQ2, true);
190  interaction->KinePtr()->ClearRunningValues();
191 
192  // set the cross section for the selected kinematics
193  evrec->SetDiffXSec(xsec * TMath::Exp(-b * gt) / tsum, kPSQ2yfE);
194 
195  return;
196  }
197  }// iterations
198 }
199 //___________________________________________________________________________
201 {
202  // Get the Primary Interacton object
203  Interaction * interaction = evrec->Summary();
204  interaction->SetBit(kISkipProcessChk);
205  interaction->SetBit(kISkipKinematicChk);
206 
207  // Initialise a random number generator
208  RandomGen * rnd = RandomGen::Instance();
209 
210  //-- For the subsequent kinematic selection with the rejection method:
211  // Calculate the max differential cross section or retrieve it from the
212  // cache. Throw an exception and quit the evg thread if a non-positive
213  // value is found.
214  //
215  // TODO: We are not offering the "fGenerateUniformly" option here.
216  double xsec_max = this->MaxXSec(evrec);
217 
218  //-- Get the kinematical limits for the generated x,y
219  const KPhaseSpace & kps = interaction->PhaseSpace();
220  Range1D_t y = kps.YLim();
222  assert(y.min>0. && y.max>0. && y.min<1. && y.max<1. && y.min<y.max);
223 
224  const double ymin = y.min + kASmallNum;
225  const double ymax = y.max - kASmallNum;
226  const double dy = ymax - ymin;
227  const double Q2min = Q2.min + kASmallNum;
228  const double Q2max = Q2.max - kASmallNum;
229  const double dQ2 = Q2max - Q2min;
230  const double tmin = kASmallNum;
231  const double tmax = fTMax - kASmallNum; // TODO: Choose realistic t bounds
232  const double dt = tmax - tmin;
233 
234  //-- Try to select a valid (Q^2,y,t) triple.
235 
236  unsigned int iter = 0;
237  bool accept=false;
238  double xsec=-1, gy=-1, gt=-1, gQ2=-1;
239 
240  while(1) {
241  iter++;
242  if(iter > kRjMaxIterations) this->throwOnTooManyIterations(iter,evrec);
243 
244  //-- Select unweighted kinematics using importance sampling method.
245  // TODO: The importance sampling envelope is not used. Currently,
246  // we just employ a standard rejection-method approach.
247 
248  gy = ymin + dy * rnd->RndKine().Rndm();
249  gt = tmin + dt * rnd->RndKine().Rndm();
250  gQ2 = Q2min + dQ2 * rnd->RndKine().Rndm();
251 
252  LOG("COHKinematics", pINFO) <<
253  "Trying: Q^2 = " << gQ2 << ", y = " << gy << ", t = " << gt;
254 
255  interaction->KinePtr()->Sety(gy);
256  interaction->KinePtr()->Sett(gt);
257  interaction->KinePtr()->SetQ2(gQ2);
258 
259  // computing cross section for the current kinematics
260  xsec = fXSecModel->XSec(interaction, kPSxyfE);
261 
262  //-- decide whether to accept the current kinematics
263  accept = (xsec_max * rnd->RndKine().Rndm() < xsec);
264 
265  //-- If the generated kinematics are accepted, finish-up module's job
266  if(accept) {
267  LOG("COHKinematics", pNOTICE)
268  << "Selected: Q^2 = " << gQ2 << ", y = " << gy << ", t = " << gt;
269 
270  // TODO: If we re-install the fGenerateUniformly option, we
271  // would compute the event weight here.
272 
273  // reset bits
274  interaction->ResetBit(kISkipProcessChk);
275  interaction->ResetBit(kISkipKinematicChk);
276 
277  // lock selected kinematics & clear running values
278  interaction->KinePtr()->SetQ2(gQ2, true);
279  interaction->KinePtr()->Sety(gy, true);
280  interaction->KinePtr()->Sett(gt, true);
281  interaction->KinePtr()->SetW(this->pionMass(interaction), true);
282  interaction->KinePtr()->ClearRunningValues();
283 
284  // set the cross section for the selected kinematics
285  evrec->SetDiffXSec(xsec, kPSxytfE);
286 
287  return;
288  }
289  }// iterations
290 }
291 //___________________________________________________________________________
293 {
294  // Get the Primary Interacton object
295  Interaction * interaction = evrec->Summary();
296  interaction->SetBit(kISkipProcessChk);
297  interaction->SetBit(kISkipKinematicChk);
298 
299  //-- Get the random number generators
300  RandomGen * rnd = RandomGen::Instance();
301 
302  //-- For the subsequent kinematic selection with the rejection method:
303  // Calculate the max differential cross section or retrieve it from the
304  // cache. Throw an exception and quit the evg thread if a non-positive
305  // value is found.
306  // If the kinematics are generated uniformly over the allowed phase
307  // space the max xsec is irrelevant
308  double xsec_max = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
309 
310  //-- Get the kinematical limits for the generated x,y
311  const KPhaseSpace & kps = interaction->PhaseSpace();
312  Range1D_t y = kps.YLim();
313  assert(y.min>0. && y.max>0. && y.min<1. && y.max<1. && y.min<y.max);
314 
315  const double xmin = kASmallNum;
316  const double xmax = 1.- kASmallNum;
317  const double ymin = y.min + kASmallNum;
318  const double ymax = y.max - kASmallNum;
319  const double dx = xmax - xmin;
320  const double dy = ymax - ymin;
321 
322  //------ Try to select a valid x,y pair
323 
324  unsigned int iter = 0;
325  bool accept=false;
326  double xsec=-1, gx=-1, gy=-1;
327 
328  while(1) {
329  iter++;
330  if(iter > kRjMaxIterations) this->throwOnTooManyIterations(iter,evrec);
331 
332  if(fGenerateUniformly) {
333  //-- Generate a x,y pair uniformly in the kinematically allowed range.
334  gx = xmin + dx * rnd->RndKine().Rndm();
335  gy = ymin + dy * rnd->RndKine().Rndm();
336 
337  } else {
338  //-- Select unweighted kinematics using importance sampling method.
339 
340  if(iter==1) {
341  LOG("COHKinematics", pNOTICE) << "Initializing the sampling envelope";
342  double Ev = interaction->InitState().ProbeE(kRfLab);
343  fEnvelope->SetRange(xmin,ymin,xmax,ymax);
344  fEnvelope->SetParameter(0, xsec_max);
345  fEnvelope->SetParameter(1, Ev);
346  }
347 
348  // Generate W,QD2 using the 2-D envelope as PDF
349  fEnvelope->GetRandom2(gx,gy);
350  }
351 
352  LOG("COHKinematics", pINFO) << "Trying: x = " << gx << ", y = " << gy;
353 
354  interaction->KinePtr()->Setx(gx);
355  interaction->KinePtr()->Sety(gy);
356 
357  // computing cross section for the current kinematics
358  xsec = fXSecModel->XSec(interaction, kPSxyfE);
359 
360  //-- decide whether to accept the current kinematics
361  if(!fGenerateUniformly) {
362  double max = fEnvelope->Eval(gx, gy);
363  double t = max * rnd->RndKine().Rndm();
364 
365  this->AssertXSecLimits(interaction, xsec, max);
366 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
367  LOG("COHKinematics", pDEBUG)
368  << "xsec= " << xsec << ", J= 1, Rnd= " << t;
369 #endif
370  accept = (t<xsec);
371  }
372  else {
373  accept = (xsec>0);
374  }
375 
376  //-- If the generated kinematics are accepted, finish-up module's job
377  if(accept) {
378  LOG("COHKinematics", pNOTICE) << "Selected: x = "<< gx << ", y = "<< gy;
379 
380  // the Rein-Sehgal COH cross section should be a triple differential cross section
381  // d^2xsec/dxdydt where t is the the square of the 4p transfer to the
382  // nucleus. The cross section used for kinematical selection should have
383  // the t-dependence integrated out. The t-dependence is of the form
384  // ~exp(-bt). Now that the x,y kinematical variables have been selected
385  // we can generate a t using the t-dependence as a PDF.
386  const InitialState & init_state = interaction->InitState();
387  double Ev = init_state.ProbeE(kRfLab);
388  double Epi = gy*Ev; // pion energy
389  double Epi2 = TMath::Power(Epi,2);
390  double pme2 = kPionMass2/Epi2;
391  double xME = kNucleonMass*gx/Epi;
392  double tA = 1. + xME - 0.5*pme2;
393  double tB = TMath::Sqrt(1.+ 2*xME) * TMath::Sqrt(1.-pme2);
394  double tmin = 2*Epi2 * (tA-tB);
395  double tmax = 2*Epi2 * (tA+tB);
396  double A = (double) init_state.Tgt().A();
397  double A13 = TMath::Power(A,1./3.);
398  double R = fRo * A13 * units::fermi; // nuclear radius
399  double R2 = TMath::Power(R,2.);
400  double b = 0.33333 * R2;
401  double tsum = (TMath::Exp(-b*tmin) - TMath::Exp(-b*tmax))/b;
402  double rt = tsum * rnd->RndKine().Rndm();
403  double gt = -1.*TMath::Log(-1.*b*rt + TMath::Exp(-1.*b*tmin))/b;
404 
405  LOG("COHKinematics", pNOTICE)
406  << "Selected: t = "<< gt << ", from ["<< tmin << ", "<< tmax << "]";
407 
408  // for uniform kinematics, compute an event weight as
409  // wght = (phase space volume)*(differential xsec)/(event total xsec)
410  if(fGenerateUniformly) {
411  double vol = y.max-y.min; // dx=1, dt: irrelevant
412  double totxsec = evrec->XSec();
413  double wght = (vol/totxsec)*xsec;
414  LOG("COHKinematics", pNOTICE) << "Kinematics wght = "<< wght;
415 
416  // apply computed weight to the current event weight
417  wght *= evrec->Weight();
418  LOG("COHKinematics", pNOTICE) << "Current event wght = " << wght;
419  evrec->SetWeight(wght);
420  }
421 
422  // reset bits
423  interaction->ResetBit(kISkipProcessChk);
424  interaction->ResetBit(kISkipKinematicChk);
425 
426  // lock selected kinematics & clear running values
427  interaction->KinePtr()->Setx(gx, true);
428  interaction->KinePtr()->Sety(gy, true);
429  interaction->KinePtr()->Sett(gt, true);
430  interaction->KinePtr()->SetW(kPionMass, true);
431  interaction->KinePtr()->SetQ2(2*kNucleonMass*gx*gy*Ev, true);
432  interaction->KinePtr()->ClearRunningValues();
433 
434  // set the cross section for the selected kinematics
435  evrec->SetDiffXSec(xsec * TMath::Exp(-b * gt) / tsum, kPSxytfE);
436 
437  return;
438  }
439  }// iterations
440 }
441 //___________________________________________________________________________
443 {
444 
445  LOG("COHKinematics", pNOTICE) << "Using AlvarezRuso Model";
446  // Get the Primary Interacton object
447  Interaction * interaction = evrec->Summary();
448  interaction->SetBit(kISkipProcessChk);
449  interaction->SetBit(kISkipKinematicChk);
450 
451  // Initialise a random number generator
452  RandomGen * rnd = RandomGen::Instance();
453 
454  //-- For the subsequent kinematic selection with the rejection method:
455  // Calculate the max differential cross section or retrieve it from the
456  // cache. Throw an exception and quit the evg thread if a non-positive
457  // value is found.
458  // If the kinematics are generated uniformly over the allowed phase
459  // space the max xsec is irrelevant
460  double xsec_max = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
461 
462  //Set up limits of integration variables
463  // Primary lepton energy
464  const double E_l_min = interaction->FSPrimLepton()->Mass();
465  const double E_l_max = interaction->InitStatePtr()->GetProbeP4(kRfLab)->E() - kPionMass;
466  // Primary lepton angle with respect to the beam axis
467  const double ctheta_l_min = 0.4;
468  const double ctheta_l_max = 1.0 - kASmallNum;
469  // Pion angle with respect to the beam axis
470  const double ctheta_pi_min = 0.4;
471  const double ctheta_pi_max = 1.0 - kASmallNum;
472  // Pion angle transverse to the beam axis
473  const double phi_min = kASmallNum;
474  const double phi_max = (2.0 * kPi) - kASmallNum;
475  //
476  const double d_E_l = E_l_max - E_l_min;
477  const double d_ctheta_l = ctheta_l_max - ctheta_l_min;
478  const double d_ctheta_pi = ctheta_pi_max - ctheta_pi_min;
479  const double d_phi = phi_max - phi_min;
480 
481  //------ Try to select a valid set of kinematics
482  unsigned int iter = 0;
483  bool accept=false;
484  double xsec=-1, g_E_l=-1, g_theta_l=-1, g_phi_l=-1, g_theta_pi=-1, g_phi_pi=-1;
485  double g_ctheta_l, g_ctheta_pi;
486 
487  while(1) {
488  iter++;
489  if(iter > kRjMaxIterations) this->throwOnTooManyIterations(iter,evrec);
490 
491  //Select kinematic point
492  g_E_l = E_l_min + d_E_l * rnd->RndKine().Rndm();
493  g_ctheta_l = ctheta_l_min + d_ctheta_l * rnd->RndKine().Rndm();
494  g_ctheta_pi = ctheta_pi_min + d_ctheta_pi * rnd->RndKine().Rndm();
495  g_phi_l = phi_min + d_phi * rnd->RndKine().Rndm();
496  // random phi is relative to phi_l
497  g_phi_pi = g_phi_l + (phi_min + d_phi * rnd->RndKine().Rndm());
498  g_theta_l = TMath::ACos(g_ctheta_l);
499  g_theta_pi = TMath::ACos(g_ctheta_pi);
500 
501  LOG("COHKinematics", pINFO) << "Trying: Lep(" <<g_E_l << ", " <<
502  g_theta_l << ", " << g_phi_l << ") Pi(" <<
503  g_theta_pi << ", " << g_phi_pi << ")";
504 
505  this->SetKinematics(g_E_l, g_theta_l, g_phi_l, g_theta_pi, g_phi_pi,
506  interaction, interaction->KinePtr());
507 
508  // computing cross section for the current kinematics
509  xsec = fXSecModel->XSec(interaction,kPSElOlOpifE) / (1E-38 * units::cm2);
510 
511  if (!fGenerateUniformly) {
512  //-- decide whether to accept the current kinematics
513  double t = xsec_max * rnd->RndKine().Rndm();
514 
515  LOG("COHKinematics", pINFO) << "Got: xsec = " << xsec << ", t = " <<
516  t << " (max_xsec = " << xsec_max << ")";
517 
518  this->AssertXSecLimits(interaction, xsec, xsec_max);
519 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
520  LOG("COHKinematics", pDEBUG)
521  << "xsec= " << xsec << ", J= 1, Rnd= " << t;
522 #endif
523  accept = (t<xsec);
524  }
525  else {
526  accept = (xsec>0);
527  }
528 
529  //-- If the generated kinematics are accepted, finish-up module's job
530  if(accept) {
531  LOG("COHKinematics", pNOTICE) << "Selected: Lepton(" <<
532  g_E_l << ", " << g_theta_l << ", " <<
533  g_phi_l << ") Pion(" << g_theta_pi << ", " << g_phi_pi << ")";
534 
535  double E_l = g_E_l;
536  double theta_l = g_theta_l;
537  double theta_pi = g_theta_pi;
538  double phi_l = g_phi_l;
539  double phi_pi = g_phi_pi;
540  const TLorentzVector P4_nu = *(interaction->InitStatePtr()->GetProbeP4(kRfLab));
541  double E_nu = P4_nu.E();
542  double E_pi= E_nu-E_l;
543  double m_l = interaction->FSPrimLepton()->Mass();
544  double m_pi = this->pionMass(interaction);
545 
546  double p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
547  TVector3 lepton_3vector = TVector3(0,0,0);
548  lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
549  TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
550 
551  double p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
552  TVector3 pion_3vector = TVector3(0,0,0);
553  pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
554  TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
555 
556  TLorentzVector q = P4_nu - P4_lep;
557  double Q2 = -q.Mag2();
558  double x = Q2/(2*E_pi*constants::kNucleonMass);
559  double y = E_pi/E_nu;
560 
561  double t = TMath::Abs( (q - P4_pion).Mag2() );
562 
563  // for uniform kinematics, compute an event weight as
564  // wght = (phase space volume)*(differential xsec)/(event total xsec)
565  if(fGenerateUniformly) {
566  // Phase space volume needs checking
567  double vol = d_E_l*d_ctheta_l*d_phi*d_ctheta_pi*d_phi;
568  double totxsec = evrec->XSec();
569  double wght = (vol/totxsec)*xsec;
570  LOG("COHKinematics", pNOTICE) << "Kinematics wght = "<< wght;
571 
572  // apply computed weight to the current event weight
573  wght *= evrec->Weight();
574  LOG("COHKinematics", pNOTICE) << "Current event wght = " << wght;
575  evrec->SetWeight(wght);
576  }
577 
578  // reset bits
579  interaction->ResetBit(kISkipProcessChk);
580  interaction->ResetBit(kISkipKinematicChk);
581  // lock selected kinematics & clear running values
582  interaction->KinePtr()->Setx(x, true);
583  interaction->KinePtr()->Sety(y, true);
584  interaction->KinePtr()->Sett(t, true);
585  interaction->KinePtr()->SetW(kPionMass, true);
586  interaction->KinePtr()->SetQ2(2*kNucleonMass*x*y*E_nu, true);
587  interaction->KinePtr()->ClearRunningValues();
588  // set the cross section for the selected kinematics
589  evrec->SetDiffXSec(xsec,kPSElOlOpifE);
590  return;
591  }
592  }//while
593 }
594 //___________________________________________________________________________
596  const double theta_l,
597  const double phi_l,
598  const double theta_pi,
599  const double phi_pi,
600  const Interaction* interaction,
601  Kinematics* kinematics) const
602 {
603  const TLorentzVector P4_nu = *(interaction->InitStatePtr()->GetProbeP4(kRfLab));
604  double E_nu = P4_nu.E();
605  double E_pi= E_nu-E_l;
606  double m_l = interaction->FSPrimLepton()->Mass();
607  double m_pi;
608  if ( interaction->ProcInfo().IsWeakCC() ) {
609  m_pi = constants::kPionMass;
610  } else {
611  m_pi = constants::kPi0Mass;
612  }
613  double p_l=0.0;
614  if (E_l > m_l) {
615  p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
616  }
617  TVector3 lepton_3vector = TVector3(0,0,0);
618  lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
619  TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
620 
621  double p_pi=0.0;
622  if (E_pi > m_pi) {
623  p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
624  }
625  TVector3 pion_3vector = TVector3(0,0,0);
626  pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
627  TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
628 
629  double Q2 = -(P4_nu-P4_lep).Mag2();
630  double x = Q2/(2*E_pi*constants::kNucleonMass);
631  double y = E_pi/E_nu;
632 
633  kinematics->Setx(x);
634  kinematics->Sety(y);
635  kinematics::UpdateWQ2FromXY(interaction);
636 
637  kinematics->SetFSLeptonP4(P4_lep );
638  kinematics->SetHadSystP4 (P4_pion); // use Hadronic System variable to store pion momentum
639 }
640 //___________________________________________________________________________
642  const double /* theta_l */ ,
643  const double /* phi_l */ ,
644  const double /* theta_pi */ ,
645  const double /* phi_pi */ ,
646  const Interaction* interaction) const
647 {
648  const TLorentzVector P4_nu = *(interaction->InitStatePtr()->GetProbeP4(kRfLab));
649  double E_nu = P4_nu.E();
650  double E_pi= E_nu-E_l;
651  double m_l = interaction->FSPrimLepton()->Mass();
652  double m_pi;
653  if ( interaction->ProcInfo().IsWeakCC() ) {
654  m_pi = constants::kPionMass;
655  }
656  else {
657  m_pi = constants::kPi0Mass;
658  }
659  if (E_l <= m_l) {
660  return false;
661  }
662  if (E_pi <= m_pi) {
663  return false;
664  }
665  return true;
666 }
667 //___________________________________________________________________________
669 {
670  // Computes the maximum differential cross section in the requested phase
671  // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
672  // method and the value is cached at a circular cache branch for retrieval
673  // during subsequent event generation.
674 
675 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
676  SLOG("COHKinematics", pDEBUG)
677  << "Scanning the allowed phase space {K} for the max(dxsec/d{K})";
678 #endif
679  double max_xsec = 0.;
680  if (fXSecModel->Id().Name() == "genie::ReinSehgalCOHPiPXSec") {
681  max_xsec = MaxXSec_ReinSehgal(in);
682  } else if ((fXSecModel->Id().Name() == "genie::BergerSehgalCOHPiPXSec2015")) {
683  max_xsec = MaxXSec_BergerSehgal(in);
684  } else if ((fXSecModel->Id().Name() == "genie::BergerSehgalFMCOHPiPXSec2015")) {
685  max_xsec = MaxXSec_BergerSehgalFM(in);
686  } else if ((fXSecModel->Id().Name() == "genie::AlvarezRusoCOHPiPXSec")) {
687  max_xsec = MaxXSec_AlvarezRuso(in);
688  }
689  else {
690  LOG("COHKinematicsGenerator",pFATAL) <<
691  "ComputeMaxXSec >> Cannot calculate max cross-section for " <<
692  fXSecModel->Id().Name();
693  }
694 
695  // Apply safety factor, since value retrieved from the cache might
696  // correspond to a slightly different energy.
697  max_xsec *= fSafetyFactor;
698 
699 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
700  SLOG("COHKinematics", pDEBUG) << in->AsString();
701  SLOG("COHKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec;
702  SLOG("COHKinematics", pDEBUG) << "Computed using alg = " << fXSecModel->Id();
703 #endif
704 
705  return max_xsec;
706 }
707 //___________________________________________________________________________
709 {
710  double max_xsec = 0;
711  const int NQ2 = 50;
712  const int Ny = 50;
713 
714  const KPhaseSpace & kps = in->PhaseSpace();
715  Range1D_t Q2r = kps.Q2Lim();
716  Q2r.max = fQ2Max;
717 
718  const double logQ2min = TMath::Log10(Q2r.min + kASmallNum);
719  const double logQ2max = TMath::Log10(Q2r.max);
720  const double dlogQ2 = (logQ2max - logQ2min) /(NQ2-1);
721 
722  for(int i=0; i<NQ2; i++) {
723  double Q2 = TMath::Power(10, logQ2min + i * dlogQ2);
724  in->KinePtr()->SetQ2(Q2);
725 
726  Range1D_t yr = kps.YLim();
727  if ((yr.max < 0) || (yr.max < yr.min) ||
728  (yr.max > 1) || (yr.min < 0)) { // forbidden kinematics
729  continue;
730  }
731  const double logymin = TMath::Log10(yr.min);
732  const double logymax = TMath::Log10(yr.max);
733  const double dlogy = (logymax - logymin) /(Ny-1);
734 
735  for(int j=0; j<Ny; j++) {
736  double gy = TMath::Power(10, logymin + j * dlogy);
737  in->KinePtr()->Sety(gy);
738 
739  /* Range1D_t tl = kps.TLim(); // TESTING! - this becomes a loop over t */
741 
742  // Note: We're not stepping through log Q^2, log y - we "unpacked"
743  double xsec = fXSecModel->XSec(in, kPSQ2yfE);
744 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
745  LOG("COHKinematics", pDEBUG)
746  << "xsec(Q2= " << Q2 << ", y= " << gy << ", t = " << gt << ") = " << xsec;
747 #endif
748  max_xsec = TMath::Max(max_xsec, xsec);
749 
750  } // y
751  } // Q2
752  return max_xsec;
753 }
754 //___________________________________________________________________________
756 {
757  double max_xsec = 0;
758  // How many sampling bins in each variable for max xsec calculation?
759  const int NQ2 = 50;
760  const int Ny = 50;
761  const int Nt = 50;
762 
763  const KPhaseSpace & kps = in->PhaseSpace();
764  Range1D_t Q2r = kps.Q2Lim();
765  Q2r.max = fQ2Max;
766 
767  const double logQ2min = TMath::Log10(Q2r.min + kASmallNum);
768  const double logQ2max = TMath::Log10(Q2r.max);
769  const double logtmin = TMath::Log10(kASmallNum);
770  const double logtmax = TMath::Log10(fTMax - kASmallNum);
771  const double dlogQ2 = (logQ2max - logQ2min) /(NQ2-1);
772  const double dlogt = (logtmax - logtmin) /(Nt-1);
773 
774  for(int i=0; i<NQ2; i++) {
775  double Q2 = TMath::Power(10, logQ2min + i * dlogQ2);
776  in->KinePtr()->SetQ2(Q2);
777 
778  Range1D_t yr = kps.YLim();
779  if ((yr.max < 0) || (yr.max < yr.min) ||
780  (yr.max > 1) || (yr.min < 0)) { // forbidden kinematics
781  continue;
782  }
783  const double logymin = TMath::Log10(yr.min);
784  const double logymax = TMath::Log10(yr.max);
785  const double dlogy = (logymax - logymin) /(Ny-1);
786 
787  for(int j=0; j<Ny; j++) {
788  double gy = TMath::Power(10, logymin + j * dlogy);
789 
790  for(int k=0; k<Nt; k++) {
791  double gt = TMath::Power(10, logtmin + k * dlogt);
792 
793  in->KinePtr()->Sety(gy);
794  in->KinePtr()->Sett(gt);
795 
796  double xsec = fXSecModel->XSec(in, kPSxyfE);
797 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
798  LOG("COHKinematics", pDEBUG)
799  << "xsec(Q2= " << Q2 << ", y= " << gy << ", t = " << gt << ") = " << xsec;
800 #endif
801  max_xsec = TMath::Max(max_xsec, xsec);
802 
803  } // t
804  } // y
805  } // Q2
806  return max_xsec;
807 }
808 //___________________________________________________________________________
810 {
811  double max_xsec = 0;
812  double Ev = in->InitState().ProbeE(kRfLab);
813 
814  const int Nx = 50;
815  const int Ny = 50;
816 
817  const KPhaseSpace & kps = in->PhaseSpace();
818  Range1D_t y = kps.YLim();
819 
820  const double logxmin = TMath::Log10(1E-5);
821  const double logxmax = TMath::Log10(1.0);
822  const double logymin = TMath::Log10(y.min);
823  const double logymax = TMath::Log10(y.max);
824 
825  const double dlogx = (logxmax - logxmin) /(Nx-1);
826  const double dlogy = (logymax - logymin) /(Ny-1);
827 
828  for(int i=0; i<Nx; i++) {
829  double gx = TMath::Power(10, logxmin + i * dlogx);
830  for(int j=0; j<Ny; j++) {
831  double gy = TMath::Power(10, logymin + j * dlogy);
832 
833  double Q2 = 2*kNucleonMass*gx*gy*Ev;
834  if(Ev>1.0 && Q2>0.01) continue;
835 
836  in->KinePtr()->Setx(gx);
837  in->KinePtr()->Sety(gy);
838 
839  double xsec = fXSecModel->XSec(in, kPSxyfE);
840 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
841  LOG("COHKinematics", pDEBUG)
842  << "xsec(x= " << gx << ", y= " << gy << ") = " << xsec;
843 #endif
844  max_xsec = TMath::Max(max_xsec, xsec);
845 
846  }//y
847  }//x
848  return max_xsec;
849 }
850 //___________________________________________________________________________
852 {
853  // Computes the maximum differential cross section in the requested phase
854  // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
855  // method and the value is cached at a circular cache branch for retrieval
856  // during subsequent event generation.
857 
858 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
859  SLOG("COHKinematics", pDEBUG)
860  << "Scanning the allowed phase space {K} for the max(dxsec/d{K})";
861 #endif
862  double max_xsec = 0.;
863  double Ev = in->InitState().ProbeE(kRfLab);
864 
865  const KPhaseSpace & kps = in->PhaseSpace();
866  Range1D_t y = kps.YLim();
867 
868  ROOT::Math::Minimizer * min = ROOT::Math::Factory::CreateMinimizer("Minuit2");
870  f.SetFactor(-1.); // Make it return negative of cross-section so we can minimize
871 
872  min->SetFunction( f );
873  min->SetMaxFunctionCalls(10000);
874  min->SetTolerance(0.05);
875 
876  const double min_el = in->FSPrimLepton()->Mass();
877  const double max_el = Ev - kPionMass;
878  const unsigned int n_el = 100;
879  const double d_el = (max_el - min_el) / double(n_el - 1);
880 
881  const double min_thetal = kASmallNum;
882  const double max_thetal = kPi / 4.0;
883  const unsigned int n_thetal = 10;
884  const double d_thetal = (max_thetal - min_thetal) / double(n_thetal - 1);
885 
886  const double min_thetapi = kASmallNum;
887  const double max_thetapi = kPi / 2.0;
888  const unsigned int n_thetapi = 10;
889  const double d_thetapi = (max_thetapi - min_thetapi) / double(n_thetapi - 1);
890 
891  //~ const double min_phipi = kPi;
892  //~ const double max_phipi = 0.5 * kPi;
893  const double min_phipi = kASmallNum;
894  const double max_phipi = 2*kPi-kASmallNum;
895  const unsigned int n_phipi = 10;
896  const double d_phipi = (max_phipi - min_phipi) / double(n_phipi - 1);
897 
898  min->SetLimitedVariable ( 0 ,"E_lep" , max_el -kASmallNum , d_el , min_el , max_el );
899  min->SetLimitedVariable ( 1 ,"theta_l" , min_thetal +kASmallNum , d_thetal , min_thetal , max_thetal );
900  min->SetLimitedVariable ( 2 ,"theta_pi" , min_thetapi+kASmallNum , d_thetapi , min_thetapi, max_thetapi );
901  min->SetLimitedVariable ( 3 ,"phi_pi" , min_phipi +kASmallNum , d_phipi , min_phipi , max_phipi );
902 
903  min->Minimize();
904  max_xsec = -min->MinValue(); //back to positive xsec
905 
906  // Apply safety factor, since value retrieved from the cache might
907  // correspond to a slightly different energy.
908  max_xsec *= fSafetyFactor;
909 
910 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
911  SLOG("COHKinematics", pDEBUG) << in->AsString();
912  SLOG("COHKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec;
913  SLOG("COHKinematics", pDEBUG) << "Computed using alg = " << fXSecModel->Id();
914 #endif
915 
916  delete min;
917 
918  return max_xsec;
919 }
920 //___________________________________________________________________________
921 double COHKinematicsGenerator::Energy(const Interaction * interaction) const
922 {
923  // Override the base class Energy() method to cache the max xsec for the
924  // neutrino energy in the LAB rather than in the hit nucleon rest frame.
925 
926  const InitialState & init_state = interaction->InitState();
927  double E = init_state.ProbeE(kRfLab);
928  return E;
929 }
930 //___________________________________________________________________________
932 {
933  double m_pi = 0.0;
934  if ( in->ProcInfo().IsWeakCC() ) {
935  m_pi = constants::kPionMass;
936  } else {
937  m_pi = constants::kPi0Mass;
938  }
939  return m_pi;
940 }
941 //___________________________________________________________________________
943  GHepRecord* evrec) const
944 {
945  LOG("COHKinematics", pWARN)
946  << "*** Could not select valid kinematics after "
947  << iters << " iterations";
948  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
950  exception.SetReason("Couldn't select kinematics");
951  exception.SwitchOnFastForward();
952  throw exception;
953 }
954 //___________________________________________________________________________
956 {
957  Algorithm::Configure(config);
958  this->LoadConfig();
959 }
960 //____________________________________________________________________________
962 {
963  Algorithm::Configure(config);
964  this->LoadConfig();
965 }
966 //____________________________________________________________________________
968 {
969  //-- COH model parameter Ro
970  GetParam( "COH-Ro", fRo );
971  //-- COH model parameter t_max for t = (q - p_pi)^2
972  GetParam( "COH-t-max", fTMax ) ;
973  //-- COH model bounds of integration for Q^2
974  GetParam( "COH-Q2-min", fQ2Min ) ;
975  GetParam( "COH-Q2-max", fQ2Max ) ;
976 
977  //-- max xsec safety factor (for rejection method) and min cached energy
978  GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 1.6 ) ;
979  GetParamDef( "Cache-MinEnergy", fEMin, -1.0 ) ;
980 
981  //-- Generate kinematics uniformly over allowed phase space and compute
982  // an event weight?
983  GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false ) ;
984 
985  //-- Maximum allowed fractional cross section deviation from maxim cross
986  // section used in rejection method
987  GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ;
988  assert(fMaxXSecDiffTolerance>=0);
989 
990  //-- Envelope employed when importance sampling is used
991  // (initialize with dummy range)
992  if(fEnvelope) delete fEnvelope;
993  fEnvelope = new TF2("CohKinEnvelope",
995  // stop ROOT from deleting this object of its own volition
996  gROOT->GetListOfFunctions()->Remove(fEnvelope);
997 }
998 //____________________________________________________________________________
double COHImportanceSamplingEnvelope(double *x, double *par)
Definition: KineUtils.cxx:1466
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
virtual void SetWeight(double wght)
Definition: GHepRecord.h:130
bool IsWeakCC(void) const
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double MaxXSec_ReinSehgal(const Interaction *in) const
double fRo
nuclear scale parameter
static const double kNucleonMass
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.
void CalculateKin_BergerSehgalFM(GHepRecord *event_rec) const
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
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
void CalculateKin_AlvarezRuso(GHepRecord *event_rec) const
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
void Configure(const Registry &config)
virtual double Weight(void) const
Definition: GHepRecord.h:124
double x(bool selected=false) const
Definition: Kinematics.cxx:99
Range1D_t YLim(void) const
y limits
void SetKinematics(const double E_l, const double theta_l, const double phi_l, const double theta_pi, const double phi_pi, const Interaction *interaction, Kinematics *kinematics) const
Range1D_t Q2Lim(void) const
Q2 limits.
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
double MaxXSec_AlvarezRuso(const Interaction *in) const
string AsString(void) const
void CalculateKin_ReinSehgal(GHepRecord *event_rec) const
static constexpr double b
Definition: Units.h:78
void ProcessEventRecord(GHepRecord *event_rec) const
Summary information for an interaction.
Definition: Interaction.h:56
double fTMax
upper bound for t = (q - p_pi)^2
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void SetFSLeptonP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:297
string Name(void) const
Definition: AlgId.h:44
static constexpr double A
Definition: Units.h:74
static constexpr double cm2
Definition: Units.h:69
double fQ2Min
lower bound of integration for Q^2 in Berger-Sehgal Model
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
Kinematical phase space.
Definition: KPhaseSpace.h:33
void Sett(double t, bool selected=false)
Definition: Kinematics.cxx:291
static const double kASmallNum
Definition: Controls.h:40
void CalculateKin_BergerSehgal(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.
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
void UpdateWQ2FromXY(const Interaction *in)
Definition: KineUtils.cxx:1290
static RunningThreadInfo * Instance(void)
double ComputeMaxXSec(const Interaction *in) const
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
double max
Definition: Range1.h:53
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
double Energy(const Interaction *in) const
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
double MaxXSec_BergerSehgal(const Interaction *in) const
void UpdateXFromQ2Y(const Interaction *in)
Definition: KineUtils.cxx:1344
bool CheckKinematics(const double E_l, const double theta_l, const double phi_l, const double theta_pi, const double phi_pi, const Interaction *interaction) const
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
double pionMass(const Interaction *in) const
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
void SetHadSystP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:307
virtual double XSec(void) const
Definition: GHepRecord.h:126
double gQ2
Definition: gtestDISSF.cxx:55
static constexpr double fermi
Definition: Units.h:55
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
TF2 * fEnvelope
2-D envelope used for importance sampling
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double min
Definition: Range1.h:52
#define pNOTICE
Definition: Messenger.h:61
void ClearRunningValues(void)
Definition: Kinematics.cxx:347
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
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
double MaxXSec_BergerSehgalFM(const Interaction *in) const
Keep info on the event generation thread currently on charge. This is used so that event generation m...
double fQ2Max
upper bound of integration for Q^2 in Berger-Sehgal Model
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
void throwOnTooManyIterations(unsigned int iters, GHepRecord *evrec) const
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63