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: G4PolarizedCompton.cc,v 1.7 2007/07/10 09:35:37 schaelic Exp $ >> 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ >> 29 // >> 30 // 26 // File name: G4PolarizedCompton 31 // File name: G4PolarizedCompton 27 // 32 // 28 // Author: Andreas Schaelicke 33 // Author: Andreas Schaelicke 29 // based on code by Michel Mair 34 // based on code by Michel Maire / Vladimir IVANTCHENKO 30 // << 31 // Class description 35 // Class description 32 // modified version respecting media and bea << 36 // 33 // using the stokes formalism << 37 // modified version respecting media and beam polarization >> 38 // using the stokes formalism >> 39 // >> 40 // Creation date: 01.05.2005 >> 41 // >> 42 // Modifications: >> 43 // >> 44 // 01-01-05, include polarization description (A.Stahl) >> 45 // 01-01-05, create asymmetry table and determine interactionlength (A.Stahl) >> 46 // 01-05-05, update handling of media polarization (A.Schalicke) >> 47 // 01-05-05, update polarized differential cross section (A.Schalicke) >> 48 // 20-05-05, added polarization transfer (A.Schalicke) >> 49 // 10-06-05, transformation between different reference frames (A.Schalicke) >> 50 // 17-10-05, correct reference frame dependence in GetMeanFreePath (A.Schalicke) >> 51 // 26-07-06, cross section recalculated (P.Starovoitov) >> 52 // 09-08-06, make it work under current geant4 release (A.Schalicke) >> 53 // 11-06-07, add PostStepGetPhysicalInteractionLength (A.Schalicke) >> 54 // ----------------------------------------------------------------------------- 34 55 35 #include "G4PolarizedCompton.hh" << 36 56 >> 57 #include "G4PolarizedCompton.hh" 37 #include "G4Electron.hh" 58 #include "G4Electron.hh" 38 #include "G4EmParameters.hh" << 59 39 #include "G4KleinNishinaCompton.hh" << 60 #include "G4StokesVector.hh" 40 #include "G4PhysicsTableHelper.hh" << 41 #include "G4PolarizationManager.hh" 61 #include "G4PolarizationManager.hh" 42 #include "G4PolarizedComptonModel.hh" 62 #include "G4PolarizedComptonModel.hh" 43 #include "G4ProductionCutsTable.hh" 63 #include "G4ProductionCutsTable.hh" 44 #include "G4StokesVector.hh" << 64 #include "G4PhysicsTableHelper.hh" 45 #include "G4SystemOfUnits.hh" << 65 #include "G4KleinNishinaCompton.hh" >> 66 #include "G4PolarizedComptonModel.hh" 46 67 47 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 48 G4PhysicsTable* G4PolarizedCompton::theAsymmet << 49 69 50 G4PolarizedCompton::G4PolarizedCompton(const G 70 G4PolarizedCompton::G4PolarizedCompton(const G4String& processName, 51 G4Proce << 71 G4ProcessType type): 52 : G4VEmProcess(processName, type) << 72 G4VEmProcess (processName, type), 53 , fType(10) << 73 buildAsymmetryTable(true), 54 , fBuildAsymmetryTable(true) << 74 useAsymmetryTable(true), 55 , fUseAsymmetryTable(true) << 75 isInitialised(false), 56 , fIsInitialised(false) << 76 selectedModel(0), 57 { << 77 mType(10), 58 SetStartFromNullFlag(true); << 78 theAsymmetryTable(NULL) 59 SetBuildTableFlag(true); << 79 { 60 SetSecondaryParticle(G4Electron::Electron()) << 80 SetLambdaBinning(90); 61 SetProcessSubType(fComptonScattering); << 81 SetMinKinEnergy(0.1*keV); 62 SetMinKinEnergyPrim(1. * MeV); << 82 SetMaxKinEnergy(100.0*GeV); 63 SetSplineFlag(true); << 64 fEmModel = nullptr; << 65 } 83 } 66 84 67 //....oooOO0OOooo........oooOO0OOooo........oo 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 68 G4PolarizedCompton::~G4PolarizedCompton() { Cl << 86 69 << 87 G4PolarizedCompton::~G4PolarizedCompton() 70 //....oooOO0OOooo........oooOO0OOooo........oo << 71 void G4PolarizedCompton::ProcessDescription(st << 72 { 88 { 73 out << "Polarized model for Compton scatteri << 89 if (theAsymmetryTable) { 74 << 75 G4VEmProcess::ProcessDescription(out); << 76 } << 77 << 78 //....oooOO0OOooo........oooOO0OOooo........oo << 79 void G4PolarizedCompton::CleanTable() << 80 { << 81 if(theAsymmetryTable) << 82 { << 83 theAsymmetryTable->clearAndDestroy(); 90 theAsymmetryTable->clearAndDestroy(); 84 delete theAsymmetryTable; 91 delete theAsymmetryTable; 85 theAsymmetryTable = nullptr; << 86 } 92 } 87 } 93 } 88 94 89 //....oooOO0OOooo........oooOO0OOooo........oo 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 90 G4bool G4PolarizedCompton::IsApplicable(const << 91 { << 92 return (&p == G4Gamma::Gamma()); << 93 } << 94 96 95 //....oooOO0OOooo........oooOO0OOooo........oo << 96 void G4PolarizedCompton::InitialiseProcess(con 97 void G4PolarizedCompton::InitialiseProcess(const G4ParticleDefinition*) 97 { 98 { 98 if(!fIsInitialised) << 99 if(!isInitialised) { 99 { << 100 isInitialised = true; 100 fIsInitialised = true; << 101 SetBuildTableFlag(true); 101 if(0 == fType) << 102 SetSecondaryParticle(G4Electron::Electron()); 102 { << 103 G4double emin = MinKinEnergy(); 103 if(nullptr == EmModel(0)) << 104 G4double emax = MaxKinEnergy(); 104 { << 105 emModel = new G4PolarizedComptonModel(); 105 SetEmModel(new G4KleinNishinaCompton() << 106 if(0 == mType) selectedModel = new G4KleinNishinaCompton(); 106 } << 107 else if(10 == mType) selectedModel = emModel; 107 } << 108 selectedModel->SetLowEnergyLimit(emin); 108 else << 109 selectedModel->SetHighEnergyLimit(emax); 109 { << 110 AddEmModel(1, selectedModel); 110 fEmModel = new G4PolarizedComptonModel() << 111 } 111 SetEmModel(fEmModel); << 112 } << 113 G4EmParameters* param = G4EmParameters::In << 114 EmModel(0)->SetLowEnergyLimit(param->MinKi << 115 EmModel(0)->SetHighEnergyLimit(param->MaxK << 116 AddEmModel(1, EmModel(0)); << 117 } << 118 } 112 } 119 113 120 //....oooOO0OOooo........oooOO0OOooo........oo 114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 121 void G4PolarizedCompton::SetModel(const G4Stri << 115 >> 116 void G4PolarizedCompton::PrintInfo() 122 { 117 { 123 if(ss == "Klein-Nishina") << 118 G4cout << " Total cross sections has a good parametrisation" 124 { << 119 << " from 10 KeV to (100/Z) GeV" 125 fType = 0; << 120 << "\n Sampling according " << selectedModel->GetName() << " model" 126 } << 121 << G4endl; 127 if(ss == "Polarized-Compton") << 122 } 128 { << 123 129 fType = 10; << 124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 130 } << 125 >> 126 void G4PolarizedCompton::SetModel(const G4String& s) >> 127 { >> 128 if(s == "Klein-Nishina") mType = 0; >> 129 if(s == "Polarized-Compton") mType = 10; 131 } 130 } 132 131 >> 132 133 //....oooOO0OOooo........oooOO0OOooo........oo 133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 134 G4double G4PolarizedCompton::GetMeanFreePath(c << 134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 135 G << 135 136 G << 136 >> 137 >> 138 G4double G4PolarizedCompton::GetMeanFreePath( >> 139 const G4Track& aTrack, >> 140 G4double previousStepSize, >> 141 G4ForceCondition* condition) 137 { 142 { 138 // *** get unploarised mean free path from l 143 // *** get unploarised mean free path from lambda table *** 139 G4double mfp = << 144 G4double mfp = G4VEmProcess::GetMeanFreePath(aTrack, previousStepSize, condition); 140 G4VEmProcess::GetMeanFreePath(aTrack, prev << 141 145 142 if(theAsymmetryTable && fUseAsymmetryTable & << 146 143 { << 147 if (theAsymmetryTable && useAsymmetryTable) { 144 mfp *= ComputeSaturationFactor(aTrack); << 148 // *** get asymmetry, if target is polarized *** 145 } << 149 const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle(); 146 if(verboseLevel >= 2) << 150 const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy(); 147 { << 151 const G4StokesVector GammaPolarization = aTrack.GetPolarization(); 148 G4cout << "G4PolarizedCompton::MeanFreePat << 152 const G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection(); 149 << G4endl; << 153 150 } << 154 G4Material* aMaterial = aTrack.GetMaterial(); 151 return mfp; << 155 G4VPhysicalVolume* aPVolume = aTrack.GetVolume(); >> 156 G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume(); >> 157 >> 158 // G4Material* bMaterial = aLVolume->GetMaterial(); >> 159 G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance(); >> 160 >> 161 const G4bool VolumeIsPolarized = polarizationManger->IsPolarized(aLVolume); >> 162 G4StokesVector ElectronPolarization = polarizationManger->GetVolumePolarization(aLVolume); >> 163 >> 164 if (!VolumeIsPolarized || mfp == DBL_MAX) return mfp; >> 165 >> 166 if (verboseLevel>=2) { >> 167 >> 168 G4cout << " Mom " << GammaDirection0 << G4endl; >> 169 G4cout << " Polarization " << GammaPolarization << G4endl; >> 170 G4cout << " MaterialPol. " << ElectronPolarization << G4endl; >> 171 G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl; >> 172 G4cout << " Log. Volume " << aLVolume->GetName() << G4endl; >> 173 G4cout << " Material " << aMaterial << G4endl; >> 174 } >> 175 >> 176 G4int midx= CurrentMaterialCutsCoupleIndex(); >> 177 G4PhysicsVector * aVector=(*theAsymmetryTable)(midx); >> 178 >> 179 G4double asymmetry=0; >> 180 if (aVector) { >> 181 G4bool isOutRange; >> 182 asymmetry = aVector->GetValue(GammaEnergy, isOutRange); >> 183 } else { >> 184 G4cout << " MaterialIndex " << midx << " is out of range \n"; >> 185 asymmetry=0; >> 186 } >> 187 >> 188 // we have to determine angle between particle motion >> 189 // and target polarisation here >> 190 // circ pol * Vec(ElectronPol)*Vec(PhotonMomentum) >> 191 // both vectors in global reference frame >> 192 >> 193 G4double pol=ElectronPolarization*GammaDirection0; >> 194 >> 195 G4double polProduct = GammaPolarization.p3() * pol; >> 196 mfp *= 1. / ( 1. + polProduct * asymmetry ); >> 197 >> 198 if (verboseLevel>=2) { >> 199 G4cout << " MeanFreePath: " << mfp / mm << " mm " << G4endl; >> 200 G4cout << " Asymmetry: " << asymmetry << G4endl; >> 201 G4cout << " PolProduct: " << polProduct << G4endl; >> 202 } >> 203 } >> 204 >> 205 return mfp; 152 } 206 } 153 207 154 //....oooOO0OOooo........oooOO0OOooo........oo << 155 G4double G4PolarizedCompton::PostStepGetPhysic 208 G4double G4PolarizedCompton::PostStepGetPhysicalInteractionLength( 156 const G4Track& aTrack, G4double previousStep << 209 const G4Track& aTrack, >> 210 G4double previousStepSize, >> 211 G4ForceCondition* condition) 157 { 212 { 158 // save previous values << 213 // *** get unploarised mean free path from lambda table *** 159 G4double nLength = theNumberOfInteractionLen << 214 G4double mfp = G4VEmProcess::PostStepGetPhysicalInteractionLength(aTrack, previousStepSize, condition); 160 G4double iLength = currentInteractionLength; << 215 161 << 216 162 // *** compute unpolarized step limit *** << 217 if (theAsymmetryTable && useAsymmetryTable) { 163 // this changes theNumberOfInteractionLength << 218 // *** get asymmetry, if target is polarized *** 164 G4double x = G4VEmProcess::PostStepGetPhysic << 219 const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle(); 165 aTrack, previousStepSize, condition); << 220 const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy(); 166 G4double x0 = x; << 221 const G4StokesVector GammaPolarization = aTrack.GetPolarization(); 167 G4double satFact = 1.0; << 222 const G4ParticleMomentum GammaDirection0 = aDynamicGamma->GetMomentumDirection(); 168 << 223 169 // *** add corrections on polarisation *** << 224 G4Material* aMaterial = aTrack.GetMaterial(); 170 if(theAsymmetryTable && fUseAsymmetryTable & << 225 G4VPhysicalVolume* aPVolume = aTrack.GetVolume(); 171 { << 226 G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume(); 172 satFact = ComputeSaturationFact << 227 173 G4double curLength = currentInteractionLen << 228 // G4Material* bMaterial = aLVolume->GetMaterial(); 174 G4double prvLength = iLength * satFact; << 229 G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance(); 175 if(nLength > 0.0) << 230 176 { << 231 const G4bool VolumeIsPolarized = polarizationManger->IsPolarized(aLVolume); 177 theNumberOfInteractionLengthLeft = << 232 G4StokesVector ElectronPolarization = polarizationManger->GetVolumePolarization(aLVolume); 178 std::max(nLength - previousStepSize / << 233 179 } << 234 if (!VolumeIsPolarized || mfp == DBL_MAX) return mfp; 180 x = theNumberOfInteractionLengthLeft * cur << 235 181 } << 236 if (verboseLevel>=2) { 182 if(verboseLevel >= 2) << 237 183 { << 238 G4cout << " Mom " << GammaDirection0 << G4endl; 184 G4cout << "G4PolarizedCompton::PostStepGPI << 239 G4cout << " Polarization " << GammaPolarization << G4endl; 185 << x / mm << " mm;" << G4endl << 240 G4cout << " MaterialPol. " << ElectronPolarization << G4endl; 186 << " unpolarized valu << 241 G4cout << " Phys. Volume " << aPVolume->GetName() << G4endl; 187 << x0 / mm << " mm." << G4endl; << 242 G4cout << " Log. Volume " << aLVolume->GetName() << G4endl; 188 } << 243 G4cout << " Material " << aMaterial << G4endl; 189 return x; << 244 } >> 245 >> 246 G4int midx= CurrentMaterialCutsCoupleIndex(); >> 247 G4PhysicsVector * aVector=(*theAsymmetryTable)(midx); >> 248 >> 249 G4double asymmetry=0; >> 250 if (aVector) { >> 251 G4bool isOutRange; >> 252 asymmetry = aVector->GetValue(GammaEnergy, isOutRange); >> 253 } else { >> 254 G4cout << " MaterialIndex " << midx << " is out of range \n"; >> 255 asymmetry=0; >> 256 } >> 257 >> 258 // we have to determine angle between particle motion >> 259 // and target polarisation here >> 260 // circ pol * Vec(ElectronPol)*Vec(PhotonMomentum) >> 261 // both vectors in global reference frame >> 262 >> 263 G4double pol=ElectronPolarization*GammaDirection0; >> 264 >> 265 G4double polProduct = GammaPolarization.p3() * pol; >> 266 mfp *= 1. / ( 1. + polProduct * asymmetry ); >> 267 >> 268 if (verboseLevel>=2) { >> 269 G4cout << " MeanFreePath: " << mfp / mm << " mm " << G4endl; >> 270 G4cout << " Asymmetry: " << asymmetry << G4endl; >> 271 G4cout << " PolProduct: " << polProduct << G4endl; >> 272 } >> 273 } >> 274 >> 275 return mfp; 190 } 276 } 191 277 192 //....oooOO0OOooo........oooOO0OOooo........oo 278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 193 G4double G4PolarizedCompton::ComputeSaturation << 194 { << 195 G4double factor = 1.0; << 196 << 197 // *** get asymmetry, if target is polarized << 198 const G4DynamicParticle* aDynamicGamma = aTr << 199 const G4double GammaEnergy = aDy << 200 const G4StokesVector GammaPolarization = << 201 G4StokesVector(aTrack.GetPolarization()); << 202 const G4ParticleMomentum GammaDirection0 = << 203 aDynamicGamma->GetMomentumDirection(); << 204 << 205 const G4Material* aMaterial = aTrack.GetMate << 206 G4VPhysicalVolume* aPVolume = aTrack.GetVolu << 207 G4LogicalVolume* aLVolume = aPVolume->GetL << 208 << 209 G4PolarizationManager* polarizationManager = << 210 G4PolarizationManager::GetInstance(); << 211 << 212 const G4bool VolumeIsPolarized = polarizatio << 213 G4StokesVector ElectronPolarization = << 214 polarizationManager->GetVolumePolarization << 215 << 216 if(VolumeIsPolarized) << 217 { << 218 if(verboseLevel >= 2) << 219 { << 220 G4cout << "G4PolarizedCompton::ComputeSa << 221 G4cout << " Mom " << GammaDirection0 << << 222 G4cout << " Polarization " << GammaPolar << 223 G4cout << " MaterialPol. " << ElectronPo << 224 G4cout << " Phys. Volume " << aPVolume-> << 225 G4cout << " Log. Volume " << aLVolume-> << 226 G4cout << " Material " << aMaterial << 227 } << 228 279 229 std::size_t midx = CurrentMa << 280 void G4PolarizedCompton::PreparePhysicsTable(const G4ParticleDefinition& part) 230 const G4PhysicsVector* aVector = nullptr; << 281 { 231 if(midx < theAsymmetryTable->size()) << 282 G4VEmProcess::PreparePhysicsTable(part); 232 { << 283 if(buildAsymmetryTable) 233 aVector = (*theAsymmetryTable)(midx); << 284 theAsymmetryTable = G4PhysicsTableHelper::PreparePhysicsTable(theAsymmetryTable); 234 } << 235 if(aVector) << 236 { << 237 G4double asymmetry = aVector->Value(Gamm << 238 << 239 // we have to determine angle between p << 240 // and target polarisation here << 241 // circ pol * Vec(ElectronPol)*Vec( << 242 // both vectors in global reference fra << 243 << 244 G4double pol = ElectronPolarizati << 245 G4double polProduct = GammaPolarization. << 246 factor /= (1. + polProduct * asymmetry); << 247 if(verboseLevel >= 2) << 248 { << 249 G4cout << " Asymmetry: " << asymme << 250 G4cout << " PolProduct: " << polPro << 251 G4cout << " Factor: " << factor << 252 } << 253 } << 254 else << 255 { << 256 G4ExceptionDescription ed; << 257 ed << "Problem with asymmetry table: mat << 258 << " is out of range or the table is << 259 G4Exception("G4PolarizedComptonModel::Co << 260 JustWarning, ed, ""); << 261 } << 262 } << 263 return factor; << 264 } 285 } 265 286 266 //....oooOO0OOooo........oooOO0OOooo........oo 287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 288 >> 289 267 void G4PolarizedCompton::BuildPhysicsTable(con 290 void G4PolarizedCompton::BuildPhysicsTable(const G4ParticleDefinition& part) 268 { 291 { 269 // *** build (unpolarized) cross section tab 292 // *** build (unpolarized) cross section tables (Lambda) 270 G4VEmProcess::BuildPhysicsTable(part); 293 G4VEmProcess::BuildPhysicsTable(part); 271 if(fBuildAsymmetryTable && fEmModel) << 294 if(buildAsymmetryTable) 272 { << 295 BuildAsymmetryTable(part); 273 G4bool isMaster = true; << 274 const G4PolarizedCompton* masterProcess = << 275 static_cast<const G4PolarizedCompton*>(G << 276 if(masterProcess && masterProcess != this) << 277 { << 278 isMaster = false; << 279 } << 280 if(isMaster) << 281 { << 282 BuildAsymmetryTable(part); << 283 } << 284 } << 285 } 296 } 286 297 287 //....oooOO0OOooo........oooOO0OOooo........oo 298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 299 >> 300 288 void G4PolarizedCompton::BuildAsymmetryTable(c 301 void G4PolarizedCompton::BuildAsymmetryTable(const G4ParticleDefinition& part) 289 { 302 { 290 // cleanup old, initialise new table << 291 CleanTable(); << 292 theAsymmetryTable = << 293 G4PhysicsTableHelper::PreparePhysicsTable( << 294 << 295 // Access to materials 303 // Access to materials 296 const G4ProductionCutsTable* theCoupleTable << 304 const G4ProductionCutsTable* theCoupleTable= 297 G4ProductionCutsTable::GetProductionCutsTa << 305 G4ProductionCutsTable::GetProductionCutsTable(); 298 G4int numOfCouples = (G4int)theCoupleTable-> << 306 size_t numOfCouples = theCoupleTable->GetTableSize(); 299 if(!theAsymmetryTable) << 307 for(size_t i=0; i<numOfCouples; ++i) { 300 { << 308 if (!theAsymmetryTable) break; 301 return; << 309 if (theAsymmetryTable->GetFlag(i)) { 302 } << 310 303 G4int nbins = LambdaBinning( << 304 G4double emin = MinKinEnergy() << 305 G4double emax = MaxKinEnergy() << 306 G4PhysicsLogVector* aVector = nullptr; << 307 G4PhysicsLogVector* bVector = nullptr; << 308 << 309 for(G4int i = 0; i < numOfCouples; ++i) << 310 { << 311 if(theAsymmetryTable->GetFlag(i)) << 312 { << 313 // create physics vector and fill it 311 // create physics vector and fill it 314 const G4MaterialCutsCouple* couple = << 312 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i); 315 theCoupleTable->GetMaterialCutsCouple( << 316 // use same parameters as for lambda 313 // use same parameters as for lambda 317 if(!aVector) << 314 G4PhysicsVector* aVector = LambdaPhysicsVector(couple); 318 { << 315 // modelManager->FillLambdaVector(aVector, couple, startFromNull); 319 aVector = new G4PhysicsLogVector(emin, << 320 bVector = aVector; << 321 } << 322 else << 323 { << 324 bVector = new G4PhysicsLogVector(*aVec << 325 } << 326 316 327 for(G4int j = 0; j <= nbins; ++j) << 317 for (G4int j = 0 ; j < LambdaBinning() ; ++j ) { 328 { << 318 G4double lowEdgeEnergy = aVector->GetLowEdgeEnergy(j); 329 G4double energy = bVector->Energy(j); << 319 G4double tasm=0.; 330 G4double tasm = 0.; << 320 G4double asym = ComputeAsymmetry(lowEdgeEnergy, couple, part, 0., tasm); 331 G4double asym = ComputeAsymmetry(ene << 321 aVector->PutValue(j,asym); 332 bVector->PutValue(j, asym); << 333 } 322 } 334 bVector->FillSecondDerivatives(); << 323 335 G4PhysicsTableHelper::SetPhysicsVector(t << 324 G4PhysicsTableHelper::SetPhysicsVector(theAsymmetryTable, i, aVector); 336 } 325 } 337 } 326 } >> 327 338 } 328 } 339 329 340 //....oooOO0OOooo........oooOO0OOooo........oo 330 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 341 G4double G4PolarizedCompton::ComputeAsymmetry( << 331 342 G4double energy, const G4MaterialCutsCouple* << 332 343 const G4ParticleDefinition& aParticle, G4dou << 333 G4double G4PolarizedCompton::ComputeAsymmetry(G4double energy, >> 334 const G4MaterialCutsCouple* couple, >> 335 const G4ParticleDefinition& aParticle, >> 336 G4double cut, >> 337 G4double & tAsymmetry) 344 { 338 { 345 G4double lAsymmetry = 0.0; 339 G4double lAsymmetry = 0.0; 346 tAsymmetry = 0; << 340 tAsymmetry=0; 347 341 >> 342 // 348 // calculate polarized cross section 343 // calculate polarized cross section 349 G4ThreeVector thePolarization = G4ThreeVecto << 344 // 350 fEmModel->SetTargetPolarization(thePolarizat << 345 G4ThreeVector thePolarization=G4ThreeVector(0.,0.,1.); 351 fEmModel->SetBeamPolarization(thePolarizatio << 346 emModel->SetTargetPolarization(thePolarization); 352 G4double sigma2 = << 347 emModel->SetBeamPolarization(thePolarization); 353 fEmModel->CrossSection(couple, &aParticle, << 348 G4double sigma2=emModel->CrossSection(couple,&aParticle,energy,cut,energy); 354 349 >> 350 // 355 // calculate unpolarized cross section 351 // calculate unpolarized cross section 356 thePolarization = G4ThreeVector(); << 352 // 357 fEmModel->SetTargetPolarization(thePolarizat << 353 thePolarization=G4ThreeVector(); 358 fEmModel->SetBeamPolarization(thePolarizatio << 354 emModel->SetTargetPolarization(thePolarization); 359 G4double sigma0 = << 355 emModel->SetBeamPolarization(thePolarization); 360 fEmModel->CrossSection(couple, &aParticle, << 356 G4double sigma0=emModel->CrossSection(couple,&aParticle,energy,cut,energy); 361 << 357 362 // determine asymmetries << 358 // determine assymmetries 363 if(sigma0 > 0.) << 359 if (sigma0>0.) { 364 { << 360 lAsymmetry=sigma2/sigma0-1.; 365 lAsymmetry = sigma2 / sigma0 - 1.; << 366 } 361 } 367 return lAsymmetry; 362 return lAsymmetry; 368 } 363 } >> 364 >> 365 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 369 366