[819] | 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 | // $Id: G4PolarizedComptonModel.cc,v 1.4 2007/05/23 08:52:20 vnivanch Exp $ |
---|
[1228] | 28 | // GEANT4 tag $Name: geant4-09-03 $ |
---|
[819] | 29 | // |
---|
| 30 | // ------------------------------------------------------------------- |
---|
| 31 | // |
---|
| 32 | // GEANT4 Class file |
---|
| 33 | // |
---|
| 34 | // |
---|
| 35 | // File name: G4PolarizedComptonModel |
---|
| 36 | // |
---|
| 37 | // Author: Andreas Schaelicke |
---|
| 38 | // |
---|
| 39 | // Creation date: 01.05.2005 |
---|
| 40 | // |
---|
| 41 | // Modifications: |
---|
| 42 | // 18-07-06 use newly calculated cross sections (P. Starovoitov) |
---|
| 43 | // 21-08-05 update interface (A. Schaelicke) |
---|
| 44 | // |
---|
| 45 | // Class Description: |
---|
| 46 | // |
---|
| 47 | // ------------------------------------------------------------------- |
---|
| 48 | // |
---|
| 49 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 50 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 51 | |
---|
| 52 | #include "G4PolarizedComptonModel.hh" |
---|
| 53 | #include "G4Electron.hh" |
---|
| 54 | #include "G4Gamma.hh" |
---|
| 55 | #include "Randomize.hh" |
---|
| 56 | #include "G4DataVector.hh" |
---|
| 57 | #include "G4ParticleChangeForGamma.hh" |
---|
| 58 | |
---|
| 59 | |
---|
| 60 | #include "G4StokesVector.hh" |
---|
| 61 | #include "G4PolarizationManager.hh" |
---|
| 62 | #include "G4PolarizationHelper.hh" |
---|
| 63 | #include "G4PolarizedComptonCrossSection.hh" |
---|
| 64 | |
---|
| 65 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 66 | |
---|
| 67 | G4PolarizedComptonModel::G4PolarizedComptonModel(const G4ParticleDefinition*, |
---|
| 68 | const G4String& nam) |
---|
| 69 | : G4KleinNishinaCompton(0,nam), |
---|
| 70 | verboseLevel(0) |
---|
| 71 | { |
---|
| 72 | crossSectionCalculator=new G4PolarizedComptonCrossSection(); |
---|
| 73 | } |
---|
| 74 | |
---|
| 75 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 76 | |
---|
| 77 | G4PolarizedComptonModel::~G4PolarizedComptonModel() |
---|
| 78 | { |
---|
| 79 | if (crossSectionCalculator) delete crossSectionCalculator; |
---|
| 80 | } |
---|
| 81 | |
---|
| 82 | |
---|
| 83 | |
---|
| 84 | G4double G4PolarizedComptonModel::ComputeAsymmetryPerAtom |
---|
| 85 | (G4double gammaEnergy, G4double /*Z*/) |
---|
| 86 | |
---|
| 87 | { |
---|
| 88 | G4double asymmetry = 0.0 ; |
---|
| 89 | |
---|
| 90 | G4double k0 = gammaEnergy / electron_mass_c2 ; |
---|
| 91 | G4double k1 = 1 + 2*k0 ; |
---|
| 92 | |
---|
| 93 | asymmetry = -k0; |
---|
| 94 | asymmetry *= (k0 + 1.)*sqr(k1)*std::log(k1) - 2.*k0*(5.*sqr(k0) + 4.*k0 + 1.); |
---|
| 95 | asymmetry /= ((k0 - 2.)*k0 -2.)*sqr(k1)*std::log(k1) + 2.*k0*(k0*(k0 + 1.)*(k0 + 8.) + 2.); |
---|
| 96 | |
---|
| 97 | // G4cout<<"energy = "<<GammaEnergy<<" asymmetry = "<<asymmetry<<"\t\t GAM = "<<k0<<G4endl; |
---|
| 98 | if (asymmetry>1.) G4cout<<"ERROR in G4PolarizedComptonModel::ComputeAsymmetryPerAtom"<<G4endl; |
---|
| 99 | |
---|
| 100 | return asymmetry; |
---|
| 101 | } |
---|
| 102 | |
---|
| 103 | |
---|
| 104 | G4double G4PolarizedComptonModel::ComputeCrossSectionPerAtom( |
---|
| 105 | const G4ParticleDefinition* pd, |
---|
| 106 | G4double kinEnergy, |
---|
| 107 | G4double Z, |
---|
| 108 | G4double A, |
---|
| 109 | G4double cut, |
---|
| 110 | G4double emax) |
---|
| 111 | { |
---|
| 112 | double xs = |
---|
| 113 | G4KleinNishinaCompton::ComputeCrossSectionPerAtom(pd,kinEnergy, |
---|
| 114 | Z,A,cut,emax); |
---|
| 115 | G4double polzz = theBeamPolarization.p3()*theTargetPolarization.z(); |
---|
| 116 | if (polzz!=0) { |
---|
| 117 | G4double asym=ComputeAsymmetryPerAtom(kinEnergy, Z); |
---|
| 118 | xs*=(1.+polzz*asym); |
---|
| 119 | } |
---|
| 120 | return xs; |
---|
| 121 | } |
---|
| 122 | |
---|
| 123 | |
---|
| 124 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 125 | |
---|
| 126 | void G4PolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, |
---|
| 127 | const G4MaterialCutsCouple*, |
---|
| 128 | const G4DynamicParticle* aDynamicGamma, |
---|
| 129 | G4double, |
---|
| 130 | G4double) |
---|
| 131 | { |
---|
| 132 | const G4Track * aTrack = fParticleChange->GetCurrentTrack(); |
---|
| 133 | G4VPhysicalVolume* aPVolume = aTrack->GetVolume(); |
---|
| 134 | G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume(); |
---|
| 135 | |
---|
| 136 | if (verboseLevel>=1) |
---|
| 137 | G4cout<<"G4PolarizedComptonModel::SampleSecondaries in " |
---|
| 138 | << aLVolume->GetName() <<G4endl; |
---|
| 139 | |
---|
| 140 | G4PolarizationManager * polarizationManager = G4PolarizationManager::GetInstance(); |
---|
| 141 | |
---|
| 142 | // obtain polarization of the beam |
---|
| 143 | theBeamPolarization = aDynamicGamma->GetPolarization(); |
---|
| 144 | theBeamPolarization.SetPhoton(); |
---|
| 145 | |
---|
| 146 | // obtain polarization of the media |
---|
| 147 | const G4bool targetIsPolarized = polarizationManager->IsPolarized(aLVolume); |
---|
| 148 | theTargetPolarization = polarizationManager->GetVolumePolarization(aLVolume); |
---|
| 149 | |
---|
| 150 | // if beam is linear polarized or target is transversely polarized |
---|
| 151 | // determine the angle to x-axis |
---|
| 152 | // (assumes same PRF as in the polarization definition) |
---|
| 153 | |
---|
| 154 | G4ThreeVector gamDirection0 = aDynamicGamma->GetMomentumDirection(); |
---|
| 155 | |
---|
| 156 | // transfere theTargetPolarization |
---|
| 157 | // into the gamma frame (problem electron is at rest) |
---|
| 158 | if (targetIsPolarized) |
---|
| 159 | theTargetPolarization.rotateUz(gamDirection0); |
---|
| 160 | |
---|
| 161 | // The scattered gamma energy is sampled according to Klein - Nishina formula. |
---|
| 162 | // The random number techniques of Butcher & Messel are used |
---|
| 163 | // (Nuc Phys 20(1960),15). |
---|
| 164 | // Note : Effects due to binding of atomic electrons are negliged. |
---|
| 165 | |
---|
| 166 | G4double gamEnergy0 = aDynamicGamma->GetKineticEnergy(); |
---|
| 167 | G4double E0_m = gamEnergy0 / electron_mass_c2 ; |
---|
| 168 | |
---|
| 169 | |
---|
| 170 | // |
---|
| 171 | // sample the energy rate of the scattered gamma |
---|
| 172 | // |
---|
| 173 | |
---|
| 174 | G4double epsilon, epsilonsq, onecost, sint2, greject ; |
---|
| 175 | |
---|
| 176 | G4double epsilon0 = 1./(1. + 2.*E0_m); |
---|
| 177 | G4double epsilon0sq = epsilon0*epsilon0; |
---|
| 178 | G4double alpha1 = - std::log(epsilon0); |
---|
| 179 | G4double alpha2 = 0.5*(1.- epsilon0sq); |
---|
| 180 | |
---|
| 181 | G4double polarization = theBeamPolarization.p3()*theTargetPolarization.p3(); |
---|
| 182 | do { |
---|
| 183 | if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) { |
---|
| 184 | epsilon = std::exp(-alpha1*G4UniformRand()); // epsilon0**r |
---|
| 185 | epsilonsq = epsilon*epsilon; |
---|
| 186 | |
---|
| 187 | } else { |
---|
| 188 | epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand(); |
---|
| 189 | epsilon = std::sqrt(epsilonsq); |
---|
| 190 | }; |
---|
| 191 | |
---|
| 192 | onecost = (1.- epsilon)/(epsilon*E0_m); |
---|
| 193 | sint2 = onecost*(2.-onecost); |
---|
| 194 | |
---|
| 195 | |
---|
| 196 | G4double gdiced = 2.*(1./epsilon+epsilon); |
---|
| 197 | G4double gdist = 1./epsilon + epsilon - sint2 |
---|
| 198 | - polarization*(1./epsilon-epsilon)*(1.-onecost); |
---|
| 199 | |
---|
| 200 | greject = gdist/gdiced; |
---|
| 201 | |
---|
| 202 | if (greject>1) G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n" |
---|
| 203 | <<" costh rejection does not work properly: "<<greject<<G4endl; |
---|
| 204 | |
---|
| 205 | } while (greject < G4UniformRand()); |
---|
| 206 | |
---|
| 207 | // |
---|
| 208 | // scattered gamma angles. ( Z - axis along the parent gamma) |
---|
| 209 | // |
---|
| 210 | |
---|
| 211 | G4double cosTeta = 1. - onecost; |
---|
| 212 | G4double sinTeta = std::sqrt (sint2); |
---|
| 213 | G4double Phi; |
---|
| 214 | do { |
---|
| 215 | Phi = twopi * G4UniformRand(); |
---|
| 216 | G4double gdiced = 1./epsilon + epsilon - sint2 |
---|
| 217 | + std::abs(theBeamPolarization.p3())* |
---|
| 218 | ( std::abs((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3()) |
---|
| 219 | +(1.-epsilon)*sinTeta*(std::sqrt(sqr(theTargetPolarization.p1()) |
---|
| 220 | + sqr(theTargetPolarization.p2())))) |
---|
| 221 | +sint2*(std::sqrt(sqr(theBeamPolarization.p1()) + sqr(theBeamPolarization.p2()))); |
---|
| 222 | |
---|
| 223 | G4double gdist = 1./epsilon + epsilon - sint2 |
---|
| 224 | + theBeamPolarization.p3()* |
---|
| 225 | ((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3() |
---|
| 226 | +(1.-epsilon)*sinTeta*(std::cos(Phi)*theTargetPolarization.p1()+ |
---|
| 227 | std::sin(Phi)*theTargetPolarization.p2())) |
---|
| 228 | -sint2*(std::cos(2.*Phi)*theBeamPolarization.p1() |
---|
| 229 | +std::sin(2.*Phi)*theBeamPolarization.p2()); |
---|
| 230 | greject = gdist/gdiced; |
---|
| 231 | |
---|
| 232 | if (greject>1.+1.e-10 || greject<0) G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n" |
---|
| 233 | <<" phi rejection does not work properly: "<<greject<<G4endl; |
---|
| 234 | |
---|
| 235 | if (greject<1.e-3) { |
---|
| 236 | G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n" |
---|
| 237 | <<" phi rejection does not work properly: "<<greject<<"\n"; |
---|
| 238 | G4cout<<" greject="<<greject<<" phi="<<Phi<<" cost="<<cosTeta<<"\n"; |
---|
| 239 | G4cout<<" gdiced="<<gdiced<<" gdist="<<gdist<<"\n"; |
---|
| 240 | G4cout<<" eps="<<epsilon<<" 1/eps="<<1./epsilon<<"\n"; |
---|
| 241 | } |
---|
| 242 | |
---|
| 243 | } while (greject < G4UniformRand()); |
---|
| 244 | G4double dirx = sinTeta*std::cos(Phi), diry = sinTeta*std::sin(Phi), dirz = cosTeta; |
---|
| 245 | |
---|
| 246 | // |
---|
| 247 | // update G4VParticleChange for the scattered gamma |
---|
| 248 | // |
---|
| 249 | |
---|
| 250 | G4ThreeVector gamDirection1 ( dirx,diry,dirz ); |
---|
| 251 | gamDirection1.rotateUz(gamDirection0); |
---|
| 252 | G4double gamEnergy1 = epsilon*gamEnergy0; |
---|
| 253 | fParticleChange->SetProposedKineticEnergy(gamEnergy1); |
---|
| 254 | |
---|
| 255 | |
---|
| 256 | if(gamEnergy1 > lowestGammaEnergy) { |
---|
| 257 | fParticleChange->ProposeMomentumDirection(gamDirection1); |
---|
| 258 | } else { |
---|
| 259 | fParticleChange->ProposeTrackStatus(fStopAndKill); |
---|
| 260 | gamEnergy1 += fParticleChange->GetLocalEnergyDeposit(); |
---|
| 261 | fParticleChange->ProposeLocalEnergyDeposit(gamEnergy1); |
---|
| 262 | } |
---|
| 263 | |
---|
| 264 | // |
---|
| 265 | // kinematic of the scattered electron |
---|
| 266 | // |
---|
| 267 | |
---|
| 268 | G4double eKinEnergy = gamEnergy0 - gamEnergy1; |
---|
| 269 | G4ThreeVector eDirection = gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1; |
---|
| 270 | eDirection = eDirection.unit(); |
---|
| 271 | |
---|
| 272 | // |
---|
| 273 | // calculate Stokesvector of final state photon and electron |
---|
| 274 | // |
---|
| 275 | G4ThreeVector nInteractionFrame; |
---|
| 276 | if((gamEnergy1 > lowestGammaEnergy) || |
---|
| 277 | (eKinEnergy > DBL_MIN)) { |
---|
| 278 | |
---|
| 279 | // determine interaction plane |
---|
| 280 | // nInteractionFrame = |
---|
| 281 | // G4PolarizationHelper::GetFrame(gamDirection1,eDirection); |
---|
| 282 | if (gamEnergy1 > lowestGammaEnergy) |
---|
| 283 | nInteractionFrame = G4PolarizationHelper::GetFrame(gamDirection1,gamDirection0); |
---|
| 284 | else |
---|
| 285 | nInteractionFrame = G4PolarizationHelper::GetFrame(gamDirection0, eDirection); |
---|
| 286 | |
---|
| 287 | // transfere theBeamPolarization and theTargetPolarization |
---|
| 288 | // into the interaction frame (note electron is in gamma frame) |
---|
| 289 | if (verboseLevel>=1) { |
---|
| 290 | G4cout << "========================================\n"; |
---|
| 291 | G4cout << " nInteractionFrame = " <<nInteractionFrame<<"\n"; |
---|
| 292 | G4cout << " GammaDirection0 = " <<gamDirection0<<"\n"; |
---|
| 293 | G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n"; |
---|
| 294 | G4cout << " electronPolarization = " <<theTargetPolarization<<"\n"; |
---|
| 295 | } |
---|
| 296 | |
---|
| 297 | theBeamPolarization.InvRotateAz(nInteractionFrame,gamDirection0); |
---|
| 298 | theTargetPolarization.InvRotateAz(nInteractionFrame,gamDirection0); |
---|
| 299 | |
---|
| 300 | if (verboseLevel>=1) { |
---|
| 301 | G4cout << "----------------------------------------\n"; |
---|
| 302 | G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n"; |
---|
| 303 | G4cout << " electronPolarization = " <<theTargetPolarization<<"\n"; |
---|
| 304 | G4cout << "----------------------------------------\n"; |
---|
| 305 | } |
---|
| 306 | |
---|
| 307 | // initialize the polarization transfer matrix |
---|
| 308 | crossSectionCalculator->Initialize(epsilon,E0_m,0., |
---|
| 309 | theBeamPolarization, |
---|
| 310 | theTargetPolarization,2); |
---|
| 311 | } |
---|
| 312 | |
---|
| 313 | // if(eKinEnergy > DBL_MIN) |
---|
| 314 | { |
---|
| 315 | // in interaction frame |
---|
| 316 | // calculate polarization transfer to the photon (in interaction plane) |
---|
| 317 | finalGammaPolarization = crossSectionCalculator->GetPol2(); |
---|
| 318 | if (verboseLevel>=1) G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n"; |
---|
| 319 | finalGammaPolarization.SetPhoton(); |
---|
| 320 | |
---|
| 321 | // translate polarization into particle reference frame |
---|
| 322 | finalGammaPolarization.RotateAz(nInteractionFrame,gamDirection1); |
---|
| 323 | //store polarization vector |
---|
| 324 | fParticleChange->ProposePolarization(finalGammaPolarization); |
---|
| 325 | if (finalGammaPolarization.mag() > 1.+1.e-8){ |
---|
| 326 | G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl; |
---|
| 327 | G4cout<<"Polarization of final photon more than 100%"<<G4endl; |
---|
| 328 | G4cout<<finalGammaPolarization<<" mag = "<<finalGammaPolarization.mag()<<G4endl; |
---|
| 329 | } |
---|
| 330 | if (verboseLevel>=1) { |
---|
| 331 | G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n"; |
---|
| 332 | G4cout << " GammaDirection1 = " <<gamDirection1<<"\n"; |
---|
| 333 | } |
---|
| 334 | } |
---|
| 335 | |
---|
| 336 | // if (ElecKineEnergy > fminimalEnergy) { |
---|
| 337 | { |
---|
| 338 | finalElectronPolarization = crossSectionCalculator->GetPol3(); |
---|
| 339 | if (verboseLevel>=1) |
---|
| 340 | G4cout << " electronPolarization1 = " <<finalElectronPolarization<<"\n"; |
---|
| 341 | |
---|
| 342 | // transfer into particle reference frame |
---|
| 343 | finalElectronPolarization.RotateAz(nInteractionFrame,eDirection); |
---|
| 344 | if (verboseLevel>=1) { |
---|
| 345 | G4cout << " electronPolarization1 = " <<finalElectronPolarization<<"\n"; |
---|
| 346 | G4cout << " ElecDirection = " <<eDirection<<"\n"; |
---|
| 347 | } |
---|
| 348 | } |
---|
| 349 | if (verboseLevel>=1) |
---|
| 350 | G4cout << "========================================\n"; |
---|
| 351 | |
---|
| 352 | |
---|
| 353 | if(eKinEnergy > DBL_MIN) { |
---|
| 354 | |
---|
| 355 | // create G4DynamicParticle object for the electron. |
---|
| 356 | G4DynamicParticle* aElectron = new G4DynamicParticle(theElectron,eDirection,eKinEnergy); |
---|
| 357 | //store polarization vector |
---|
| 358 | if (finalElectronPolarization.mag() > 1.+1.e-8){ |
---|
| 359 | G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl; |
---|
| 360 | G4cout<<"Polarization of final electron more than 100%"<<G4endl; |
---|
| 361 | G4cout<<finalElectronPolarization<<" mag = "<<finalElectronPolarization.mag()<<G4endl; |
---|
| 362 | } |
---|
| 363 | aElectron->SetPolarization(finalElectronPolarization.p1(), |
---|
| 364 | finalElectronPolarization.p2(), |
---|
| 365 | finalElectronPolarization.p3()); |
---|
| 366 | fvect->push_back(aElectron); |
---|
| 367 | } |
---|
| 368 | } |
---|
| 369 | |
---|
| 370 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... |
---|
| 371 | |
---|
| 372 | |
---|