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