GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CEvNSEventGenerator.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 <TMath.h>
14 #include <Math/IFunction.h>
15 #include <Math/GSLMinimizer1D.h>
16 #include <Math/BrentMinimizer1D.h>
17 
37 
38 using namespace genie;
39 using namespace genie::constants;
40 using namespace genie::controls;
41 using namespace genie::utils;
42 
43 //___________________________________________________________________________
45 EventRecordVisitorI("genie::CEvNSEventGenerator")
46 {
47 
48 }
49 //___________________________________________________________________________
51 EventRecordVisitorI("genie::COHElKinematicsGenerator", config)
52 {
53 
54 }
55 //___________________________________________________________________________
57 {
58 
59 }
60 //___________________________________________________________________________
62 {
63  this -> GenerateKinematics (event);
64  this -> AddFinalStateNeutrino (event);
65  this -> AddRecoilNucleus (event);
66 }
67 //___________________________________________________________________________
69 {
70  Interaction * interaction = event->Summary();
71  interaction->SetBit(kISkipProcessChk);
72  interaction->SetBit(kISkipKinematicChk);
73 
74  // Get the random number generators
76 
77  // Get the kinematical limits
78  const KPhaseSpace & kps = interaction->PhaseSpace();
79  Range1D_t Q2 = kps.Q2Lim();
80  assert(Q2.min > 0. && Q2.min < Q2.max);
81  const InitialState & init_state = interaction -> InitState();
82  double E = init_state.ProbeE(kRfLab);
83  const double Q2min = Q2.min;
84  const double Q2max = Q2.max;
85  const double dQ2 = Q2max - Q2min;
86 
87  // Access cross section algorithm for running thread
89  const EventGeneratorI * evg = rtinfo->RunningThread();
90  fXSecModel = evg->CrossSectionAlg();
91 
92  double gQ2 = -1; // generated Q2
93  double gxsec = -1; // dsig/dQ2 at generated Q2
94 
95  // Generate kinematics
96  if(fGenerateUniformly) {
97  LOG("CEvNS", pFATAL)
98  << "Option to generate kinematics uniformly not supported";
99  exit(1);
100 /*
101  gQ2 = Q2min + dQ2 * rnd->RndKine().Rndm();
102  LOG("CEvNS", pINFO) << "Trying: Q2 = " << gQ2;
103  interaction->KinePtr()->SetQ2(gQ2);
104 
105  // Computing cross section for the current kinematics
106  gxsec = fXSecModel->XSec(interaction, kPSQ2fE);
107  if(gxsec<=0){
108  LOG("CEvNS", pWARN)
109  << "Non-positive x-section for selected Q2 = " << gQ2 << "GeV^2";
110  }
111 
112  double weight = 1; // to implement if fGenerateUniformly option is enabled
113  LOG("CEvNS", pDEBUG) << "Kinematics wght = "<< weight;
114  // apply computed weight to the current event weight
115  weight *= event->Weight();
116  LOG("CEvNS", pDEBUG) << "Current event wght = " << weight;
117  event->SetWeight(weight);
118 */
119  }
120  else {
121  // For the subsequent kinematic selection with the rejection method:
122  // Calculate the max differential cross section.
123  // Always at Q^2 = 0 for energies and model tested,
124  // but go ahead and do the calculation nevertheless.
125  ROOT::Math::IBaseFunctionOneDim * xsec_func =
126  new utils::gsl::dXSec_dQ2_E(fXSecModel, interaction,-1.);
127  ROOT::Math::BrentMinimizer1D minimizer;
128  minimizer.SetFunction(*xsec_func,Q2min,Q2max);
129  minimizer.Minimize(1000,1,1E-5);
130  double Q2_for_xsec_max = minimizer.XMinimum();
131  interaction->KinePtr()->SetQ2(Q2_for_xsec_max);
132  double xsec_max = fXSecModel->XSec(interaction, kPSQ2fE);
133  delete xsec_func;
134  LOG("CEvNS", pNOTICE)
135  << "Maximizing dsig(Q2;E = " << E << "GeV)/dQ2 gave a value of "
136  << xsec_max/(units::cm2) << " cm2/GeV^2 at Q2 = "
137  << Q2_for_xsec_max << " GeV^2";
138 
139  // Try to select a valid Q2
140  unsigned int iter = 0;
141  while(1) {
142  iter++;
143  if(iter > kRjMaxIterations) {
144  LOG("CEvNS", pWARN)
145  << "*** Could not select a valid Q2 after " << iter << " iterations";
146  event->EventFlags()->SetBitNumber(kKineGenErr, true);
148  exception.SetReason("Couldn't select kinematics");
149  exception.SwitchOnFastForward();
150  throw exception;
151  } // max iterations
152 
153  gQ2 = Q2min + dQ2 * rnd->RndKine().Rndm();
154  LOG("CEvNS", pINFO) << "Trying: Q2 = " << gQ2;
155  interaction->KinePtr()->SetQ2(gQ2);
156 
157  // Computing cross section for the current kinematics
158  gxsec = fXSecModel->XSec(interaction, kPSQ2fE);
159 
160  if(gxsec > xsec_max) {
161  double frac = TMath::Abs(gxsec-xsec_max)/xsec_max;
162  if(frac > fMaxXSecDiffTolerance) {
163  LOG("CEvNS", pWARN)
164  << "Current computed cross-section (" << gxsec/(units::cm2)
165  << " cm2/GeV^2) exceeds the maximum cross-section ("
166  << xsec_max/(units::cm2) << " beyond the specified tolerance";
167  }
168  }
169 
170  // Decide whether to accept the current kinematic point
171  double t = fSafetyFactor * xsec_max * rnd->RndKine().Rndm();
172  //this->AssertXSecLimits(interaction, gxsec, xsec_max);
173  LOG("CEvNS", pINFO)
174  << "dxsec/dQ2 = " << gxsec/(units::cm2) << " cm2/GeV^2"
175  << "J = 1, rnd = " << t;
176  bool accept = (t<gxsec);
177  if(accept) break; // exit loop
178  } // 1
179  } // generate uniformly
180 
181  LOG("CEvNS", pNOTICE) << "Selected Q2 = " << gQ2 << " GeV^2";
182 
183  // reset bits
184  interaction->ResetBit(kISkipProcessChk);
185  interaction->ResetBit(kISkipKinematicChk);
186 
187  // lock selected kinematics & clear running values
188  interaction->KinePtr()->SetQ2(gQ2, true);
189  interaction->KinePtr()->ClearRunningValues();
190 
191  // Set the cross section for the selected kinematics
192  event->SetDiffXSec(gxsec,kPSQ2fE);
193 }
194 //___________________________________________________________________________
196 {
197  GHepParticle * probe = event->Probe();
198  GHepParticle * target = event->TargetNucleus();
199 
200  int target_pdgc = target->Pdg();
201  double M = PDGLibrary::Instance()->Find(target_pdgc)->Mass(); // units: GeV
202 
203  double Ev = probe->E(); // neutrino energy, units: GeV
204  double Q2 = event->Summary()->Kine().Q2(true); // selected momentum transfer, units: GeV^2
205 
206  const TLorentzVector & vtx = *(probe->X4());
207  TLorentzVector x4l(vtx); // position 4-vector
208 
209  // Compute the final state neutino energy
210 
211  double y = Q2/(2*M*Ev);
212  double El = (1-y)*Ev;
213 
214  LOG("CEvNS", pNOTICE)
215  << "Final state neutrino energy: E = " << El << " GeV";
216 
217  // Compute the final state neutrino momentum components
218  // along and perpendicular to the incoming neutrino direction
219 
220  double ml2 = 0;
221  double plp = El - 0.5*(Q2+ml2)/Ev; // p(//)
222  double plt = TMath::Sqrt(TMath::Max(0.,El*El-plp*plp-ml2)); // p(-|)
223 
224  LOG("CEvNS", pNOTICE)
225  << "Final state neutrino momentum components: |p//| = "
226  << plp << " GeV, [pT] = " << plt << " GeV";
227 
228  // Randomize transverse components
229  RandomGen * rnd = RandomGen::Instance();
230  double phi = 2*kPi * rnd->RndLep().Rndm();
231  double pltx = plt * TMath::Cos(phi);
232  double plty = plt * TMath::Sin(phi);
233 
234  // Take a unit vector along the neutrino direction @ the LAB
235  TVector3 unit_nudir = probe->P4()->Vect().Unit();
236 
237  // Rotate lepton momentum vector from the reference frame (x'y'z') where
238  // {z':(neutrino direction), z'x':(theta plane)} to the LAB
239  TVector3 p3l(pltx,plty,plp);
240  p3l.RotateUz(unit_nudir);
241 
242  // Lepton 4-momentum in the LAB
243  TLorentzVector p4l(p3l,El);
244 
245  LOG("CEvNS", pNOTICE)
246  << "Final state neutrino 4-momentum: " << utils::print::P4AsString(&p4l);
247 
248  event->AddParticle(
249  probe->Pdg(), kIStStableFinalState, event->ProbePosition(),
250  -1,-1,-1, p4l, x4l);
251 
252  event->Summary()->KinePtr()->SetFSLeptonP4(p4l);
253 }
254 //___________________________________________________________________________
256 {
257  GHepParticle * probe = event->Probe();
258  GHepParticle * target = event->TargetNucleus();
259  GHepParticle * fsl = event->Particle(probe->FirstDaughter());
260 
261  const TLorentzVector & p4probe = * ( probe -> P4() );
262  const TLorentzVector & p4target = * ( target -> P4() );
263  const TLorentzVector & p4fsl = * ( fsl -> P4() );
264 
265  const TLorentzVector & p4recoil = p4probe + p4target - p4fsl;
266 
267  LOG("CEvNS", pNOTICE)
268  << "Recoil 4-momentum: " << utils::print::P4AsString(&p4recoil);
269 
270  const TLorentzVector & vtx = *(probe->X4());
271 
272  event->AddParticle(
273  event->TargetNucleus()->Pdg(),
275  event->TargetNucleusPosition(),
276  -1,-1,-1, p4recoil, vtx);
277 }
278 //___________________________________________________________________________
280 {
281  Algorithm::Configure(config);
282  this->LoadConfig();
283 }
284 //____________________________________________________________________________
286 {
287  Algorithm::Configure(config);
288  this->LoadConfig();
289 }
290 //____________________________________________________________________________
292 {
293  fXSecModel = 0;
294 
295  // max xsec safety factor (for rejection method) and min cached energy
296  this->GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 1.05 ) ;
297 
298  // Generate kinematics uniformly over allowed phase space and compute
299  // an event weight?
300  this->GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false ) ;
301 
302  // Maximum allowed fractional cross section deviation from maxim cross
303  // section used in rejection method
304  this->GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ;
305  assert(fMaxXSecDiffTolerance>=0);
306 }
307 //____________________________________________________________________________
const XSecAlgorithmI * fXSecModel
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
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
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
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...
int FirstDaughter(void) const
Definition: GHepParticle.h:68
A simple [min,max] interval for doubles.
Definition: Range1.h:42
#define pFATAL
Definition: Messenger.h:56
Defines the EventGeneratorI interface.
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:27
Range1D_t Q2Lim(void) const
Q2 limits.
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
void AddRecoilNucleus(GHepRecord *event) const
int Pdg(void) const
Definition: GHepParticle.h:63
Summary information for an interaction.
Definition: Interaction.h:56
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 Configure(const Registry &config)
static constexpr double cm2
Definition: Units.h:69
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
Kinematical phase space.
Definition: KPhaseSpace.h:33
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:293
#define pINFO
Definition: Messenger.h:62
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:50
void AddFinalStateNeutrino(GHepRecord *event) const
#define pWARN
Definition: Messenger.h:60
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
void GenerateKinematics(GHepRecord *event) const
static RunningThreadInfo * Instance(void)
double max
Definition: Range1.h:53
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
void ProcessEventRecord(GHepRecord *event) const
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
double gQ2
Definition: gtestDISSF.cxx:55
double min
Definition: Range1.h:52
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
#define pNOTICE
Definition: Messenger.h:61
void ClearRunningValues(void)
Definition: Kinematics.cxx:347
bool GetParamDef(const RgKey &name, T &p, const T &def) const
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
const EventGeneratorI * RunningThread(void)
double ProbeE(RefFrame_t rf) const
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
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 UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:48