QwGeant4
QweakSimEPEvent.hh
Go to the documentation of this file.
1 //=============================================================================
2 //
3 // ---------------------------
4 // | Doxygen File Information |
5 // ---------------------------
6 //
7 /**
8 
9  \file QweakSimEPEvent.hh
10 
11  \author Jie Pan
12 
13 */
14 //=============================================================================
15 
16 #ifndef QWEAKSIMEPEVENT_H
17 #define QWEAKSIMEPEVENT_H
18 
19 // geant4 includes
20 #include "Randomize.hh"
21 #include "G4Types.hh"
22 #include "G4ThreeVector.hh"
23 
24 // user includes
25 #include "QweakSimSystemOfUnits.hh"
26 #include "QweakSimFieldMap.hh"
27 
28 // user classes
31 
32 //template <class value_t, unsigned int value_n>
33 //class QweakSimFieldMap<value_t,value_n>;
34 
36 
37  private:
38  //define target geometry parameter,
39  //TODO: these parameters will be acquired from QweakSimTarget
40 
41  static const G4double TargetLength;
42 
43  // definition of a mil = inch/1000
44  static const G4double mil;
45  static const G4double TargetWindowThickness;
46 
47  static const G4double M_n; //neutron mass in MeV/c^2
48  static const G4double M_p; // proton mass in MeV/c^2
49 
50  static const G4double Theta_Min;
51 
52  G4int TypeSetting;
53  G4int ReactionType; ///< \ref reaction_type used for event generation
54  G4int ReactionRegion; ///< \ref reaction_region used for event generation
55  G4int kActiveOctantNumber; ///< Active octant number in the simulation, 0 will enable all octants
56 
57  G4int Isotropy; ///< \ref isotropy used for event generation
58  //G4double PhiAngle_Min;
59  //G4double PhiAngle_Max;
60  //G4double ThetaAngle_Min;
61  //G4double ThetaAngle_Max;
62  //G4double EPrime_Min;
63  //G4double EPrime_Max;
65  G4double SchwingerDeltaE;
66  G4double BeamEnergy;
67 
68  G4double myPositionZ;
69 
70  // Lookup Table storage type
71  typedef float value_t; // type of all values stored in the field
72 
73  // Lookup Table field values
74  static const G4int value_n = 15; // number of values at each point in the field
75  static const G4int value_d = 4; // number of dimensions of the coordinates
76  //std::vector< G4double > fMin;
77  //std::vector< G4double > fMax;
78  //std::vector< G4double > fStep;
79 
81 
82  G4ThreeVector GetMomentumDirection();
83  G4int SuperElasticCheck(G4double E_in, G4double E_out, G4double theta, G4double &xsec);
84  G4double ResMod507(G4int sf,G4double w2,G4double q2,G4double *xval);
85  G4double Sigma_EEPrime(G4double eni,G4double eprime,G4double theta, G4double &q2);
86 
89 
90  public:
91 
93  //QweakSimEPEvent();
94  virtual ~QweakSimEPEvent();
95 
96  void SetActiveOctantNumber(G4int kaot) { kActiveOctantNumber = kaot; };
98 
99  void SetIsotropy(G4int isot) { Isotropy = isot; };
100  G4int GetIsotropy( ) {return Isotropy;};
101  //------------------------------------------------------------------------
102  /*
103  void SetEPrime_Min(G4double energy) {myUserInfo->SetEPrime_Min(energy); CheckLookupTableBounds();};
104  G4double GetEPrime_Min() {return myUserInfo->GetEPrime_Min();};
105 
106  void SetEPrime_Max(G4double energy) {myUserInfo->SetEPrime_Max(energy); CheckLookupTableBounds();};
107  G4double GetEPrime_Max() {return myUserInfo->GetEPrime_Max();};
108 
109  void SetThetaAngle_Min(G4double ang) {myUserInfo->SetThetaAngle_Min(ang); CheckLookupTableBounds();};
110  G4double GetThetaAngle_Min() {return myUserInfo->GetThetaAngle_Min();};
111 
112  void SetThetaAngle_Max(G4double ang) {myUserInfo->SetThetaAngle_Max(ang); CheckLookupTableBounds();};
113  G4double GetThetaAngle_Max() {return myUserInfo->GetThetaAngle_Max();};
114 
115  void SetPhiAngle_Min(G4double ang) {myUserInfo->SetPhiAngle_Min(ang);};
116  G4double GetPhiAngle_Min() {return myUserInfo->GetPhiAngle_Min();};
117 
118  void SetPhiAngle_Max(G4double ang) {myUserInfo->SetPhiAngle_Max(ang);};
119  G4double GetPhiAngle_Max() {return myUserInfo->GetPhiAngle_Max();};
120 
121  void SetLuminosity(G4double lum) {myUserInfo->SetLuminosity(lum);};
122  G4double GetLuminosity() {return myUserInfo->GetLuminosity();};
123 
124  void SetPhaseSpace(G4double ps) {myUserInfo->SetPhaseSpace(ps);}
125  G4double GetPhaseSpace() {return myUserInfo->GetPhaseSpace();};
126 
127  void SetBeamEnergy(G4double energy = 1.160*GeV);
128  G4double GetBeamEnergy() {return myUserInfo->GetBeamEnergy();}
129  */
130  //----------------------------------------------------------------------
131  void SetEPrime_Min(G4double energy);
132  G4double GetEPrime_Min();
133 
134  void SetEPrime_Max(G4double energy);
135  G4double GetEPrime_Max();
136 
137  void SetThetaAngle_Min(G4double ang);
138  G4double GetThetaAngle_Min();
139 
140  void SetThetaAngle_Max(G4double ang);
141  G4double GetThetaAngle_Max();
142 
143  void SetPhiAngle_Min(G4double ang);
144  G4double GetPhiAngle_Min();
145 
146  void SetPhiAngle_Max(G4double ang);
147  G4double GetPhiAngle_Max();
148 
149  void SetLuminosity(G4double lum);
150  G4double GetLuminosity();
151 
152  void SetPhaseSpace(G4double ps);
153  G4double GetPhaseSpace();
154 
155  void SetBeamEnergy(G4double energy = 1.160*GeV);
156  G4double GetBeamEnergy();
157  //-----------------------------------------------------------------
158  // Both these values should be the same unless SchwingerDeltaE < 15 MeV.
159 
160  // Sets Radiative Cutoff for Generator 7 Lookuptables.
161  // Default and minumum value is 15 MeV.
162  void SetElasticPeakDeltaE(G4double energy = 15*MeV) { if (energy < 15*MeV) { energy = 15*MeV; G4cout << "***** Value Out of Range: ElasticPeakDeltaE set to : 15 MeV *****" << G4endl; } ElasticPeakDeltaE = energy; } ;
163  G4double GetElasticPeakDeltaE() { return ElasticPeakDeltaE; } ;
164 
165  // Sets DeltaE in the Schwinger Correction to elastic cross sections
166  // Default Value is 15 MeV
167  void SetSchwingerDeltaE(G4double energy = 15*MeV) { SchwingerDeltaE = energy; } ;
168  G4double GetSchwingerDeltaE() { return SchwingerDeltaE; } ;
169 
170  //-----------------------------------------------------------------
171  void SetReactionType(G4int rt) { ReactionType = rt; TypeSetting = rt;};
172  G4int GetReactionType() {return ReactionType; };
173 
175  G4int GetReactionRegion() {return ReactionRegion; };
176 
177  G4double GetVertexZ();
178  void GetanEvent(G4double E_in,
179  std::vector< G4double > &CrossSection,
180  G4double &weight_n,
181  G4double &Q2,
182  G4double &E_out,
183  G4ThreeVector &MomentumDirection,
184  G4double &theta,
185  G4double &phi,
186  G4double &Asymmetry);
187 
188  G4double Elastic_Cross_Section_Proton( G4double E_in,
189  G4double Theta,
190  G4double &fWeightN,
191  G4double &Q2,
192  G4double &E_out);
193 
194  G4double Elastic_Cross_Section_Aluminum( G4double E_in,
195  G4double Theta,
196  G4double &fWeightN,
197  G4double &Q2,
198  G4double &E_out);
199 
200  G4double Quasi_Elastic_Neutron( G4double E_in,
201  G4double Theta,
202  G4double &fWeightN,
203  G4double &Q2,
204  G4double &E_out);
205 
206  G4double Delta_Resonance(G4double E_in,
207  G4double Theta,
208  G4double &fWeightN,
209  G4double &Q2,
210  G4double &E_out);
211 
212  G4double Moller_Scattering(G4double E_in, G4double theta1,
213  G4double &E_out1, G4double &E_out2, G4double &theta2,
214  G4double &q2, G4double &fWeightN, G4double &asymmetry);
215 
216  const std::vector< G4double > Radiative_Cross_Section_Lookup( G4double E_in,
217  G4double Theta,
218  G4double &fWeightN,
219  G4double &Q2,
220  G4double &E_out);
221 
222  void CreateLookupTable();
223 
224  void CheckLookupTableBounds();
225 
226  G4double AlNuclInel(G4double E_in,
227  G4double Theta,
228  G4double &fWeightN,
229  G4double &Q2,
230  G4double &E_out);
231 
232  G4double AlGDR( G4double E_in,
233  G4double Theta,
234  G4double &fWeightN,
235  G4double &Q2,
236  G4double &E_out);
237 
238  G4double Elastic_Cross_Section_Carbon( G4double E_in,
239  G4double Theta,
240  G4double &fWeightN,
241  G4double &Q2,
242  G4double &E_out);
243 
244  G4double CGDR( G4double E_in,
245  G4double Theta,
246  G4double &fWeightN,
247  G4double &Q2,
248  G4double &E_out );
249 
250  G4double CNuclInel( G4double E_in,
251  G4double Theta,
252  G4int nState,
253  G4double &fWeightN,
254  G4double &Q2,
255  G4double &E_out);
256 
257 
258  G4double Pion_PhotoProduction(G4double E_in,
259  G4double Theta,
260  G4double &fWeightN,
261  G4double &Q2,
262  G4double &E_out);
263 
264  G4double Pion_PhotoProductionAl(G4double E_in,
265  G4double Theta,
266  G4double &fWeightN,
267  G4double &Q2,
268  G4double &E_out);
269 
270  G4double Pion_PhotoProductionCarbon(G4double E_in,
271  G4double Theta,
272  G4double &fWeightN,
273  G4double &Q2,
274  G4double &E_out);
275 
276  G4double Quasi_Elastic_Bosted(G4double E_in, G4double Theta, G4int Zin,
277  G4int Ain, G4double &fWeightN, G4double &Q2, G4double &E_out);
278  G4double NuclearInelastic_Bosted(G4double E_in, G4double Theta,
279  G4int Zin, G4int Ain,
280  G4double &fWeightN, G4double &Q2,
281  G4double &E_out);
282 
283  G4double Aluminum_Excited_State(G4double E_in,
284  G4double Theta,
285  G4double E_lvl,
286  G4double fit_c,
287  G4double fit_mean,
288  G4double fit_sigma,
289  G4double &fWeightN,
290  G4double &Q2,
291  G4double &E_out);
292 
293  G4double Horowitz_DW_Xsect(G4double E_in,
294  G4double theta,
295  const G4String path,
296  G4double &fWeightN,
297  G4double &Q2,
298  G4double &E_out);
299 
300  G4double Horowitz_DW_Asym(G4double theta,
301  const G4String path);
302 
303  void F1F2QE09(G4int Z, G4int A, G4double QSQ,
304  G4double wsq, G4double &F1, G4double &F2);
305  void F1F2IN09(G4int Z, G4int A, G4double QSQ,
306  G4double wsq, G4double &F1, G4double &F2);
307 
308  void christy507(G4double wsq, G4double Q2, G4double &F1,
309  G4double &R, G4double &sigt, G4double &sigl);
310  void resmodd(G4double w2, G4double q2, G4double xval[50], G4double &sig);
311  G4double resmod507_v2(G4double sf,G4double w2, G4double q2,
312  G4double xval[50]);
313  G4double MEC2009(G4int a, G4double q2,G4double w2);
314  G4double fitemc(G4double X, G4int A);
315 
316  G4double AlloyScattering(G4double E_in, G4double Theta, G4int Zin,
317  G4int Ain, G4double &fWeightN, G4double &Q2,
318  G4double &E_out);
319  G4double Fshell(G4int Z, G4int A, G4double q2);
320  G4double Fgauss(G4int Z, G4int A, G4double q2);
321  G4double FF_BESSEL(G4int Z, G4int A, G4double q2, G4bool &ofr);
322 
323  G4double GetAsymmetry_EP(G4double theta, G4double energy);
324  G4double GetAsymmetry_EN(G4double theta, G4double energy);
325  G4double GetAsymmetry_AL(G4double theta, G4double energy);
326  G4double GetAsymmetry_Be(G4double theta, G4double energy);
327  G4double GetAsymmetry_Pi(G4double Q2_pi);
328  G4double Born_Approx_Nuclei_Asym(G4double Z,
329  G4double A,
330  G4double energy,
331  G4double eprime,
332  G4double qq);
333 
334 };
335 
336 #endif
void GetanEvent(G4double E_in, std::vector< G4double > &CrossSection, G4double &weight_n, G4double &Q2, G4double &E_out, G4ThreeVector &MomentumDirection, G4double &theta, G4double &phi, G4double &Asymmetry)
G4int ReactionRegion
reaction region used for event generation
QweakSimEPEventMessenger * EventGen_Messenger
G4double SchwingerDeltaE
G4double Quasi_Elastic_Bosted(G4double E_in, G4double Theta, G4int Zin, G4int Ain, G4double &fWeightN, G4double &Q2, G4double &E_out)
G4int ReactionType
reaction type used for event generation
void SetEPrime_Max(G4double energy)
G4double GetThetaAngle_Max()
QweakSimFieldMap< value_t, value_n > * fLookupTable
void SetEPrime_Min(G4double energy)
G4double GetSchwingerDeltaE()
void SetReactionType(G4int rt)
void SetPhaseSpace(G4double ps)
G4ThreeVector GetMomentumDirection()
G4double GetThetaAngle_Min()
G4double Fshell(G4int Z, G4int A, G4double q2)
G4double fitemc(G4double X, G4int A)
G4double Moller_Scattering(G4double E_in, G4double theta1, G4double &E_out1, G4double &E_out2, G4double &theta2, G4double &q2, G4double &fWeightN, G4double &asymmetry)
static const G4int value_n
UI control for the event generator.
G4double Pion_PhotoProductionAl(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
G4double AlNuclInel(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
G4int SuperElasticCheck(G4double E_in, G4double E_out, G4double theta, G4double &xsec)
G4double GetPhiAngle_Min()
A multi-dimensional grid of values with interpolation methods.
G4double CGDR(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
void SetReactionRegion(G4int rr)
G4int GetActiveOctantNumber()
G4double Pion_PhotoProductionCarbon(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
G4double GetAsymmetry_Pi(G4double Q2_pi)
void SetPhiAngle_Max(G4double ang)
void F1F2QE09(G4int Z, G4int A, G4double QSQ, G4double wsq, G4double &F1, G4double &F2)
G4double Aluminum_Excited_State(G4double E_in, G4double Theta, G4double E_lvl, G4double fit_c, G4double fit_mean, G4double fit_sigma, G4double &fWeightN, G4double &Q2, G4double &E_out)
static const G4double TargetWindowThickness
G4double AlloyScattering(G4double E_in, G4double Theta, G4int Zin, G4int Ain, G4double &fWeightN, G4double &Q2, G4double &E_out)
void SetIsotropy(G4int isot)
G4double Horowitz_DW_Xsect(G4double E_in, G4double theta, const G4String path, G4double &fWeightN, G4double &Q2, G4double &E_out)
static const G4double M_p
G4double GetAsymmetry_AL(G4double theta, G4double energy)
void SetThetaAngle_Min(G4double ang)
G4double ElasticPeakDeltaE
G4double CNuclInel(G4double E_in, G4double Theta, G4int nState, G4double &fWeightN, G4double &Q2, G4double &E_out)
void SetSchwingerDeltaE(G4double energy=15 *MeV)
G4double GetVertexZ()
G4double Elastic_Cross_Section_Carbon(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
G4double ResMod507(G4int sf, G4double w2, G4double q2, G4double *xval)
void SetActiveOctantNumber(G4int kaot)
G4double GetBeamEnergy()
void SetElasticPeakDeltaE(G4double energy=15 *MeV)
G4double GetAsymmetry_EP(G4double theta, G4double energy)
G4double Quasi_Elastic_Neutron(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
static const G4double Theta_Min
void SetBeamEnergy(G4double energy=1.160 *GeV)
G4double AlGDR(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
G4double GetLuminosity()
static const G4double mil
QweakSimEPEvent(QweakSimUserInformation *myUI)
G4double Delta_Resonance(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
void SetThetaAngle_Max(G4double ang)
void F1F2IN09(G4int Z, G4int A, G4double QSQ, G4double wsq, G4double &F1, G4double &F2)
void christy507(G4double wsq, G4double Q2, G4double &F1, G4double &R, G4double &sigt, G4double &sigl)
G4double Elastic_Cross_Section_Proton(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
G4double Born_Approx_Nuclei_Asym(G4double Z, G4double A, G4double energy, G4double eprime, G4double qq)
void SetLuminosity(G4double lum)
G4double NuclearInelastic_Bosted(G4double E_in, G4double Theta, G4int Zin, G4int Ain, G4double &fWeightN, G4double &Q2, G4double &E_out)
static const G4double M_n
G4int Isotropy
isotropy used for event generation
G4int kActiveOctantNumber
Active octant number in the simulation, 0 will enable all octants.
G4double GetAsymmetry_Be(G4double theta, G4double energy)
G4double MEC2009(G4int a, G4double q2, G4double w2)
G4double GetPhaseSpace()
G4double GetAsymmetry_EN(G4double theta, G4double energy)
G4double GetElasticPeakDeltaE()
G4double Horowitz_DW_Asym(G4double theta, const G4String path)
G4double GetPhiAngle_Max()
static const G4int value_d
static const G4double TargetLength
virtual ~QweakSimEPEvent()
G4double Elastic_Cross_Section_Aluminum(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
G4double GetEPrime_Max()
G4double FF_BESSEL(G4int Z, G4int A, G4double q2, G4bool &ofr)
const std::vector< G4double > Radiative_Cross_Section_Lookup(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
G4double resmod507_v2(G4double sf, G4double w2, G4double q2, G4double xval[50])
QweakSimUserInformation * myUserInfo
G4double GetEPrime_Min()
G4double Sigma_EEPrime(G4double eni, G4double eprime, G4double theta, G4double &q2)
G4double Fgauss(G4int Z, G4int A, G4double q2)
G4double Pion_PhotoProduction(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
void resmodd(G4double w2, G4double q2, G4double xval[50], G4double &sig)
void SetPhiAngle_Min(G4double ang)