Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // Author: Sebastien Incerti 27 // 22 January 2012 28 // on base of G4BoldyshevTripletModel (original version) 29 // and G4LivermoreRayleighModel (MT version) 30 31 #include "G4BoldyshevTripletModel.hh" 32 #include "G4PhysicalConstants.hh" 33 #include "G4SystemOfUnits.hh" 34 #include "G4Log.hh" 35 #include "G4Exp.hh" 36 37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 38 39 using namespace std; 40 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 42 43 const G4int G4BoldyshevTripletModel::maxZ; 44 G4PhysicsFreeVector* G4BoldyshevTripletModel::data[] = {nullptr}; 45 46 G4BoldyshevTripletModel::G4BoldyshevTripletModel(const G4ParticleDefinition*, const G4String& nam) 47 :G4VEmModel(nam),smallEnergy(4.*MeV) 48 { 49 fParticleChange = nullptr; 50 51 lowEnergyLimit = 4.0*electron_mass_c2; 52 momentumThreshold_c = energyThreshold = xb = xn = lowEnergyLimit; 53 54 verboseLevel= 0; 55 // Verbosity scale for debugging purposes: 56 // 0 = nothing 57 // 1 = calculation of cross sections, file openings... 58 // 2 = entering in methods 59 60 if(verboseLevel > 0) 61 { 62 G4cout << "G4BoldyshevTripletModel is constructed " << G4endl; 63 } 64 } 65 66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 67 68 G4BoldyshevTripletModel::~G4BoldyshevTripletModel() 69 { 70 if(IsMaster()) { 71 for(G4int i=0; i<maxZ; ++i) { 72 if(data[i]) { 73 delete data[i]; 74 data[i] = nullptr; 75 } 76 } 77 } 78 } 79 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 81 82 void G4BoldyshevTripletModel::Initialise(const G4ParticleDefinition*, 83 const G4DataVector&) 84 { 85 if (verboseLevel > 1) 86 { 87 G4cout << "Calling Initialise() of G4BoldyshevTripletModel." 88 << G4endl 89 << "Energy range: " 90 << LowEnergyLimit() / MeV << " MeV - " 91 << HighEnergyLimit() / GeV << " GeV isMaster: " << IsMaster() 92 << G4endl; 93 } 94 // compute values only once 95 energyThreshold = 1.1*electron_mass_c2; 96 momentumThreshold_c = std::sqrt(energyThreshold * energyThreshold 97 - electron_mass_c2*electron_mass_c2); 98 G4double momentumThreshold_N = momentumThreshold_c/electron_mass_c2; 99 G4double t = 0.5*G4Log(momentumThreshold_N + 100 std::sqrt(momentumThreshold_N*momentumThreshold_N + 1.0)); 101 //G4cout << 0.5*asinh(momentumThreshold_N) << " " << t << G4endl; 102 G4double sinht = std::sinh(t); 103 G4double cosht = std::cosh(t); 104 G4double logsinht = G4Log(2.*sinht); 105 G4double J1 = 0.5*(t*cosht/sinht - logsinht); 106 G4double J2 = (-2./3.)*logsinht + t*cosht/sinht 107 + (sinht - t*cosht*cosht*cosht)/(3.*sinht*sinht*sinht); 108 109 xb = 2.*(J1-J2)/J1; 110 xn = 1. - xb/6.; 111 112 if(IsMaster()) 113 { 114 // Access to elements 115 const char* path = G4FindDataDir("G4LEDATA"); 116 117 G4ProductionCutsTable* theCoupleTable = 118 G4ProductionCutsTable::GetProductionCutsTable(); 119 120 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize(); 121 122 for(G4int i=0; i<numOfCouples; ++i) 123 { 124 const G4Material* material = 125 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial(); 126 const G4ElementVector* theElementVector = material->GetElementVector(); 127 std::size_t nelm = material->GetNumberOfElements(); 128 129 for (std::size_t j=0; j<nelm; ++j) 130 { 131 G4int Z = std::min((*theElementVector)[j]->GetZasInt(), maxZ); 132 if(!data[Z]) { ReadData(Z, path); } 133 } 134 } 135 } 136 if(!fParticleChange) { 137 fParticleChange = GetParticleChangeForGamma(); 138 } 139 } 140 141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 142 143 G4double 144 G4BoldyshevTripletModel::MinPrimaryEnergy(const G4Material*, 145 const G4ParticleDefinition*, 146 G4double) 147 { 148 return lowEnergyLimit; 149 } 150 151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 152 153 void G4BoldyshevTripletModel::ReadData(size_t Z, const char* path) 154 { 155 if (verboseLevel > 1) 156 { 157 G4cout << "Calling ReadData() of G4BoldyshevTripletModel" 158 << G4endl; 159 } 160 161 if(data[Z]) { return; } 162 163 const char* datadir = path; 164 165 if(!datadir) 166 { 167 datadir = G4FindDataDir("G4LEDATA"); 168 if(!datadir) 169 { 170 G4Exception("G4BoldyshevTripletModel::ReadData()", 171 "em0006",FatalException, 172 "Environment variable G4LEDATA not defined"); 173 return; 174 } 175 } 176 177 data[Z] = new G4PhysicsFreeVector(0,/*spline=*/true); 178 std::ostringstream ost; 179 ost << datadir << "/livermore/tripdata/pp-trip-cs-" << Z <<".dat"; 180 std::ifstream fin(ost.str().c_str()); 181 182 if( !fin.is_open()) 183 { 184 G4ExceptionDescription ed; 185 ed << "G4BoldyshevTripletModel data file <" << ost.str().c_str() 186 << "> is not opened!" << G4endl; 187 G4Exception("G4BoldyshevTripletModel::ReadData()", 188 "em0003",FatalException, 189 ed,"G4LEDATA version should be G4EMLOW6.27 or later."); 190 return; 191 } 192 193 else 194 { 195 196 if(verboseLevel > 3) { G4cout << "File " << ost.str() 197 << " is opened by G4BoldyshevTripletModel" << G4endl;} 198 199 data[Z]->Retrieve(fin, true); 200 } 201 202 // Activation of spline interpolation 203 data[Z]->FillSecondDerivatives(); 204 } 205 206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 207 208 G4double G4BoldyshevTripletModel::ComputeCrossSectionPerAtom( 209 const G4ParticleDefinition* part, 210 G4double GammaEnergy, G4double Z, G4double, G4double, G4double) 211 { 212 if (verboseLevel > 1) 213 { 214 G4cout << "Calling ComputeCrossSectionPerAtom() of G4BoldyshevTripletModel" 215 << G4endl; 216 } 217 218 if (GammaEnergy < lowEnergyLimit) { return 0.0; } 219 220 G4double xs = 0.0; 221 G4int intZ = std::max(1, std::min(G4lrint(Z), maxZ)); 222 G4PhysicsFreeVector* pv = data[intZ]; 223 224 // if element was not initialised 225 // do initialisation safely for MT mode 226 if(!pv) 227 { 228 InitialiseForElement(part, intZ); 229 pv = data[intZ]; 230 if(!pv) { return xs; } 231 } 232 // x-section is taken from the table 233 xs = pv->Value(GammaEnergy); 234 235 if(verboseLevel > 1) 236 { 237 G4cout << "*** Triplet conversion xs for Z=" << Z << " at energy E(MeV)=" 238 << GammaEnergy/MeV << " cs=" << xs/millibarn << " mb" << G4endl; 239 } 240 return xs; 241 } 242 243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 244 245 void G4BoldyshevTripletModel::SampleSecondaries( 246 std::vector<G4DynamicParticle*>* fvect, 247 const G4MaterialCutsCouple* /*couple*/, 248 const G4DynamicParticle* aDynamicGamma, 249 G4double, G4double) 250 { 251 252 // The energies of the secondary particles are sampled using 253 // a modified Wheeler-Lamb model (see PhysRevD 7 (1973), 26) 254 if (verboseLevel > 1) { 255 G4cout << "Calling SampleSecondaries() of G4BoldyshevTripletModel" 256 << G4endl; 257 } 258 259 G4double photonEnergy = aDynamicGamma->GetKineticEnergy(); 260 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection(); 261 262 G4double epsilon; 263 264 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine(); 265 266 // recoil electron thould be 3d particle 267 G4DynamicParticle* particle3 = nullptr; 268 static const G4double costlim = std::cos(4.47*CLHEP::pi/180.); 269 270 G4double loga, f1_re, greject, cost; 271 G4double cosThetaMax = (energyThreshold - electron_mass_c2 272 + electron_mass_c2*(energyThreshold + electron_mass_c2)/photonEnergy ) 273 /momentumThreshold_c; 274 if (cosThetaMax > 1.) { 275 //G4cout << "G4BoldyshevTripletModel::SampleSecondaries: ERROR cosThetaMax= " 276 // << cosThetaMax << G4endl; 277 cosThetaMax = 1.0; 278 } 279 280 G4double logcostm = G4Log(cosThetaMax); 281 do { 282 cost = G4Exp(logcostm*rndmEngine->flat()); 283 G4double are = 1./(14.*cost*cost); 284 G4double bre = (1.-5.*cost*cost)/(2.*cost); 285 loga = G4Log((1.+ cost)/(1.- cost)); 286 f1_re = 1. - bre*loga; 287 greject = (cost < costlim) ? are*f1_re : 1.0; 288 } while(greject < rndmEngine->flat()); 289 290 // Calculo de phi - elecron de recoil 291 G4double sint2 = (1. - cost)*(1. + cost); 292 G4double fp = 1. - sint2*loga/(2.*cost) ; 293 G4double rt, phi_re; 294 do { 295 phi_re = twopi*rndmEngine->flat(); 296 rt = (1. - std::cos(2.*phi_re)*fp/f1_re)/twopi; 297 } while(rt < rndmEngine->flat()); 298 299 // Calculo de la energia - elecron de recoil - relacion momento maximo <-> angulo 300 G4double S = electron_mass_c2*(2.* photonEnergy + electron_mass_c2); 301 G4double P2 = S - electron_mass_c2*electron_mass_c2; 302 303 G4double D2 = 4.*S * electron_mass_c2*electron_mass_c2 + P2*P2*sint2; 304 G4double ener_re = electron_mass_c2 * (S + electron_mass_c2*electron_mass_c2)/sqrt(D2); 305 306 if(ener_re >= energyThreshold) 307 { 308 G4double electronRKineEnergy = ener_re - electron_mass_c2; 309 G4double sint = std::sqrt(sint2); 310 G4ThreeVector electronRDirection (sint*std::cos(phi_re), sint*std::sin(phi_re), cost); 311 electronRDirection.rotateUz(photonDirection); 312 particle3 = new G4DynamicParticle (G4Electron::Electron(), 313 electronRDirection, 314 electronRKineEnergy); 315 } 316 else 317 { 318 // deposito la energia ener_re - electron_mass_c2 319 // G4cout << "electron de retroceso " << ener_re << G4endl; 320 fParticleChange->ProposeLocalEnergyDeposit(std::max(0.0, ener_re - electron_mass_c2)); 321 ener_re = 0.0; 322 } 323 324 // Depaola (2004) suggested distribution for e+e- energy 325 // VI: very suspect that 1 random number is not enough 326 // and sampling below is not correct - should be fixed 327 G4double re = rndmEngine->flat(); 328 329 G4double a = std::sqrt(16./xb - 3. - 36.*re*xn + 36.*re*re*xn*xn + 6.*xb*re*xn); 330 G4double c1 = G4Exp(G4Log((-6. + 12.*re*xn + xb + 2*a)*xb*xb)/3.); 331 epsilon = c1/(2.*xb) + (xb - 4.)/(2.*c1) + 0.5; 332 333 G4double photonEnergy1 = photonEnergy - ener_re ; 334 // resto al foton la energia del electron de retro. 335 G4double positronTotEnergy = std::max(epsilon*photonEnergy1, electron_mass_c2); 336 G4double electronTotEnergy = std::max(photonEnergy1 - positronTotEnergy, electron_mass_c2); 337 338 static const G4double a1 = 1.6; 339 static const G4double a2 = 0.5333333333; 340 G4double uu = -G4Log(rndmEngine->flat()*rndmEngine->flat()); 341 G4double u = (0.25 > rndmEngine->flat()) ? uu*a1 : uu*a2; 342 343 G4double thetaEle = u*electron_mass_c2/electronTotEnergy; 344 G4double sinte = std::sin(thetaEle); 345 G4double coste = std::cos(thetaEle); 346 347 G4double thetaPos = u*electron_mass_c2/positronTotEnergy; 348 G4double sintp = std::sin(thetaPos); 349 G4double costp = std::cos(thetaPos); 350 351 G4double phi = twopi * rndmEngine->flat(); 352 G4double sinp = std::sin(phi); 353 G4double cosp = std::cos(phi); 354 355 // Kinematics of the created pair: 356 // the electron and positron are assumed to have a symetric angular 357 // distribution with respect to the Z axis along the parent photon 358 359 G4double electronKineEnergy = electronTotEnergy - electron_mass_c2; 360 361 G4ThreeVector electronDirection (sinte*cosp, sinte*sinp, coste); 362 electronDirection.rotateUz(photonDirection); 363 364 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(), 365 electronDirection, 366 electronKineEnergy); 367 368 G4double positronKineEnergy = positronTotEnergy - electron_mass_c2; 369 370 G4ThreeVector positronDirection (-sintp*cosp, -sintp*sinp, costp); 371 positronDirection.rotateUz(photonDirection); 372 373 // Create G4DynamicParticle object for the particle2 374 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(), 375 positronDirection, positronKineEnergy); 376 // Fill output vector 377 378 fvect->push_back(particle1); 379 fvect->push_back(particle2); 380 381 if(particle3) { fvect->push_back(particle3); } 382 383 // kill incident photon 384 fParticleChange->SetProposedKineticEnergy(0.); 385 fParticleChange->ProposeTrackStatus(fStopAndKill); 386 } 387 388 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 389 390 #include "G4AutoLock.hh" 391 namespace { G4Mutex BoldyshevTripletModelMutex = G4MUTEX_INITIALIZER; } 392 393 void G4BoldyshevTripletModel::InitialiseForElement( 394 const G4ParticleDefinition*, G4int Z) 395 { 396 G4AutoLock l(&BoldyshevTripletModelMutex); 397 // G4cout << "G4BoldyshevTripletModel::InitialiseForElement Z= " 398 // << Z << G4endl; 399 if(!data[Z]) { ReadData(Z); } 400 l.unlock(); 401 } 402 403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 404