GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
COHDNuEventGenerator.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  Author: Iker de Icaza <i.de-icaza-astiz \at sussex.ac.uk>
7  University of Sussex
8 
9  Costas Andreopoulos <c.andreopoulos \at cern.ch>
10  University of Liverpool
11 */
12 //____________________________________________________________________________
13 
14 #include <TMath.h>
15 #include <Math/IFunction.h>
16 #include <Math/GSLMinimizer1D.h>
17 #include <Math/BrentMinimizer1D.h>
18 
39 
40 using namespace genie;
41 using namespace genie::constants;
42 using namespace genie::controls;
43 using namespace genie::utils;
44 
45 //___________________________________________________________________________
47  EventRecordVisitorI("genie::COHDNuEventGenerator")
48 {
49 
50 }
51 //___________________________________________________________________________
53  EventRecordVisitorI("genie::COHDNuEventGenerator", config)
54 {
55 
56 }
57 //___________________________________________________________________________
59 {
60 
61 }
62 //___________________________________________________________________________
64 {
65  this -> GenerateKinematics (event);
66  this -> AddFinalStateDarkNeutrino (event);
67  this -> AddRecoilNucleus (event);
68 }
69 //___________________________________________________________________________
71 {
72  Interaction * interaction = event->Summary();
73  interaction->SetBit(kISkipProcessChk);
74  interaction->SetBit(kISkipKinematicChk);
75 
76  // Get the random number generators
78 
79  // Get the kinematical limits
80  const InitialState & init_state = interaction -> InitState();
81  double E_nu = init_state.ProbeE(kRfLab);
82 
83  // Access cross section algorithm for running thread
85  const EventGeneratorI * evg = rtinfo->RunningThread();
86  fXSecModel = evg->CrossSectionAlg();
87 
88  double gDNuE = -1; // generated Dark Neutrino Energy
89  double gxsec = -1; // dsig/dDENu at generated DNuE
90 
91  // Generate kinematics
92  if(fGenerateUniformly) {
93  LOG("COHDNu", pFATAL)
94  << "Option to generate kinematics uniformly not supported";
95  exit(1);
96  }
97  else {
98  // For the subsequent kinematic selection with the rejection method:
99  // Calculate the max differential cross section.
100  // Always at Q^2 = 0 for energies and model tested,
101  // but go ahead and do the calculation nevertheless.
102  utils::gsl::dXSec_dEDNu_E xsec_func(fXSecModel, interaction, fDNuMass, -1.);
103  Range1D_t DNuEnergy = xsec_func.IntegrationRange();
104  double dDNuE = DNuEnergy.max - DNuEnergy.min;
105  ROOT::Math::BrentMinimizer1D minimizer;
106  minimizer.SetFunction( xsec_func, DNuEnergy.min, DNuEnergy.max);
107  minimizer.Minimize(1000, 1, 1E-5);
108  double DNuE_for_xsec_max = minimizer.XMinimum();
109  double xsec_max = -1. * xsec_func(DNuE_for_xsec_max); // xsec in units of 1E-38 cm2/GeV
110 
111  // Try to select a valid E_N
112  unsigned int iter = 0;
113  while(1) {
114  iter++;
115  if(iter > kRjMaxIterations) {
116  LOG("COHDNu", pWARN)
117  << "*** Could not select a valid DNuE after " << iter << " iterations";
118  event->EventFlags()->SetBitNumber(kKineGenErr, true);
120  exception.SetReason("Couldn't select kinematics");
121  exception.SwitchOnFastForward();
122  throw exception;
123  } // max iterations
124 
125  gDNuE = DNuEnergy.min + dDNuE * rnd->RndKine().Rndm();
126  LOG("COHDNu", pINFO) << "Trying: E_N = " << gDNuE;
127 
128  // Computing cross section for the current kinematics
129  gxsec = -1. * xsec_func(gDNuE);
130 
131  if(gxsec > xsec_max) {
132  double frac = TMath::Abs(gxsec-xsec_max)/xsec_max;
133  if(frac > fMaxXSecDiffTolerance) {
134  LOG("COHDNu", pWARN)
135  << "Current computed cross-section (" << gxsec/(units::cm2)
136  << " cm2/GeV^2) exceeds the maximum cross-section ("
137  << xsec_max/(units::cm2) << " beyond the specified tolerance";
138  }
139  }
140 
141  // Decide whether to accept the current kinematic point
142  double t = fSafetyFactor * xsec_max * rnd->RndKine().Rndm();
143  //this->AssertXSecLimits(interaction, gxsec, xsec_max);
144  LOG("COHDNu", pINFO)
145  << "dxsec/dQ2 = " << gxsec/(units::cm2) << " cm2/GeV^2"
146  << "J = 1, rnd = " << t;
147  bool accept = (t<gxsec);
148  if(accept) break; // exit loop
149  } // 1
150  } // generate uniformly
151 
152  // reset bits
153  interaction->ResetBit(kISkipProcessChk);
154  interaction->ResetBit(kISkipKinematicChk);
155 
156  double M_target = interaction->InitState().Tgt().Mass();
157 
158  double ETimesM = E_nu * M_target;
159  double EPlusM = E_nu + M_target;
160 
161  double p_DNu = TMath::Sqrt(gDNuE*gDNuE - fDNuMass2);
162  double cos_theta_DNu = (gDNuE*EPlusM - ETimesM - 0.5*fDNuMass2) / (E_nu * p_DNu);
163  double theta_DNu = TMath::ACos(cos_theta_DNu);
164 
165  // Take a unit vector along the neutrino direction @ the LAB
166  GHepParticle * probe = event->Probe();
167  TVector3 unit_nudir = probe->P4()->Vect().Unit();
168 
169  TVector3 DNu_3vector = TVector3(0,0,0);
170  double phi = 2.*kPi * rnd->RndKine().Rndm();
171  DNu_3vector.SetMagThetaPhi(p_DNu, theta_DNu, phi);
172  DNu_3vector.RotateUz(unit_nudir);
173  TLorentzVector P4_DNu = TLorentzVector(DNu_3vector, gDNuE);
174  interaction->KinePtr()->SetFSLeptonP4(P4_DNu);
175 
176  // lock selected kinematics & clear running values
177  double gQ2 = -(P4_DNu - *(probe->P4())).M2();
178  interaction->KinePtr()->SetQ2(gQ2, true);
179  interaction->KinePtr()->ClearRunningValues();
180 
181  // Set the cross section for the selected kinematics
182  event->SetDiffXSec(gxsec * (1E-38*units::cm2), kPSEDNufE);
183 }
184 //___________________________________________________________________________
186 {
187  GHepParticle * probe = event->Probe();
188 
189  const TLorentzVector & vtx = *(probe->X4());
190  TLorentzVector x4l(vtx); // position 4-vector
191 
192  event->AddParticle(probe -> Pdg() > 0 ? kPdgDarkNeutrino : kPdgAntiDarkNeutrino,
194  event->ProbePosition(), event->TargetNucleusPosition(),
195  -1,-1, event->Summary()->Kine().FSLeptonP4(), x4l);
196 }
197 //___________________________________________________________________________
199 {
200  GHepParticle * probe = event->Probe();
201  GHepParticle * target = event->TargetNucleus();
202  GHepParticle * fsl = event->Particle(probe->FirstDaughter());
203 
204  const TLorentzVector & p4probe = * ( probe -> P4() );
205  const TLorentzVector & p4target = * ( target -> P4() );
206  const TLorentzVector & p4fsl = * ( fsl -> P4() );
207 
208  const TLorentzVector & p4recoil = p4probe + p4target - p4fsl;
209 
210  LOG("COHDNu", pNOTICE)
211  << "Recoil 4-momentum: " << utils::print::P4AsString(&p4recoil);
212 
213  const TLorentzVector & vtx = *(probe->X4());
214 
215  event->AddParticle(event->TargetNucleus()->Pdg(),
217  event->TargetNucleusPosition(), event->ProbePosition(),
218  -1,-1, p4recoil, vtx);
219 }
220 //___________________________________________________________________________
222 {
223  Algorithm::Configure(config);
224  this->LoadConfig();
225 }
226 //____________________________________________________________________________
228 {
229  Algorithm::Configure(config);
230  this->LoadConfig();
231 }
232 //____________________________________________________________________________
234 {
235  fXSecModel = 0;
236 
237  // max xsec safety factor (for rejection method) and min cached energy
238  this->GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 1.05 ) ;
239 
240  // Generate kinematics uniformly over allowed phase space and compute
241  // an event weight?
242  this->GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false ) ;
243 
244  // Maximum allowed fractional cross section deviation from maxim cross
245  // section used in rejection method
246  this->GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ;
247  assert(fMaxXSecDiffTolerance>=0);
248 
249  fDNuMass = 0.;
250  this->GetParam("Dark-NeutrinoMass", fDNuMass);
252 
253 }
254 //____________________________________________________________________________
void ProcessEventRecord(GHepRecord *event) const
Range1D_t IntegrationRange(void) const
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
void GenerateKinematics(GHepRecord *event) const
#define pFATAL
Definition: Messenger.h:56
Defines the EventGeneratorI interface.
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:27
const XSecAlgorithmI * fXSecModel
double Mass(void) const
Definition: Target.cxx:224
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
const int kPdgAntiDarkNeutrino
Definition: PDGCodes.h:222
int Pdg(void) const
Definition: GHepParticle.h:63
Summary information for an interaction.
Definition: Interaction.h:56
void AddRecoilNucleus(GHepRecord *event) const
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void SetFSLeptonP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:297
static constexpr double cm2
Definition: Units.h:69
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:293
#define pINFO
Definition: Messenger.h:62
void AddFinalStateDarkNeutrino(GHepRecord *event) const
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:50
#define pWARN
Definition: Messenger.h:60
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
static RunningThreadInfo * Instance(void)
double max
Definition: Range1.h:53
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
double gQ2
Definition: gtestDISSF.cxx:55
const InitialState & InitState(void) const
Definition: Interaction.h:69
double min
Definition: Range1.h:52
#define pNOTICE
Definition: Messenger.h:61
void ClearRunningValues(void)
Definition: Kinematics.cxx:347
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
virtual int ProbePosition(void) const
Definition: GHepRecord.cxx:352
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
const int kPdgDarkNeutrino
Definition: PDGCodes.h:221
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 UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
void Configure(const Registry &config)
Initial State information.
Definition: InitialState.h:48