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 #include "G4BoldyshevTripletModel.hh"
00038 #include "G4PhysicalConstants.hh"
00039 #include "G4SystemOfUnits.hh"
00040
00041
00042
00043 using namespace std;
00044
00045
00046
00047 G4BoldyshevTripletModel::G4BoldyshevTripletModel(const G4ParticleDefinition*,
00048 const G4String& nam)
00049 :G4VEmModel(nam),fParticleChange(0),smallEnergy(4.*MeV),isInitialised(false),
00050 crossSectionHandler(0),meanFreePathTable(0)
00051 {
00052 lowEnergyLimit = 4.0*electron_mass_c2;
00053 highEnergyLimit = 100 * GeV;
00054 SetHighEnergyLimit(highEnergyLimit);
00055
00056 verboseLevel= 0;
00057
00058
00059
00060
00061
00062
00063
00064 if(verboseLevel > 0) {
00065 G4cout << "Triplet Gamma conversion is constructed " << G4endl
00066 << "Energy range: "
00067 << lowEnergyLimit / MeV << " MeV - "
00068 << highEnergyLimit / GeV << " GeV"
00069 << G4endl;
00070 }
00071 }
00072
00073
00074
00075 G4BoldyshevTripletModel::~G4BoldyshevTripletModel()
00076 {
00077 if (crossSectionHandler) delete crossSectionHandler;
00078 }
00079
00080
00081
00082 void
00083 G4BoldyshevTripletModel::Initialise(const G4ParticleDefinition*,
00084 const G4DataVector&)
00085 {
00086 if (verboseLevel > 3)
00087 G4cout << "Calling G4BoldyshevTripletModel::Initialise()" << G4endl;
00088
00089 if (crossSectionHandler)
00090 {
00091 crossSectionHandler->Clear();
00092 delete crossSectionHandler;
00093 }
00094
00095
00096
00097 crossSectionHandler = new G4CrossSectionHandler();
00098 crossSectionHandler->Initialise(0,lowEnergyLimit,100.*GeV,400);
00099 G4String crossSectionFile = "tripdata/pp-trip-cs-";
00100 crossSectionHandler->LoadData(crossSectionFile);
00101
00102
00103
00104 if (verboseLevel > 0) {
00105 G4cout << "Loaded cross section files for Livermore GammaConversion" << G4endl;
00106 G4cout << "To obtain the total cross section this should be used only " << G4endl
00107 << "in connection with G4NuclearGammaConversion " << G4endl;
00108 }
00109
00110 if (verboseLevel > 0) {
00111 G4cout << "Livermore Electron Gamma Conversion model is initialized " << G4endl
00112 << "Energy range: "
00113 << LowEnergyLimit() / MeV << " MeV - "
00114 << HighEnergyLimit() / GeV << " GeV"
00115 << G4endl;
00116 }
00117
00118 if(isInitialised) return;
00119 fParticleChange = GetParticleChangeForGamma();
00120 isInitialised = true;
00121 }
00122
00123
00124
00125 G4double
00126 G4BoldyshevTripletModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00127 G4double GammaEnergy,
00128 G4double Z, G4double,
00129 G4double, G4double)
00130 {
00131 if (verboseLevel > 3) {
00132 G4cout << "Calling ComputeCrossSectionPerAtom() of G4BoldyshevTripletModel"
00133 << G4endl;
00134 }
00135 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0;
00136
00137 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
00138 return cs;
00139 }
00140
00141
00142
00143 void G4BoldyshevTripletModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00144 const G4MaterialCutsCouple* ,
00145 const G4DynamicParticle* aDynamicGamma,
00146 G4double,
00147 G4double)
00148 {
00149
00150
00151
00152
00153 if (verboseLevel > 3)
00154 G4cout << "Calling SampleSecondaries() of G4BoldyshevTripletModel" << G4endl;
00155
00156 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
00157 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
00158
00159 G4double epsilon ;
00160 G4double p0 = electron_mass_c2;
00161
00162 G4double positronTotEnergy, electronTotEnergy, thetaEle, thetaPos;
00163 G4double ener_re=0., theta_re, phi_re, phi;
00164
00165
00166
00167 G4double energyThreshold = sqrt(2.)*electron_mass_c2;
00168 energyThreshold = 1.1*electron_mass_c2;
00169
00170
00171 G4double momentumThreshold_c = sqrt(energyThreshold * energyThreshold - electron_mass_c2*electron_mass_c2);
00172 G4double momentumThreshold_N = momentumThreshold_c/electron_mass_c2;
00173
00174
00175
00176 G4double SigmaTot = (28./9.) * std::log ( 2.* photonEnergy / electron_mass_c2 ) - 218. / 27. ;
00177 G4double X_0 = 2. * ( sqrt(momentumThreshold_N*momentumThreshold_N + 1) -1 );
00178 G4double SigmaQ = (82./27. - (14./9.) * log (X_0) + 4./15.*X_0 - 0.0348 * X_0 * X_0);
00179 G4double recoilProb = G4UniformRand();
00180
00181
00182 if (recoilProb >= SigmaQ/SigmaTot)
00183 {
00184
00185 G4double cosThetaMax = ( ( energyThreshold - electron_mass_c2 ) / (momentumThreshold_c) + electron_mass_c2*
00186 ( energyThreshold + electron_mass_c2 ) / (photonEnergy*momentumThreshold_c) );
00187
00188 if (cosThetaMax > 1) G4cout << "ERRORE " << G4endl;
00189
00190 G4double r1;
00191 G4double r2;
00192 G4double are, bre, loga, f1_re, greject, cost;
00193
00194 do {
00195 r1 = G4UniformRand();
00196 r2 = G4UniformRand();
00197
00198 cost = pow(cosThetaMax,r1);
00199 theta_re = acos(cost);
00200 are = 1./(14.*cost*cost);
00201 bre = (1.-5.*cost*cost)/(2.*cost);
00202 loga = log((1.+ cost)/(1.- cost));
00203 f1_re = 1. - bre*loga;
00204
00205 if ( theta_re >= 4.47*CLHEP::pi/180.)
00206 {
00207 greject = are*f1_re;
00208 } else {
00209 greject = 1. ;
00210 }
00211 } while(greject < r2);
00212
00213
00214
00215 G4double r3, r4, rt;
00216
00217 do {
00218
00219 r3 = G4UniformRand();
00220 r4 = G4UniformRand();
00221 phi_re = twopi*r3 ;
00222 G4double sint2 = 1. - cost*cost ;
00223 G4double fp = 1. - sint2*loga/(2.*cost) ;
00224 rt = (1.-cos(2.*phi_re)*fp/f1_re)/(2.*pi) ;
00225
00226 } while(rt < r4);
00227
00228
00229
00230 G4double S = electron_mass_c2*(2.* photonEnergy + electron_mass_c2);
00231 G4double D2 = 4.*S * electron_mass_c2*electron_mass_c2
00232 + (S - electron_mass_c2*electron_mass_c2)
00233 *(S - electron_mass_c2*electron_mass_c2)*sin(theta_re)*sin(theta_re);
00234 ener_re = electron_mass_c2 * (S + electron_mass_c2*electron_mass_c2)/sqrt(D2);
00235
00236
00237
00238
00239 G4double dxEle_re=sin(theta_re)*std::cos(phi_re),dyEle_re=sin(theta_re)*std::sin(phi_re), dzEle_re=cos(theta_re);
00240
00241 G4double electronRKineEnergy = std::max(0.,ener_re - electron_mass_c2) ;
00242
00243 G4ThreeVector electronRDirection (dxEle_re, dyEle_re, dzEle_re);
00244 electronRDirection.rotateUz(photonDirection);
00245
00246 G4DynamicParticle* particle3 = new G4DynamicParticle (G4Electron::Electron(),
00247 electronRDirection,
00248 electronRKineEnergy);
00249 fvect->push_back(particle3);
00250
00251 }
00252 else
00253 {
00254
00255
00256 fParticleChange->ProposeLocalEnergyDeposit(ener_re - electron_mass_c2);
00257 }
00258
00259
00260
00261
00262 G4double t = 0.5*log(momentumThreshold_N + sqrt(momentumThreshold_N*momentumThreshold_N+1));
00263
00264 G4cout << 0.5*asinh(momentumThreshold_N) << " " << t << G4endl;
00265
00266 G4double J1 = 0.5*(t*cosh(t)/sinh(t) - log(2.*sinh(t)));
00267 G4double J2 = (-2./3.)*log(2.*sinh(t)) + t*cosh(t)/sinh(t) + (sinh(t)-t*pow(cosh(t),3))/(3.*pow(sinh(t),3));
00268 G4double b = 2.*(J1-J2)/J1;
00269
00270 G4double n = 1 - b/6.;
00271 G4double re=0.;
00272 re = G4UniformRand();
00273 G4double a = 0.;
00274 G4double b1 = 16. - 3.*b - 36.*b*re*n + 36.*b*pow(re,2.)*pow(n,2.) +
00275 6.*pow(b,2.)*re*n;
00276 a = pow((b1/b),0.5);
00277 G4double c1 = (-6. + 12.*re*n + b + 2*a)*pow(b,2.);
00278 epsilon = (pow(c1,1./3.))/(2.*b) + (b-4.)/(2.*pow(c1,1./3.))+0.5;
00279
00280 G4double photonEnergy1 = photonEnergy - ener_re ;
00281 positronTotEnergy = epsilon*photonEnergy1;
00282 electronTotEnergy = photonEnergy1 - positronTotEnergy;
00283
00284 G4double momento_e = sqrt(electronTotEnergy*electronTotEnergy -
00285 electron_mass_c2*electron_mass_c2) ;
00286 G4double momento_p = sqrt(positronTotEnergy*positronTotEnergy -
00287 electron_mass_c2*electron_mass_c2) ;
00288
00289 thetaEle = acos((sqrt(p0*p0/(momento_e*momento_e) +1.)- p0/momento_e)) ;
00290 thetaPos = acos((sqrt(p0*p0/(momento_p*momento_p) +1.)- p0/momento_p)) ;
00291 phi = twopi * G4UniformRand();
00292
00293 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
00294 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
00295
00296
00297
00298
00299
00300
00301 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
00302
00303
00304
00305 G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
00306 electronDirection.rotateUz(photonDirection);
00307
00308 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
00309 electronDirection,
00310 electronKineEnergy);
00311
00312
00313 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
00314
00315
00316
00317 G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
00318 positronDirection.rotateUz(photonDirection);
00319
00320
00321 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
00322 positronDirection, positronKineEnergy);
00323
00324
00325
00326 fvect->push_back(particle1);
00327 fvect->push_back(particle2);
00328
00329
00330
00331
00332
00333 fParticleChange->SetProposedKineticEnergy(0.);
00334 fParticleChange->ProposeTrackStatus(fStopAndKill);
00335
00336 }
00337
00338
00339
00340
00341
00342