GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
QELEventGeneratorSuSA.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  or see $GENIE/LICENSE
6 
7  For the class documentation see the corresponding header file.
8 
9 */
10 //____________________________________________________________________________
11 
12 #include <TMath.h>
13 
16 #include "Framework/Conventions/GBuild.h"
34 
39 
41 
42 using namespace genie;
43 using namespace genie::controls;
44 using namespace genie::constants;
45 using namespace genie::utils;
46 
47 namespace { // anonymous namespace (file only visibility)
48  const double eps = std::numeric_limits<double>::epsilon();
49 }
50 //___________________________________________________________________________
52 KineGeneratorWithCache("genie::QELEventGeneratorSuSA")
53 {
54 
55 }
56 //___________________________________________________________________________
58 KineGeneratorWithCache("genie::QELEventGeneratorSuSA", config)
59 {
60 
61 }
62 //___________________________________________________________________________
64 {
65 
66 }
67 //___________________________________________________________________________
69 {
70  // If we're working with a free nucleon target, then the SuSAv2 calculation
71  // isn't set up to handle it. In these cases, we delegate the work to another
72  // EventRecordVisitorI object (likely QELEventGenerator) configured to run as
73  // a sub-algorithm.
74  if ( !event->Summary()->InitState().Tgt().IsNucleus() ) {
76  }
77 
78  //-- Access cross section algorithm for running thread
80  const EventGeneratorI * evg = rtinfo->RunningThread();
81  fXSecModel = evg->CrossSectionAlg();
82  //event->Print();
83  this -> SelectLeptonKinematics(event);
84  //event->Print();
85  this -> AddTargetNucleusRemnant(event);
86  //event->Print();
87  this -> GenerateNucleon(event);
88  //event->Print();
89 }
90 //___________________________________________________________________________
92 {
93 
94  // -- Event Properties -----------------------------//
95  Interaction * interaction = event->Summary();
96  Kinematics * kinematics = interaction->KinePtr();
97 
98  // Choose the appropriate minimum Q^2 value based on the interaction
99  // mode (this is important for EM interactions since the differential
100  // cross section blows up as Q^2 --> 0)
101  double Q2min = genie::controls::kMinQ2Limit; // CC/NC limit
102  if ( interaction->ProcInfo().IsEM() ) Q2min = genie::utils::kinematics
103  ::electromagnetic::kMinQ2Limit; // EM limit
104 
105  // The SuSA 1p1h model kinematics works in a system where
106  // the whole nuclear target system has no momentum.
107  double Enu = interaction->InitState().ProbeE(kRfLab);
108 
109  int NuPDG = interaction->InitState().ProbePdg();
110  int TgtPDG = interaction->InitState().TgtPdg();
111  // interacton vtx
112  TLorentzVector v4(*event->Probe()->X4());
113  TLorentzVector tempp4(0.,0.,0.,0.);
114 
115  GHepParticle * nucleus = event->TargetNucleus();
116  bool have_nucleus = nucleus != 0;
117 
118  // -- Lepton Kinematic Limits ----------------------------------------- //
119  double Costh = 0.0; // lepton angle
120  double CosthMax = 1.0;
121  double CosthMin = -1.0;
122 
123  double T = 0.0; // lepton kinetic energy
124  double TMax = std::numeric_limits<double>::max();
125  double TMin = 0.0;
126 
127  double Plep = 0.0; // lepton 3 momentum
128  double Elep = 0.0; // lepton energy
129  double LepMass = interaction->FSPrimLepton()->Mass();
130 
131  double Q0 = 0.0; // energy component of q four vector
132  double Q3 = 0.0; // magnitude of transfered 3 momentum
133  double Q2 = 0.0; // properly Q^2 (Q squared) - transfered 4 momentum.
134 
135  // Set lepton KE TMax for for throwing rndm in the accept/reject loop.
136  // We can accidentally set it too high, because the xsec will return zero.
137  // This way if someone reuses this code, they are not tripped up by it.
138  TMax = Enu - LepMass;
139 
140  // Set Tmin for throwing rndm in the accept/reject loop
141  // the hadron tensors we expect will be limited in q3
142  // therefore also the outgoing lepton KE can't be too low or costheta too backward
143  // make the accept/reject loop more efficient by using Min values.
144  if(Enu < fQ3Max){
145  TMin = 0 ;
146  CosthMin = -1 ;
147  } else {
148  TMin = TMath::Sqrt(TMath::Power(LepMass, 2) + TMath::Power((Enu - fQ3Max), 2)) - LepMass;
149  CosthMin = TMath::Sqrt(1 - TMath::Power((fQ3Max / Enu ), 2));
150  }
151 
152 
153  // -- Generate and Test the Kinematics----------------------------------//
154 
155  RandomGen * rnd = RandomGen::Instance();
156  bool accept = false;
157  unsigned int iter = 0;
158  unsigned int maxIter = kRjMaxIterations * 1000;
159 
160  //e-scat xsecs blow up close to theta=0, MC methods won't work so well...
161  // NOTE: SuSAv2 1p1h e-scatting has not been validated yet, use with caution
162  if ( NuPDG == 11 ) maxIter *= 1000;
163 
164  // Get Max XSec:
165  double XSecMax = this->MaxXSec( event );
166 
167  LOG("Kinematics", pDEBUG) << "Max XSec = " << XSecMax;
168 
169  // loop over different (randomly) selected T and Costh
170  while (!accept) {
171  iter++;
172  if(iter > maxIter) {
173  // error if try too many times
174  LOG("QELEvent", pERROR)
175  << "Couldn't select a valid Tmu, CosTheta pair after "
176  << iter << " iterations";
177  event->EventFlags()->SetBitNumber(kKineGenErr, true);
179  exception.SetReason("Couldn't select lepton kinematics");
180  exception.SwitchOnFastForward();
181  throw exception;
182  }
183 
184  // generate random kinetic energy T and Costh
185  T = TMin + (TMax-TMin)*rnd->RndKine().Rndm();
186  Costh = CosthMin + (CosthMax-CosthMin)*rnd->RndKine().Rndm();
187 
188  // Calculate useful values for judging this choice
189  Plep = TMath::Sqrt( T * (T + (2.0 * LepMass))); // ok is sqrt(E2 - m2)
190  Q3 = TMath::Sqrt(Plep*Plep + Enu*Enu - 2.0 * Plep * Enu * Costh);
191 
192  // Calculate what we need to enforce the minimum Q^2 cut
193  Q0 = Enu - (T + LepMass);
194  Q2 = Q3*Q3 - Q0*Q0;
195 
196  // Anti neutrino elastic scattering case - include delta in xsec
197  if (!have_nucleus){
198  genie::utils::mec::Getq0q3FromTlCostl(T, Costh, Enu, LepMass, Q0, Q3);
199  Q3 = sqrt(Q0*Q0+2*kNucleonMass*Q0);
200  genie::utils::mec::GetTlCostlFromq0q3(Q0, Q3, Enu, LepMass, T, Costh);
201  LOG("QELEvent", pDEBUG) << " Anu elastic case. T, Costh: " << T << ", " << Costh ;
202  LOG("QELEvent", pDEBUG) << " Anu elastic case. Q0, Q3: " << Q0 << ", " << Q3 ;
203  }
204 
205  // Don't bother doing hard work if the selected Q3 is greater than Q3Max
206  // or if Q^2 is less than the minimum allowed value
207  if ( Q3 < fQ3Max && Q2 >= Q2min ){
208 
209  kinematics->SetKV(kKVTl, T);
210  kinematics->SetKV(kKVctl, Costh);
211  LOG("QELEvent", pDEBUG) << " T, Costh, Q2: " << T << ", " << Costh << ", " << Q2;
212 
213  double XSec = fXSecModel->XSec(interaction, kPSTlctl);
214 
215  // Some debugging if things go wrong here...
216  if (XSec > XSecMax) {
217  LOG("QELEvent", pDEBUG) << "XSec in cm2 is " << XSec/(units::cm2);
218  LOG("QELEvent", pDEBUG) << "XSec in cm2 /neutron is " << XSec/(units::cm2*pdg::IonPdgCodeToZ(TgtPDG));
219  LOG("QELEvent", pDEBUG) << "XSecMax in cm2 /neutron is " << XSecMax/(units::cm2*pdg::IonPdgCodeToZ(TgtPDG));
220  LOG("QELEvent", pERROR) << " T, Costh, Q2: " << T << ", " << Costh << ", " << Q2;
221  LOG("QELEvent", pERROR) << "XSec is > XSecMax for nucleus " << TgtPDG << " "
222  << XSec << " > " << XSecMax
223  << " don't let this happen.";
224  }
225  // decide whether to accept or reject these kinematics
226  this->AssertXSecLimits( interaction, XSec, XSecMax );
227  accept = XSec > XSecMax*rnd->RndKine().Rndm();
228  LOG("QELEvent", pINFO) << "Xsec, Max, Accept: " << XSec << ", "
229  << XSecMax << ", " << accept;
230  LOG("QELEvent", pDEBUG) << "XSec in cm2 /neutron is " << XSec/(units::cm2*pdg::IonPdgCodeToZ(TgtPDG));
231  LOG("QELEvent", pDEBUG) << "XSecMax in cm2 /neutron is " << XSecMax/(units::cm2*pdg::IonPdgCodeToZ(TgtPDG));
232 
233 
234  }// end if passes q3 test
235  } // end main accept-reject loop
236 
237  // -- finish lepton kinematics
238  // If the code got here, then we accepted some kinematics
239  // and we can proceed to generate the final state.
240 
241  // define coordinate system wrt neutrino: z along neutrino, xy perp
242 
243  // Cos theta gives us z, the rest in xy:
244  double PlepZ = Plep * Costh;
245  double PlepXY = Plep * TMath::Sqrt(1. - TMath::Power(Costh,2));
246 
247  // random rotation about unit vector for phi direction
248  double phi= 2 * kPi * rnd->RndLep().Rndm();
249  // now fill x and y from PlepXY
250  double PlepX = PlepXY * TMath::Cos(phi);
251  double PlepY = PlepXY * TMath::Sin(phi);
252 
253  // Rotate lepton momentum vector from the reference frame (x'y'z') where
254  // {z':(neutrino direction), z'x':(theta plane)} to the LAB
255  TVector3 unit_nudir = event->Probe()->P4()->Vect().Unit();
256  TVector3 p3l(PlepX, PlepY, PlepZ);
257  p3l.RotateUz(unit_nudir);
258 
259  // Lepton 4-momentum in LAB
260  Elep = TMath::Sqrt(LepMass*LepMass + PlepX*PlepX + PlepY*PlepY + PlepZ*PlepZ);
261  TLorentzVector p4l(p3l,Elep);
262 
263  // Figure out the final-state primary lepton PDG code
264  int pdgc = interaction->FSPrimLepton()->PdgCode();
265  int momidx = event->ProbePosition();
266 
267  // -- Store Values ------------------------------------------//
268  // -- Interaction: Q2
269  Q0 = Enu - Elep;
270  Q2 = Q3*Q3 - Q0*Q0;
271  double gy = Q0 / Enu;
272  double gx = kinematics::Q2YtoX(Enu, 2 * kNucleonMass, Q2, gy);
273  double gW = kinematics::XYtoW(Enu, 2 * kNucleonMass, gx, gy);
274 
275  interaction->KinePtr()->SetQ2(Q2, true);
276  interaction->KinePtr()->Sety(gy, true);
277  interaction->KinePtr()->Setx(gx, true);
278  interaction->KinePtr()->SetW(gW, true);
279  interaction->KinePtr()->SetFSLeptonP4(p4l);
280 
281  // -- Lepton
282  event->AddParticle( pdgc, kIStStableFinalState, momidx, -1, -1, -1, p4l, v4);
283 
284  LOG("QELEvent",pDEBUG) << "~~~ LEPTON DONE ~~~";
285 }
286 //___________________________________________________________________________
288 {
289  // add the remnant nuclear target at the GHEP record
290 
291  LOG("QELEvent", pINFO) << "Adding final state nucleus";
292 
293  double Px = 0;
294  double Py = 0;
295  double Pz = 0;
296  double E = 0;
297 
298  GHepParticle * nucleus = event->TargetNucleus();
299  bool have_nucleus = nucleus != 0;
300  if (!have_nucleus) return;
301 
302  int A = nucleus->A();
303  int Z = nucleus->Z();
304 
305  int fd = nucleus->FirstDaughter();
306  int ld = nucleus->LastDaughter();
307 
308  for(int id = fd; id <= ld; id++) {
309 
310  // compute A,Z for final state nucleus & get its PDG code and its mass
311  GHepParticle * particle = event->Particle(id);
312  assert(particle);
313  int pdgc = particle->Pdg();
314  bool is_p = pdg::IsProton (pdgc);
315  bool is_n = pdg::IsNeutron(pdgc);
316 
317  if (is_p) Z--;
318  if (is_p || is_n) A--;
319 
320  Px += particle->Px();
321  Py += particle->Py();
322  Pz += particle->Pz();
323  E += particle->E();
324 
325  }//daughters
326 
327  TParticlePDG * remn = 0;
328  int ipdgc = pdg::IonPdgCode(A, Z);
329  remn = PDGLibrary::Instance()->Find(ipdgc);
330  if(!remn) {
331  LOG("HadronicVtx", pFATAL)
332  << "No particle with [A = " << A << ", Z = " << Z
333  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
334  assert(remn);
335  }
336 
337  double Mi = nucleus->Mass();
338  Px *= -1;
339  Py *= -1;
340  Pz *= -1;
341  E = Mi-E;
342 
343  // Add the nucleus to the event record
344  LOG("QELEvent", pINFO)
345  << "Adding nucleus [A = " << A << ", Z = " << Z
346  << ", pdgc = " << ipdgc << "]";
347 
348  int imom = event->TargetNucleusPosition();
349  event->AddParticle(
350  ipdgc,kIStStableFinalState, imom,-1,-1,-1, Px,Py,Pz,E, 0,0,0,0);
351 
352  LOG("QELEvent", pINFO) << "Done";
353  LOG("QELEvent", pINFO) << *event;
354 }
355 //___________________________________________________________________________
357 {
358  // We need a kinematic limits accept/reject loop here, so generating the
359  // initial hadrons is combined with generating the recoil hadrons...
360 
361  LOG("QELEvent",pDEBUG) << "Generate Nucleon - Start";
362 
363  Interaction * interaction = event->Summary();
364  GHepParticle * neutrino = event->Probe();
365  assert(neutrino);
366  TLorentzVector p4nu(*neutrino->P4());
367 
368  // get final state primary lepton & its 4-momentum
369  GHepParticle * fsl = event->FinalStatePrimaryLepton();
370  assert(fsl);
371  TLorentzVector p4l(*fsl->P4());
372 
373  // calculate 4-momentum transfer from these
374  TLorentzVector Q4 = p4nu - p4l;
375  //Q4.Print();
376 
377  // get the target nucleon and nucleus.
378  // the remnant nucleus is apparently set, except for its momentum.
379  GHepParticle * target_nucleus = event->TargetNucleus();
380  bool have_nucleus = (target_nucleus != 0);
381  //assert(target_nucleus);
382  GHepParticle * initial_nucleon = event->HitNucleon();
383  assert(initial_nucleon);
384  GHepParticle * remnant_nucleus = event->RemnantNucleus();
385  //assert(remnant_nucleus);
386 
387  // instantiate an empty local target nucleus, so I can use existing methods
388  // to get a momentum from the prevailing Fermi-motion distribution
389  int tgtpdg;
390  if(have_nucleus) tgtpdg = target_nucleus->Pdg();
391  else tgtpdg = kPdgTgtFreeP;
392  Target tgt(tgtpdg);
393 
394  // These things need to be saved through to the end of the accept loop.
395  bool accept = false;
396  TVector3 p3i;
397  unsigned int iter = 0;
398 
399  int initial_nucleon_pdg = initial_nucleon->Pdg();
400  int final_nucleon_pdg = interaction->RecoilNucleonPdg();
401 
402  TLorentzVector p4initial_nucleon;
403  TLorentzVector p4final_nucleon;
404  double removalenergy;
405 
406  //remnant kinematic alterations
407  double pxb = 0;
408  double pyb = 0;
409  double pzb = 0;
410 
411  //===========================================================================
412  // Choose nucleons from the prevailing fermi-motion distribution.
413  // Possible to produce kinematically unallowed (Pauli blocked) nucleon.
414  // Find out, and if so choose them again with this accept/reject loop.
415  // Pauli blocking was included in the SuSAv2 tensor tables, so it should not
416  // be allowed to affect the inclusive xsec.
417  while(!accept){
418  iter++;
419  if(iter > kRjMaxIterations) {
420  // error if try too many times
421  LOG("QELEvent", pWARN)
422  << "Couldn't select a valid nucleon after "
423  << iter << " iterations";
424  event->EventFlags()->SetBitNumber(kKineGenErr, true);
426  exception.SetReason("Couldn't select initial hadron kinematics");
427  exception.SwitchOnFastForward();
428  throw exception;
429  }
430 
431  // generate nucleons
432  // tgt is a Target object for local use, just waiting to be filled.
433  // this sets the pdg of each nucleon and its momentum from user chosen nuclear model
434 
435  double hitNucPos = initial_nucleon->X4()->Vect().Mag();
436  tgt.SetHitNucPdg(initial_nucleon_pdg);
437  fNuclModel->GenerateNucleon(tgt,hitNucPos);
438  p3i = fNuclModel->Momentum3();
439 
440  // Default: Calculate the removal energy as in Guille's thesis - this is a simplicification of
441  // a fairly complex aproach employed in SuSAv2, but we expect it to work pretty well.
442  // We should write something about this in the implementation technical paper ...
443  // IMPORTANT CAVEAT: By default we choose to allow the binding energy to depend on the interaction
444  // (as it should), but this means we don't corrolate the chosen Eb with the intial nucleon
445  // momentum. Therefore we can sometimes have initial state nucleons with KE>Eb. This isn't
446  // great, but fot the moment it's what we have (working on improvments).
447 
448  double q3 = Q4.Vect().Mag();
449 
450  if(!have_nucleus){
451  // For elastic case no Fermi momentum and no Eb
452  removalenergy = 0.0;
453  p3i.SetXYZ(0.0,0.0,0.0);
454  }
455  else if(fForceEbFromModel) removalenergy = fNuclModel->RemovalEnergy();
456  else if(fForceFixEb) removalenergy = fEbOR;
457  else{
458  if(q3<0.827){
459  removalenergy = -0.017687 + 0.0564*q3;
460  }
461  else{
462  removalenergy = 0.0289558;
463  }
464  // The above is for Carbon, should shift for different targets
465  removalenergy += (fNuclModel->RemovalEnergy()) - fEbC;
466  if(removalenergy<0.005) removalenergy=0.005;
467  }
468 
469  //removalenergy=0.0;
470  //initial_nucleon->SetRemovalEnergy(removalenergy);
471 
472  LOG("QELEvent",pDEBUG) << "q3 for this event:" << q3;
473  LOG("QELEvent",pDEBUG) << "Removal energy:" << removalenergy;
474 
475  // Now write down the initial nucleon four-vector for this choice
476  double mass = PDGLibrary::Instance()->Find( initial_nucleon_pdg )->Mass();
477  double mass2 = mass*mass;
478  double energy = TMath::Sqrt(p3i.Mag2() + mass2);
479  p4initial_nucleon.SetPxPyPzE(p3i.Px(),p3i.Py(),p3i.Pz(),energy);
480 
481  //One rather unsubtle option for making sure nucleons remain bound.
482  //This will give us bound nucleons with a sensible missing eneergy
483  //but the distribution of Fermi motion will look crazy.
484  // Anything we do is wrong (without semi-inclusive inputs) you just
485  // have to decide what is less wrong!
486  if(fForceBound && (energy-mass>removalenergy)) continue;
487 
488  // cast the removal energy as the energy component of a 4-vector for later.
489  TLorentzVector tLVebind(0., 0., 0., -1.0 * (removalenergy));
490 
491  // Taken from MEC event generator:
492  // calculate recoil nucleon cluster 4-momentum (tLVebind is negative)
493  p4final_nucleon = p4initial_nucleon + Q4 + tLVebind;
494 
495  // Put on shell as in the Aggregator
496  // This is a bit of a horrible approximation but it is hard to think
497  // of anything simple and better without semi-inclusive model predictions.
498  // However, we are working on an improvments.
499  if(have_nucleus){
500  double En = p4final_nucleon.E();
501  double M = PDGLibrary::Instance()->Find( final_nucleon_pdg )->Mass();
502  double pmag_old = p4final_nucleon.Vect().Mag();
503 
504  double pmag_new = TMath::Sqrt(utils::math::NonNegative(En*En-M*M));
505  double scale = pmag_new / pmag_old;
506 
507  double pxn = scale * p4final_nucleon.Px();
508  double pyn = scale * p4final_nucleon.Py();
509  double pzn = scale * p4final_nucleon.Pz();
510 
511  p4final_nucleon.SetPxPyPzE(pxn,pyn,pzn,En);
512  }
513 
514  // Extra momentum xfer to keep nucleon on shell - I'll give this to the remnant
515 
516  pxb = p4nu.Px()-p4l.Px()-p4final_nucleon.Px();
517  pyb = p4nu.Py()-p4l.Py()-p4final_nucleon.Py();
518  pzb = p4nu.Pz()-p4l.Pz()-p4final_nucleon.Pz();
519 
520  LOG("QELEvent",pDEBUG) << "Remnant momentum is: (" << pxb << ", " << pyb << ", " << pzb << ")";
521 
522 
523  // Pauli blocking check:
524  // Test if the resulting four-vector corresponds to a high-enough energy.
525  // Fail the accept if we couldn't put this thing on-shell.
526  // Basically: is energy of the nucleon positive after we subtracting Eb
527  if (p4final_nucleon.E() < PDGLibrary::Instance()->Find(final_nucleon_pdg)->Mass()) {
528  accept = false;
529  LOG("QELEvent",pDEBUG) << "Rejected nucleon, can't be put on-shell";
530  LOG("QELEvent",pDEBUG) << "Nucleon invariant mass:" << p4final_nucleon.M();
531  LOG("QELEvent",pDEBUG) << "Nucleon real mass:" << PDGLibrary::Instance()->Find(final_nucleon_pdg)->Mass();
532  LOG("QELEvent",pDEBUG) << "Nucleon 4 momenutum:";
533  //p4final_nucleon.Print();
534  LOG("QELEvent",pDEBUG) << "Removal energy:" << removalenergy;
535  LOG("QELEvent",pDEBUG) << "Q4 transfer:";
536  //Q4.Print();
537  }
538  else {
539  accept = true;
540  LOG("QELEvent",pDEBUG) << "Nucleon accepted, Q4 is";
541  //Q4.Print();
542  LOG("QELEvent",pDEBUG) << "Initial nucleon mass is" << sqrt((p4initial_nucleon.E()*p4initial_nucleon.E())-(p4initial_nucleon.Vect().Mag()*p4initial_nucleon.Vect().Mag()));
543  LOG("QELEvent",pDEBUG) << "Final nucleon mass is" << sqrt((p4final_nucleon.E()*p4final_nucleon.E())-(p4final_nucleon.Vect().Mag()*p4final_nucleon.Vect().Mag()));
544  LOG("QELEvent",pDEBUG) << "Nucleon invariant mass:" << p4final_nucleon.M();
545  LOG("QELEvent",pDEBUG) << "Nucleon real mass:" << PDGLibrary::Instance()->Find(final_nucleon_pdg)->Mass();
546  LOG("QELEvent",pDEBUG) << "Nucleon 4 momenutum:";
547  //p4final_nucleon.Print();
548  LOG("QELEvent",pDEBUG) << "Removal energy:" << removalenergy;
549  LOG("QELEvent",pDEBUG) << "Q4 transfer:";
550  }
551 
552  } // end accept loop
553 
554  // we got here if we accepted the final state hadron system
555  // so now we need to write everything to ghep
556 
557  // First the initial state nucleon.
558  initial_nucleon->SetMomentum(p4initial_nucleon);
559 
560  // and the remnant nucleus
561  if(have_nucleus){
562  double Mi = PDGLibrary::Instance()->Find(target_nucleus->Pdg() )-> Mass();
563  remnant_nucleus->SetMomentum(pxb,pyb,pzb, Mi - p4initial_nucleon.E() + removalenergy);
564  }
565 
566  // Getting this v4 is a position, i.e. a position within the nucleus (?)
567  // possibly it takes on a non-trivial value only for local Fermi gas
568  // or for sophisticated treatments of intranuclear rescattering.
569  TLorentzVector v4(*neutrino->X4());
570 
571  // Now add the final nucleon
572 
573  interaction->KinePtr()->SetHadSystP4(p4final_nucleon);
574 
576  event->AddParticle(interaction->RecoilNucleonPdg(), ist, event->HitNucleonPosition(),-1,-1,-1, interaction->KinePtr()->HadSystP4(), v4);
577 
578  LOG("QELEvent",pDEBUG) << "Generate Nucleon - End";
579 
580 }
581 //___________________________________________________________________________
583 {
584  Algorithm::Configure(config);
585  this->LoadConfig();
586 }
587 //___________________________________________________________________________
589 {
590  Algorithm::Configure(config);
591  this->LoadConfig();
592 }
593 //___________________________________________________________________________
595 {
596  fNuclModel = 0;
597  RgKey nuclkey = "NuclearModel";
598  fNuclModel = dynamic_cast<const NuclearModelI *> ( this->SubAlg(nuclkey) );
599  assert( fNuclModel );
600 
601  // Sub-algorithm for generating events on free nucleon targets
602  // (not handled by the SuSAv2 calculation)
603  fFreeNucleonEventGenerator = dynamic_cast<const EventRecordVisitorI* >(
604  this->SubAlg("FreeNucleonEventGenerator") );
605  assert( fFreeNucleonEventGenerator );
606 
607  //-- Maximum q3 in input hadron tensors
608  GetParam( "QEL-Q3Max", fQ3Max ) ;
609 
610  //-- Whether to force nucleons to be bound
611  GetParam( "QEL-ForceBound", fForceBound) ;
612 
613  //-- Whether to force Eb to come from the nuclear model
614  GetParam( "QEL-ForceEbFromModel", fForceEbFromModel) ;
615 
616  //-- Whether to force some fixed Eb
617  GetParam( "QEL-ForceFixedEb", fForceFixEb) ;
618  GetParam( "QEL-EbOR", fEbOR) ;
619 
620  //-- Carbon Eb - needed for scaling
621  this->GetParam( "RFG-NucRemovalE@Pdg=1000060120", fEbC);
622 
623  //-- Safety factor for the maximum differential cross section
624  GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor , 5.00 ) ;
625 
626  //-- Minimum energy for which max xsec would be cached, forcing explicit
627  // calculation for lower energies
628  //-- I've set this to an extremely high value to avoid problems with the
629  // SuSAv2 model seen during testing. Everything seems to work OK when
630  // the cache is disabled. - S. Gardiner, 1 July 2020
631  GetParamDef( "Cache-MinEnergy", fEMin, 1000.00 ) ;
632 
633  //-- Maximum allowed fractional cross section deviation from maxim cross
634  // section used in rejection method
635  GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ;
636  assert( fMaxXSecDiffTolerance >= 0. );
637 
638  //-- Generate kinematics uniformly over allowed phase space and compute
639  // an event weight? NOT IMPLEMENTED FOR SUSA YET!
640  GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false ) ;
641 }
642 //____________________________________________________________________________
644  const Interaction * interaction ) const
645 {
646  // Computes the maximum differential cross section in the requested phase
647  // space. This method overloads KineGeneratorWithCache::ComputeMaxXSec
648  // method and the value is cached at a circular cache branch for retrieval
649  // during subsequent event generation.
650  // The computed max differential cross section does not need to be the exact
651  // maximum. The number used in the rejection method will be scaled up by a
652  // safety factor. But it needs to be fast - do not use very small steps.
653 
654  double max_xsec = utils::mec::GetMaxXSecTlctl( *fXSecModel, *interaction );
655 
656  // Apply safety factor, since value retrieved from the cache might
657  // correspond to a slightly different energy
658  max_xsec *= fSafetyFactor;
659 
660  #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
661  SLOG("QELKinematics", pDEBUG) << interaction->AsString();
662  SLOG("QELKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec;
663  SLOG("QELKinematics", pDEBUG) << "Computed using alg = " << *fXSecModel;
664  #endif
665 
666  return max_xsec;
667 }
668 //___________________________________________________________________________
int Z(void) const
bool GetTlCostlFromq0q3(double q0, double q3, double Enu, double ml, double &Tl, double &costl)
Definition: MECUtils.cxx:82
double RemovalEnergy(void) const
Definition: NuclearModelI.h:65
bool fGenerateUniformly
uniform over allowed phase space + event weight?
#define pERROR
Definition: Messenger.h:59
double E(void) const
Get energy.
Definition: GHepParticle.h:91
TRandom3 & RndLep(void) const
rnd number generator used by final state primary lepton generators
Definition: RandomGen.h:62
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
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
double fSafetyFactor
ComputeMaxXSec -&gt; ComputeMaxXSec * fSafetyFactor.
int FirstDaughter(void) const
Definition: GHepParticle.h:68
int RecoilNucleonPdg(void) const
recoil nucleon pdg
virtual void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
#define pFATAL
Definition: Messenger.h:56
bool IsNucleus(void) const
Definition: Target.cxx:272
double GetMaxXSecTlctl(const XSecAlgorithmI &xsec_model, const Interaction &inter, const double tolerance=0.01, const double safety_factor=1.2, const int max_n_layers=100)
Definition: MECUtils.cxx:334
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
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
void SelectLeptonKinematics(GHepRecord *event) const
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:66
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:46
double Mass(Resonance_t res)
resonance mass (GeV)
void SetMomentum(const TLorentzVector &p4)
double Mass(void) const
Mass that corresponds to the PDG code.
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 Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
void AddTargetNucleusRemnant(GHepRecord *event) const
string AsString(void) const
double XYtoW(double Ev, double M, double x, double y)
Definition: KineUtils.cxx:1192
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:75
static const double kMinQ2Limit
Definition: Controls.h:41
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
const double epsilon
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:284
int Pdg(void) const
Definition: GHepParticle.h:63
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
Summary information for an interaction.
Definition: Interaction.h:56
double ComputeMaxXSec(const Interaction *in) const
int LastDaughter(void) const
Definition: GHepParticle.h:69
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
virtual void ProcessEventRecord(GHepRecord *event_rec) const =0
void SetFSLeptonP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:297
static constexpr double A
Definition: Units.h:74
static constexpr double cm2
Definition: Units.h:69
void GenerateNucleon(GHepRecord *event) const
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
int ProbePdg(void) const
Definition: InitialState.h:64
const int kPdgTgtFreeP
Definition: PDGCodes.h:199
bool Getq0q3FromTlCostl(double Tl, double costl, double Enu, double ml, double &q0, double &q3)
Definition: MECUtils.cxx:121
#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
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
bool IsEM(void) const
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
void SetKV(KineVar_t kv, double value)
Definition: Kinematics.cxx:335
static RunningThreadInfo * Instance(void)
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
const EventRecordVisitorI * fFreeNucleonEventGenerator
double Q2YtoX(double Ev, double M, double Q2, double y)
Definition: KineUtils.cxx:1222
string RgKey
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
void ProcessEventRecord(GHepRecord *event) const
void Configure(const Registry &config)
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
void SetHitNucPdg(int pdgc)
Definition: Target.cxx:171
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
int TgtPdg(void) const
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
void SetHadSystP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:307
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double NonNegative(double x)
Definition: MathUtils.cxx:273
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:55
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
int A(void) const
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
virtual bool GenerateNucleon(const Target &) const =0
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...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
Keep info on the event generation thread currently on charge. This is used so that event generation m...
const NuclearModelI * fNuclModel
enum genie::EGHepStatus GHepStatus_t
#define pDEBUG
Definition: Messenger.h:63
double Py(void) const
Get Py.
Definition: GHepParticle.h:89
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345