QwGeant4
wiser_pion.h File Reference
#include "TF1.h"
+ Include dependency graph for wiser_pion.h:
+ This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Macros

#define pi_w   3.14159
 
#define me_w   0.511E-3
 
#define alpha_w   0.007299
 
#define WISER_VERBOSE   0
 

Functions

G4double wiser_all_fit (G4double *x, G4double *par)
 
G4double wiser_sigma (G4double Ebeam, G4double pf, G4double thf, G4double rad_len, G4int type)
 

Macro Definition Documentation

#define alpha_w   0.007299

Definition at line 7 of file wiser_pion.h.

#define me_w   0.511E-3

Definition at line 6 of file wiser_pion.h.

#define pi_w   3.14159

Definition at line 5 of file wiser_pion.h.

#define WISER_VERBOSE   0

Definition at line 8 of file wiser_pion.h.

Referenced by wiser_sigma().

Function Documentation

G4double wiser_all_fit ( G4double *  x,
G4double *  par 
)

Definition at line 161 of file wiser_pion.h.

Referenced by wiser_sigma().

161  {
162  // Primary variable x[0] is photon energy in [GeV]
163 
164  // Parameters are:
165  // par[0] Beam energy [GeV]
166  //G4double Ebeam = par[0];
167  // par[1] Final particle momentum [GeV/c]
168  G4double pf = par[1];
169  // par[2] Final particle angle [rad]
170  G4double thf = par[2];
171  // par[3] Type (as defined in wiser_all_sig)
172  G4int type = (G4int) par[3];
173  // par[4] Minimum invariant mass of the residual system [GeV]
174  G4double M_X = par[4];
175 
176 
177  const G4double mass_p = 0.9383;
178  const G4double mass_p2 = mass_p*mass_p;
179  const G4double mass_pi = 0.13957;
180  const G4double mass_K = 0.4973;
181 
182  G4double mass[] = {mass_pi, mass_pi, mass_K, mass_K, mass_p, mass_p};
183 
184  G4double E_gamma = x[0];
185 
186  double s = mass_p2 + 2.0*E_gamma*mass_p;
187 
188  /* Wiser's fit pi+ pi- k+ k- p+ p- */
189  G4double A1[] = {566., 486., 368., 18.2, 1.33E5, 1.63E3 };
190  G4double A2[] = {829., 115., 1.91, 307., 5.69E4, -4.30E3};
191  G4double A3[] = {1.79, 1.77, 1.91, 0.98, 1.41, 1.79 };
192  G4double A4[] = {2.10, 2.18, 1.15, 1.83, .72, 2.24 };
193  G4double A6 = 1.90;
194  G4double A7 = -.0117;
195 
196 
197  // Boost to CoM
198  double beta_cm = E_gamma/(E_gamma+mass_p);
199  double gamma_cm = 1.0/sqrt(1.0 - beta_cm*beta_cm);
200 
201  double p_cm_z = -gamma_cm*beta_cm*sqrt(pf*pf+mass[type]*mass[type])
202  +gamma_cm*pf*cos(thf);
203 
204  double pT = pf*sin(thf);
205  double p_cm = sqrt( pT*pT + p_cm_z*p_cm_z );
206  double Ef = sqrt( p_cm*p_cm + mass[type]*mass[type] );
207 
208  double p_cm_max = sqrt(s +pow(M_X*M_X - mass[type]*mass[type],2.0)/s -
209  2.0*(M_X*M_X + mass[type]*mass[type]) )/2.0;
210 
211  double X_R = p_cm/p_cm_max;
212 
213  if( X_R > 1.0 ){ return 0.0; } // Kinematically forbidden
214 
215 
216  if( type != 4 ){ // Everything but proton
217  return ( A1[type] + A2[type]/sqrt(s) )*
218  pow(1.0 - X_R + A3[type]*A3[type]/s, A4[type])/E_gamma;
219  } else {
220  G4double U_MAN = fabs(2.0*mass_p2 - 2.0*mass_p*Ef);
221 
222  return ( A1[type] + A2[type]/sqrt(s) )*
223  pow(1.0 - X_R + A3[type]*A3[type]/s, A4[type])/pow(1.0 + U_MAN,A6+A7*s)/E_gamma;
224  }
225 
226  return 0.0;
227 }

+ Here is the caller graph for this function:

G4double wiser_sigma ( G4double  Ebeam,
G4double  pf,
G4double  thf,
G4double  rad_len,
G4int  type 
)

Definition at line 12 of file wiser_pion.h.

References wiser_all_fit(), and WISER_VERBOSE.

Referenced by QweakSimEPEvent::Pion_PhotoProduction(), QweakSimEPEvent::Pion_PhotoProductionAl(), and QweakSimEPEvent::Pion_PhotoProductionCarbon().

