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 // ------------------------------------------------------------------- 27 // 28 // Geant4 Class file 29 // 30 // File name: G4PolarizedAnnihilationModel 31 // 32 // Author: Andreas Schaelicke 33 // 34 // Class Description: 35 // Implementation of polarized gamma Annihilation scattering on free electron 36 37 #include "G4PolarizedAnnihilationModel.hh" 38 39 #include "G4Gamma.hh" 40 #include "G4ParticleChangeForGamma.hh" 41 #include "G4PhysicalConstants.hh" 42 #include "G4PolarizationHelper.hh" 43 #include "G4PolarizationManager.hh" 44 #include "G4PolarizedAnnihilationXS.hh" 45 #include "G4StokesVector.hh" 46 #include "G4TrackStatus.hh" 47 48 G4PolarizedAnnihilationModel::G4PolarizedAnnihilationModel( 49 const G4ParticleDefinition* p, const G4String& nam) 50 : G4eeToTwoGammaModel(p, nam) 51 , fCrossSectionCalculator(nullptr) 52 , fParticleChange(nullptr) 53 , fVerboseLevel(0) 54 { 55 fCrossSectionCalculator = new G4PolarizedAnnihilationXS(); 56 fBeamPolarization = G4StokesVector::ZERO; 57 fTargetPolarization = G4StokesVector::ZERO; 58 fFinalGamma1Polarization = G4StokesVector::ZERO; 59 fFinalGamma2Polarization = G4StokesVector::ZERO; 60 } 61 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 63 G4PolarizedAnnihilationModel::~G4PolarizedAnnihilationModel() 64 { 65 delete fCrossSectionCalculator; 66 } 67 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 69 void G4PolarizedAnnihilationModel::Initialise(const G4ParticleDefinition* part, 70 const G4DataVector& dv) 71 { 72 G4eeToTwoGammaModel::Initialise(part, dv); 73 if(fParticleChange) 74 { 75 return; 76 } 77 fParticleChange = GetParticleChangeForGamma(); 78 } 79 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 81 G4double G4PolarizedAnnihilationModel::ComputeCrossSectionPerElectron( 82 G4double kinEnergy) 83 { 84 // cross section from base model 85 G4double xs = G4eeToTwoGammaModel::ComputeCrossSectionPerElectron(kinEnergy); 86 87 G4double polzz = fBeamPolarization.z() * fTargetPolarization.z(); 88 G4double poltt = fBeamPolarization.x() * fTargetPolarization.x() + 89 fBeamPolarization.y() * fTargetPolarization.y(); 90 if(polzz != 0 || poltt != 0) 91 { 92 G4double xval, lasym, tasym; 93 ComputeAsymmetriesPerElectron(kinEnergy, xval, lasym, tasym); 94 xs *= (1. + polzz * lasym + poltt * tasym); 95 } 96 97 return xs; 98 } 99 100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 101 void G4PolarizedAnnihilationModel::ComputeAsymmetriesPerElectron( 102 G4double ene, G4double& valueX, G4double& valueA, G4double& valueT) 103 { 104 // *** calculate asymmetries 105 G4double gam = 1. + ene / electron_mass_c2; 106 G4double xs0 = fCrossSectionCalculator->TotalXSection( 107 0., 1., gam, G4StokesVector::ZERO, G4StokesVector::ZERO); 108 G4double xsA = fCrossSectionCalculator->TotalXSection( 109 0., 1., gam, G4StokesVector::P3, G4StokesVector::P3); 110 G4double xsT1 = fCrossSectionCalculator->TotalXSection( 111 0., 1., gam, G4StokesVector::P1, G4StokesVector::P1); 112 G4double xsT2 = fCrossSectionCalculator->TotalXSection( 113 0., 1., gam, G4StokesVector::P2, G4StokesVector::P2); 114 G4double xsT = 0.5 * (xsT1 + xsT2); 115 116 valueX = xs0; 117 valueA = xsA / xs0 - 1.; 118 valueT = xsT / xs0 - 1.; 119 120 if((valueA < -1) || (1 < valueA)) 121 { 122 G4ExceptionDescription ed; 123 ed << " ERROR PolarizedAnnihilationPS::ComputeAsymmetries \n"; 124 ed << " something wrong in total cross section calculation (valueA)\n"; 125 ed << " LONG: " << valueX << "\t" << valueA << "\t" << valueT 126 << " energy = " << gam << G4endl; 127 G4Exception("G4PolarizedAnnihilationModel::ComputeAsymmetriesPerElectron", 128 "pol004", JustWarning, ed); 129 } 130 if((valueT < -1) || (1 < valueT)) 131 { 132 G4ExceptionDescription ed; 133 ed << " ERROR PolarizedAnnihilationPS::ComputeAsymmetries \n"; 134 ed << " something wrong in total cross section calculation (valueT)\n"; 135 ed << " TRAN: " << valueX << "\t" << valueA << "\t" << valueT 136 << " energy = " << gam << G4endl; 137 G4Exception("G4PolarizedAnnihilationModel::ComputeAsymmetriesPerElectron", 138 "pol005", JustWarning, ed); 139 } 140 } 141 142 void G4PolarizedAnnihilationModel::SampleSecondaries( 143 std::vector<G4DynamicParticle*>* fvect, const G4MaterialCutsCouple*, 144 const G4DynamicParticle* dp, G4double, G4double) 145 { 146 const G4Track* aTrack = fParticleChange->GetCurrentTrack(); 147 148 // kill primary 149 fParticleChange->SetProposedKineticEnergy(0.); 150 fParticleChange->ProposeTrackStatus(fStopAndKill); 151 152 // V.Ivanchenko add protection against zero kin energy 153 G4double PositKinEnergy = dp->GetKineticEnergy(); 154 155 if(PositKinEnergy == 0.0) 156 { 157 G4double cosTeta = 2. * G4UniformRand() - 1.; 158 G4double sinTeta = std::sqrt((1.0 - cosTeta) * (1.0 + cosTeta)); 159 G4double phi = twopi * G4UniformRand(); 160 G4ThreeVector dir(sinTeta * std::cos(phi), sinTeta * std::sin(phi), 161 cosTeta); 162 fvect->push_back( 163 new G4DynamicParticle(G4Gamma::Gamma(), dir, electron_mass_c2)); 164 fvect->push_back( 165 new G4DynamicParticle(G4Gamma::Gamma(), -dir, electron_mass_c2)); 166 return; 167 } 168 169 // *** obtain and save target and beam polarization *** 170 G4PolarizationManager* polarizationManager = 171 G4PolarizationManager::GetInstance(); 172 173 // obtain polarization of the beam 174 fBeamPolarization = G4StokesVector(aTrack->GetPolarization()); 175 176 // obtain polarization of the media 177 G4VPhysicalVolume* aPVolume = aTrack->GetVolume(); 178 G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume(); 179 const G4bool targetIsPolarized = polarizationManager->IsPolarized(aLVolume); 180 fTargetPolarization = polarizationManager->GetVolumePolarization(aLVolume); 181 182 if(fVerboseLevel >= 1) 183 { 184 G4cout << "G4PolarizedComptonModel::SampleSecondaries in " 185 << aLVolume->GetName() << G4endl; 186 } 187 188 // transfer target electron polarization in frame of positron 189 if(targetIsPolarized) 190 fTargetPolarization.rotateUz(dp->GetMomentumDirection()); 191 192 G4ParticleMomentum PositDirection = dp->GetMomentumDirection(); 193 194 // polar asymmetry: 195 G4double polarization = fBeamPolarization.p3() * fTargetPolarization.p3(); 196 197 G4double gamam1 = PositKinEnergy / electron_mass_c2; 198 G4double gama = gamam1 + 1., gamap1 = gamam1 + 2.; 199 G4double sqgrate = std::sqrt(gamam1 / gamap1) / 2., 200 sqg2m1 = std::sqrt(gamam1 * gamap1); 201 202 // limits of the energy sampling 203 G4double epsilmin = 0.5 - sqgrate, epsilmax = 0.5 + sqgrate; 204 G4double epsilqot = epsilmax / epsilmin; 205 206 // sample the energy rate of the created gammas 207 // note: for polarized partices, the actual dicing strategy 208 // will depend on the energy, and the degree of polarization !! 209 G4double epsil; 210 G4double gmax = 1. + std::fabs(polarization); // crude estimate 211 212 fCrossSectionCalculator->Initialize(epsilmin, gama, 0., fBeamPolarization, 213 fTargetPolarization); 214 if(fCrossSectionCalculator->DiceEpsilon() < 0.) 215 { 216 G4ExceptionDescription ed; 217 ed << "ERROR in PolarizedAnnihilationPS::PostStepDoIt\n" 218 << "epsilmin DiceRoutine not appropriate ! " 219 << fCrossSectionCalculator->DiceEpsilon() << G4endl; 220 G4Exception("G4PolarizedAnnihilationModel::SampleSecondaries", "pol006", 221 JustWarning, ed); 222 } 223 224 fCrossSectionCalculator->Initialize(epsilmax, gama, 0., fBeamPolarization, 225 fTargetPolarization); 226 if(fCrossSectionCalculator->DiceEpsilon() < 0) 227 { 228 G4ExceptionDescription ed; 229 ed << "ERROR in PolarizedAnnihilationPS::PostStepDoIt\n" 230 << "epsilmax DiceRoutine not appropriate ! " 231 << fCrossSectionCalculator->DiceEpsilon() << G4endl; 232 G4Exception("G4PolarizedAnnihilationModel::SampleSecondaries", "pol007", 233 JustWarning, ed); 234 } 235 236 G4int ncount = 0; 237 G4double trejectmax = 0.; 238 G4double treject; 239 240 do 241 { 242 epsil = epsilmin * std::pow(epsilqot, G4UniformRand()); 243 244 fCrossSectionCalculator->Initialize(epsil, gama, 0., fBeamPolarization, 245 fTargetPolarization, 1); 246 247 treject = fCrossSectionCalculator->DiceEpsilon(); 248 treject *= epsil; 249 250 if(treject > gmax || treject < 0.) 251 { 252 G4ExceptionDescription ed; 253 ed << "ERROR in PolarizedAnnihilationPS::PostStepDoIt\n" 254 << " eps (" << epsil 255 << ") rejection does not work properly: " << treject << G4endl; 256 G4Exception("G4PolarizedAnnihilationModel::SampleSecondaries", "pol008", 257 JustWarning, ed); 258 } 259 ++ncount; 260 if(treject > trejectmax) 261 trejectmax = treject; 262 if(ncount > 1000) 263 { 264 G4ExceptionDescription ed; 265 ed << "WARNING in PolarizedAnnihilationPS::PostStepDoIt\n" 266 << "eps dicing very inefficient =" << trejectmax / gmax << ", " 267 << treject / gmax << ". For secondary energy = " << epsil << " " 268 << ncount << G4endl; 269 G4Exception("G4PolarizedAnnihilationModel::SampleSecondaries", "pol009", 270 JustWarning, ed); 271 break; 272 } 273 274 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko 275 } while(treject < gmax * G4UniformRand()); 276 277 // scattered Gamma angles. ( Z - axis along the parent positron) 278 G4double cost = (epsil * gamap1 - 1.) / (epsil * sqg2m1); 279 G4double sint = std::sqrt((1. + cost) * (1. - cost)); 280 G4double phi = 0.; 281 G4double beamTrans = 282 std::sqrt(sqr(fBeamPolarization.p1()) + sqr(fBeamPolarization.p2())); 283 G4double targetTrans = 284 std::sqrt(sqr(fTargetPolarization.p1()) + sqr(fTargetPolarization.p2())); 285 286 do 287 { 288 phi = twopi * G4UniformRand(); 289 fCrossSectionCalculator->Initialize(epsil, gama, 0., fBeamPolarization, 290 fTargetPolarization, 2); 291 292 G4double gdiced = fCrossSectionCalculator->getVar(0); 293 gdiced += fCrossSectionCalculator->getVar(3) * fBeamPolarization.p3() * 294 fTargetPolarization.p3(); 295 gdiced += 1. * 296 (std::fabs(fCrossSectionCalculator->getVar(1)) + 297 std::fabs(fCrossSectionCalculator->getVar(2))) * 298 beamTrans * targetTrans; 299 gdiced += 1. * std::fabs(fCrossSectionCalculator->getVar(4)) * 300 (std::fabs(fBeamPolarization.p3()) * targetTrans + 301 std::fabs(fTargetPolarization.p3()) * beamTrans); 302 303 G4double gdist = fCrossSectionCalculator->getVar(0); 304 gdist += fCrossSectionCalculator->getVar(3) * fBeamPolarization.p3() * 305 fTargetPolarization.p3(); 306 gdist += fCrossSectionCalculator->getVar(1) * 307 (std::cos(phi) * fBeamPolarization.p1() + 308 std::sin(phi) * fBeamPolarization.p2()) * 309 (std::cos(phi) * fTargetPolarization.p1() + 310 std::sin(phi) * fTargetPolarization.p2()); 311 gdist += fCrossSectionCalculator->getVar(2) * 312 (std::cos(phi) * fBeamPolarization.p2() - 313 std::sin(phi) * fBeamPolarization.p1()) * 314 (std::cos(phi) * fTargetPolarization.p2() - 315 std::sin(phi) * fTargetPolarization.p1()); 316 gdist += 317 fCrossSectionCalculator->getVar(4) * 318 (std::cos(phi) * fBeamPolarization.p3() * fTargetPolarization.p1() + 319 std::cos(phi) * fBeamPolarization.p1() * fTargetPolarization.p3() + 320 std::sin(phi) * fBeamPolarization.p3() * fTargetPolarization.p2() + 321 std::sin(phi) * fBeamPolarization.p2() * fTargetPolarization.p3()); 322 323 treject = gdist / gdiced; 324 if(treject > 1. + 1.e-10 || treject < 0) 325 { 326 G4ExceptionDescription ed; 327 ed << "!!!ERROR in PolarizedAnnihilationPS::PostStepDoIt\n" 328 << " phi rejection does not work properly: " << treject << G4endl; 329 G4cout << " gdiced = " << gdiced << G4endl; 330 G4cout << " gdist = " << gdist << G4endl; 331 G4cout << " epsil = " << epsil << G4endl; 332 G4Exception("G4PolarizedAnnihilationModel::SampleSecondaries", "pol009", 333 JustWarning, ed); 334 } 335 336 if(treject < 1.e-3) 337 { 338 G4ExceptionDescription ed; 339 ed << "!!!ERROR in PolarizedAnnihilationPS::PostStepDoIt\n" 340 << " phi rejection does not work properly: " << treject << "\n"; 341 G4cout << " gdiced=" << gdiced << " gdist=" << gdist << "\n"; 342 G4cout << " epsil = " << epsil << G4endl; 343 G4Exception("G4PolarizedAnnihilationModel::SampleSecondaries", "pol010", 344 JustWarning, ed); 345 } 346 347 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko 348 } while(treject < G4UniformRand()); 349 350 G4double dirx = sint * std::cos(phi); 351 G4double diry = sint * std::sin(phi); 352 G4double dirz = cost; 353 354 // kinematic of the created pair 355 G4double TotalAvailableEnergy = PositKinEnergy + 2 * electron_mass_c2; 356 G4double Phot1Energy = epsil * TotalAvailableEnergy; 357 G4double Phot2Energy = (1. - epsil) * TotalAvailableEnergy; 358 359 // *** prepare calculation of polarization transfer *** 360 G4ThreeVector Phot1Direction(dirx, diry, dirz); 361 362 // get interaction frame 363 G4ThreeVector nInteractionFrame = 364 G4PolarizationHelper::GetFrame(PositDirection, Phot1Direction); 365 366 // define proper in-plane and out-of-plane component of initial spins 367 fBeamPolarization.InvRotateAz(nInteractionFrame, PositDirection); 368 fTargetPolarization.InvRotateAz(nInteractionFrame, PositDirection); 369 370 // calculate spin transfere matrix 371 372 fCrossSectionCalculator->Initialize(epsil, gama, phi, fBeamPolarization, 373 fTargetPolarization, 2); 374 375 Phot1Direction.rotateUz(PositDirection); 376 // create G4DynamicParticle object for the particle1 377 G4DynamicParticle* aParticle1 = 378 new G4DynamicParticle(G4Gamma::Gamma(), Phot1Direction, Phot1Energy); 379 fFinalGamma1Polarization = fCrossSectionCalculator->GetPol2(); 380 G4double n1 = fFinalGamma1Polarization.mag2(); 381 if(n1 > 1.) 382 { 383 G4ExceptionDescription ed; 384 ed << "ERROR: PolarizedAnnihilation Polarization Vector at epsil = " 385 << epsil << " is too large!!! \n" 386 << "annihi pol1= " << fFinalGamma1Polarization << ", (" << n1 << ")\n"; 387 fFinalGamma1Polarization *= 1. / std::sqrt(n1); 388 G4Exception("G4PolarizedAnnihilationModel::SampleSecondaries", "pol011", 389 JustWarning, ed); 390 } 391 392 // define polarization of first final state photon 393 fFinalGamma1Polarization.SetPhoton(); 394 fFinalGamma1Polarization.RotateAz(nInteractionFrame, Phot1Direction); 395 aParticle1->SetPolarization(fFinalGamma1Polarization.p1(), 396 fFinalGamma1Polarization.p2(), 397 fFinalGamma1Polarization.p3()); 398 399 fvect->push_back(aParticle1); 400 401 // ********************************************************************** 402 403 G4double Eratio = Phot1Energy / Phot2Energy; 404 G4double PositP = 405 std::sqrt(PositKinEnergy * (PositKinEnergy + 2. * electron_mass_c2)); 406 G4ThreeVector Phot2Direction(-dirx * Eratio, -diry * Eratio, 407 (PositP - dirz * Phot1Energy) / Phot2Energy); 408 Phot2Direction.rotateUz(PositDirection); 409 // create G4DynamicParticle object for the particle2 410 G4DynamicParticle* aParticle2 = 411 new G4DynamicParticle(G4Gamma::Gamma(), Phot2Direction, Phot2Energy); 412 413 // define polarization of second final state photon 414 fFinalGamma2Polarization = fCrossSectionCalculator->GetPol3(); 415 G4double n2 = fFinalGamma2Polarization.mag2(); 416 if(n2 > 1.) 417 { 418 G4ExceptionDescription ed; 419 ed << "ERROR: PolarizedAnnihilation Polarization Vector at epsil = " 420 << epsil << " is too large!!! \n"; 421 ed << "annihi pol2= " << fFinalGamma2Polarization << ", (" << n2 << ")\n"; 422 423 G4Exception("G4PolarizedAnnihilationModel::SampleSecondaries", "pol012", 424 JustWarning, ed); 425 fFinalGamma2Polarization *= 1. / std::sqrt(n2); 426 } 427 fFinalGamma2Polarization.SetPhoton(); 428 fFinalGamma2Polarization.RotateAz(nInteractionFrame, Phot2Direction); 429 aParticle2->SetPolarization(fFinalGamma2Polarization.p1(), 430 fFinalGamma2Polarization.p2(), 431 fFinalGamma2Polarization.p3()); 432 433 fvect->push_back(aParticle2); 434 } 435