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: G4LivermorePolarizedGammaConversionModel.hh,v 1.2 2010-11-23 16:42:15 flongo Exp $ 26 // 27 // 27 // Authors: G.Depaola & F.Longo 28 // Authors: G.Depaola & F.Longo 28 // 29 // 29 // 30 // 30 31 31 #include "G4LivermorePolarizedGammaConversionM 32 #include "G4LivermorePolarizedGammaConversionModel.hh" 32 #include "G4PhysicalConstants.hh" 33 #include "G4PhysicalConstants.hh" 33 #include "G4SystemOfUnits.hh" 34 #include "G4SystemOfUnits.hh" 34 #include "G4Electron.hh" 35 #include "G4Electron.hh" 35 #include "G4Positron.hh" 36 #include "G4Positron.hh" 36 #include "G4ParticleChangeForGamma.hh" 37 #include "G4ParticleChangeForGamma.hh" 37 #include "G4Log.hh" 38 #include "G4Log.hh" 38 #include "G4AutoLock.hh" << 39 #include "G4Exp.hh" 39 #include "G4Exp.hh" 40 #include "G4ProductionCutsTable.hh" 40 #include "G4ProductionCutsTable.hh" 41 41 42 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 43 43 44 using namespace std; 44 using namespace std; 45 namespace { G4Mutex LivermorePolarizedGammaCon << 46 45 47 G4PhysicsFreeVector* G4LivermorePolarizedGamma << 46 G4int G4LivermorePolarizedGammaConversionModel::maxZ = 99; >> 47 G4LPhysicsFreeVector* G4LivermorePolarizedGammaConversionModel::data[] = {0}; 48 48 49 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 50 50 51 G4LivermorePolarizedGammaConversionModel::G4Li 51 G4LivermorePolarizedGammaConversionModel::G4LivermorePolarizedGammaConversionModel( 52 const G4ParticleDefinition*, const G4String 52 const G4ParticleDefinition*, const G4String& nam) 53 :G4VEmModel(nam), smallEnergy(2.*MeV), isIni << 53 :G4VEmModel(nam), isInitialised(false),smallEnergy(2.*MeV) 54 { 54 { 55 fParticleChange = nullptr; 55 fParticleChange = nullptr; 56 lowEnergyLimit = 2*electron_mass_c2; 56 lowEnergyLimit = 2*electron_mass_c2; 57 57 58 Phi=0.; 58 Phi=0.; 59 Psi=0.; 59 Psi=0.; 60 60 61 verboseLevel= 0; 61 verboseLevel= 0; 62 // Verbosity scale: 62 // Verbosity scale: 63 // 0 = nothing 63 // 0 = nothing 64 // 1 = calculation of cross sections, file o 64 // 1 = calculation of cross sections, file openings, samping of atoms 65 // 2 = entering in methods 65 // 2 = entering in methods 66 66 67 if(verboseLevel > 0) { 67 if(verboseLevel > 0) { 68 G4cout << "Livermore Polarized GammaConver 68 G4cout << "Livermore Polarized GammaConversion is constructed " 69 << G4endl; 69 << G4endl; 70 } 70 } 71 71 72 } 72 } 73 73 74 //....oooOO0OOooo........oooOO0OOooo........oo 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 75 75 76 G4LivermorePolarizedGammaConversionModel::~G4L 76 G4LivermorePolarizedGammaConversionModel::~G4LivermorePolarizedGammaConversionModel() 77 { 77 { 78 if(IsMaster()) { 78 if(IsMaster()) { 79 for(G4int i=0; i<maxZ; ++i) { 79 for(G4int i=0; i<maxZ; ++i) { 80 if(data[i]) { 80 if(data[i]) { 81 delete data[i]; 81 delete data[i]; 82 data[i] = nullptr; << 82 data[i] = 0; 83 } 83 } 84 } 84 } 85 } 85 } 86 } 86 } 87 87 88 //....oooOO0OOooo........oooOO0OOooo........oo 88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 89 89 90 void G4LivermorePolarizedGammaConversionModel: 90 void G4LivermorePolarizedGammaConversionModel::Initialise(const G4ParticleDefinition* particle, 91 const G 91 const G4DataVector& cuts) 92 { 92 { 93 if (verboseLevel > 1) 93 if (verboseLevel > 1) 94 { 94 { 95 G4cout << "Calling1 G4LivermorePolarized 95 G4cout << "Calling1 G4LivermorePolarizedGammaConversionModel::Initialise()" 96 << G4endl 96 << G4endl 97 << "Energy range: " 97 << "Energy range: " 98 << LowEnergyLimit() / MeV << " MeV - " 98 << LowEnergyLimit() / MeV << " MeV - " 99 << HighEnergyLimit() / GeV << " G 99 << HighEnergyLimit() / GeV << " GeV" 100 << G4endl; 100 << G4endl; 101 } 101 } 102 102 >> 103 103 if(IsMaster()) 104 if(IsMaster()) 104 { 105 { 105 // Initialise element selector << 106 >> 107 // Initialise element selector >> 108 106 InitialiseElementSelectors(particle, cut 109 InitialiseElementSelectors(particle, cuts); 107 110 108 // Access to elements 111 // Access to elements 109 const char* path = G4FindDataDir("G4LEDA << 112 >> 113 char* path = getenv("G4LEDATA"); 110 114 111 G4ProductionCutsTable* theCoupleTable = 115 G4ProductionCutsTable* theCoupleTable = 112 G4ProductionCutsTable::GetProductionCutsTabl 116 G4ProductionCutsTable::GetProductionCutsTable(); 113 117 114 G4int numOfCouples = (G4int)theCoupleTab << 118 G4int numOfCouples = theCoupleTable->GetTableSize(); 115 119 116 for(G4int i=0; i<numOfCouples; ++i) 120 for(G4int i=0; i<numOfCouples; ++i) 117 { 121 { 118 const G4Material* material = 122 const G4Material* material = 119 theCoupleTable->GetMaterialCutsCouple(i) 123 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial(); 120 const G4ElementVector* theElementVector = 124 const G4ElementVector* theElementVector = material->GetElementVector(); 121 std::size_t nelm = material->GetNumberOfEl << 125 G4int nelm = material->GetNumberOfElements(); 122 126 123 for (std::size_t j=0; j<nelm; ++j) << 127 for (G4int j=0; j<nelm; ++j) 124 { 128 { 125 G4int Z = (G4int)(*theElementVector)[j 129 G4int Z = (G4int)(*theElementVector)[j]->GetZ(); 126 if(Z < 1) { Z = 1; } 130 if(Z < 1) { Z = 1; } 127 else if(Z > maxZ) { Z = maxZ; } 131 else if(Z > maxZ) { Z = maxZ; } 128 if(!data[Z]) { ReadData(Z, path); } 132 if(!data[Z]) { ReadData(Z, path); } 129 } 133 } 130 } 134 } 131 } 135 } 132 if(isInitialised) { return; } 136 if(isInitialised) { return; } 133 fParticleChange = GetParticleChangeForGamma( 137 fParticleChange = GetParticleChangeForGamma(); 134 isInitialised = true; 138 isInitialised = true; >> 139 135 } 140 } 136 141 >> 142 137 //....oooOO0OOooo........oooOO0OOooo........oo 143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 138 144 139 void G4LivermorePolarizedGammaConversionModel: 145 void G4LivermorePolarizedGammaConversionModel::InitialiseLocal( 140 const G4ParticleDefinition*, G4VEmModel* 146 const G4ParticleDefinition*, G4VEmModel* masterModel) 141 { 147 { 142 SetElementSelectors(masterModel->GetElementS 148 SetElementSelectors(masterModel->GetElementSelectors()); 143 } 149 } 144 150 145 //....oooOO0OOooo........oooOO0OOooo........oo 151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 146 152 147 G4double G4LivermorePolarizedGammaConversionMo 153 G4double G4LivermorePolarizedGammaConversionModel::MinPrimaryEnergy(const G4Material*, 148 const G4ParticleDefinition*, G4doub 154 const G4ParticleDefinition*, G4double) 149 { 155 { 150 return lowEnergyLimit; 156 return lowEnergyLimit; 151 } 157 } 152 158 153 //....oooOO0OOooo........oooOO0OOooo........oo 159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 154 160 155 void G4LivermorePolarizedGammaConversionModel: << 161 void G4LivermorePolarizedGammaConversionModel::ReadData(size_t Z, const char* path) 156 { 162 { 157 if (verboseLevel > 1) 163 if (verboseLevel > 1) 158 { 164 { 159 G4cout << "Calling ReadData() of G4Liver 165 G4cout << "Calling ReadData() of G4LivermorePolarizedGammaConversionModel" 160 << G4endl; 166 << G4endl; 161 } 167 } 162 168 163 if(data[Z]) { return; } 169 if(data[Z]) { return; } 164 170 165 const char* datadir = path; 171 const char* datadir = path; 166 172 167 if(!datadir) 173 if(!datadir) 168 { 174 { 169 datadir = G4FindDataDir("G4LEDATA"); << 175 datadir = getenv("G4LEDATA"); 170 if(!datadir) 176 if(!datadir) 171 { 177 { 172 G4Exception("G4LivermorePolarizedGammaConv 178 G4Exception("G4LivermorePolarizedGammaConversionModel::ReadData()", 173 "em0006",FatalException, 179 "em0006",FatalException, 174 "Environment variable G4LEDATA not d 180 "Environment variable G4LEDATA not defined"); 175 return; 181 return; 176 } 182 } 177 } 183 } 178 // << 184 179 data[Z] = new G4PhysicsFreeVector(0,/*spline << 185 // >> 186 >> 187 data[Z] = new G4LPhysicsFreeVector(); >> 188 180 // 189 // >> 190 181 std::ostringstream ost; 191 std::ostringstream ost; 182 ost << datadir << "/livermore/pair/pp-cs-" < 192 ost << datadir << "/livermore/pair/pp-cs-" << Z <<".dat"; 183 std::ifstream fin(ost.str().c_str()); 193 std::ifstream fin(ost.str().c_str()); 184 194 185 if( !fin.is_open()) 195 if( !fin.is_open()) 186 { 196 { 187 G4ExceptionDescription ed; 197 G4ExceptionDescription ed; 188 ed << "G4LivermorePolarizedGammaConversi 198 ed << "G4LivermorePolarizedGammaConversionModel data file <" << ost.str().c_str() 189 << "> is not opened!" << G4endl; 199 << "> is not opened!" << G4endl; 190 G4Exception("G4LivermorePolarizedGammaCo 200 G4Exception("G4LivermorePolarizedGammaConversionModel::ReadData()", 191 "em0003",FatalException, 201 "em0003",FatalException, 192 ed,"G4LEDATA version should be G4EMLOW6. 202 ed,"G4LEDATA version should be G4EMLOW6.27 or later."); 193 return; 203 return; 194 } 204 } 195 else 205 else 196 { 206 { 197 207 198 if(verboseLevel > 3) { G4cout << "File " 208 if(verboseLevel > 3) { G4cout << "File " << ost.str() 199 << " is opened by G4LivermorePolar 209 << " is opened by G4LivermorePolarizedGammaConversionModel" << G4endl;} 200 210 201 data[Z]->Retrieve(fin, true); 211 data[Z]->Retrieve(fin, true); 202 } 212 } 203 213 204 // Activation of spline interpolation 214 // Activation of spline interpolation 205 data[Z]->FillSecondDerivatives(); << 215 data[Z] ->SetSpline(true); 206 216 207 } 217 } 208 218 209 //....oooOO0OOooo........oooOO0OOooo........oo 219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 210 220 211 G4double G4LivermorePolarizedGammaConversionMo 221 G4double G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom( 212 const G 222 const G4ParticleDefinition*, 213 G4double GammaEnergy, 223 G4double GammaEnergy, 214 G4double Z, G4double, 224 G4double Z, G4double, 215 G4double, G4double) 225 G4double, G4double) 216 { 226 { 217 if (verboseLevel > 1) { 227 if (verboseLevel > 1) { 218 G4cout << "G4LivermorePolarizedGammaConver 228 G4cout << "G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom()" 219 << G4endl; 229 << G4endl; 220 } 230 } 221 if (GammaEnergy < lowEnergyLimit) { return 0 231 if (GammaEnergy < lowEnergyLimit) { return 0.0; } 222 232 223 G4double xs = 0.0; 233 G4double xs = 0.0; 224 234 225 G4int intZ=G4int(Z); 235 G4int intZ=G4int(Z); 226 236 227 if(intZ < 1 || intZ > maxZ) { return xs; } 237 if(intZ < 1 || intZ > maxZ) { return xs; } 228 238 229 G4PhysicsFreeVector* pv = data[intZ]; << 239 G4LPhysicsFreeVector* pv = data[intZ]; 230 240 231 // if element was not initialised 241 // if element was not initialised 232 // do initialisation safely for MT mode 242 // do initialisation safely for MT mode 233 if(!pv) 243 if(!pv) 234 { 244 { 235 InitialiseForElement(0, intZ); 245 InitialiseForElement(0, intZ); 236 pv = data[intZ]; 246 pv = data[intZ]; 237 if(!pv) { return xs; } 247 if(!pv) { return xs; } 238 } 248 } 239 // x-section is taken from the table 249 // x-section is taken from the table 240 xs = pv->Value(GammaEnergy); 250 xs = pv->Value(GammaEnergy); 241 251 242 if(verboseLevel > 0) 252 if(verboseLevel > 0) 243 { 253 { 244 G4int n = G4int(pv->GetVectorLength() - << 254 G4int n = pv->GetVectorLength() - 1; 245 G4cout << "****** DEBUG: tcs value for 255 G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)=" 246 << GammaEnergy/MeV << G4endl; 256 << GammaEnergy/MeV << G4endl; 247 G4cout << " cs (Geant4 internal unit) 257 G4cout << " cs (Geant4 internal unit)=" << xs << G4endl; 248 G4cout << " -> first cs value in EA 258 G4cout << " -> first cs value in EADL data file (iu) =" << (*pv)[0] << G4endl; 249 G4cout << " -> last cs value in EA 259 G4cout << " -> last cs value in EADL data file (iu) =" << (*pv)[n] << G4endl; 250 G4cout << "*************************** 260 G4cout << "*********************************************************" << G4endl; 251 } 261 } 252 262 253 return xs; 263 return xs; 254 } 264 } 255 265 256 //....oooOO0OOooo........oooOO0OOooo........oo 266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 257 267 258 void 268 void 259 G4LivermorePolarizedGammaConversionModel::Samp 269 G4LivermorePolarizedGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, 260 const G4MaterialCutsCouple* 270 const G4MaterialCutsCouple* couple, 261 const G4DynamicParticle* aDy 271 const G4DynamicParticle* aDynamicGamma, 262 G4double, 272 G4double, 263 G4double) 273 G4double) 264 { 274 { 265 275 266 // Fluorescence generated according to: 276 // Fluorescence generated according to: 267 // J. Stepanek ,"A program to determine the 277 // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic 268 // subshell ionisation by a particle or due 278 // subshell ionisation by a particle or due to deexcitation or decay of radionuclides", 269 // Comp. Phys. Comm. 1206 pp 1-1-9 (1997) 279 // Comp. Phys. Comm. 1206 pp 1-1-9 (1997) >> 280 270 if (verboseLevel > 3) 281 if (verboseLevel > 3) 271 G4cout << "Calling SampleSecondaries() of 282 G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedGammaConversionModel" << G4endl; 272 283 273 G4double photonEnergy = aDynamicGamma->GetKi 284 G4double photonEnergy = aDynamicGamma->GetKineticEnergy(); >> 285 // Within energy limit? 274 286 275 if(photonEnergy <= lowEnergyLimit) 287 if(photonEnergy <= lowEnergyLimit) 276 { 288 { 277 fParticleChange->ProposeTrackStatus(fSto 289 fParticleChange->ProposeTrackStatus(fStopAndKill); 278 fParticleChange->SetProposedKineticEnerg 290 fParticleChange->SetProposedKineticEnergy(0.); 279 return; 291 return; 280 } 292 } 281 293 >> 294 282 G4ThreeVector gammaPolarization0 = aDynamicG 295 G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization(); 283 G4ThreeVector gammaDirection0 = aDynamicGamm 296 G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection(); 284 297 285 // Make sure that the polarization vector is 298 // Make sure that the polarization vector is perpendicular to the 286 // gamma direction. If not 299 // gamma direction. If not >> 300 287 if(!(gammaPolarization0.isOrthogonal(gammaDi 301 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0)) 288 { // only for testing now 302 { // only for testing now 289 gammaPolarization0 = GetRandomPolarizati 303 gammaPolarization0 = GetRandomPolarization(gammaDirection0); 290 } 304 } 291 else 305 else 292 { 306 { 293 if ( gammaPolarization0.howOrthogonal(ga 307 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0) 294 { 308 { 295 gammaPolarization0 = GetPerpendicularPolar 309 gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0); 296 } 310 } 297 } 311 } 298 312 299 // End of Protection 313 // End of Protection 300 314 >> 315 301 G4double epsilon ; 316 G4double epsilon ; 302 G4double epsilon0Local = electron_mass_c2 / 317 G4double epsilon0Local = electron_mass_c2 / photonEnergy ; 303 318 304 // Do it fast if photon energy < 2. MeV 319 // Do it fast if photon energy < 2. MeV 305 320 306 if (photonEnergy < smallEnergy ) 321 if (photonEnergy < smallEnergy ) 307 { 322 { 308 epsilon = epsilon0Local + (0.5 - epsilon 323 epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand(); 309 } 324 } 310 else 325 else 311 { 326 { 312 // Select randomly one element in the cu << 327 >> 328 // Select randomly one element in the current material >> 329 313 const G4ParticleDefinition* particle = 330 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition(); 314 const G4Element* element = SelectRandomA 331 const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy); 315 332 316 if (element == nullptr) << 333 >> 334 if (element == 0) 317 { 335 { 318 G4cout << "G4LivermorePolarizedGamma 336 G4cout << "G4LivermorePolarizedGammaConversionModel::SampleSecondaries - element = 0" << G4endl; 319 return; 337 return; 320 } 338 } 321 339 322 340 323 G4IonisParamElm* ionisation = element->G << 341 G4IonisParamElm* ionisation = element->GetIonisation(); 324 if (ionisation == nullptr) << 342 >> 343 if (ionisation == 0) 325 { 344 { 326 G4cout << "G4LivermorePolarizedGamma 345 G4cout << "G4LivermorePolarizedGammaConversionModel::SampleSecondaries - ionisation = 0" << G4endl; 327 return; 346 return; 328 } 347 } 329 348 330 // Extract Coulomb factor for this Eleme 349 // Extract Coulomb factor for this Element >> 350 331 G4double fZ = 8. * (ionisation->GetlogZ3 351 G4double fZ = 8. * (ionisation->GetlogZ3()); 332 if (photonEnergy > 50. * MeV) fZ += 8. * 352 if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb()); 333 353 334 // Limits of the screening variable 354 // Limits of the screening variable 335 G4double screenFactor = 136. * epsilon0L 355 G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ; 336 G4double screenMax = G4Exp ((42.24 - fZ) 356 G4double screenMax = G4Exp ((42.24 - fZ)/8.368) - 0.952 ; 337 G4double screenMin = std::min(4.*screenF 357 G4double screenMin = std::min(4.*screenFactor,screenMax) ; 338 358 339 // Limits of the energy sampling 359 // Limits of the energy sampling 340 G4double epsilon1 = 0.5 - 0.5 * sqrt(1. 360 G4double epsilon1 = 0.5 - 0.5 * sqrt(1. - screenMin / screenMax) ; 341 G4double epsilonMin = std::max(epsilon0L 361 G4double epsilonMin = std::max(epsilon0Local,epsilon1); 342 G4double epsilonRange = 0.5 - epsilonMin 362 G4double epsilonRange = 0.5 - epsilonMin ; 343 363 344 // Sample the energy rate of the created 364 // Sample the energy rate of the created electron (or positron) 345 G4double screen; 365 G4double screen; 346 G4double gReject ; 366 G4double gReject ; 347 367 348 G4double f10 = ScreenFunction1(screenMin 368 G4double f10 = ScreenFunction1(screenMin) - fZ; 349 G4double f20 = ScreenFunction2(screenMin 369 G4double f20 = ScreenFunction2(screenMin) - fZ; 350 G4double normF1 = std::max(f10 * epsilon 370 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.); 351 G4double normF2 = std::max(1.5 * f20,0.) 371 G4double normF2 = std::max(1.5 * f20,0.); 352 372 353 do { 373 do { 354 if (normF1 / (normF1 + normF2) > G4Uni 374 if (normF1 / (normF1 + normF2) > G4UniformRand() ) 355 { 375 { 356 epsilon = 0.5 - epsilonRange * pow 376 epsilon = 0.5 - epsilonRange * pow(G4UniformRand(), 0.3333) ; 357 screen = screenFactor / (epsilon * 377 screen = screenFactor / (epsilon * (1. - epsilon)); 358 gReject = (ScreenFunction1(screen) 378 gReject = (ScreenFunction1(screen) - fZ) / f10 ; 359 } 379 } 360 else 380 else 361 { 381 { 362 epsilon = epsilonMin + epsilonRang 382 epsilon = epsilonMin + epsilonRange * G4UniformRand(); 363 screen = screenFactor / (epsilon * 383 screen = screenFactor / (epsilon * (1 - epsilon)); 364 gReject = (ScreenFunction2(screen) 384 gReject = (ScreenFunction2(screen) - fZ) / f20 ; >> 385 >> 386 365 } 387 } 366 } while ( gReject < G4UniformRand() ); 388 } while ( gReject < G4UniformRand() ); >> 389 367 } // End of epsilon sampling 390 } // End of epsilon sampling 368 391 369 // Fix charges randomly 392 // Fix charges randomly >> 393 370 G4double electronTotEnergy; 394 G4double electronTotEnergy; 371 G4double positronTotEnergy; 395 G4double positronTotEnergy; 372 396 >> 397 >> 398 // if (G4int(2*G4UniformRand())) 373 if (G4UniformRand() > 0.5) 399 if (G4UniformRand() > 0.5) 374 { 400 { 375 electronTotEnergy = (1. - epsilon) * pho 401 electronTotEnergy = (1. - epsilon) * photonEnergy; 376 positronTotEnergy = epsilon * photonEner 402 positronTotEnergy = epsilon * photonEnergy; 377 } 403 } 378 else 404 else 379 { 405 { 380 positronTotEnergy = (1. - epsilon) * pho 406 positronTotEnergy = (1. - epsilon) * photonEnergy; 381 electronTotEnergy = epsilon * photonEner 407 electronTotEnergy = epsilon * photonEnergy; 382 } 408 } 383 409 384 // Scattered electron (positron) angles. ( Z 410 // Scattered electron (positron) angles. ( Z - axis along the parent photon) 385 // Universal distribution suggested by L. Ur 411 // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211), 386 // derived from Tsai distribution (Rev. Mod. 412 // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977) >> 413 >> 414 /* >> 415 G4double u; >> 416 const G4double a1 = 0.625; >> 417 G4double a2 = 3. * a1; >> 418 >> 419 if (0.25 > G4UniformRand()) >> 420 { >> 421 u = - log(G4UniformRand() * G4UniformRand()) / a1 ; >> 422 } >> 423 else >> 424 { >> 425 u = - log(G4UniformRand() * G4UniformRand()) / a2 ; >> 426 } >> 427 */ >> 428 387 G4double Ene = electronTotEnergy/electron_ma 429 G4double Ene = electronTotEnergy/electron_mass_c2; // Normalized energy 388 430 389 G4double cosTheta = 0.; 431 G4double cosTheta = 0.; 390 G4double sinTheta = 0.; 432 G4double sinTheta = 0.; 391 433 392 SetTheta(&cosTheta,&sinTheta,Ene); 434 SetTheta(&cosTheta,&sinTheta,Ene); >> 435 >> 436 // G4double theta = u * electron_mass_c2 / photonEnergy ; >> 437 // G4double phi = twopi * G4UniformRand() ; >> 438 393 G4double phi,psi=0.; 439 G4double phi,psi=0.; 394 440 395 //corrected e+ e- angular angular distributi 441 //corrected e+ e- angular angular distribution //preliminary! >> 442 >> 443 // if(photonEnergy>50*MeV) >> 444 // { 396 phi = SetPhi(photonEnergy); 445 phi = SetPhi(photonEnergy); 397 psi = SetPsi(photonEnergy,phi); 446 psi = SetPsi(photonEnergy,phi); >> 447 // } >> 448 //else >> 449 // { >> 450 //psi = G4UniformRand()*2.*pi; >> 451 //phi = pi; // coplanar >> 452 // } >> 453 398 Psi = psi; 454 Psi = psi; 399 Phi = phi; 455 Phi = phi; >> 456 //G4cout << "PHI " << phi << G4endl; >> 457 //G4cout << "PSI " << psi << G4endl; 400 458 401 G4double phie, phip; 459 G4double phie, phip; 402 G4double choice, choice2; 460 G4double choice, choice2; 403 choice = G4UniformRand(); 461 choice = G4UniformRand(); 404 choice2 = G4UniformRand(); 462 choice2 = G4UniformRand(); 405 463 406 if (choice2 <= 0.5) 464 if (choice2 <= 0.5) 407 { 465 { 408 // do nothing 466 // do nothing 409 // phi = phi; 467 // phi = phi; 410 } 468 } 411 else 469 else 412 { 470 { 413 phi = -phi; 471 phi = -phi; 414 } 472 } 415 473 416 if (choice <= 0.5) 474 if (choice <= 0.5) 417 { 475 { 418 phie = psi; //azimuthal angle for the el 476 phie = psi; //azimuthal angle for the electron 419 phip = phie+phi; //azimuthal angle for t 477 phip = phie+phi; //azimuthal angle for the positron 420 } 478 } 421 else 479 else 422 { 480 { 423 // opzione 1 phie / phip equivalenti 481 // opzione 1 phie / phip equivalenti >> 482 424 phip = psi; //azimuthal angle for the po 483 phip = psi; //azimuthal angle for the positron 425 phie = phip + phi; //azimuthal angle for 484 phie = phip + phi; //azimuthal angle for the electron 426 } 485 } 427 486 428 487 429 // Electron Kinematics 488 // Electron Kinematics >> 489 430 G4double dirX = sinTheta*cos(phie); 490 G4double dirX = sinTheta*cos(phie); 431 G4double dirY = sinTheta*sin(phie); 491 G4double dirY = sinTheta*sin(phie); 432 G4double dirZ = cosTheta; 492 G4double dirZ = cosTheta; 433 G4ThreeVector electronDirection(dirX,dirY,di 493 G4ThreeVector electronDirection(dirX,dirY,dirZ); 434 494 435 // Kinematics of the created pair: 495 // Kinematics of the created pair: 436 // the electron and positron are assumed to 496 // the electron and positron are assumed to have a symetric angular 437 // distribution with respect to the Z axis a 497 // distribution with respect to the Z axis along the parent photon 438 498 >> 499 //G4double localEnergyDeposit = 0. ; >> 500 439 G4double electronKineEnergy = std::max(0.,el 501 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ; 440 502 441 SystemOfRefChange(gammaDirection0,electronDi 503 SystemOfRefChange(gammaDirection0,electronDirection, 442 gammaPolarization0); 504 gammaPolarization0); 443 505 444 G4DynamicParticle* particle1 = new G4Dynamic 506 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(), 445 electronDirection, 507 electronDirection, 446 electronKineEnergy); 508 electronKineEnergy); 447 509 448 // The e+ is always created (even with kinet 510 // The e+ is always created (even with kinetic energy = 0) for further annihilation >> 511 449 Ene = positronTotEnergy/electron_mass_c2; // 512 Ene = positronTotEnergy/electron_mass_c2; // Normalized energy 450 513 451 cosTheta = 0.; 514 cosTheta = 0.; 452 sinTheta = 0.; 515 sinTheta = 0.; 453 516 454 SetTheta(&cosTheta,&sinTheta,Ene); 517 SetTheta(&cosTheta,&sinTheta,Ene); 455 518 456 // Positron Kinematics 519 // Positron Kinematics >> 520 457 dirX = sinTheta*cos(phip); 521 dirX = sinTheta*cos(phip); 458 dirY = sinTheta*sin(phip); 522 dirY = sinTheta*sin(phip); 459 dirZ = cosTheta; 523 dirZ = cosTheta; 460 G4ThreeVector positronDirection(dirX,dirY,di 524 G4ThreeVector positronDirection(dirX,dirY,dirZ); 461 525 462 G4double positronKineEnergy = std::max(0.,po 526 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ; 463 SystemOfRefChange(gammaDirection0,positronDi 527 SystemOfRefChange(gammaDirection0,positronDirection, 464 gammaPolarization0); 528 gammaPolarization0); 465 529 466 // Create G4DynamicParticle object for the p 530 // Create G4DynamicParticle object for the particle2 467 G4DynamicParticle* particle2 = new G4Dynamic 531 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(), 468 532 positronDirection, positronKineEnergy); >> 533 >> 534 469 fvect->push_back(particle1); 535 fvect->push_back(particle1); 470 fvect->push_back(particle2); 536 fvect->push_back(particle2); 471 537 472 // Kill the incident photon 538 // Kill the incident photon 473 fParticleChange->SetProposedKineticEnergy(0. 539 fParticleChange->SetProposedKineticEnergy(0.); 474 fParticleChange->ProposeTrackStatus(fStopAnd 540 fParticleChange->ProposeTrackStatus(fStopAndKill); >> 541 475 } 542 } 476 543 477 //....oooOO0OOooo........oooOO0OOooo........oo 544 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 478 545 479 G4double G4LivermorePolarizedGammaConversionMo 546 G4double G4LivermorePolarizedGammaConversionModel::ScreenFunction1(G4double screenVariable) 480 { 547 { 481 // Compute the value of the screening functi 548 // Compute the value of the screening function 3*phi1 - phi2 >> 549 482 G4double value; 550 G4double value; >> 551 483 if (screenVariable > 1.) 552 if (screenVariable > 1.) 484 value = 42.24 - 8.368 * log(screenVariable 553 value = 42.24 - 8.368 * log(screenVariable + 0.952); 485 else 554 else 486 value = 42.392 - screenVariable * (7.796 - 555 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable); 487 556 488 return value; 557 return value; 489 } 558 } 490 559 491 560 492 561 493 G4double G4LivermorePolarizedGammaConversionMo 562 G4double G4LivermorePolarizedGammaConversionModel::ScreenFunction2(G4double screenVariable) 494 { 563 { 495 // Compute the value of the screening functi 564 // Compute the value of the screening function 1.5*phi1 - 0.5*phi2 >> 565 496 G4double value; 566 G4double value; 497 567 498 if (screenVariable > 1.) 568 if (screenVariable > 1.) 499 value = 42.24 - 8.368 * log(screenVariable 569 value = 42.24 - 8.368 * log(screenVariable + 0.952); 500 else 570 else 501 value = 41.405 - screenVariable * (5.828 - 571 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable); 502 572 503 return value; 573 return value; 504 } 574 } 505 575 506 576 507 void G4LivermorePolarizedGammaConversionModel: 577 void G4LivermorePolarizedGammaConversionModel::SetTheta(G4double* p_cosTheta, G4double* p_sinTheta, G4double Energy) 508 { 578 { >> 579 509 // to avoid computational errors since Theta 580 // to avoid computational errors since Theta could be very small 510 // Energy in Normalized Units (!) 581 // Energy in Normalized Units (!) 511 582 512 G4double Momentum = sqrt(Energy*Energy -1); 583 G4double Momentum = sqrt(Energy*Energy -1); 513 G4double Rand = G4UniformRand(); 584 G4double Rand = G4UniformRand(); 514 585 515 *p_cosTheta = (Energy*((2*Rand)- 1) + Moment 586 *p_cosTheta = (Energy*((2*Rand)- 1) + Momentum)/((Momentum*(2*Rand-1))+Energy); 516 *p_sinTheta = (2*sqrt(Rand*(1-Rand)))/(Momen 587 *p_sinTheta = (2*sqrt(Rand*(1-Rand)))/(Momentum*(2*Rand-1)+Energy); 517 } 588 } 518 589 519 590 520 591 521 G4double G4LivermorePolarizedGammaConversionMo 592 G4double G4LivermorePolarizedGammaConversionModel::SetPhi(G4double Energy) 522 { 593 { >> 594 >> 595 523 G4double value = 0.; 596 G4double value = 0.; 524 G4double Ene = Energy/MeV; 597 G4double Ene = Energy/MeV; 525 598 526 G4double pl[4]; 599 G4double pl[4]; >> 600 >> 601 527 G4double pt[2]; 602 G4double pt[2]; 528 G4double xi = 0; 603 G4double xi = 0; 529 G4double xe = 0.; 604 G4double xe = 0.; 530 G4double n1=0.; 605 G4double n1=0.; 531 G4double n2=0.; 606 G4double n2=0.; 532 607 >> 608 533 if (Ene>=50.) 609 if (Ene>=50.) 534 { 610 { 535 const G4double ay0=5.6, by0=18.6, aa0=2. 611 const G4double ay0=5.6, by0=18.6, aa0=2.9, ba0 = 8.16E-3; 536 const G4double aw = 0.0151, bw = 10.7, c 612 const G4double aw = 0.0151, bw = 10.7, cw = -410.; 537 613 538 const G4double axc = 3.1455, bxc = -1.11 614 const G4double axc = 3.1455, bxc = -1.11, cxc = 310.; 539 615 540 pl[0] = Fln(ay0,by0,Ene); 616 pl[0] = Fln(ay0,by0,Ene); 541 pl[1] = aa0 + ba0*(Ene); 617 pl[1] = aa0 + ba0*(Ene); 542 pl[2] = Poli(aw,bw,cw,Ene); 618 pl[2] = Poli(aw,bw,cw,Ene); 543 pl[3] = Poli(axc,bxc,cxc,Ene); 619 pl[3] = Poli(axc,bxc,cxc,Ene); 544 620 545 const G4double abf = 3.1216, bbf = 2.68; 621 const G4double abf = 3.1216, bbf = 2.68; 546 pt[0] = -1.4; 622 pt[0] = -1.4; 547 pt[1] = abf + bbf/Ene; 623 pt[1] = abf + bbf/Ene; 548 624 >> 625 >> 626 >> 627 //G4cout << "PL > 50. "<< pl[0] << " " << pl[1] << " " << pl[2] << " " <<pl[3] << " " << G4endl; >> 628 549 xi = 3.0; 629 xi = 3.0; 550 xe = Encu(pl,pt,xi); 630 xe = Encu(pl,pt,xi); >> 631 //G4cout << "ENCU "<< xe << G4endl; 551 n1 = Fintlor(pl,pi) - Fintlor(pl,xe); 632 n1 = Fintlor(pl,pi) - Fintlor(pl,xe); 552 n2 = Finttan(pt,xe) - Finttan(pt,0.); 633 n2 = Finttan(pt,xe) - Finttan(pt,0.); 553 } 634 } 554 else 635 else 555 { 636 { 556 const G4double ay0=0.144, by0=0.11; 637 const G4double ay0=0.144, by0=0.11; 557 const G4double aa0=2.7, ba0 = 2.74; 638 const G4double aa0=2.7, ba0 = 2.74; 558 const G4double aw = 0.21, bw = 10.8, cw 639 const G4double aw = 0.21, bw = 10.8, cw = -58.; 559 const G4double axc = 3.17, bxc = -0.87, 640 const G4double axc = 3.17, bxc = -0.87, cxc = -6.; 560 641 561 pl[0] = Fln(ay0, by0, Ene); 642 pl[0] = Fln(ay0, by0, Ene); 562 pl[1] = Fln(aa0, ba0, Ene); 643 pl[1] = Fln(aa0, ba0, Ene); 563 pl[2] = Poli(aw,bw,cw,Ene); 644 pl[2] = Poli(aw,bw,cw,Ene); 564 pl[3] = Poli(axc,bxc,cxc,Ene); 645 pl[3] = Poli(axc,bxc,cxc,Ene); 565 646 >> 647 //G4cout << "PL < 50."<< pl[0] << " " << pl[1] << " " << pl[2] << " " <<pl[3] << " " << G4endl; >> 648 //G4cout << "ENCU "<< xe << G4endl; 566 n1 = Fintlor(pl,pi) - Fintlor(pl,xe); 649 n1 = Fintlor(pl,pi) - Fintlor(pl,xe); >> 650 567 } 651 } 568 652 569 653 570 G4double n=0.; 654 G4double n=0.; 571 n = n1+n2; 655 n = n1+n2; 572 656 573 G4double c1 = 0.; 657 G4double c1 = 0.; 574 c1 = Glor(pl, xe); 658 c1 = Glor(pl, xe); 575 659 >> 660 /* >> 661 G4double xm = 0.; >> 662 xm = Flor(pl,pl[3])*Glor(pl,pl[3]); >> 663 */ >> 664 576 G4double r1,r2,r3; 665 G4double r1,r2,r3; 577 G4double xco=0.; 666 G4double xco=0.; 578 667 579 if (Ene>=50.) 668 if (Ene>=50.) 580 { 669 { 581 r1= G4UniformRand(); 670 r1= G4UniformRand(); 582 if( r1>=n2/n) 671 if( r1>=n2/n) 583 { 672 { 584 do 673 do 585 { 674 { 586 r2 = G4UniformRand(); 675 r2 = G4UniformRand(); 587 value = Finvlor(pl,xe,r2); 676 value = Finvlor(pl,xe,r2); 588 xco = Glor(pl,value)/c1; 677 xco = Glor(pl,value)/c1; 589 r3 = G4UniformRand(); 678 r3 = G4UniformRand(); 590 } while(r3>=xco); 679 } while(r3>=xco); 591 } 680 } 592 else 681 else 593 { 682 { 594 value = Finvtan(pt,n,r1); 683 value = Finvtan(pt,n,r1); 595 } 684 } 596 } 685 } 597 else 686 else 598 { 687 { 599 do 688 do 600 { 689 { 601 r2 = G4UniformRand(); 690 r2 = G4UniformRand(); 602 value = Finvlor(pl,xe,r2); 691 value = Finvlor(pl,xe,r2); 603 xco = Glor(pl,value)/c1; 692 xco = Glor(pl,value)/c1; 604 r3 = G4UniformRand(); 693 r3 = G4UniformRand(); 605 } while(r3>=xco); 694 } while(r3>=xco); 606 } 695 } >> 696 >> 697 // G4cout << "PHI = " <<value << G4endl; 607 return value; 698 return value; 608 } 699 } 609 << 610 //....oooOO0OOooo........oooOO0OOooo........oo << 611 << 612 G4double G4LivermorePolarizedGammaConversionMo 700 G4double G4LivermorePolarizedGammaConversionModel::SetPsi(G4double Energy, G4double PhiLocal) 613 { 701 { >> 702 614 G4double value = 0.; 703 G4double value = 0.; 615 G4double Ene = Energy/MeV; 704 G4double Ene = Energy/MeV; 616 705 617 G4double p0l[4]; 706 G4double p0l[4]; 618 G4double ppml[4]; 707 G4double ppml[4]; 619 G4double p0t[2]; 708 G4double p0t[2]; 620 G4double ppmt[2]; 709 G4double ppmt[2]; 621 710 622 G4double xi = 0.; 711 G4double xi = 0.; 623 G4double xe0 = 0.; 712 G4double xe0 = 0.; 624 G4double xepm = 0.; 713 G4double xepm = 0.; 625 714 626 if (Ene>=50.) 715 if (Ene>=50.) 627 { 716 { 628 const G4double ay00 = 3.4, by00 = 9.8, a 717 const G4double ay00 = 3.4, by00 = 9.8, aa00 = 1.34, ba00 = 5.3; 629 const G4double aw0 = 0.014, bw0 = 9.7, c 718 const G4double aw0 = 0.014, bw0 = 9.7, cw0 = -2.E4; 630 const G4double axc0 = 3.1423, bxc0 = -2. 719 const G4double axc0 = 3.1423, bxc0 = -2.35, cxc0 = 0.; 631 const G4double ay0p = 1.53, by0p = 3.2, 720 const G4double ay0p = 1.53, by0p = 3.2, aap = 0.67, bap = 8.5E-3; 632 const G4double awp = 6.9E-3, bwp = 12.6, 721 const G4double awp = 6.9E-3, bwp = 12.6, cwp = -3.8E4; 633 const G4double axcp = 2.8E-3,bxcp = -3.1 722 const G4double axcp = 2.8E-3,bxcp = -3.133; 634 const G4double abf0 = 3.1213, bbf0 = 2.6 723 const G4double abf0 = 3.1213, bbf0 = 2.61; 635 const G4double abfpm = 3.1231, bbfpm = 2 724 const G4double abfpm = 3.1231, bbfpm = 2.84; 636 725 637 p0l[0] = Fln(ay00, by00, Ene); 726 p0l[0] = Fln(ay00, by00, Ene); 638 p0l[1] = Fln(aa00, ba00, Ene); 727 p0l[1] = Fln(aa00, ba00, Ene); 639 p0l[2] = Poli(aw0, bw0, cw0, Ene); 728 p0l[2] = Poli(aw0, bw0, cw0, Ene); 640 p0l[3] = Poli(axc0, bxc0, cxc0, Ene); 729 p0l[3] = Poli(axc0, bxc0, cxc0, Ene); 641 730 642 ppml[0] = Fln(ay0p, by0p, Ene); 731 ppml[0] = Fln(ay0p, by0p, Ene); 643 ppml[1] = aap + bap*(Ene); 732 ppml[1] = aap + bap*(Ene); 644 ppml[2] = Poli(awp, bwp, cwp, Ene); 733 ppml[2] = Poli(awp, bwp, cwp, Ene); 645 ppml[3] = Fln(axcp,bxcp,Ene); 734 ppml[3] = Fln(axcp,bxcp,Ene); 646 735 647 p0t[0] = -0.81; 736 p0t[0] = -0.81; 648 p0t[1] = abf0 + bbf0/Ene; 737 p0t[1] = abf0 + bbf0/Ene; 649 ppmt[0] = -0.6; 738 ppmt[0] = -0.6; 650 ppmt[1] = abfpm + bbfpm/Ene; 739 ppmt[1] = abfpm + bbfpm/Ene; 651 740 >> 741 //G4cout << "P0L > 50"<< p0l[0] << " " << p0l[1] << " " << p0l[2] << " " <<p0l[3] << " " << G4endl; >> 742 //G4cout << "PPML > 50"<< ppml[0] << " " << ppml[1] << " " << ppml[2] << " " <<ppml[3] << " " << G4endl; >> 743 652 xi = 3.0; 744 xi = 3.0; 653 xe0 = Encu(p0l, p0t, xi); 745 xe0 = Encu(p0l, p0t, xi); >> 746 //G4cout << "ENCU1 "<< xe0 << G4endl; 654 xepm = Encu(ppml, ppmt, xi); 747 xepm = Encu(ppml, ppmt, xi); >> 748 //G4cout << "ENCU2 "<< xepm << G4endl; 655 } 749 } 656 else 750 else 657 { 751 { 658 const G4double ay00 = 2.82, by00 = 6.35; 752 const G4double ay00 = 2.82, by00 = 6.35; 659 const G4double aa00 = -1.75, ba00 = 0.25 753 const G4double aa00 = -1.75, ba00 = 0.25; 660 754 661 const G4double aw0 = 0.028, bw0 = 5., cw 755 const G4double aw0 = 0.028, bw0 = 5., cw0 = -50.; 662 const G4double axc0 = 3.14213, bxc0 = -2 756 const G4double axc0 = 3.14213, bxc0 = -2.3, cxc0 = 5.7; 663 const G4double ay0p = 1.56, by0p = 3.6; 757 const G4double ay0p = 1.56, by0p = 3.6; 664 const G4double aap = 0.86, bap = 8.3E-3; 758 const G4double aap = 0.86, bap = 8.3E-3; 665 const G4double awp = 0.022, bwp = 7.4, c 759 const G4double awp = 0.022, bwp = 7.4, cwp = -51.; 666 const G4double xcp = 3.1486; 760 const G4double xcp = 3.1486; 667 761 668 p0l[0] = Fln(ay00, by00, Ene); 762 p0l[0] = Fln(ay00, by00, Ene); 669 p0l[1] = aa00+pow(Ene, ba00); 763 p0l[1] = aa00+pow(Ene, ba00); 670 p0l[2] = Poli(aw0, bw0, cw0, Ene); 764 p0l[2] = Poli(aw0, bw0, cw0, Ene); 671 p0l[3] = Poli(axc0, bxc0, cxc0, Ene); 765 p0l[3] = Poli(axc0, bxc0, cxc0, Ene); 672 ppml[0] = Fln(ay0p, by0p, Ene); 766 ppml[0] = Fln(ay0p, by0p, Ene); 673 ppml[1] = aap + bap*(Ene); 767 ppml[1] = aap + bap*(Ene); 674 ppml[2] = Poli(awp, bwp, cwp, Ene); 768 ppml[2] = Poli(awp, bwp, cwp, Ene); 675 ppml[3] = xcp; 769 ppml[3] = xcp; >> 770 676 } 771 } 677 772 678 G4double a,b=0.; 773 G4double a,b=0.; 679 774 680 if (Ene>=50.) 775 if (Ene>=50.) 681 { 776 { 682 if (PhiLocal>xepm) 777 if (PhiLocal>xepm) 683 { 778 { 684 b = (ppml[0]+2*ppml[1]*ppml[2]*Flor( 779 b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml,PhiLocal)); 685 } 780 } 686 else 781 else 687 { 782 { 688 b = Ftan(ppmt,PhiLocal); 783 b = Ftan(ppmt,PhiLocal); 689 } 784 } 690 if (PhiLocal>xe0) 785 if (PhiLocal>xe0) 691 { 786 { 692 a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l 787 a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,PhiLocal)); 693 } 788 } 694 else 789 else 695 { 790 { 696 a = Ftan(p0t,PhiLocal); 791 a = Ftan(p0t,PhiLocal); 697 } 792 } 698 } 793 } 699 else 794 else 700 { 795 { 701 b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml 796 b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml,PhiLocal)); 702 a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,Phi 797 a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,PhiLocal)); 703 } 798 } 704 G4double nr =0.; 799 G4double nr =0.; 705 800 706 if (b>a) 801 if (b>a) 707 { 802 { 708 nr = 1./b; 803 nr = 1./b; 709 } 804 } 710 else 805 else 711 { 806 { 712 nr = 1./a; 807 nr = 1./a; 713 } 808 } 714 809 715 G4double r1,r2=0.; 810 G4double r1,r2=0.; 716 G4double r3 =-1.; 811 G4double r3 =-1.; 717 do 812 do 718 { 813 { 719 r1 = G4UniformRand(); 814 r1 = G4UniformRand(); 720 r2 = G4UniformRand(); 815 r2 = G4UniformRand(); 721 //value = r2*pi; 816 //value = r2*pi; 722 value = 2.*r2*pi; 817 value = 2.*r2*pi; 723 r3 = nr*(a*cos(value)*cos(value) + b*sin 818 r3 = nr*(a*cos(value)*cos(value) + b*sin(value)*sin(value)); 724 }while(r1>r3); 819 }while(r1>r3); 725 820 726 return value; 821 return value; 727 } 822 } 728 823 729 //....oooOO0OOooo........oooOO0OOooo........oo << 730 824 731 G4double G4LivermorePolarizedGammaConversionMo 825 G4double G4LivermorePolarizedGammaConversionModel::Poli 732 (G4double a, G4double b, G4double c, G4double 826 (G4double a, G4double b, G4double c, G4double x) 733 { 827 { 734 G4double value=0.; 828 G4double value=0.; 735 if(x>0.) 829 if(x>0.) 736 { 830 { 737 value =(a + b/x + c/(x*x*x)); 831 value =(a + b/x + c/(x*x*x)); 738 } 832 } 739 else 833 else 740 { 834 { 741 //G4cout << "ERROR in Poli! " << G4endl; 835 //G4cout << "ERROR in Poli! " << G4endl; 742 } 836 } 743 return value; 837 return value; 744 } 838 } 745 << 746 //....oooOO0OOooo........oooOO0OOooo........oo << 747 << 748 G4double G4LivermorePolarizedGammaConversionMo 839 G4double G4LivermorePolarizedGammaConversionModel::Fln 749 (G4double a, G4double b, G4double x) 840 (G4double a, G4double b, G4double x) 750 { 841 { 751 G4double value=0.; 842 G4double value=0.; 752 if(x>0.) 843 if(x>0.) 753 { 844 { 754 value =(a*log(x)-b); 845 value =(a*log(x)-b); 755 } 846 } 756 else 847 else 757 { 848 { 758 //G4cout << "ERROR in Fln! " << G4endl; 849 //G4cout << "ERROR in Fln! " << G4endl; 759 } 850 } 760 return value; 851 return value; 761 } 852 } 762 853 763 //....oooOO0OOooo........oooOO0OOooo........oo << 764 854 765 G4double G4LivermorePolarizedGammaConversionMo 855 G4double G4LivermorePolarizedGammaConversionModel::Encu 766 (G4double* p_p1, G4double* p_p2, G4double x0) 856 (G4double* p_p1, G4double* p_p2, G4double x0) 767 { 857 { 768 G4int i=0; 858 G4int i=0; 769 G4double fx = 1.; 859 G4double fx = 1.; 770 G4double x = x0; 860 G4double x = x0; 771 const G4double xmax = 3.0; 861 const G4double xmax = 3.0; 772 862 773 for(i=0; i<100; ++i) 863 for(i=0; i<100; ++i) 774 { 864 { 775 fx = (Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p 865 fx = (Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p_p2, x))/ 776 (Fdlor(p_p1,x) - Fdtan(p_p2,x)); 866 (Fdlor(p_p1,x) - Fdtan(p_p2,x)); 777 x -= fx; 867 x -= fx; 778 if(x > xmax) { return xmax; } 868 if(x > xmax) { return xmax; } >> 869 // x -= (Flor(p_p1, x)*Glor(p_p1,x) - Ftan(p_p2, x))/ >> 870 // (Fdlor(p_p1,x) - Fdtan(p_p2,x)); >> 871 // fx = Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p_p2, x); >> 872 // G4cout << std::fabs(fx) << " " << i << " " << x << "dentro ENCU " << G4endl; 779 if(std::fabs(fx) <= x*1.0e-6) { break; } 873 if(std::fabs(fx) <= x*1.0e-6) { break; } 780 } 874 } 781 875 782 if(x < 0.0) { x = 0.0; } 876 if(x < 0.0) { x = 0.0; } 783 return x; 877 return x; 784 } 878 } 785 879 786 //....oooOO0OOooo........oooOO0OOooo........oo << 787 880 788 G4double G4LivermorePolarizedGammaConversionMo 881 G4double G4LivermorePolarizedGammaConversionModel::Flor(G4double* p_p1, G4double x) 789 { 882 { 790 G4double value =0.; 883 G4double value =0.; >> 884 // G4double y0 = p_p1[0]; >> 885 // G4double A = p_p1[1]; 791 G4double w = p_p1[2]; 886 G4double w = p_p1[2]; 792 G4double xc = p_p1[3]; 887 G4double xc = p_p1[3]; 793 888 794 value = 1./(pi*(w*w + 4.*(x-xc)*(x-xc))); 889 value = 1./(pi*(w*w + 4.*(x-xc)*(x-xc))); 795 return value; 890 return value; 796 } 891 } 797 892 798 //....oooOO0OOooo........oooOO0OOooo........oo << 799 893 800 G4double G4LivermorePolarizedGammaConversionMo 894 G4double G4LivermorePolarizedGammaConversionModel::Glor(G4double* p_p1, G4double x) 801 { 895 { 802 G4double value =0.; 896 G4double value =0.; 803 G4double y0 = p_p1[0]; 897 G4double y0 = p_p1[0]; 804 G4double A = p_p1[1]; 898 G4double A = p_p1[1]; 805 G4double w = p_p1[2]; 899 G4double w = p_p1[2]; 806 G4double xc = p_p1[3]; 900 G4double xc = p_p1[3]; 807 901 808 value = (y0 *pi*(w*w + 4.*(x-xc)*(x-xc)) + 902 value = (y0 *pi*(w*w + 4.*(x-xc)*(x-xc)) + 2.*A*w); 809 return value; 903 return value; 810 } 904 } 811 905 812 //....oooOO0OOooo........oooOO0OOooo........oo << 813 906 814 G4double G4LivermorePolarizedGammaConversionMo 907 G4double G4LivermorePolarizedGammaConversionModel::Fdlor(G4double* p_p1, G4double x) 815 { 908 { 816 G4double value =0.; 909 G4double value =0.; >> 910 //G4double y0 = p_p1[0]; 817 G4double A = p_p1[1]; 911 G4double A = p_p1[1]; 818 G4double w = p_p1[2]; 912 G4double w = p_p1[2]; 819 G4double xc = p_p1[3]; 913 G4double xc = p_p1[3]; 820 914 821 value = (-16.*A*w*(x-xc))/ 915 value = (-16.*A*w*(x-xc))/ 822 (pi*(w*w+4.*(x-xc)*(x-xc))*(w*w+4.*(x-xc)* 916 (pi*(w*w+4.*(x-xc)*(x-xc))*(w*w+4.*(x-xc)*(x-xc))); 823 return value; 917 return value; 824 } 918 } 825 919 826 //....oooOO0OOooo........oooOO0OOooo........oo << 827 920 828 G4double G4LivermorePolarizedGammaConversionMo 921 G4double G4LivermorePolarizedGammaConversionModel::Fintlor(G4double* p_p1, G4double x) 829 { 922 { 830 G4double value =0.; 923 G4double value =0.; 831 G4double y0 = p_p1[0]; 924 G4double y0 = p_p1[0]; 832 G4double A = p_p1[1]; 925 G4double A = p_p1[1]; 833 G4double w = p_p1[2]; 926 G4double w = p_p1[2]; 834 G4double xc = p_p1[3]; 927 G4double xc = p_p1[3]; 835 928 836 value = y0*x + A*atan( 2*(x-xc)/w) / pi; 929 value = y0*x + A*atan( 2*(x-xc)/w) / pi; 837 return value; 930 return value; 838 } 931 } 839 932 840 933 841 G4double G4LivermorePolarizedGammaConversionMo 934 G4double G4LivermorePolarizedGammaConversionModel::Finvlor(G4double* p_p1, G4double x, G4double r) 842 { 935 { 843 G4double value = 0.; 936 G4double value = 0.; 844 G4double nor = 0.; 937 G4double nor = 0.; >> 938 //G4double y0 = p_p1[0]; >> 939 // G4double A = p_p1[1]; 845 G4double w = p_p1[2]; 940 G4double w = p_p1[2]; 846 G4double xc = p_p1[3]; 941 G4double xc = p_p1[3]; 847 942 848 nor = atan(2.*(pi-xc)/w)/(2.*pi*w) - atan(2. 943 nor = atan(2.*(pi-xc)/w)/(2.*pi*w) - atan(2.*(x-xc)/w)/(2.*pi*w); 849 value = xc - (w/2.)*tan(-2.*r*nor*pi*w+atan( 944 value = xc - (w/2.)*tan(-2.*r*nor*pi*w+atan(2*(xc-x)/w)); 850 945 851 return value; 946 return value; 852 } 947 } 853 948 854 //....oooOO0OOooo........oooOO0OOooo........oo << 855 949 856 G4double G4LivermorePolarizedGammaConversionMo 950 G4double G4LivermorePolarizedGammaConversionModel::Ftan(G4double* p_p1, G4double x) 857 { 951 { 858 G4double value =0.; 952 G4double value =0.; 859 G4double a = p_p1[0]; 953 G4double a = p_p1[0]; 860 G4double b = p_p1[1]; 954 G4double b = p_p1[1]; 861 955 862 value = a /(x-b); 956 value = a /(x-b); 863 return value; 957 return value; 864 } 958 } 865 959 866 //....oooOO0OOooo........oooOO0OOooo........oo << 867 960 868 G4double G4LivermorePolarizedGammaConversionMo 961 G4double G4LivermorePolarizedGammaConversionModel::Fdtan(G4double* p_p1, G4double x) 869 { 962 { 870 G4double value =0.; 963 G4double value =0.; 871 G4double a = p_p1[0]; 964 G4double a = p_p1[0]; 872 G4double b = p_p1[1]; 965 G4double b = p_p1[1]; 873 966 874 value = -1.*a / ((x-b)*(x-b)); 967 value = -1.*a / ((x-b)*(x-b)); 875 return value; 968 return value; 876 } 969 } 877 970 878 //....oooOO0OOooo........oooOO0OOooo........oo << 879 971 880 G4double G4LivermorePolarizedGammaConversionMo 972 G4double G4LivermorePolarizedGammaConversionModel::Finttan(G4double* p_p1, G4double x) 881 { 973 { 882 G4double value =0.; 974 G4double value =0.; 883 G4double a = p_p1[0]; 975 G4double a = p_p1[0]; 884 G4double b = p_p1[1]; 976 G4double b = p_p1[1]; 885 977 >> 978 886 value = a*log(b-x); 979 value = a*log(b-x); 887 return value; 980 return value; 888 } 981 } 889 982 890 //....oooOO0OOooo........oooOO0OOooo........oo << 891 983 892 G4double G4LivermorePolarizedGammaConversionMo 984 G4double G4LivermorePolarizedGammaConversionModel::Finvtan(G4double* p_p1, G4double cnor, G4double r) 893 { 985 { 894 G4double value =0.; 986 G4double value =0.; 895 G4double a = p_p1[0]; 987 G4double a = p_p1[0]; 896 G4double b = p_p1[1]; 988 G4double b = p_p1[1]; 897 989 898 value = b*(1-G4Exp(r*cnor/a)); 990 value = b*(1-G4Exp(r*cnor/a)); 899 991 900 return value; 992 return value; 901 } 993 } 902 994 >> 995 >> 996 >> 997 903 //....oooOO0OOooo........oooOO0OOooo........oo 998 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 904 999 905 G4ThreeVector G4LivermorePolarizedGammaConvers 1000 G4ThreeVector G4LivermorePolarizedGammaConversionModel::SetPerpendicularVector(G4ThreeVector& a) 906 { 1001 { 907 G4double dx = a.x(); 1002 G4double dx = a.x(); 908 G4double dy = a.y(); 1003 G4double dy = a.y(); 909 G4double dz = a.z(); 1004 G4double dz = a.z(); 910 G4double x = dx < 0.0 ? -dx : dx; 1005 G4double x = dx < 0.0 ? -dx : dx; 911 G4double y = dy < 0.0 ? -dy : dy; 1006 G4double y = dy < 0.0 ? -dy : dy; 912 G4double z = dz < 0.0 ? -dz : dz; 1007 G4double z = dz < 0.0 ? -dz : dz; 913 if (x < y) { 1008 if (x < y) { 914 return x < z ? G4ThreeVector(-dy,dx,0) : G 1009 return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy); 915 }else{ 1010 }else{ 916 return y < z ? G4ThreeVector(dz,0,-dx) : G 1011 return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0); 917 } 1012 } 918 } 1013 } 919 1014 920 //....oooOO0OOooo........oooOO0OOooo........oo 1015 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 921 1016 922 G4ThreeVector G4LivermorePolarizedGammaConvers 1017 G4ThreeVector G4LivermorePolarizedGammaConversionModel::GetRandomPolarization(G4ThreeVector& direction0) 923 { 1018 { 924 G4ThreeVector d0 = direction0.unit(); 1019 G4ThreeVector d0 = direction0.unit(); 925 G4ThreeVector a1 = SetPerpendicularVector(d0 1020 G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal 926 G4ThreeVector a0 = a1.unit(); // unit vector 1021 G4ThreeVector a0 = a1.unit(); // unit vector 927 1022 928 G4double rand1 = G4UniformRand(); 1023 G4double rand1 = G4UniformRand(); 929 1024 930 G4double angle = twopi*rand1; // random pola 1025 G4double angle = twopi*rand1; // random polar angle 931 G4ThreeVector b0 = d0.cross(a0); // cross pr 1026 G4ThreeVector b0 = d0.cross(a0); // cross product 932 1027 933 G4ThreeVector c; 1028 G4ThreeVector c; 934 1029 935 c.setX(std::cos(angle)*(a0.x())+std::sin(ang 1030 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x()); 936 c.setY(std::cos(angle)*(a0.y())+std::sin(ang 1031 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y()); 937 c.setZ(std::cos(angle)*(a0.z())+std::sin(ang 1032 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z()); 938 1033 939 G4ThreeVector c0 = c.unit(); 1034 G4ThreeVector c0 = c.unit(); 940 1035 941 return c0; << 1036 return c0; >> 1037 942 } 1038 } 943 1039 944 //....oooOO0OOooo........oooOO0OOooo........oo 1040 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 945 1041 946 G4ThreeVector G4LivermorePolarizedGammaConvers 1042 G4ThreeVector G4LivermorePolarizedGammaConversionModel::GetPerpendicularPolarization 947 (const G4ThreeVector& gammaDirection, const G4 1043 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const 948 { 1044 { >> 1045 949 // 1046 // 950 // The polarization of a photon is always pe 1047 // The polarization of a photon is always perpendicular to its momentum direction. 951 // Therefore this function removes those vec 1048 // Therefore this function removes those vector component of gammaPolarization, which 952 // points in direction of gammaDirection 1049 // points in direction of gammaDirection 953 // 1050 // 954 // Mathematically we search the projection o 1051 // Mathematically we search the projection of the vector a on the plane E, where n is the 955 // plains normal vector. 1052 // plains normal vector. 956 // The basic equation can be found in each g 1053 // The basic equation can be found in each geometry book (e.g. Bronstein): 957 // p = a - (a o n)/(n o n)*n 1054 // p = a - (a o n)/(n o n)*n 958 1055 959 return gammaPolarization - gammaPolarization 1056 return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection; 960 } 1057 } 961 1058 962 //....oooOO0OOooo........oooOO0OOooo........oo 1059 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 963 1060 >> 1061 964 void G4LivermorePolarizedGammaConversionModel: 1062 void G4LivermorePolarizedGammaConversionModel::SystemOfRefChange 965 (G4ThreeVector& direction0,G4ThreeVector& 1063 (G4ThreeVector& direction0,G4ThreeVector& direction1, 966 G4ThreeVector& polarization0) 1064 G4ThreeVector& polarization0) 967 { 1065 { 968 // direction0 is the original photon directi 1066 // direction0 is the original photon direction ---> z 969 // polarization0 is the original photon pola 1067 // polarization0 is the original photon polarization ---> x 970 // need to specify y axis in the real refere 1068 // need to specify y axis in the real reference frame ---> y 971 G4ThreeVector Axis_Z0 = direction0.unit(); 1069 G4ThreeVector Axis_Z0 = direction0.unit(); 972 G4ThreeVector Axis_X0 = polarization0.unit() 1070 G4ThreeVector Axis_X0 = polarization0.unit(); 973 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_ 1071 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed; 974 1072 975 G4double direction_x = direction1.getX(); 1073 G4double direction_x = direction1.getX(); 976 G4double direction_y = direction1.getY(); 1074 G4double direction_y = direction1.getY(); 977 G4double direction_z = direction1.getZ(); 1075 G4double direction_z = direction1.getZ(); 978 1076 979 direction1 = (direction_x*Axis_X0 + directio << 1077 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit(); >> 1078 980 } 1079 } 981 1080 >> 1081 982 //....oooOO0OOooo........oooOO0OOooo........oo 1082 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 983 1083 >> 1084 #include "G4AutoLock.hh" >> 1085 namespace { G4Mutex LivermorePolarizedGammaConversionModelMutex = G4MUTEX_INITIALIZER; } >> 1086 984 void G4LivermorePolarizedGammaConversionModel: 1087 void G4LivermorePolarizedGammaConversionModel::InitialiseForElement( 985 const G4ParticleDefiniti 1088 const G4ParticleDefinition*, 986 G4int Z) 1089 G4int Z) 987 { 1090 { 988 G4AutoLock l(&LivermorePolarizedGammaConvers 1091 G4AutoLock l(&LivermorePolarizedGammaConversionModelMutex); >> 1092 // G4cout << "G4LivermorePolarizedGammaConversionModel::InitialiseForElement Z= " >> 1093 // << Z << G4endl; 989 if(!data[Z]) { ReadData(Z); } 1094 if(!data[Z]) { ReadData(Z); } 990 l.unlock(); 1095 l.unlock(); 991 } 1096 } 992 1097