12  {
13  /*
14  Adapted from the infamous Wiser code from Peter Bosted
15 
16  Seamus Riordan
17  riordan@jlab.org
18  September 4, 2012
19 
20  type = 0 pi+
21  1 pi-
22  2 K+
23  3 K-
24  4 p
25  5 p-bar
26 
27  type defines the type of particle produced
28  Ebeam and pf are the beam energy in GeV and final particle momentum
29  in GeV/c respectively thf is the lab angle relative to the beam of
30  the produced particle in radians
31 
32  rad_len is the *effective* *total* radiator. [not in percent, typically
33  called bt in the literature] This means put in the effective radiator for
34  internal radiation BEFORE and AFTER the vertex (total 0.05 for JLab kinematics)
35  on top of real radiation lengths from ext brem times 4/3
36 
37  Returns the cross section in units of nanobarns/GeV/str
38  */
39 
40  const int ntype = 6;
41 
42  G4double A5[] = {-5.49, -5.23, -5.91, -4.45, -6.77, -6.53 };
43  G4double A6[] = {-1.73, -1.82, -1.74, -3.23, 1.90, -2.45 };
44 
45  const G4double mass_p = 0.9383;
46  const G4double mass_p2 = mass_p*mass_p;
47  const G4double mass_pi = 0.13957;
48  const G4double mass_K = 0.4973;
49  const G4double mass_Lambda = 1.116;
50 
51  G4double mass[] = {mass_pi, mass_pi, mass_K, mass_K, mass_p, mass_p};
52  G4double *mass2 = new G4double[ntype];
53 
54  int i;
55  for( i = 0; i < 6; i++ ){
56  mass2[i] = mass[i]*mass[i];
57  }
58 
59  // Calculate minimum photon energy required to produce such a system
60  // with the particle of transverse momentum pf*sin(thf)
61  // This is the case where in the CoM frame, the particle and the
62  // residual system are going back to back, transverse to the incoming particles
63  // and all components of the residual system are at rest relative to one another
64 
65  double M_X;// invariant mass of the minimum residual system
66 
67  switch( type ){
68  // pi+ can just be N
69  case 0:
70  M_X = mass_p;
71  break;
72  // pi- needs a pi+ N
73  case 1:
74  M_X = mass_p + mass_pi;
75  break;
76  // K+ can be created with just a Lambda, final state
77  // NOTE: This is different from the original Wiser
78  // code! They had just proton mass, but you still have
79  // to have an additional strange quark lying around
80  // The lowest energy configuration you can have like this
81  // is the Lambda
82  case 2:
83  M_X = mass_Lambda;
84  break;
85  // p/p-bar production is twice the mass of the proton
86  case 4:
87  case 5:
88  M_X = 2.0*mass_p;
89  break;
90  default:
91  fprintf(stderr, "%s: %s line %d - Forbidden type passed to Wiser parameterization\n",
92  __PRETTY_FUNCTION__, __FILE__, __LINE__ );
93  exit(1);
94  break;
95  }
96 
97  double MX2 = M_X*M_X;
98  double Ef = sqrt(pf*pf + mass2[type]);
99 
100  G4double E_gamma_min = ( MX2 - mass2[type] - mass_p2 + 2.0*mass_p*Ef )/
101  ( 2.0*( mass_p - Ef + pf*cos(thf) ) );
102 
103  // Parameterization parameters
104  double pT = pf*sin(thf);
105  double ML = sqrt(pT*pT + mass2[type]);
106 
107  double sigma, sig_e;
108  double fitres;
109 
110  if(WISER_VERBOSE)
111  G4cout <<"In wiser_pion.h: \n"<< Ebeam <<"\t" << pf << "\t" << thf <<"\t" << rad_len <<"\t" << type << G4endl;
112 
113  if(WISER_VERBOSE){
114  G4cout << "E_gamma_min:: " << E_gamma_min << G4endl;
115  G4cout << "Ebeam:: " << Ebeam << G4endl;
116  }
117 
118  if( E_gamma_min > 0.0 && E_gamma_min < Ebeam ){
119 
120  TF1 *wiserfit = new TF1("wiserfit", wiser_all_fit, E_gamma_min, Ebeam, 5);
121  wiserfit->SetParameter(0, Ebeam);
122  wiserfit->SetParameter(1, pf);
123  wiserfit->SetParameter(2, thf);
124  wiserfit->SetParameter(3, (G4double) type);
125  wiserfit->SetParameter(4, M_X);
126 
127  fitres = wiserfit->Integral(E_gamma_min, Ebeam);
128 
129  delete wiserfit;
130 
131  if( type != 4 ){
132  sig_e = fitres*exp(A5[type]*ML)*exp(A6[type]*pT*pT/Ef);
133  } else {
134  sig_e = fitres*exp(A5[type]*ML);
135  }
136 
137  // Factor of 1000 is required for units
138  sigma = pf*pf*sig_e*rad_len*1000.0/Ef;
139 
140  if(WISER_VERBOSE)
141  G4cout << "Pion crsec :: " << sigma << G4endl;
142 
143  delete mass2;
144  return sigma;
145  } else {
146  // Kinematically forbidden
147  if(WISER_VERBOSE)
148  G4cout << "--** Creation of Pions kinematically Forbidden **--" << G4endl;
149  delete mass2;
150  return 0.0;
151  }
152 
153  if(WISER_VERBOSE)
154  G4cout << "--** No pions. E_gamma is outside the valid range **--" << G4endl;
155 
156  delete mass2;
157  return 0.0;
158 }
#define WISER_VERBOSE
Definition: wiser_pion.h:8
G4double wiser_all_fit(G4double *x, G4double *par)
Definition: wiser_pion.h:161

+ Here is the call graph for this function:

+ Here is the caller graph for this function: