00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #ifndef G4ElasticHadrNucleusHE_h
00046 #define G4ElasticHadrNucleusHE_h 1
00047
00048 #include <vector>
00049
00050 #include "globals.hh"
00051 #include "G4ParticleDefinition.hh"
00052 #include "G4ParticleChange.hh"
00053 #include "G4Nucleus.hh"
00054 #include "G4HadronElastic.hh"
00055
00056 class G4NistManager;
00057
00058 static const G4int NHADRONS = 26;
00059 static const G4int ONQ0 = 5;
00060 static const G4int ONQ2 = 100;
00061 static const G4int NENERGY = 30;
00062 static const G4int NQTABLE = NENERGY*ONQ2;
00063
00064
00066
00067
00068
00069 class G4ElasticData
00070 {
00071 public:
00072
00073 G4ElasticData(const G4ParticleDefinition* h,
00074 G4int Z, G4double A, G4double* eGeV);
00075
00076 ~G4ElasticData(){}
00077
00078 const G4ParticleDefinition* Hadron() {return hadr;}
00079
00080 private:
00081 void DefineNucleusParameters(G4double A);
00082 const G4ParticleDefinition* hadr;
00083
00084
00085 G4ElasticData & operator=(const G4ElasticData &right);
00086 G4ElasticData(const G4ElasticData&);
00087
00088 public:
00089 G4int AtomicWeight;
00090 G4double R1, R2, Pnucl, Aeff;
00091 G4double limitQ2;
00092 G4double massGeV;
00093 G4double mass2GeV2;
00094 G4double massA;
00095 G4double massA2;
00096 G4int dnkE[NENERGY];
00097 G4double maxQ2[NENERGY];
00098
00099
00100 G4double TableQ2[ONQ2];
00101 G4double TableCrossSec[NQTABLE];
00102 };
00103
00105
00106
00107
00108 class G4ElasticHadrNucleusHE : public G4HadronElastic
00109 {
00110 public:
00111
00112 G4ElasticHadrNucleusHE(const G4String& name = "hElasticGlauber");
00113
00114 virtual ~G4ElasticHadrNucleusHE();
00115
00116 virtual G4double SampleInvariantT(const G4ParticleDefinition* p,
00117 G4double plab,
00118 G4int Z, G4int A);
00119
00120 virtual void Description() const;
00121
00122 G4double SampleT(const G4ParticleDefinition* p,
00123 G4double plab,
00124 G4int Z, G4int A);
00125
00126 G4double HadronNucleusQ2_2(G4ElasticData * pElD, G4int Z,
00127 G4double plabGeV, G4double tmax);
00128
00129 void DefineHadronValues(G4int Z);
00130
00131 G4double GetLightFq2(G4int Z, G4int A, G4double Q);
00132 G4double GetHeavyFq2(G4int Z, G4int Nucleus, G4double *LineFq2);
00133
00134 G4double GetQ2_2(G4int N, G4double * Q,
00135 G4double * F, G4double R);
00136
00137 G4double LineInterpol(G4double p0, G4double p2,
00138 G4double c1, G4double c2,
00139 G4double p);
00140
00141 G4double HadrNucDifferCrSec(G4int Z, G4int Nucleus, G4double Q2);
00142
00143 void InterpolateHN(G4int n, const G4double EnP[],
00144 const G4double C0P[], const G4double C1P[],
00145 const G4double B0P[], const G4double B1P[]);
00146
00147
00148 G4ElasticHadrNucleusHE & operator=(const G4ElasticHadrNucleusHE &right);
00149 G4ElasticHadrNucleusHE(const G4ElasticHadrNucleusHE&);
00150
00151 G4double GetBinomCof( G4int n, G4int m );
00152
00153 G4double GetFt(G4double Q2);
00154
00155 G4double GetDistrFun(G4double Q2);
00156
00157 G4double GetQ2(G4double Ran);
00158
00159 G4double HadronProtonQ2(const G4ParticleDefinition * aHadron,
00160 G4double inLabMom);
00161
00162 void GetKinematics(const G4ParticleDefinition * aHadron,
00163 G4double MomentumH);
00164 private:
00165
00166 void Binom();
00167
00168
00169
00170 G4int iHadrCode;
00171 G4int iHadron;
00172 G4int HadronCode[NHADRONS];
00173 G4int HadronType[NHADRONS];
00174 G4int HadronType1[NHADRONS];
00175
00176
00177
00178 G4double lowestEnergyLimit;
00179 G4double plabLowLimit;
00180 G4double dQ2;
00181
00182
00183
00184 G4double MbToGeV2;
00185 G4double sqMbToGeV;
00186 G4double Fm2ToGeV2;
00187 G4double GeV2;
00188 G4double protonM;
00189 G4double protonM2;
00190
00191
00192
00193 G4double hMass;
00194 G4double hMass2;
00195 G4double hLabMomentum;
00196 G4double hLabMomentum2;
00197 G4double MomentumCM;
00198 G4double HadrEnergy;
00199
00200
00201
00202 G4double R1, R2, Pnucl, Aeff;
00203 G4int NumbN;
00204
00205
00206
00207 G4double HadrTot, HadrSlope, HadrReIm, TotP,
00208 DDSect2, DDSect3, ConstU, FmaxT;
00209
00210
00211 G4double BoundaryP[7], BoundaryTL[7], BoundaryTG[7];
00212
00213
00214
00215 G4double Slope1, Slope2, Coeff1, Coeff2, MaxTR;
00216 G4double Slope0, Coeff0;
00217
00218 G4double aAIm, aDIm, Dtot11;
00219
00220 G4double Energy[NENERGY];
00221 G4double LowEdgeEnergy[NENERGY];
00222
00223 G4double SetBinom[240][240];
00224
00225 G4ElasticData* SetOfElasticData[NHADRONS][93];
00226 G4NistManager* nistManager;
00227
00228 };
00229
00231
00232 inline
00233 G4double G4ElasticHadrNucleusHE::LineInterpol(G4double p1, G4double p2,
00234 G4double c1, G4double c2,
00235 G4double p)
00236 {
00237
00238
00239
00240
00241 return c1+(p-p1)*(c2-c1)/(p2-p1);
00242 }
00243
00245
00246 inline
00247 void G4ElasticHadrNucleusHE::InterpolateHN(G4int n, const G4double EnP[],
00248 const G4double C0P[], const G4double C1P[],
00249 const G4double B0P[], const G4double B1P[])
00250 {
00251 G4int i;
00252
00253 for(i=1; i<n; i++) if(hLabMomentum <= EnP[i]) break;
00254
00255 if(i == n) i = n - 1;
00256
00257 Coeff0 = LineInterpol(EnP[i], EnP[i-1], C0P[i], C0P[i-1], hLabMomentum);
00258 Coeff1 = LineInterpol(EnP[i], EnP[i-1], C1P[i], C1P[i-1], hLabMomentum);
00259 Slope0 = LineInterpol(EnP[i], EnP[i-1], B0P[i], B0P[i-1], hLabMomentum);
00260 Slope1 = LineInterpol(EnP[i], EnP[i-1], B1P[i], B1P[i-1], hLabMomentum);
00261
00262
00263
00264 }
00265
00267
00268 inline
00269 G4double G4ElasticHadrNucleusHE::GetBinomCof( G4int numN, G4int numM )
00270 {
00271 if ( numN >= numM && numN <= 240) return SetBinom[numN][numM];
00272 else return 0.;
00273 }
00274
00276
00277 inline
00278 G4double G4ElasticHadrNucleusHE::GetDistrFun(G4double Q2)
00279 {
00280 return GetFt(Q2)/FmaxT;
00281 }
00282
00283 #endif