Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // >> 27 // $Id: G4PolarizedComptonModel.cc 69847 2013-05-16 09:36:18Z gcosmo $ >> 28 // 26 // ------------------------------------------- 29 // ------------------------------------------------------------------- 27 // 30 // 28 // Geant4 Class file << 31 // GEANT4 Class file >> 32 // 29 // 33 // 30 // File name: G4PolarizedComptonModel 34 // File name: G4PolarizedComptonModel 31 // 35 // 32 // Author: Andreas Schaelicke 36 // Author: Andreas Schaelicke >> 37 // >> 38 // Creation date: 01.05.2005 >> 39 // >> 40 // Modifications: >> 41 // 18-07-06 use newly calculated cross sections (P. Starovoitov) >> 42 // 21-08-05 update interface (A. Schaelicke) >> 43 // >> 44 // Class Description: >> 45 // >> 46 // ------------------------------------------------------------------- >> 47 // >> 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 33 50 34 #include "G4PolarizedComptonModel.hh" 51 #include "G4PolarizedComptonModel.hh" 35 << 36 #include "G4Exp.hh" << 37 #include "G4Log.hh" << 38 #include "G4ParticleChangeForGamma.hh" << 39 #include "G4PhysicalConstants.hh" 52 #include "G4PhysicalConstants.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" 40 #include "G4PolarizationManager.hh" 61 #include "G4PolarizationManager.hh" 41 #include "G4PolarizationHelper.hh" 62 #include "G4PolarizationHelper.hh" 42 #include "G4PolarizedComptonXS.hh" << 63 #include "G4PolarizedComptonCrossSection.hh" 43 #include "G4StokesVector.hh" << 44 #include "G4SystemOfUnits.hh" << 45 64 46 //....oooOO0OOooo........oooOO0OOooo........oo 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 66 47 G4PolarizedComptonModel::G4PolarizedComptonMod 67 G4PolarizedComptonModel::G4PolarizedComptonModel(const G4ParticleDefinition*, 48 << 68 const G4String& nam) 49 : G4KleinNishinaCompton(nullptr, nam) << 69 : G4KleinNishinaCompton(0,nam), 50 , fVerboseLevel(0) << 70 verboseLevel(0) 51 { 71 { 52 fCrossSectionCalculator = new G4PolarizedCom << 72 crossSectionCalculator=new G4PolarizedComptonCrossSection(); 53 fBeamPolarization = G4StokesVector::ZE << 54 fTargetPolarization = G4StokesVector::ZE << 55 } 73 } 56 74 57 //....oooOO0OOooo........oooOO0OOooo........oo 75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 76 58 G4PolarizedComptonModel::~G4PolarizedComptonMo 77 G4PolarizedComptonModel::~G4PolarizedComptonModel() 59 { 78 { 60 delete fCrossSectionCalculator; << 79 if (crossSectionCalculator) delete crossSectionCalculator; 61 } 80 } 62 81 63 //....oooOO0OOooo........oooOO0OOooo........oo << 82 64 G4double G4PolarizedComptonModel::ComputeAsymm << 83 65 << 84 G4double G4PolarizedComptonModel::ComputeAsymmetryPerAtom >> 85 (G4double gammaEnergy, G4double /*Z*/) >> 86 66 { 87 { 67 G4double asymmetry = 0.0; << 88 G4double asymmetry = 0.0 ; 68 89 69 G4double k0 = gammaEnergy / electron_mass_c2 << 90 G4double k0 = gammaEnergy / electron_mass_c2 ; 70 G4double k1 = 1. + 2. * k0; << 91 G4double k1 = 1 + 2*k0 ; 71 92 72 asymmetry = -k0; << 93 asymmetry = -k0; 73 asymmetry *= << 94 asymmetry *= (k0 + 1.)*sqr(k1)*std::log(k1) - 2.*k0*(5.*sqr(k0) + 4.*k0 + 1.); 74 (k0 + 1.) * sqr(k1) * G4Log(k1) - 2. * k0 << 95 asymmetry /= ((k0 - 2.)*k0 -2.)*sqr(k1)*std::log(k1) + 2.*k0*(k0*(k0 + 1.)*(k0 + 8.) + 2.); 75 asymmetry /= ((k0 - 2.) * k0 - 2.) * sqr(k1) << 76 2. * k0 * (k0 * (k0 + 1.) * (k0 << 77 96 78 if(asymmetry > 1.) << 97 // G4cout<<"energy = "<<GammaEnergy<<" asymmetry = "<<asymmetry<<"\t\t GAM = "<<k0<<G4endl; 79 { << 98 if (asymmetry>1.) G4cout<<"ERROR in G4PolarizedComptonModel::ComputeAsymmetryPerAtom"<<G4endl; 80 G4ExceptionDescription ed; << 81 ed << "ERROR in G4PolarizedComptonModel::C << 82 << " asymmetry = " << asymmetry << "\n" << 83 G4Exception("G4PolarizedComptonModel::Comp << 84 JustWarning, ed); << 85 } << 86 99 87 return asymmetry; << 100 return asymmetry; 88 } 101 } 89 102 90 //....oooOO0OOooo........oooOO0OOooo........oo << 103 91 G4double G4PolarizedComptonModel::ComputeCross 104 G4double G4PolarizedComptonModel::ComputeCrossSectionPerAtom( 92 const G4ParticleDefinition* pd, G4double kin << 105 const G4ParticleDefinition* pd, 93 G4double cut, G4double emax) << 106 G4double kinEnergy, >> 107 G4double Z, >> 108 G4double A, >> 109 G4double cut, >> 110 G4double emax) 94 { 111 { 95 G4double xs = G4KleinNishinaCompton::Compute << 112 double xs = 96 pd, kinEnergy, Z, A, cut, emax); << 113 G4KleinNishinaCompton::ComputeCrossSectionPerAtom(pd,kinEnergy, 97 G4double polzz = fBeamPolarization.p3() * fT << 114 Z,A,cut,emax); 98 if(polzz > 0.0) << 115 G4double polzz = theBeamPolarization.p3()*theTargetPolarization.z(); 99 { << 116 if (polzz!=0) { 100 G4double asym = ComputeAsymmetryPerAtom(ki << 117 G4double asym=ComputeAsymmetryPerAtom(kinEnergy, Z); 101 xs *= (1. + polzz * asym); << 118 xs*=(1.+polzz*asym); 102 } 119 } 103 return xs; 120 return xs; 104 } 121 } 105 122 >> 123 106 //....oooOO0OOooo........oooOO0OOooo........oo 124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 107 void G4PolarizedComptonModel::SampleSecondarie << 125 108 std::vector<G4DynamicParticle*>* fvect, cons << 126 void G4PolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, 109 const G4DynamicParticle* aDynamicGamma, G4do << 127 const G4MaterialCutsCouple*, >> 128 const G4DynamicParticle* aDynamicGamma, >> 129 G4double, >> 130 G4double) 110 { 131 { 111 // do nothing below the threshold << 132 const G4Track * aTrack = fParticleChange->GetCurrentTrack(); 112 if(aDynamicGamma->GetKineticEnergy() <= LowE << 133 G4VPhysicalVolume* aPVolume = aTrack->GetVolume(); 113 { << 134 G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume(); 114 return; << 135 115 } << 136 if (verboseLevel>=1) >> 137 G4cout<<"G4PolarizedComptonModel::SampleSecondaries in " >> 138 << aLVolume->GetName() <<G4endl; 116 139 117 const G4Track* aTrack = fParticleChang << 140 G4PolarizationManager * polarizationManager = G4PolarizationManager::GetInstance(); 118 G4VPhysicalVolume* aPVolume = aTrack->GetVol << 119 G4LogicalVolume* aLVolume = aPVolume->GetL << 120 << 121 if(fVerboseLevel >= 1) << 122 { << 123 G4cout << "G4PolarizedComptonModel::Sample << 124 << aLVolume->GetName() << G4endl; << 125 } << 126 G4PolarizationManager* polarizationManager = << 127 G4PolarizationManager::GetInstance(); << 128 141 129 // obtain polarization of the beam 142 // obtain polarization of the beam 130 fBeamPolarization = G4StokesVector(aDynamicG << 143 theBeamPolarization = aDynamicGamma->GetPolarization(); 131 fBeamPolarization.SetPhoton(); << 144 theBeamPolarization.SetPhoton(); 132 145 133 // obtain polarization of the media 146 // obtain polarization of the media 134 G4bool targetIsPolarized = polarizationManag << 147 const G4bool targetIsPolarized = polarizationManager->IsPolarized(aLVolume); 135 fTargetPolarization = polarizationManager->G << 148 theTargetPolarization = polarizationManager->GetVolumePolarization(aLVolume); 136 149 137 // if beam is linear polarized or target is << 150 // if beam is linear polarized or target is transversely polarized 138 // determine the angle to x-axis 151 // determine the angle to x-axis 139 // (assumes same PRF as in the polarization 152 // (assumes same PRF as in the polarization definition) >> 153 140 G4ThreeVector gamDirection0 = aDynamicGamma- 154 G4ThreeVector gamDirection0 = aDynamicGamma->GetMomentumDirection(); 141 155 142 // transfer fTargetPolarization << 156 // transfere theTargetPolarization 143 // into the gamma frame (problem electron is 157 // into the gamma frame (problem electron is at rest) 144 if(targetIsPolarized) << 158 if (targetIsPolarized) 145 { << 159 theTargetPolarization.rotateUz(gamDirection0); 146 fTargetPolarization.rotateUz(gamDirection0 << 147 } << 148 // The scattered gamma energy is sampled acc << 149 // Klein - Nishina formula. << 150 // The random number techniques of Butcher & << 151 // (Nuc Phys 20(1960),15). << 152 // Note : Effects due to binding of atomic e << 153 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 154 G4double gamEnergy0 = aDynamicGamma->GetKine 166 G4double gamEnergy0 = aDynamicGamma->GetKineticEnergy(); 155 G4double E0_m = gamEnergy0 / electron_ << 167 G4double E0_m = gamEnergy0 / electron_mass_c2 ; 156 168 157 // sample the energy rate of the scattered g << 158 G4double epsilon, sint2; << 159 G4double onecost = 0.0; << 160 G4double Phi = 0.0; << 161 G4double greject = 1.0; << 162 G4double cosTeta = 1.0; << 163 G4double sinTeta = 0.0; << 164 << 165 G4double eps0 = 1. / (1. + 2. * E0_m); << 166 G4double epsilon0sq = eps0 * eps0; << 167 G4double alpha1 = -G4Log(eps0); << 168 G4double alpha2 = alpha1 + 0.5 * (1. - e << 169 << 170 G4double polarization = fBeamPolarization.p3 << 171 << 172 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra << 173 G4int nloop = 0; << 174 G4bool end = fals << 175 169 176 G4double rndm[3]; << 170 // 177 << 171 // sample the energy rate of the scattered gamma 178 do << 172 // 179 { << 173 180 do << 174 G4double epsilon, epsilonsq, onecost, sint2, greject ; 181 { << 175 182 ++nloop; << 176 G4double eps0 = 1./(1. + 2.*E0_m); 183 // false interaction if too many iterati << 177 G4double epsilon0sq = eps0*eps0; 184 if(nloop > fLoopLim) << 178 G4double alpha1 = - std::log(eps0); 185 { << 179 G4double alpha2 = 0.5*(1.- epsilon0sq); 186 PrintWarning(aDynamicGamma, nloop, gre << 180 187 "too many iterations"); << 181 G4double polarization = theBeamPolarization.p3()*theTargetPolarization.p3(); 188 return; << 182 do { 189 } << 183 if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) { 190 << 184 epsilon = std::exp(-alpha1*G4UniformRand()); // epsilon0**r 191 // 3 random numbers to sample scattering << 185 epsilonsq = epsilon*epsilon; 192 rndmEngineMod->flatArray(3, rndm); << 186 193 << 187 } else { 194 if(alpha1 > alpha2 * rndm[0]) << 188 epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand(); 195 { << 189 epsilon = std::sqrt(epsilonsq); 196 epsilon = G4Exp(-alpha1 * rndm[1]); << 190 }; 197 } << 191 198 else << 192 onecost = (1.- epsilon)/(epsilon*E0_m); 199 { << 193 sint2 = onecost*(2.-onecost); 200 epsilon = std::sqrt(epsilon0sq + (1. - << 194 201 } << 195 202 << 196 G4double gdiced = 2.*(1./epsilon+epsilon); 203 onecost = (1. - epsilon) / (epsilon * E0 << 197 G4double gdist = 1./epsilon + epsilon - sint2 204 sint2 = onecost * (2. - onecost); << 198 - polarization*(1./epsilon-epsilon)*(1.-onecost); 205 << 199 206 G4double gdiced = 2. * (1. / epsilon + e << 200 greject = gdist/gdiced; 207 G4double gdist = 1. / epsilon + epsilon << 201 208 polarization * (1. / ep << 202 if (greject>1) G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n" 209 << 203 <<" costh rejection does not work properly: "<<greject<<G4endl; 210 greject = gdist / gdiced; << 204 211 << 205 } while (greject < G4UniformRand()); 212 if(greject > 1.0) << 206 213 { << 207 // 214 PrintWarning(aDynamicGamma, nloop, gre << 208 // scattered gamma angles. ( Z - axis along the parent gamma) 215 "theta majoranta wrong"); << 209 // 216 } << 210 217 // Loop checking, 03-Aug-2015, Vladimir << 211 G4double cosTeta = 1. - onecost; 218 } while(greject < rndm[2]); << 212 G4double sinTeta = std::sqrt (sint2); 219 << 213 G4double Phi; 220 // assuming phi loop successful << 214 do { 221 end = true; << 215 Phi = twopi * G4UniformRand(); 222 << 216 G4double gdiced = 1./epsilon + epsilon - sint2 223 // scattered gamma angles. ( Z - axis alon << 217 + std::abs(theBeamPolarization.p3())* 224 cosTeta = 1. - onecost; << 218 ( std::abs((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3()) 225 sinTeta = std::sqrt(sint2); << 219 +(1.-epsilon)*sinTeta*(std::sqrt(sqr(theTargetPolarization.p1()) 226 do << 220 + sqr(theTargetPolarization.p2())))) 227 { << 221 +sint2*(std::sqrt(sqr(theBeamPolarization.p1()) + sqr(theBeamPolarization.p2()))); 228 ++nloop; << 222 229 << 223 G4double gdist = 1./epsilon + epsilon - sint2 230 // 2 random numbers to sample scattering << 224 + theBeamPolarization.p3()* 231 rndmEngineMod->flatArray(2, rndm); << 225 ((1./epsilon-epsilon)*cosTeta*theTargetPolarization.p3() 232 << 226 +(1.-epsilon)*sinTeta*(std::cos(Phi)*theTargetPolarization.p1()+ 233 // false interaction if too many iterati << 227 std::sin(Phi)*theTargetPolarization.p2())) 234 Phi = twopi * rndm[0]; << 228 -sint2*(std::cos(2.*Phi)*theBeamPolarization.p1() 235 if(nloop > fLoopLim) << 229 +std::sin(2.*Phi)*theBeamPolarization.p2()); 236 { << 230 greject = gdist/gdiced; 237 PrintWarning(aDynamicGamma, nloop, gre << 231 238 "too many iterations"); << 232 if (greject>1.+1.e-10 || greject<0) G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n" 239 return; << 233 <<" phi rejection does not work properly: "<<greject<<G4endl; 240 } << 234 241 << 235 if (greject<1.e-3) { 242 G4double gdiced = 1. / epsilon + epsilon << 236 G4cout<<"ERROR in PolarizedComptonScattering::PostStepDoIt\n" 243 std::abs(fBeamPolariza << 237 <<" phi rejection does not work properly: "<<greject<<"\n"; 244 (std::abs((1. / epsi << 238 G4cout<<" greject="<<greject<<" phi="<<Phi<<" cost="<<cosTeta<<"\n"; 245 fTargetPol << 239 G4cout<<" gdiced="<<gdiced<<" gdist="<<gdist<<"\n"; 246 (1. - epsilon) * si << 240 G4cout<<" eps="<<epsilon<<" 1/eps="<<1./epsilon<<"\n"; 247 (std::sqrt(sqr(fT << 241 } 248 sqr(fT << 242 249 sint2 * (std::sqrt(sqr << 243 } while (greject < G4UniformRand()); 250 sqr << 244 G4double dirx = sinTeta*std::cos(Phi), diry = sinTeta*std::sin(Phi), dirz = cosTeta; 251 << 252 G4double gdist = << 253 1. / epsilon + epsilon - sint2 + << 254 fBeamPolarization.p3() * << 255 ((1. / epsilon - epsilon) * cosTeta << 256 (1. - epsilon) * sinTeta * << 257 (std::cos(Phi) * fTargetPolarizat << 258 std::sin(Phi) * fTargetPolarizat << 259 sint2 * (std::cos(2. * Phi) * fBeamPol << 260 std::sin(2. * Phi) * fBeamPol << 261 greject = gdist / gdiced; << 262 << 263 if(greject > 1.0) << 264 { << 265 PrintWarning(aDynamicGamma, nloop, gre << 266 "phi majoranta wrong"); << 267 } << 268 << 269 if(greject < 1.e-3) << 270 { << 271 PrintWarning(aDynamicGamma, nloop, gre << 272 "phi loop ineffective"); << 273 // restart theta loop << 274 end = false; << 275 break; << 276 } << 277 << 278 // Loop checking, 03-Aug-2015, Vladimir << 279 } while(greject < rndm[1]); << 280 } while(!end); << 281 G4double dirx = sinTeta * std::cos(Phi); << 282 G4double diry = sinTeta * std::sin(Phi); << 283 G4double dirz = cosTeta; << 284 245 >> 246 // 285 // update G4VParticleChange for the scattere 247 // update G4VParticleChange for the scattered gamma 286 G4ThreeVector gamDirection1(dirx, diry, dirz << 248 // >> 249 >> 250 G4ThreeVector gamDirection1 ( dirx,diry,dirz ); 287 gamDirection1.rotateUz(gamDirection0); 251 gamDirection1.rotateUz(gamDirection0); 288 G4double gamEnergy1 = epsilon * gamEnergy0; << 252 G4double gamEnergy1 = epsilon*gamEnergy0; >> 253 fParticleChange->SetProposedKineticEnergy(gamEnergy1); 289 254 290 G4double edep = 0.0; << 255 291 if(gamEnergy1 > lowestSecondaryEnergy) << 256 if(gamEnergy1 > lowestGammaEnergy) { 292 { << 293 fParticleChange->ProposeMomentumDirection( 257 fParticleChange->ProposeMomentumDirection(gamDirection1); 294 fParticleChange->SetProposedKineticEnergy( << 258 } else { 295 } << 296 else << 297 { << 298 fParticleChange->ProposeTrackStatus(fStopA 259 fParticleChange->ProposeTrackStatus(fStopAndKill); 299 fParticleChange->SetProposedKineticEnergy( << 260 gamEnergy1 += fParticleChange->GetLocalEnergyDeposit(); 300 edep = gamEnergy1; << 261 fParticleChange->ProposeLocalEnergyDeposit(gamEnergy1); 301 } 262 } >> 263 >> 264 // >> 265 // kinematic of the scattered electron >> 266 // 302 267 303 // calculate Stokes vector of final state ph << 268 G4double eKinEnergy = gamEnergy0 - gamEnergy1; 304 G4ThreeVector nInteractionFrame = << 269 G4ThreeVector eDirection = gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1; 305 G4PolarizationHelper::GetFrame(gamDirectio << 270 eDirection = eDirection.unit(); 306 << 307 // transfer fBeamPolarization and fTargetPol << 308 // into the interaction frame (note electron << 309 if(fVerboseLevel >= 1) << 310 { << 311 G4cout << "=============================== << 312 G4cout << " nInteractionFrame = " << nInte << 313 G4cout << " GammaDirection0 = " << gamDire << 314 G4cout << " gammaPolarization = " << fBeam << 315 G4cout << " electronPolarization = " << fT << 316 } << 317 271 318 fBeamPolarization.InvRotateAz(nInteractionFr << 272 // 319 fTargetPolarization.InvRotateAz(nInteraction << 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 } 320 296 321 if(fVerboseLevel >= 1) << 297 theBeamPolarization.InvRotateAz(nInteractionFrame,gamDirection0); 322 { << 298 theTargetPolarization.InvRotateAz(nInteractionFrame,gamDirection0); 323 G4cout << "------------------------------- << 299 324 G4cout << " gammaPolarization = " << fBeam << 300 if (verboseLevel>=1) { 325 G4cout << " electronPolarization = " << fT << 301 G4cout << "----------------------------------------\n"; 326 G4cout << "------------------------------- << 302 G4cout << " gammaPolarization = " <<theBeamPolarization<<"\n"; 327 } << 303 G4cout << " electronPolarization = " <<theTargetPolarization<<"\n"; >> 304 G4cout << "----------------------------------------\n"; >> 305 } 328 306 329 // initialize the polarization transfer matr << 307 // initialize the polarization transfer matrix 330 fCrossSectionCalculator->Initialize(epsilon, << 308 crossSectionCalculator->Initialize(epsilon,E0_m,0., 331 fTargetP << 309 theBeamPolarization, >> 310 theTargetPolarization,2); >> 311 } 332 312 333 if(gamEnergy1 > lowestSecondaryEnergy) << 313 // if(eKinEnergy > DBL_MIN) 334 { 314 { 335 // in interaction frame 315 // in interaction frame 336 // calculate polarization transfer to the 316 // calculate polarization transfer to the photon (in interaction plane) 337 fFinalGammaPolarization = fCrossSectionCal << 317 finalGammaPolarization = crossSectionCalculator->GetPol2(); 338 if(fVerboseLevel >= 1) << 318 if (verboseLevel>=1) G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n"; 339 { << 319 finalGammaPolarization.SetPhoton(); 340 G4cout << " gammaPolarization1 = " << fF << 341 } << 342 fFinalGammaPolarization.SetPhoton(); << 343 320 344 // translate polarization into particle re 321 // translate polarization into particle reference frame 345 fFinalGammaPolarization.RotateAz(nInteract << 322 finalGammaPolarization.RotateAz(nInteractionFrame,gamDirection1); 346 if(fFinalGammaPolarization.mag() > 1. + 1. << 323 //store polarization vector 347 { << 324 fParticleChange->ProposePolarization(finalGammaPolarization); 348 G4ExceptionDescription ed; << 325 if (finalGammaPolarization.mag() > 1.+1.e-8){ 349 ed << "ERROR in Polarizaed Compton Scatt << 326 G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl; 350 ed << "Polarization of final photon more << 327 G4cout<<"Polarization of final photon more than 100%"<<G4endl; 351 ed << fFinalGammaPolarization << 328 G4cout<<finalGammaPolarization<<" mag = "<<finalGammaPolarization.mag()<<G4endl; 352 << " mag = " << fFinalGammaPolarizati << 353 G4Exception("G4PolarizedComptonModel::Sa << 354 FatalException, ed); << 355 } 329 } 356 // store polarization vector << 330 if (verboseLevel>=1) { 357 fParticleChange->ProposePolarization(fFina << 331 G4cout << " gammaPolarization1 = " <<finalGammaPolarization<<"\n"; 358 if(fVerboseLevel >= 1) << 332 G4cout << " GammaDirection1 = " <<gamDirection1<<"\n"; 359 { << 360 G4cout << " gammaPolarization1 = " << fF << 361 G4cout << " GammaDirection1 = " << gamDi << 362 } 333 } 363 } 334 } 364 335 365 // kinematic of the scattered electron << 336 // if (ElecKineEnergy > fminimalEnergy) { 366 G4double eKinEnergy = gamEnergy0 - gamEnergy << 367 << 368 if(eKinEnergy > lowestSecondaryEnergy) << 369 { 337 { 370 G4ThreeVector eDirection = << 338 finalElectronPolarization = crossSectionCalculator->GetPol3(); 371 gamEnergy0 * gamDirection0 - gamEnergy1 << 339 if (verboseLevel>=1) 372 eDirection = eDirection.unit(); << 340 G4cout << " electronPolarization1 = " <<finalElectronPolarization<<"\n"; 373 << 341 374 finalElectronPolarization = fCrossSectionC << 375 if(fVerboseLevel >= 1) << 376 { << 377 G4cout << " electronPolarization1 = " << << 378 << G4endl; << 379 } << 380 // transfer into particle reference frame 342 // transfer into particle reference frame 381 finalElectronPolarization.RotateAz(nIntera << 343 finalElectronPolarization.RotateAz(nInteractionFrame,eDirection); 382 if(fVerboseLevel >= 1) << 344 if (verboseLevel>=1) { 383 { << 345 G4cout << " electronPolarization1 = " <<finalElectronPolarization<<"\n"; 384 G4cout << " electronPolarization1 = " << << 346 G4cout << " ElecDirection = " <<eDirection<<"\n"; 385 << G4endl << " ElecDirection = " << 386 } 347 } >> 348 } >> 349 if (verboseLevel>=1) >> 350 G4cout << "========================================\n"; >> 351 >> 352 >> 353 if(eKinEnergy > DBL_MIN) { 387 354 388 // create G4DynamicParticle object for the 355 // create G4DynamicParticle object for the electron. 389 G4DynamicParticle* aElectron = << 356 G4DynamicParticle* aElectron = new G4DynamicParticle(theElectron,eDirection,eKinEnergy); 390 new G4DynamicParticle(theElectron, eDire << 357 //store polarization vector 391 // store polarization vector << 358 if (finalElectronPolarization.mag() > 1.+1.e-8){ 392 if(finalElectronPolarization.mag() > 1. + << 359 G4cout<<"ERROR in Polarizaed Compton Scattering !"<<G4endl; 393 { << 360 G4cout<<"Polarization of final electron more than 100%"<<G4endl; 394 G4ExceptionDescription ed; << 361 G4cout<<finalElectronPolarization<<" mag = "<<finalElectronPolarization.mag()<<G4endl; 395 ed << "ERROR in Polarized Compton Scatte << 396 ed << "Polarization of final electron mo << 397 ed << finalElectronPolarization << 398 << " mag = " << finalElectronPolariza << 399 G4Exception("G4PolarizedComptonModel::Sa << 400 FatalException, ed); << 401 } 362 } 402 aElectron->SetPolarization(finalElectronPo 363 aElectron->SetPolarization(finalElectronPolarization.p1(), 403 finalElectronPo << 364 finalElectronPolarization.p2(), 404 finalElectronPo << 365 finalElectronPolarization.p3()); 405 fvect->push_back(aElectron); 366 fvect->push_back(aElectron); 406 } 367 } 407 else << 408 { << 409 edep += eKinEnergy; << 410 } << 411 // energy balance << 412 if(edep > 0.0) << 413 { << 414 fParticleChange->ProposeLocalEnergyDeposit << 415 } << 416 } 368 } 417 369 418 //....oooOO0OOooo........oooOO0OOooo........oo 370 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 419 void G4PolarizedComptonModel::PrintWarning(con << 371 420 G4i << 372 421 G4d << 422 con << 423 { << 424 G4ExceptionDescription ed; << 425 ed << "Problem of scattering sampling: " << << 426 << "Niter= " << nloop << " grej= " << gre << 427 << " cos(theta)= " << 1.0 - onecos << " p << 428 << "Gamma E(MeV)= " << dp->GetKineticEner << 429 << " dir= " << dp->GetMomentumDirection() << 430 << " pol= " << dp->GetPolarization(); << 431 G4Exception("G4PolarizedComptonModel::Sample << 432 JustWarning, ed, ""); << 433 } << 434 373