Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // ------------------------------------------- 28 // 29 // GEANT4 Class 30 // File name: G4PAIPhotModel.cc 31 // 32 // Author: Vladimir.Grichine@cern.ch on base o 33 // 34 // Creation date: 07.10.2013 35 // 36 // Modifications: 37 // 38 // 39 40 #include "G4PAIPhotModel.hh" 41 42 #include "G4SystemOfUnits.hh" 43 #include "G4PhysicalConstants.hh" 44 #include "G4Region.hh" 45 #include "G4ProductionCutsTable.hh" 46 #include "G4MaterialCutsCouple.hh" 47 #include "G4MaterialTable.hh" 48 #include "G4RegionStore.hh" 49 50 #include "Randomize.hh" 51 #include "G4Electron.hh" 52 #include "G4Positron.hh" 53 #include "G4Gamma.hh" 54 #include "G4Poisson.hh" 55 #include "G4Step.hh" 56 #include "G4Material.hh" 57 #include "G4DynamicParticle.hh" 58 #include "G4ParticleDefinition.hh" 59 #include "G4ParticleChangeForLoss.hh" 60 #include "G4PAIPhotData.hh" 61 #include "G4DeltaAngle.hh" 62 63 ////////////////////////////////////////////// 64 65 using namespace std; 66 67 G4PAIPhotModel::G4PAIPhotModel(const G4Particl 68 : G4VEmModel(nam),G4VEmFluctuationModel(nam) 69 fVerbose(0), 70 fModelData(nullptr), 71 fParticle(nullptr) 72 { 73 fElectron = G4Electron::Electron(); 74 fPositron = G4Positron::Positron(); 75 76 fParticleChange = nullptr; 77 78 if(p) { SetParticle(p); } 79 else { SetParticle(fElectron); } 80 81 // default generator 82 SetAngularDistribution(new G4DeltaAngle()); 83 fLowestTcut = 12.5*CLHEP::eV; 84 } 85 86 ////////////////////////////////////////////// 87 88 G4PAIPhotModel::~G4PAIPhotModel() 89 { 90 if(IsMaster()) { delete fModelData; fModelDa 91 } 92 93 ////////////////////////////////////////////// 94 95 void G4PAIPhotModel::Initialise(const G4Partic 96 const G4DataVector& cuts) 97 { 98 if(fVerbose > 1) 99 { 100 G4cout<<"G4PAIPhotModel::Initialise for "< 101 } 102 SetParticle(p); 103 fParticleChange = GetParticleChangeForLoss() 104 105 if( IsMaster() ) 106 { 107 delete fModelData; 108 fMaterialCutsCoupleVector.clear(); 109 110 G4double tmin = LowEnergyLimit()*fRatio; 111 G4double tmax = HighEnergyLimit()*fRatio; 112 fModelData = new G4PAIPhotData(tmin, tmax, 113 114 // Prepare initialization 115 const G4MaterialTable* theMaterialTable = 116 size_t numOfMat = G4Material::GetNumberO 117 size_t numRegions = fPAIRegionVector.size( 118 119 // protect for unit tests 120 if(0 == numRegions) { 121 G4Exception("G4PAIModel::Initialise()"," 122 "no G4Regions are registered 123 fPAIRegionVector.push_back(G4RegionStore 124 ->GetRegion("DefaultRegionForTheWorld 125 numRegions = 1; 126 } 127 128 for( size_t iReg = 0; iReg < numRegions; + 129 { 130 const G4Region* curReg = fPAIRegionVecto 131 G4Region* reg = const_cast<G4Region*>(cu 132 133 for(size_t jMat = 0; jMat < numOfMat; ++ 134 { 135 G4Material* mat = (*theMaterialTable)[jMat]; 136 const G4MaterialCutsCouple* cutCouple = reg- 137 if(nullptr != cutCouple) 138 { 139 if(fVerbose > 1) 140 { 141 G4cout << "Reg <" <<curReg->GetName() << 142 << mat->GetName() << "> fCouple= " 143 << cutCouple << ", idx= " << cutCouple-> 144 <<" " << p->GetParticleName() 145 <<", cuts.size() = " << cuts.size() << G 146 } 147 // check if this couple is not already ini 148 149 size_t n = fMaterialCutsCoupleVector.size( 150 151 G4bool isnew = true; 152 if( 0 < n ) 153 { 154 for(size_t i=0; i<fMaterialCutsCoupleVec 155 { 156 if(cutCouple == fMaterialCutsCoupleVec 157 isnew = false; 158 break; 159 } 160 } 161 } 162 // initialise data banks 163 if(isnew) { 164 fMaterialCutsCoupleVector.push_back(cutC 165 G4double deltaCutInKinEnergy = cuts[cutC 166 fModelData->Initialise(cutCouple, deltaC 167 } 168 } 169 } 170 } 171 InitialiseElementSelectors(p, cuts); 172 } 173 } 174 175 ////////////////////////////////////////////// 176 177 void G4PAIPhotModel::InitialiseLocal(const G4P 178 G4VEmModel* masterModel) 179 { 180 fModelData = static_cast<G4PAIPhotModel*>(ma 181 fMaterialCutsCoupleVector = static_cast<G4PA 182 SetElementSelectors( masterModel->GetElement 183 } 184 185 ////////////////////////////////////////////// 186 187 G4double G4PAIPhotModel::MinEnergyCut(const G4 188 const G4MaterialCutsCouple*) 189 { 190 return fLowestTcut; 191 } 192 193 ////////////////////////////////////////////// 194 195 G4double G4PAIPhotModel::ComputeDEDXPerVolume( 196 const G4ParticleDefinition* p, 197 G4double kineticEnergy, 198 G4double cutEnergy) 199 { 200 G4int coupleIndex = FindCoupleIndex(CurrentC 201 if(0 > coupleIndex) { return 0.0; } 202 203 G4double cut = std::min(MaxSecondaryEnergy(p 204 G4double scaledTkin = kineticEnergy*fRatio; 205 G4double dedx = fChargeSquare*fModelData->DE 206 return dedx; 207 } 208 209 ////////////////////////////////////////////// 210 211 G4double G4PAIPhotModel::CrossSectionPerVolume 212 const G4ParticleDefinition* p, 213 G4double kineticEnergy, 214 G4double cutEnergy, 215 G4double maxEnergy ) 216 { 217 G4int coupleIndex = FindCoupleIndex(CurrentC 218 if(0 > coupleIndex) { return 0.0; } 219 220 G4double tmax = std::min(MaxSecondaryEnergy( 221 if(tmax <= cutEnergy) { return 0.0; } 222 223 G4double scaledTkin = kineticEnergy*fRatio; 224 G4double xs = fChargeSquare*fModelData->Cros 225 scal 226 return xs; 227 } 228 229 ////////////////////////////////////////////// 230 // 231 // It is analog of PostStepDoIt in terms of se 232 // 233 234 void G4PAIPhotModel::SampleSecondaries(std::ve 235 const G4MaterialCutsCouple* matCC, 236 const G4DynamicParticle* dp, 237 G4double tmin, 238 G4double maxEnergy) 239 { 240 G4int coupleIndex = FindCoupleIndex(matCC); 241 if(0 > coupleIndex) { return; } 242 243 SetParticle(dp->GetDefinition()); 244 245 G4double kineticEnergy = dp->GetKineticEnerg 246 247 G4double tmax = MaxSecondaryEnergy(fParticle 248 249 if( maxEnergy < tmax) tmax = maxEnergy; 250 if( tmin >= tmax) return; 251 252 G4ThreeVector direction = dp->GetMomentumDir 253 G4double scaledTkin = kineticEnergy*fRat 254 G4double totalEnergy = kineticEnergy + fM 255 G4double totalMomentum = sqrt(kineticEnergy 256 G4double plRatio = fModelData->GetPla 257 258 if( G4UniformRand() <= plRatio ) 259 { 260 G4double deltaTkin = fModelData->SamplePos 261 262 // G4cout<<"G4PAIPhotModel::SampleSecondar 263 // <<"; Tkin = "<<kineticEnergy/keV<<" keV 264 265 if( deltaTkin <= 0. && fVerbose > 0) 266 { 267 G4cout<<"G4PAIPhotModel::SampleSecondary 268 } 269 if( deltaTkin <= 0.) { return; } 270 271 if( deltaTkin > tmax) { deltaTkin = tmax; 272 273 const G4Element* anElement = SelectTargetA 274 275 G4int Z = anElement->GetZasInt(); 276 277 auto deltaRay = new G4DynamicParticle(fEle 278 GetAngularDistribution()->SampleDirectio 279 Z, matCC->GetMaterial()), 280 deltaTkin); 281 282 // primary change 283 284 kineticEnergy -= deltaTkin; 285 286 if( kineticEnergy <= 0. ) // kill primary 287 { 288 fParticleChange->SetProposedKineticEnerg 289 fParticleChange->ProposeLocalEnergyDepos 290 return; 291 } 292 else 293 { 294 G4ThreeVector dir = totalMomentum*direct 295 direction = dir.unit(); 296 fParticleChange->SetProposedKineticEnerg 297 fParticleChange->SetProposedMomentumDire 298 vdp->push_back(deltaRay); 299 } 300 } 301 else // secondary X-ray CR photon 302 { 303 G4double deltaTkin = fModelData->Sampl 304 305 // G4cout<<"PAIPhotonModel PhotonPostStep 306 307 if( deltaTkin <= 0. ) 308 { 309 G4cout<<"G4PAIPhotonModel::SampleSeconda 310 } 311 if( deltaTkin <= 0.) return; 312 313 if( deltaTkin >= kineticEnergy ) // stop p 314 { 315 deltaTkin = kineticEnergy; 316 kineticEnergy = 0.0; 317 } 318 G4double costheta = 0.; // G4UniformRand() 319 G4double sintheta = sqrt((1.+costheta)*(1. 320 321 // direction of the 'Cherenkov' photon 322 G4double phi = twopi*G4UniformRand(); 323 G4double dirx = sintheta*cos(phi), diry = 324 325 G4ThreeVector deltaDirection(dirx,diry,dir 326 deltaDirection.rotateUz(direction); 327 328 if( kineticEnergy > 0.) // primary change 329 { 330 kineticEnergy -= deltaTkin; 331 fParticleChange->SetProposedKineticEnerg 332 } 333 else // stop primary, but pass X-ray CR 334 { 335 // fParticleChange->ProposeLocalEnergyDe 336 fParticleChange->SetProposedKineticEnerg 337 } 338 // create G4DynamicParticle object for pho 339 340 auto photonRay = new G4DynamicParticle; 341 photonRay->SetDefinition( G4Gamma::Gamma() 342 photonRay->SetKineticEnergy( deltaTkin ); 343 photonRay->SetMomentumDirection(deltaDirec 344 345 vdp->push_back(photonRay); 346 } 347 return; 348 } 349 350 ////////////////////////////////////////////// 351 352 G4double G4PAIPhotModel::SampleFluctuations( 353 const G4MaterialCutsC 354 const G4DynamicPartic 355 const G4double, const 356 const G4double step, 357 { 358 // return 0.; 359 G4int coupleIndex = FindCoupleIndex(matCC); 360 if(0 > coupleIndex) { return eloss; } 361 362 SetParticle(aParticle->GetDefinition()); 363 364 // G4cout << "G4PAIPhotModel::SampleFluctuat 365 // << " Eloss(keV)= " << eloss/keV << " in 366 // << matCC->GetMaterial()->GetName() << G4e 367 368 G4double Tkin = aParticle->GetKineticE 369 G4double scaledTkin = Tkin*fRatio; 370 371 G4double loss = fModelData->SampleAlongStepP 372 scaledTkin, 373 step*fChargeSquare); 374 loss += fModelData->SampleAlongStepPlasmonTr 375 376 377 // G4cout<<" PAIPhotModel::SampleFluctuatio 378 // <<step/mm<<" mm"<<G4endl; 379 return loss; 380 381 } 382 383 ////////////////////////////////////////////// 384 // 385 // Returns the statistical estimation of the e 386 // 387 388 389 G4double G4PAIPhotModel::Dispersion(const G4Ma 390 const G4Dy 391 const G4double tcut, 392 const G4double tmax, 393 const G4double step) 394 { 395 G4double particleMass = aParticle->GetMass( 396 G4double electronDensity = material->GetElec 397 G4double kineticEnergy = aParticle->GetKinet 398 G4double q = aParticle->GetCharge()/eplus; 399 G4double etot = kineticEnergy + particleMass 400 G4double beta2 = kineticEnergy*(kineticEnerg 401 G4double siga = (tmax/beta2 - 0.5*tcut) * t 402 * electronDensity * q * q; 403 404 return siga; 405 } 406 407 ////////////////////////////////////////////// 408 409 G4double G4PAIPhotModel::MaxSecondaryEnergy( c 410 G4double kinEnergy) 411 { 412 SetParticle(p); 413 G4double tmax = kinEnergy; 414 if(p == fElectron) { tmax *= 0.5; } 415 else if(p != fPositron) { 416 G4double ratio= electron_mass_c2/fMass; 417 G4double gamma= kinEnergy/fMass + 1.0; 418 tmax = 2.0*electron_mass_c2*(gamma*gamma - 419 (1. + 2.0*gamma*ratio + rati 420 } 421 return tmax; 422 } 423 424 ////////////////////////////////////////////// 425 426 void G4PAIPhotModel::DefineForRegion(const G4R 427 { 428 fPAIRegionVector.push_back(r); 429 } 430 431 // 432 // 433 ////////////////////////////////////////////// 434