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 // Geant4 Class file 29 // 30 // File name: G4PolarizedIonisationModel 31 // 32 // Author: A.Schaelicke on base of Vlad 33 // 34 // Class Description: 35 // Implementation of energy loss and delta-e 36 // (including polarization effects) 37 38 #include "G4PolarizedIonisationModel.hh" 39 40 #include "G4ParticleChangeForLoss.hh" 41 #include "G4PhysicalConstants.hh" 42 #include "G4PolarizationHelper.hh" 43 #include "G4PolarizationManager.hh" 44 #include "G4PolarizedIonisationBhabhaXS.hh" 45 #include "G4PolarizedIonisationMollerXS.hh" 46 #include "Randomize.hh" 47 48 //....oooOO0OOooo........oooOO0OOooo........oo 49 G4PolarizedIonisationModel::G4PolarizedIonisat 50 const G4ParticleDefinition* p, const G4Strin 51 : G4MollerBhabhaModel(p, nam) 52 , fCrossSectionCalculator(nullptr) 53 { 54 fBeamPolarization = G4StokesVector::ZERO; 55 fTargetPolarization = G4StokesVector::ZERO; 56 57 fPositronPolarization = G4StokesVector::ZERO 58 fElectronPolarization = G4StokesVector::ZERO 59 60 isElectron = (p == theElectron); // necessa 61 // G4Molle 62 63 if(!isElectron) 64 { 65 G4cout << " buildBhabha cross section " << 66 fCrossSectionCalculator = new G4PolarizedI 67 } 68 else 69 { 70 G4cout << " buildMoller cross section " << 71 fCrossSectionCalculator = new G4PolarizedI 72 } 73 } 74 75 //....oooOO0OOooo........oooOO0OOooo........oo 76 G4PolarizedIonisationModel::~G4PolarizedIonisa 77 { 78 if(fCrossSectionCalculator) 79 { 80 delete fCrossSectionCalculator; 81 } 82 } 83 84 //....oooOO0OOooo........oooOO0OOooo........oo 85 G4double G4PolarizedIonisationModel::ComputeCr 86 const G4ParticleDefinition* pd, G4double kin 87 G4double emax) 88 { 89 G4double xs = G4MollerBhabhaModel::ComputeCr 90 pd, kinEnergy, cut, emax); 91 G4double factor = 1.; 92 if(xs != 0.) 93 { 94 G4double tmax = MaxSecondaryEnergy(pd, kin 95 tmax = std::min(emax, tmax); 96 97 if(std::fabs(cut / emax - 1.) < 1.e-10) 98 return xs; 99 100 if(cut < tmax) 101 { 102 G4double xmin = cut / kinEnergy; 103 G4double xmax = tmax / kinEnergy; 104 G4double gam = kinEnergy / electron 105 G4double crossPol = fCrossSectionCalcula 106 xmin, xmax, gam, fBeamPolarization, fT 107 G4double crossUnpol = fCrossSectionCalcu 108 xmin, xmax, gam, G4StokesVector::ZERO, 109 if(crossUnpol > 0.) 110 factor = crossPol / crossUnpol; 111 } 112 } 113 return xs * factor; 114 } 115 116 //....oooOO0OOooo........oooOO0OOooo........oo 117 void G4PolarizedIonisationModel::SampleSeconda 118 std::vector<G4DynamicParticle*>* vdp, const 119 const G4DynamicParticle* dp, G4double tmin, 120 { 121 // *** obtain and save target and beam polar 122 G4PolarizationManager* polarizationManager = 123 G4PolarizationManager::GetInstance(); 124 125 const G4Track* aTrack = fParticleChange->Get 126 127 // obtain polarization of the beam 128 fBeamPolarization = G4StokesVector(dp->GetPo 129 130 // obtain polarization of the media 131 G4VPhysicalVolume* aPVolume = aTrack->Get 132 G4LogicalVolume* aLVolume = aPVolume->G 133 const G4bool targetIsPolarized = polarizatio 134 fTargetPolarization = polarizationManager->G 135 136 // transfer target polarization in interacti 137 if(targetIsPolarized) 138 fTargetPolarization.rotateUz(dp->GetMoment 139 140 G4double tmax = std::min(maxEnergy, MaxSecon 141 if(tmin >= tmax) 142 return; 143 144 G4double polL = fBeamPolarization.z() * fTar 145 polL = std::fabs(polL); 146 G4double polT = fBeamPolarization.x() * fTar 147 fBeamPolarization.y() * fTar 148 polT = std::fabs(polT); 149 150 G4double kineticEnergy = dp->GetKineticEnerg 151 G4double energy = kineticEnergy + ele 152 G4double totalMomentum = 153 std::sqrt(kineticEnergy * (energy + electr 154 G4double xmin = tmin / kineticEnergy; 155 G4double xmax = tmax / kineticEnergy; 156 G4double gam = energy / electron_mass_c2; 157 G4double gamma2 = gam * gam; 158 G4double gmo = gam - 1.; 159 G4double gmo2 = gmo * gmo; 160 G4double gmo3 = gmo2 * gmo; 161 G4double gpo = gam + 1.; 162 G4double gpo2 = gpo * gpo; 163 G4double gpo3 = gpo2 * gpo; 164 G4double x, y, q, grej, grej2; 165 G4double z = 0.; 166 G4double xs = 0., phi = 0.; 167 G4ThreeVector direction = dp->GetMomentumDir 168 169 //(Polarized) Moller (e-e-) scattering 170 if(isElectron) 171 { 172 // *** dice according to polarized cross s 173 G4double G = ((2.0 * gam - 1.0) / gamma2) 174 G4double H = (sqr(gam - 1.0) / gamma2) * 175 (1. + polT + polL * ((gam + 3 176 177 y = 1.0 - xmax; 178 grej = 1.0 - G * xmax + xmax * xmax * (H 179 grej2 = 1.0 - G * xmin + xmin * xmin * (H 180 if(grej2 > grej) 181 grej = grej2; 182 G4double prefM = gamma2 * classic_electr_r 183 (gmo2 * (gam + 1.0)); 184 grej *= prefM; 185 do 186 { 187 q = G4UniformRand(); 188 x = xmin * xmax / (xmin * (1.0 - q) + xm 189 if(fCrossSectionCalculator) 190 { 191 fCrossSectionCalculator->Initialize(x, 192 fT 193 xs = fCrossSectionCalculator->XSection 194 195 z = xs * sqr(x) * 4.; 196 if(grej < z) 197 { 198 G4ExceptionDescription ed; 199 ed << "WARNING : error in Moller rej 200 << " z = " << z << " grej=" << gr 201 G4Exception("G4PolarizedIonisationMo 202 JustWarning, ed); 203 } 204 } 205 else 206 { 207 G4cout << "No calculator in Moller sca 208 } 209 // Loop checking, 03-Aug-2015, Vladimir 210 } while(grej * G4UniformRand() > z); 211 // Bhabha (e+e-) scattering 212 } 213 else 214 { 215 // *** dice according to polarized cross s 216 y = xmax * xmax; 217 grej = 0.; 218 grej += y * y * gmo3 * (1. + (polL + polT) 219 grej += -2. * xmin * xmin * xmin * gam * g 220 (1. - (polL + polT) * (gam + 3.) / 221 grej += y * y * gmo * (3. * gamma2 + 6. * 222 (1. + (polL * (3. * gam + 1.) * (g 223 polT * ((gam + 2.) * gamma2 224 (gmo * (3. * gam * (gam + 225 grej /= gpo3; 226 grej += -xmin * (2. * gamma2 + 4. * gam + 227 (1. - gam * (polL * (2. * gam + 1. 228 (2. * gam * (gam + 2.) + 1 229 gpo2; 230 grej += gamma2 / (gamma2 - 1.); 231 G4double prefB = 232 classic_electr_radius * classic_electr_r 233 grej *= prefB; 234 235 do 236 { 237 q = G4UniformRand(); 238 x = xmin * xmax / (xmin * (1.0 - q) + xm 239 if(fCrossSectionCalculator) 240 { 241 fCrossSectionCalculator->Initialize(x, 242 fT 243 xs = fCrossSectionCalculator->XSection 244 245 z = xs * sqr(x) * 4.; 246 } 247 else 248 { 249 G4cout << "No calculator in Bhabha sca 250 } 251 252 if(z > grej) 253 { 254 G4ExceptionDescription ed; 255 ed << "G4PolarizedIonisationModel::Sam 256 << "Majorant " << grej << " < " << 257 << " e+e- (Bhabha) scattering" 258 << " at KinEnergy " << kineticEnerg 259 G4Exception("G4PolarizedIonisationMode 260 JustWarning, ed); 261 } 262 // Loop checking, 03-Aug-2015, Vladimir 263 } while(grej * G4UniformRand() > z); 264 } 265 266 // polar asymmetries (due to transverse pola 267 if(fCrossSectionCalculator) 268 { 269 grej = xs * 2.; 270 do 271 { 272 phi = twopi * G4UniformRand(); 273 fCrossSectionCalculator->Initialize(x, g 274 fTar 275 xs = fCrossSectionCalculator->XSection(G 276 G 277 if(xs > grej) 278 { 279 if(isElectron) 280 { 281 G4ExceptionDescription ed; 282 ed << "Majorant " << grej << " < " < 283 << "\n" 284 << " e-e- (Moller) scattering\n" 285 << "PHI DICING\n"; 286 G4Exception("G4PolarizedIonisationMo 287 JustWarning, ed); 288 } 289 else 290 { 291 G4ExceptionDescription ed; 292 ed << "Majorant " << grej << " < " < 293 << "\n" 294 << " e+e- (Bhabha) scattering\n" 295 << "PHI DICING\n"; 296 G4Exception("G4PolarizedIonisationMo 297 JustWarning, ed); 298 } 299 } 300 // Loop checking, 03-Aug-2015, Vladimir 301 } while(grej * G4UniformRand() > xs); 302 } 303 304 // fix kinematics of delta electron 305 G4double deltaKinEnergy = x * kineticEnergy; 306 G4double deltaMomentum = 307 std::sqrt(deltaKinEnergy * (deltaKinEnergy 308 G4double cost = deltaKinEnergy * (energy + e 309 (deltaMomentum * totalMoment 310 G4double sint = 1.0 - cost * cost; 311 if(sint > 0.0) 312 sint = std::sqrt(sint); 313 314 G4ThreeVector deltaDirection(-sint * std::co 315 cost); 316 deltaDirection.rotateUz(direction); 317 318 // primary change 319 kineticEnergy -= deltaKinEnergy; 320 fParticleChange->SetProposedKineticEnergy(ki 321 322 if(kineticEnergy > DBL_MIN) 323 { 324 G4ThreeVector dir = 325 totalMomentum * direction - deltaMomentu 326 direction = dir.unit(); 327 fParticleChange->SetProposedMomentumDirect 328 } 329 330 // create G4DynamicParticle object for delta 331 G4DynamicParticle* delta = 332 new G4DynamicParticle(theElectron, deltaDi 333 vdp->push_back(delta); 334 335 // get interaction frame 336 G4ThreeVector nInteractionFrame = 337 G4PolarizationHelper::GetFrame(direction, 338 339 if(fCrossSectionCalculator) 340 { 341 // calculate mean final state polarization 342 fBeamPolarization.InvRotateAz(nInteraction 343 fTargetPolarization.InvRotateAz(nInteracti 344 fCrossSectionCalculator->Initialize(x, gam 345 fTarge 346 347 // electron/positron 348 fPositronPolarization = fCrossSectionCalcu 349 fPositronPolarization.RotateAz(nInteractio 350 351 fParticleChange->ProposePolarization(fPosi 352 353 // electron 354 fElectronPolarization = fCrossSectionCalcu 355 fElectronPolarization.RotateAz(nInteractio 356 delta->SetPolarization(fElectronPolarizati 357 fElectronPolarizati 358 } 359 else 360 { 361 fPositronPolarization = G4StokesVector::ZE 362 fElectronPolarization = G4StokesVector::ZE 363 } 364 } 365