MECGenerator.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  Steve Dytman <dytman+ \at pitt.edu>
11  Pittsburgh University
12 
13  For the class documentation see the corresponding header file.
14 
15  Important revisions after version 2.0.0 :
16  @ Sep 22, 2008 - CA
17  Skeleton was first added in version 2.5.1
18  @ Nov 24-30, 2010 - CA
19  Major development leading to the first complete version of the generator.
20  @ Nov 20, 2015 - CA, SD
21  Add proper exception handling for failure of phase space decay.
22 */
23 //____________________________________________________________________________
24 
25 #include <TMath.h>
26 
40 
49 
50 using namespace genie;
51 using namespace genie::utils;
52 using namespace genie::constants;
53 using namespace genie::controls;
54 
55 //___________________________________________________________________________
57 EventRecordVisitorI("genie::MECGenerator")
58 {
59 
60 }
61 //___________________________________________________________________________
63 EventRecordVisitorI("genie::MECGenerator", config)
64 {
65 
66 }
67 //___________________________________________________________________________
69 {
70 
71 }
72 //___________________________________________________________________________
74 {
75  //-- Access cross section algorithm for running thread
77  const EventGeneratorI * evg = rtinfo->RunningThread();
78  fXSecModel = evg->CrossSectionAlg();
79  if (fXSecModel->Id().Name() == "genie::EmpiricalMECPXSec2015") {
80  this -> AddTargetRemnant (event); /// shortly, this will be handled by the InitialStateAppender module
81  this -> GenerateFermiMomentum(event);
82  this -> SelectEmpiricalKinematics(event);
83  // TODO: test removing `AddFinalStateLepton` and replacing it with
84  // PrimaryLeptonGenerator::ProcessEventRecord(evrec);
85  this -> AddFinalStateLepton(event);
86  // TODO: maybe put `RecoilNucleonCluster` in a `MECHadronicSystemGenerator` class,
87  // name it `GenerateEmpiricalDiNucleonCluster` or something...
88  this -> RecoilNucleonCluster(event);
89  // TODO: `DecayNucleonCluster` should probably be in `MECHadronicSystemGenerator`,
90  // if we make that...
91  this -> DecayNucleonCluster(event);
92  } else if (fXSecModel->Id().Name() == "genie::NievesSimoVacasMECPXSec2016") {
93  this -> SelectNSVLeptonKinematics(event);
94  this -> AddTargetRemnant(event);
95  this -> GenerateNSVInitialHadrons(event);
96  // Note: this method in `MECTensor/MECTensorGenerator.cxx` appeared to be a straight
97  // copy of an earlier version of the `DecayNucleonCluster` method here - but, watch
98  // for this...
99  this -> DecayNucleonCluster(event);
100  }
101  else {
102  LOG("MECGenerator",pFATAL) <<
103  "ProcessEventRecord >> Cannot calculate kinematics for " <<
104  fXSecModel->Id().Name();
105  }
106 
107 
108 }
109 //___________________________________________________________________________
111 {
112 // Add the remnant nucleus (= initial nucleus - nucleon cluster) in the
113 // event record.
114 
115  GHepParticle * target = event->TargetNucleus();
116  GHepParticle * cluster = event->HitNucleon();
117 
118  int Z = target->Z();
119  int A = target->A();
120 
121  if(cluster->Pdg() == kPdgClusterNN) { A-=2; ; }
122  if(cluster->Pdg() == kPdgClusterNP) { A-=2; Z-=1; }
123  if(cluster->Pdg() == kPdgClusterPP) { A-=2; Z-=2; }
124 
125  int ipdgc = pdg::IonPdgCode(A, Z);
126 
127  const TLorentzVector & p4cluster = *(cluster->P4());
128  const TLorentzVector & p4tgt = *(target ->P4());
129 
130  const TLorentzVector p4 = p4tgt - p4cluster;
131  const TLorentzVector v4(0.,0.,0., 0.);
132 
133  int momidx = event->TargetNucleusPosition();
134 
135  /*
136  if( A == 2 && Z == 2){
137  // residual nucleus was just two protons, not a nucleus we know.
138  // this might not conserve energy...
139  event->AddParticle(kPdgProton,kIStStableFinalState, momidx,-1,-1,-1, p4,v4);
140  // v4,v4 because I'm lazy, give the four momentum to one of the protons, not the other
141  event->AddParticle(kPdgProton,kIStStableFinalState, momidx,-1,-1,-1, v4,v4);
142  } else if ( A == 2 && Z == 0){
143  // residual nucleus was just two neutrons, not a nucleus we know.
144  // this might not conserve energy...
145  event->AddParticle(kPdgNeutron,kIStStableFinalState, momidx,-1,-1,-1, p4,v4);
146  // v4,v4 because I'm lazy, give the four momentum to one of the protons, not the other
147  event->AddParticle(kPdgNeutron,kIStStableFinalState, momidx,-1,-1,-1, v4,v4);
148  } else {
149  // regular nucleus, including deuterium
150  event->AddParticle(ipdgc,kIStStableFinalState, momidx,-1,-1,-1, p4,v4);
151  }
152  */
153 
154  event->AddParticle(ipdgc,kIStStableFinalState, momidx,-1,-1,-1, p4,v4);
155 
156 }
157 //___________________________________________________________________________
159 {
160 // Generate the initial state di-nucleon cluster 4-momentum.
161 // Draw Fermi momenta for each of the two nucleons.
162 // Sum the two Fermi momenta to calculate the di-nucleon momentum.
163 // For simplicity, keep the di-nucleon cluster on the mass shell.
164 //
165  GHepParticle * target_nucleus = event->TargetNucleus();
166  assert(target_nucleus);
167  GHepParticle * nucleon_cluster = event->HitNucleon();
168  assert(nucleon_cluster);
169  GHepParticle * remnant_nucleus = event->RemnantNucleus();
170  assert(remnant_nucleus);
171 
172  // generate a Fermi momentum for each nucleon
173 
174  Target tgt(target_nucleus->Pdg());
175  PDGCodeList pdgv = this->NucleonClusterConstituents(nucleon_cluster->Pdg());
176  assert(pdgv.size()==2);
177  tgt.SetHitNucPdg(pdgv[0]);
179  TVector3 p3a = fNuclModel->Momentum3();
180  tgt.SetHitNucPdg(pdgv[1]);
182  TVector3 p3b = fNuclModel->Momentum3();
183 
184  LOG("FermiMover", pINFO)
185  << "1st nucleon (code = " << pdgv[0] << ") generated momentum: ("
186  << p3a.Px() << ", " << p3a.Py() << ", " << p3a.Pz() << "), "
187  << "|p| = " << p3a.Mag();
188  LOG("FermiMover", pINFO)
189  << "2nd nucleon (code = " << pdgv[1] << ") generated momentum: ("
190  << p3b.Px() << ", " << p3b.Py() << ", " << p3b.Pz() << "), "
191  << "|p| = " << p3b.Mag();
192 
193  // calcute nucleon cluster momentum
194 
195  TVector3 p3 = p3a + p3b;
196 
197  LOG("FermiMover", pINFO)
198  << "di-nucleon cluster momentum: ("
199  << p3.Px() << ", " << p3.Py() << ", " << p3.Pz() << "), "
200  << "|p| = " << p3.Mag();
201 
202  // target (initial) nucleus and nucleon cluster mass
203 
204  double Mi = PDGLibrary::Instance()->Find(target_nucleus->Pdg() )-> Mass(); // initial nucleus mass
205  double M2n = PDGLibrary::Instance()->Find(nucleon_cluster->Pdg())-> Mass(); // nucleon cluster mass
206 
207  // nucleon cluster energy
208 
209  double EN = TMath::Sqrt(p3.Mag2() + M2n*M2n);
210 
211  // set the remnant nucleus and nucleon cluster 4-momenta at GHEP record
212 
213  TLorentzVector p4nclust ( p3.Px(), p3.Py(), p3.Pz(), EN );
214  TLorentzVector p4remnant (-1*p3.Px(), -1*p3.Py(), -1*p3.Pz(), Mi-EN);
215 
216  nucleon_cluster->SetMomentum(p4nclust);
217  remnant_nucleus->SetMomentum(p4remnant);
218 
219  // set the nucleon cluster 4-momentum at the interaction summary
220 
221  event->Summary()->InitStatePtr()->TgtPtr()->SetHitNucP4(p4nclust);
222 }
223 //___________________________________________________________________________
225 {
226 // Select interaction kinematics using the rejection method.
227 //
228 
229  // Access cross section algorithm for running thread
231  const EventGeneratorI * evg = rtinfo->RunningThread();
232  fXSecModel = evg->CrossSectionAlg();
233 
234  Interaction * interaction = event->Summary();
235  double Ev = interaction->InitState().ProbeE(kRfHitNucRest);
236 
237  // **** NOTE / TODO:
238  // **** Hardcode bogus limits for the time-being
239  // **** Should be able to get limits via Interaction::KPhaseSpace
240  double Q2min = 0.01;
241  double Q2max = 8.00;
242  double Wmin = 1.88;
243  double Wmax = 3.00;
244 
245  // Scan phase-space for the maximum differential cross section
246  // at the current neutrino energy
247  const int nq=30;
248  const int nw=20;
249  double dQ2 = (Q2max-Q2min) / (nq-1);
250  double dW = (Wmax-Wmin ) / (nw-1);
251  double xsec_max = 0;
252  for(int iw=0; iw<nw; iw++) {
253  for(int iq=0; iq<nq; iq++) {
254  double Q2 = Q2min + iq*dQ2;
255  double W = Wmin + iw*dW;
256  interaction->KinePtr()->SetQ2(Q2);
257  interaction->KinePtr()->SetW (W);
258  double xsec = fXSecModel->XSec(interaction, kPSWQ2fE);
259  xsec_max = TMath::Max(xsec, xsec_max);
260  }
261  }
262  LOG("MEC", pNOTICE) << "xsec_max (E = " << Ev << " GeV) = " << xsec_max;
263 
264  // Select kinematics
266  unsigned int iter = 0;
267  bool accept = false;
268  while(1) {
269  iter++;
270  if(iter > kRjMaxIterations) {
271  LOG("MEC", pWARN)
272  << "Couldn't select a valid W, Q^2 pair after "
273  << iter << " iterations";
274  event->EventFlags()->SetBitNumber(kKineGenErr, true);
276  exception.SetReason("Couldn't select kinematics");
277  exception.SwitchOnFastForward();
278  throw exception;
279  }
280 
281  // Generate next pair
282  double gQ2 = Q2min + (Q2max-Q2min) * rnd->RndKine().Rndm();
283  double gW = Wmin + (Wmax -Wmin ) * rnd->RndKine().Rndm();
284 
285  // Calculate d2sigma/dQ2dW
286  interaction->KinePtr()->SetQ2(gQ2);
287  interaction->KinePtr()->SetW (gW);
288  double xsec = fXSecModel->XSec(interaction, kPSWQ2fE);
289 
290  // Decide whether to accept the current kinematics
291  double t = xsec_max * rnd->RndKine().Rndm();
292  double J = 1; // jacobean
293  accept = (t < J*xsec);
294 
295  // If the generated kinematics are accepted, finish-up module's job
296  if(accept) {
297  LOG("MEC", pINFO) << "Selected: Q^2 = " << gQ2 << ", W = " << gW;
298  double gx = 0;
299  double gy = 0;
300  // More accurate calculation of the mass of the cluster than 2*Mnucl
301  int nucleon_cluster_pdg = interaction->InitState().Tgt().HitNucPdg();
302  double M2n = PDGLibrary::Instance()->Find(nucleon_cluster_pdg)->Mass();
303  bool is_em = interaction->ProcInfo().IsEM();
304  kinematics::WQ2toXY(Ev,M2n,gW,gQ2,gx,gy);
305 
306  LOG("MEC", pINFO) << "x = " << gx << ", y = " << gy;
307  // lock selected kinematics & clear running values
308  interaction->KinePtr()->SetQ2(gQ2, true);
309  interaction->KinePtr()->SetW (gW, true);
310  interaction->KinePtr()->Setx (gx, true);
311  interaction->KinePtr()->Sety (gy, true);
312  interaction->KinePtr()->ClearRunningValues();
313 
314  return;
315  }//accepted?
316  }//iter
317 }
318 //___________________________________________________________________________
320 {
321 // Add the final-state primary lepton in the event record.
322 // Compute its 4-momentum based on the selected interaction kinematics.
323 //
324  Interaction * interaction = event->Summary();
325 
326  // The boost back to the lab frame was missing, that is included now with the introduction of the beta factor
327  const InitialState & init_state = interaction->InitState();
328  const TLorentzVector & pnuc4 = init_state.Tgt().HitNucP4(); //[@LAB]
329  TVector3 beta = pnuc4.BoostVector();
330 
331  // Boosting the incoming neutrino to the NN-cluster rest frame
332  // Neutrino 4p
333  TLorentzVector * p4v = event->Probe()->GetP4(); // v 4p @ LAB
334  p4v->Boost(-1.*beta); // v 4p @ NN-cluster rest frame
335 
336  // Look-up selected kinematics
337  double Q2 = interaction->Kine().Q2(true);
338  double y = interaction->Kine().y(true);
339 
340  // Auxiliary params
341  double Ev = interaction->InitState().ProbeE(kRfHitNucRest);
342  LOG("MEC", pNOTICE) << "neutrino energy = " << Ev;
343  double ml = interaction->FSPrimLepton()->Mass();
344  double ml2 = TMath::Power(ml,2);
345 
346  // Compute the final state primary lepton energy and momentum components
347  // along and perpendicular the neutrino direction
348  double El = (1-y)*Ev;
349  double plp = El - 0.5*(Q2+ml2)/Ev; // p(//)
350  double plt = TMath::Sqrt(TMath::Max(0.,El*El-plp*plp-ml2)); // p(-|)
351 
352  LOG("MEC", pNOTICE)
353  << "fsl: E = " << El << ", |p//| = " << plp << ", |pT| = " << plt;
354 
355  // Randomize transverse components
357  double phi = 2*kPi * rnd->RndLep().Rndm();
358  double pltx = plt * TMath::Cos(phi);
359  double plty = plt * TMath::Sin(phi);
360 
361  // Take a unit vector along the neutrino direction
362  // WE NEED THE UNIT VECTOR ALONG THE NEUTRINO DIRECTION IN THE NN-CLUSTER REST FRAME, NOT IN THE LAB FRAME
363  TVector3 unit_nudir = p4v->Vect().Unit(); //We use this one, which is in the NN-cluster rest frame
364  // Rotate lepton momentum vector from the reference frame (x'y'z') where
365  // {z':(neutrino direction), z'x':(theta plane)} to the LAB
366  TVector3 p3l(pltx,plty,plp);
367  p3l.RotateUz(unit_nudir);
368 
369  // Lepton 4-momentum in LAB
370  TLorentzVector p4l(p3l,El);
371 
372  // Boost final state primary lepton to the lab frame
373  p4l.Boost(beta); // active Lorentz transform
374 
375  // Figure out the final-state primary lepton PDG code
376  int pdgc = interaction->FSPrimLepton()->PdgCode();
377 
378  // Lepton 4-position (= interacton vtx)
379  TLorentzVector v4(*event->Probe()->X4());
380 
381  int momidx = event->ProbePosition();
382  event->AddParticle(
383  pdgc, kIStStableFinalState, momidx, -1, -1, -1, p4l, v4);
384 }
385 //___________________________________________________________________________
387 {
388  // get di-nucleon cluster & its 4-momentum
389  GHepParticle * nucleon_cluster = event->HitNucleon();
390  assert(nucleon_cluster);
391  TLorentzVector * tmp=nucleon_cluster->GetP4();
392  TLorentzVector p4cluster(*tmp);
393  delete tmp;
394 
395  // get neutrino & its 4-momentum
396  GHepParticle * neutrino = event->Probe();
397  assert(neutrino);
398  TLorentzVector p4v(*neutrino->P4());
399 
400  // get final state primary lepton & its 4-momentum
401  GHepParticle * fsl = event->FinalStatePrimaryLepton();
402  assert(fsl);
403  TLorentzVector p4l(*fsl->P4());
404 
405  // calculate 4-momentum transfer
406  TLorentzVector q = p4v - p4l;
407 
408  // calculate recoil nucleon cluster 4-momentum
409  TLorentzVector p4cluster_recoil = p4cluster + q;
410 
411  // work-out recoil nucleon cluster code
412  LOG("MEC", pINFO) << "Interaction summary";
413  LOG("MEC", pINFO) << *event->Summary();
414  int recoil_nucleon_cluster_pdg = event->Summary()->RecoilNucleonPdg();
415 
416  // vtx position in nucleus coord system
417  TLorentzVector v4(*neutrino->X4());
418 
419  // add to the event record
420  event->AddParticle(
421  recoil_nucleon_cluster_pdg, kIStDecayedState,
422  2, -1, -1, -1, p4cluster_recoil, v4);
423 }
424 //___________________________________________________________________________
426 {
427 // Perform a phase-space decay of the nucleon cluster and add its decay
428 // products in the event record
429 //
430  LOG("MEC", pINFO) << "Decaying nucleon cluster...";
431 
432  // get di-nucleon cluster
433  int nucleon_cluster_id = 5;
434  GHepParticle * nucleon_cluster = event->Particle(nucleon_cluster_id);
435  assert(nucleon_cluster);
436 
437  // get decay products
438  PDGCodeList pdgv = this->NucleonClusterConstituents(nucleon_cluster->Pdg());
439  LOG("MEC", pINFO) << "Decay product IDs: " << pdgv;
440 
441  // Get the decay product masses
442  vector<int>::const_iterator pdg_iter;
443  int i = 0;
444  double * mass = new double[pdgv.size()];
445  double sum = 0;
446  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
447  int pdgc = *pdg_iter;
448  double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
449  mass[i++] = m;
450  sum += m;
451  }
452 
453  LOG("MEC", pINFO)
454  << "Performing a phase space decay to "
455  << pdgv.size() << " particles / total mass = " << sum;
456 
457  TLorentzVector * p4d = nucleon_cluster->GetP4();
458  TLorentzVector * v4d = nucleon_cluster->GetX4();
459 
460  LOG("MEC", pINFO)
461  << "Decaying system p4 = " << utils::print::P4AsString(p4d);
462 
463  // Set the decay
464  bool permitted = fPhaseSpaceGenerator.SetDecay(*p4d, pdgv.size(), mass);
465  if(!permitted) {
466  LOG("MEC", pERROR)
467  << " *** Phase space decay is not permitted \n"
468  << " Total particle mass = " << sum << "\n"
469  << " Decaying system p4 = " << utils::print::P4AsString(p4d);
470  // clean-up
471  delete [] mass;
472  delete p4d;
473  delete v4d;
474  // throw exception
475  event->EventFlags()->SetBitNumber(kHadroSysGenErr, true);
477  exception.SetReason("Decay not permitted kinematically");
478  exception.SwitchOnStepBack();
479  exception.SetReturnStep(0);
480  throw exception;
481  }
482 
483  // Get the maximum weight
484  double wmax = -1;
485  for(int idec=0; idec<200; idec++) {
486  double w = fPhaseSpaceGenerator.Generate();
487  wmax = TMath::Max(wmax,w);
488  }
489  assert(wmax>0);
490  wmax *= 2;
491 
492  LOG("MEC", pNOTICE)
493  << "Max phase space gen. weight = " << wmax;
494 
496  bool accept_decay=false;
497  unsigned int itry=0;
498  while(!accept_decay)
499  {
500  itry++;
501 
503  // report, clean-up and return
504  LOG("MEC", pWARN)
505  << "Couldn't generate an unweighted phase space decay after "
506  << itry << " attempts";
507  // clean up
508  delete [] mass;
509  delete p4d;
510  delete v4d;
511  // throw exception
512  event->EventFlags()->SetBitNumber(kHadroSysGenErr, true);
514  exception.SetReason("Couldn't select decay after N attempts");
515  exception.SwitchOnStepBack();
516  exception.SetReturnStep(0);
517  throw exception;
518  }
519  double w = fPhaseSpaceGenerator.Generate();
520  if(w > wmax) {
521  LOG("MEC", pWARN)
522  << "Decay weight = " << w << " > max decay weight = " << wmax;
523  }
524  double gw = wmax * rnd->RndDec().Rndm();
525  accept_decay = (gw<=w);
526 
527  LOG("MEC", pINFO)
528  << "Decay weight = " << w << " / R = " << gw
529  << " - accepted: " << accept_decay;
530 
531  } //!accept_decay
532 
533  // Insert the decay products in the event record
534  TLorentzVector v4(*v4d);
536  int idp = 0;
537  for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
538  int pdgc = *pdg_iter;
539  TLorentzVector * p4fin = fPhaseSpaceGenerator.GetDecay(idp);
540  event->AddParticle(pdgc, ist, nucleon_cluster_id,-1,-1,-1, *p4fin, v4);
541  idp++;
542  }
543 
544  // Clean-up
545  delete [] mass;
546  delete p4d;
547  delete v4d;
548 }
549 //___________________________________________________________________________
551 {
552  bool allowdup = true;
553  PDGCodeList pdgv(allowdup);
554 
555  if(pdgc == kPdgClusterNN) {
556  pdgv.push_back(kPdgNeutron);
557  pdgv.push_back(kPdgNeutron);
558  }
559  else
560  if(pdgc == kPdgClusterNP) {
561  pdgv.push_back(kPdgNeutron);
562  pdgv.push_back(kPdgProton);
563  }
564  else
565  if(pdgc == kPdgClusterPP) {
566  pdgv.push_back(kPdgProton);
567  pdgv.push_back(kPdgProton);
568  }
569  else
570  {
571  LOG("MEC", pERROR)
572  << "Unknown di-nucleon cluster PDG code (" << pdgc << ")";
573  }
574 
575  return pdgv;
576 }
577 //___________________________________________________________________________
579 {
580  // -- implementation -- //
581  // The IFIC Valencia model can provide three different hadron tensors.
582  // The user probably wants all CC QE-like 2p2h events
583  // But could in principle get the no-delta component if they want (deactivated incode)
584  int FullDeltaNodelta = 1; // 1: full, 2: only delta, 3: zero delta
585 
586  // -- limit the maximum XS for the accept/reject loop -- //
587  //
588  // MaxXSec parameters. This whole calculation could be in it's own function?
589  // these need to lead to a number that is safely large enough, or crash the run.
590  double XSecMaxPar1 = 2.2504;
591  double XSecMaxPar2 = 9.41158;
592 
593  // -- Event Properties -----------------------------//
594  Interaction * interaction = event->Summary();
595  Kinematics * kinematics = interaction->KinePtr();
596 
597  double Enu = interaction->InitState().ProbeE(kRfHitNucRest);
598 
599  int NuPDG = interaction->InitState().ProbePdg();
600  int TgtPDG = interaction->InitState().TgtPdg();
601  // interacton vtx
602  TLorentzVector v4(*event->Probe()->X4());
603  TLorentzVector tempp4(0.,0.,0.,0.);
604 
605  // -- Lepton Kinematic Limits ----------------------------------------- //
606  double Costh = 0.0; // lepton angle
607  double CosthMax = 1.0;
608  double CosthMin = -1.0;
609 
610  double T = 0.0; // lepton kinetic energy
611  double TMax = std::numeric_limits<double>::max();
612  double TMin = 0.0;
613 
614  double Plep = 0.0; // lepton 3 momentum
615  double Elep = 0.0; // lepton energy
616  double LepMass = interaction->FSPrimLepton()->Mass();
617 
618  double Q0 = 0.0; // energy component of q four vector
619  double Q3 = 0.0; // magnitude of transfered 3 momentum
620  double Q2 = 0.0; // properly Q^2 (Q squared) - transfered 4 momentum.
621 
622  // Set lepton KE TMax for for throwing rndm in the accept/reject loop.
623  // We can accidentally set it too high, because the xsec will return zero.
624  // This way if someone reuses this code, they are not tripped up by it.
625  TMax = Enu - LepMass;
626 
627  // Set Tmin for throwing rndm in the accept/reject loop
628  // the hadron tensors we expect will be limited in q3
629  // therefore also the outgoing lepton KE can't be too low or costheta too backward
630  // make the accept/reject loop more efficient by using Min values.
631  if(Enu < fQ3Max){
632  TMin = 0 ;
633  CosthMin = -1 ;
634  } else {
635  TMin = TMath::Sqrt(TMath::Power(LepMass, 2) + TMath::Power((Enu - fQ3Max), 2)) - LepMass;
636  CosthMin = TMath::Sqrt(1 - TMath::Power((fQ3Max / Enu ), 2));
637  }
638 
639  // The accept/reject loop tests a rand against a maxxsec - must scale with A.
640  int NuclearA = 12;
641  double NuclearAfactorXSecMax = 1.0;
642  if (TgtPDG != kPdgTgtC12) {
643  if (TgtPDG > kPdgTgtFreeN && TgtPDG) {
644  NuclearA = pdg::IonPdgCodeToA(TgtPDG);
645  // The QE-like portion scales as A, but the Delta portion increases faster, not simple.
646  // so this gives additional safety factor. Remember, we need a safe max, not precise max.
647  if (NuclearA < 12) NuclearAfactorXSecMax *= NuclearA / 12.0;
648  else NuclearAfactorXSecMax *= TMath::Power(NuclearA/12.0, 1.4);
649  }
650  else {
651  LOG("MEC", pERROR) << "Trying to scale XSecMax for larger nuclei, but "
652  << TgtPDG << " isn't a nucleus?";
653  assert(false);
654  }
655  }
656 
657  // -- Generate and Test the Kinematics----------------------------------//
658 
660  bool accept = false;
661  unsigned int iter = 0;
662 
663  // loop over different (randomly) selected T and Costh
664  while (!accept) {
665  iter++;
666  if(iter > kRjMaxIterations) {
667  // error if try too many times
668  LOG("MEC", pWARN)
669  << "Couldn't select a valid Tmu, CosTheta pair after "
670  << iter << " iterations";
671  event->EventFlags()->SetBitNumber(kKineGenErr, true);
673  exception.SetReason("Couldn't select lepton kinematics");
674  exception.SwitchOnFastForward();
675  throw exception;
676  }
677 
678  // generate random kinetic energy T and Costh
679  T = TMin + (TMax-TMin)*rnd->RndKine().Rndm();
680  Costh = CosthMin + (CosthMax-CosthMin)*rnd->RndKine().Rndm();
681 
682  // Calculate useful values for judging this choice
683  Plep = TMath::Sqrt( T * (T + (2.0 * LepMass))); // ok is sqrt(E2 - m2)
684  Q3 = TMath::Sqrt(Plep*Plep + Enu*Enu - 2.0 * Plep * Enu * Costh);
685 
686  // Don't bother doing hard work if the selected Q3 is greater than Q3Max
687  if (Q3 < fQ3Max){
688 
689  kinematics->SetKV(kKVTl, T);
690  kinematics->SetKV(kKVctl, Costh);
691 
692  // decide whether to accept or reject these kinematics
693  // AND set the chosen two-nucleon system
694 
695  // to save time, use a pre-calculated max cross-section XSecMax
696  // it doesn't matter what it is, as long as it is big enough.
697  // RIK asks can XSecMax can be pushed to the q0q3 part of the calculation
698  // where the XS doesn't depend much on Enu.
699  // instead, this implementation uses a rough dependence on log10(Enu).
700  // starting around 0.5 GeV, the log10(max) is linear vs. log10(Enu)
701  // 1.10*TMath::Power(10.0, 2.2504 * TMath::Log10(Enu) - 9.41158)
702 
703  if (FullDeltaNodelta == 1){
704  // this block for the user who wants all CC QE-like 2p2h events
705 
706  // extract xsecmax from the spline making process for C12 and other nuclei.
707  // plot Log10(E) on horizontal and Log10(xsecmax) vertical
708  // and fit a line. Use that plus 1.35 safety factors to limit the accept/reject loop.
709  double XSecMax = 1.35 * TMath::Power(10.0, XSecMaxPar1 * TMath::Log10(Enu) - XSecMaxPar2);
710  if (NuclearA > 12) XSecMax *= NuclearAfactorXSecMax; // Scale it by A, precomputed above.
711 
712  LOG("MEC", pDEBUG) << " T, Costh: " << T << ", " << Costh ;
713 
714 
715  // We need four different cross sections. Right now, pursue the
716  // inelegant method of calling XSec four times - there is
717  // definitely some runtime inefficiency here, but it is not awful
718 
719  // first, get delta-less all
720  if (NuPDG > 0) {
721  interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterNN);
722  }
723  else {
724  interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterPP);
725  }
726  double XSec = fXSecModel->XSec(interaction, kPSTlctl);
727  // now get all with delta
728  interaction->ExclTagPtr()->SetResonance(genie::kP33_1232);
729  double XSecDelta = fXSecModel->XSec(interaction, kPSTlctl);
730  // get PN with delta
731  interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterNP);
732  double XSecDeltaPN = fXSecModel->XSec(interaction, kPSTlctl);
733  // now get delta-less PN
735  double XSecPN = fXSecModel->XSec(interaction, kPSTlctl);
736 
737  if (XSec > XSecMax) {
738  LOG("MEC", pERROR) << "XSec is > XSecMax for nucleus " << TgtPDG << " "
739  << XSec << " > " << XSecMax
740  << " don't let this happen.";
741  }
742  assert(XSec <= XSecMax);
743  accept = XSec > XSecMax*rnd->RndKine().Rndm();
744  LOG("MEC", pINFO) << "Xsec, Max, Accept: " << XSec << ", "
745  << XSecMax << ", " << accept;
746 
747  if(accept){
748  // If it passes the All cross section we still need to do two things:
749  // * Was the initial state pn or not?
750  // * Do we assign the reaction to have had a Delta on the inside?
751 
752  // PDD means from the part of the XSec with an internal Delta line
753  // that (at the diagram level) did not produce a pion in the final state.
754 
755  bool isPDD = false;
756 
757  // Find out if we should use a pn initial state
758  double myrand = rnd->RndKine().Rndm();
759  double pnFraction = XSecPN / XSec;
760  LOG("MEC", pDEBUG) << "Test for pn: xsec_pn = " << XSecPN
761  << "; xsec = " << XSec
762  << "; pn_fraction = " << pnFraction
763  << "; random number val = " << myrand;
764 
765  if (myrand <= pnFraction) {
766  // yes it is, add a PN initial state to event record
767  event->AddParticle(kPdgClusterNP, kIStNucleonTarget,
768  1, -1, -1, -1, tempp4, v4);
769  interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterNP);
770 
771  // Its a pn, so test for Delta by comparing DeltaPN/PN
772  if (rnd->RndKine().Rndm() <= XSecDeltaPN / XSecPN) {
773  isPDD = true;
774  }
775  }
776  else {
777  // no it is not a PN, add either NN or PP initial state to event record.
778  if (NuPDG > 0) {
779  event->AddParticle(kPdgClusterNN, kIStNucleonTarget,
780  1, -1, -1, -1, tempp4, v4);
781  interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterNN);
782  }
783  else {
784  event->AddParticle(kPdgClusterPP, kIStNucleonTarget,
785  1, -1, -1, -1, tempp4, v4);
786  interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(kPdgClusterPP);
787  }
788  // its not pn, so test for Delta (XSecDelta-XSecDeltaPN)/(XSec-XSecPN)
789  // right, both numerator and denominator are total not pn.
790  if (rnd->RndKine().Rndm() <=
791  (XSecDelta - XSecDeltaPN) / (XSec - XSecPN)) {
792  isPDD = true;
793  }
794  }
795 
796  // now test whether we tagged this as a pion event
797  // and assign that fact to the Exclusive State tag
798  // later, we can query const XclsTag & xcls = interaction->ExclTag()
799  if (isPDD){
800  interaction->ExclTagPtr()->SetResonance(genie::kP33_1232);
801  }
802 
803 
804  } // end if accept
805  } // end if delta == 1
806 
807  /* One can make simpler versions of the above for the
808  FullDeltaNodelta == 2 (only delta)
809  or
810  FullDeltaNodelta == 3 (set Delta FF = 1, lose interference effect).
811  but I (Rik) don't see what the use-case is for these, genratorly speaking.
812  */
813 
814  }// end if passes q3 test
815  } // end while
816 
817  // -- finish lepton kinematics
818  // If the code got here, then we accepted some kinematics
819  // and we can proceed to generate the final state.
820 
821  // define coordinate system wrt neutrino: z along neutrino, xy perp
822 
823  // Cos theta gives us z, the rest in xy:
824  double PlepZ = Plep * Costh;
825  double PlepXY = Plep * TMath::Sqrt(1. - TMath::Power(Costh,2));
826 
827  // random rotation about unit vector for phi direction
828  double phi= 2 * kPi * rnd->RndLep().Rndm();
829  // now fill x and y from PlepXY
830  double PlepX = PlepXY * TMath::Cos(phi);
831  double PlepY = PlepXY * TMath::Sin(phi);
832 
833  // Rotate lepton momentum vector from the reference frame (x'y'z') where
834  // {z':(neutrino direction), z'x':(theta plane)} to the LAB
835  TVector3 unit_nudir = event->Probe()->P4()->Vect().Unit();
836  TVector3 p3l(PlepX, PlepY, PlepZ);
837  p3l.RotateUz(unit_nudir);
838 
839  // Lepton 4-momentum in LAB
840  Elep = TMath::Sqrt(LepMass*LepMass + PlepX*PlepX + PlepY*PlepY + PlepZ*PlepZ);
841  TLorentzVector p4l(p3l,Elep);
842 
843  // Figure out the final-state primary lepton PDG code
844  int pdgc = interaction->FSPrimLepton()->PdgCode();
845  int momidx = event->ProbePosition();
846 
847  // -- Store Values ------------------------------------------//
848  // -- Interaction: Q2
849  Q0 = Enu - Elep;
850  Q2 = Q3*Q3 - Q0*Q0;
851  double gy = Q0 / Enu;
852  double gx = kinematics::Q2YtoX(Enu, 2 * kNucleonMass, Q2, gy);
853  double gW = kinematics::XYtoW(Enu, 2 * kNucleonMass, gx, gy);
854 
855  interaction->KinePtr()->SetQ2(Q2, true);
856  interaction->KinePtr()->Sety(gy, true);
857  interaction->KinePtr()->Setx(gx, true);
858  interaction->KinePtr()->SetW(gW, true);
859  interaction->KinePtr()->SetFSLeptonP4(p4l);
860  // in later methods
861  // will also set the four-momentum and W^2 of the hadron system.
862 
863  // -- Lepton
864  event->AddParticle( pdgc, kIStStableFinalState, momidx, -1, -1, -1, p4l, v4);
865 
866  LOG("MEC",pDEBUG) << "~~~ LEPTON DONE ~~~";
867 }
868 //___________________________________________________________________________
870 {
871  // We need a kinematic limits accept/reject loop here, so generating the
872  // initial hadrons is combined with generating the recoil hadrons...
873 
874  LOG("MEC",pDEBUG) << "Generate Initial Hadrons - Start";
875 
876  Interaction * interaction = event->Summary();
877  GHepParticle * neutrino = event->Probe();
878  assert(neutrino);
879  TLorentzVector p4nu(*neutrino->P4());
880 
881  // get final state primary lepton & its 4-momentum
882  GHepParticle * fsl = event->FinalStatePrimaryLepton();
883  assert(fsl);
884  TLorentzVector p4l(*fsl->P4());
885 
886  // calculate 4-momentum transfer from these
887  TLorentzVector Q4 = p4nu - p4l;
888 
889  // get the target two-nucleon cluster and nucleus.
890  // the remnant nucleus is apparently set, except for its momentum.
891  GHepParticle * target_nucleus = event->TargetNucleus();
892  assert(target_nucleus);
893  GHepParticle * initial_nucleon_cluster = event->HitNucleon();
894  assert(initial_nucleon_cluster);
895  GHepParticle * remnant_nucleus = event->RemnantNucleus();
896  assert(remnant_nucleus);
897 
898  // -- make a two-nucleon system, then give them some momenta.
899 
900  // instantiate an empty local target nucleus, so I can use existing methods
901  // to get a momentum from the prevailing Fermi-motion distribution
902  Target tgt(target_nucleus->Pdg());
903 
904  // NucleonClusterConstituents is an implementation within this class, called with this
905  // It using the nucleon cluster from the earlier tests for a pn state,
906  // the method returns a vector of pdgs, which hopefully will be of size two.
907 
908  PDGCodeList pdgv = this->NucleonClusterConstituents(initial_nucleon_cluster->Pdg());
909  assert(pdgv.size()==2);
910 
911  // These things need to be saved through to the end of the accept loop.
912  bool accept = false;
913  TVector3 p31i;
914  TVector3 p32i;
915  unsigned int iter = 0;
916 
917  int initial_nucleon_cluster_pdg = initial_nucleon_cluster->Pdg();
918  int final_nucleon_cluster_pdg = 0;
919  if (neutrino->Pdg() > 0) {
920  if (initial_nucleon_cluster->Pdg() == kPdgClusterNP) {
921  final_nucleon_cluster_pdg = kPdgClusterPP;
922  }
923  else if (initial_nucleon_cluster->Pdg() == kPdgClusterNN) {
924  final_nucleon_cluster_pdg = kPdgClusterNP;
925  }
926  else {
927  LOG("MEC", pERROR) << "Wrong pdg for a CC neutrino MEC interaction"
928  << initial_nucleon_cluster->Pdg();
929  }
930  }
931  else if (neutrino->Pdg() < 0) {
932  if (initial_nucleon_cluster->Pdg() == kPdgClusterNP) {
933  final_nucleon_cluster_pdg = kPdgClusterNN;
934  }
935  else if (initial_nucleon_cluster->Pdg() == kPdgClusterPP) {
936  final_nucleon_cluster_pdg = kPdgClusterNP;
937  }
938  else {
939  LOG("MEC", pERROR) << "Wrong pdg for a CC anti-neutrino MEC interaction"
940  << initial_nucleon_cluster->Pdg();
941  }
942  }
943 
944  TLorentzVector p4initial_cluster;
945  TLorentzVector p4final_cluster;
946  TLorentzVector p4remnant_nucleus;
947  double removalenergy1;
948  double removalenergy2;
949 
950  //===========================================================================
951  // Choose two nucleons from the prevailing fermi-motion distribution.
952  // Some produce kinematically unallowed combinations initial cluster and Q2
953  // Find out, and if so choose them again with this accept/reject loop.
954  // Some kinematics are especially tough
955  while(!accept){
956  iter++;
957  if(iter > kRjMaxIterations*1000) {
958  // error if try too many times
959  LOG("MEC", pWARN)
960  << "Couldn't select a valid W, Q^2 pair after "
961  << iter << " iterations";
962  event->EventFlags()->SetBitNumber(kKineGenErr, true);
964  exception.SetReason("Couldn't select initial hadron kinematics");
965  exception.SwitchOnFastForward();
966  throw exception;
967  }
968 
969  // generate nucleons
970  // tgt is a Target object for local use, just waiting to be filled.
971  // this sets the pdg of each nucleon and its momentum from a Fermi-gas
972  // Nieves et al. would use a local Fermi gas here, not this, but ok.
973  // so momentum from global Fermi gas, local Fermi gas, or spectral function
974  // and removal energy ~0.025 GeV, correlated with density, or from SF distribution
975  tgt.SetHitNucPdg(pdgv[0]);
977  p31i = fNuclModel->Momentum3();
978  removalenergy1 = fNuclModel->RemovalEnergy();
979  tgt.SetHitNucPdg(pdgv[1]);
981  p32i = fNuclModel->Momentum3();
982  removalenergy2 = fNuclModel->RemovalEnergy();
983 
984  // not sure -- could give option to use Nieves q-value here.
985 
986  // Now write down the initial cluster four-vector for this choice
987  TVector3 p3i = p31i + p32i;
988  double mass2 = PDGLibrary::Instance()->Find( initial_nucleon_cluster_pdg )->Mass();
989  mass2 *= mass2;
990  double energy = TMath::Sqrt(p3i.Mag2() + mass2);
991  p4initial_cluster.SetPxPyPzE(p3i.Px(),p3i.Py(),p3i.Pz(),energy);
992 
993  // cast the removal energy as the energy component of a 4-vector for later.
994  TLorentzVector tLVebind(0., 0., 0., -1.0 * (removalenergy1 + removalenergy2));
995 
996  // RIK: You might ask why is this the right place to subtract ebind?
997  // It is okay. Physically, I'm subtracting it from q. The energy
998  // transfer to the nucleon is 50 MeV less. the energy cost to put this
999  // cluster on-shell. What Jan says he does in PRC.86.015504 is this:
1000  // The nucleons are assumed to be in a potential well
1001  // V = Efermi + 8 MeV.
1002  // The Fermi energy is subtracted from each initial-state nucleon
1003  // (I guess he does it dynamically based on Ef = p2/2M or so which
1004  // is what we are doing above, on average). Then after the reaction,
1005  // another 8 MeV is subtracted at that point; a small adjustment to the
1006  // momentum is needed to keep the nucleon on shell.
1007 
1008  // calculate recoil nucleon cluster 4-momentum (tLVebind is negative)
1009  p4final_cluster = p4initial_cluster + Q4 + tLVebind;
1010 
1011  // Test if the resulting four-vector corresponds to a high-enough invariant mass.
1012  // Fail the accept if we couldn't put this thing on-shell.
1013  if (p4final_cluster.M() <
1014  PDGLibrary::Instance()->Find(final_nucleon_cluster_pdg )->Mass()) {
1015  accept = false;
1016  } else {
1017  accept = true;
1018  }
1019 
1020  } // end accept loop
1021 
1022  // we got here if we accepted the final state two-nucleon system
1023  // so now we need to write everything to ghep
1024 
1025  // First the initial state nucleons.
1026  initial_nucleon_cluster->SetMomentum(p4initial_cluster);
1027 
1028  // and the remnant nucleus
1029  double Mi = PDGLibrary::Instance()->Find(target_nucleus->Pdg() )-> Mass();
1030  remnant_nucleus->SetMomentum(-1.0*p4initial_cluster.Px(),
1031  -1.0*p4initial_cluster.Py(),
1032  -1.0*p4initial_cluster.Pz(),
1033  Mi - p4initial_cluster.E() + removalenergy1 + removalenergy2);
1034 
1035  // Now the final nucleon cluster.
1036 
1037  // Getting this v4 is a position, i.e. a position within the nucleus (?)
1038  // possibly it takes on a non-trivial value only for local Fermi gas
1039  // or for sophisticated treatments of intranuclear rescattering.
1040  TLorentzVector v4(*neutrino->X4());
1041 
1042  // Now write the (undecayed) final two-nucleon system
1043  GHepParticle p1(final_nucleon_cluster_pdg, kIStDecayedState,
1044  2, -1, -1, -1, p4final_cluster, v4);
1045 
1046  //p1.SetRemovalEnergy(removalenergy1 + removalenergy2);
1047  // The "bound particle" concept applies only to p or n.
1048  // Instead, add this directly to the remnant nucleon a few lines above.
1049 
1050  // actually, this is not an status1 particle, so it is not picked up
1051  // by the aggregator. anyway, the aggregator does not run until the very end.
1052 
1053  event->AddParticle(p1);
1054 
1055  interaction->KinePtr()->SetHadSystP4(p4final_cluster);
1056 }
1057 //___________________________________________________________________________
1059 {
1060  Algorithm::Configure(config);
1061  this->LoadConfig();
1062 }
1063 //___________________________________________________________________________
1065 {
1066  Algorithm::Configure(config);
1067  this->LoadConfig();
1068 }
1069 //___________________________________________________________________________
1071 {
1072  fNuclModel = 0;
1073  RgKey nuclkey = "NuclearModel";
1074  fNuclModel = dynamic_cast<const NuclearModelI *> (this->SubAlg(nuclkey));
1075  assert(fNuclModel);
1076 
1077  GetParam( "NSV-Q3Max", fQ3Max ) ;
1078 }
1079 //___________________________________________________________________________
1080 
const double kPi
int Z(void) const
void ProcessEventRecord(GHepRecord *event) const
Basic constants.
TLorentzVector * GetX4(void) const
const XML_Char * target
Definition: expat.h:268
void SelectNSVLeptonKinematics(GHepRecord *event) const
double RemovalEnergy(void) const
Definition: NuclearModelI.h:57
const NuclearModelI * fNuclModel
Definition: MECGenerator.h:66
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:141
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
TRandom3 & RndLep(void) const
rnd number generator used by final state primary lepton generators
Definition: RandomGen.h:63
static const double kNucleonMass
Definition: Constants.h:78
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
int HitNucPdg(void) const
Definition: Target.cxx:321
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:265
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
const int kPdgClusterNP
Definition: PDGCodes.h:192
enum genie::EGHepStatus GHepStatus_t
void Configure(const Registry &config)
#define pFATAL
Definition: Messenger.h:57
Defines the EventGeneratorI interface.
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:61
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:34
static const unsigned int kMaxUnweightDecayIterations
Definition: Controls.h:61
double Wmax[kNWBins]
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:42
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
double Mass(Resonance_t res)
resonance mass (GeV)
void SetMomentum(const TLorentzVector &p4)
Definition: config.py:1
Float_t tmp
Definition: plot.C:36
Double_t beta
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
void SetResonance(Resonance_t res)
Definition: XclsTag.cxx:123
double XYtoW(double Ev, double M, double x, double y)
Definition: KineUtils.cxx:1087
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:67
A list of PDG codes.
Definition: PDGCodeList.h:33
const int kPdgClusterNN
Definition: PDGCodes.h:191
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:321
double y(bool selected=false) const
Definition: Kinematics.cxx:122
Float_t Z
Definition: plot.C:38
int Pdg(void) const
Definition: GHepParticle.h:64
void GenerateNSVInitialHadrons(GHepRecord *event) const
Summary information for an interaction.
Definition: Interaction.h:56
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:92
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:97
void SetFSLeptonP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:307
string Name(void) const
Definition: AlgId.h:45
void RecoilNucleonCluster(GHepRecord *event) const
void DecayNucleonCluster(GHepRecord *event) const
TLorentzVector * GetP4(void) const
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
Definition: KineUtils.cxx:1046
const Kinematics & Kine(void) const
Definition: Interaction.h:71
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
int ProbePdg(void) const
Definition: InitialState.h:65
double energy
Definition: plottest35.C:25
const int kPdgTgtFreeN
Definition: PDGCodes.h:177
#define pINFO
Definition: Messenger.h:63
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:51
Misc GENIE control constants.
PDGCodeList NucleonClusterConstituents(int pdgc) const
Double_t xsec[nknots]
Definition: testXsec.C:47
void SelectEmpiricalKinematics(GHepRecord *event) const
#define pWARN
Definition: Messenger.h:61
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
bool IsEM(void) const
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:241
XclsTag * ExclTagPtr(void) const
Definition: Interaction.h:77
void SetKV(KineVar_t kv, double value)
Definition: Kinematics.cxx:345
static RunningThreadInfo * Instance(void)
static const double A
Definition: Units.h:82
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:289
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
double Q2YtoX(double Ev, double M, double Q2, double y)
Definition: KineUtils.cxx:1117
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
string RgKey
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
void SetHitNucPdg(int pdgc)
Definition: Target.cxx:188
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:253
int TgtPdg(void) const
Target * TgtPtr(void) const
Definition: InitialState.h:68
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:69
const XSecAlgorithmI * fXSecModel
Definition: MECGenerator.h:64
assert(nhit_max >=nhit_nbins)
void SetHadSystP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:317
const int kPdgTgtC12
Definition: PDGCodes.h:179
double gQ2
Definition: gtestDISSF.cxx:55
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
virtual bool GenerateNucleon(const Target &) const =0
const InitialState & InitState(void) const
Definition: Interaction.h:69
void AddTargetRemnant(GHepRecord *event) const
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
void GenerateFermiMomentum(GHepRecord *event) const
double Wmin[kNWBins]
const int kPdgProton
Definition: PDGCodes.h:65
double T
Definition: Xdiff_gwt.C:5
int A(void) const
double Q2(bool selected=false) const
Definition: Kinematics.cxx:135
#define pNOTICE
Definition: Messenger.h:62
void ClearRunningValues(void)
Definition: Kinematics.cxx:357
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:67
static const unsigned int kRjMaxIterations
Definition: Controls.h:27
Double_t sum
Definition: plot.C:31
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
const EventGeneratorI * RunningThread(void)
void kinematics()
Definition: kinematics.C:10
Float_t w
Definition: plot.C:20
#define W(x)
double ProbeE(RefFrame_t rf) const
const int kPdgNeutron
Definition: PDGCodes.h:67
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
TRandom3 & RndDec(void) const
rnd number generator used by decay models
Definition: RandomGen.h:57
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
Keep info on the event generation thread currently on charge. This is used so that event generation m...
void AddFinalStateLepton(GHepRecord *event) const
Root of GENIE utility namespaces.
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
TGenPhaseSpace fPhaseSpaceGenerator
Definition: MECGenerator.h:65
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
const int kPdgClusterPP
Definition: PDGCodes.h:193
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353