VertexGenerator.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  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Sep 19, 2007 - CA
14  Added code to generate the interaction according to a realistic nuclear
15  density and made that the default setting.
16  @ Dec 01, 2007 - CA
17  For COH and ve- interactions setting the vertex on the nuclear boundary
18  rather than inside the nucleus.
19  @ Sep 15, 2009 - CA
20  IsFake() and IsNucleus() are no longer available in GHepParticle.
21  Use pdg::IsPseudoParticle() and pdg::IsIon().
22  @ Feb 12, 2013 - CA (code from Rosen Matev)
23  Handle the IMD annihilation channel.
24  @ Jun 03, 2016 - JJ (SD)
25  Move code to generate the position to a public method, so this code
26  can easily be reused by other classes. (Specifically LwlynSmithQELCCPXSec
27  and NievesQELCCPXSec to generate a position before calculating the xsec
28  when making splines).
29 */
30 //____________________________________________________________________________
31 
32 #include <TMath.h>
33 
48 
49 using namespace genie;
50 using namespace genie::utils;
51 using namespace genie::controls;
52 using namespace genie::constants;
53 
54 //___________________________________________________________________________
56 EventRecordVisitorI("genie::VertexGenerator")
57 {
58 
59 }
60 //___________________________________________________________________________
62 EventRecordVisitorI("genie::VertexGenerator", config)
63 {
64 
65 }
66 //___________________________________________________________________________
68 {
69 
70 }
71 //___________________________________________________________________________
73 {
74 // generate a vtx and set it to all GHEP physical particles
75  Interaction * interaction = evrec->Summary();
76  GHepParticle * nucltgt = evrec->TargetNucleus();
77  TVector3 vtx(9999999.,999999.,999999.);
78  if(!nucltgt){
79  vtx.SetXYZ(0.,0.,0.);
80  }else{
81  double A = nucltgt->A();
82  vtx = GenerateVertex(interaction,A);
83  }
84 
85  // Copy the vertex info to the particles already in the event record
86  //
87  TObjArrayIter piter(evrec);
88  GHepParticle * p = 0;
89  while( (p = (GHepParticle *) piter.Next()) )
90  {
91  if(pdg::IsPseudoParticle(p->Pdg())) continue;
92  if(pdg::IsIon (p->Pdg())) continue;
93 
94  LOG("Vtx", pINFO) << "Setting vertex position for: " << p->Name();
95  p->SetPosition(vtx.x(), vtx.y(), vtx.z(), 0.);
96  }
97 }
98 //___________________________________________________________________________
100 {
101  Algorithm::Configure(config);
102  this->LoadConfig();
103 }
104 //____________________________________________________________________________
106 {
107  Algorithm::Configure(config);
108  this->LoadConfig();
109 }
110 //____________________________________________________________________________
112 {
113 
114  GetParam( "VtxGenerationMethod", fVtxGenMethod ) ;
115  GetParam( "NUCL-R0", fR0 ) ; //fm
116 
117 }
118 //____________________________________________________________________________
119 TVector3 VertexGenerator::GenerateVertex(const Interaction * interaction,
120  double A) const{
122  TVector3 vtx(999999.,999999.,999999.);
123 
124  //GHepParticle * nucltgt = evrec->TargetNucleus();
125 
126  bool uniform = fVtxGenMethod==0;
127  bool realistic = !uniform;
128 
129  //if(!nucltgt) {
130  //vtx.SetXYZ(0.,0.,0.);
131  //}
132  //else {
133  //double A = nucltgt->A();
134  double R = fR0 * TMath::Power(A, 1./3.);
135 
136  //Interaction * interaction = evrec->Summary();
137  const ProcessInfo & proc_info = interaction->ProcInfo();
138  bool is_coh = proc_info.IsCoherent() || proc_info.IsCoherentElas();
139  bool is_ve = proc_info.IsInverseMuDecay() ||
140  proc_info.IsIMDAnnihilation() ||
141  proc_info.IsNuElectronElastic();
142 
143  if(is_coh||is_ve) {
144  // ** For COH or ve- set a vertex positon on the nuclear boundary
145  //
146  LOG("Vtx", pINFO) << "Setting vertex on the nuclear boundary";
147  double phi = 2*kPi * rnd->RndFsi().Rndm();
148  double cosphi = TMath::Cos(phi);
149  double sinphi = TMath::Sin(phi);
150  double costheta = -1 + 2 * rnd->RndFsi().Rndm();
151  double sintheta = TMath::Sqrt(1-costheta*costheta);
152  vtx.SetX(R*sintheta*cosphi);
153  vtx.SetY(R*sintheta*sinphi);
154  vtx.SetZ(R*costheta);
155  }
156  else {
157  // ** For other events on nuclear targets set the interaction vertex
158  // ** using the requested method: either using a realistic nuclear
159  // ** density or by setting it uniformly within the nucleus
160  //
161  if(realistic) {
162  // Generate the vertex using a realistic nuclear density
163  //
164  LOG("Vtx", pINFO)
165  << "Generating vertex according to a realistic nuclear density profile";
166  // get inputs to the rejection method
167  double ymax = -1;
168  double rmax = 3*R;
169  double dr = R/40.;
170  for(double r = 0; r < rmax; r+=dr) {
171  ymax = TMath::Max(ymax, r*r * utils::nuclear::Density(r,(int)A));
172  }
173  ymax *= 1.2;
174 
175  // select a vertex using the rejection method
176  unsigned int iter = 0;
177  while(1) {
178  iter++;
179 
180  // throw an exception if it hasn't find a solution after many attempts
181  if(iter > kRjMaxIterations) {
182  LOG("Vtx", pWARN)
183  << "Couldn't generate a vertex position after " << iter << " iterations";
184  //evrec->EventFlags()->SetBitNumber(kGenericErr, true);
186  exception.SetReason("Couldn't generate vertex");
187  exception.SwitchOnFastForward();
188  throw exception;
189  }
190 
191  double r = rmax * rnd->RndFsi().Rndm();
192  double t = ymax * rnd->RndFsi().Rndm();
193  double y = r*r * utils::nuclear::Density(r,(int)A);
194  if(y > ymax) {
195  LOG("Vtx", pERROR)
196  << "y = " << y << " > ymax = " << ymax
197  << " for r = " << r << ", A = " << A;
198  }
199  bool accept = (t < y);
200  if(accept) {
201  double phi = 2*kPi * rnd->RndFsi().Rndm();
202  double cosphi = TMath::Cos(phi);
203  double sinphi = TMath::Sin(phi);
204  double costheta = -1 + 2 * rnd->RndFsi().Rndm();
205  double sintheta = TMath::Sqrt(1-costheta*costheta);
206  vtx.SetX(r*sintheta*cosphi);
207  vtx.SetY(r*sintheta*sinphi);
208  vtx.SetZ(r*costheta);
209  break;
210  }
211  }//w(1)
212  } //use density?
213 
214  if(uniform) {
215  // Generate the vertex uniformly within the nucleus
216  //
217  LOG("Vtx", pINFO)
218  << "Generating intranuclear vertex uniformly in volume";
219  while(vtx.Mag() > R) {
220  vtx.SetX(-R + 2*R * rnd->RndFsi().Rndm());
221  vtx.SetY(-R + 2*R * rnd->RndFsi().Rndm());
222  vtx.SetZ(-R + 2*R * rnd->RndFsi().Rndm());
223  }
224  }// uniform?
225 
226  } // coh or ve-?
227  //} // nuclear target ?
228 
229  LOG("Vtx", pINFO)
230  << "Generated vtx @ r = " << vtx.Mag() << " fm / "
231  << print::Vec3AsString(&vtx);
232  return vtx;
233 }
const double kPi
TVector3 GenerateVertex(const Interaction *in, double A) const
int fVtxGenMethod
vtx generation method (0: uniform, 1: according to nuclear density [def])
double fR0
parameter controlling nuclear sizes
Basic constants.
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
Definition: RandomGen.h:60
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:135
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
double Density(double r, int A, double ring=0.)
bool IsInverseMuDecay(void) const
const char * p
Definition: xmltok.h:285
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
Definition: config.py:1
bool IsIMDAnnihilation(void) const
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
int Pdg(void) const
Definition: GHepParticle.h:64
Double_t ymax
Definition: plot.C:25
string Name(void) const
Name that corresponds to the PDG code.
Summary information for an interaction.
Definition: Interaction.h:56
void SetPosition(const TLorentzVector &v4)
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
bool IsNuElectronElastic(void) const
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
#define R(x)
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:330
#define pINFO
Definition: Messenger.h:63
Misc GENIE control constants.
#define pWARN
Definition: Messenger.h:61
void Configure(const Registry &config)
static const double A
Definition: Units.h:82
bool IsIon(int pdgc)
Definition: PDGUtils.cxx:40
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
bool IsPseudoParticle(int pdgc)
Definition: PDGUtils.cxx:25
TRandom3 r(0)
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
int A(void) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
bool IsCoherentElas(void) const
static const unsigned int kRjMaxIterations
Definition: Controls.h:27
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
string Vec3AsString(const TVector3 *vec)
Definition: PrintUtils.cxx:87
Root of GENIE utility namespaces.
void ProcessEventRecord(GHepRecord *event_rec) const
bool IsCoherent(void) const
Definition: ProcessInfo.cxx:97