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
00046
00047 #include "G4QProbability.hh"
00048 #include "G4SystemOfUnits.hh"
00049
00050 G4QProbability::G4QProbability(G4int PDG)
00051 {
00052 S0 = 1.*GeV*GeV;
00053 pom_Alpha = 1.0808;
00054 pom_Alphaprime = 0.25/GeV/GeV;
00055 qex_Gamma = 9./GeV/GeV;
00056 qex_R2 = 27./GeV/GeV;
00057 qex_Alphaprime = 1.5/GeV/GeV;
00058
00059 G4int aP = std::abs(PDG);
00060 if (PDG==2212 || PDG==2112) InitForNucleon();
00061 else if(PDG==111 || aP==211) InitForPion();
00062 else if(PDG==130 || PDG==310 || aP==311 || aP==321) InitForKaon();
00063 else if(PDG==22) InitForGamma();
00064 else if(PDG > 3000) InitForHyperon();
00065 else if(PDG <-2000) InitForAntiBaryon();
00066 else
00067 {
00068 G4cout<<"-Warning-G4QProbability is initialized for PDGCode="<<PDG<<" as Pion"<<G4endl;
00069 InitForPion();
00070 }
00071 pom_sqC=std::sqrt(pom_C);
00072 }
00073
00074 G4double G4QProbability::GetCutPomProbability(const G4double s_value, const G4double imp2,
00075 const G4int nPom)
00076 {
00077 static const G4int nft=11;
00078 static const G4int nf1=nft-1;
00079 static const G4double ft[nft]={1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880.,3628800.};
00080 if(nPom<0) return 0.;
00081 G4double f=ft[nf1];
00082 if(nPom<nft) f=ft[nPom];
00083 else for(G4int i=nft; i<= nPom; i++) f*=i;
00084 G4double e=PomEikonal(s_value,imp2); e+=e;
00085 return std::exp(-e)*std::pow(e,nPom)/pom_C/f;
00086 }
00087
00088 G4double G4QProbability::GetCutQexProbability(const G4double s_value, const G4double imp2,
00089 const G4int nQex)
00090 {
00091 static const G4int nft=11;
00092 static const G4int nf1=nft-1;
00093 static const G4double ft[nft]={1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880.,3628800.};
00094 if(nQex<0) return 0.;
00095 G4double f=ft[nf1];
00096 if(nQex<nft) f=ft[nQex];
00097 else for(G4int i=nft; i<= nQex; i++) f*=i;
00098 G4double e=QexEikonal(s_value,imp2); e+=e;
00099 return std::exp(-e)*std::pow(e,nQex)/f;
00100 }
00101
00102 void G4QProbability::InitForNucleon()
00103 {
00104 pom_Gamma = 2.16/GeV/GeV;
00105 pom_C = 1.4;
00106 pom_R2 = 3.30/GeV/GeV;
00107 }
00108
00109 void G4QProbability::InitForHyperon()
00110 {
00111 pom_Gamma = 2.16/GeV/GeV;
00112 pom_C = 1.4;
00113 pom_R2 = 3.30/GeV/GeV;
00114 }
00115 void G4QProbability::InitForAntiBaryon()
00116 {
00117 pom_Gamma = 2.16/GeV/GeV;
00118 pom_C = 1.4;
00119 pom_R2 = 3.30/GeV/GeV;
00120 }
00121
00122 void G4QProbability::InitForPion()
00123 {
00124 pom_Gamma = 2.16/GeV/GeV;
00125 pom_C = 1.6;
00126 pom_R2 = 2.36/GeV/GeV;
00127 }
00128
00129 void G4QProbability::InitForKaon()
00130 {
00131 pom_Gamma = 1.92/GeV/GeV;
00132 pom_C = 1.8;
00133 pom_R2 = 1.96/GeV/GeV;
00134 }
00135
00136 void G4QProbability::InitForGamma()
00137 {
00138 pom_Gamma = 2.16/GeV/GeV;
00139 pom_C = 1.7;
00140 pom_R2 = 2.16/GeV/GeV;
00141 }
00142
00143 G4double G4QProbability::Expand(G4double z)
00144 {
00145 G4double sum=1.;
00146 G4double current=1.;
00147 for(G4int j=2; j<21; j++)
00148 {
00149 current *= -z*(j-1)/j/j;
00150 sum+=current;
00151 }
00152 return sum;
00153 }
00154
00155 G4double G4QProbability::GetQexTotProbability(const G4double s_value, const G4double imp2)
00156 {
00157 G4double ExpPom=std::exp(-PomEikonal(s_value,imp2));
00158 G4double ExpQex=std::exp(-QexEikonal(s_value,imp2));
00159 G4double Amp=(ExpQex*(1.-ExpPom) + sqr(pom_sqC-1.)*ExpPom*(1.-ExpQex))/pom_C;
00160 return Amp+Amp;
00161 }
00162
00163 G4double G4QProbability::GetQexElProbability(const G4double s_value, const G4double imp2)
00164 {
00165 G4double ExpPom=std::exp(-PomEikonal(s_value,imp2));
00166 G4double ExpQex=std::exp(-QexEikonal(s_value,imp2));
00167 G4double Amp=(ExpQex*(1.-ExpPom) + sqr(pom_sqC-1.)*ExpPom*(1.-ExpQex))/pom_C;
00168 return Amp*Amp;
00169 }
00170
00171 G4double G4QProbability::GetQexDubDiffProbability(const G4double s_value, const G4double imp2)
00172 {
00173 G4double ExpPom=std::exp(-PomEikonal(s_value,imp2));
00174 G4double ExpQex=std::exp(-QexEikonal(s_value,imp2));
00175 G4double Amp=sqr(pom_sqC-1.)*(ExpQex*(1.-ExpPom) + ExpPom*(1.-ExpQex))/pom_C;
00176 return Amp*Amp;
00177 }
00178
00179 G4double G4QProbability::GetQexSinDiffProbability(const G4double s_value, const G4double imp2)
00180 {
00181 G4double ExpPom=std::exp(-PomEikonal(s_value,imp2));
00182 G4double ExpQex=std::exp(-QexEikonal(s_value,imp2));
00183 G4double Amp=(pom_sqC-1.)*(ExpQex*(1.-ExpPom) - (pom_sqC-1.)*ExpPom*(1.-ExpQex))/pom_C;
00184 return Amp*Amp;
00185 }
00186
00187 G4double G4QProbability::GetQexDiffProbability(const G4double s_value, const G4double imp2)
00188 {
00189 return GetQexDubDiffProbability(s_value,imp2)+2*GetQexSinDiffProbability(s_value,imp2);
00190 }
00191
00192 G4double G4QProbability::GetQexInelProbability(const G4double s_value, const G4double imp2)
00193 {
00194 G4double ExpPom=std::exp(-PomEikonal(s_value,imp2));
00195 G4double ExpQex=std::exp(-QexEikonal(s_value,imp2));
00196 G4double Amp=sqr(pom_sqC-1.)*(ExpQex*(1.-ExpPom) + ExpPom*(1.-ExpQex))/pom_C;
00197 return Amp+Amp-Amp*Amp;
00198 }