QwGeant4
QweakSimEPEvent.cc
Go to the documentation of this file.
1 //=============================================================================
2 //
3 // ---------------------------
4 // | Doxygen File Information |
5 // ---------------------------
6 //
7 /**
8 
9  \file QweakSimEPEvent.cc
10 
11  \author Jie Pan
12 
13 */
14 //=============================================================================
15 
16 #include "QweakSimEPEvent.hh"
17 
18 // geant4 includes
19 #include "Randomize.hh"
20 
21 // ROOT include?
22 #include "TGraph.h"
23 #include "TMath.h"
24 
25 // user includes
26 #include "QweakSimEPEventMessenger.hh"
27 #include "QweakSimUserInformation.hh"
28 #include "QweakSimPhysicalConstants.hh"
29 
30 #include "wiser_pion.h"
31 #include "Alinel_crsec_gen.h"
32 
33 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
34 
35  const G4double QweakSimEPEvent::M_n = 939.5656*MeV; //neutron mass in MeV/c^2
36  const G4double QweakSimEPEvent::M_p = 938.2796*MeV; // proton mass in MeV/c^2
37 
38  const G4double QweakSimEPEvent::Theta_Min = 1.745329E-4;
39 
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
41 
42 QweakSimEPEvent::QweakSimEPEvent( QweakSimUserInformation* myUI)
43 : fLookupTable(0)
44 {
45  G4cout << "###### Calling QweakSimEPEvent::QweakSimEPEvent () " << G4endl;
46 
47  Isotropy = 1;
48 
49  //SetPhiAngle_Min(-16.0*degree);
50  //SetPhiAngle_Max(16.0*degree);
51 
52  //SetThetaAngle_Min(4.0*degree);
53  //SetThetaAngle_Max(13.5*degree);
54 
55  //SetEPrime_Max(1.159*GeV);
56  //SetEPrime_Min(0.059*GeV);
57 
58  //SetBeamEnergy(1.16*GeV);
59 
60  ElasticPeakDeltaE = 15*MeV;
61  SchwingerDeltaE = 15*MeV;
62 
63  TypeSetting = 1;
64  ReactionType = 1;
65  ReactionRegion = 1;
66  kActiveOctantNumber = 0; //default octant all, choose from [0,8]
67 
68  myUserInfo = myUI;
69  myUserInfo->SetEPEvent(this);
70 
71  EventGen_Messenger = new QweakSimEPEventMessenger(this);
72 
73  G4cout << "###### Leaving QweakSimEPEvent::QweakSimEPEvent () " << G4endl;
74 }
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
77 
78 QweakSimEPEvent::~QweakSimEPEvent()
79 {
80 
81  G4cout << "###### Calling QweakSimEPEvent::~QweakSimEPEvent () " << G4endl;
82 
83  if( EventGen_Messenger) delete EventGen_Messenger;
84 
85  G4cout << "Created a great many events: " << myUserInfo->GetNumberOfEventToBeProcessed() << G4endl;
86 
87  G4cout << "###### Leaving QweakSimEPEvent::~QweakSimEPEvent () " << G4endl;
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91 
92 G4double QweakSimEPEvent::GetVertexZ()
93 {
94  /** \page reaction_region reaction region
95  * There are a number of pre-set reaction regions. These do not automatically set the
96  * reaction type, so don't expect aluminum exit window simulations by setting reaction
97  * region to 3...
98  * \li 1: LH2 Target
99  * \li 2: front entrance window
100  * \li 3: back exit window
101  * \li 4: Dummy Target (1% US Al)
102  * \li 5: Dummy Target (2% US Al)
103  * \li 6: Dummy Target (4% US Al)
104  * \li 7: Dummy Target (2% DS Al)
105  * \li 8: Dummy Target (4% DS Al)
106  * \li 9: Dummy Target (8% DS Al)
107  * \li 10: Dummy Target (US Carbon)
108  * \li 11: Dummy Target (DS Carbon)
109  * \li else: target
110  */
111 
112  if (ReactionRegion == 1) // target
113  myPositionZ = myUserInfo->TargetCenterPositionZ + (G4UniformRand()-0.5)*(myUserInfo->TargetLength);
114 
115  else if (ReactionRegion == 2) // front entrance window
116  myPositionZ = myUserInfo->TargetCenterPositionZ - 0.5*(myUserInfo->TargetLength)
117  - (myUserInfo->TargetEntranceWindowThickness)*G4UniformRand();
118 
119  else if(ReactionRegion == 3) // back exit window
120  myPositionZ = myUserInfo->TargetCenterPositionZ + 0.5*(myUserInfo->TargetLength)
121  + (myUserInfo->TargetExitWindowNippleThickness)*G4UniformRand();
122 
123  else if(ReactionType == 7 && (ReactionRegion == 4 || ReactionRegion == 5 || ReactionRegion == 6 || ReactionRegion == 10)) // Generator 7 US Dummy Targets
124  myPositionZ = myUserInfo->TargetCenterPositionZ + myUserInfo->TargetUSDummyPositionOffsetZ;
125 
126  else if(ReactionType == 7 && (ReactionRegion == 7 || ReactionRegion == 8 || ReactionRegion == 9 || ReactionRegion == 11)) // Generator 7 DS Dummy Targets
127  myPositionZ = myUserInfo->TargetCenterPositionZ + myUserInfo->TargetDSDummyPositionOffsetZ;
128 
129  else if(ReactionRegion == 4 || ReactionRegion == 5 || ReactionRegion == 6 || ReactionRegion == 10) // Al & C US Dummy Targets
130  myPositionZ = myUserInfo->TargetCenterPositionZ - 0.5*(myUserInfo->TargetLength)
131  - (myUserInfo->TargetEntranceWindowThickness)*G4UniformRand();
132 
133  else if(ReactionRegion == 7 || ReactionRegion == 8 || ReactionRegion == 9 || ReactionRegion == 11) // Al & C DS Dummy Targets
134  myPositionZ = myUserInfo->TargetCenterPositionZ + 0.5*(myUserInfo->TargetLength)
135  + (myUserInfo->TargetExitWindowNippleThickness)*G4UniformRand();
136 
137 
138 
139  else
140  myPositionZ = myUserInfo->TargetCenterPositionZ + (G4UniformRand()-0.5)*(myUserInfo->TargetLength); //default region
141 
142  return myPositionZ;
143 }
144 
145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
146 
147 G4ThreeVector QweakSimEPEvent::GetMomentumDirection()
148 {
149 
150  // Generate flat phi distribution
151  G4double PhiAngle = GetPhiAngle_Min() + G4UniformRand()*(GetPhiAngle_Max() - GetPhiAngle_Min());
152  // If active octant = 0, all octants are used
153  if (kActiveOctantNumber == 0) PhiAngle += (G4RandFlat::shootInt(8) * 45.0 * degree);
154  G4double cosPhi = cos(PhiAngle);
155  G4double sinPhi = sin(PhiAngle);
156 
157 
158  /** \page isotropy isotropy
159  * The cross section used for weighting Monte Carlo events depends on how those
160  * events are thrown. There are two options in the Qweak Monte Carlo:
161  * \li isotropy = 0: flat theta distribution, requires a modified weight for the cross section -> QweakSimUserPrimaryEvent::CrossSectionWeight
162  * \li isotropy = 1: uniform spherical theta/phi distribution, use regular cross section -> QweakSimUserPrimaryEvent::CrossSection
163  * \li else: not defined
164  */
165 
166  G4double cosTheta = 1.0;
167  G4double sinTheta = 0.0;
168  G4double ThetaAngle = 0.0;
169  G4double E_out = 0.0;
170  //G4double xsec = 0.0;
171  //G4int SuperElastic = 1;
172 
173  //while (SuperElastic) {
174  if (Isotropy == 0) {
175  // Generate flat theta distribution
176  ThetaAngle = GetThetaAngle_Min() + G4UniformRand()*(GetThetaAngle_Max() - GetThetaAngle_Min());
177  cosTheta = cos(ThetaAngle);
178  sinTheta = sin(ThetaAngle);
179 
180  } else if (Isotropy == 1) {
181  // Generate uniform distribution on spherical surface. See for example
182  // http://hypernews.slac.stanford.edu/HyperNews/geant4/get/particles/31/2.html
183  // or more generally http://mathworld.wolfram.com/SpherePointPicking.html
184  G4double cosThetaMax = cos(GetThetaAngle_Max());
185  G4double cosThetaMin = cos(GetThetaAngle_Min());
186  cosTheta = cosThetaMin + G4UniformRand()*(cosThetaMax - cosThetaMin);
187  sinTheta = sqrt(1. - cosTheta * cosTheta);
188  ThetaAngle = acos(cosTheta);
189 
190  } else {
191  G4cerr << "Warning: unkown isotropy type. Pick 0 or 1." << G4endl;
192  }
193  E_out = (G4UniformRand()*(GetEPrime_Max()-GetEPrime_Min()) + GetEPrime_Min());
194  //SuperElastic = SuperElasticCheck(myUserInfo->GetBeamEnergy(), E_out, ThetaAngle, xsec);
195  //}
196  myUserInfo->SetEPrime(E_out);
197 
198  G4double ux = sinTheta * cosPhi;
199  G4double uy = sinTheta * sinPhi;
200  G4double uz = cosTheta;
201  G4ThreeVector myNormMomentum(ux,uy,uz);
202 
203  // Rotate the momentum to the active octant (if octant = 0 all octants are used)
204  if (kActiveOctantNumber > 8 )
205  {
206  myNormMomentum.rotateZ( (kActiveOctantNumber-9 +4*G4RandFlat::shootInt(2))*45.0*degree);
207  }else if (kActiveOctantNumber > 0){
208  myNormMomentum.rotateZ( (kActiveOctantNumber-1)*45.0*degree);
209  }
210 
211  return myNormMomentum;
212 }
213 
214 
215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
216 
217 void QweakSimEPEvent::GetanEvent(G4double E_in,
218  std::vector< G4double > &fCrossSection,
219  G4double &fWeightN,
220  G4double &Q2,
221  G4double &E_out,
222  G4ThreeVector &OutgoingMomentumDirection,
223  G4double &theta,
224  G4double &phi,
225  G4double &Asymmetry)
226 {
227  // incoming momentum
228  G4ThreeVector IncomingMomentumDirection = myUserInfo->GetNormMomentum();
229 
230  G4double IncomingMomentumDirectionMag = IncomingMomentumDirection.mag();
231 
232 
233  // outgoing momentum
234  OutgoingMomentumDirection = GetMomentumDirection();
235  theta = OutgoingMomentumDirection.theta()/degree;
236  phi = OutgoingMomentumDirection.phi()/degree;
237  // relative theta angle
238 
239  G4double dot = OutgoingMomentumDirection.dot(IncomingMomentumDirection);
240 
241  G4double RelativeThetaAngle = acos(dot/IncomingMomentumDirectionMag);
242 
243 
244  if(ReactionType==0 || TypeSetting==0) // combination of all process, Random int [1,8]
245  {
246  TypeSetting = 0;
247  ReactionType = G4int(G4UniformRand()*7.0+1.0+0.5);
248  }
249  else
250  {
251  TypeSetting = ReactionType;
252  }
253 
254  //std::cout<<"ReactionType: "<<ReactionType<<", TypeSetting: "<<TypeSetting<<std::endl;
255 
256  /** \page reaction_type reaction type
257  * There are a number of reaction types implemented in the Qweak Monte Carlo:
258  * \li 0: random selection of some of these reaction types
259  * \li 1: elastic scattering from hydrogen (see QweakSimEPEvent::Elastic_Cross_Section_Proton)
260  * \li 2: elastic scattering from aluminum (see QweakSimEPEvent::Elastic_Cross_Section_Aluminum)
261  * \li 3: quasi-elastic scattering from proton in aluminum (see QweakSimEPEvent::Elastic_Cross_Section_Proton)
262  * \li 4: quasi-elastic scattering from neutron in aluminum (see QweakSimEPEvent::Quasi_Elastic_Neutron)
263  * \li 5: Delta resonance (see QweakSimEPEvent::Delta_Resonance)
264  * \li 6: Moller scattering (see QweakSimEPEvent::Moller_Scattering)
265  * \li 7: Radiative Scattering (see QweakSimEPEvent::Radiative_Cross_Section_Proton)
266  * \li 8: Quasielastic scattering from aluminum using the Bosted fit
267  * \li 88: Pion photoprodcution on LH2
268  * \li 89: Pion photoproduction on aluminum
269  * \li 90: Pion photoproduction on carbon
270  * \li 9: Aluminum Nuclear inelastic Single Particle States
271  * \li 10: Aluminum GDR
272  * \li 11: Aluminum Nuclear Radiated Inelastic from Bosted fit
273  * \li 12: Zn scattering approximation
274  * \li 13: Mg scattering approximation
275  * \li 14: Cu scattering approximation
276  * \li 15: Cr scattering approximation
277  * \li 16: Fe scattering approximation
278  * \li 17: Si scattering approximation
279  * \li 20: Elastic Scattering from Carbon (see QweakSimEPEvent::Elastic_Cross_Section_Carbon)
280  * \li 21: Alternate Elastic FF from FB fits (see QweakSimEPEvent::AlloyScattering) **USE THIS ONE**
281  * \li 22: Carbon GDR (see QweakSimEPEvent::CGDR)
282  * \li 23: Carbon-12 4.439 MeV, 2+ excited state scattering
283  * \li 24: Carbon-12 7.654 MeV, 0+ excited state scattering
284  * \li 25: Carbon-12 9.641 MeV, 3- excited state scattering
285  * \li 100: Isotropic and with Delta E in a uniform Distribution over 500 MeV
286  * \li 2700: Horowitz Distorted Wave Al Elastic Cross-section/Asym Calculation
287  * \li 2701: Aluminum 0.84376 MeV Jp=1/2+ dT=0 Excited State Scattering
288  * \li 2702: Aluminum 1.01445 MeV Jp=3/2+ dT=0 Excited State Scattering
289  * \li 2703: Aluminum 2.2111 MeV Jp=7/2+ dT=0 Excited State Scattering
290  * \li 2704: Aluminum 2.7349 MeV Jp=5/2+ dT=0 Excited State Scattering
291  * \li 2705: Aluminum 2.9931 MeV Jp=(3/2+)/(9/2+) dT=0 doublet Excited State Scattering
292  * \li 2706: Aluminum 4.580 MeV Jp=7/2+ dT=0 Excited State Scattering
293  * \li 2707: Aluminum 4.8116 MeV Jp=5/2+ dT=0 Excited State Scattering
294  * \li 2708: Aluminum 5.4328 MeV Jp=7/2 dT=0 Excited State Scattering
295  * \li 2709: Aluminum 5.6673 MeV Jp=9/2+ dT=0 Excited State Scattering
296  * \li 2710: Aluminum 7.2272 MeV Jp=9/2- dT=0? Excited State Scattering
297  * \li 2711: Aluminum 7.4771 MeV Jp=7/2- dT=0? Excited State Scattering
298  */
299 
300  if(ReactionType==1) //LH2 target
301  {
302  //G4double A = 1.0;
303  //G4double Z = 1.0;
304  //G4double Mass = Z*M_p+(A-Z)*M_n;
305  fCrossSection[0] = Elastic_Cross_Section_Proton(E_in, RelativeThetaAngle, fWeightN, Q2, E_out);
306  Asymmetry = GetAsymmetry_EP(RelativeThetaAngle, E_in);
307 
308  }
309 
310  else if(ReactionType==2) //Elastic Aluminum windows/dummy
311  {
312  G4double Z = 13.0;
313  G4double A = 27.0;
314  //G4double Mass = Z*M_p+(A-Z)*M_n;
315  fCrossSection[0] = Elastic_Cross_Section_Aluminum(E_in, RelativeThetaAngle, fWeightN, Q2, E_out);
316  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
317  //Asymmetry = GetAsymmetry_AL(RelativeThetaAngle, E_in);
318  }
319 
320  else if(ReactionType==3) // Aluminum window quasi-elastic proton (assume free proton)
321  {
322  //G4double A = 1.0;
323  //G4double Z = 1.0;
324  //G4double Mass = M_p;
325  fCrossSection[0] = Elastic_Cross_Section_Proton(E_in, RelativeThetaAngle, fWeightN, Q2, E_out);
326  Asymmetry = GetAsymmetry_EP(RelativeThetaAngle, E_in);
327  }
328 
329  else if(ReactionType==4) // Aluminum window quasi-elastic neutron (assume free neutron)
330  {
331  //G4double A = 1.0;
332  //G4double Z = 1.0; // Z needs to be set to 1 for neutron quasi elastic scattering
333  //G4double Mass = M_n;
334  fCrossSection[0] = Quasi_Elastic_Neutron(E_in, RelativeThetaAngle, fWeightN, Q2, E_out);
335  Asymmetry = GetAsymmetry_EN(RelativeThetaAngle, E_in);
336  }
337 
338  else if(ReactionType==5) // Delta resonances
339  {
340  fCrossSection[0] = Delta_Resonance(E_in,RelativeThetaAngle, fWeightN, Q2, E_out);
341  //std::cout<<E_in<<" "<<ThetaAngle/degree<<" "<<fWeightN<<" "<<Q2<<" "<<E_out<<std::endl;
342  Asymmetry = GetAsymmetry_Pi(Q2);
343  }
344  else if(ReactionType==6) // Moller scattering
345  {
346  // Peiqing, Dec 12, 2011
347  // Small angle recoil electrons are directly dumped for now.
348  G4double E_recoil;
349  G4double ThetaRecoil;
350  fCrossSection[0] = Moller_Scattering(E_in, RelativeThetaAngle, E_out,
351  E_recoil, ThetaRecoil,
352  Q2, fWeightN, Asymmetry);
353  }
354  else if(ReactionType==7) // Radiative Cross Section Lookup Table
355  {
356  fCrossSection = Radiative_Cross_Section_Lookup(myUserInfo->GetBeamEnergy(), RelativeThetaAngle, fWeightN, Q2, E_out);
357  Asymmetry = GetAsymmetry_EP(RelativeThetaAngle, E_in);
358  }
359  else if(ReactionType==8) // Quasielastic Bosted - Aluminum
360  {
361  G4int Z = 13;
362  G4int A = 27;
363  fCrossSection[0] = Quasi_Elastic_Bosted(E_in, RelativeThetaAngle, Z, A, fWeightN, Q2, E_out);
364  }
365  else if(ReactionType==88) //--- LH2 target, pion photo-production cross section 3.35 GeV
366  {
367  fCrossSection[0] = Pion_PhotoProduction(E_in, RelativeThetaAngle, fWeightN, Q2, E_out);
368  Asymmetry = GetAsymmetry_EP(RelativeThetaAngle, E_in); //--- use the elastic asymmetry for now
369  }
370  else if(ReactionType==89) //--- Al target, pion photo-production cross section 3.35 GeV
371  {
372  fCrossSection[0] = Pion_PhotoProductionAl(E_in, RelativeThetaAngle, fWeightN, Q2, E_out);
373  Asymmetry = GetAsymmetry_EP(RelativeThetaAngle, E_in); //--- use the elastic asymmetry for now
374  }
375  else if(ReactionType==90) //--- Carbon target, pion photo-production cross section 3.35 GeV
376  {
377  fCrossSection[0] = Pion_PhotoProductionCarbon(E_in, RelativeThetaAngle, fWeightN, Q2, E_out);
378  Asymmetry = GetAsymmetry_EP(RelativeThetaAngle, E_in); //--- use the elastic asymmetry for now
379  }
380  else if(ReactionType==9) // Nuclear Inelastics Single Particle States - Aluminum
381  {
382  fCrossSection[0] = AlNuclInel(E_in, RelativeThetaAngle, fWeightN, Q2, E_out);
383  }
384  else if(ReactionType==10) // GDR - Aluminum
385  {
386  G4double Z = 13.;
387  G4double A = 27.;
388  fCrossSection[0] = AlGDR(E_in, RelativeThetaAngle, fWeightN, Q2, E_out);
389  Asymmetry = -1.0*Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
390  // Factor of -1.0, for GDR being isovector 2017/8/25 - KDB
391  }
392  else if(ReactionType==11) // Nuclear Inelastic Bosted - Aluminum
393  {
394  G4double Z = 13;
395  G4double A = 27;
396  fCrossSection[0] = NuclearInelastic_Bosted(E_in, RelativeThetaAngle, Z, A, fWeightN, Q2, E_out);
397  }
398  else if(ReactionType==12) // Alloy: Zinc elastic scattering approximation
399  {
400  G4double Z = 30;
401  G4double A = 64;
402  fCrossSection[0] = AlloyScattering(E_in, RelativeThetaAngle, Z, A, fWeightN, Q2, E_out);
403  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
404  }
405  else if(ReactionType==13) // Alloy: Magnesium elastic scattering approximation
406  {
407  G4double Z = 12;
408  //G4double A = 24; //24Mg currently used Gaussian approx - 2016/10/25 KDB
409  G4double A = 26; //Will use FB coefficients
410  fCrossSection[0] = AlloyScattering(E_in, RelativeThetaAngle, Z, A, fWeightN, Q2, E_out);
411  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
412  }
413  else if(ReactionType==14) // Alloy: Copper elastic scattering approximation
414  {
415  G4double Z = 29;
416  G4double A = 63;//Use 63Cu, largest abundance of stable isotope.
417  fCrossSection[0] = AlloyScattering(E_in, RelativeThetaAngle, Z, A, fWeightN, Q2, E_out);
418  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
419  }
420  else if(ReactionType==15) // Alloy: Chromium elastic scattering approximation
421  {
422  G4double Z = 24;
423  G4double A = 52;
424  fCrossSection[0] = AlloyScattering(E_in, RelativeThetaAngle, Z, A, fWeightN, Q2, E_out);
425  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
426  }
427  else if(ReactionType==16) // Alloy: Iron elastic scattering approximation
428  {
429  G4double Z = 26;
430  G4double A = 56;
431  fCrossSection[0] = AlloyScattering(E_in, RelativeThetaAngle, Z, A, fWeightN, Q2, E_out);
432  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
433  }
434  else if(ReactionType==17) // Alloy: Silicon elastic scattering approximation
435  {
436  G4double Z = 14;
437  G4double A = 28;
438  fCrossSection[0] = AlloyScattering(E_in, RelativeThetaAngle, Z, A, fWeightN, Q2, E_out);
439  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
440  }
441  else if(ReactionType==18) // Alloy: Manganese elastic scattering approximation
442  {
443  G4double Z = 25;
444  G4double A = 55;
445  fCrossSection[0] = AlloyScattering(E_in, RelativeThetaAngle, Z, A, fWeightN, Q2, E_out);
446  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
447  }
448  else if(ReactionType==19) // Alloy: Titanium elastic scattering approximation
449  {
450  G4double Z = 22;
451  G4double A = 48;
452  fCrossSection[0] = AlloyScattering(E_in, RelativeThetaAngle, Z, A, fWeightN, Q2, E_out);
453  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
454  }
455  else if(ReactionType==20) // Elastic C
456  {
457  fCrossSection[0] = Elastic_Cross_Section_Carbon(E_in, RelativeThetaAngle, fWeightN, Q2, E_out);
458  }
459  else if(ReactionType==21) // Elastic C Alternate
460  {
461  G4double Z = 6;
462  G4double A = 12;
463  fCrossSection[0] = AlloyScattering(E_in, RelativeThetaAngle, Z, A, fWeightN, Q2, E_out);
464  }
465  else if(ReactionType==22) // GDR - Carbon
466  {
467  fCrossSection[0] = CGDR(E_in, RelativeThetaAngle, fWeightN, Q2, E_out);
468  }
469  else if(ReactionType==23) // Inelastic Nuclear Excited States - Carbon-12 4.439 MeV, 2+
470  {
471  fCrossSection[0] = CNuclInel(E_in, RelativeThetaAngle, 1, fWeightN, Q2, E_out);
472  }
473  else if(ReactionType==24) // Inelastic Nuclear Excited States - Carbon-12 7.654 MeV, 0+
474  {
475  fCrossSection[0] = CNuclInel(E_in, RelativeThetaAngle, 2, fWeightN, Q2, E_out);
476  }
477  else if(ReactionType==25) // Inelastic Nuclear Excited States - Carbon-12 9.641 MeV, 3-
478  {
479  fCrossSection[0] = CNuclInel(E_in, RelativeThetaAngle, 3, fWeightN, Q2, E_out);
480  }
481  else if(ReactionType==100)
482  {
483  fCrossSection[0] = 0;
484  E_out = E_in*G4UniformRand();
485  if(E_out < 0.0) E_out = E_in;
486  G4double STH = sin(RelativeThetaAngle/2.0);
487  Q2 = 4.0*E_in*E_out*STH*STH;
488 
489  // G4cout << E_in << "\t" << RelativeThetaAngle << "\t" << E_out << "\t" << Q2 << G4endl;
490  }
491  else if(ReactionType==2700) // Horowitz Distorted Wave Al Elastic Cross-section/Asym Calculation
492  {
493  fCrossSection[0] = Horowitz_DW_Xsect(E_in,
494  RelativeThetaAngle,
495  "horowitz_alloy_tables/elasticAl27.csv",
496  fWeightN,
497  Q2,
498  E_out);
499  Asymmetry = Horowitz_DW_Asym(RelativeThetaAngle,
500  "horowitz_alloy_tables/elasticAl27.csv");
501  }
502  else if(ReactionType==2701) // Al Nuclear Excited State 0.84376 MeV Jp = 1/2+, dT=0
503  {
504  G4double Z = 13.0; // # of protons
505  G4double A = 27.0; // Atomic Mass [amu]
506  fCrossSection[0] = Aluminum_Excited_State(E_in,
507  RelativeThetaAngle,
508  0.84376, //[MeV]
509  0.000478, //[unitless]
510  0.9599, //[fm^-1]
511  0.3072, //[fm^-1]
512  fWeightN,
513  Q2,
514  E_out);
515  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
516  }
517  else if(ReactionType==2702) // Al Nuclear Excited State 1.01445 MeV Jp = 3/2+, dT=0
518  {
519  G4double Z = 13.0; // # of protons
520  G4double A = 27.0; // Atomic Mass [amu]
521  fCrossSection[0] = Aluminum_Excited_State(E_in,
522  RelativeThetaAngle,
523  1.01445, //[MeV]
524  0.000719, //[unitless]
525  0.9049, //[fm^-1]
526  0.2815, //[fm^-1]
527  fWeightN,
528  Q2,
529  E_out);
530  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
531  }
532  else if(ReactionType==2703) // Al Nuclear Excited State 2.2111 MeV Jp = 7/2+, dT=0
533  {
534  G4double Z = 13.0; // # of protons
535  G4double A = 27.0; // Atomic Mass [amu]
536  fCrossSection[0] = Aluminum_Excited_State(E_in,
537  RelativeThetaAngle,
538  2.2111, //[MeV]
539  0.002351, //[unitless]
540  0.8895, //[fm^-1]
541  0.2647, //[fm^-1]
542  fWeightN,
543  Q2,
544  E_out);
545  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
546  }
547  else if(ReactionType==2704) // Al Nuclear Excited State 2.7349 MeV Jp = 5/2+, dT=0
548  {
549  G4double Z = 13.0; // # of protons
550  G4double A = 27.0; // Atomic Mass [amu]
551  fCrossSection[0] = Aluminum_Excited_State(E_in,
552  RelativeThetaAngle,
553  2.7349, //[MeV]
554  0.000318, //[unitless]
555  0.9615, //[fm^-1]
556  0.4139, //[fm^-1]
557  fWeightN,
558  Q2,
559  E_out);
560  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
561  }
562  else if(ReactionType==2705) // Al Nuclear Excited State 2.9931 MeV
563  //(2.98200 MeV + 3.0042 MeV doublet) Jp = (3/2+)/(9/2+), dT=0
564  {
565  G4double Z = 13.0; // # of protons
566  G4double A = 27.0; // Atomic Mass [amu]
567  fCrossSection[0] = Aluminum_Excited_State(E_in,
568  RelativeThetaAngle,
569  2.9931, //[MeV]
570  0.001665, //[unitless]
571  0.9577, //[fm^-1]
572  0.2999, //[fm^-1]
573  fWeightN,
574  Q2,
575  E_out);
576  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
577  }
578  else if(ReactionType==2706) // Al Nuclear Excited State 4.580 MeV Jp = 7/2+, dT=0
579  {
580  G4double Z = 13.0; // # of protons
581  G4double A = 27.0; // Atomic Mass [amu]
582  fCrossSection[0] = Aluminum_Excited_State(E_in,
583  RelativeThetaAngle,
584  4.580, //[MeV]
585  0.000159, //[unitless]
586  1.3896, //[fm^-1]
587  0.4825, //[fm^-1]
588  fWeightN,
589  Q2,
590  E_out);
591  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
592  }
593  else if(ReactionType==2707) // Al Nuclear Excited State 4.8116 MeV Jp = 5/2+, dT=0
594  {
595  G4double Z = 13.0; // # of protons
596  G4double A = 27.0; // Atomic Mass [amu]
597  fCrossSection[0] = Aluminum_Excited_State(E_in,
598  RelativeThetaAngle,
599  4.8116, //[MeV]
600  0.000230, //[unitless]
601  0.2536, //[fm^-1]
602  0.5944, //[fm^-1]
603  fWeightN,
604  Q2,
605  E_out);
606  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
607  }
608  else if(ReactionType==2708) // Al Nuclear Excited State 5.4328 MeV Jp = 7/2, dT=0
609  {
610  G4double Z = 13.0; // # of protons
611  G4double A = 27.0; // Atomic Mass [amu]
612  fCrossSection[0] = Aluminum_Excited_State(E_in,
613  RelativeThetaAngle,
614  5.4328, //[MeV]
615  0.000269, //[unitless]
616  0.7793, //[fm^-1]
617  0.4466, //[fm^-1]
618  fWeightN,
619  Q2,
620  E_out);
621  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
622  }
623  else if(ReactionType==2709) // Al Nuclear Excited State 5.6673 MeV Jp = 9/2+, dT=0
624  {
625  G4double Z = 13.0; // # of protons
626  G4double A = 27.0; // Atomic Mass [amu]
627  fCrossSection[0] = Aluminum_Excited_State(E_in,
628  RelativeThetaAngle,
629  5.6673, //[MeV]
630  0.000137, //[unitless]
631  0.8912, //[fm^-1]
632  0.3636, //[fm^-1]
633  fWeightN,
634  Q2,
635  E_out);
636  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
637  }
638  else if(ReactionType==2710) // Al Nuclear Excited State 7.2272 MeV Jp = 9/2-, dT=0?
639  {
640  G4double Z = 13.0; // # of protons
641  G4double A = 27.0; // Atomic Mass [amu]
642  fCrossSection[0] = Aluminum_Excited_State(E_in,
643  RelativeThetaAngle,
644  7.2272, //[MeV]
645  0.000356, //[unitless]
646  1.0592, //[fm^-1]
647  0.3218, //[fm^-1]
648  fWeightN,
649  Q2,
650  E_out);
651  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
652  }
653  else if(ReactionType==2711) // Al Nuclear Excited State 7.4771 MeV Jp = 1/2+, dT=0?
654  {
655  G4double Z = 13.0; // # of protons
656  G4double A = 27.0; // Atomic Mass [amu]
657  fCrossSection[0] = Aluminum_Excited_State(E_in,
658  RelativeThetaAngle,
659  7.4771, //[MeV]
660  0.000210, //[unitless]
661  1.0785, //[fm^-1]
662  0.2994, //[fm^-1]
663  fWeightN,
664  Q2,
665  E_out);
666  Asymmetry = Born_Approx_Nuclei_Asym(Z, A, E_in, E_out, Q2);//[ppm]
667  }
668  else if(ReactionType==2750) // Horowitz Distorted Wave Zn Elastic Cross-section/Asym Calculation
669  {
670  fCrossSection[0] = Horowitz_DW_Xsect(E_in,
671  RelativeThetaAngle,
672  "horowitz_alloy_tables/elasticZn64.csv",
673  fWeightN,
674  Q2,
675  E_out);
676  Asymmetry = Horowitz_DW_Asym(RelativeThetaAngle,
677  "horowitz_alloy_tables/elasticZn64.csv");
678  }
679  else if(ReactionType==2751) // Horowitz Distorted Wave Mg Elastic Cross-section/Asym Calculation
680  {
681  fCrossSection[0] = Horowitz_DW_Xsect(E_in,
682  RelativeThetaAngle,
683  "horowitz_alloy_tables/elasticMg24.csv",
684  fWeightN,
685  Q2,
686  E_out);
687  Asymmetry = Horowitz_DW_Asym(RelativeThetaAngle,
688  "horowitz_alloy_tables/elasticMg24.csv");
689  }
690  else if(ReactionType==2752) // Horowitz Distorted Wave Cu Elastic Cross-section/Asym Calculation
691  {
692  fCrossSection[0] = Horowitz_DW_Xsect(E_in,
693  RelativeThetaAngle,
694  "horowitz_alloy_tables/elasticCu63.csv",
695  fWeightN,
696  Q2,
697  E_out);
698  Asymmetry = Horowitz_DW_Asym(RelativeThetaAngle,
699  "horowitz_alloy_tables/elasticCu63.csv");
700  }
701  else if(ReactionType==2753) // Horowitz Distorted Wave Cr Elastic Cross-section/Asym Calculation
702  {
703  fCrossSection[0] = Horowitz_DW_Xsect(E_in,
704  RelativeThetaAngle,
705  "horowitz_alloy_tables/elasticCr52.csv",
706  fWeightN,
707  Q2,
708  E_out);
709  Asymmetry = Horowitz_DW_Asym(RelativeThetaAngle,
710  "horowitz_alloy_tables/elasticCr52.csv");
711  }
712  else if(ReactionType==2754) // Horowitz Distorted Wave Fe Elastic Cross-section/Asym Calculation
713  {
714  fCrossSection[0] = Horowitz_DW_Xsect(E_in,
715  RelativeThetaAngle,
716  "horowitz_alloy_tables/elasticFe56.csv",
717  fWeightN,
718  Q2,
719  E_out);
720  Asymmetry = Horowitz_DW_Asym(RelativeThetaAngle,
721  "horowitz_alloy_tables/elasticFe56.csv");
722  }
723  else if(ReactionType==2755) // Horowitz Distorted Wave Si Elastic Cross-section/Asym Calculation
724  {
725  fCrossSection[0] = Horowitz_DW_Xsect(E_in,
726  RelativeThetaAngle,
727  "horowitz_alloy_tables/elasticSi28.csv",
728  fWeightN,
729  Q2,
730  E_out);
731  Asymmetry = Horowitz_DW_Asym(RelativeThetaAngle,
732  "horowitz_alloy_tables/elasticSi28.csv");
733  }
734  else
735  {
736  G4cout << "Warning! ReactionType Not SELECTED!!!!" << G4endl;
737  }
738 
739  if (ReactionRegion == 1) SetLuminosity(myUserInfo->TargetLuminosityLH2);
740  else if (ReactionRegion == 2) SetLuminosity(myUserInfo->TargetLuminosityUSALWindow);
741  else if (ReactionRegion == 3) SetLuminosity(myUserInfo->TargetLuminosityDSALWindow);
742  else if (ReactionRegion == 4) SetLuminosity(myUserInfo->TargetLuminosityUSALDummy1);
743  else if (ReactionRegion == 5) SetLuminosity(myUserInfo->TargetLuminosityUSALDummy2);
744  else if (ReactionRegion == 6) SetLuminosity(myUserInfo->TargetLuminosityUSALDummy4);
745  else if (ReactionRegion == 7) SetLuminosity(myUserInfo->TargetLuminosityDSALDummy2);
746  else if (ReactionRegion == 8) SetLuminosity(myUserInfo->TargetLuminosityDSALDummy4);
747  else if (ReactionRegion == 9) SetLuminosity(myUserInfo->TargetLuminosityDSALDummy8);
748  else if (ReactionRegion == 10) SetLuminosity(myUserInfo->TargetLuminosityUSCDummy);
749  else if (ReactionRegion == 11) SetLuminosity(myUserInfo->TargetLuminosityDSCDummy);
750  else SetLuminosity(myUserInfo->TargetLuminosityLH2);
751 
752  // Allow calculation of the luminosity that includes fractional weights for
753  // other elements in the Aluminum alloy for DS 4% dummy. K. Bartlett 2016/9/13
754  if (ReactionRegion == 8 && (ReactionType==2 || ReactionType==7 || ReactionType==8 ||
755  ReactionType==9 || ReactionType==10 || ReactionType==11 || ReactionType==89 ||
756  ReactionType==2700 || ReactionType==2701 || ReactionType==2702 || ReactionType==2703 ||
757  ReactionType==2704 || ReactionType==2705 || ReactionType==2706 || ReactionType==2707 ||
758  ReactionType==2708 || ReactionType==2709 || ReactionType==2710 || ReactionType==2711)){
759  SetLuminosity(myUserInfo->TargetLuminosityDSALDummy4_Al);
760  }
761  else if (ReactionRegion == 8 && (ReactionType==12 || ReactionType==2750)) SetLuminosity(myUserInfo->TargetLuminosityDSALDummy4_Zn);
762  else if (ReactionRegion == 8 && (ReactionType==13 || ReactionType==2751)) SetLuminosity(myUserInfo->TargetLuminosityDSALDummy4_Mg);
763  else if (ReactionRegion == 8 && (ReactionType==14 || ReactionType==2752)) SetLuminosity(myUserInfo->TargetLuminosityDSALDummy4_Cu);
764  else if (ReactionRegion == 8 && (ReactionType==15 || ReactionType==2753)) SetLuminosity(myUserInfo->TargetLuminosityDSALDummy4_Cr);
765  else if (ReactionRegion == 8 && (ReactionType==16 || ReactionType==2754)) SetLuminosity(myUserInfo->TargetLuminosityDSALDummy4_Fe);
766  else if (ReactionRegion == 8 && (ReactionType==17 || ReactionType==2755)) SetLuminosity(myUserInfo->TargetLuminosityDSALDummy4_Si);
767  else if (ReactionRegion == 8 && ReactionType==18) SetLuminosity(myUserInfo->TargetLuminosityDSALDummy4_Mn);
768  else if (ReactionRegion == 8 && ReactionType==19) SetLuminosity(myUserInfo->TargetLuminosityDSALDummy4_Ti);
769 
770 
771  //G4double phase_space = (GetPhiAngle_Max()-GetPhiAngle_Min())*(cos(GetThetaAngle_Min())-cos(GetThetaAngle_Max())); // units [Sr]
772 
773  G4double phase_space = (GetPhiAngle_Max()-GetPhiAngle_Min());
774 
775  if (GetIsotropy() == 0) phase_space *= (GetThetaAngle_Max()-GetThetaAngle_Min());// * sin(RelativeThetaAngle);
776  else if (GetIsotropy() == 1) phase_space *= (cos(GetThetaAngle_Min())-cos(GetThetaAngle_Max()));
777  else phase_space *= (cos(GetThetaAngle_Min())-cos(GetThetaAngle_Max()));
778 
779  // For ReactionTypes that are differential in E'
780  if (ReactionType == 7 || ReactionType == 88 || ReactionType == 89 || ReactionType == 90) phase_space *= (GetEPrime_Max()-GetEPrime_Min())/1000.0; // units [Sr*GeV]
781  // Delta_Resonance (reac 5) and Quasi_Elastic_Bosted (reac 8) are also
782  // differential in Eprime, but the xsect returned by these functions
783  // already includes the additional phase space factor. - K.Mesick 13/May/14
784 
785 
786  // If events are being thrown in all octants, this increases the phase space accordingly
787  if (!GetActiveOctantNumber()) phase_space *= 8.0;
788 
789 
790  SetPhaseSpace(phase_space);
791 
792 }
793 
794 
795 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
796 // jpan@nuclear.uwinnipeg.ca
797 // calculate proton cross sections using the dipole fits to the form factors.
798 // Angles are restricted to be greater than .01 deg to avoid
799 // division by 0 when evaluating the Mott cross section.
800 
801 G4double QweakSimEPEvent::Elastic_Cross_Section_Proton(G4double E_in,
802  G4double Theta,
803  G4double &fWeightN,
804  G4double &Q2,
805  G4double &E_out)
806 {
807  G4double Lamda_2 = 0.710;
808  G4double M_p = 938.2796 * MeV; // proton mass in MeV
809  G4double mu = 2.793;
810  G4double Z = 1.0;
811  G4double A = 1.0;
812  G4double M = M_p*A;
813  G4double myhbarc = hbarc / MeV / fermi; // 197.3269718 in MeV fm
814  G4double alpha = 1.0/137.035999074;
815  G4double CC = myhbarc*alpha/2.0; // 0.719982242379
816  G4double Electron_Mass = 0.511 * MeV; //electron mass in MeV
817 
818 // E_in units is MeV
819 
820  const G4double theta_min = 0.01 * degree;
821  if (Theta < theta_min) {
822  Theta = theta_min;
823  G4cout << "Warning: Elastic_Cross_Section_Proton: theta less than " << theta_min << G4endl;
824  G4cout << "Warning: Elastic_Cross_Section_Proton: theta was set to " << theta_min << G4endl;
825  }
826 
827  G4double CTH = cos(Theta/2.);
828  G4double STH = sin(Theta/2.);
829  G4double T2THE = STH*STH/CTH/CTH;
830  G4double ETA = 1.0+2.0*E_in*STH*STH/M;
831  E_out = E_in/ETA;
832  Q2 = 4.0*E_in*E_out*STH*STH; // MeV^2
833  G4double tau = Q2/4.0/M/M;
834 // Mott scatering cross-section, including recoil correction
835  G4double CrossSection = (Z*CC/E_in*CTH/STH/STH)*(Z*CC/E_in*CTH/STH/STH)/ETA;
836 // Units: ub/sr
837  G4double Mott = CrossSection*10000.0;
838 // Cross section
839  G4double GEP_DIPOLE = 1.0/(1.0+Q2/1.E6/Lamda_2)/(1.0+Q2/1.E6/Lamda_2);
840  G4double GMP_DIPOLE = GEP_DIPOLE*mu;
841  G4double FAC = 1.0/(1.0+tau);
842 
843  //The next two line is to add schwinger
844  G4double FunctionofTheta = log (STH*STH) * log (CTH*CTH);
845  G4double delta_Schwinger = (-2.0*alpha/pi) * ((log(E_in/15.0) - 13.0/12.0) * (log(Q2/(Electron_Mass*Electron_Mass)) - 1.0) + 17.0/36.0 + FunctionofTheta/2.0);
846 // G4double omega_Sch = E_in - 15;
847  G4double Sigma_Dipole;
848  Sigma_Dipole = Mott*(GEP_DIPOLE*GEP_DIPOLE*FAC+tau*GMP_DIPOLE*GMP_DIPOLE*(FAC+2.*T2THE));
849  Sigma_Dipole *= (1.0 + delta_Schwinger);
850  fWeightN = Sigma_Dipole*sin(Theta);
851  return Sigma_Dipole;
852 }
853 
854 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
855 // jpan@nuclear.uwinnipeg.ca
856 // calculate Aluminum elastic cross sections
857 
858 G4double QweakSimEPEvent::Elastic_Cross_Section_Aluminum(G4double E_in,
859  G4double Theta,
860  G4double &fWeightN,
861  G4double &Q2,
862  G4double &E_out)
863 {
864  G4double Theta_Min = 1.745329E-4;
865  G4double M_p = 938.2796; // proton mass in MeV
866  G4double Z = 13.0;
867  G4double A = 27.0;
868  G4double M = M_p*A;
869 
870 // harmonic oscillator well parameter a0 ~1.76 fm
871  G4double a = 2.98; //unit: fm
872  G4double ap = sqrt(0.427); //unit: fm
873  G4double a0 = sqrt((a*a-1.5*ap*ap)/(3.5-10/Z-1.5/A));
874 
875  G4double Q = 14.6; //unit fm^(-2)
876 // G4double J = 5.0/2.0;
877 
878 // G4double gamma = 2.792847351; //This is the magnetic moment
879 // G4double Omega = -3/2*(1+2*gamma)*a0*a0;
880 // G4double Gamma = 15/4*gamma*pow(a0,4);
881 // G4double mu = 3.69; // This is NOT the magnetic moment
882 // G4double Omega_mu = -6.36; // Omega/mu
883 // G4double Gamma_mu = 20.7; // Gamma/mu
884 
885  if (Theta<Theta_Min)
886  Theta = Theta_Min;
887 
888 // E_in unit is MeV, q2 unit is fm^(-2)
889  G4double CTH = cos(Theta/2.0);
890  G4double STH = sin(Theta/2.0);
891 // G4double T2THE = STH*STH/CTH/CTH;
892  G4double ETA = 1.0+2.0*E_in*STH*STH/M;
893  E_out = E_in/ETA;
894 
895  Q2 = 4*E_in*E_out*STH*STH; //unit: MeV^2
896  G4double q2 = Q2/1000000*(1.0/0.197)*(1.0/0.197); //convert MeV^2 into fm^(-2)
897  G4double x = (1.0/4.0)*q2*a0*a0;
898 
899 // std::ofstream EventDataFile("Event.dat", std::ios::app);
900 // EventDataFile<<"E_f="<<E_f<<" q2="<<q2<<" "<<"x="<<x<<G4endl;
901 
902 // Electric form factor (data fit)
903  G4double F0 = (1.0/Z)*( Z-4.0/3.0*(Z-5.0)*x+4.0/15.0*(Z-8.0)*x*x)*exp(-x);
904  G4double F2 = (1.0-2.0/7.0*x)*exp(-x);
905  G4double Fe = sqrt( F0*F0+(7.0/450.0)*q2*q2*(Q*Q/Z/Z)*F2*F2 );
906  Fe=Fe*exp(-0.25*q2*ap*ap); //correction for finite proton size
907  Fe=Fe*exp(x/A); //correction for center-of-well motion
908  G4double Fe_2 = pow(Fe,2);
909 
910 
911 // EventDataFile<<"F0="<<F0<<" "<<"F2="<<F2<<" "<<"Fe_2="<<Fe_2<<G4endl;
912 
913 // Magnetic form factor (theoretical)
914 // G4double Fm1 = (1.0-2.0/5.0*(1.0+2.0*gamma)/(1.0+gamma/2.0)*x+6.0/35.0*gamma*x*x/(1.0+gamma/2.0))*exp(-x);
915 // G4double Fm3 = (1.0-2.0/3.0*gamma*x/(1.0+2.0*gamma))*exp(-x);
916 // G4double Fm5 = exp(-x);
917 // G4double Fm_2 = Fm1*Fm1+(4.0/525.0)*pow(q2*Omega_mu*Fm3,2)+(2.0/33075.0)*pow(q2,4)*pow(Gamma_mu*Fm5,2);
918 
919 // EventDataFile<<"Fm1="<<Fm1<<" "<<"Fm3="<<Fm3<<" "<<"Fm5="<<Fm5<<" "<<"Fm_2="<<Fm_2<<G4endl;
920 
921 // form factor square
922 // G4double F_2 = Z*Z*Fe_2+(1.0+2.0*T2THE)*mu*mu/3.0*(J+1.0)/J*q2/4.0/pow((M_p/0.197),2)*Fm_2;
923  G4double F_2 = Fe_2;//+(1.0+2.0*T2THE)*mu*mu/3.0*(J+1.0)/J*q2/4.0/pow((M_p/0.197),2)*Fm_2;
924 
925 // EventDataFile<<"F_2="<<F_2<<" "<<G4endl;
926 
927  //The next lines are to add schwinger
928  G4double Electron_Mass = 0.511 * MeV; //electorn mass in MeV
929  G4double alpha = 1.0/137.035999074;
930 
931  G4double FunctionofTheta = log (STH*STH) * log (CTH*CTH);
932  G4double delta_Schwinger = (-2.0*alpha/pi) * ((log(E_in/15.0) - 13.0/12.0) * (log(Q2/(Electron_Mass*Electron_Mass)) - 1.0) + 17.0/36.0 + FunctionofTheta/2.0);
933 
934 // cross section units: ub/sr
935 // G4double SigmaMott = ((0.72/E_in)*cos(Theta/2)/(sin(Theta/2)^2))^2 /( 1+2*E_in/M*(sin(Theta/2)^2))*10000 ;
936  G4double SigmaMott = pow(((0.72/E_in)*CTH/(STH*STH)),2)/(1+2*E_in/M*STH*STH)*10000 ;
937  SigmaMott *= (Z*Z);
938 
939 // EventDataFile<<"SigmaMott="<<SigmaMott<<" "<<"SigmaMott*F_2="<<SigmaMott*F_2<<" "<<"weight_n="<<SigmaMott*F_2*sin(Theta)<<G4endl;
940 
941  G4double Xsect = SigmaMott*F_2;
942  Xsect = Xsect * (1.0 + delta_Schwinger);
943  fWeightN = Xsect*sin(Theta);
944 
945  // G4cout<< "q(1/fm)"<<sqrt(q2) <<" SigmaMott="<<SigmaMott<<" F_2="<<F_2<<" SigmaMott*F_2="<<Xsect<<G4endl;
946  return Xsect;
947 }
948 
949 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
950 //jpan@nuclear.uwinnipeg.ca Sat Apr 18 11:28:18 CDT 2009
951 // In this subroutine units are MeV
952 
953 G4double QweakSimEPEvent::Quasi_Elastic_Neutron(G4double E_in,
954  G4double Theta,
955  G4double &fWeightN,
956  G4double &Q2,
957  G4double &E_out)
958 {
959 // cout<<"===>>>>Calling Quasi_Elastic_Neutron"<<endl;
960  G4double Theta_Min = 1.745329E-4;
961  G4double Lamda_2 = 0.710;
962  G4double M_p = 938.2796; // proton mass in MeV
963 // G4double mu = 2.793;
964  G4double mu_n = -1.91;
965  G4double Z = 1.0;
966 // G4double A = 1.0;
967  G4double M = M_p;
968 
969 // E_in units is MeV
970 
971  if (Theta<Theta_Min)
972  Theta = Theta_Min;
973 
974  G4double CTH = cos(Theta/2.0);
975  G4double STH = sin(Theta/2.0);
976  G4double T2THE = STH*STH/CTH/CTH;
977  G4double ETA = 1.0+2.0*E_in*STH*STH/M;
978  E_out = E_in/ETA;
979 // MeV^2
980  Q2 = 4.*E_in*E_out*STH*STH;
981  G4double tau = Q2/4./M/M;
982 // Mott scatering
983  G4double CrossSection = (Z*0.72/E_in*CTH/STH/STH)*(Z*0.72/E_in*CTH/STH/STH)/ETA;
984 // Units: ub/sr
985  G4double Mott = CrossSection*10000.0;
986 
987 // Cross section
988  G4double GEP_DIPOLE = 1.0/(1.0+Q2/1.E6/Lamda_2)*(1.0+Q2/1.E6/Lamda_2);
989 // G4double GMP_DIPOLE = GEP_DIPOLE*mu;
990  G4double GEn_DIPOLE = 0.0;
991  G4double GMn_DIPOLE = GEP_DIPOLE*mu_n;
992  G4double FAC = 1.0/(1.0+tau);
993 
994 // G4double Sigma_Dipole_p = Mott*(GEP_DIPOLE*GEP_DIPOLE*FAC+tau*GMP_DIPOLE*GMP_DIPOLE*(FAC+2.*T2THE));
995  G4double Sigma_Dipole_n = Mott*(GEn_DIPOLE*GEn_DIPOLE*FAC+tau*GMn_DIPOLE*GMn_DIPOLE*(FAC+2.*T2THE));
996 // G4double Sigma_Dipole = Z*Sigma_Dipole_p + (A-Z)*Sigma_Dipole_n;
997 
998 // cout<<"===>>>>leaving Quasi_Elastic_Neutron"<<endl;
999  fWeightN = Sigma_Dipole_n*sin(Theta);
1000  return Sigma_Dipole_n;
1001 }
1002 
1003 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1004 
1005 // Inelastic generator interface
1006 // Peiqing, Nov 30,2011
1007 
1008 G4double QweakSimEPEvent::Delta_Resonance(G4double E_in,
1009  G4double Theta,
1010  G4double &fWeightN,
1011  G4double &Q2,
1012  G4double &E_out)
1013 {
1014  G4double E_beam = myUserInfo->GetBeamEnergy(); // maximum energy (beam energy) in MeV
1015  G4double M_electron = 0.511; // minimum energy (electron mass) in MeV
1016 
1017  // G4double Theta_Min = 1.745329E-4;
1018  if (Theta<Theta_Min)
1019  Theta = Theta_Min;
1020 
1021  // Generate flat energy distribution of outgoing electron
1022  E_out = M_electron + G4UniformRand()*(E_beam - M_electron);
1023 
1024  // TODO: total energy phase space should be reduced to improve the efficiency.
1025  G4double xsect = Sigma_EEPrime(E_in/1000.0, E_out/1000.0, Theta, Q2); // ub/sr/GeV
1026  Q2 = Q2*1e6; // convert to MeV^2 for Q2
1027 
1028  //std::cout<<"Q2="<<Q2<<" xsect="<<xsect<<std::endl;
1029  xsect = xsect*(E_beam - M_electron)/1000.0;
1030  //std::cout<<"xsect*1.165 = "<<xsect<<std::endl;
1031  fWeightN = xsect*sin(Theta);
1032 
1033  if(xsect == 0) // if E > E_max, reject the event
1034  {
1035  E_out = M_electron;
1036  Q2 = 0.0;
1037  }
1038 
1039  return xsect;
1040 }
1041 
1042 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1043 
1044 G4double QweakSimEPEvent::Sigma_EEPrime(G4double eni, G4double eprime, G4double theta, G4double &q2)
1045 {
1046 
1047 // Peiqing, Nov 28, 2011
1048 // This is adopted from QWGEANT3 SUBROUTINE SIGMA_EEPRIME (ENI,EPRIME,THETA,WVAL,SIGMA_EEP)
1049 //
1050 // Calculates the differential cross section for inelastic
1051 // ep scattering into unit solid angle and scattered
1052 // electron energy EPRIME. The cross section is
1053 // returned in microbarns/sr/GeV.
1054 
1055  G4double w2,xval1[51],xvalL[51];
1056  G4double sigT,sigL,epsilon;
1057  G4int i;
1058 
1059  G4double t2,epmax,xneu,gamma;
1060 
1061  G4double mp = 0.9382727;
1062  G4double mp2 = mp*mp;
1063  G4double alpha = 1.0/137.036;
1064  G4double mpion = 0.135;
1065 
1066  G4double xval[101]={ 0,
1067  0.12286E+01,0.15200E+01,0.15053E+01,0.16960E+01,0.16466E+01,
1068  0.14055E+01,0.12247E+00,0.22000E+00,0.84591E-01,0.90400E-01,
1069  0.16956E+00,0.29002E+00,0.74046E+01,0.65122E+01,0.58506E+01,
1070  0.16990E+01,0.45729E+01,0.53546E+03,-.97628E+03,0.82539E+05,
1071  0.29494E+01,0.20352E+01,0.12422E+00,0.19560E+01,0.22333E+01,
1072  -.32484E+00,0.24212E+00,0.57737E-01,0.30497E+01,0.63111E+01,
1073  0.37579E+00,0.41100E+01,0.20668E-01,0.45490E+03,0.54493E-01,
1074  0.22565E+01,0.40369E+01,0.43734E+00,0.65625E+00,0.36182E+00,
1075  0.55216E-01,0.41789E+00,0.18104E+00,0.91306E+00,0.18116E+01,
1076  0.24152E+03,0.19329E+01,0.38000E+00,0.28732E+01,-.53116E+00,
1077  0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,
1078  0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,
1079  0.00000E+00,0.00000E+00,0.99495E+01,0.70370E-02,0.16172E+01,
1080  0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,
1081  0.20782E+05,0.79523E+04,0.54933E+00,0.00000E+00,0.68629E+01,
1082  0.10369E+01,0.88112E+00,0.00000E+00,0.35659E+03,0.20281E-02,
1083  0.40336E+02,0.00000E+00,0.45242E+00,0.20000E-01,0.28691E+00,
1084  0.00000E+00,0.25115E+03,0.10663E+01,0.55422E+01,0.30350E+00,
1085  0.10541E+01,0.20527E+01,0.00000E+00,0.00000E+00,0.13055E+02,
1086  0.27997E+01,0.39546E+00,0.46372E+00,0.49972E+01,0.00000E+00};
1087 
1088 // Check that we are not above maximum allowed E assuming threshold W=Mp+Mpi
1089  t2 = theta/2.0;
1090  G4double sin_t2 = sin(t2);
1091  G4double tan_t2 = tan(t2);
1092  G4double sin2_t2 = sin_t2*sin_t2;
1093  G4double tan2_t2 = tan_t2*tan_t2;
1094  G4double sigma_eep = 0.0;
1095  epmax = eni-mpion*(1.0+0.5*mpion/mp);
1096  epmax = epmax/(1.0+(2.0*eni/mp)*sin2_t2);
1097  if (eprime>=epmax)
1098  return 0;
1099 
1100 // Calculate Q**2, W**2, epsilon, gamma
1101  q2 = 4.0*eni*eprime*sin2_t2;
1102  xneu = eni-eprime;
1103  w2 = mp2 + 2.0*mp*xneu - q2;
1104  epsilon = 1.0/(1.0+2.0*(1.0+xneu*xneu/q2)*tan2_t2);
1105  gamma = alpha*eprime*(w2-mp2)/((2.0*pi)*(2.0*pi));
1106  gamma = gamma/(q2*mp*eni*(1.0-epsilon));
1107  // wval = sqrt(w2);
1108 
1109  for (i=1;i<=50;i++)
1110  {
1111  xval1[i] = xval[i];
1112  xvalL[i] = xval[50+i];
1113  if(i<=12)
1114  xvalL[i] = xval1[i];
1115  }
1116 
1117  xvalL[43] = xval1[47];
1118  xvalL[44] = xval1[48];
1119  xvalL[50] = xval1[50];
1120 
1121  sigT = ResMod507(1,w2,q2,&xval1[0]);
1122  sigL = ResMod507(2,w2,q2,&xvalL[0]);
1123  sigma_eep=sigT+epsilon*sigL;
1124  sigma_eep=sigma_eep*gamma;
1125 
1126  return sigma_eep;
1127 }
1128 
1129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1130 
1131 
1132 G4double QweakSimEPEvent::ResMod507(G4int sf, G4double w2, G4double q2, G4double* xval)
1133 {
1134  G4double xb,mass[8],width[8];
1135  G4double height[8],rescoef[7][5];
1136  G4double nr_coef[4][5],sigr[8],sig_nr;
1137  G4double intwidth[8],k,kcm,kr[8],kcmr[8],ppicm,ppi2cm;
1138  G4double petacm,ppicmr[8],ppi2cmr[8],petacmr[8],epicmr[8],epi2cmr[8];
1139  G4double eetacmr[8],epicm,epi2cm,eetacm,br[8][4],ang[8];
1140  G4double pgam[8],pwid[8][4],x0[8],q20;
1141  G4double sig_res,t,xpr[3],m0;
1142  G4int i,j,num;
1143 
1144  G4double xt,p1,p2,p3;
1145 
1146  G4double mp = 0.9382727;
1147  G4double mpi = 0.135;
1148  G4double meta = 0.547;
1149  G4double mp2 = mp*mp;
1150  G4double w = sqrt(w2);
1151  //wdif[1] = w - (mp + mpi);
1152  //wdif[2] = w - (mp + 2.*mpi);
1153 
1154  m0 = 0.125;
1155  if(sf==2)
1156  m0 = xval[49];
1157 
1158  if(sf==1)
1159  {
1160  q20 = 0.05;
1161  }
1162  else
1163  {
1164  q20 = 0.300;
1165  }
1166 
1167 // single pion branching ratios
1168  br[1][1] = 1.0; // P33(1232)
1169  br[2][1] = 0.45; // S11(1535)
1170  br[3][1] = 0.65; // D13(1520)
1171  br[4][1] = 0.65; // F15(1680)
1172  br[5][1] = 0.4; // S11(1650)
1173  br[6][1] = 0.65; // P11(1440) roper
1174  br[7][1] = 0.50; // F37(1950)
1175 
1176 // eta branching ratios
1177  br[1][3] = 0.0; // !!! P33(1232)
1178  br[2][3] = 0.45; // !!! S11(1535)
1179  br[3][3] = 0.0; // !!! D13(1520)
1180  br[4][3] = 0.0; // !!! F15(1680)
1181  br[5][3] = 0.1; // !!! S11(1650)
1182  br[6][3] = 0.0; // !!! P11(1440) roper
1183  br[7][3] = 0.0; // !!! F37(1950)
1184 
1185 // 2-pion branching ratios
1186  for (i=1;i<=7;i++)
1187  {
1188  br[i][2] = 1.-br[i][1]-br[i][3];
1189  }
1190 
1191 // Meson angular momentum
1192  ang[1] = 1.; // P33(1232)
1193  ang[2] = 0.; // S11(1535)
1194  ang[3] = 2.; // D13(1520)
1195  ang[4] = 3.; // F15(1680)
1196  ang[5] = 0.; // S15(1650)
1197  ang[6] = 1.; // P11(1440) roper
1198  ang[7] = 3.; // F37(1950)
1199 
1200  for( i=1;i<=7;i++) // resonance damping parameter
1201  {
1202  x0[i] = 0.145;
1203  }
1204 
1205  x0[1] = 0.130;
1206 
1207  for (i=1;i<=7;i++)
1208  {
1209  br[i][2] = 1.-br[i][1]-br[i][3];
1210  }
1211 
1212  // dip = 1./((1.+q2/0.71)*(1.+q2/0.71)); // Dipole parameterization
1213  // mon = 1./(1.+q2/0.71);
1214 
1215  xb = q2/(q2+w2-mp2);
1216  xpr[1] = 1.+(w2-(mp+mpi)*(mp+mpi))/(q2+q20);
1217  xpr[1] = 1./xpr[1];
1218  xpr[2] = 1.+(w2-(mp+mpi+mpi)*(mp+mpi+mpi))/(q2+q20);
1219  xpr[2] = 1./xpr[2];
1220 
1221  t = log(log((q2+m0)/0.330*0.330)/log(m0/0.330*0.330));
1222 
1223 // Calculate kinematics needed for threshold Relativistic B-W
1224 
1225  k = (w2 - mp2)/2./mp;
1226  kcm = (w2-mp2)/2./w;
1227 
1228  epicm = (w2 + mpi*mpi -mp2 )/2./w;
1229  ppicm = sqrt(fmax(0.0,(epicm*epicm - mpi*mpi)));
1230  epi2cm = (w2 + (2.*mpi)*(2.*mpi) -mp2 )/2./w;
1231  ppi2cm = sqrt(fmax(0.0,(epi2cm*epi2cm - (2.*mpi)*(2.*mpi))));
1232  eetacm = (w2 + meta*meta -mp2 )/2./w;
1233  petacm = sqrt(fmax(0.0,(eetacm*eetacm - meta*meta)));
1234 
1235  num = 0;
1236 
1237  for (i=1;i<=6;i++) // Read in resonance masses
1238  {
1239  num = num + 1;
1240  mass[i] = xval[i];
1241  }
1242 
1243  for (i=1;i<=6;i++) // Read in resonance widths
1244  {
1245  num = num + 1;
1246  intwidth[i] = xval[num];
1247  width[i] = intwidth[i];
1248  }
1249 
1250  if(sf==2) // Put in 4th resonance region
1251  {
1252  mass[7] = xval[43];
1253  intwidth[7] = xval[44];
1254  width[7] = intwidth[7];
1255  }
1256  else
1257  {
1258  mass[7] = xval[47];
1259  intwidth[7] = xval[48];
1260  width[7] = intwidth[7];
1261  }
1262 
1263  for (i=1;i<=7;i++)
1264  {
1265  kr[i] = (mass[i]*mass[i]-mp2)/2./mp;
1266  kcmr[i] = (mass[i]*mass[i]-mp2)/2./mass[i];
1267  epicmr[i] = (mass[i]*mass[i] + mpi*mpi -mp2 )/2./mass[i];
1268  ppicmr[i] = sqrt(fmax(0.0,(epicmr[i]*epicmr[i] - mpi*mpi)));
1269  epi2cmr[i] = (mass[i]*mass[i] + (2.*mpi)*(2.*mpi) -mp2 )/2./mass[i];
1270  ppi2cmr[i] = sqrt(fmax(0.0,(epi2cmr[i]*epi2cmr[i] - (2.*mpi)*(2.*mpi))));
1271  eetacmr[i] = (mass[i]*mass[i] + meta*meta -mp2 )/2./mass[i];
1272  petacmr[i] = sqrt(fmax(0.0,(eetacmr[i]*eetacmr[i] - meta*meta)));
1273 
1274  // Calculate partial widths
1275  pwid[i][1] = intwidth[i]*pow(ppicm/ppicmr[i],(2.*ang[i]+1.))
1276  *pow((ppicmr[i]*ppicmr[i]+x0[i]*x0[i])/(ppicm*ppicm+x0[i]*x0[i]),ang[i]);
1277  // 1-pion decay mode
1278 
1279  pwid[i][2] = intwidth[i]*pow(ppi2cm/ppi2cmr[i],(2.*ang[i]+4.))
1280  *pow((ppi2cmr[i]*ppi2cmr[i]+x0[i]*x0[i])/(ppi2cm*ppi2cm+x0[i]*x0[i]),(ang[i]+2));
1281  // 2-pion decay mode
1282 
1283  pwid[i][2] = w/mass[i]*pwid[i][2];
1284  pwid[i][3] = 0.; // !!! eta decay mode
1285 
1286  if(i==2 || i==5)
1287  {
1288  pwid[i][3]=intwidth[i]*pow(petacm/petacmr[i],(2.*ang[i]+1.))
1289  *pow((petacmr[i]*petacmr[i]+x0[i]*x0[i])/(petacm*petacm+x0[i]*x0[i]),ang[i]);
1290  // eta decay only for S11's
1291  }
1292 
1293  pgam[i] = (kcm/kcmr[i])*(kcm/kcmr[i])*(kcmr[i]*kcmr[i]+x0[i]*x0[i])/(kcm*kcm+x0[i]*x0[i]);
1294  pgam[i] = intwidth[i]*pgam[i];
1295  width[i] = br[i][1]*pwid[i][1]+br[i][2]*pwid[i][2]+br[i][3]*pwid[i][3];
1296  }
1297  // End resonance kinematics and Widths calculations
1298 
1299  // Begin resonance Q^2 dependence calculations
1300  for (i=1;i<=6;i++)
1301  {
1302  for (j=1;j<=4;j++)
1303  {
1304  num = num + 1;
1305  rescoef[i][j]=xval[num];
1306  }
1307 
1308  if(sf==1)
1309  {
1310  height[i] = rescoef[i][1]*(1.+rescoef[i][2]*q2/(1.+rescoef[i][3]*q2))/ pow(1.+q2/0.91,rescoef[i][4]);
1311  }
1312  else
1313  {
1314  height[i] = rescoef[i][1]*q2/(1.+rescoef[i][2]*q2)*exp(-1.*rescoef[i][3]*q2);
1315  }
1316 
1317  height[i] = height[i]*height[i];
1318  }
1319 
1320 
1321  if(sf==2) // 4th resonance region
1322  {
1323  height[7] = xval[45]*q2/(1.+xval[46]*q2)*exp(-1.*xval[47]*q2);
1324  }
1325  else
1326  {
1327  height[7] = xval[49]/(1.+q2/0.91);
1328  }
1329 
1330  height[7] = height[7]*height[7];
1331 
1332 // End resonance Q^2 dependence calculations
1333 
1334  for (i=1;i<=3;i++) // Non-Res coefficients
1335  {
1336  for (j=1;j<=4;j++)
1337  {
1338  num = num + 1;
1339  nr_coef[i][j]=xval[num];
1340  }
1341  }
1342 
1343 // Calculate Breit-Wigners for all resonances
1344 
1345  sig_res = 0.0;
1346  for (i=1;i<=7;i++)
1347  {
1348  sigr[i] = width[i]*pgam[i]/((w2 - mass[i]*mass[i])*(w2 - mass[i]*mass[i])
1349  + (mass[i]*width[i])*(mass[i]*width[i]));
1350  sigr[i] = height[i]*kr[i]/k*kcmr[i]/kcm*sigr[i]/intwidth[i];
1351  sig_res = sig_res + sigr[i];
1352  }
1353  sig_res = sig_res*w;
1354 
1355 // Finish resonances / start non-res background calculation
1356 
1357  sig_nr = 0.;
1358  if(sf==1)
1359  {
1360  xt = (q2+xval[37])/(w2-1.15+q2+xval[37]);
1361  p2 = xval[38]*exp(-1.*xval[39]*q2)+xval[40]*log(q2/xval[41]+1.);
1362  p3 = xval[42]*exp(-1.*xval[43]*q2)+xval[44]*log(q2/xval[45]+1.);
1363  p1 = xval[46]*(1.+xval[50]*q2/(1.+0.534*q2));
1364 
1365  sig_nr = p1*pow(xt,p2)*pow(1.-xt,p3);
1366  }
1367  else if(sf==2)
1368  {
1369 
1370  for (i=1;i<=1;i++)
1371  {
1372  sig_nr = sig_nr + nr_coef[i][1]*pow(1.-xpr[i],(nr_coef[i][3]+nr_coef[i][2]*t)) /(1.-xb)/(q2+q20)
1373  *pow(q2/(q2+q20), nr_coef[i][4]) *pow(xpr[i],(xval[41]+xval[42]*t));
1374  }
1375  }
1376 
1377  G4double sig = sig_res + sig_nr;
1378  return sig;
1379 }
1380 
1381 ////....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1382 //
1383 // Peiqing, Dec. 11, 2011, Moller elastic scattering
1384 //
1385 // kinematic relations and cross section formulae
1386 // are from Wagner, et al. NIMA294 (1990) 541-548
1387 // and Arrington, et al. NIMA311 (1992) 39-48, etc.
1388 
1389 G4double QweakSimEPEvent::Moller_Scattering(G4double E_in, G4double theta1,
1390  G4double &E_out1, G4double &E_out2, G4double &theta2,
1391  G4double &q2, G4double &fWeightN, G4double &asymmetry)
1392 {
1393  G4double alpha = 1.0/137.036;
1394  G4double G_F = 1.166e-11; // Fermi constant in MeV^(-2)
1395  G4double M_electron = 0.511; // Incident Electron mass in MeV
1396  G4double sin2_theta_W = 0.02381; // sin2_theta_W at 0.026 GeV, taken from E158
1397 
1398  G4double theta_CM = 2.0*atan(tan(theta1)*sqrt((E_in+M_electron)/(2.*M_electron)));
1399  E_out1 = M_electron+(E_in-M_electron)*(cos(theta_CM/2.0)*cos(theta_CM/2.0));
1400  E_out2 = E_in - E_out1 + M_electron;
1401  theta2 = atan(2.0*M_electron/(E_in+M_electron)/tan(theta1));
1402  //G4double momentum = sqrt(E_out1*E_out1 - M_electron*M_electron);
1403  // compute q2 in MeV^2
1404  //q2 = 4.0*E_in*momentum*sin(theta1/2.0)*sin(theta1/2.0);
1405  q2 = 2*M_electron*(M_electron+E_in)*sin(theta_CM/2.0)*sin(theta_CM/2.0);
1406 
1407  // calculate Moller cross-section (units are ub/sr)
1408  G4double Xsect = pow((1+cos(theta_CM))*(3+cos(theta_CM)*cos(theta_CM)),2);
1409  // Xsect = pow(alpha/2.0/M_electron,2)*Xsect/pow(sin(theta_CM),4); // alpha^2/(4m^2) = 0.199 b/Sr
1410  Xsect = (1.99e4)*Xsect/pow(sin(theta_CM),4);
1411  fWeightN = Xsect*sin(theta1);
1412  asymmetry = -M_electron*E_in*(G_F/(sqrt(2.0)*alpha*3.1415927))
1413  *16*sin(theta_CM)*sin(theta_CM)/(3+cos(theta_CM)*cos(theta_CM))*(1-4*sin2_theta_W);
1414  asymmetry = asymmetry*1e6; // in [ppm]
1415 
1416  //std::cout<<"E_in E1 E2 theta1 theta2 Xsect Weight Q2 asym"<<std::endl;
1417  //std::cout<<E_in<<" "<<E_out1<<" "<<E_out2<<" "<<theta1*180/3.14<<" "
1418  // <<theta2*180/3.14<<" "<<Xsect<<" "<<fWeightN<<" "<<q2*1e-6<<" "<<asymmetry<<std::endl;
1419 
1420  return Xsect;
1421 }
1422 
1423 ////....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1424 // Jim Dowd
1425 // ---------------------------------------------------------
1426 // Calculates the Cross Section weighting factor for
1427 // radiated scattering from a lookup table.
1428 // ---------------------------------------------------------
1429 //
1430 // Beam Energy must be 3.35 GeV, 1.16 GeV, or 877 MeV
1431 
1432 const std::vector< G4double > QweakSimEPEvent::Radiative_Cross_Section_Lookup(G4double E_in,
1433  G4double Theta,
1434  G4double &fWeightN,
1435  G4double &Q2,
1436  G4double &E_out)
1437 {
1438  //G4int SuperElastic = 1;
1439  //G4double xsec = 0.0;
1440  Double_t value[value_n] = {0.0};
1441  std::vector< G4double > CrossSection;
1442  G4double xsec = 0.0;
1443  if (Theta<GetThetaAngle_Min())
1444  Theta = GetThetaAngle_Min();
1445 
1446  //while (SuperElastic) {
1447  //E_out = (G4UniformRand()*(GetEPrime_Max()-GetEPrime_Min()) + GetEPrime_Min());
1448  //SuperElastic = SuperElasticCheck(E_in, E_out, Theta, xsec);
1449  //G4cout << "El XSec: " << xsec << G4endl;
1450  //}
1451 
1452  E_out = myUserInfo->GetEPrime();
1453 
1454  Double_t Zpos;
1455 
1456  if ( ReactionRegion == 1 ) {
1457  Zpos = myUserInfo->GetOriginVertexPositionZ();
1458  }
1459  else {
1460  Zpos = 0.0;
1461  }
1462  const Double_t pos[value_d] = {Zpos,E_in,E_out,Theta};
1463  //Double_t value[value_n] = {0.0};
1464  //std::vector< G4double > CrossSection;
1465 
1466  fLookupTable->GetValue(pos,value);
1467  Q2 = value[1];
1468  fWeightN = value[6]*sin(Theta*degree);
1469  //fWeightN = value[6];
1470  //G4cout << "Rad XSec: " << value[12]*0.001 << G4endl;
1471  //}
1472 
1473  CrossSection.push_back(value[6]*0.001); // 0 total radiated cross section
1474  CrossSection.push_back(value[2]*0.001); // 1 total born cross section
1475  CrossSection.push_back(value[3]*0.001); // 2 inelastic born cross section
1476  CrossSection.push_back(value[4]*0.001); // 3 quasi-elastic born cross section
1477  CrossSection.push_back(value[6]*0.001); // 4 total radiated cross section
1478  CrossSection.push_back(value[7]*0.001); // 5 elastic radiated cross section
1479  CrossSection.push_back(value[8]*0.001); // 6 quasi-elastic radiated cross section
1480  CrossSection.push_back(value[9]*0.001); // 7 deep-inelastic radiated cross section
1481  CrossSection.push_back(value[11]*0.001); // 8 total radiated cross section internal only
1482  CrossSection.push_back(value[12]*0.001); // 9 elastic radiated cross section internal only
1483  CrossSection.push_back(value[14]*0.001); // 10 quasi-elastic radiated cross section internal only
1484  CrossSection.push_back(value[13]*0.001); // 11 deep-inelastic radiated cross section internal only
1485  CrossSection.push_back(0.0); // 12 elastic peak cross section
1486 
1487  if (SuperElasticCheck(myUserInfo->GetBeamEnergy(), E_out, Theta, xsec)) CrossSection[5] = 0.0;
1488  else if (SuperElasticCheck(myUserInfo->GetBeamEnergy(), E_out+ElasticPeakDeltaE, Theta, xsec)) {
1489  CrossSection[12] = xsec/(ElasticPeakDeltaE*0.001);
1490  CrossSection[5] = 0.0;
1491  }
1492  CrossSection[4] = CrossSection[5] + CrossSection[6] + CrossSection[7] + CrossSection[12];
1493  CrossSection[0] = CrossSection[4];
1494 
1495  return CrossSection; // units [ub/sr/GeV]
1496 }
1497 
1498 ////....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1499 // Jim Dowd
1500 // ---------------------------------------------------------
1501 // Checks the event to see if it is Super-Elastic or not.
1502 // ---------------------------------------------------------
1503 
1504 
1505 G4int QweakSimEPEvent::SuperElasticCheck(G4double E_in, G4double E_out, G4double theta, G4double &xsec)
1506 {
1507  if (myUserInfo->GetReactionType()!=7) return 0; // Aborts if not using lookup table
1508 
1509  G4double Weight; // dummy
1510  G4double Q2; // dummy
1511 
1512  G4double E_elastic;
1513 
1514  xsec = Elastic_Cross_Section_Proton(E_in, theta, Weight, Q2, E_elastic);
1515 
1516  if (E_out > E_elastic) {
1517  //G4cout << "###################################################" << G4endl;
1518  //G4cout << "---------------------------------------------------" << G4endl;
1519  //G4cout << "E_in: " << E_in << G4endl;
1520  //G4cout << "E_out: " << E_out << G4endl;
1521  //G4cout << "E_elastic:" << E_elastic << G4endl;
1522  //G4cout << "Theta: " << theta*180.0/3.14159 << G4endl;
1523  //G4cout << "Event is SuperElestic! Generating new event." << G4endl;
1524  return 1;
1525  }
1526  else {
1527  //G4cout << "###################################################" << G4endl;
1528  //G4cout << "E_in: " << E_in << G4endl;
1529  //G4cout << "E_out: " << E_out << G4endl;
1530  //G4cout << "E_elastic:" << E_elastic << G4endl;
1531  //G4cout << "Theta: " << theta*180.0/3.14159 << G4endl;
1532  //G4cout << "Event is NOT SuperElastic" << G4endl;
1533  return 0;
1534  }
1535 }
1536 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1537 
1538 
1539 
1540 
1541 ////....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1542 // Jim Dowd
1543 // ---------------------------------------------------------
1544 // Creates lookup table for calculating elastic radiative
1545 // cross sections on hydrogen.
1546 // ---------------------------------------------------------
1547 
1548 void QweakSimEPEvent::CreateLookupTable()
1549 {
1550 
1551  G4cout << "---> Calling CreateLookupTable() <---" << G4endl;
1552 
1553  G4double energy = GetBeamEnergy();
1554  if (energy != 3350*MeV && energy != 1160*MeV && energy != 877*MeV) {
1555  G4cout << "#### Current beam energy is not a valid choice for lookup table" << G4endl;
1556  G4cout << "#### Setting beam energy to default value for lookup table (3.35 GeV)" << G4endl;
1557  energy = 3350*MeV;
1558  SetBeamEnergy(energy);
1559  }
1560 
1561  TString target;
1562  if (ReactionRegion == 1) target = "lh2";
1563  else if (ReactionRegion == 2) target = "USALWindow";
1564  else if (ReactionRegion == 3) target = "DSALWindow";
1565  else if (ReactionRegion == 4) target = "USALDummy1";
1566  else if (ReactionRegion == 5) target = "USALDummy2";
1567  else if (ReactionRegion == 6) target = "USALDummy4";
1568  else if (ReactionRegion == 7) target = "DSALDummy2";
1569  else if (ReactionRegion == 8) target = "DSALDummy4";
1570  else if (ReactionRegion == 9) target = "DSALDummy8";
1571  else if (ReactionRegion == 10) target = "USCDummy";
1572  else if (ReactionRegion == 11) target = "DSCDummy";
1573 
1574  TString filename = "radiative_cross_section_";
1575  filename += target;
1576  filename += "_";
1577  filename += (Int_t)energy;
1578  filename += "MeV";
1579  //if (ReactionRegion == 1 && energy != 3350*MeV) filename += "_v2";
1580  if (energy != 3350*MeV) filename += "_v2";
1581  filename += ".dat";
1582 
1583  TString filepath = "./radiative_lookup_tables/";
1584 
1585  std::ifstream in;
1586  in.open(filepath + filename);
1587  if (!in.is_open()) {
1588  G4cout << "#### Failed to open lookup table data file \"" << filepath << filename << "\"" << G4endl;
1589  return;
1590  }
1591  else
1592  {
1593  in.close();
1594  G4cout << "#### Found lookup table data file \"" << filepath << filename << "\"" << G4endl;
1595  if (fLookupTable != 0) delete fLookupTable;
1596  fLookupTable = new QweakSimFieldMap<value_t,value_n>(value_d);
1597  fLookupTable->ReadTextFile((filepath +filename).Data());
1598  }
1599  G4cout << "+++++++++++++++++++++++++++++" << G4endl;
1600  G4cout << "++ Lookup Table Parameters ++" << G4endl;
1601  G4cout << "+++++++++++++++++++++++++++++" << G4endl;
1602  G4cout << "Beam Energy: " << fLookupTable->GetMinimum(1)/1000 << " GeV" << G4endl;
1603  G4cout << "Target: " << target << G4endl;
1604  G4cout << "E\' Lower Bound: " << fLookupTable->GetMinimum(2)/1000 << " GeV" << G4endl;
1605  G4cout << "E\' Upper Bound: " << fLookupTable->GetMaximum(2)/1000 << " GeV" << G4endl;
1606  G4cout << "Theta Lower Bound: " << fLookupTable->GetMinimum(3)/degree << " degrees" << G4endl;
1607  G4cout << "Theta Upper Bound: " << fLookupTable->GetMaximum(3)/degree << " degrees" << G4endl;
1608  G4cout << "+++++++++++++++++++++++++++++" << G4endl;
1609 
1610  if (energy != GetBeamEnergy()) SetBeamEnergy(energy);
1611  else CheckLookupTableBounds();
1612 
1613 
1614  G4cout << "---> Leaving CreateLookupTable() <---" << G4endl;
1615 
1616 
1617  //********************************************************************
1618 
1619  /*
1620  // **********************************************************************
1621  // **********************************************************************
1622  // The code below is used to convert the raw output of the Bosted code
1623  // into a format the GEANT4 can read.
1624  // **********************************************************************
1625 
1626 
1627  //G4int entries = 93610;
1628  G4int entries = 1;
1629  //G4int entries[3] = { 1295, 1665, 1073};
1630 
1631  coord_t coord[value_d] = {0.0};
1632  value_t field[value_n] = {0.0};
1633 
1634  std::vector< G4double > fMin;
1635  std::vector< G4double > fMax;
1636  std::vector< G4double > fStep;
1637 
1638  //std::vector< G4double > fMin[3];
1639  //std::vector< G4double > fMax[3];
1640  //std::vector< G4double > fStep[3];
1641 
1642  std::string target[10] = { "DSALDummy2",
1643  "DSALDummy4",
1644  "DSALDummy8",
1645  "USALDummy1",
1646  "USALDummy2",
1647  "USALDummy4",
1648  "DSALWindow",
1649  "USALWindow",
1650  "DSCDummy",
1651  "USCDummy" };
1652 
1653  std::string energy_in[3] = { "877mev",
1654  "1gev",
1655  "3gev" };
1656 
1657  std::string energy_out[3] = { "877MeV",
1658  "1160MeV",
1659  "3350MeV" };
1660 
1661  std::string filepath = "/group/qweak/jdowd/externals_qweak/OUT/";
1662  std::string name, name1, name2, name3;
1663 
1664  // Set min, max, and step size for each coordinate in the lookup table.
1665  // Z position in radiation lengths
1666  fMin.push_back(myUserInfo->TargetCenterPositionZ - 0.5*myUserInfo->TargetLength);
1667  fMax.push_back(myUserInfo->TargetCenterPositionZ + 0.5*myUserInfo->TargetLength);
1668  fStep.push_back(myUserInfo->TargetLength/10.0);
1669 
1670  //fMin.push_back(0.0);
1671  //fMax.push_back(0.0);
1672  //fStep.push_back(1.0);
1673 
1674  //fMin[0].push_back(0.0);
1675  //fMax[0].push_back(0.0);
1676  //fStep[0].push_back(1.0);
1677 
1678  //fMin[1].push_back(0.0);
1679  //fMax[1].push_back(0.0);
1680  //fStep[1].push_back(1.0);
1681 
1682  //fMin[2].push_back(0.0);
1683  //fMax[2].push_back(0.0);
1684  //fStep[2].push_back(1.0);
1685 
1686 
1687  // Beam Energy in GeV
1688  fMin.push_back(1.16*GeV);
1689  fMax.push_back(1.16*GeV);
1690  fStep.push_back(0.005*GeV);
1691 
1692  //fMin[0].push_back(0.877*GeV);
1693  //fMax[0].push_back(0.877*GeV);
1694  //fStep[0].push_back(0.05*GeV);
1695 
1696  //fMin[1].push_back(1.16*GeV);
1697  //fMax[1].push_back(1.16*GeV);
1698  //fStep[1].push_back(0.05*GeV);
1699 
1700  //fMin[2].push_back(3.35*GeV);
1701  //fMax[2].push_back(3.35*GeV);
1702  //fStep[2].push_back(0.05*GeV);
1703 
1704  // E prime in GeV
1705  fMin.push_back(0.015*GeV);
1706  fMax.push_back(1.16*GeV);
1707  fStep.push_back(0.005*GeV);
1708 
1709  //fMin[0].push_back(0.026*GeV);
1710  //fMax[0].push_back(0.876*GeV);
1711  //fStep[0].push_back(0.025*GeV);
1712 
1713  //fMin[1].push_back(0.059*GeV);
1714  //fMax[1].push_back(1.159*GeV);
1715  //fStep[1].push_back(0.025*GeV);
1716 
1717  //fMin[2].push_back(0.150*GeV);
1718  //fMax[2].push_back(1.550*GeV);
1719  //fStep[2].push_back(0.050*GeV);
1720 
1721  // Theta angle
1722  fMin.push_back(2.00*degree);
1723  fMax.push_back(20.0*degree);
1724  fStep.push_back(0.5*degree);
1725 
1726  //fMin[0].push_back(2.00*degree);
1727  //fMax[0].push_back(20.0*degree);
1728  //fStep[0].push_back(0.5*degree);
1729 
1730  //fMin[1].push_back(2.00*degree);
1731  //fMax[1].push_back(20.0*degree);
1732  //fStep[1].push_back(0.5*degree);
1733 
1734  //fMin[2].push_back(2.00*degree);
1735  //fMax[2].push_back(20.0*degree);
1736  //fStep[2].push_back(0.5*degree);
1737 
1738  //G4cout << "Target Min: " << fMin[0] << G4endl;
1739  //G4cout << "Target Max: " << fMax[0] << G4endl;
1740 
1741  for (Int_t n = 0; n < (Int_t)fStep.size(); n++) {
1742  entries *= (G4int)( (fMax[n]-fMin[n])/fStep[n] + 1.5 );
1743  }
1744  std::ifstream in;
1745 
1746  G4cout << "!!!!!!!!!!!!!!!!!!!!!! Debug !!!!!!!!!!!!!!!!!!!!!!" << G4endl;
1747  G4cout << " fMin: " << fMin[0] << " " << fMin[1] << " " << fMin[2] << " " << fMin[3] << G4endl;
1748  G4cout << " fMax: " << fMax[0] << " " << fMax[1] << " " << fMax[2] << " " << fMax[3] << G4endl;
1749  G4cout << " fStep: " << fStep[0]<< " " << fStep[1]<< " " << fStep[2]<< " " << fStep[3]<< G4endl;
1750  G4cout << "===================================================================" << G4endl;
1751 
1752  //for (Int_t target_num = 0; target_num < 10; target_num++) {
1753  for (Int_t target_num = 0; target_num < 1; target_num++) {
1754  //for (Int_t energy_num = 0; energy_num < 3; energy_num++) {
1755  for (Int_t energy_num = 0; energy_num < 1; energy_num++) {
1756  //in.open("./radiative_lookup.dat");
1757 
1758  name = filepath + "lookup_" + target[target_num] + "_" + energy_out[energy_num] + ".out";
1759 
1760  in.open("./test_v2.out");
1761  //in.open(name.c_str());
1762  if (!in.is_open())
1763  G4cout << "#### Failed to open data file for lookup table" << G4endl;
1764  else
1765  {
1766  G4cout << "#### Found lookup table data file" << G4endl;
1767  if (fLookupTable != 0) delete fLookupTable;
1768  fLookupTable = new QweakSimFieldMap<value_t,value_n>(value_d);
1769  //fLookupTable->SetMinimumMaximumStep(fMin[energy_num],fMax[energy_num],fStep[energy_num]);
1770  fLookupTable->SetMinimumMaximumStep(fMin,fMax,fStep);
1771 
1772 
1773  //for (Int_t n = 0; n < (Int_t)fStep[energy_num].size(); n++) {
1774  // entries *= (G4int)( (fMax[energy_num][n]-fMin[energy_num][n])/fStep[energy_num][n] + 1.5 );
1775  //}
1776 
1777  // Filling Lookup Table
1778  //for (Int_t line = 0; line < entries[energy_num]; line++) {
1779  for (Int_t line = 0; line < entries; line++) {
1780  if (in.peek() == EOF) {
1781  G4cout << "#### Error reading \'elastic_lookup.dat\': File contains only "
1782  << line << " of " << entries << " expected lines. ####" << G4endl;
1783  break;
1784  }
1785 
1786  for (G4int i = 0; i < value_d; i++) in >> coord[i];
1787  for (G4int j = 0; j < value_n; j++) in >> field[j];
1788 
1789  //G4cout << "Initial: " << coord[0] << " " << coord[1] << " " << coord[2] << " " << coord[3] << G4endl;
1790 
1791  // Add units
1792  // Z position converted from rad lengths to GEANT coords
1793  coord[0] = myUserInfo->TargetCenterPositionZ - 0.5*myUserInfo->TargetLength
1794  + coord[0]*myUserInfo->TargetLength/0.0396;
1795  coord[1] *= GeV; // Beam Energy
1796  coord[2] *= GeV; // E prime
1797  coord[3] *= degree; // Theta
1798 
1799  //G4cout << "Final: " << coord[0] << " " << coord[1] << " " << coord[2] << " " << coord[3] << G4endl;
1800  //G4cout << "-----------------------------------------------------------------" << G4endl;
1801 
1802  //field[0] // Bjorken x (unitless) // This value has too few sig figs
1803  field[1] *= GeV*GeV; // Q2 // This value has too few sig figs
1804  //field[2] // total born cross section (ub/Sr)
1805  //field[3] // inelastic born cross section(ub/Sr)
1806  //field[4] // quasi-elastic born cross section (ub/Sr)
1807  //field[5] // quasi-elastic factor
1808  //field[6] // total radiated cross section (ub/Sr)
1809  //field[7] // elastic radiated cross section (ub/Sr)
1810  //field[8] // quasi-elastic radiated cross section (ub/Sr)
1811  //field[9] // deep-inelastic radiated cross section (ub/Sr)
1812  //field[10] // charge correction
1813  //field[11] // total radiated cross section internal only (ub/Sr)
1814  //field[12] // elastic radiated cross section internal only (ub/Sr)
1815  //field[13] // deep-inelastic radiated cross section internal only (ub/Sr)
1816  //field[14] // quasi-elastic radiated cross section internal only (ub/Sr)
1817 
1818  fLookupTable->Set(coord,field);
1819  }
1820  G4cout << "===== Filling of Lookup Table complete! =====" << G4endl;
1821  }
1822  in.close();
1823  fLookupTable->WriteTextFile("./radiative_lookup_1160_v2.out");
1824 
1825  name = filepath + "final/radiative_lookup_" + target[target_num] + "_" + energy_out[energy_num] + ".dat";
1826 
1827  //fLookupTable->WriteTextFile("./radiative_lookup_test.out");
1828  //fLookupTable->WriteTextFile(name.c_str());
1829 
1830  G4cout << "Wrote: " << name << G4endl;
1831 
1832  }
1833  }
1834  */
1835 }
1836 
1837 
1838 ////....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1839 // Jim Dowd
1840 // ---------------------------------------------------------
1841 // Checks the ranges of input parameters (E, E', and theta)
1842 // and insures they are within the lookup table bounds.
1843 // ---------------------------------------------------------
1844 
1845 void QweakSimEPEvent::CheckLookupTableBounds()
1846 {
1847  G4cout << "---> Calling CheckLookupTableBounds() <---" << G4endl;
1848  if (fLookupTable == 0) return;
1849  if (GetBeamEnergy() != (G4double)fLookupTable->GetMinimum(1) ) {
1850  G4cout << "!!!! Beam Energy is out of bounds" << G4endl;
1851  G4cout << "---> Changing Beam Energy to default value (" << fLookupTable->GetMinimum(1)/1000 << " GeV)" << G4endl;
1852  SetBeamEnergy( (G4double)fLookupTable->GetMinimum(1) );
1853  CreateLookupTable();
1854  return;
1855  }
1856 
1857  if (GetEPrime_Min() < fLookupTable->GetMinimum(2) || GetEPrime_Min() > fLookupTable->GetMaximum(2)) {
1858  SetEPrime_Min(fLookupTable->GetMinimum(2));
1859  G4cout << "!!!! Minimum E\' is out of bounds" << G4endl;
1860  G4cout << "---> Changing Minimum E\' to lower bound (" << fLookupTable->GetMinimum(2)/1000 << " GeV)" << G4endl;
1861  }
1862  if (GetEPrime_Max() > fLookupTable->GetMaximum(2) || GetEPrime_Max() < fLookupTable->GetMinimum(2)) {
1863  SetEPrime_Max(fLookupTable->GetMaximum(2));
1864  G4cout << "!!!! Maximum E\' is out of bounds" << G4endl;
1865  G4cout << "---> Changing Maximum E\' to upper bound (" << fLookupTable->GetMaximum(2)/1000 << " GeV)" << G4endl;
1866  }
1867  if (GetEPrime_Min() > GetEPrime_Max()) {
1868  SetEPrime_Min(GetEPrime_Max());
1869  G4cout << "!!!! Minimum E\' greater than maximum E\'" << G4endl;
1870  G4cout << "---> Changing Minimum E\' to " << GetEPrime_Max()/1000 << " GeV" << G4endl;
1871  }
1872 
1873  if (GetThetaAngle_Min() < fLookupTable->GetMinimum(3) || GetThetaAngle_Min() > fLookupTable->GetMaximum(3)) {
1874  SetThetaAngle_Min(fLookupTable->GetMinimum(3));
1875  G4cout << "!!!! Minimum Theta is out of bounds" << G4endl;
1876  G4cout << "---> Changing Minimum Theta to lower bound (" << fLookupTable->GetMinimum(3)/degree << " degrees)" << G4endl;
1877  }
1878  if (GetThetaAngle_Max() > fLookupTable->GetMaximum(3) || GetThetaAngle_Max() < fLookupTable->GetMinimum(3)) {
1879  SetThetaAngle_Max(fLookupTable->GetMaximum(3));
1880  G4cout << "!!!! Maximum Theta is out of bounds" << G4endl;
1881  G4cout << "---> Changing Maximum Theta to upper bound (" << fLookupTable->GetMaximum(3)/degree << " degrees)" << G4endl;
1882  }
1883  if (GetThetaAngle_Min() > GetThetaAngle_Max()) {
1884  SetThetaAngle_Min(GetThetaAngle_Max());
1885  G4cout << "!!!! Minimum Theta greater than maximum Theta" << G4endl;
1886  G4cout << "---> Changing Minimum Theta to " << GetThetaAngle_Max() << " degrees" << G4endl;
1887  }
1888  G4cout << "---> Leaving CheckLookupTableBounds() <---" << G4endl;
1889 }
1890 
1891 
1892 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1893 
1894 G4double QweakSimEPEvent::AlNuclInel(G4double E_in, // MeV
1895  G4double Theta, //radians
1896  G4double &fWeightN,
1897  G4double &Q2,
1898  G4double &E_out)
1899 {
1900  ///
1901  G4double M_p = 938.2796; // proton mass in MeV
1902 //G4double Z_Al = 13.0;
1903  G4double A_Al = 27.0;
1904  G4double M_Al = M_p*A_Al;
1905 
1906 //G4double CTH = cos(Theta/2.0);
1907  G4double STH = sin(Theta/2.0);
1908 
1909  // now get qsq
1910  G4double ETA = 1.0+2.0*E_in*STH*STH/M_Al;
1911  E_out = E_in/ETA;
1912 
1913  Q2 = 4*E_in*E_out*STH*STH; //unit: MeV^2
1914 
1915  Theta *= 180/3.14159; //deg
1916  // in ub/sr
1917  G4double xsect = Alinel_crsec_gen(E_in,Theta);
1918 
1919  Theta *= 3.14159/180; // bk to radians
1920 
1921  fWeightN = xsect*sin(Theta);
1922 
1923  return xsect;
1924 }
1925 
1926 //=============================================================================
1927 // Aluminum Giant Dipole Resonance Cross-Section
1928 // 2017/8/25 Original: Rupesh Silwal, Updated: Kurtis Bartlett
1929 // Makes use of Inelastic cross-section parameterization from J. Goldemberg et.
1930 // al, Nucl. Phys. 43 (1963) pg.242, also Al total photo-absorption GDR
1931 // measurement tabulated in "Atlas of Giant Dipole Resoances" IAEA Report
1932 // A.V. Varlamov et. al. January 1999
1933 //============================================================================
1934 G4double QweakSimEPEvent::AlGDR(G4double E_in, //MeV
1935  G4double Theta, // radians
1936  G4double &fWeightN,
1937  G4double &Q2,
1938  G4double &E_out){
1939 
1940  G4double M_p = 938.272081; // Proton mass [MeV] - 2016 PDG
1941  G4double Z_Al = 13.0;
1942  G4double N_Al = 14.0;
1943  G4double A_Al = 27.0;
1944  G4double M_Al = M_p*A_Al; // [MeV]
1945 
1946  G4double CTH = cos(Theta/2.0);
1947  G4double STH = sin(Theta/2.0);
1948 
1949  // now get qsq
1950  G4double ETA = 1.0+2.0*E_in*STH*STH/M_Al;
1951  E_out = E_in/ETA;
1952  Q2 = 4.0*E_in*E_out*STH*STH; // [MeV^2]
1953 
1954  // GDR peak energy
1955  G4double E_peak = 20.8; // [MeV]
1956  G4double fac = 9.0; // FWHM [MeV]
1957  //Use fac to scale xsection to reflect energy integration.
1958 
1959  G4double mu = ((N_Al*Z_Al)/A_Al)*M_Al; //Reduced mass [MeV]
1960 
1961  G4double long_ratio = pow(N_Al/A_Al,2.0)/(2.0*mu)*Q2/E_peak;
1962  G4double tran_ratio = 0.5*pow(N_Al/A_Al,2.0)*(E_peak/mu)*((1.0+STH*STH)/(CTH*CTH));
1963 
1964  G4double elasticAl = Horowitz_DW_Xsect(E_in,
1965  Theta,
1966  "horowitz_alloy_tables/elasticAl27.csv",
1967  fWeightN,
1968  Q2,
1969  E_out);
1970 
1971  G4double xsect = elasticAl*fac*(long_ratio + tran_ratio); // [ub/sr]
1972 
1973  fWeightN = xsect*sin(Theta);
1974 
1975  //G4cerr << "-----CrossSection: " << xsect << G4endl;
1976 
1977  return xsect;
1978 }
1979 
1980 //////////////////////////
1981 // Elastic e-12C Scattering
1982 // Marty McHugh
1983 // mjmchugh@jlab.org
1984 // Calculate the Elastic Carbon cross-section using PWBA Form Factors (FF)
1985 // deterimined using the method described in HANS F. EHRBNBERG et. all Phys. Rev. 113,2 (1959)
1986 G4double QweakSimEPEvent::Elastic_Cross_Section_Carbon( G4double E_in, // Incoming Energy [MeV]
1987  G4double Theta, // Scattering Angle [rad]
1988  G4double &fWeightN,
1989  G4double &Q2,
1990  G4double &E_out ) {
1991 
1992  // Minumum Scattering angle allowed
1993  G4double Theta_Min = 0.5*pi/180;
1994  if (Theta<Theta_Min) Theta = Theta_Min;
1995 
1996  // Physical Constants
1997  G4double M_A = 931.494 * MeV; // a.m.u. in MeV
1998  G4double Z = 6.0; // # of protons
1999  G4double A = 12.0; // Atomic Weight
2000  G4double M = M_A*A; // Nuclear Mass
2001  G4double a = 1.65; // well parameter [fm]
2002  G4double alpha_FS = 1.0/137.035999074; // Fine Structure Constant
2003  G4double m_e = 0.511 * MeV; // electron mass in MeV
2004 
2005  // Kinematic variables
2006  G4double CTH = cos(Theta/2.0);
2007  G4double STH = sin(Theta/2.0);
2008  G4double ETA = 1.0+2.0*E_in*STH*STH/M;
2009  E_out = E_in/ETA;
2010  Q2 = 4*E_in*E_out*STH*STH; // [MeV^2]
2011  G4double myhbarc = hbarc / MeV / fermi; // 197.3269631 MeV fm
2012  G4double q2 = Q2/(myhbarc*myhbarc); //convert MeV^2 into fm^(-2)
2013 
2014  // Dynamics
2015  // Mott Cross Section
2016  G4double SigmaMott = pow(alpha_FS*myhbarc*Z/(2.0*E_in),2); // [fm^2]
2017  SigmaMott *= (CTH*CTH)/pow(STH,4)*(E_out/E_in)*10000; // [ub] 10000 is the conversion factor
2018  G4double FF = (1.0 - (q2*a*a)/9.0)*exp(-(q2*a*a)/4.0); // Form Factor
2019  G4double sigma = SigmaMott*FF*FF;
2020  // Schwinger correction from Mo and Tsai Rev.Mod.Phys. 41 205 (1969) eq. (II.2)
2021  G4double delta_Schwinger = 0.0;
2022  if(SchwingerDeltaE!=0) {
2023  G4double FunctionofTheta = log (STH*STH) * log (CTH*CTH);
2024  delta_Schwinger = (-2.0*alpha_FS/pi) * ((log(E_in/SchwingerDeltaE) - 13.0/12.0)
2025  * (log(Q2/(m_e*m_e)) - 1.0) + 17.0/36.0 + FunctionofTheta/2.0);
2026  }
2027  sigma = sigma*(1.0 + delta_Schwinger);
2028  fWeightN = sigma*sin(Theta);
2029 
2030  return sigma;
2031 }
2032 
2033 //////////////////////////
2034 // 12C Giant Dipole Resonance (GDR)
2035 // Marty McHugh
2036 // Uses theory from J. Goldemberg et al, Nucl. Phys. 43 (1963) 242.
2037 // with peak energy from J. Goldemberg and W. C. Barber, Phys. Rev. 134, B963 (1964).
2038 G4double QweakSimEPEvent::CGDR( G4double E_in, // Incoming Energy [MeV]
2039  G4double Theta, // Scattering Angle [rad]
2040  G4double &fWeightN,
2041  G4double &Q2,
2042  G4double &E_out ) {
2043 
2044  // Nuclear Properties (Carbon-12)
2045  G4double M_A = 931.494 * MeV; // a.m.u. in [MeV]
2046  G4double Z = 6.0; // # of Protons
2047  G4double N = 6.0; // # of Neutrons
2048  G4double A = 12.0; // Atomic Weight
2049  G4double M = M_A*A; // Atomic Mass [MeV]
2050  G4double E_peak = 24.0*MeV; // GDR Peak Energy [MeV]
2051  G4double Energy_Integration_Factor = 6.5; // [MeV] FWHM of the GDR
2052  // G4double alpha_FS = 1.0/137.035999074; // Fine Structure Constant
2053  // G4double m_e = 0.511 * MeV; // electron mass in MeV
2054 
2055  // Kinematic Calculations
2056  G4double CTH = cos(Theta/2.0);
2057  G4double STH = sin(Theta/2.0);
2058  G4double ETA = 1.0+2.0*E_in*STH*STH/M;
2059  E_out = E_in/ETA-E_peak; // Outgoing electron energy [MeV]
2060  Q2 = 4*E_in*E_out*STH*STH; // [MeV^2]
2061 
2062  // Dynamics
2063  G4double LongRatio = (N*Q2)/(2*A*Z*M*E_peak); // Longitudinal component
2064  G4double TranRatio = (N*E_peak)/(2*A*Z*M)*(1+STH*STH)/(CTH*CTH); // Transverse component
2065  G4double sigma_el = AlloyScattering(E_in, Theta, Z, A, fWeightN, Q2, E_out);
2066  G4double sigma = (LongRatio+TranRatio)*sigma_el*Energy_Integration_Factor; // [ub/sr]
2067 
2068  // Commented out because this is possibly double counting.
2069  // There is allready a Schwinger correction in AlloyScattering.
2070  /* "Schwinger"-like correction for soft photon radiation: from Meister and Griffy PhysRev.133.B1032 (1964)
2071  G4double delta_Schwinger = 0.0;
2072  if(SchwingerDeltaE!=0) {
2073  delta_Schwinger = (alpha_FS/pi)*( log(SchwingerDeltaE*SchwingerDeltaE/(E_in*E_out) ) * ( log(Q2/(m_e*m_e)) - 1.0 )
2074  -0.5*log(E_in/E_out)*log(E_in/E_out) + 13.0/6.0*log(Q2/(m_e*m_e))
2075  - 28.0/9.0);
2076  }
2077  sigma = sigma*(1.0 + delta_Schwinger);
2078  */
2079  fWeightN = sigma*sin(Theta);
2080 
2081  return sigma;
2082 }
2083 
2084 //////////////////////////
2085 // 12C Inelastic Nuclear Excited States
2086 // Marty McHugh
2087 // States included:
2088 // 1 - 4.439 MeV, 2+
2089 // 2 - 7.543 MeV, 0+
2090 // 3 - 9.641 MeV, 3-
2091 // FF^2 calculated using gaussian fits to published FF^2 data.
2092 G4double QweakSimEPEvent::CNuclInel(G4double E_in, // MeV
2093  G4double Theta, //radians
2094  G4int nState,
2095  G4double &fWeightN,
2096  G4double &Q2,
2097  G4double &E_out)
2098 {
2099  // Physical Constants
2100  G4double M_A = 931.494 * MeV; // a.m.u. in MeV
2101  G4double Z = 6.0; // # of protons
2102  G4double A = 12.0; // Atomic Weight
2103  G4double M = M_A*A; // Nuclear Mass
2104  G4double alpha_FS = 1.0/137.035999074; // Fine Structure Constant
2105  G4double myhbarc = hbarc / MeV / fermi; // 197.3269631 MeV fm
2106  G4double m_e = 0.511 * MeV; // electron mass in MeV
2107  G4double EState[3] = { 4.439*MeV, // Energy Levels
2108  7.543*MeV,
2109  9.641*MeV};
2110 
2111  // Kinematic Calculations
2112  G4double CTH = cos(Theta/2.0);
2113  G4double STH = sin(Theta/2.0);
2114  G4double ETA = 1.0+2.0*E_in*STH*STH/M;
2115  E_out = E_in/ETA-EState[nState-1]; // Outgoing electron energy [MeV] subtracted excitation energy
2116  Q2 = 4*E_in*E_out*STH*STH; // [MeV^2]
2117  G4double q = sqrt(Q2/(myhbarc*myhbarc)); // [fm^-1]
2118 
2119  // Mott Cross Section
2120  G4double SigmaMott = pow(alpha_FS*myhbarc*Z/(2.0*E_in),2); // [fm^2]
2121  SigmaMott *= (CTH*CTH)/pow(STH,4)*(E_out/E_in)*10000; // [ub] 10000 is the conversion factor
2122 
2123  // FF^2 only from first three excited states (others are at least 1 order of magnitude smaller)
2124  G4double FF2;
2125  // Fits to published data. Descrbed in Ancillary Elog 257
2126  if(nState==1) { // 4.44 MeV, 2+: PhysRevLett27.745 and PhysRev.136.B1580
2127  TF1* ff2 = new TF1("4439","0.0147*exp(-3.25*(x-1.15)**2)",0.5,1.8);
2128  FF2 = ff2->Eval(q);
2129  delete ff2;
2130  } else if(nState==2) { // 7.65 Mev, 0+: PhysRev.136.B1580
2131  TF1* ff2 = new TF1("7654","0.0358*exp(-5.00*(x-0.96)**2)",0.5,1.8);
2132  FF2 = ff2->Eval(q);
2133  delete ff2;
2134  } else if(nState==3) { // 9.64 MeV, 3-: PhysRev.136.B1580
2135  TF1* ff2 = new TF1("9641","0.0526*exp(-2.45*(x-1.25)**2)",0.5,1.8);
2136  FF2 = ff2->Eval(q);
2137  delete ff2;
2138  } else {
2139  G4cout << "Please select from the first three states. nState = {1,2,3} only" << G4endl;
2140  return 0;
2141  }
2142 
2143  G4double sigma = SigmaMott*FF2;
2144 
2145  // "Schwinger"-like correction for soft photon radiation: from Meister and Griffy PhysRev.133.B1032 (1964)
2146  G4double delta_Schwinger = 0.0;
2147  if(SchwingerDeltaE!=0) {
2148  delta_Schwinger = (alpha_FS/pi)*( log(SchwingerDeltaE*SchwingerDeltaE/(E_in*E_out) ) * ( log(Q2/(m_e*m_e)) - 1.0 )
2149  -0.5*log(E_in/E_out)*log(E_in/E_out) + 13.0/6.0*log(Q2/(m_e*m_e))
2150  - 28.0/9.0);
2151  }
2152  sigma = sigma*(1.0 + delta_Schwinger);
2153 
2154  fWeightN = sigma*sin(Theta);
2155  return sigma;
2156 }
2157 
2158 /////// --------------------------------------------------------------------
2159 //
2160 //--- Fang Guo
2161 //
2162 //--- Calculates the Cross Section weighting factor for
2163 //--- pion photo-production with the wiser code.
2164 //--- LH2 target
2165 //
2166 // EPrime_Max, EPrime_Min in GeV
2167 // E_in is in MeV,
2168 //
2169 G4double QweakSimEPEvent::Pion_PhotoProduction(G4double E_in, //MeV
2170  G4double Theta, //radians
2171  G4double &fWeightN,
2172  G4double &Q2,
2173  G4double &E_out)
2174 {
2175  const G4double Mpi = 139.57018*MeV;
2176 
2177  if (Theta<GetThetaAngle_Min())
2178  {
2179  Theta = GetThetaAngle_Min();
2180  }
2181 
2182  if(E_in < Mpi) E_in = Mpi;
2183  E_out = Mpi + G4UniformRand()*(E_in - Mpi);
2184  G4double pf = sqrt(pow(E_out,2.0) - pow(Mpi,2.0)); //--- final momentum in MeV
2185 
2186  // Calculate Q2
2187  G4double STH = sin(Theta/2.0);
2188  Q2 = 4.0*E_in*E_out*STH*STH; // [MeV^2]
2189 
2190  //---
2191  //--- radiation length from page 12 of "The Qweak target design and safety document"
2192  //--- http://qweak.jlab.org/DocDB/0010/001041/002/Qweak%20Target%20PDR.pdf
2193  //---
2194  G4double lh2_length = myPositionZ - (myUserInfo->TargetCenterPositionZ) + 0.5*(myUserInfo->TargetLength);
2195  G4double intern_rad_len = 0.0204; // internal radiator from Mo/Tsai
2196  G4double rad_len = myUserInfo->TargetEntranceWindowThickness/(8.896*cm) + lh2_length/(871.9*cm) + intern_rad_len; //--- radiation length
2197  //G4cout << "radiation length: " << rad_len << G4endl;
2198 
2199  G4int type = 1; //--- 0: pi+, 1: pi-
2200 
2201  // wiser_sigma() needs E_in & pf in GeV, Theta in radians,
2202  // rad_len in fraction (not %)
2203  // wiser_sigma returns crsec in nb/GeV/sr
2204  G4double xsect = (1.0/1000.00) * wiser_sigma(E_in/GeV, pf/GeV, Theta, rad_len, type); //--- nanobarns/GeV/str --> ub/GeV/sr
2205 
2206  fWeightN = xsect*sin(Theta);
2207 
2208  // G4cout << E_in/GeV <<"\t" << pf/GeV << "\t" << Theta <<"\t" << rad_len <<"\t" << type << G4endl;
2209  // G4cout << "-- ** Pion crsec :: " << xsect << " ** --" << G4endl;
2210 
2211  return xsect;
2212 }
2213 
2214 //--- Al Dummy Target
2215 
2216 G4double QweakSimEPEvent::Pion_PhotoProductionAl(G4double E_in, //MeV
2217  G4double Theta, //radians
2218  G4double &fWeightN,
2219  G4double &Q2,
2220  G4double &E_out)
2221 {
2222  // Pi+- mass - PDG2014
2223  const G4double Mpi = 139.57018*MeV;
2224 
2225  if (Theta<GetThetaAngle_Min())
2226  {
2227  Theta = GetThetaAngle_Min();
2228  }
2229 
2230  //G4double tmp_EpMax = GetEPrime_Max();
2231  //if(GetEPrime_Max() > E_in) tmp_EpMax = E_in;
2232 
2233  //E_out = (G4UniformRand()*(tmp_EpMax - GetEPrime_Min()) + GetEPrime_Min()); //--- final total energy in MeV
2234 
2235  if(E_in < Mpi) E_in = Mpi;
2236  E_out = Mpi + G4UniformRand()*(E_in - Mpi);
2237  G4double pf = sqrt(pow(E_out,2.0) - pow(Mpi,2.0)); //--- final momentum in MeV
2238 
2239  // Calculate Q2
2240  G4double STH = sin(Theta/2.0);
2241  Q2 = 4.0*E_in*E_out*STH*STH; // [MeV^2]
2242 
2243  //---
2244  //--- radiation length from page 12 of "The Qweak target design and safety document"
2245  //--- http://qweak.jlab.org/DocDB/0010/001041/002/Qweak%20Target%20PDR.pdf
2246  //---
2247 
2248  // The below was original code from Fang, which is incorrect, I believe (D. Armstrong)
2249  /*
2250  G4double Al_length = 0.0*mm;
2251  if (ReactionRegion == 2 || ReactionRegion == 4 || ReactionRegion == 5 || ReactionRegion == 6 || ReactionRegion == 10) //-- Entrance & US Al dummy
2252  Al_length = myPositionZ - (myUserInfo->TargetCenterPositionZ) + 0.5*(myUserInfo->TargetLength) + (myUserInfo->TargetEntranceWindowThickness);
2253  else if (ReactionRegion == 3 || ReactionRegion == 7 || ReactionRegion == 8 || ReactionRegion == 9 || ReactionRegion == 11) //-- Exit & DS Al dummy
2254  Al_length = myPositionZ - (myUserInfo->TargetCenterPositionZ) - 0.5*(myUserInfo->TargetLength);
2255 
2256  G4double intern_rad_len = 0.0204; // internal radiator from Mo/Tsai
2257  G4double rad_len = Al_length/(8.896*cm) + intern_rad_len; //--- radiation length
2258  */
2259 
2260  // Need to calculate radiation length upstream of the scattering vertex, in order for Wiser code to estimate bremss spectrum for
2261  // the pion photoproduction. Since the Al windows are relatively thin, we just assume scattering takes place halfway along
2262  // the Al material, on average.
2263 
2264  G4double Al_length = 0.0*mm;
2265  G4double lh2_length = 0.0*mm;
2266  if (ReactionRegion == 2 ) //-- Entrance window
2267  {
2268  Al_length = 0.5*(myUserInfo->TargetEntranceWindowThickness);
2269  }
2270  else if (ReactionRegion == 3 ) //-- Exit window
2271  {
2272  Al_length = (myUserInfo->TargetEntranceWindowThickness) + 0.5*(myUserInfo->TargetExitWindowNippleThickness);
2273  lh2_length = myUserInfo->TargetLength;
2274  }
2275  else if (ReactionRegion == 8 ) //-- 4% DS window No entrance window, no hydrogen here...
2276  {
2277  Al_length = 0.5*(myUserInfo->TargetThicknessDSALDummy4);
2278  }
2279 
2280  G4double intern_rad_len = 0.0204; // internal radiator from Mo/Tsai
2281  G4double rad_len = Al_length/(8.896*cm) + intern_rad_len + lh2_length/(866.0*cm); //--- radiation length in fractional terms
2282 
2283  //G4cout << "ReactionRegion: " << ReactionRegion << G4endl;
2284  //G4cout << "myPositionZ: " << myPositionZ << G4endl;
2285  //G4cout << "myUserInfo->TargetCenterPositionZ: " << myUserInfo->TargetCenterPositionZ << G4endl;
2286  //G4cout << "myUserInfo->TargetLength: " << myUserInfo->TargetLength << G4endl;
2287  //G4cout << "Al_length: " << Al_length << G4endl;
2288  //G4cout << "lh2_length: " << lh2_length << G4endl;
2289  //G4cout << "rad_len: " << rad_len << G4endl;
2290  //G4cout << "Theta: " << Theta << G4endl;
2291 
2292  G4int type = 1; //--- 0: pi+ 1:pi-
2293 
2294  // wiser_sigma() needs E_in & pf in GeV, Theta in radians,
2295  // rad_len in fraction (not %)
2296  // wiser_sigma returns crsec in nb/GeV/sr
2297  //
2298  // For a nucleus, need to get cross section from the neutrons and from the protons.
2299  // To generate a pi^- from the proton, need two-pion production (can generate pi^- from neutron via single pion production).
2300  // assume exact isospin symmetry, so to get pi^- cross section for neutron, use pi^+ cross section for proton, which is what
2301  // is in the Wiser code fits.
2302  // Ignore FSI or other nuclear effects (treat 27Al as a bag of 13 free protons and 14 free neutrons)
2303 
2304  G4double xsect_proton = (1.0/1000.00) * wiser_sigma(E_in/GeV, pf/GeV, Theta, rad_len, type); //--- nanobarns/GeV/str --> ub/GeV/sr
2305 
2306  // below: switch type to 0 to get the pi+ cross section, which we will use for the pi^- cross section from neutrons
2307  G4double xsect_neutron = (1.0/1000.00) * wiser_sigma(E_in/GeV, pf/GeV, Theta, rad_len, 0); //--- nanobarns/GeV/str --> ub/GeV/sr
2308  G4double xsect = 13.0*xsect_proton + 14.0*xsect_neutron;
2309  fWeightN = xsect*sin(Theta);
2310 
2311  //G4cout << GetEPrime_Max() << "\t" << GetEPrime_Min() << "\t" << E_in/GeV <<"\t" << E_out/GeV << "\t" << pf/GeV << "\t" << Q2/GeV/GeV << "\t" << Theta <<"\t" << rad_len <<"\t" << type << G4endl;
2312  //G4cout << "------ Pion crsec :: " << xsect << " ------" << G4endl;
2313 
2314  return xsect;
2315 }
2316 
2317 
2318 //--- Carbon Target
2319 
2320 G4double QweakSimEPEvent::Pion_PhotoProductionCarbon(G4double E_in, //MeV
2321  G4double Theta, //radians
2322  G4double &fWeightN,
2323  G4double &Q2,
2324  G4double &E_out)
2325 {
2326  const G4double Mpi = 0.13957*GeV; ///--- GeV
2327 
2328  if (Theta<GetThetaAngle_Min())
2329  Theta = GetThetaAngle_Min();
2330 
2331  // need to evaluate Q2 properply, in order to evaluate
2332  // the internal radiation length more accurately below
2333  Q2 = 0.026;
2334 
2335  G4double tmp_EpMax = GetEPrime_Max();
2336  if(GetEPrime_Max() > E_in) tmp_EpMax = E_in;
2337 
2338  E_out = (G4UniformRand()*(tmp_EpMax - GetEPrime_Min()) + GetEPrime_Min()); //--- final total energy in GeV
2339  G4double pf = sqrt(pow(E_out,2) - pow(Mpi,2)); //--- final momentum in GeV
2340 
2341  //---
2342  //--- radiation length calculated from http://pdg.lbl.gov/2013/reviews/rpp2013-rev-atomic-nuclear-prop.pdf
2343  //---
2344  //G4double C_length = myPositionZ - (myUserInfo->TargetCenterPositionZ) + 0.5*(myUserInfo->TargetLength);
2345  G4double C_length = 0.0*mm;
2346  if (ReactionRegion == 2 || ReactionRegion == 4 || ReactionRegion == 5 || ReactionRegion == 6 || ReactionRegion == 10) //-- Entrance & US C dummy
2347  C_length = myPositionZ - (myUserInfo->TargetCenterPositionZ) + 0.5*(myUserInfo->TargetLength) + (myUserInfo->TargetEntranceWindowThickness);
2348  else if (ReactionRegion == 3 || ReactionRegion == 7 || ReactionRegion == 8 || ReactionRegion == 9 || ReactionRegion == 11) //-- Exit & DS C dummy
2349  C_length = myPositionZ - (myUserInfo->TargetCenterPositionZ) - 0.5*(myUserInfo->TargetLength);
2350 
2351  G4double intern_rad_len = 0.0204; // internal radiator from Mo/Tsai
2352  G4double rad_len = C_length/(19.32*cm) + intern_rad_len; //--- radiation length
2353 
2354  G4int type = 1; //--- 0: pi+ 1:pi-
2355 
2356  // wiser_sigma() needs E_in & pf in GeV, Theta in radians,
2357  // rad_len in fraction (not %)
2358  // wiser_sigma returns crsec in nb/GeV/sr
2359  G4double xsect = (1.0/1000.00) * wiser_sigma(E_in/GeV, pf/GeV, Theta, rad_len, type); //--- nanobarns/GeV/str --> ub/GeV/sr
2360 
2361  fWeightN = xsect*sin(Theta);
2362 
2363  //G4cout << E_in/GeV <<"\t" << pf/GeV << "\t" << Theta <<"\t" << rad_len <<"\t" << type << G4endl;
2364  //G4cout << "------ Pion crsec :: " << xsect << " ------" << G4endl;
2365 
2366  return xsect;
2367 }
2368 
2369 //==============================================================================
2370 // Aluminum Nuclear Excited State Generator
2371 // Added 2017/5/16
2372 // Kurtis Bartlett kbartlet@jlab.org
2373 // Fit to MIT Bates and Glasgow Electroexcitation of discrete excited states.
2374 // See Qweak ELOG 1671
2375 //==============================================================================
2376 
2377 G4double QweakSimEPEvent::Aluminum_Excited_State(G4double E_in, //[MeV]
2378  G4double Theta, //[radians]
2379  G4double E_lvl, //[MeV]
2380  G4double fit_c, //[unitless]
2381  G4double fit_mean, //[fm^-1]
2382  G4double fit_sigma, //[fm^-1]
2383  G4double &fWeightN,
2384  G4double &Q2,
2385  G4double &E_out){
2386  G4double xsect = 0.0;
2387  G4double FF2 = 0.0;
2388  const G4double alpha = 1/137.035999139; // Fine Structure Constant 2016-PDG
2389  const G4double M_p = 938.272081; // Proton Mass [MeV] 2016-PDG
2390  const G4double M_n = 939.565413; // Neutron Mass [MeV] 2016-PDG
2391  const G4double m_e = 0.5109989461*MeV; // Electron mass 2016 PDG
2392  const G4double Z_Al = 13.0; // # of Protons
2393  const G4double A_Al = 27.0; // Atomic Mass [amu]
2394  const G4double N_Al = A_Al-Z_Al; // # of Neutrons
2395  const G4double M_Al = M_p*Z_Al + M_n*N_Al; // Atomic Mass [MeV]
2396  const G4double hbar_c = 197.3269788; // Conversion Constant [MeV fm]
2397  const G4double R = 3.94;// Al Uniform-density charge radius [fm]
2398  const G4double deltaE = 20.0*MeV; // Schwinger correction delta E parameter
2399 
2400  G4double STH = sin(Theta/2.0);
2401  G4double CTH = cos(Theta/2.0);
2402 
2403  //Calculate Recoil Factor Eta
2404  G4double ETA = 1.0 + 2.0*(E_in/M_Al)*STH*STH;
2405 
2406  //Calculate Eprime
2407  //See H. Uberall, "Electron Scattering from Complex Nuclei",
2408  //(Academic, NY, 1971) eq.2-10a and eq.2-10d
2409  E_out = (E_in - E_lvl - ((E_lvl*E_lvl)/(2.0*M_Al)))/ETA; //[MeV]
2410 
2411  //Calculate/Return Q^2 and convert to q [fm^-1]
2412  Q2 = 4.0*E_in*E_out*STH*STH; //[MeV^2]
2413  G4double q = sqrt(Q2/(hbar_c*hbar_c)); //[fm^-1]
2414 
2415  //Calculate Q_eff (Coulomb Corrections)
2416  //See Ryan et al. Phys.Rev.C27, 2515 (1983)
2417  G4double q_eff = q*(1.0 + (3.0*Z_Al*alpha*hbar_c)/(2.0*E_in*R)); //[fm^-1]
2418 
2419  //Calculate Mott Cross-section [fm^2/sr]
2420  G4double sigma_mott = pow((Z_Al*alpha*hbar_c*CTH)/(2.0*E_in*STH*STH),2.0)/ETA;
2421  //Convert from [fm^2/sr] to [ubarn/sr]
2422  sigma_mott = 10000.0*sigma_mott;
2423 
2424  //Calculate Gaussian Parameterization of Form Factor Squared
2425  //using q_eff in range 0.3-1.8 [fm^-1]
2426  if((q_eff >= 0.3) && (q_eff <= 1.8)){
2427  FF2 = fit_c*exp(-0.5*pow((q_eff-fit_mean)/(fit_sigma),2.0));
2428  } else{
2429  FF2 = 0.0;
2430  }
2431 
2432  //Calculate Cross-section
2433  xsect = sigma_mott*FF2;
2434 
2435  // Calculate Schwinger Radiative Correction
2436  // Using eq.6-8d pg.488 from "Electron Scattering from Complex Nuclei"
2437  // Part B Herber Uberall - 2017/8/24 K. Bartlett
2438  G4double delta_Schwinger = 0.0;
2439  delta_Schwinger = (TMath::Log(TMath::Sqrt(E_in*E_out)/deltaE)-(13.0/12.0))*(TMath::Log(Q2/(m_e*m_e))-1.0);
2440  delta_Schwinger += (17.0/36.0) + 0.25*TMath::Log(E_in/E_out)*TMath::Log(E_in/E_out) + 0.5*((1.0/6.0)*TMath::Pi()*TMath::Pi()-TMath::DiLog(TMath::Cos(Theta/2.0)*TMath::Cos(Theta/2.0)));
2441  delta_Schwinger *= (2.0*alpha/TMath::Pi());
2442 
2443  xsect *= TMath::Exp(-delta_Schwinger);
2444 
2445 
2446  //Return weighted Cross-Section
2447  fWeightN = xsect*sin(Theta);
2448 
2449  //Debug statement
2450  //G4cout << "E_out: " << E_out << " Q2: " << Q2 << " q_eff: " << q_eff << " FF2: " << FF2 << " Xsect: " << xsect << G4endl;
2451 
2452  return xsect;
2453 }
2454 
2455 //==============================================================================
2456 // Horowitz Distorted Wave Aluminum Elastic Cross-section and Asymmetry
2457 // calculation. Uses Eval method from Cern ROOT TGraph to interpolate from
2458 // calculated cross-section and asymmetry data set from his 2014 publication.
2459 // See Phys. Rev. C89 045503 (2014) for more info about his calculation.
2460 // Note this is for a fixed beam energy of 1.160 GeV.
2461 // Added 2017/7/19
2462 // Kurtis Bartlett kbartlet@jlab.org
2463 //==============================================================================
2464 G4double QweakSimEPEvent::Horowitz_DW_Xsect(G4double E_in, //[MeV]
2465  G4double theta, //[radians]
2466  const G4String path,
2467  G4double &fWeightN,
2468  G4double &Q2,
2469  G4double &E_out)
2470 {
2471  const G4double M = 13.0*938.272081 + 14.0*939.565413; //2016 PDG - Al mass [MeV]
2472  const G4double alpha = 1.0/137.035999139; // Fine Structure Const 2016 PDG
2473  const G4double deltaE = 20.0*MeV; // Schwinger correction delta E parameter
2474  const G4double m_e = 0.5109989461*MeV; // Electron mass 2016 PDG
2475  // Calculate kinematic variables
2476  G4double STH = sin(theta/2.0);
2477  G4double ETA = 1.0+2.0*(E_in/M)*STH*STH;
2478  E_out = E_in/ETA; //E prime [MeV]
2479  Q2 = 4.0*E_in*E_out*STH*STH; //[MeV^2]
2480 
2481  // Interpolate Horowitz Xsect calculation
2482  G4double radtodeg = (180.0/TMath::Pi());
2483  G4double xsect = 0.0;
2484  TGraph *graph = new TGraph(path, "%lg %lg %*lg", ",");
2485  graph->SetBit(graph->kIsSortedX); // Faster queires on TGraph object.
2486  // Check theta value inside of table range (4.0-13.5 [deg])
2487  if((theta*radtodeg) >= 4.0 && (theta*radtodeg) <= 13.5)
2488  {
2489  xsect = graph->Eval((theta*radtodeg), 0, "S"); // Spline Interpolation
2490  }
2491  // Calculate Schwinger Radiative Correction
2492  // Using eq.6-8c pg.487 from "Electron Scattering from Complex Nuclei"
2493  // Part B Herbert Uberall - 2017/7/26 K. Bartlett
2494  G4double delta_Schwinger = 0.0;
2495  delta_Schwinger = (TMath::Log(E_in/deltaE)-(13.0/12.0))*(TMath::Log(Q2/(m_e*m_e))-1.0);
2496  delta_Schwinger += (17.0/36.0) + 0.5*((1.0/6.0)*TMath::Pi()*TMath::Pi()-TMath::DiLog(TMath::Cos(theta/2.0)*TMath::Cos(theta/2.0)));
2497  delta_Schwinger *= (2.0*alpha/TMath::Pi());
2498 
2499  xsect *= TMath::Exp(-delta_Schwinger)*1000.0; //Convert from [mb/sr] to [ub/sr]
2500 
2501  fWeightN = xsect*sin(theta);
2502  delete graph; // Delete TGraph object.
2503  return xsect;
2504 }
2505 
2506 G4double QweakSimEPEvent::Horowitz_DW_Asym(G4double theta, // [radians]
2507  const G4String path)
2508 {
2509  G4double radtodeg = (180.0/3.1415927);
2510  G4double asym = 0.0; // Asymmetry
2511  TGraph *graph = new TGraph(path, "%lg %*lg %lg", ",");
2512  graph->SetBit(graph->kIsSortedX); // Faster queires on TGraph object.
2513  if((theta*radtodeg) >= 4.0 && (theta*radtodeg) <= 13.5)
2514  {
2515  asym = graph->Eval((theta*radtodeg), 0, "S"); // Spline Interpolation
2516  }
2517  delete graph; // Delete TGraph object.
2518  return asym*1e6; // [ppm]
2519 }
2520 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
2521 
2522 // Quasielastic fit from Bosted et al
2523 // Added 19/March kamyers@jlab.org
2524 
2525 G4double QweakSimEPEvent::Quasi_Elastic_Bosted(G4double E_in, //MeV
2526  G4double Theta,
2527  G4int Zin,
2528  G4int Ain,
2529  G4double &fWeightN,
2530  G4double &Q2,
2531  G4double &E_out)
2532 {
2533 
2534  G4double F1 = 0.0;
2535  G4double F2 = 0.0;
2536  G4double W1 = 0.0;
2537  G4double W2 = 0.0;
2538  G4double Wsq = 0.0;
2539  G4double xsect = 0.0;
2540  G4double Ein = 0.0;
2541  G4double Eout = 0.0;
2542 
2543  // In this subroutine units are in GeV, so
2544  Ein = E_in/1000.; // maximum energy in GeV
2545  G4double M_electron = 0.000511; // minimum energy (electron mass) in GeV
2546  G4double Mp = 0.93828;
2547 
2548  // G4double Theta_Min = 1.745329E-4;
2549  if (Theta<GetThetaAngle_Min())
2550  Theta = GetThetaAngle_Min();
2551 
2552  // Generate flat energy distribution of outgoing electron
2553  Eout = M_electron + G4UniformRand()*(Ein - M_electron);
2554 
2555  G4double CTH = cos(Theta/2.0);
2556  G4double STH = sin(Theta/2.0);
2557  G4double T2THE = STH*STH/CTH/CTH;
2558  G4double Nu = Ein - Eout;
2559  Q2 = 4.0*Ein*Eout*STH*STH;
2560  Wsq = Mp*Mp + 2.0*Mp*Nu - Q2;
2561 
2562  // Mott scattering
2563  G4double MOTT = pow((0.00072/Ein*CTH/STH/STH),2);
2564  MOTT = MOTT*1.0E4; // Units: ub/sr/GeV
2565 
2566  //G4cout << "Ein, Eout, Nu, Theta: " << Ein << ", " << Eout << ", " << Nu << ", " << Theta/degree << G4endl;
2567  //G4cout << "Q2, W2, Mott: " << Q2 << ", " << Wsq << ", " << MOTT << G4endl;
2568 
2569  F1F2QE09(Zin, Ain, Q2, Wsq, F1, F2);
2570  W1 = F1/Mp;
2571  W2 = F2/Nu;
2572 
2573  xsect = MOTT*(W2 + 2.0*T2THE*W1)*(Ein - M_electron);
2574  // G4cout << "F1, F2, W1, W2: " << F1 << ", " << F2 << ", " << W1 << ", " << W2 << G4endl;
2575 
2576  // In some cases a negative F2 is returned giving a negative cross section
2577  if (xsect <= 0) xsect = 0.0;
2578 
2579  // G4cout << "xsect: " << xsect << G4endl;
2580 
2581  fWeightN = xsect*sin(Theta);
2582 
2583  E_out = Eout*1000.; // Convert back to MeV before leaving
2584  Q2 = Q2*1.0E6; // Seems like this needs to be returned in MeV**2
2585  return xsect;
2586 }
2587 
2588 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
2589 
2590 // Nuclear inelastic fit from Bosted et al
2591 // Added 14/May/2014 kamyers@jlab.org
2592 
2593 G4double QweakSimEPEvent::NuclearInelastic_Bosted(G4double E_in, //MeV
2594  G4double Theta, //radians
2595  G4int Zin,
2596  G4int Ain,
2597  G4double &fWeightN,
2598  G4double &Q2,
2599  G4double &E_out)
2600 {
2601 
2602  G4double F1 = 0.0;
2603  G4double F2 = 0.0;
2604  G4double W1 = 0.0;
2605  G4double W2 = 0.0;
2606  G4double Wsq = 0.0;
2607  G4double xsect = 0.0;
2608  G4double Ein = 0.0;
2609  G4double Eout = 0.0;
2610 
2611  // In this subroutine units are in GeV, so
2612  Ein = E_in/1000.; // maximum energy in GeV
2613  G4double M_electron = 0.000511; // minimum energy (electron mass) in GeV
2614  G4double Mp = 0.93828;
2615 
2616  // G4double Theta_Min = 1.745329E-4;
2617  if (Theta<GetThetaAngle_Min())
2618  Theta = GetThetaAngle_Min();
2619 
2620  // Generate flat energy distribution of outgoing electron
2621  Eout = M_electron + G4UniformRand()*(Ein - M_electron);
2622 
2623  G4double CTH = cos(Theta/2.0);
2624  G4double STH = sin(Theta/2.0);
2625  G4double T2THE = STH*STH/CTH/CTH;
2626  G4double Nu = Ein - Eout;
2627  Q2 = 4.0*Ein*Eout*STH*STH;
2628  Wsq = Mp*Mp + 2.0*Mp*Nu - Q2;
2629 
2630  // Mott scattering
2631  G4double MOTT = pow((0.00072/Ein*CTH/STH/STH),2);
2632  MOTT = MOTT*1.0E4; // Units: ub/sr/GeV
2633 
2634  //G4cout << "Ein, Eout, Nu, Theta: " << Ein << ", " << Eout << ", " << Nu << ", " << Theta/degree << G4endl;
2635  //G4cout << "Q2, W2, Mott: " << Q2 << ", " << Wsq << ", " << MOTT << G4endl;
2636 
2637  F1F2IN09(Zin, Ain, Q2, Wsq, F1, F2);
2638  W1 = F1/Mp;
2639  W2 = F2/Nu;
2640 
2641  xsect = MOTT*(W2 + 2.0*T2THE*W1)*(Ein - M_electron);
2642  //G4cout << "F1, F2, W1, W2: " << F1 << ", " << F2 << ", " << W1 << ", " << W2 << G4endl;
2643 
2644  //G4cout << E_in << "\t" << Theta << "\t" << MOTT*(Ein - M_electron)*Zin*Zin << "\t" << (W2 + 2.0*T2THE*W1)/(Zin*Zin) << G4endl;
2645 
2646  // In some cases a negative F2 is returned giving a negative cross section
2647  if (xsect <= 0) xsect = 0.0;
2648 
2649  fWeightN = xsect*sin(Theta);
2650  //G4cout << "xsect, fweight: " << xsect << ", "<< fWeightN << G4endl;
2651 
2652  E_out = Eout*1000.; // Convert back to MeV before leaving
2653  Q2 = Q2*1.0E6; // Seems like this needs to be returned in MeV**2
2654  return xsect;
2655 }
2656 
2657 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
2658 
2659 void QweakSimEPEvent::F1F2QE09(G4int Z, G4int IA, G4double QSQ,
2660  G4double wsq, G4double &F1, G4double &F2)
2661 {
2662 //=======================================================================
2663 // SUBROUTINE F1F2QE09(Z, A, QSQ, wsq, F1, F2)
2664 //
2665 // Calculates quasielastic A(e,e')X structure functions F1 and F2 PER NUCLEUS
2666 // for A>2 uses superscaling from Sick, Donnelly, Maieron, nucl-th/0109032
2667 // for A=2 uses pre-integrated Paris wave function (see ~bosted/smear.f)
2668 // coded by P. Bosted August to October, 2006
2669 //
2670 // input: Z, A (real*8) Z and A of nucleus (shoud be 2.0D0 for deueron)
2671 // Qsq (real*8) is 4-vector momentum transfer squared (positive in
2672 // chosen metric)
2673 // Wsq (real*8) is invarinat mass squared of final state calculated
2674 // assuming electron scattered from a free proton
2675 //
2676 // outputs: F1, F2 (real*8) are structure functions per nucleus
2677 //
2678 // Note: Deuteron agrees well with Laget (see ~bosted/eg1b/laget.f) for
2679 // a) Q2<1 gev**2 and dsig > 1% of peak: doesnt describe tail at high W
2680 // b) Q2>1 gev**2 on wings of q.e. peak. But, this model is up
2681 // to 50% too big at top of q.e. peak. BUT, F2 DOES agree very
2682 // nicely with Osipenko et al data from CLAS, up to 5 GeV**2
2683 
2684 // G4cout << "Z, A, Q2, W2: " << Z << ", " << IA << ", " << QSQ << ", " << wsq << G4endl;
2685 
2686  G4double avgN, Pauli_sup1, Pauli_sup2, GEP, GEN, GMP, GMN, Q, Q3, Q4;
2687  G4double amp = 0.93828;
2688  G4double amd = 1.8756;
2689  G4double RMUP = 2.792782;
2690  G4double RMUN = -1.913148;
2691  G4double pz, Nu, dpz, pznom, pzmin;
2692  G4double QV, TAU, FY, dwmin, w2p;
2693  G4double kappa, lam, lamp, taup, squigglef, psi, psip, nuL, nuT;
2694  G4double kf, Es, GM2bar, GE2bar, W2bar, Delta, GL, GT;
2695  G4int izz, izzmin, izp, izznom, izdif;
2696 
2697  // Look up tables for deuteron case
2698  G4double fyd[200] = {
2699  0.00001,0.00002,0.00003,0.00005,0.00006,0.00009,0.00010,0.00013,
2700  0.00015,0.00019,0.00021,0.00026,0.00029,0.00034,0.00038,0.00044,
2701  0.00049,0.00057,0.00062,0.00071,0.00078,0.00089,0.00097,0.00109,
2702  0.00119,0.00134,0.00146,0.00161,0.00176,0.00195,0.00211,0.00232,
2703  0.00252,0.00276,0.00299,0.00326,0.00352,0.00383,0.00412,0.00447,
2704  0.00482,0.00521,0.00560,0.00603,0.00648,0.00698,0.00747,0.00803,
2705  0.00859,0.00921,0.00985,0.01056,0.01126,0.01205,0.01286,0.01376,
2706  0.01467,0.01569,0.01671,0.01793,0.01912,0.02049,0.02196,0.02356,
2707  0.02525,0.02723,0.02939,0.03179,0.03453,0.03764,0.04116,0.04533,
2708  0.05004,0.05565,0.06232,0.07015,0.07965,0.09093,0.10486,0.12185,
2709  0.14268,0.16860,0.20074,0.24129,0.29201,0.35713,0.44012,0.54757,
2710  0.68665,0.86965,1.11199,1.43242,1.86532,2.44703,3.22681,4.24972,
2711  5.54382,7.04016,8.48123,9.40627,9.40627,8.48123,7.04016,5.54382,
2712  4.24972,3.22681,2.44703,1.86532,1.43242,1.11199,0.86965,0.68665,
2713  0.54757,0.44012,0.35713,0.29201,0.24129,0.20074,0.16860,0.14268,
2714  0.12185,0.10486,0.09093,0.07965,0.07015,0.06232,0.05565,0.05004,
2715  0.04533,0.04116,0.03764,0.03453,0.03179,0.02939,0.02723,0.02525,
2716  0.02356,0.02196,0.02049,0.01912,0.01793,0.01671,0.01569,0.01467,
2717  0.01376,0.01286,0.01205,0.01126,0.01056,0.00985,0.00921,0.00859,
2718  0.00803,0.00747,0.00698,0.00648,0.00603,0.00560,0.00521,0.00482,
2719  0.00447,0.00412,0.00383,0.00352,0.00326,0.00299,0.00276,0.00252,
2720  0.00232,0.00211,0.00195,0.00176,0.00161,0.00146,0.00134,0.00119,
2721  0.00109,0.00097,0.00089,0.00078,0.00071,0.00062,0.00057,0.00049,
2722  0.00044,0.00038,0.00034,0.00029,0.00026,0.00021,0.00019,0.00015,
2723  0.00013,0.00010,0.00009,0.00006,0.00005,0.00003,0.00002,0.00001};
2724 
2725  G4double avp2[200] = {
2726  1.0,0.98974,0.96975,0.96768,0.94782,0.94450,0.92494,0.92047,
2727  0.90090,0.89563,0.87644,0.87018,0.85145,0.84434,0.82593,0.81841,
2728  0.80021,0.79212,0.77444,0.76553,0.74866,0.73945,0.72264,0.71343,
2729  0.69703,0.68740,0.67149,0.66182,0.64631,0.63630,0.62125,0.61154,
2730  0.59671,0.58686,0.57241,0.56283,0.54866,0.53889,0.52528,0.51581,
2731  0.50236,0.49291,0.47997,0.47063,0.45803,0.44867,0.43665,0.42744,
2732  0.41554,0.40656,0.39511,0.38589,0.37488,0.36611,0.35516,0.34647,
2733  0.33571,0.32704,0.31656,0.30783,0.29741,0.28870,0.27820,0.26945,
2734  0.25898,0.25010,0.23945,0.23023,0.21943,0.20999,0.19891,0.18911,
2735  0.17795,0.16793,0.15669,0.14667,0.13553,0.12569,0.11504,0.10550,
2736  0.09557,0.08674,0.07774,0.06974,0.06184,0.05484,0.04802,0.04203,
2737  0.03629,0.03129,0.02654,0.02247,0.01867,0.01545,0.01251,0.01015,
2738  0.00810,0.00664,0.00541,0.00512,0.00512,0.00541,0.00664,0.00810,
2739  0.01015,0.01251,0.01545,0.01867,0.02247,0.02654,0.03129,0.03629,
2740  0.04203,0.04802,0.05484,0.06184,0.06974,0.07774,0.08674,0.09557,
2741  0.10550,0.11504,0.12569,0.13553,0.14667,0.15669,0.16793,0.17795,
2742  0.18911,0.19891,0.20999,0.21943,0.23023,0.23945,0.25010,0.25898,
2743  0.26945,0.27820,0.28870,0.29741,0.30783,0.31656,0.32704,0.33571,
2744  0.34647,0.35516,0.36611,0.37488,0.38589,0.39511,0.40656,0.41554,
2745  0.42744,0.43665,0.44867,0.45803,0.47063,0.47997,0.49291,0.50236,
2746  0.51581,0.52528,0.53889,0.54866,0.56283,0.57241,0.58686,0.59671,
2747  0.61154,0.62125,0.63630,0.64631,0.66182,0.67149,0.68740,0.69703,
2748  0.71343,0.72264,0.73945,0.74866,0.76553,0.77444,0.79212,0.80021,
2749  0.81841,0.82593,0.84434,0.85145,0.87018,0.87644,0.89563,0.90090,
2750  0.92047,0.92494,0.94450,0.94782,0.96768,0.96975,0.98974,1.0};
2751 
2752  // Peter Bosted's correction params
2753  /*
2754  G4double pb[20] = {
2755  0.1023E+02, 0.1052E+01, 0.2485E-01, 0.1455E+01,
2756  0.5650E+01,-0.2889E+00, 0.4943E-01,-0.8183E-01,
2757  -0.7495E+00, 0.8426E+00,-0.2829E+01, 0.1607E+01,
2758  0.1733E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00,
2759  0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00};
2760  */
2761  G4double y,R;
2762 
2763  G4double P[24] = {
2764  5.1377e-03, 9.8071e-01, 4.6379e-02, 1.6433e+00,
2765  6.9826e+00, -2.2655e-01, 1.1095e-01, 2.7945e-02,
2766  4.0643e-01, 1.6076e+00, -7.5460e+00, 4.4418e+00,
2767  -3.7464e-01, 1.0414e-01, -2.6852e-01, 9.6653e-01,
2768  -1.9055e+00, 9.8965e-01, 2.0613e+02, -4.5536e-02,
2769  2.4902e-01, -1.3728e-01, 2.9201e+01, 4.9280e-03};
2770 
2771  // return if proton: future change this to allow for
2772  // equivalent W resolution
2773  F1 = 0.;
2774  F2 = 0.;
2775  avgN = IA - Z;
2776  if (IA==1) return;
2777 
2778  // some kinematic factors. Return if Nu or QSQ is negative
2779  Nu = (wsq - amp*amp + QSQ) / 2. / amp;
2780 
2781  //G4cout << "In call... IA, Nu, QSQ = " << IA << ", " << Nu << ", " << QSQ << G4endl;
2782  if(Nu <= 0.0 || QSQ < 0.) return;
2783  TAU = QSQ / 4.0 / amp / amp;
2784  QV = sqrt(Nu*Nu + QSQ);
2785 
2786  // Bosted fit for nucleon form factors Phys. Rev. C 51, p. 409 (1995)
2787  Q = sqrt(QSQ);
2788  Q3 = QSQ * Q;
2789  Q4 = QSQ*QSQ;
2790  GEP = 1./ (1. + 0.14 * Q + 3.01 * QSQ + 0.02 * Q3 + 1.20 * Q4 + 0.32 * pow(Q,5));
2791  GMP = RMUP * GEP;
2792  GMN = RMUN / (1.- 1.74 * Q + 9.29 * QSQ - 7.63 * Q3 + 4.63 * Q4);
2793  GEN = 1.25 * RMUN * TAU / (1. + 18.3 * TAU) / pow((1. + QSQ / 0.71),2);
2794 
2795  //G4cout << "Form Factors: " << GEP << ", " << GMP << ", " << GEN << ", " << GMN << G4endl;
2796 
2797  // Get kf and Es from superscaling from Sick, Donnelly, Maieron,
2798  // nucl-th/0109032
2799  if(IA==2) kf=0.085;
2800  if(IA==2) Es=0.0022;
2801  // changed 4/09
2802  if(IA==3) kf=0.115;
2803  if(IA==3) Es=0.001 ;
2804  // changed 4/09
2805  if(IA>3) kf=0.19;
2806  if(IA>3) Es=0.017;
2807  if(IA>7) kf=0.228;
2808  if(IA>7) Es=0.020;
2809  // changed 5/09
2810  if(IA>7) Es=0.0165;
2811  if(IA>16) kf=0.230;
2812  if(IA>16) Es=0.025;
2813  if(IA>25) kf=0.236;
2814  if(IA>25) Es=0.018;
2815  if(IA>38) kf=0.241;
2816  if(IA>38) Es=0.028;
2817  if(IA>55) kf=0.241;
2818  if(IA>55) Es=0.023;
2819  if(IA>60) kf=0.245;
2820  if(IA>60) Es=0.028;
2821  // changed 5/09
2822  if(IA>55) Es=0.018;
2823 
2824  // Pauli suppression model from Tsai RMP 46,816(74) eq.B54
2825  if ((QV > 2.* kf) || (IA == 1)) {
2826  Pauli_sup2 =1.0;
2827  } else {
2828  Pauli_sup2 = 0.75 * (QV / kf) * (1.0 - (pow((QV / kf),2))/12.);
2829  }
2830  Pauli_sup1 = Pauli_sup2;
2831 
2832  //G4cout << "kf, Es, Paulisup1,2: " << kf << ", " << Es << ", " << Pauli_sup1 << ", " << Pauli_sup2 << G4endl;
2833 
2834  // structure functions with off shell factors
2835  kappa = QV / 2. / amp;
2836  lam = Nu / 2. / amp;
2837  lamp = lam - Es / 2. / amp;
2838  taup = kappa*kappa - lamp*lamp;
2839  squigglef = sqrt(1. + pow((kf/amp),2)) -1.;
2840 
2841  // Very close to treshold, could have a problem
2842  if(1.+lamp <= 0.) return;
2843  if(taup * (1. + taup) <= 0.) return;
2844 
2845  psi = (lam - TAU ) / sqrt(squigglef) / sqrt((1.+lam )* TAU + kappa * sqrt(TAU * (1. + TAU)));
2846  psip = (lamp - taup) / sqrt(squigglef) / sqrt((1.+lamp)*taup + kappa * sqrt(taup * (1. + taup)));
2847  nuL = pow((TAU / kappa / kappa),2);
2848 
2849  // changed definition of nuT from
2850  // nuT = TAU / 2. / kappa**2 + tan(thr/2.)**2
2851  // to this, in order to separate out F1 and F2 (F1 prop. to tan2 term)
2852  nuT = TAU / 2. / kappa / kappa;
2853 
2854  GM2bar = Pauli_sup1 * (Z * GMP*GMP + avgN * GMN*GMN);
2855  GE2bar = Pauli_sup2 * (Z * GEP*GEP + avgN * GEN*GEN);
2856  //G4double W1bar = TAU * GM2bar;
2857  W2bar = (GE2bar + TAU * GM2bar) / (1. + TAU);
2858 
2859  Delta = squigglef * (1. - psi*psi) * (sqrt(TAU * (1.+TAU)) / kappa + squigglef/3. *
2860  (1. - psi*psi) * TAU / kappa / kappa);
2861  GL = kappa*kappa / TAU * (GE2bar + Delta * W2bar) / 2. / kappa / (1. + squigglef *
2862  (1. + psi*psi) / 2.);
2863  GT = (2. * TAU * GM2bar + Delta * W2bar) / 2. / kappa / (1. + squigglef *
2864  (1. + psi*psi) / 2.);
2865 
2866  //G4cout << "nuL, nuT, GL, GT: " << nuL << ", " << nuT << ", " << GL << ", " << GT << G4endl;
2867 
2868  // added to prevent negative xsections:
2869  if (GT < 0) {
2870  GT = 0;
2871  //G4cout << "Reset GT to zero" << G4endl;
2872  }
2873 
2874  // from Maria Barbaro: see Amaro et al., PRC71,015501(2005).
2875  FY = 1.5576 / (1. + 1.7720*1.7720 * pow((psip + 0.3014),2)) / (1. + exp(-2.4291 * psip)) / kf;
2876 
2877 // == Note: Below is for deuteron and haven't == //
2878 // == verified if code conversion is correct == //
2879  // Use PWIA and Paris W.F. for deuteron to get better FY
2880  if (IA == 2) {
2881  // value assuming average p2=0.
2882  pz = (QSQ - 2. * amp * Nu ) / 2. / QV;
2883  izz = int((pz + 1.0) / 0.01) + 1;
2884  if (izz < 1) izz = 1;
2885  if (izz > 200) izz = 200;
2886  izznom = izz;
2887 
2888  // ignoring energy term, estimate change in pz to compensate
2889  //! for avp2 term
2890  dpz = avp2[izznom] / 2. / QV;
2891  izdif = dpz * 150.;
2892  dwmin=1.E6;
2893  izzmin=0;
2894 
2895  G4int izpmax;
2896  if ((izznom + izdif) < 1) izpmax = 1;
2897  if (izpmax > 200) izpmax = 200;
2898  for (izp = izznom; izp <= izpmax; izp++) {
2899  pz = -1. + 0.01 * (izp-0.5);
2900  // *** this version gives worse agreement with laget than
2901  // w2p = (amd + Nu - sqrt(amp**2 + avp2(izp)))**2 -
2902  // > QV**2 + 2. * QV * pz - avp2(izp)
2903  // this version!
2904  w2p = pow((amd + Nu - amp ),2) - QV*QV + 2. * QV * pz - avp2[izp];
2905 
2906  // if passed first minimum, quit looking so don't find second one
2907  if (abs(w2p - amp*amp) > dwmin) {
2908  if (izzmin < 2) izzmin = 2;
2909  if (izzmin > 199) izzmin = 199;
2910  izz = izzmin;
2911  } else if (abs(w2p - amp*amp) < dwmin) {
2912  dwmin = abs(w2p - amp*amp);
2913  izzmin = izp;
2914  }
2915  }
2916 
2917  // search for minimum in 1/10th bins locally
2918  pznom = -1. + 0.01 * (izz-0.5);
2919  dwmin=1.E6;
2920  for (izp = 1; izp <= 19; izp++) {
2921  pz = pznom - 0.01 + 0.001 * izp;
2922  // *** this version gives worse agreement with laget than
2923  // w2p = (amd + Nu - sqrt(amp**2 + avp2(izz)))**2 -
2924  // > QV**2 + 2. * QV * pz - avp2(izz)
2925  // this version!
2926  w2p = pow((amd + Nu - amp ),2) - QV*QV + 2. * QV * pz - avp2[izz];
2927  if (abs(w2p - amp*amp) < dwmin) {
2928  dwmin = abs(w2p - amp*amp);
2929  pzmin = pz;
2930  }
2931  }
2932 
2933  if (dwmin >= 1.e6 || abs(pznom-pzmin) > 0.01) {
2934  // > write(6,'(1x,''error in dwmin,pzmin'',3i4,6f7.3)')
2935  // > izznom,izzmin,izz,QSQ,wsq,w2p,dwmin/1.e6,pzmin,pznom
2936  if (pzmin < pznom) FY = fyd[izz] - (fyd[izz-1] - fyd[izz]) * (pzmin - pznom) / 0.01;
2937  else FY = fyd[izz] + (fyd[izz+1] - fyd[izz]) * (pzmin - pznom) / 0.01;
2938  }
2939  }
2940 // == To here == //
2941 
2942  // final results
2943  F2 = Nu * FY * (nuL * GL + nuT * GT);
2944  F1 = amp * FY * GT / 2.;
2945 
2946  //G4cout << "nu, Fy, nuL, GL, nuT, GT, amp: " << G4endl;
2947  //G4cout << Nu << ", " << FY << ", " << nuL << ", " << GL << ", " << nuT << ", " << GT << ", " << amp << G4endl;
2948 
2949  if (F1 < 0.0) F1 = 0.;
2950  if (Nu > 0. && F1 > 0.) R = (F2 / Nu) / (F1 / amp) * (1. + Nu*Nu / QSQ) - 1.0;
2951  else R = 0.4/QSQ;
2952 
2953 
2954  // apply correction factors
2955  if ( IA > 2 ) {
2956  y = (wsq -amp*amp) / QV;
2957  // F1 = F1 * (1. + pb(8) + pb(9) * y +
2958  // > pb(10)*y**2 + pb(11)*y**3 + pb(12)*y**4 )
2959  // R = R * (1. + pb(13))
2960  // F2 = Nu * F1/amp * (1. + R) / (1. + Nu**2/QSQ)
2961 
2962  // correction to correction Vahe
2963  if (wsq > 0.0) {
2964  F1=F1*(1.0+P[7]+P[8]*y+P[9]*y*y +P[10]*pow(y,3) +P[11]*pow(y,4));
2965  R = R * ( 1.0 + P[12] );
2966  F2 = Nu * F1/amp * (1. + R) / (1. + Nu*Nu/QSQ);
2967  if (F1 < 0.0) F1=0.0;
2968  }
2969  }
2970 
2971  return;
2972 }
2973 
2974 void QweakSimEPEvent::F1F2IN09(G4int Z, G4int IA, G4double qsq,
2975  G4double wsq, G4double &F1, G4double &F2)
2976 {
2977 /*--------------------------------------------------------------------
2978  Fit to inelastic cross sections for A(e,e')X
2979  valid for all W<3 GeV and all Q2<10 GeV2
2980 
2981  Inputs: Z, A (real*8) are Z and A of nucleus
2982  (use Z=0., A=1. to get free neutron)
2983  Qsq (real*8) is 4-vector momentum transfer squared (positive in
2984  chosen metric)
2985  Wsq (real*8) is invarinat mass squared of final state calculated
2986  assuming electron scattered from a free proton
2987 
2988  outputs: F1, F2 (real*8) are structure functions per nucleus
2989  Version of 10/20/2006 P. Bosted
2990 --------------------------------------------------------------------*/
2991 
2992  G4double W, x;
2993  G4double nu,siglp,sigtp,F1pp,F1dp;
2994  G4double W1,W2,sigt,Rc,sigl,F1d, F1p,qv;
2995  G4double DW2DPF,Wsqp,pf,kf,Es,dw2des,Fyuse;
2996  G4double x4, emcfac;
2997  G4int ism;
2998  G4double PM = 0.93828;
2999 
3000  // This is for exp(-xx**2/2.), from teste.f
3001  G4double XXp[99] = {
3002  -3.000,-2.939,-2.878,-2.816,-2.755,-2.694,-2.633,-2.571,-2.510,
3003  -2.449,-2.388,-2.327,-2.265,-2.204,-2.143,-2.082,-2.020,-1.959,
3004  -1.898,-1.837,-1.776,-1.714,-1.653,-1.592,-1.531,-1.469,-1.408,
3005  -1.347,-1.286,-1.224,-1.163,-1.102,-1.041,-0.980,-0.918,-0.857,
3006  -0.796,-0.735,-0.673,-0.612,-0.551,-0.490,-0.429,-0.367,-0.306,
3007  -0.245,-0.184,-0.122,-0.061, 0.000, 0.061, 0.122, 0.184, 0.245,
3008  0.306, 0.367, 0.429, 0.490, 0.551, 0.612, 0.673, 0.735, 0.796,
3009  0.857, 0.918, 0.980, 1.041, 1.102, 1.163, 1.224, 1.286, 1.347,
3010  1.408, 1.469, 1.531, 1.592, 1.653, 1.714, 1.776, 1.837, 1.898,
3011  1.959, 2.020, 2.082, 2.143, 2.204, 2.265, 2.327, 2.388, 2.449,
3012  2.510, 2.571, 2.633, 2.694, 2.755, 2.816, 2.878, 2.939, 3.000};
3013  // these are 100x bigger for convenience
3014  G4double fyp[99] = {
3015  0.0272,0.0326,0.0390,0.0464,0.0551,0.0651,0.0766,0.0898,0.1049,
3016  0.1221,0.1416,0.1636,0.1883,0.2159,0.2466,0.2807,0.3182,0.3595,
3017  0.4045,0.4535,0.5066,0.5637,0.6249,0.6901,0.7593,0.8324,0.9090,
3018  0.9890,1.0720,1.1577,1.2454,1.3349,1.4254,1.5163,1.6070,1.6968,
3019  1.7849,1.8705,1.9529,2.0313,2.1049,2.1731,2.2350,2.2901,2.3379,
3020  2.3776,2.4090,2.4317,2.4454,2.4500,2.4454,2.4317,2.4090,2.3776,
3021  2.3379,2.2901,2.2350,2.1731,2.1049,2.0313,1.9529,1.8705,1.7849,
3022  1.6968,1.6070,1.5163,1.4254,1.3349,1.2454,1.1577,1.0720,0.9890,
3023  0.9090,0.8324,0.7593,0.6901,0.6249,0.5637,0.5066,0.4535,0.4045,
3024  0.3595,0.3182,0.2807,0.2466,0.2159,0.1883,0.1636,0.1416,0.1221,
3025  0.1049,0.0898,0.0766,0.0651,0.0551,0.0464,0.0390,0.0326,0.0272};
3026 
3027  // deuteron fit parameters
3028  G4double xvald0[50] = {
3029  0.1964E+01, 0.1086E+01, 0.5313E-02, 0.1265E+01, 0.8000E+01,
3030  0.2979E+00, 0.1354E+00, 0.2200E+00, 0.8296E-01, 0.9578E-01,
3031  0.1094E+00, 0.3794E+00, 0.8122E+01, 0.5189E+01, 0.3290E+01,
3032  0.1870E+01, 0.6110E+01,-0.3464E+02, 0.9000E+03, 0.1717E+01,
3033  0.4335E-01, 0.1915E+03, 0.2232E+00, 0.2119E+01, 0.2088E+01,
3034  -0.3029E+00, 0.2012E+00, 0.1104E-02, 0.2276E-01,-0.4562E+00,
3035  0.2397E+00, 0.1204E+01, 0.2321E-01, 0.5419E+03, 0.2247E+00,
3036  0.2168E+01, 0.2266E+03, 0.7649E-01, 0.1457E+01, 0.1318E+00,
3037  -0.7534E+02, 0.1776E+00, 0.1636E+01, 0.1350E+00,-0.5596E-02,
3038  0.5883E-02, 0.1934E+01, 0.3800E+00, 0.3319E+01, 0.1446E+00};
3039 
3040  G4double P[24] = {
3041  5.1377e-03, 9.8071e-01, 4.6379e-02, 1.6433e+00,
3042  6.9826e+00, -2.2655e-01, 1.1095e-01, 2.7945e-02,
3043  4.0643e-01, 1.6076e+00, -7.5460e+00, 4.4418e+00,
3044  -3.7464e-01, 1.0414e-01, -2.6852e-01, 9.6653e-01,
3045  -1.9055e+00, 9.8965e-01, 2.0613e+02, -4.5536e-02,
3046  2.4902e-01, -1.3728e-01, 2.9201e+01, 4.9280e-03};
3047 
3048  G4double F1M;
3049  nu = (wsq - PM*PM + qsq) / 2. / PM;
3050  qv = sqrt(nu*nu + qsq);
3051 
3052  if (wsq <= 0.0) {
3053  W = 0.0;
3054  } else {
3055  W = sqrt(wsq);
3056  x = qsq / (2.0 * PM * nu);
3057  }
3058  if (wsq <= 0.0) x = 0.0;
3059 
3060  if (IA <= 2) return;
3061 
3062  // For nuclei
3063  if (IA > 2) {
3064  sigt = 0.;
3065  sigl = 0.;
3066  F1d = 0.;
3067  F1p = 0.;
3068  }
3069 
3070  // Modifed to use Superscaling from Sick, Donnelly, Maieron,
3071  // nucl-th/0109032
3072  if(IA==2) kf=0.085;
3073  if(IA==2) Es=0.0022;
3074  // changed 4/09
3075  if(IA==3) kf=0.115;
3076  if(IA==3) Es=0.001 ;
3077  // changed 4/09
3078  if(IA>3) kf=0.19;
3079  if(IA>3) Es=0.017;
3080  if(IA>7) kf=0.228;
3081  if(IA>7) Es=0.020;
3082  // changed 5/09
3083  if(IA>7) Es=0.0165;
3084  if(IA>16) kf=0.230;
3085  if(IA>16) Es=0.025;
3086  if(IA>25) kf=0.236;
3087  if(IA>25) Es=0.018;
3088  if(IA>38) kf=0.241;
3089  if(IA>38) Es=0.028;
3090  if(IA>55) kf=0.241;
3091  if(IA>55) Es=0.023;
3092  if(IA>60) kf=0.245;
3093  if(IA>60) Es=0.028;
3094  // changed 5/09
3095  if(IA>55) Es=0.018;
3096 
3097  // adjust pf to give right width based on kf
3098  pf = 0.5 * kf;
3099  // assume this is 2 * pf * qv
3100  DW2DPF = 2. * qv;
3101  dw2des = 2. * (nu + PM) ;
3102  // switched to using 99 bins!
3103  for (ism = 0; ism < 99; ism++) {
3104  Fyuse = fyp[ism]/100.;
3105  Wsqp = wsq + XXp[ism] * pf * DW2DPF - Es * dw2des;
3106  if (Wsqp > 1.159) {
3107  christy507(Wsqp,qsq,F1pp,Rc,sigtp,siglp);
3108  resmodd(Wsqp,qsq,xvald0,F1dp);
3109  F1d = F1d + F1dp * Fyuse;
3110  F1p = F1p + F1pp * Fyuse;
3111  sigt = sigt + sigtp * Fyuse;
3112  sigl = sigl + siglp * Fyuse;
3113  }
3114  }
3115 
3116  Rc = 0.;
3117  if(sigt > 0.) Rc = sigl / sigt;
3118  W1 = (2. * Z * F1d + (IA - 2. * Z) * (2. * F1d - F1p)) / PM;
3119  W1= W1*(1.0+P[13]*x+P[14]*pow(x,2)+P[15]*pow(x,3)+P[16]*pow(x,4)+P[17]*pow(x,5));
3120 
3121  if(W > 0.0) W1=W1*pow((1.0+(P[20]*W+P[21]*pow(W,2))/(1.0+P[22]*qsq)),2);
3122  F1M = MEC2009(IA,qsq,wsq);
3123  W1 = W1 + F1M;
3124  if(wsq > 0.0 ) Rc = Rc * ( 1.0 + P[6] + P[23]*IA );
3125  W2 = W1 * (1. + Rc) / (1. + nu*nu / qsq);
3126 
3127  x4 = qsq / 2. / PM / nu;
3128  emcfac = fitemc(x4, IA);
3129  F1 = PM * W1 * emcfac;
3130  F2 = nu * W2 * emcfac;
3131 
3132  return;
3133 
3134 }
3135 
3136 //-----------------------------------------------------------------//
3137 void QweakSimEPEvent::christy507(G4double w2,G4double q2,G4double &F1,
3138  G4double &R, G4double &sigT, G4double &sigL)
3139 {
3140 // M.E. Christy and P.E. Bosted, ``Empirical Fit to Precision
3141 // Inclusive Electron-Proton Cross Sections in the Resonance Region'',
3142 // (arXiv:0712.3731). To be submitted to Phys. Rev. C.
3143 
3144  G4double xval1[50],xvalL[50];
3145 
3146  G4double mp = .9382727;
3147  G4double mp2 = mp*mp;
3148  G4double alpha = 1./137.036;
3149 
3150  G4double xval[100] = {
3151  0.12298E+01,0.15304E+01,0.15057E+01,0.16980E+01,0.16650E+01,
3152  0.14333E+01,0.13573E+00,0.22000E+00,0.82956E-01,0.95782E-01,
3153  0.10936E+00,0.37944E+00,0.77805E+01,0.42291E+01,0.12598E+01,
3154  0.21242E+01,0.63351E+01,0.68232E+04,0.33521E+05,0.25686E+01,
3155  0.60347E+00,0.21240E+02,0.55746E-01,0.24886E+01,0.23305E+01,
3156  -.28789E+00,0.18607E+00,0.63534E-01,0.19790E+01,-.56175E+00,
3157  0.38964E+00,0.54883E+00,0.22506E-01,0.46213E+03,0.19221E+00,
3158  0.19141E+01,0.24606E+03,0.67469E-01,0.13501E+01,0.12054E+00,
3159  -.89360E+02,0.20977E+00,0.15715E+01,0.90736E-01,-.38495E-02,
3160  0.10362E-01,0.19341E+01,0.38000E+00,0.34187E+01,0.14462E+00,
3161  0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,
3162  0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,
3163  0.00000E+00,0.00000E+00,0.29414E+02,0.19910E+02,0.22587E+00,
3164  0.00000E+00,0.00000E+00,0.38565E+04,0.65717E+00,0.00000E+00,
3165  0.15792E+03,0.97046E+02,0.31042E+00,0.00000E+00,0.42160E+01,
3166  0.38200E-01,0.12182E+01,0.00000E+00,0.13764E+02,0.31393E+00,
3167  0.29997E+01,0.00000E+00,0.55124E+01,0.53743E-01,0.13091E+01,
3168  0.00000E+00,0.86746E+02,0.40864E-04,0.40294E+01,0.31285E+01,
3169  0.33403E+00,0.49623E+01,0.00000E+00,0.00000E+00,0.11000E+02,
3170  0.18951E+01,0.51376E+00,0.00000E+00,0.42802E+01,0.00000E+00};
3171 
3172  for (G4int i = 0; i<50; i++) {
3173  xval1[i] = xval[i];
3174  xvalL[i] = xval[50+i] ;
3175  if(i < 12) xvalL[i] = xval1[i];
3176  }
3177  xvalL[42] = xval1[46];
3178  xvalL[43] = xval1[47];
3179  xvalL[49] = xval1[49];
3180 
3181  //G4double xb = q2/(w2+q2-mp2);
3182  sigT = resmod507_v2(1,w2,q2,xval1);
3183  sigL = resmod507_v2(2,w2,q2,xvalL);
3184 
3185  F1 = sigT*(w2-mp2)/8./pi/pi/alpha/0.3894e3;
3186  //G4double FL = sigL*2.*xb*(w2-mp2)/8./pi/pi/alpha/0.3894e3;
3187  R = sigL/sigT;
3188 
3189  return;
3190 }
3191 
3192 // -------------------------------------------------------------------------//
3193 void QweakSimEPEvent::resmodd(G4double w2, G4double q2,
3194  G4double xval[50], G4double &sig)
3195 {
3196  //! returns F1 for average of free proton and neutron
3197  //! for given W2, Q2
3198 
3199  G4double W,mp,mp2,mass[7],width[7];
3200  G4double height[7],rescoef[6][4];
3201  G4double nr_coef[3][4],wdif,sig_nr;
3202  G4double sigr[7];
3203  G4double mpi,meta,intwidth[7],k,kcm,kr[7],kcmr[7],ppicm,ppi2cm;
3204  G4double petacm,ppicmr[7],ppi2cmr[7],petacmr[7],epicmr[7],epi2cmr[7];
3205  G4double eetacmr[7],epicm,epi2cm,eetacm,br[7][2],ang[7],sigrsv[7];
3206  G4double pgam[7],pwid[7][2],x0[7],dip;
3207  G4double sig_res,xpr,alpha,F1;
3208  G4int i,j,num,iw;
3209 
3210  G4double xval0[12] = {
3211  0.12298E+01,0.15304E+01,0.15057E+01,0.16980E+01,0.16650E+01,
3212  0.14333E+01,0.13573E+00,0.22000E+00,0.82956E-01,0.95782E-01,
3213  0.10936E+00,0.37944E+00
3214  };
3215 
3216  G4double w2sv;
3217 
3218  mp = 0.9382727;
3219  mpi = 0.135;
3220  meta = 0.547;
3221  mp2 = mp*mp;
3222  alpha = 1./137.036;
3223 
3224  sig = 0.;
3225 
3226  if(w2 < 1.07327*1.07327 || w2 > 25 || q2 < 0.0 || q2 > 11.0) {
3227  G4cerr << "ERROR: " << __FILE__ << " line " << __LINE__ << G4endl;
3228  G4cerr << " W2/Q2 check failed: should be 1.07327^2<W2<25 and 0 < Q2 < 11" << G4endl;
3229  G4cerr << " but are w2 q2: " << w2 << " " << q2 << G4endl;
3230  exit(1);
3231  }
3232 
3233  // ! branching ratios
3234  br[0][0] = 1.0;
3235  br[1][0] = 0.5;
3236  br[2][0] = 0.65;
3237  br[3][0] = 0.65;
3238  br[4][0] = 0.4;
3239  br[5][0] = 0.65;
3240  br[6][0] = 0.6;
3241 
3242  // ! angular momenta
3243  ang[0] = 1.; // !!! P33(1232)
3244  ang[1] = 0.; // !!! S11(1535)
3245  ang[2] = 2.; // !!! D13(1520)
3246  ang[3] = 3.; // !!! F15(1680)
3247  ang[4] = 0.; // !!! S15(1650)
3248  ang[5] = 1.; // !!! P11(1440) roper
3249  ang[6] = 3.; // !!! ? 4th resonance region
3250 
3251  // ! x0 parameter
3252  for (i = 0; i < 7; i++) x0[i] = 0.215;
3253  x0[0] = 0.1446;
3254 
3255  // ! out branching ratio
3256  for (i = 0; i < 7; i++) br[i][1] = 1.-br[i][0];
3257 
3258  // ! remember w2
3259  w2sv = w2;
3260 
3261  // ! uses xvals of 1-12, 47, and 48
3262  // ! move masses, wdiths into local variables
3263  // ! pyb changed to be fixed
3264  num = 0;
3265  for (i = 0; i < 6; i++) {
3266  num = num + 1;
3267  mass[i] = xval0[i];
3268  }
3269  for (i = 0; i < 6; i++) {
3270  num = num + 1;
3271  intwidth[i] = xval0[num-1];
3272  }
3273 
3274  //! changed to allow delta width, mass to vary
3275  //! taken out again since xval[1] used in MEC
3276  mass[6] = 1.9341;
3277  intwidth[6] = 0.380;
3278 
3279  iw = G4int(1000.*sqrt(w2));
3280  W = 0.001 * (iw+0.5);
3281  w2 = W*W;
3282  wdif = W - (mp + mpi);
3283 
3284  // ! Calculate kinematics needed for threshold Relativistic B-W
3285  k = (w2 - mp2) / 2. / mp;
3286  kcm = (w2 - mp2) / 2. / W;
3287  epicm = (w2 + mpi*mpi -mp2 ) / 2. / W;
3288  ppicm = pow(std::max(0.0,(epicm*epicm - mpi*mpi)),0.5);
3289  epi2cm = (w2 + pow((2.*mpi),2) -mp2 ) / 2. / W;
3290  ppi2cm = pow(std::max(0.0,(epi2cm*epi2cm - pow((2.*mpi),2))),0.5);
3291  eetacm = (w2 + meta*meta -mp2 ) / 2. / W;
3292  petacm = pow(std::max(0.0,(eetacm*eetacm - meta*meta)),0.5);
3293  for (i = 0; i < 7; i++) {
3294  kr[i] = (mass[i]*mass[i]-mp2)/2./mp;
3295  kcmr[i] = (mass[i]*mass[i]-mp2)/2./mass[i];
3296  epicmr[i] = (mass[i]*mass[i] + mpi*mpi -mp2 )/2./mass[i];
3297  ppicmr[i] = pow(std::max(0.0,(epicmr[i]*epicmr[i] - mpi*mpi)),0.5);
3298  epi2cmr[i] = (mass[i]*mass[i] + pow(2.*mpi,2) -mp2 )/2./mass[i];
3299  ppi2cmr[i] = pow(std::max(0.0,(epi2cmr[i]*epi2cmr[i] - pow((2.*mpi),2))),0.5);
3300  eetacmr[i] = (mass[i]*mass[i] + meta*meta -mp2 )/2./mass[i];
3301  petacmr[i] = pow(std::max(0.0,(eetacmr[i]*eetacmr[i] - meta*meta)),0.5);
3302 
3303  // ! Calculate partial widths
3304  pwid[i][0] = intwidth[i]*pow((ppicm/ppicmr[i]),(2.*ang[i]+1.))*pow(((ppicmr[i]*ppicmr[i]+x0[i]*x0[i])/(ppicm*ppicm+x0[i]*x0[i])),ang[i]);
3305 
3306  if(i != 1)
3307  pwid[i][1] = intwidth[i]*pow((ppi2cm/ppi2cmr[i]),(2.*ang[i]+4.))*
3308  pow(((ppi2cmr[i]*ppi2cmr[i]+x0[i]*x0[i])/(ppi2cm*ppi2cm+x0[i]*x0[i])),(ang[i]+2.))* W / mass[i];
3309  else
3310  pwid[i][1] = intwidth[1]*pow((petacm/petacmr[i]),(2.*ang[i]+1.))*
3311  pow(((petacmr[i]*petacmr[i]+x0[i]*x0[i])/(petacm*petacm+x0[i]*x0[i])),ang[i]);
3312 
3313  pgam[i] = pow((kcm/kcmr[i]),2)*(kcmr[i]*kcmr[i]+x0[i]*x0[i])/(kcm*kcm+x0[i]*x0[i]);
3314  pgam[i] = intwidth[i]*pgam[i];
3315  width[i] = br[i][0]*pwid[i][0]+br[i][1]*pwid[i][1];
3316  sigr[i] = width[i] * pgam[i] / (pow(w2 - mass[i]*mass[i],2) + pow(mass[i]*width[i],2)) * kr[i] / k * kcmr[i] / kcm / intwidth[i];
3317  }
3318 
3319  w2 = w2sv;
3320 
3321  // ! get parameters into local variables
3322  num = 12;
3323  //! resonance height coefficients. xvals of 13-36
3324  for (i = 0; i < 6; i++) {
3325  for (j = 0; j < 4; j++) {
3326  num = num + 1;
3327  rescoef[i][j]=xval[num-1];
3328  }
3329  }
3330  // ! Non-Res coefficients xvals of 37-44
3331  for (i = 0; i < 2; i++) {
3332  for (j = 0; j < 4; j++) {
3333  num = num + 1;
3334  nr_coef[i][j]=xval[num-1];
3335  }
3336  }
3337 
3338  // ! Begin resonance Q^2 dependence calculations CCC
3339  // ! uses xvals 49
3340  for (i = 0; i < 6; i++) {
3341  height[i] = rescoef[i][0]*(1.+rescoef[i][1] * q2 / (1. + rescoef[i][2] * q2))/pow((1. + q2/0.91),rescoef[i][3]);
3342  }
3343  dip = 1./pow((1. + q2 / 0.91),2);
3344  height[6] = xval[48]*dip;
3345  sig_res = 0.;
3346  for (i = 0; i < 7; i++) {
3347  sigrsv[i] = height[i]*height[i] * sigr[i];
3348  sig_res = sig_res + sigrsv[i];
3349  }
3350 
3351  sig_res = sig_res * sqrt(w2);
3352 
3353  //! Begin non-resonant part uses xvals 45, 46, 50
3354  //! Depends on both W2 and Q2 so can't easily precalculate
3355  sig_nr = 0.;
3356  xpr = 1.+(w2-pow(mp+mpi,2))/(q2+0.05);
3357  xpr = 1./xpr;
3358  W = pow(w2,0.5);
3359  wdif = W - (mp + mpi);
3360  for (i = 0; i < 2; i++) {
3361  sig_nr = sig_nr +(nr_coef[i][0]*pow(wdif,((2*(i+1)+1)/2.)))/pow(q2+nr_coef[i][1],(nr_coef[i][2]+nr_coef[i][3]*q2+xval[44+i]*q2*q2));
3362  }
3363  sig_nr = sig_nr * xpr;
3364  sig = sig_res + sig_nr;
3365 
3366  F1 = sig * (w2-mp2)/8./pi/pi/alpha/0.3894e3;
3367  sig = F1;
3368 }
3369 
3370 
3371 //------------------------------------------------------------------//
3372 // Used for Aluminum inelastic...
3373 G4double QweakSimEPEvent::resmod507_v2(G4double sf,G4double w2,
3374  G4double q2,G4double xval[50])
3375 {
3376 
3377  G4double W,mp,mp2,xb;
3378  G4double mass[7] = {0,0,0,0,0,0,0};
3379  G4double width[7] = {0,0,0,0,0,0,0};
3380  G4double height[7] = {0,0,0,0,0,0,0};
3381  G4double intwidth[7] = {0,0,0,0,0,0,0};
3382  G4double rescoef[6][4];
3383  G4double nr_coef[3][4],sigr[7],wdif[2],sig_nr,h_nr[3];
3384  G4double mpi,meta,k,kcm,kr[7],kcmr[7],ppicm,ppi2cm;
3385  G4double petacm,ppicmr[7],ppi2cmr[7],petacmr[7],epicmr[7],epi2cmr[7];
3386  G4double eetacmr[7],epicm,epi2cm,eetacm,br[7][3],ang[7];
3387  G4double pgam[7],pwid[7][3],x0[7],q20;
3388  G4double sig_res,t,xpr[2],m0,sig;
3389  G4int i,j,num;
3390 
3391  mp = 0.9382727;
3392  mpi = 0.135;
3393  meta = 0.547;
3394  mp2 = mp*mp;
3395  W = sqrt(w2);
3396  wdif[0] = W - (mp + mpi);
3397  wdif[1] = W - (mp + 2.*mpi);
3398 
3399  m0 = 0.125;
3400  if(sf == 2) m0 = xval[48];
3401  if(sf == 1) {
3402  q20 = 0.05;
3403  } else {
3404  q20 = 0.125;
3405  }
3406 
3407  // single pion branching ratios CCCC
3408 
3409  br[0][0] = 1.0; //!!! P33(1232)
3410  br[1][0] = 0.45; //!!! S11(1535)
3411  br[2][0] = 0.65; //!!! D13(1520)
3412  br[3][0] = 0.65; //!!! F15(1680)
3413  br[4][0] = 0.4; //!!! S11(1650)
3414  br[5][0] = 0.65; //!!! P11(1440) roper
3415  br[6][0] = 0.50 ; //!!! F37(1950)
3416 
3417  br[0][2] = 0.0; //!!! P33(1232)
3418  br[1][2] = 0.45; //!!! S11(1535)
3419  br[2][2] = 0.0; //!!! D13(1520)
3420  br[3][2] = 0.0; //!!! F15(1680)
3421  br[4][2] = 0.1; //!!! S11(1650)
3422  br[5][2] = 0.0; //!!! P11(1440) roper
3423  br[6][2] = 0.0; //!!! F37(1950)
3424 
3425  // 2-pion branching ratios CCCC
3426 
3427  for (i = 0; i < 7; i++) {
3428  br[i][1] = 1.-br[i][0]-br[i][2];
3429  }
3430 
3431  // Meson angular momentum CCCC
3432  ang[0] = 1.; //!!! P33(1232)
3433  ang[1] = 0.; //!!! S11(1535)
3434  ang[2] = 2.; //!!! D13(1520)
3435  ang[3] = 3.; //!!! F15(1680)
3436  ang[4] = 0.; //!!! S15(1650)
3437  ang[5] = 1.; //!!! P11(1440) roper
3438  ang[6] = 3.; //!!! F37(1950)
3439 
3440  for (i = 0; i < 7; i++) { //!!! resonance damping parameter !!!
3441  x0[i] = 0.215;
3442  }
3443  x0[0] = 0.15;
3444  x0[0] = xval[49];
3445 
3446  xb = q2/(q2+w2-mp2);
3447  xpr[0] = 1.+(w2-pow((mp+mpi),2))/(q2+q20);
3448  xpr[0] = 1./xpr[0];
3449  xpr[1] = 1.+(w2-pow((mp+mpi+mpi),2))/(q2+q20);
3450  xpr[1] = 1./xpr[1];
3451 
3452  t = log(log((q2+m0)/0.330/0.330)/log(m0/0.330/0.330));
3453 
3454  // Calculate kinematics needed for threshold Relativistic B-W
3455 
3456  k = (w2 - mp2)/2./mp;
3457  kcm = (w2-mp2)/2./W;
3458 
3459  epicm = (w2 + mpi*mpi -mp2 )/2./W;
3460  ppicm = pow(std::max(0.0,(epicm*epicm - mpi*mpi)),0.5);
3461  epi2cm = (w2 + pow((2.*mpi),2) -mp2 )/2./W;
3462  ppi2cm = pow(std::max(0.0,(epi2cm*epi2cm - pow((2.*mpi),2))),0.5);
3463  eetacm = (w2 + meta*meta -mp2 )/2./W;
3464  petacm = pow(std::max(0.0,(eetacm*eetacm - meta*meta)),0.5);
3465 
3466  num = 0;
3467 
3468  for (i = 0; i < 6; i++) { // !!! Read in resonance masses !!!
3469  num = num + 1;
3470  mass[i] = xval[i];
3471  }
3472  for (i = 0; i < 6; i++) { // !!! Read in resonance widths !!!
3473  num = num + 1;
3474  intwidth[i] = xval[num-1];
3475  width[i] = intwidth[i];
3476  }
3477  if (sf == 2) { // !!! Put in 4th resonance region !!!
3478  mass[6] = xval[42];
3479  intwidth[6] = xval[43];
3480  width[6] = intwidth[6];
3481  } else {
3482  mass[6] = xval[46];
3483  intwidth[6] = xval[47];
3484  width[6] = intwidth[6];
3485  }
3486 
3487  for (i = 0; i < 7; i++) {
3488  kr[i] = (mass[i]*mass[i]-mp2)/2./mp;
3489  kcmr[i] = (mass[i]*mass[i]-mp2)/2./mass[i];
3490  epicmr[i] = (mass[i]*mass[i] + mpi*mpi -mp2 )/2./mass[i];
3491  ppicmr[i] = pow(std::max(0.0,(epicmr[i]*epicmr[i] - mpi*mpi)),0.5);
3492  epi2cmr[i] = (mass[i]*mass[i] + pow((2.*mpi),2) -mp2 )/2./mass[i];
3493  ppi2cmr[i] = pow(std::max(0.0,(epi2cmr[i]*epi2cmr[i] - pow((2.*mpi),2))),0.5);
3494  eetacmr[i] = (mass[i]*mass[i] + meta*meta -mp2 )/2./mass[i];
3495  petacmr[i] = pow(std::max(0.0,(eetacmr[i]*eetacmr[i] - meta*meta)),0.5);
3496 
3497  // CCC Calculate partial widths CCC
3498 
3499  pwid[i][0] = intwidth[i]*pow((ppicm/ppicmr[i]),(2.*ang[i]+1.))*pow(((ppicmr[i]*ppicmr[i]+x0[i]*x0[i])/(ppicm*ppicm+x0[i]*x0[i])),ang[i]); // !!! 1-pion decay mode
3500 
3501  pwid[i][1] = intwidth[i]*pow((ppi2cm/ppi2cmr[i]),(2.*ang[i]+4.))*pow(((ppi2cmr[i]*ppi2cmr[i]+x0[i]*x0[i])/(ppi2cm*ppi2cm+x0[i]*x0[i])),(ang[i]+2)); // !!! 2-pion decay mode
3502 
3503  pwid[i][1] = W/mass[i]*pwid[i][1];
3504  pwid[i][2] = 0.; // !!! eta decay mode
3505 
3506  if(i == 1 || i == 4) {
3507  pwid[i][2] = intwidth[i]*pow((petacm/petacmr[i]),(2.*ang[i]+1.))*pow(((petacmr[i]*petacmr[i]+x0[i]*x0[i])/(petacm*petacm+x0[i]*x0[i])),ang[i]); // !!! eta decay only for S11's
3508  }
3509 
3510  pgam[i] = pow((kcm/kcmr[i]),2)*(pow(kcmr[i],2)+x0[i]*x0[i])/(kcm*kcm+x0[i]*x0[i]);
3511  pgam[i] = intwidth[i]*pgam[i];
3512  width[i] = br[i][0]*pwid[i][0]+br[i][1]*pwid[i][1]+br[i][2]*pwid[i][2];
3513  }
3514 
3515  //CCC End resonance kinematics and Widths calculations CCC
3516 
3517  // CCC Begin resonance Q^2 dependence calculations CCC
3518 
3519  for (i = 0; i < 6; i++) {
3520  for (j = 0; j < 4; j++) {
3521  num = num + 1;
3522  rescoef[i][j]=xval[num-1];
3523  }
3524 
3525  if(sf == 1) {
3526  height[i] = rescoef[i][0]*(1.+rescoef[i][1]*q2/(1.+rescoef[i][2]*q2))/pow((1.+q2/0.91),rescoef[i][3]);
3527  } else {
3528  height[i] = rescoef[i][0]*q2/(1.+rescoef[i][1]*q2)*exp(-1.*rescoef[i][2]*q2);
3529  }
3530 
3531  height[i] = height[i]*height[i];
3532  }
3533 
3534  if(sf == 2) { // !!! 4th resonance region !!!
3535  height[6] = xval[44]*q2/(1.+xval[45]*q2)*exp(-1.*xval[46]*q2);
3536  } else {
3537  height[6] = xval[48]/(1.+q2/0.91);
3538  }
3539 
3540  height[6] = height[6]*height[6];
3541 
3542  // CCC End resonance Q^2 dependence calculations CCC
3543 
3544  for (i = 0; i < 3; i++) { // !!! Non-Res coefficients !!!
3545  for (j = 0; j < 4; j++) {
3546  num = num + 1;
3547  nr_coef[i][j]=xval[num-1];
3548  }
3549  }
3550 
3551  // CCC Calculate Breit-Wigners for all resonances CCC
3552 
3553  sig_res = 0.0;
3554  for (i = 0; i < 7; i++) {
3555  sigr[i] = width[i]*pgam[i]/(pow(w2 - mass[i]*mass[i],2.) + pow(mass[i]*width[i],2.));
3556  sigr[i] = height[i]*kr[i]/k*kcmr[i]/kcm*sigr[i]/intwidth[i];
3557  sig_res = sig_res + sigr[i];
3558  }
3559  sig_res = sig_res*W;
3560  // CCC Finish resonances / start non-res background calculation CCC
3561  sig_nr = 0.;
3562  if(sf == 1) {
3563  for (i = 0; i < 2; i++) {
3564  h_nr[i] = nr_coef[i][0]/pow(q2+nr_coef[i][1],nr_coef[i][2] + nr_coef[i][3]*q2+xval[44+i]*q2*q2);
3565  sig_nr = sig_nr +h_nr[i]*(pow(wdif[0],(2*(i+1)+1)/2.));
3566  }
3567 
3568  sig_nr = sig_nr*xpr[0];
3569 
3570  } else if (sf == 2) {
3571  for (i = 0; i < 1; i++) {
3572  sig_nr = sig_nr + nr_coef[i][0]*pow((1.-xpr[i]),(nr_coef[i][2]+nr_coef[i][1]*t))/(1.-xb)/(q2+q20);
3573  sig_nr = sig_nr*pow(q2/(q2+q20),nr_coef[i][3])*pow(xpr[i],(xval[40]+xval[41]*t));
3574  }
3575  }
3576  sig = sig_res + sig_nr;
3577 
3578  return sig;
3579 }
3580 
3581 G4double QweakSimEPEvent::MEC2009(G4int a, G4double q2,G4double w2)
3582 {
3583  //! fit to low q2 dip region: purefly empirical
3584  //! assume contribution is pure transverse
3585  G4double f1 = 0.0;
3586  G4double am = 0.9383;
3587  G4double w,nu;
3588 
3589  G4double P[24] = {
3590  5.1377e-03, 9.8071e-01, 4.6379e-02, 1.6433e+00,
3591  6.9826e+00, -2.2655e-01, 1.1095e-01, 2.7945e-02,
3592  4.0643e-01, 1.6076e+00, -7.5460e+00, 4.4418e+00,
3593  -3.7464e-01, 1.0414e-01, -2.6852e-01, 9.6653e-01,
3594  -1.9055e+00, 9.8965e-01, 2.0613e+02, -4.5536e-02,
3595  2.4902e-01, -1.3728e-01, 2.9201e+01, 4.9280e-03};
3596 
3597  G4double p18, x, f1corr;
3598 
3599  if (w2 <= 0.0) return 0;
3600  w = pow(w2,0.5);
3601  nu = (w2 - am*am + q2) / 2. / am;
3602  x = q2 / (2.0 * am * nu );
3603 
3604  if (a < 2.5) return 0;
3605 
3606  p18 = P[18];
3607  // ! special case for 3He
3608  if (a > 2.5 && a < 3.5) p18 = 70;
3609  // ! special case for 4He
3610  if (a > 3.5 && a < 4.5) p18 = 170;
3611  // ! new values for C, Al, Cu
3612  if(a > 4.5) p18 = 215;
3613  if(a > 20.) p18 = 235;
3614  if(a > 50.) p18 = 230;
3615 
3616 
3617  f1corr = P[0]*exp(-pow(w-P[1],2)/P[2])/(pow(1.0 + std::max(0.3,q2)/P[3],P[4]))*pow(nu,P[5])*(1.0 + p18 * pow(a,(1.0 + P[19] * x)));
3618 
3619  f1 = f1corr;
3620 
3621  if (f1 <= 1.0E-9 ) f1 = 0.0;
3622 
3623  return f1;
3624 }
3625 
3626 G4double QweakSimEPEvent::fitemc(G4double X, G4int A)
3627 {
3628  /*!---------------------------------------------------------------------
3629  ! Fit to EMC effect. Steve Rock 8/3/94
3630  ! Funciton returns value of sigma(A)/sigma(d)
3631  ! with no isoscalerity correction
3632  ! A= atomic number
3633  ! x = Bjorken x.
3634  !
3635  ! Fit of sigma(A)/sigma(d) to form C*A**alpha where A is atomic number
3636  ! First data at each x was fit to form C*A**alpha. The point A=2(d)
3637  ! was includded with a value of 1 and an error of about 2%.
3638  ! For x>=.125 Javier Gomez fit of 7/93 to E139 data was used.
3639  ! For .09 >=x>=.0085 NMC data from Amaudruz et al Z. Phys C. 51,387(91)
3640  ! Steve did the fit for alpha and C to the He/d. C/d and Ca/d NMC data.
3641  ! Alpha(x) was fit to a 9 term polynomial a0 +a1*x +a2*x**2 + a3*x**3 ..
3642  ! C(x) was fit to a 3 term polynomial in natural logs as
3643  ! Ln(C) = c0 + c1*Ln(x) + c2*[Ln(x)]**2.
3644 
3645  ! 6/2/98 ***** Bug (which set x= .00885 if x was out of range) fixed
3646  ! also gave value at x=.0085 if x>.88
3647  ! 11/05 PYB PEB modified to use value at x=0.7 if x>0.7, because
3648  ! beyond that assume Fermi motion is giving the rise, and we
3649  ! already are taking that into account with the y-smearing of
3650  ! the inelastic
3651  !-----------------------------------------------------------------------
3652  */
3653  G4double ALPHA,C,LN_C,X_U;
3654  G4double fitval;
3655 
3656  // !Chisq= 19. for 30 points
3657  // !Term Coeficient
3658  G4double ALPHA_COEF[9] = {
3659  -6.98871401E-02,
3660  2.18888887E+00,
3661  -2.46673765E+01,
3662  1.45290967E+02,
3663  -4.97236711E+02,
3664  1.01312929E+03,
3665  -1.20839250E+03,
3666  7.75766802E+02,
3667  -2.05872410E+02};
3668 
3669  // !Chisq= 22. for 30 points
3670  // !Term Coeficient
3671  G4double C_COEF[3] = {
3672  1.69029097E-02,
3673  1.80889367E-02,
3674  5.04268396E-03};
3675 
3676  fitval = 1.;
3677  if (A < 2.5) return 0;
3678 
3679  if ( (X > 0.70) || (X < 0.0085) ) {
3680  // !Out of range of fit
3681  if (X < 0.0085) X_U =.0085;
3682  if (X > 0.70) X_U = 0.70;
3683  } else {
3684  X_U = X;
3685  }
3686 
3687  LN_C = C_COEF[0];
3688  for (G4int i = 1; i <= 2; i++) {
3689  LN_C = LN_C + C_COEF[i] * pow(log(X_U),i);
3690  }
3691 
3692  C = exp(LN_C);
3693 
3694  ALPHA = ALPHA_COEF[0];
3695  for (G4int i = 1; i <= 8; i++) {
3696  ALPHA = ALPHA + ALPHA_COEF[i] * pow(X_U,i);
3697  }
3698 
3699  fitval = C*pow(A,ALPHA);
3700  return fitval;
3701 }
3702 
3703 
3704 
3705 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
3706 G4double QweakSimEPEvent::AlloyScattering(G4double E_in, //MeV
3707  G4double Theta,
3708  G4int Zin,
3709  G4int Ain,
3710  G4double &fWeightN,
3711  G4double &Q2,
3712  G4double &E_out)
3713 {
3714  G4double M_p = 938.2796 * MeV;
3715  G4double M = M_p * Ain;
3716  G4double alpha_FS = 1.0/137.035999074; // Fine Structure Constant
3717  G4double m_e = 0.511 * MeV; // electron mass in MeV
3718 
3719 
3720  G4double STH = sin(Theta/2.);
3721  G4double CTH = cos(Theta/2.);
3722 
3723  E_out = E_in/(1+2*E_in/M*STH*STH);
3724  Q2 = 4.*E_in*E_out*STH*STH; // MeV^2
3725 
3726 
3727  // Mott cross section
3728  G4double sigMott = pow(Zin*0.719982/E_in*CTH/(STH*STH),2);
3729  sigMott = sigMott/(1+2*E_in/M*STH*STH)*10000; // 10000 is conversion from fm^2 to ub
3730 
3731  // Form factor
3732  G4int NUC_MODEL = 1; // USE Fourier-Bessel if available
3733  G4bool OUT_OF_RANGE = 0;
3734  //G4double eng = E_in/1e-3; // GeV // unused
3735  G4double qsq = Q2/1e6; // GeV^2
3736 
3737  G4double FF = 1;
3738 
3739  if (Ain < 3) {
3740  G4cout << "** No generator fo A < 3 ** " << G4endl;
3741  return 0;
3742  }
3743  else if (Ain >= 3) {
3744  OUT_OF_RANGE = 0;
3745  if (NUC_MODEL == 1) {
3746  FF = FF_BESSEL(Zin,Ain,qsq,OUT_OF_RANGE);// !only for some Nuclei,for limited Q2
3747  }
3748  if (OUT_OF_RANGE || (NUC_MODEL == 0)) { //!use if FF_BESSEL out of range
3749  if (Ain == 3) { // !3HE ! added 5/30/96 SER
3750  G4cout << "** No generator for He3 **" << G4endl;
3751  return 0;
3752  }
3753  else if (Ain <= 20) {
3754  FF = Fshell(Zin,Ain,qsq);
3755  }
3756  else { // !ia >20
3757  FF = Fgauss(Zin,Ain,qsq);
3758  }
3759  }
3760  }
3761 
3762  G4double xsec = sigMott*FF*FF;
3763  // Schwinger correction from Mo and Tsai Rev.Mod.Phys. 41 205 (1969) eq. (II.2)
3764  G4double delta_Schwinger = 0.0;
3765  if(SchwingerDeltaE!=0) {
3766  G4double FunctionofTheta = log (STH*STH) * log (CTH*CTH);
3767  delta_Schwinger = (-2.0*alpha_FS/pi) * ((log(E_in/SchwingerDeltaE) - 13.0/12.0)
3768  * (log(Q2/(m_e*m_e)) - 1.0) + 17.0/36.0 + FunctionofTheta/2.0);
3769  }
3770  xsec = xsec*(1.0 + delta_Schwinger);
3771  fWeightN = xsec*sin(Theta);
3772 
3773  //G4cout << E_in << "\t" << Theta << "\t" << Q2 << "\t"
3774  // << SchwingerDeltaE << "\t" << delta_Schwinger << G4endl;
3775 
3776  return xsec;
3777 }
3778 
3779 G4double QweakSimEPEvent::Fgauss(G4int /* iZ */ /* unused */, G4int iA, G4double T)
3780 {
3781  // cout << "** Calling Fgauss **" << endl;
3782 
3783  G4double fgauss = 0.;
3784 
3785  G4double avgA = iA;
3786  G4double radius = 1.07*pow(avgA,1./3.);
3787 
3788  //! from H. de Vries: Nuclear Charge Density Distributions from Elastic Electron
3789  //! Scattering. in Atomic Data and Nuclear Data Tables 36, 495(1987)
3790  //! 8/9/96
3791  if(iA==205)radius=5.470;
3792  if(iA==56) radius=3.729;
3793  if(iA==28) radius=3.085;
3794  if(iA==27) radius=3.035;
3795  if(iA==24) radius=3.08;// Mg-24
3796 
3797  G4double x2 = (T/pow(.197328,2))*pow(radius,2);
3798  G4double CHARR = (T/pow(.197328,2))*(2.4*2.4)/6.;
3799 
3800  //cout << " T (Q^2), CHARR: " << T << ", " << CHARR << endl;
3801  // cout << "** Fgauss CHARR:: " << CHARR << endl;
3802  if (CHARR<80) fgauss = exp(-CHARR)/(1.+x2/6.);
3803 
3804  return fgauss;
3805 }
3806 
3807 G4double QweakSimEPEvent::Fshell(G4int iZ, G4int iA, G4double T)
3808 {
3809  // cout << "** Calling Fshell **" << endl;
3810 
3811  G4double fshell = 0.;
3812  G4double avgA = iA;
3813  G4double radius = 1.07*pow(avgA,1./3.);
3814  //! from H. de Vries: Nuclear Charge Density Distributions from Elastic Electron
3815  //! Scattering. in Atomic Data and Nuclear Data Tables 36, 495(1987)
3816  //!8/9/96
3817  if(iA==16) radius=2.737;
3818  if(iA==15) radius=2.611;
3819  if(iA==14) radius=2.540;
3820  if(iA==12) radius=2.464;
3821  if(iA== 9) radius=2.519;
3822  if(iA== 7) radius=2.39;
3823  if(iA== 6) radius=2.56;
3824  if(iA== 4) radius=1.676;
3825  G4double x2 = (T/pow(.197328,2))*pow(radius,2);
3826  G4double alp = (iZ-2.)/3.;
3827  G4double CHARR = x2*(2.+3.*alp)/(12.+30.*alp);
3828 
3829  // cout << "** Fshell CHARR:: " << CHARR << endl;
3830  if (CHARR<80) fshell = exp(-CHARR)*(1.-alp*x2/(6.+15.*alp));
3831 
3832  return fshell;
3833 }
3834 
3835 G4double QweakSimEPEvent::FF_BESSEL (G4int iZ, G4int iA, G4double T,
3836  G4bool &OUT_OF_RANGE){
3837 
3838  // Calculates PWBA form factor for various nuclei using
3839  // Fourier-Bessel coefficients. (Inverse Fourier transform)
3840  //
3841  // Sources:
3842  // Technique-
3843  // 1)B. Dreher et al., Nuc. Phys. A235 219-248 (1974)
3844  // Data-
3845  // 1)H. de Vries et al., Atomic Data and Nuclear Data Tables 36,
3846  // 495-536 (1987) - Table IV pg. 520
3847  // 2)G. Fricke et al., Atomic Data and Nuclear Data Tables 60,
3848  // 177-285 (1995) - Table IX pg. 264
3849  //
3850  // Note: Can be used for nuclei with coeffients listed below
3851  // for a limited Q^2 range.
3852  //----------------------------------------------------------------
3853  // 2003/ 3/28 - Corrected for divide 0/0 when qm =0 - Steve Rock
3854  // 2016/10/25 - Updated references, included addition FB coeffients,
3855  // for more nuclei and isotopes of existing nuclei. - K. Bartlett
3856  // 2016/10/26 - Included a few more FB coefficients for various
3857  // isotopes.
3858 
3859  G4double a[20];
3860 
3861  //3He2 Helium-3
3862  G4double a_3He[] = {0.20020E-01, 0.41934E-01, 0.36254E-01, 0.17941E-01,
3863  0.46608E-02, 0.46834E-02, 0.52042E-02, 0.38280E-02,
3864  0.25661E-02, 0.14182E-02, 0.61390E-03, 0.22929E-03};
3865  //12C6 Carbon-12
3866  G4double a_12C[] = {0.15721E-01, 0.38732E-01, 0.36808E-01, 0.14671E-01,
3867  -0.43277E-02,-0.97752E-02,-0.68908E-02,-0.27631E-02,
3868  -0.63568E-03, 0.71809E-04, 0.18441E-03, 0.75066E-04,
3869  0.51069E-04, 0.14308E-04, 0.23170E-05, 0.68465E-06};
3870  //15N7 Nitrogen-15
3871  G4double a_15N[] = {0.25491E-01, 0.50618E-01, 0.29822E-01,-0.55196E-02,
3872  -0.15913E-01,-0.76184E-02,-0.23992E-02,-0.47940E-03};
3873 
3874  //16O8 Oxygen-16
3875  G4double a_16O[] = {0.20238E-01, 0.44793E-01, 0.33533E-01, 0.35030E-02,
3876  -0.12293E-01,-0.10329E-01,-0.34036E-02,-0.41627E-03,
3877  -0.94435E-03,-0.25571E-03, 0.23759E-03,-0.10603E-03,
3878  0.41480E-04, 0.00000E-03};
3879  //26Mg12 Magnesium-26
3880  G4double a_26Mg[] = {2.9280e-2, 5.5520e-2, 2.3910e-2,-1.3940e-2,-1.5840e-2,
3881  -3.5680e-3, 1.9670e-3, 1.4470e-3,-8.6820e-5,-1.4000e-4,
3882  7.0000e-5,-2.0000e-5, 3.0000e-6, 1.0000e-7,-2.0000e-7};
3883 
3884  //27Al13 Aluminum-27
3885  G4double a_27Al[] = {0.43418E-01, 0.60298E-01, 0.28950E-02,-0.23522E-01,
3886  -0.79791E-02, 0.23010E-02, 0.10794E-02, 0.12574E-03,
3887  -0.13021E-03, 0.56563E-04,-0.18011E-04, 0.42869E-05};
3888 
3889  //28Si14 Silicon-28
3890  G4double a_28Si[] = {0.33495E-01, 0.59533E-01, 0.20979E-01,-0.16900E-01,
3891  -0.14998E-01,-0.93248E-03, 0.33266E-02, 0.59244E-03,
3892  -0.40013E-03, 0.12242E-03,-0.12994E-04,-0.92784E-05,
3893  0.72595E-05,-0.42096E-05};
3894  //29Si14 Silicon-29
3895  G4double a_29Si[] = {3.3521e-02, 5.9679e-02, 2.0593e-02,-1.8646e-02,
3896  -1.6550e-02,-1.1922e-03, 2.8025e-03,-6.7353e-05,
3897  -3.4619e-04, 1.7611e-04,-5.7173e-06, 1.2371e-05};
3898  //30Si14 Silicon-30
3899  G4double a_30Si[] = {2.8397e-02, 5.4163e-02, 2.5167e-02,-1.2858e-02,
3900  -1.7592e-02,-4.6722e-03, 2.4804e-03, 1.4760e-03,
3901  -3.0168e-04, 4.8346e-05, 0.0000e+00,-5.1570e-06,
3902  3.0261e-06};
3903  //48Ti22 Titanium-48
3904  G4double a_48Ti[] = {2.7850e-02, 5.5432e-02, 2.6369e-02,-1.7091e-02,
3905  -2.1798e-02,-2.4889e-03, 7.6631e-03, 3.4554e-03,
3906  -6.7477e-04, 1.0764e-04,-1.6564e-06,-5.5566e-06};
3907  //50Ti22 Titanium-50
3908  G4double a_50Ti[] = {3.1818e-02, 5.8556e-02, 1.9637e-02,-2.4309e-02,
3909  -1.8748e-02, 3.3741e-03, 8.9961e-03, 3.7954e-03,
3910  -4.1238e-04, 1.2540e-04};
3911  //50Cr24 Chromium-50
3912  G4double a_50Cr[] = {3.9174e-02, 6.1822e-02, 6.8550e-03,-3.0170e-02,
3913  -9.8745e-03, 8.7944e-03, 6.8502e-03,-9.3609e-04,
3914  -2.4962e-03,-1.5361e-03,-7.3687e-04};
3915  //52Cr24 Chromium-52
3916  G4double a_52Cr[] = {0.39287e-1, 0.62477e-1, 0.62482e-2,-0.32885e-1,
3917  -0.10648e-1, 0.10520e-1, 0.85478e-2,-0.24003e-3,
3918  -0.20499e-2,-0.12001e-2,-0.56649e-3};
3919  //54Cr24 Chromium-54
3920  G4double a_54Cr[] = {3.9002e-02, 6.0305e-02, 4.5845e-03,-3.0723e-02,
3921  -9.1355e-03, 9.3251e-03, 6.0583e-03,-1.5602e-03,
3922  -7.6809e-04, 7.6809e-04,-3.4804e-04};
3923  //54Fe26 Iron-54
3924  G4double a_54Fe[] = {4.2339e-02, 6.4428e-02, 1.5840e-03,-3.5171e-02,
3925  -1.0116e-02, 1.2069e-02, 6.2230e-03,-1.2045e-03,
3926  1.3561e-04, 1.0428e-05,-1.6980e-05, 9.1817e-06,
3927  -3.9988e-06, 1.5731e-06,-5.7862e-07, 2.0186e-07,
3928  -6.7892e-08};
3929  //56Fe26 Iron-56
3930  G4double a_56Fe[] = {0.42018E-1, 0.62337E-1, 0.23995e-3,-0.32776E-1,
3931  -0.79941E-2, 0.10844E-1, 0.49123e-2,-0.22144e-2,
3932  -0.18146E-3, 0.37261E-3,-0.23296E-3, 0.11494E-3,
3933  -0.50596E-4, 0.20652E-4,-0.79428E-5, 0.28986E-5,
3934  -0.10075E-5};
3935  //58Fe26 Iron-58
3936  G4double a_58Fe[] = {4.1791e-02, 6.0524e-02,-1.4978e-03,-3.1183e-02,
3937  -5.8013e-03, 1.0611e-02, 4.1629e-03,-2.9045e-06,
3938  5.4106e-04,-3.8689e-04, 2.0514e-04,-9.5237e-05,
3939  4.0707e-05,-1.6346e-05, 6.2233e-06,-2.2568e-06,
3940  7.8077e-07};
3941  //63Cu29 Copper-63
3942  G4double a_63Cu[] = {4.5598e-02, 6.0706e-02,-7.8616e-03,-3.1638e-02,
3943  -1.4447e-03, 1.0953e-02, 4.2578e-03,-2.4224e-04,
3944  -3.0067e-04, 2.3903e-04,-1.2910e-04, 6.0195e-05,
3945  -2.5755e-05, 1.0332e-05,-3.9330e-06, 1.4254e-06,
3946  -4.9221e-07};
3947  //65Cu29 Copper-65
3948  G4double a_65Cu[] = {0.45444E-01, 0.59544E-01,-0.94968E-02,-0.31561e-01,
3949  0.22898E-03, 0.11189E-01, 0.37360E-02,-0.64873E-03,
3950  -0.51133E-03, 0.43765E-03,-0.24276E-03, 0.11507E-03,
3951  -0.49761E-04, 0.20140E-04,-0.76945E-05, 0.28055E-05,
3952  -0.97411E-06};
3953  //64Zn30 Zinc-64
3954  G4double a_64Zn[] = {0.47038E-01, 0.61536E-01,-0.90045E-02,-0.30669E-01,
3955  -0.78705E-03, 0.10034E-01, 0.14053E-02,-0.20640E-02,
3956  0.35105E-03, 0.27303E-04,-0.63811E-04, 0.40893E-04,
3957  -0.20311E-04, 0.88986E-05,-0.35849E-05, 0.13522E-05,
3958  -0.38635E-06};
3959  //66Zn30 Zinc-66
3960  G4double a_66Zn[] = {4.6991e-02, 6.0995e-02,-9.6693e-03,-3.0457e-02,
3961  -5.3435e-04, 9.7083e-03, 1.4091e-03,-7.0813e-04,
3962  2.0809e-04,-4.8275e-05, 7.2680e-06, 9.1369e-07,
3963  -1.4874e-06, 8.8831e-07,-4.1689e-07, 1.7283e-07,
3964  -6.5968e-08};
3965  //68Zn30 Zinc-68
3966  G4double a_68Zn[] = {4.6654e-02, 5.8827e-02,-1.2283e-02,-2.9865e-02,
3967  2.5669e-03, 1.0235e-02, 3.1861e-03,-1.7351e-04,
3968  -4.2979e-04, 3.3700e-04,-1.8435e-04, 8.7043e-05,
3969  -3.7612e-05, 1.5220e-05,-5.8282e-06, 2.1230e-06,
3970  -7.3709e-07};
3971  //70Zn30 Zinc-70
3972  G4double a_70Zn[] = {4.6362e-02, 5.7130e-02,-1.3877e-02,-3.0030e-02,
3973  3.5341e-03, 1.0113e-02, 4.1029e-03, 7.6469e-04,
3974  -1.0138e-03, 6.0837e-04,-2.9929e-04, 1.3329e-04,
3975  -5.5502e-05, 2.1893e-05,-8.2286e-06, 2.9559e-06,
3976  -1.0148e-06};
3977 
3978  //196Pt78 Platinum-196
3979  G4double a_196Pt[] = {0.50218e-1, 0.53722e-1,-0.35015e-1,-0.34588e-1,
3980  0.23564e-1, 0.14340e-1,-0.13270e-1,-0.51212e-2,
3981  0.56088e-2, 0.14890e-2,-0.10928e-2, 0.55662e-3,
3982  -0.50557e-4,-0.19708e-3, 0.24016e-3};
3983 
3984 
3985  // Change to units of fm**(-1)
3986  G4double qq = sqrt ( T ) / 0.197328; // used to be q
3987  // cout << "q(1/fm):: " << qq << endl;
3988 
3989  G4double ff_bessel=0;
3990 
3991  if( qq <= 0.005 ){
3992  OUT_OF_RANGE=false; //! line added 4/29/98
3993  ff_bessel = 1.00000;
3994  return 0;
3995  }
3996 
3997  ff_bessel=0.0;
3998  OUT_OF_RANGE=false;
3999  Float_t R_max=0.;
4000  G4int n=0;
4001 
4002  if(iA == 3 && iZ == 2){ // 3He2
4003  if(qq>10.10 ) OUT_OF_RANGE=true;
4004  R_max=5.;
4005  n = 12;
4006  for(G4int i=0;i<n;i++)
4007  a[i] = a_3He[i];
4008  }
4009  else if(iA == 12 && iZ == 6){ // 12C6
4010  if(qq>4.01 ) OUT_OF_RANGE=true;
4011  R_max=8.;
4012  n = 16;
4013  for(G4int i=0;i<n;i++)
4014  a[i] = a_12C[i];
4015  }
4016  else if(iA == 15 && iZ == 7){ // 15N7
4017  if(qq>3.17) OUT_OF_RANGE=true;
4018  R_max=7.;
4019  n = 8;
4020  for(G4int i=0;i<n;i++)
4021  a[i] = a_15N[i];
4022  }
4023  else if(iA == 16 && iZ == 8){ // 16O8
4024  if(qq>2.77) OUT_OF_RANGE=true;
4025  R_max=8.;
4026  n = 13;
4027  for(G4int i=0;i<n;i++)
4028  a[i] = a_16O[i];
4029  }
4030  else if(iA == 26 && iZ == 26){ // 26Mg12
4031  if(qq>3.0) OUT_OF_RANGE=true;
4032  //if(qq<0.60) OUT_OF_RANGE=true;
4033  R_max=8.;
4034  n = 15;
4035  for(G4int i=0;i<n;i++)
4036  a[i] = a_26Mg[i];
4037  }
4038  else if(iA == 27 && iZ == 13){ // 27Al13
4039  if(qq>2.70) OUT_OF_RANGE=true;
4040  //if(qq<0.47) OUT_OF_RANGE=true;
4041  R_max=7.;
4042  n = 12;
4043  for(G4int i=0;i<n;i++)
4044  a[i] = a_27Al[i];
4045  }
4046  else if(iA == 28 && iZ == 14){ // 28Si14
4047  if(qq>2.64) OUT_OF_RANGE=true;
4048  //if(qq<0.25) OUT_OF_RANGE=true;
4049  R_max = 8.;
4050  n = 14;
4051  for(G4int i=0;i<n;i++)
4052  a[i] = a_28Si[i];
4053  }
4054  else if(iA == 29 && iZ == 14){ // 29Si14
4055  if(qq>2.64) OUT_OF_RANGE=true;
4056  //if(qq<0.25) OUT_OF_RANGE=true;
4057  R_max = 8.;
4058  n = 12;
4059  for(G4int i=0;i<n;i++)
4060  a[i] = a_29Si[i];
4061  }
4062  else if(iA == 30 && iZ == 14){ // 30Si14
4063  if(qq>2.64) OUT_OF_RANGE=true;
4064  //if(qq<0.25) OUT_OF_RANGE=true;
4065  R_max = 8.5;
4066  n = 13;
4067  for(G4int i=0;i<n;i++)
4068  a[i] = a_30Si[i];
4069  }
4070  else if (iA == 48 && iZ == 22) { // 48Ti22
4071  if (qq>2.20) OUT_OF_RANGE=true;
4072  //if (qq<0.15) OUT_OF_RANGE=true;
4073  R_max=10.0;
4074  n=12;
4075  for(G4int i=0;i<n;i++)
4076  a[i] = a_48Ti[i];
4077  }
4078  else if (iA == 50 && iZ == 22) { // 50Ti22
4079  if (qq>2.20) OUT_OF_RANGE=true;
4080  //if (qq<0.15) OUT_OF_RANGE=true;
4081  R_max=9.5;
4082  n=10;
4083  for(G4int i=0;i<n;i++)
4084  a[i] = a_50Ti[i];
4085  }
4086  else if (iA == 50 && iZ == 24) { // 50Cr24
4087  if (qq>2.59) OUT_OF_RANGE=true;
4088  //if (qq<0.15) OUT_OF_RANGE=true;
4089  R_max=9.0;
4090  n=11;
4091  for(G4int i=0;i<n;i++)
4092  a[i] = a_50Cr[i];
4093  }
4094  else if (iA == 52 && iZ == 24) { // 52Cr24
4095  if (qq>2.59) OUT_OF_RANGE=true;
4096  //if (qq<0.15) OUT_OF_RANGE=true;
4097  R_max=9.0;
4098  n=11;
4099  for(G4int i=0;i<n;i++)
4100  a[i] = a_52Cr[i];
4101  }
4102  else if (iA == 54 && iZ == 24) { // 54Cr24
4103  if (qq>2.59) OUT_OF_RANGE=true;
4104  //if (qq<0.15) OUT_OF_RANGE=true;
4105  R_max=9.0;
4106  n=11;
4107  for(G4int i=0;i<n;i++)
4108  a[i] = a_54Cr[i];
4109  }
4110  else if(iA == 54 && iZ == 26){ // 54Fe26
4111  if(qq>2.22) OUT_OF_RANGE=true;
4112  //if(qq<0.51) OUT_OF_RANGE=true;
4113  R_max=9.0;
4114  n = 17;
4115  for(G4int i=0;i<n;i++)
4116  a[i] = a_54Fe[i];
4117  }
4118  else if(iA == 56 && iZ == 26){ // 56Fe26
4119  if(qq>2.22) OUT_OF_RANGE=true;
4120  //if(qq<0.51) OUT_OF_RANGE=true;
4121  R_max=9.0;
4122  n = 17;
4123  for(G4int i=0;i<n;i++)
4124  a[i] = a_56Fe[i];
4125  }
4126  else if(iA == 58 && iZ == 26){ // 58Fe26
4127  if(qq>2.22) OUT_OF_RANGE=true;
4128  //if(qq<0.51) OUT_OF_RANGE=true;
4129  R_max=9.0;
4130  n = 17;
4131  for(G4int i=0;i<n;i++)
4132  a[i] = a_58Fe[i];
4133  }
4134  else if(iA == 63 && iZ == 29){ // 63Cu29
4135  if(qq>2.22) OUT_OF_RANGE=true;
4136  //if(qq<0.51) OUT_OF_RANGE=true;
4137  R_max=9.;
4138  n = 17;
4139  for(G4int i=0;i<n;i++)
4140  a[i] = a_63Cu[i];
4141  }
4142  else if(iA == 65 && iZ == 29){ // 65Cu29
4143  if(qq>2.22) OUT_OF_RANGE=true;
4144  //if(qq<0.51) OUT_OF_RANGE=true;
4145  R_max=9.;
4146  n = 17;
4147  for(G4int i=0;i<n;i++)
4148  a[i] = a_65Cu[i];
4149  }
4150  else if(iA == 64 && iZ == 30){ // 64Zn30
4151  if(qq>2.22) OUT_OF_RANGE=true;
4152  //if(qq<0.51) OUT_OF_RANGE=true;
4153  R_max=9.;
4154  n = 17;
4155  for(G4int i=0;i<n;i++)
4156  a[i] = a_64Zn[i];
4157  }
4158  else if(iA == 66 && iZ == 30){ // 66Zn30
4159  if(qq>2.22) OUT_OF_RANGE=true;
4160  //if(qq<0.51) OUT_OF_RANGE=true;
4161  R_max=9.;
4162  n = 17;
4163  for(G4int i=0;i<n;i++)
4164  a[i] = a_66Zn[i];
4165  }
4166  else if(iA == 68 && iZ == 30){ // 68Zn30
4167  if(qq>2.22) OUT_OF_RANGE=true;
4168  //if(qq<0.51) OUT_OF_RANGE=true;
4169  R_max=9.;
4170  n = 17;
4171  for(G4int i=0;i<n;i++)
4172  a[i] = a_68Zn[i];
4173  }
4174  else if(iA == 70 && iZ == 30){ // 70Zn30
4175  if(qq>2.22) OUT_OF_RANGE=true;
4176  //if(qq<0.51) OUT_OF_RANGE=true;
4177  R_max=9.;
4178  n = 17;
4179  for(G4int i=0;i<n;i++)
4180  a[i] = a_70Zn[i];
4181  }
4182  else if(iA == 196 && iZ == 78){ // 196Pt78
4183  if(qq>2.28) OUT_OF_RANGE=true;
4184  if(qq<0.34) OUT_OF_RANGE=true;
4185  R_max=12.;
4186  n = 15;
4187  for(G4int i=0;i<n;i++)
4188  a[i] = a_196Pt[i];
4189  }
4190  else
4191  OUT_OF_RANGE=true;
4192 
4193  if(OUT_OF_RANGE || (R_max == 0.)){
4194  ff_bessel=0.;
4195  return 0;
4196  }
4197 
4198  G4double qmu = 3.14159265 / R_max;
4199  G4double ff=0.;
4200  G4double sinqR = sin(qq*R_max);
4201  G4double qmux,qm,qp;
4202 
4203  G4double PI4=12.56637062; //! 4*PI
4204  G4double PI2=6.28318531; //! 2*PI
4205 
4206  for(G4int j=0;j<n;j++){
4207  qmux = qmu * Float_t(j+1);
4208  qm = qq - qmux;
4209  qp = qq + qmux;
4210  if(fabs(qm)>1.E-6)
4211  ff = ff + a[j]*(pow((-1.),(j+1)))*sinqR/(qm*qp);
4212  else
4213  ff = ff + a[j]*pow(R_max,2)/(PI2*(j+1));
4214  }
4215 
4216  if((qq*R_max>1.E-20) && (ff<1.E20))
4217  ff_bessel = PI4/Float_t(iZ)/qq *ff;
4218  else
4219  ff_bessel = 0.;
4220 
4221  if (ff_bessel < 0) ff_bessel*=-1;
4222 
4223  return ff_bessel;
4224 }
4225 
4226 ////....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
4227 //
4228 // ---------------------------------------------------------
4229 // Calculates the asymmetry weighting factor for elastic
4230 // scattering from the proton. Asymmetry returned in ppm
4231 // ---------------------------------------------------------
4232 
4233 G4double QweakSimEPEvent::GetAsymmetry_EP(G4double theta, G4double energy)
4234 {
4235 
4236  const G4double Mp=0.938;
4237  const G4double gf=0.0000116637;
4238  const G4double alpha=1./137.;
4239  //const G4double s2tw=0.231;
4240  const G4double s2tw_low=0.2387;
4241  const G4double qwp=0.0713;
4242  const G4double qwn=-0.988;
4243  //const G4double NT=14.;
4244  //const G4double ZT=13.;
4245 
4246 // Radiative correction factors [from Musolf's Phys. Rep. 1994]
4247  //const G4double rpv=-0.054;
4248  //const G4double rnv=-0.0143;
4249  //const G4double rt0a=0.0;
4250  const G4double rt1a=-0.23;
4251  //const G4double r0a=0.023;
4252 
4253 // Check for minimum theta
4254  const G4double theta_min = 1.745329E-4;
4255  if (theta<theta_min)
4256  theta = theta_min;
4257 
4258 // The energy is in MeV after the cross section subroutine -
4259 // change back to GeV
4260  energy=energy/1000.0;
4261 
4262 // Calculate Q2
4263  G4double Q2_ep=4*energy*energy*sin(theta/2.0)*sin(theta/2.0);
4264  Q2_ep=Q2_ep/(1.0+2.0*energy/Mp*sin(theta/2.0)*sin(theta/2.0));
4265 
4266 // Proton Asymmetry
4267 // Kinematic Factors
4268  G4double tau=Q2_ep/4.0/(Mp*Mp);
4269  G4double epsilon=1.0/(1.0+2.0*(1.0+tau)*tan(theta/2.0)*tan(theta/2.0));
4270  G4double epsilonp=sqrt(tau*(1.0+tau)*(1.0-epsilon*epsilon));
4271 
4272 // Form factor calculations
4273 // Form Factors: Electromagetic
4274  G4double gvd=1.0/( (1.0+Q2_ep/0.71)*(1.0+Q2_ep/0.71) );
4275  G4double gepg=gvd;
4276  G4double gmpg=2.793*gvd;
4277  G4double geng=1.91*tau*gvd/(1.+5.6*tau);
4278  G4double gmng=-1.91*gvd;
4279 
4280 // Form Factors: Neutral-weak, Axial
4281 // Assume: Gs_E,M=0, G8_A,Gs_A=0
4282  //G4double gad=1/((1.0+Q2_ep/1.001/1.001)*(1.0+Q2_ep/1.001/1.001));
4283  //G4double gsa=-0.12/((1.+Q2_ep/1.06/1.06)*(1.+Q2_ep/1.06/1.06));
4284  G4double gt1a=1.2695/((1+Q2_ep/1.001/1.001)*(1+Q2_ep/1.001/1.001));
4285  //G4double g8a=0.0;
4286  //gsa=0.0;
4287 
4288 // Use SM radiative correction factors
4289 // Since we use SM values for Qw(p) and Qw(n)
4290 // we no longer need the radiative corrections
4291 // on those quantities.
4292 // gepz=(1.-4.*s2tw)*(1.+rpv)*gepg-(1.+rnv)*geng
4293 // gmpz=(1.-4.*s2tw)*(1.+rpv)*gmpg-(1.+rnv)*gmng
4294 
4295  G4double gepz=qwp*gepg+qwn*geng;
4296  G4double gmpz=qwp*gmpg+qwn*gmng;
4297  G4double gapz=-(1.0+rt1a)*gt1a;
4298 
4299  G4double asym=-gf/(4.0*pi*alpha*sqrt(2.0))*Q2_ep*1e6;
4300  asym=asym/(epsilon*gepg*gepg+tau*gmpg*gmpg);
4301  asym=asym*(epsilon*gepg*gepz+tau*gmpg*gmpz-(1.0-4.0*s2tw_low)*epsilonp*gmpg*gapz);
4302 
4303  return asym;
4304 }
4305 
4306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
4307 
4308 
4309 //---------------------------------------------------------
4310 // Calculates the asymmetry weighting factor for elastic
4311 // e-Al scattering using a formula from Donnelly. Returns
4312 // the asymmetry in units of ppm.
4313 //
4314 // Qw(p) value from J. Erler, Phys. Rev. D72, 073003 (2005)
4315 // Qw(n) value from using C_1u and C_1d values from PDG2008
4316 //
4317 //---------------------------------------------------------
4318 
4319 G4double QweakSimEPEvent::GetAsymmetry_AL (G4double theta, G4double energy)
4320 {
4321 // Needed constants
4322  const G4double gf=0.0000116637;
4323  const G4double alpha=1./137.;
4324  const G4double qwp=0.0713;
4325  const G4double qwn=-0.988;
4326  const G4double NT=14.;
4327  const G4double ZT=13.;
4328 
4329 // The energy is in MeV after the cross section subroutine -
4330 // change back to GeV
4331  energy = energy/1000.0;
4332 
4333 // Claculate Q2
4334  G4double Q2=4*energy*energy*sin(theta/2.)*sin(theta/2.);
4335 
4336 // Correction for Coulomb distortion for Al
4337 // Q2=Q2*(1.+(3.*ZT*0.197/137./(2.*energy*2.98)))*(1.+(3.*ZT*0.197/137./(2.*energy*2.98)));
4338 
4339 // Aluminum Asymmetry
4340  G4double asym=-gf/(4.*pi*alpha*sqrt(2.))*1e6*Q2*(qwp+qwn*NT/ZT);
4341 
4342  return asym;
4343 }
4344 
4345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
4346 
4347 //---------------------------------------------------------
4348 // Calculates the asymmetry weighting factor for elastic
4349 // e-Be scattering using a formula from Donnelly. Returns
4350 // the asymmetry in units of ppm.
4351 //
4352 // Qw(p) value from J. Erler, Phys. Rev. D72, 073003 (2005)
4353 // Qw(n) value from using C_1u and C_1d values from PDG2008
4354 //
4355 //---------------------------------------------------------
4356 
4357 G4double QweakSimEPEvent::GetAsymmetry_Be (G4double theta, G4double energy)
4358 {
4359  const G4double Mp=0.938;
4360  const G4double gf=0.0000116637;
4361  const G4double alpha=1./137.;
4362  const G4double qwp=0.0713;
4363  const G4double qwn=-0.988;
4364  const G4double NT=5.;
4365  const G4double ZT=4.;
4366 
4367 // The energy is in MeV after the cross section subroutine -
4368 // change back to GeV
4369  energy = energy/1000.;
4370 
4371 // Claculate Q2
4372  G4double Q2=4*energy*energy*sin(theta/2.)*sin(theta/2.);
4373  Q2=Q2/(1.+2.*energy/(9.*Mp)*sin(theta/2.)*sin(theta/2.));
4374 
4375 // Correction for Coulomb distortion
4376  Q2=Q2*(1.+(3./2.*sqrt(3./5.)*ZT*0.197/137./(energy*1.77)))*(1.+(3./2.*sqrt(3./5.)*ZT*0.197/137./(energy*1.77)));
4377 
4378 // Beryllium Asymmetry
4379  G4double asym=-gf/(4.*pi*alpha*sqrt(2.))*1e6*Q2*(qwp+qwn*NT/ZT);
4380 
4381  return asym;
4382 }
4383 
4384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
4385 //==============================================================================
4386 // Born_Approx_Nuclei_Asym
4387 //
4388 // K. Bartlett 2016/9/12
4389 //
4390 // Member function that does the first order Born approximation of the
4391 // parity-violating eletron scattering asymmetry. Essentially a copy of the
4392 // QweakSimEPEvent::GetAsymmetry_Al & QweakSimEPEvent::GetAsymmetry_Be.
4393 // However this member function is a bit more general, so that it can be used
4394 // with general Z/A from the other elements in the Aluminum window/dummy
4395 // targets. Note that is very much an approximation. Most of these elements
4396 // have diffraction minima that cause the asymmetry to deviate away from the
4397 // Born approximation. Large uncertainties will have to be used. Further
4398 // investigation of the uncertaintes for each element in the alloy will have
4399 // to be made.
4400 //
4401 // Formula form taken from:
4402 // Horowitz, "Parity violating elastic scattering from Al-27 and the Qweak
4403 // Measurement", Physical Review C 89, 045503 (2014)
4404 //
4405 // Original from:
4406 // Donnelly "Parity-Violating Electron Scattering", Prog. Particle Nuclear
4407 // Physics 24, pg.179
4408 //
4409 // Returns asymmetry in units of parts-per-million [ppm].
4410 //==============================================================================
4411 G4double QweakSimEPEvent::Born_Approx_Nuclei_Asym(G4double Z, G4double A, G4double energy, G4double eprime, G4double qq){
4412 
4413  //Define constants
4414  const G4double G_F = 1.1663787e-5; //Fermi coupling constant [GeV^-2] PDG-2014
4415  const G4double alpha = 1/137.035999074; //Fine structure constant PDG-2014
4416  const G4double Qw_p = 0.0721; //Horowitz paper
4417  const G4double Qw_n = -0.9878; //Horowitz paper
4418 
4419  //Convert beam energy/outgoing energy from MeV to GeV
4420  energy = energy/1000.0;
4421  eprime = eprime/1000.0;
4422  //Convert Q^2 from MeV^2 to GeV^2
4423  qq = qq/1000000.0;
4424 
4425  G4double A_0 = -G_F/(4*pi*alpha*sqrt(2.0));//[GeV^2]
4426  G4double born_asym = 1.0e6*A_0*qq*(Qw_p+Qw_n*((A-Z)/Z));//[ppm]
4427 
4428  return born_asym;
4429 }
4430 
4431 
4432 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
4433 // Calculates the asymmetry weighting factor for elastic
4434 // scattering from the neutron.
4435 
4436 G4double QweakSimEPEvent::GetAsymmetry_EN(G4double theta, G4double energy)
4437 {
4438 
4439  const G4double Mn= 0.939565; // Mp=0.938;
4440  const G4double gf=0.0000116637;
4441  const G4double alpha=1./137.;
4442  //const G4double s2tw=0.231;
4443  const G4double s2tw_low=0.2387;
4444  const G4double qwp=0.0713;
4445  const G4double qwn=-0.988;
4446  //const G4double NT=14.;
4447  //const G4double ZT=13.;
4448 
4449 // Radiative correction factors [from Musolf's Phys. Rep. 1994]
4450  //const G4double rpv=-0.054;
4451  //const G4double rnv=-0.0143;
4452  //const G4double rt0a=0.0;
4453  const G4double rt1a=-0.23;
4454  //const G4double r0a=0.023;
4455 
4456 // Check for minimum theta
4457  const G4double theta_min = 1.745329E-4;
4458  if (theta<theta_min)
4459  theta = theta_min;
4460 
4461 // The energy is in MeV after the cross section subroutine -
4462 // change back to GeV
4463  energy=energy/1000.0;
4464 
4465 // Calculate Q2
4466  G4double Q2_ep=4*energy*energy*sin(theta/2.0)*sin(theta/2.0);
4467  Q2_ep=Q2_ep/(1.0+2.0*energy/Mn*sin(theta/2.0)*sin(theta/2.0));
4468 
4469 // Proton Asymmetry
4470 // Kinematic Factors
4471  G4double tau=Q2_ep/4.0/(Mn*Mn);
4472  G4double epsilon=1.0/(1.0+2.0*(1.0+tau)*tan(theta/2.0)*tan(theta/2.0));
4473  G4double epsilonp=sqrt(tau*(1.0+tau)*(1.0-epsilon*epsilon));
4474 
4475 // Form factor calculations
4476 // Form Factors: Electromagetic
4477  G4double gvd=1.0/( (1.0+Q2_ep/0.71)*(1.0+Q2_ep/0.71) );
4478  G4double gepg=gvd;
4479  G4double gmpg=2.793*gvd;
4480  G4double geng=1.91*tau*gvd/(1.+5.6*tau);
4481  G4double gmng=-1.91*gvd;
4482  //G4double gesg= 0;
4483  //G4double gmsg= 0;
4484 
4485 // Form Factors: Neutral-weak, Axial
4486 // Assume: Gs_E,M=0, G8_A,Gs_A=0
4487  //G4double gad=1/((1.0+Q2_ep/1.001/1.001)*(1.0+Q2_ep/1.001/1.001));
4488  //G4double gsa=-0.12/((1.+Q2_ep/1.06/1.06)*(1.+Q2_ep/1.06/1.06));
4489  G4double gt1a=1.2695/((1+Q2_ep/1.001/1.001)*(1+Q2_ep/1.001/1.001));
4490  //G4double g8a=0.0;
4491  //gsa=0.0;
4492 
4493 // Use SM radiative correction factors
4494 // Since we use SM values for Qw(p) and Qw(n)
4495 // we no longer need the radiative corrections
4496 // on those quantities.
4497 // gepz=(1.-4.*s2tw)*(1.+rpv)*gepg-(1.+rnv)*geng
4498 // gmpz=(1.-4.*s2tw)*(1.+rpv)*gmpg-(1.+rnv)*gmng
4499 
4500  G4double gapz=-(1.0+rt1a)*gt1a;
4501 
4502 // G4double gepz=qwp*gepg+qwn*geng;
4503 // G4double gmpz=qwp*gmpg+qwn*gmng;
4504 // G4double asym=-gf/(4.0*pi*alpha*sqrt(2.0))*Q2_ep*1e6;
4505 // asym=asym/(epsilon*gepg*gepg+tau*gmpg*gmpg);
4506 // asym=asym*(epsilon*gepg*gepz+tau*gmpg*gmpz-(1.0-4.0*s2tw_low)*epsilonp*gmpg*gapz);
4507 // //std::cout<<"asym_proton="<<asym<<std::endl;
4508 
4509  G4double genz=qwp*geng+qwn*gepg;
4510  G4double gmnz=qwp*gmng+qwn*gmpg;
4511  G4double asym_neutron = -gf/(4.0*pi*alpha*sqrt(2.0))*Q2_ep*1e6;
4512  asym_neutron =asym_neutron /(epsilon*geng*geng+tau*gmng*gmng);
4513  asym_neutron =asym_neutron *(epsilon*geng*genz+tau*gmng*gmnz-(1.0-4.0*s2tw_low)*epsilonp*gmng*gapz);
4514  //std::cout<<"asym_neutron="<<asym_neutron<<std::endl<<std::endl;
4515 
4516  return asym_neutron;
4517 }
4518 
4519 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
4520 //
4521 //---------------------------------------------------------
4522 // Calculates the asymmetry weighting factor for inelastic
4523 // scattering from the proton. The simple approximation
4524 // A_LR = -100 (ppm/GeV^2) * Q^2
4525 // is used. Returns asymmetry in units of ppm.
4526 // Reference:
4527 // H.-W. Hammer and D. Drechsel, Z. Phys. A353, 321-331 (1995)
4528 //
4529 //---------------------------------------------------------
4530 
4531 G4double QweakSimEPEvent::GetAsymmetry_Pi(G4double Q2_pi)
4532 {
4533 // Inelastic Asymmetry
4534  Q2_pi = Q2_pi/1e6;
4535  G4double asym=-100*Q2_pi;
4536  return asym;
4537 }
4538 
4539 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
4540 
4541 
4542 ////....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
4543 // Jim Dowd
4544 // ---------------------------------------------------------
4545 // Setters and getters need to be defined here because
4546 // they call myUserInfo methods
4547 // ---------------------------------------------------------
4548 
4549 void QweakSimEPEvent::SetEPrime_Min(G4double energy) {myUserInfo->SetEPrime_Min(energy); CheckLookupTableBounds();}
4550 G4double QweakSimEPEvent::GetEPrime_Min() {return myUserInfo->GetEPrime_Min();}
4551 
4552 void QweakSimEPEvent::SetEPrime_Max(G4double energy) {myUserInfo->SetEPrime_Max(energy); CheckLookupTableBounds();}
4553 G4double QweakSimEPEvent::GetEPrime_Max() {return myUserInfo->GetEPrime_Max();}
4554 
4555 void QweakSimEPEvent::SetThetaAngle_Min(G4double ang) {myUserInfo->SetThetaAngle_Min(ang); CheckLookupTableBounds();}
4556 G4double QweakSimEPEvent::GetThetaAngle_Min() {return myUserInfo->GetThetaAngle_Min();}
4557 
4558 void QweakSimEPEvent::SetThetaAngle_Max(G4double ang) {myUserInfo->SetThetaAngle_Max(ang); CheckLookupTableBounds();}
4559 G4double QweakSimEPEvent::GetThetaAngle_Max() {return myUserInfo->GetThetaAngle_Max();}
4560 
4561 void QweakSimEPEvent::SetPhiAngle_Min(G4double ang) {myUserInfo->SetPhiAngle_Min(ang);}
4562 G4double QweakSimEPEvent::GetPhiAngle_Min() {return myUserInfo->GetPhiAngle_Min();}
4563 
4564 void QweakSimEPEvent::SetPhiAngle_Max(G4double ang) {myUserInfo->SetPhiAngle_Max(ang);}
4565 G4double QweakSimEPEvent::GetPhiAngle_Max() {return myUserInfo->GetPhiAngle_Max();}
4566 
4567 void QweakSimEPEvent::SetLuminosity(G4double lum) {myUserInfo->SetLuminosity(lum);}
4568 G4double QweakSimEPEvent::GetLuminosity() {return myUserInfo->GetLuminosity();}
4569 
4570 void QweakSimEPEvent::SetPhaseSpace(G4double ps) {myUserInfo->SetPhaseSpace(ps);}
4571 G4double QweakSimEPEvent::GetPhaseSpace() {return myUserInfo->GetPhaseSpace();}
4572 
4573 void QweakSimEPEvent::SetBeamEnergy(G4double energy) {
4574  if (energy>0) {
4575  CheckLookupTableBounds();
4576  myUserInfo->SetBeamEnergy(energy);
4577  }
4578  else {
4579  G4cout << G4endl << "##### Beam Energy must be greater than zero" << G4endl << G4endl;
4580  }
4581 }
4582 G4double QweakSimEPEvent::GetBeamEnergy() {return myUserInfo->GetBeamEnergy();}
4583 
4584 
4585 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
4586 //jpan@nuclear.uwinnipeg.ca
4587 // To determine the energy of the electron incident on the scattering
4588 // vertex inside the target, the energy loss up to the scattering
4589 // event (inside the target material and walls) must be simulated.
4590 // The way this is done here is to first determine the scattering vertex
4591 // position, and find the path length of the electron, then calculate the
4592 // collision loss and radiation loss, which are added up to total energy
4593 // loss.
4594 
4595 /*
4596 G4double QweakSimPrimaryGeneratorAction::PreTracking_Energy_Loss(G4double E_in, G4double Vertex_Z, G4double Target_A)
4597 {
4598  //Incident energy E_in in units of MeV
4599  G4double E_Loss_Collision, E_Loss_Radiation, E_Loss;
4600 
4601  //part 1. energy loss by collision
4602  // Bethe-Bloch formula for electron collision energy loss, taken from Leo's eqn 2.66 and 2.26
4603 
4604  //define constant factor: Avagadro's Na=6.02e23; r=e*e/m/c/c is the classical electron radius;
4605  G4double 2_Pi_Na_re2_me_c2 = 0.1535; //units:[MeV cm^2/g]
4606 
4607  G4double m = 0.511; //electron mass in [MeV/c^2]
4608  G4double rho_H2 = 0.071; //density of H2, [g/cm^3]
4609  G4double rho_Al = 2.70; //density of Aluminum
4610  G4double I_H2 = 19.0; // mean excitation potential for H2, in units of [eV],
4611  // I/Z = 12+7/Z [eV] if Z<13
4612  G4double I_Al = 166.0; // mean excitation potential for Aluminum, in units of [eV]
4613  // I/Z = 9.76+58.8*Z^(-1.19) [eV] if Z>=13
4614 
4615  // For relativistic electron: E = gamma*m*c^2
4616  // P = gamma*m*beta*c, E^2 = P^2*c^2 + m^2*c^4
4617  // => beta = sqrt(1-m^2/E^2)
4618  G4double beta = sqrt(1.0 - m*m/E_in/E_in) ;
4619  G4double gamma = E_in/m;
4620  G4double eta = beta*gamma;
4621 
4622  G4double tau = E_in/m; //kinetic energy of electron in units of [mc^2]
4623 
4624  //calculation for the density correction "delta", taken from Leo's 2.30
4625  // hnup = h*nu_p is the plasma frequency of the material
4626  G4double hnup_H2 = 28.816e-9 *sqrt(rho_H2*Z/A);
4627  G4double hnup_Al = 28.816e-9 *sqrt(rho_Al*Z/A);
4628 
4629  G4double C0_H2 = -(2.0*log(I_H2/(hnup_H2))+1.0);
4630  G4double delta_H2 = 4.6052*log10(eta)+C0_H2;
4631 
4632  G4double C0_Al = -(2.0*log(I_Al/(hnup_Al))+1.0); //C0_Al = -4.24
4633  G4double delta_Al = 4.6052*log10(eta)+C0_Al;
4634 
4635  //calculation for the shell correction "C", taken from Leo's 2.33
4636  G4double C_H2 = (0.422377/eta/eta +0.0304043/eta/eta/eta/eta -0.00038106/eta/eta/eta/eta/eta/eta)*I_H2*I_H2*1.0e-6
4637  + (3.850190/eta/eta -0.1667989/eta/eta/eta/eta +0.00157955/eta/eta/eta/eta/eta/eta)*I_H2*I_H2*1.0e-9;
4638 
4639  G4double C_Al = (0.422377/eta/eta +0.0304043/eta/eta/eta/eta -0.00038106/eta/eta/eta/eta/eta/eta)*I_Al*I_Al*1.0e-6
4640  + (3.850190/eta/eta -0.1667989/eta/eta/eta/eta +0.00157955/eta/eta/eta/eta/eta/eta)*I_Al*I_Al*1.0e-9;
4641 
4642  G4double F_tau = 1-beta*beta+(tau*tau/8.0-(2*r+1)*log(2.0))/(tau+1.0)/(tau+1.0);
4643  G4double dE_dx_H2 = 2_Pi_Na_re2_me_c2*rho_H2*Z*A/beta/beta*(log(tau*tau*(tau+2.0)/2/(I_H2/m/c/c)/(I_H2/m/c/c))+F_tau-delta-2.0*C_H2/Z);
4644  G4double dE_dx_Al = 2_Pi_Na_re2_me_c2*rho_Al*Z*A/beta/beta*(log(tau*tau*(tau+2.0)/2/(I_Al/m/c/c)/(I_Al/m/c/c))+F_tau-delta-2.0*C_Al/Z);
4645 
4646  //Part 2. energy loss by radiation: Bremsstrahlung
4647 
4648  E_Loss = E_Loss_Collision + E_Loss_Radiation;
4649  return E_Loss;
4650 
4651 }
4652 
4653 */
4654 
4655 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
QweakSimEPEvent::GetanEvent
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)
QweakSimEPEvent::ReactionRegion
G4int ReactionRegion
reaction region used for event generation
Definition: QweakSimEPEvent.hh:54
QweakSimEPEvent::EventGen_Messenger
QweakSimEPEventMessenger * EventGen_Messenger
Definition: QweakSimEPEvent.hh:87
QweakSimUserInformation::GetNormMomentum
G4ThreeVector GetNormMomentum() const
Definition: QweakSimUserInformation.hh:93
QweakSimEPEvent::SchwingerDeltaE
G4double SchwingerDeltaE
Definition: QweakSimEPEvent.hh:65
QweakSimEPEvent::Quasi_Elastic_Bosted
G4double Quasi_Elastic_Bosted(G4double E_in, G4double Theta, G4int Zin, G4int Ain, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:2525
QweakSimEPEvent::ReactionType
G4int ReactionType
reaction type used for event generation
Definition: QweakSimEPEvent.hh:53
QweakSimEPEvent::SetEPrime_Max
void SetEPrime_Max(G4double energy)
Definition: QweakSimEPEvent.cc:4552
QweakSimEPEvent::GetThetaAngle_Max
G4double GetThetaAngle_Max()
Definition: QweakSimEPEvent.cc:4559
QweakSimEPEvent::fLookupTable
QweakSimFieldMap< value_t, value_n > * fLookupTable
Definition: QweakSimEPEvent.hh:80
QweakSimEPEvent::SetEPrime_Min
void SetEPrime_Min(G4double energy)
Definition: QweakSimEPEvent.cc:4549
QweakSimUserInformation::SetThetaAngle_Max
void SetThetaAngle_Max(G4double ang)
Definition: QweakSimUserInformation.hh:278
QweakSimEPEvent::SetPhaseSpace
void SetPhaseSpace(G4double ps)
Definition: QweakSimEPEvent.cc:4570
QweakSimEPEvent::GetMomentumDirection
G4ThreeVector GetMomentumDirection()
QweakSimEPEvent::GetThetaAngle_Min
G4double GetThetaAngle_Min()
Definition: QweakSimEPEvent.cc:4556
Alinel_crsec_gen
G4double Alinel_crsec_gen(G4double E_in=1160, G4double th=8.0)
Definition: Alinel_crsec_gen.h:18
QweakSimEPEvent::Fshell
G4double Fshell(G4int Z, G4int A, G4double q2)
Definition: QweakSimEPEvent.cc:3807
QweakSimEPEvent::fitemc
G4double fitemc(G4double X, G4int A)
Definition: QweakSimEPEvent.cc:3626
QweakSimEPEvent::Moller_Scattering
G4double Moller_Scattering(G4double E_in, G4double theta1, G4double &E_out1, G4double &E_out2, G4double &theta2, G4double &q2, G4double &fWeightN, G4double &asymmetry)
Definition: QweakSimEPEvent.cc:1389
QweakSimEPEvent::value_n
static const G4int value_n
Definition: QweakSimEPEvent.hh:74
QweakSimEPEventMessenger
UI control for the event generator.
Definition: QweakSimEPEventMessenger.hh:40
QweakSimEPEvent::Pion_PhotoProductionAl
G4double Pion_PhotoProductionAl(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:2216
QweakSimEPEvent::AlNuclInel
G4double AlNuclInel(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:1894
QweakSimEPEvent::SuperElasticCheck
G4int SuperElasticCheck(G4double E_in, G4double E_out, G4double theta, G4double &xsec)
Definition: QweakSimEPEvent.cc:1505
QweakSimEPEvent::GetPhiAngle_Min
G4double GetPhiAngle_Min()
Definition: QweakSimEPEvent.cc:4562
QweakSimFieldMap
A multi-dimensional grid of values with interpolation methods.
Definition: QweakSimFieldMap.hh:51
QweakSimEPEvent::CGDR
G4double CGDR(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:2038
QweakSimUserInformation::SetEPrime_Max
void SetEPrime_Max(G4double energy)
Definition: QweakSimUserInformation.hh:272
QweakSimEPEvent::GetActiveOctantNumber
G4int GetActiveOctantNumber()
Definition: QweakSimEPEvent.hh:97
QweakSimEPEvent::Pion_PhotoProductionCarbon
G4double Pion_PhotoProductionCarbon(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:2320
QweakSimEPEvent::GetAsymmetry_Pi
G4double GetAsymmetry_Pi(G4double Q2_pi)
Definition: QweakSimEPEvent.cc:4531
QweakSimEPEvent::SetPhiAngle_Max
void SetPhiAngle_Max(G4double ang)
Definition: QweakSimEPEvent.cc:4564
QweakSimUserInformation::SetEPEvent
void SetEPEvent(QweakSimEPEvent *EP)
Definition: QweakSimUserInformation.hh:203
QweakSimEPEvent::F1F2QE09
void F1F2QE09(G4int Z, G4int A, G4double QSQ, G4double wsq, G4double &F1, G4double &F2)
Definition: QweakSimEPEvent.cc:2659
QweakSimEPEvent::Aluminum_Excited_State
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)
Definition: QweakSimEPEvent.cc:2377
QweakSimEPEvent::AlloyScattering
G4double AlloyScattering(G4double E_in, G4double Theta, G4int Zin, G4int Ain, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:3706
QweakSimEPEvent::Horowitz_DW_Xsect
G4double Horowitz_DW_Xsect(G4double E_in, G4double theta, const G4String path, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:2464
QweakSimEPEvent::M_p
static const G4double M_p
Definition: QweakSimEPEvent.hh:48
QweakSimEPEvent::GetAsymmetry_AL
G4double GetAsymmetry_AL(G4double theta, G4double energy)
Definition: QweakSimEPEvent.cc:4319
QweakSimEPEvent::SetThetaAngle_Min
void SetThetaAngle_Min(G4double ang)
Definition: QweakSimEPEvent.cc:4555
QweakSimEPEvent::ElasticPeakDeltaE
G4double ElasticPeakDeltaE
Definition: QweakSimEPEvent.hh:64
QweakSimEPEvent::CNuclInel
G4double CNuclInel(G4double E_in, G4double Theta, G4int nState, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:2092
QweakSimEPEvent::GetVertexZ
G4double GetVertexZ()
QweakSimEPEvent::Elastic_Cross_Section_Carbon
G4double Elastic_Cross_Section_Carbon(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:1986
QweakSimEPEvent::ResMod507
G4double ResMod507(G4int sf, G4double w2, G4double q2, G4double *xval)
Definition: QweakSimEPEvent.cc:1132
QweakSimUserInformation::GetOriginVertexPositionZ
G4double GetOriginVertexPositionZ() const
Definition: QweakSimUserInformation.hh:228
QweakSimEPEvent::GetBeamEnergy
G4double GetBeamEnergy()
Definition: QweakSimEPEvent.cc:4582
QweakSimEPEvent::GetAsymmetry_EP
G4double GetAsymmetry_EP(G4double theta, G4double energy)
Definition: QweakSimEPEvent.cc:4233
QweakSimEPEvent::Quasi_Elastic_Neutron
G4double Quasi_Elastic_Neutron(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:953
QweakSimEPEvent::Theta_Min
static const G4double Theta_Min
Definition: QweakSimEPEvent.hh:50
QweakSimEPEvent::SetBeamEnergy
void SetBeamEnergy(G4double energy=1.160 *GeV)
Definition: QweakSimEPEvent.cc:4573
QweakSimEPEvent::AlGDR
G4double AlGDR(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:1934
wiser_sigma
G4double wiser_sigma(G4double Ebeam, G4double pf, G4double thf, G4double rad_len, G4int type)
Definition: wiser_pion.h:12
QweakSimEPEvent::GetLuminosity
G4double GetLuminosity()
Definition: QweakSimEPEvent.cc:4568
QweakSimEPEvent::QweakSimEPEvent
QweakSimEPEvent(QweakSimUserInformation *myUI)
Definition: QweakSimEPEvent.cc:42
QweakSimEPEvent::Delta_Resonance
G4double Delta_Resonance(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:1008
QweakSimEPEvent::SetThetaAngle_Max
void SetThetaAngle_Max(G4double ang)
Definition: QweakSimEPEvent.cc:4558
QweakSimEPEvent::F1F2IN09
void F1F2IN09(G4int Z, G4int A, G4double QSQ, G4double wsq, G4double &F1, G4double &F2)
Definition: QweakSimEPEvent.cc:2974
QweakSimEPEvent::christy507
void christy507(G4double wsq, G4double Q2, G4double &F1, G4double &R, G4double &sigt, G4double &sigl)
Definition: QweakSimEPEvent.cc:3137
QweakSimUserInformation::SetEPrime
void SetEPrime(G4double energy)
Definition: QweakSimUserInformation.hh:292
QweakSimUserInformation::SetBeamEnergy
void SetBeamEnergy(G4double energy)
Definition: QweakSimUserInformation.hh:289
QweakSimEPEvent::Elastic_Cross_Section_Proton
G4double Elastic_Cross_Section_Proton(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:801
QweakSimEPEvent::Born_Approx_Nuclei_Asym
G4double Born_Approx_Nuclei_Asym(G4double Z, G4double A, G4double energy, G4double eprime, G4double qq)
Definition: QweakSimEPEvent.cc:4411
QweakSimEPEvent::SetLuminosity
void SetLuminosity(G4double lum)
Definition: QweakSimEPEvent.cc:4567
QweakSimEPEvent::NuclearInelastic_Bosted
G4double NuclearInelastic_Bosted(G4double E_in, G4double Theta, G4int Zin, G4int Ain, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:2593
QweakSimEPEvent::M_n
static const G4double M_n
Definition: QweakSimEPEvent.hh:47
QweakSimEPEvent::Isotropy
G4int Isotropy
isotropy used for event generation
Definition: QweakSimEPEvent.hh:57
QweakSimEPEvent::kActiveOctantNumber
G4int kActiveOctantNumber
Active octant number in the simulation, 0 will enable all octants.
Definition: QweakSimEPEvent.hh:55
QweakSimEPEvent::GetAsymmetry_Be
G4double GetAsymmetry_Be(G4double theta, G4double energy)
Definition: QweakSimEPEvent.cc:4357
QweakSimEPEvent::MEC2009
G4double MEC2009(G4int a, G4double q2, G4double w2)
Definition: QweakSimEPEvent.cc:3581
QweakSimEPEvent::GetPhaseSpace
G4double GetPhaseSpace()
Definition: QweakSimEPEvent.cc:4571
QweakSimUserInformation::GetNumberOfEventToBeProcessed
G4int GetNumberOfEventToBeProcessed() const
Definition: QweakSimUserInformation.hh:100
QweakSimEPEvent::GetAsymmetry_EN
G4double GetAsymmetry_EN(G4double theta, G4double energy)
Definition: QweakSimEPEvent.cc:4436
QweakSimUserInformation::SetEPrime_Min
void SetEPrime_Min(G4double energy)
Definition: QweakSimUserInformation.hh:269
QweakSimEPEvent::Horowitz_DW_Asym
G4double Horowitz_DW_Asym(G4double theta, const G4String path)
Definition: QweakSimEPEvent.cc:2506
QweakSimEPEvent::value_d
static const G4int value_d
Definition: QweakSimEPEvent.hh:75
QweakSimEPEvent::GetPhiAngle_Max
G4double GetPhiAngle_Max()
Definition: QweakSimEPEvent.cc:4565
QweakSimEPEvent::~QweakSimEPEvent
virtual ~QweakSimEPEvent()
Definition: QweakSimEPEvent.cc:78
QweakSimEPEvent::Elastic_Cross_Section_Aluminum
G4double Elastic_Cross_Section_Aluminum(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:858
QweakSimEPEvent::GetEPrime_Max
G4double GetEPrime_Max()
Definition: QweakSimEPEvent.cc:4553
QweakSimEPEvent::FF_BESSEL
G4double FF_BESSEL(G4int Z, G4int A, G4double q2, G4bool &ofr)
Definition: QweakSimEPEvent.cc:3835
QweakSimEPEvent::Radiative_Cross_Section_Lookup
const std::vector< G4double > Radiative_Cross_Section_Lookup(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:1432
QweakSimEPEvent::resmod507_v2
G4double resmod507_v2(G4double sf, G4double w2, G4double q2, G4double xval[50])
Definition: QweakSimEPEvent.cc:3373
QweakSimUserInformation::SetThetaAngle_Min
void SetThetaAngle_Min(G4double ang)
Definition: QweakSimUserInformation.hh:275
QweakSimEPEvent::myUserInfo
QweakSimUserInformation * myUserInfo
Definition: QweakSimEPEvent.hh:88
QweakSimEPEvent::GetEPrime_Min
G4double GetEPrime_Min()
Definition: QweakSimEPEvent.cc:4550
QweakSimEPEvent::Sigma_EEPrime
G4double Sigma_EEPrime(G4double eni, G4double eprime, G4double theta, G4double &q2)
Definition: QweakSimEPEvent.cc:1044
QweakSimEPEvent::Fgauss
G4double Fgauss(G4int Z, G4int A, G4double q2)
Definition: QweakSimEPEvent.cc:3779
QweakSimEPEvent::Pion_PhotoProduction
G4double Pion_PhotoProduction(G4double E_in, G4double Theta, G4double &fWeightN, G4double &Q2, G4double &E_out)
Definition: QweakSimEPEvent.cc:2169
QweakSimEPEvent::resmodd
void resmodd(G4double w2, G4double q2, G4double xval[50], G4double &sig)
Definition: QweakSimEPEvent.cc:3193
QweakSimEPEvent::SetPhiAngle_Min
void SetPhiAngle_Min(G4double ang)
Definition: QweakSimEPEvent.cc:4561