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: G4PolarizedAnnihilation 31 // 32 // Author: A. Schaelicke on base of Vladimir 33 // 34 // Class Description: 35 // Polarized process of e+ annihilation into 36 37 #include "G4PolarizedAnnihilation.hh" 38 39 #include "G4DynamicParticle.hh" 40 #include "G4MaterialCutsCouple.hh" 41 #include "G4PhysicsTableHelper.hh" 42 #include "G4PhysicsVector.hh" 43 #include "G4PolarizationHelper.hh" 44 #include "G4PolarizationManager.hh" 45 #include "G4PolarizedAnnihilationModel.hh" 46 #include "G4ProductionCutsTable.hh" 47 #include "G4SystemOfUnits.hh" 48 #include "G4StokesVector.hh" 49 50 //....oooOO0OOooo........oooOO0OOooo........oo 51 G4PolarizedAnnihilation::G4PolarizedAnnihilati 52 : G4eplusAnnihilation(name) 53 , fAsymmetryTable(nullptr) 54 , fTransverseAsymmetryTable(nullptr) 55 { 56 fEmModel = new G4PolarizedAnnihilationModel( 57 SetEmModel(fEmModel); 58 } 59 60 //....oooOO0OOooo........oooOO0OOooo........oo 61 G4PolarizedAnnihilation::~G4PolarizedAnnihilat 62 63 //....oooOO0OOooo........oooOO0OOooo........oo 64 void G4PolarizedAnnihilation::CleanTables() 65 { 66 if(fAsymmetryTable) 67 { 68 fAsymmetryTable->clearAndDestroy(); 69 delete fAsymmetryTable; 70 fAsymmetryTable = nullptr; 71 } 72 if(fTransverseAsymmetryTable) 73 { 74 fTransverseAsymmetryTable->clearAndDestroy 75 delete fTransverseAsymmetryTable; 76 fTransverseAsymmetryTable = nullptr; 77 } 78 } 79 80 //....oooOO0OOooo........oooOO0OOooo........oo 81 G4double G4PolarizedAnnihilation::GetMeanFreeP 82 83 84 { 85 G4double mfp = 86 G4VEmProcess::GetMeanFreePath(track, previ 87 88 if(nullptr != fAsymmetryTable && nullptr != 89 { 90 mfp *= ComputeSaturationFactor(track); 91 } 92 if(verboseLevel >= 2) 93 { 94 G4cout << "G4PolarizedAnnihilation::MeanFr 95 << G4endl; 96 } 97 return mfp; 98 } 99 100 //....oooOO0OOooo........oooOO0OOooo........oo 101 G4double G4PolarizedAnnihilation::PostStepGetP 102 const G4Track& track, G4double previousStepS 103 { 104 // save previous values 105 G4double nLength = theNumberOfInteractionLen 106 G4double iLength = currentInteractionLength; 107 108 // *** compute unpolarized step limit *** 109 // this changes theNumberOfInteractionLength 110 G4double x = G4VEmProcess::PostStepGetPhysic 111 track, previousStepSize, condition); 112 G4double x0 = x; 113 G4double satFact = 1.0; 114 115 // *** add corrections on polarisation *** 116 if(nullptr != fAsymmetryTable && nullptr != 117 { 118 satFact = ComputeSaturationFact 119 G4double curLength = currentInteractionLen 120 G4double prvLength = iLength * satFact; 121 if(nLength > 0.0) 122 { 123 theNumberOfInteractionLengthLeft = 124 std::max(nLength - previousStepSize / 125 } 126 x = theNumberOfInteractionLengthLeft * cur 127 } 128 if(verboseLevel >= 2) 129 { 130 G4cout << "G4PolarizedAnnihilation::PostSt 131 << x / mm << " mm;" << G4endl 132 << " unpola 133 << std::setprecision(8) << x0 / mm 134 } 135 return x; 136 } 137 138 //....oooOO0OOooo........oooOO0OOooo........oo 139 G4double G4PolarizedAnnihilation::ComputeSatur 140 { 141 const G4Material* aMaterial = track.GetMater 142 G4VPhysicalVolume* aPVolume = track.GetVolum 143 G4LogicalVolume* aLVolume = aPVolume->GetL 144 145 G4PolarizationManager* polarizationManager = 146 G4PolarizationManager::GetInstance(); 147 148 const G4bool volumeIsPolarized = polarizatio 149 G4StokesVector electronPolarization = 150 polarizationManager->GetVolumePolarization 151 152 G4double factor = 1.0; 153 154 if(volumeIsPolarized) 155 { 156 // *** get asymmetry, if target is polariz 157 const G4DynamicParticle* aDynamicPositron 158 const G4double positronEnergy = aDynamicPo 159 const G4StokesVector positronPolarization 160 G4StokesVector(track.GetPolarization()); 161 const G4ParticleMomentum positronDirection 162 aDynamicPositron->GetMomentumDirection() 163 164 if(verboseLevel >= 2) 165 { 166 G4cout << "G4PolarizedAnnihilation::Comp 167 G4cout << " Mom " << positronDirection0 168 G4cout << " Polarization " << positronPo 169 G4cout << " MaterialPol. " << electronPo 170 G4cout << " Phys. Volume " << aPVolume-> 171 G4cout << " Log. Volume " << aLVolume-> 172 G4cout << " Material " << aMaterial 173 } 174 175 std::size_t midx = CurrentMa 176 const G4PhysicsVector* aVector = nullptr; 177 const G4PhysicsVector* bVector = nullptr; 178 if(midx < fAsymmetryTable->size()) 179 { 180 aVector = (*fAsymmetryTable)(midx); 181 } 182 if(midx < fTransverseAsymmetryTable->size( 183 { 184 bVector = (*fTransverseAsymmetryTable)(m 185 } 186 if(aVector && bVector) 187 { 188 G4double lAsymmetry = aVector->Value(pos 189 G4double tAsymmetry = bVector->Value(pos 190 G4double polZZ = 191 positronPolarization.z() * (electronPo 192 G4double polXX = 193 positronPolarization.x() * 194 (electronPolarization * 195 G4PolarizationHelper::GetParticleFram 196 G4double polYY = 197 positronPolarization.y() * 198 (electronPolarization * 199 G4PolarizationHelper::GetParticleFram 200 201 factor /= (1. + polZZ * lAsymmetry + (po 202 203 if(verboseLevel >= 2) 204 { 205 G4cout << " Asymmetry: " << lAsymm 206 << G4endl; 207 G4cout << " PolProduct: " << polXX 208 << G4endl; 209 G4cout << " Factor: " << factor 210 } 211 } 212 else 213 { 214 G4ExceptionDescription ed; 215 ed << "Problem with asymmetry tables: ma 216 << " is out of range or tables are no 217 G4Exception("G4PolarizedAnnihilation::Co 218 JustWarning, ed, ""); 219 } 220 } 221 return factor; 222 } 223 224 //....oooOO0OOooo........oooOO0OOooo........oo 225 void G4PolarizedAnnihilation::BuildPhysicsTabl 226 const G4ParticleDefinition& part) 227 { 228 G4VEmProcess::BuildPhysicsTable(part); 229 if(isTheMaster) 230 { 231 BuildAsymmetryTables(part); 232 } 233 } 234 235 //....oooOO0OOooo........oooOO0OOooo........oo 236 void G4PolarizedAnnihilation::BuildAsymmetryTa 237 const G4ParticleDefinition& part) 238 { 239 // cleanup old, initialise new table 240 CleanTables(); 241 fAsymmetryTable = G4PhysicsTableHelper::Prep 242 fTransverseAsymmetryTable = 243 G4PhysicsTableHelper::PreparePhysicsTable( 244 if(nullptr == fAsymmetryTable) return; 245 246 // Access to materials 247 const G4ProductionCutsTable* theCoupleTable 248 G4ProductionCutsTable::GetProductionCutsTa 249 G4int numOfCouples = (G4int)theCoupleTable-> 250 for(G4int i = 0; i < numOfCouples; ++i) 251 { 252 if(fAsymmetryTable->GetFlag(i)) 253 { 254 // create physics vector and fill it 255 const G4MaterialCutsCouple* couple = 256 theCoupleTable->GetMaterialCutsCouple( 257 258 // use same parameters as for lambda 259 G4PhysicsVector* aVector = LambdaPhysics 260 G4PhysicsVector* tVector = LambdaPhysics 261 G4int nn = (G4int)aVector->GetVectorLeng 262 for(G4int j = 0; j < nn; ++j) 263 { 264 G4double energy = aVector->Energy(j); 265 G4double tasm = 0.; 266 G4double asym = ComputeAsymmetry(energ 267 aVector->PutValue(j, asym); 268 tVector->PutValue(j, tasm); 269 } 270 if(aVector->GetSpline()) { 271 aVector->FillSecondDerivatives(); 272 tVector->FillSecondDerivatives(); 273 } 274 G4PhysicsTableHelper::SetPhysicsVector(f 275 G4PhysicsTableHelper::SetPhysicsVector(f 276 t 277 } 278 } 279 } 280 281 //....oooOO0OOooo........oooOO0OOooo........oo 282 G4double G4PolarizedAnnihilation::ComputeAsymm 283 G4double energy, const G4MaterialCutsCouple* 284 const G4ParticleDefinition& aParticle, G4dou 285 { 286 G4double lAsymmetry = 0.0; 287 tAsymmetry = 0.0; 288 289 // calculate polarized cross section 290 G4ThreeVector targetPolarization = G4ThreeVe 291 fEmModel->SetTargetPolarization(targetPolari 292 fEmModel->SetBeamPolarization(targetPolariza 293 G4double sigma2 = 294 fEmModel->CrossSection(couple, &aParticle, 295 296 // calculate transversely polarized cross se 297 targetPolarization = G4ThreeVector(1., 0., 0 298 fEmModel->SetTargetPolarization(targetPolari 299 fEmModel->SetBeamPolarization(targetPolariza 300 G4double sigma3 = 301 fEmModel->CrossSection(couple, &aParticle, 302 303 // calculate unpolarized cross section 304 targetPolarization = G4ThreeVector(); 305 fEmModel->SetTargetPolarization(targetPolari 306 fEmModel->SetBeamPolarization(targetPolariza 307 G4double sigma0 = 308 fEmModel->CrossSection(couple, &aParticle, 309 310 // determine asymmetries 311 if(sigma0 > 0.) 312 { 313 lAsymmetry = sigma2 / sigma0 - 1.; 314 tAsymmetry = sigma3 / sigma0 - 1.; 315 } 316 return lAsymmetry; 317 } 318 319 //....oooOO0OOooo........oooOO0OOooo........oo 320 void G4PolarizedAnnihilation::ProcessDescripti 321 { 322 out << "Polarized model for positron annihil 323 } 324