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 /* 26 /* 27 Authors: 27 Authors: 28 M. Omer and R. Hajima on 15 November 2019 28 M. Omer and R. Hajima on 15 November 2019 29 contact: 29 contact: 30 omer.mohamed@jaea.go.jp and hajima.ryoichi@q 30 omer.mohamed@jaea.go.jp and hajima.ryoichi@qst.go.jp 31 Publication Information: 31 Publication Information: 32 1- M. Omer, R. Hajima, Validating polarizati 32 1- M. Omer, R. Hajima, Validating polarization effects in gamma-rays elastic scattering by Monte 33 Carlo simulation, New J. Phys., vol. 21, 201 33 Carlo simulation, New J. Phys., vol. 21, 2019, pp. 113006 (1-10), 34 https://doi.org/10.1088/1367-2630/ab4d8a 34 https://doi.org/10.1088/1367-2630/ab4d8a 35 */ 35 */ 36 36 37 #include "G4JAEAPolarizedElasticScatteringMode 37 #include "G4JAEAPolarizedElasticScatteringModel.hh" 38 #include "G4SystemOfUnits.hh" 38 #include "G4SystemOfUnits.hh" 39 #include "G4AutoLock.hh" 39 #include "G4AutoLock.hh" 40 namespace { G4Mutex G4JAEAPolarizedElasticScat 40 namespace { G4Mutex G4JAEAPolarizedElasticScatteringModelMutex = G4MUTEX_INITIALIZER; } 41 using namespace std; 41 using namespace std; 42 42 43 //....oooOO0OOooo........oooOO0OOooo........oo 43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 44 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 45 45 46 G4PhysicsFreeVector* G4JAEAPolarizedElasticSca 46 G4PhysicsFreeVector* G4JAEAPolarizedElasticScatteringModel::dataCS[] = {nullptr}; 47 G4DataVector* G4JAEAPolarizedElasticScattering 47 G4DataVector* G4JAEAPolarizedElasticScatteringModel::Polarized_ES_Data[] = {nullptr}; 48 48 49 G4JAEAPolarizedElasticScatteringModel::G4JAEAP 49 G4JAEAPolarizedElasticScatteringModel::G4JAEAPolarizedElasticScatteringModel() 50 :G4VEmModel("G4JAEAPolarizedElasticScatterin 50 :G4VEmModel("G4JAEAPolarizedElasticScatteringModel"),isInitialised(false) 51 { 51 { 52 fParticleChange = 0; 52 fParticleChange = 0; 53 lowEnergyLimit = 100 * keV; //low energy 53 lowEnergyLimit = 100 * keV; //low energy limit for JAEAElasticScattering cross section data 54 fLinearPolarizationSensitvity1=1; 54 fLinearPolarizationSensitvity1=1; 55 fLinearPolarizationSensitvity2=1; 55 fLinearPolarizationSensitvity2=1; 56 fCircularPolarizationSensitvity=1; 56 fCircularPolarizationSensitvity=1; 57 57 58 verboseLevel= 0; 58 verboseLevel= 0; 59 // Verbosity scale for debugging purposes: 59 // Verbosity scale for debugging purposes: 60 // 0 = nothing 60 // 0 = nothing 61 // 1 = calculation of cross sections, file o 61 // 1 = calculation of cross sections, file openings... 62 // 2 = entering in methods 62 // 2 = entering in methods 63 63 64 if(verboseLevel > 0) 64 if(verboseLevel > 0) 65 { 65 { 66 G4cout << "G4JAEAPolarizedElasticScatter 66 G4cout << "G4JAEAPolarizedElasticScatteringModel is constructed " << G4endl; 67 } 67 } 68 68 69 } 69 } 70 70 71 //....oooOO0OOooo........oooOO0OOooo........oo 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 72 72 73 G4JAEAPolarizedElasticScatteringModel::~G4JAEA 73 G4JAEAPolarizedElasticScatteringModel::~G4JAEAPolarizedElasticScatteringModel() 74 { 74 { 75 if(IsMaster()) { 75 if(IsMaster()) { 76 for(G4int i=0; i<=maxZ; ++i) { 76 for(G4int i=0; i<=maxZ; ++i) { 77 if(dataCS[i]) { 77 if(dataCS[i]) { 78 delete dataCS[i]; 78 delete dataCS[i]; 79 dataCS[i] = nullptr; 79 dataCS[i] = nullptr; 80 } 80 } 81 if (Polarized_ES_Data[i]){ 81 if (Polarized_ES_Data[i]){ 82 delete Polarized_ES_Data[i]; 82 delete Polarized_ES_Data[i]; 83 Polarized_ES_Data[i] = nullptr; 83 Polarized_ES_Data[i] = nullptr; 84 } 84 } 85 } 85 } 86 } 86 } 87 } 87 } 88 88 89 //....oooOO0OOooo........oooOO0OOooo........oo 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 90 90 91 void G4JAEAPolarizedElasticScatteringModel::In 91 void G4JAEAPolarizedElasticScatteringModel::Initialise(const G4ParticleDefinition* particle, 92 const G4DataVector& cuts) 92 const G4DataVector& cuts) 93 { 93 { 94 if (verboseLevel > 1) 94 if (verboseLevel > 1) 95 { 95 { 96 G4cout << "Calling Initialise() of G4JAE 96 G4cout << "Calling Initialise() of G4JAEAPolarizedElasticScatteringModel." << G4endl 97 << "Energy range: " 97 << "Energy range: " 98 << LowEnergyLimit() / eV << " eV - " 98 << LowEnergyLimit() / eV << " eV - " 99 << HighEnergyLimit() / GeV << " GeV" 99 << HighEnergyLimit() / GeV << " GeV" 100 << G4endl; 100 << G4endl; 101 } 101 } 102 102 103 if(IsMaster()) { 103 if(IsMaster()) { 104 104 105 // Initialise element selector 105 // Initialise element selector 106 InitialiseElementSelectors(particle, cuts) 106 InitialiseElementSelectors(particle, cuts); 107 107 108 // Access to elements 108 // Access to elements 109 const char* path = G4FindDataDir("G4LEDATA 109 const char* path = G4FindDataDir("G4LEDATA"); 110 G4ProductionCutsTable* theCoupleTable = 110 G4ProductionCutsTable* theCoupleTable = 111 G4ProductionCutsTable::GetProductionCuts 111 G4ProductionCutsTable::GetProductionCutsTable(); 112 G4int numOfCouples = (G4int)theCoupleTable 112 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize(); 113 113 114 for(G4int i=0; i<numOfCouples; ++i) 114 for(G4int i=0; i<numOfCouples; ++i) 115 { 115 { 116 const G4MaterialCutsCouple* couple = 116 const G4MaterialCutsCouple* couple = 117 theCoupleTable->GetMaterialCutsCouple(i); 117 theCoupleTable->GetMaterialCutsCouple(i); 118 const G4Material* material = couple->GetMate 118 const G4Material* material = couple->GetMaterial(); 119 const G4ElementVector* theElementVector = ma 119 const G4ElementVector* theElementVector = material->GetElementVector(); 120 std::size_t nelm = material->GetNumberOfElem 120 std::size_t nelm = material->GetNumberOfElements(); 121 121 122 for (std::size_t j=0; j<nelm; ++j) 122 for (std::size_t j=0; j<nelm; ++j) 123 { 123 { 124 G4int Z = G4lrint((*theElementVector)[j] 124 G4int Z = G4lrint((*theElementVector)[j]->GetZ()); 125 if(Z < 1) { Z = 1; } 125 if(Z < 1) { Z = 1; } 126 else if(Z > maxZ) { Z = maxZ; } 126 else if(Z > maxZ) { Z = maxZ; } 127 if( (!dataCS[Z]) ) { ReadData(Z, path); 127 if( (!dataCS[Z]) ) { ReadData(Z, path); } 128 } 128 } 129 } 129 } 130 } 130 } 131 131 132 if(isInitialised) { return; } 132 if(isInitialised) { return; } 133 fParticleChange = GetParticleChangeForGamma( 133 fParticleChange = GetParticleChangeForGamma(); 134 isInitialised = true; 134 isInitialised = true; 135 135 136 } 136 } 137 137 138 //....oooOO0OOooo........oooOO0OOooo........oo 138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 139 139 140 void G4JAEAPolarizedElasticScatteringModel::In 140 void G4JAEAPolarizedElasticScatteringModel::InitialiseLocal(const G4ParticleDefinition*, 141 G4VEmModel* masterModel) 141 G4VEmModel* masterModel) 142 { 142 { 143 SetElementSelectors(masterModel->GetElementS 143 SetElementSelectors(masterModel->GetElementSelectors()); 144 } 144 } 145 145 146 //....oooOO0OOooo........oooOO0OOooo........oo 146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 147 147 148 void G4JAEAPolarizedElasticScatteringModel::Re 148 void G4JAEAPolarizedElasticScatteringModel::ReadData(std::size_t Z, const char* path) 149 { 149 { 150 if (verboseLevel > 1) 150 if (verboseLevel > 1) 151 { 151 { 152 G4cout << "Calling ReadData() of G4JAEAP 152 G4cout << "Calling ReadData() of G4JAEAPolarizedElasticScatteringModel" 153 << G4endl; 153 << G4endl; 154 } 154 } 155 155 156 if(dataCS[Z]) { return; } 156 if(dataCS[Z]) { return; } 157 157 158 const char* datadir = path; 158 const char* datadir = path; 159 if(!datadir) 159 if(!datadir) 160 { 160 { 161 datadir = G4FindDataDir("G4LEDATA"); 161 datadir = G4FindDataDir("G4LEDATA"); 162 if(!datadir) 162 if(!datadir) 163 { 163 { 164 G4Exception("G4JAEAPolarizedElasticScatter 164 G4Exception("G4JAEAPolarizedElasticScatteringModel::ReadData()","em0006", 165 FatalException, 165 FatalException, 166 "Environment variable G4LEDATA not d 166 "Environment variable G4LEDATA not defined"); 167 return; 167 return; 168 } 168 } 169 } 169 } 170 170 171 std::ostringstream ostCS; 171 std::ostringstream ostCS; 172 ostCS << datadir << "/JAEAESData/amp_Z_" << 172 ostCS << datadir << "/JAEAESData/amp_Z_" << Z ; 173 std::ifstream ES_Data_Buffer(ostCS.str().c_s 173 std::ifstream ES_Data_Buffer(ostCS.str().c_str(),ios::binary); 174 if( !ES_Data_Buffer.is_open() ) 174 if( !ES_Data_Buffer.is_open() ) 175 { 175 { 176 G4ExceptionDescription ed; 176 G4ExceptionDescription ed; 177 ed << "G4JAEAPolarizedElasticScattering 177 ed << "G4JAEAPolarizedElasticScattering Model data file <" << ostCS.str().c_str() 178 << "> is not opened!" << G4endl; 178 << "> is not opened!" << G4endl; 179 G4Exception("G4JAEAPolarizedElasticScatt 179 G4Exception("G4JAEAPolarizedElasticScatteringModel::ReadData()","em0003",FatalException, 180 ed,"G4LEDATA version should be G4EMLOW7. 180 ed,"G4LEDATA version should be G4EMLOW7.11 or later. Polarized Elastic Scattering Data are not loaded"); 181 return; 181 return; 182 } 182 } 183 else 183 else 184 { 184 { 185 if(verboseLevel > 3) { 185 if(verboseLevel > 3) { 186 G4cout << "File " << ostCS.str() 186 G4cout << "File " << ostCS.str() 187 << " is opened by G4JAEAPolarizedElas 187 << " is opened by G4JAEAPolarizedElasticScatteringModel" << G4endl; 188 } 188 } 189 } 189 } 190 190 191 191 192 if (!Polarized_ES_Data[Z]) 192 if (!Polarized_ES_Data[Z]) 193 Polarized_ES_Data[Z] = new G4DataVector(); 193 Polarized_ES_Data[Z] = new G4DataVector(); 194 194 195 G4float buffer_var; 195 G4float buffer_var; 196 while (ES_Data_Buffer.read(reinterpret_cast< 196 while (ES_Data_Buffer.read(reinterpret_cast<char*>(&buffer_var),sizeof(float))) 197 { 197 { 198 Polarized_ES_Data[Z]->push_back(buffer_v 198 Polarized_ES_Data[Z]->push_back(buffer_var); 199 } 199 } 200 200 201 dataCS[Z] = new G4PhysicsFreeVector(300,0.01 201 dataCS[Z] = new G4PhysicsFreeVector(300,0.01,3.,/*spline=*/true); 202 202 203 for (G4int i=0;i<300;++i) 203 for (G4int i=0;i<300;++i) 204 dataCS[Z]->PutValues(i,10.*i*1e-3,Polarize 204 dataCS[Z]->PutValues(i,10.*i*1e-3,Polarized_ES_Data[Z]->at(i)*1e-22); 205 205 206 // Activation of spline interpolation 206 // Activation of spline interpolation 207 dataCS[Z] ->FillSecondDerivatives(); 207 dataCS[Z] ->FillSecondDerivatives(); 208 } 208 } 209 209 210 //....oooOO0OOooo........oooOO0OOooo........oo 210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 211 211 212 G4double G4JAEAPolarizedElasticScatteringModel 212 G4double G4JAEAPolarizedElasticScatteringModel::ComputeCrossSectionPerAtom( 213 const G4ParticleDefinitio 213 const G4ParticleDefinition*, 214 G4double GammaEnergy, 214 G4double GammaEnergy, 215 G4double Z, G4double, 215 G4double Z, G4double, 216 G4double, G4double) 216 G4double, G4double) 217 { 217 { 218 //Select the energy-grid point closest to th 218 //Select the energy-grid point closest to the photon energy 219 // G4double *whichenergy = lower_bound(ES 219 // G4double *whichenergy = lower_bound(ESdata[0],ESdata[0]+300,GammaEnergy); 220 // int energyindex = max(0,(int)(whichene 220 // int energyindex = max(0,(int)(whichenergy-ESdata[0]-1)); 221 221 222 if (verboseLevel > 1) 222 if (verboseLevel > 1) 223 { 223 { 224 G4cout << "G4JAEAPolarizedElasticScatter 224 G4cout << "G4JAEAPolarizedElasticScatteringModel::ComputeCrossSectionPerAtom()" 225 << G4endl; 225 << G4endl; 226 } 226 } 227 227 228 if(GammaEnergy < lowEnergyLimit) { return 0. 228 if(GammaEnergy < lowEnergyLimit) { return 0.0; } 229 229 230 G4double xs = 0.0; 230 G4double xs = 0.0; 231 231 232 G4int intZ = G4lrint(Z); 232 G4int intZ = G4lrint(Z); 233 233 234 if(intZ < 1 || intZ > maxZ) { return xs; } 234 if(intZ < 1 || intZ > maxZ) { return xs; } 235 235 236 G4PhysicsFreeVector* pv = dataCS[intZ]; 236 G4PhysicsFreeVector* pv = dataCS[intZ]; 237 237 238 // if element was not initialised 238 // if element was not initialised 239 // do initialisation safely for MT mode 239 // do initialisation safely for MT mode 240 if(!pv) { 240 if(!pv) { 241 InitialiseForElement(0, intZ); 241 InitialiseForElement(0, intZ); 242 pv = dataCS[intZ]; 242 pv = dataCS[intZ]; 243 if(!pv) { return xs; } 243 if(!pv) { return xs; } 244 } 244 } 245 245 246 std::size_t n = pv->GetVectorLength() - 1; 246 std::size_t n = pv->GetVectorLength() - 1; 247 247 248 G4double e = GammaEnergy; 248 G4double e = GammaEnergy; 249 if(e >= pv->Energy(n)) { 249 if(e >= pv->Energy(n)) { 250 xs = (*pv)[n]; 250 xs = (*pv)[n]; 251 } else if(e >= pv->Energy(0)) { 251 } else if(e >= pv->Energy(0)) { 252 xs = pv->Value(e); 252 xs = pv->Value(e); 253 } 253 } 254 254 255 if(verboseLevel > 0) 255 if(verboseLevel > 0) 256 { 256 { 257 G4cout << "****** DEBUG: tcs value for 257 G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)=" 258 << e << G4endl; 258 << e << G4endl; 259 G4cout << " cs (Geant4 internal unit) 259 G4cout << " cs (Geant4 internal unit)=" << xs << G4endl; 260 G4cout << " -> first E*E*cs value i 260 G4cout << " -> first E*E*cs value in CS data file (iu) =" << (*pv)[0] 261 << G4endl; 261 << G4endl; 262 G4cout << " -> last E*E*cs value i 262 G4cout << " -> last E*E*cs value in CS data file (iu) =" << (*pv)[n] 263 << G4endl; 263 << G4endl; 264 G4cout << "*************************** 264 G4cout << "*********************************************************" 265 << G4endl; 265 << G4endl; 266 } 266 } 267 267 268 return (xs); 268 return (xs); 269 } 269 } 270 270 271 //....oooOO0OOooo........oooOO0OOooo........oo 271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 272 272 273 void G4JAEAPolarizedElasticScatteringModel::Sa 273 void G4JAEAPolarizedElasticScatteringModel::SampleSecondaries( 274 std::vector<G4DynamicParti 274 std::vector<G4DynamicParticle*>*, 275 const G4MaterialCutsCouple 275 const G4MaterialCutsCouple* couple, 276 const G4DynamicParticle* a 276 const G4DynamicParticle* aDynamicGamma, 277 G4double, G4double) 277 G4double, G4double) 278 { 278 { 279 if (verboseLevel > 1) { 279 if (verboseLevel > 1) { 280 280 281 G4cout << "Calling SampleSecondaries() of 281 G4cout << "Calling SampleSecondaries() of G4JAEAPolarizedElasticScatteringModel." 282 << G4endl; 282 << G4endl; 283 } 283 } 284 G4double photonEnergy0 = aDynamicGamma->GetK 284 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy(); 285 285 286 // absorption of low-energy gamma 286 // absorption of low-energy gamma 287 if (photonEnergy0 <= lowEnergyLimit) 287 if (photonEnergy0 <= lowEnergyLimit) 288 { 288 { 289 fParticleChange->ProposeTrackStatus(fSto 289 fParticleChange->ProposeTrackStatus(fStopAndKill); 290 fParticleChange->SetProposedKineticEnerg 290 fParticleChange->SetProposedKineticEnergy(0.); 291 fParticleChange->ProposeLocalEnergyDepos 291 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0); 292 return ; 292 return ; 293 } 293 } 294 294 295 const G4ParticleDefinition* particle = aDyn 295 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition(); 296 const G4Element* elm = SelectRandomAtom(coup 296 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0); 297 G4int Z = G4lrint(elm->GetZ()); 297 G4int Z = G4lrint(elm->GetZ()); 298 298 299 //Getting the corresponding distrbution 299 //Getting the corresponding distrbution 300 G4int energyindex=round(100*photonEnergy0)-1 300 G4int energyindex=round(100*photonEnergy0)-1; 301 G4double a1=0, a2=0, a3=0,a4=0; 301 G4double a1=0, a2=0, a3=0,a4=0; 302 for (G4int i=0;i<=180;++i) 302 for (G4int i=0;i<=180;++i) 303 { 303 { 304 a1=Polarized_ES_Data[Z]->at(4*i+300+181* 304 a1=Polarized_ES_Data[Z]->at(4*i+300+181*4*(energyindex)); 305 a2=Polarized_ES_Data[Z]->at(4*i+1+300+18 305 a2=Polarized_ES_Data[Z]->at(4*i+1+300+181*4*(energyindex)); 306 a3=Polarized_ES_Data[Z]->at(4*i+2+300+18 306 a3=Polarized_ES_Data[Z]->at(4*i+2+300+181*4*(energyindex)); 307 a4=Polarized_ES_Data[Z]->at(4*i+3+300+18 307 a4=Polarized_ES_Data[Z]->at(4*i+3+300+181*4*(energyindex)); 308 distribution[i]=a1*a1+a2*a2+a3*a3+a4*a4; 308 distribution[i]=a1*a1+a2*a2+a3*a3+a4*a4; 309 } 309 } 310 310 311 CLHEP::RandGeneral GenThetaDist(distribution 311 CLHEP::RandGeneral GenThetaDist(distribution,180); 312 //Intial sampling of the scattering angle. T 312 //Intial sampling of the scattering angle. To be updated for the circular polarization 313 G4double theta = CLHEP::pi*GenThetaDist.shoo 313 G4double theta = CLHEP::pi*GenThetaDist.shoot(); 314 //G4double theta =45.*CLHEP::pi/180.; 314 //G4double theta =45.*CLHEP::pi/180.; 315 //Theta is in degree to call scattering ampl 315 //Theta is in degree to call scattering amplitudes 316 G4int theta_in_degree =round(theta*180./CLHE 316 G4int theta_in_degree =round(theta*180./CLHEP::pi); 317 317 318 //theta_in_degree=45; 318 //theta_in_degree=45; 319 319 320 G4double am1=0,am2=0,am3=0,am4=0,aparaSquare 320 G4double am1=0,am2=0,am3=0,am4=0,aparaSquare=0,aperpSquare=0, 321 apara_aper_Asterisk=0,img_apara_aper_Aster 321 apara_aper_Asterisk=0,img_apara_aper_Asterisk=0; 322 am1=Polarized_ES_Data[Z]->at(4*theta_in_degr 322 am1=Polarized_ES_Data[Z]->at(4*theta_in_degree+300+181*4*(energyindex)); 323 am2=Polarized_ES_Data[Z]->at(4*theta_in_degr 323 am2=Polarized_ES_Data[Z]->at(4*theta_in_degree+1+300+181*4*(energyindex)); 324 am3=Polarized_ES_Data[Z]->at(4*theta_in_degr 324 am3=Polarized_ES_Data[Z]->at(4*theta_in_degree+2+300+181*4*(energyindex)); 325 am4=Polarized_ES_Data[Z]->at(4*theta_in_degr 325 am4=Polarized_ES_Data[Z]->at(4*theta_in_degree+3+300+181*4*(energyindex)); 326 aparaSquare=am1*am1+am2*am2; 326 aparaSquare=am1*am1+am2*am2; 327 aperpSquare=am3*am3+am4*am4; 327 aperpSquare=am3*am3+am4*am4; 328 apara_aper_Asterisk=2*a1*a3+2*a2*a4; 328 apara_aper_Asterisk=2*a1*a3+2*a2*a4; 329 img_apara_aper_Asterisk=2*a1*a4-2*a2*a3; 329 img_apara_aper_Asterisk=2*a1*a4-2*a2*a3; 330 330 331 G4ThreeVector Direction_Unpolarized(0.,0.,0. 331 G4ThreeVector Direction_Unpolarized(0.,0.,0.); 332 G4ThreeVector Direction_Linear1(0.,0.,0.); 332 G4ThreeVector Direction_Linear1(0.,0.,0.); 333 G4ThreeVector Direction_Linear2(0.,0.,0.); 333 G4ThreeVector Direction_Linear2(0.,0.,0.); 334 G4ThreeVector Direction_Circular(0.,0.,0.); 334 G4ThreeVector Direction_Circular(0.,0.,0.); 335 G4ThreeVector Polarization_Unpolarized(0.,0. 335 G4ThreeVector Polarization_Unpolarized(0.,0.,0.); 336 G4ThreeVector Polarization_Linear1(0.,0.,0.) 336 G4ThreeVector Polarization_Linear1(0.,0.,0.); 337 G4ThreeVector Polarization_Linear2(0.,0.,0.) 337 G4ThreeVector Polarization_Linear2(0.,0.,0.); 338 G4ThreeVector Polarization_Circular(0.,0.,0. 338 G4ThreeVector Polarization_Circular(0.,0.,0.); 339 339 340 //Stokes parameters for the incoming and out 340 //Stokes parameters for the incoming and outgoing photon 341 G4double Xi1=0, Xi2=0, Xi3=0, Xi1_Prime=0,Xi 341 G4double Xi1=0, Xi2=0, Xi3=0, Xi1_Prime=0,Xi2_Prime=0,Xi3_Prime=0; 342 342 343 //Getting the Stokes parameters for the inco 343 //Getting the Stokes parameters for the incoming photon 344 G4ThreeVector gammaPolarization0 = aDynamicG 344 G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization(); 345 Xi1=gammaPolarization0.x(); 345 Xi1=gammaPolarization0.x(); 346 Xi2=gammaPolarization0.y(); 346 Xi2=gammaPolarization0.y(); 347 Xi3=gammaPolarization0.z(); 347 Xi3=gammaPolarization0.z(); 348 348 349 //Polarization vector must be unit vector (5 349 //Polarization vector must be unit vector (5% tolerance) 350 if ((gammaPolarization0.mag())>1.05 || (Xi1* 350 if ((gammaPolarization0.mag())>1.05 || (Xi1*Xi1>1.05) || (Xi2*Xi2>1.05) || (Xi3*Xi3>1.05)) 351 { 351 { 352 G4Exception("G4JAEAPolarizedElasticScatt 352 G4Exception("G4JAEAPolarizedElasticScatteringModel::SampleSecondaries()","em1006", 353 JustWarning, 353 JustWarning, 354 "WARNING: G4JAEAPolarizedElasticScatteri 354 "WARNING: G4JAEAPolarizedElasticScatteringModel is only compatible with a unit polarization vector."); 355 return; 355 return; 356 } 356 } 357 //Unpolarized gamma rays 357 //Unpolarized gamma rays 358 if (Xi1==0 && Xi2==0 && Xi3==0) 358 if (Xi1==0 && Xi2==0 && Xi3==0) 359 { 359 { 360 G4double Phi_Unpolarized=0; 360 G4double Phi_Unpolarized=0; 361 if (fLinearPolarizationSensitvity1) 361 if (fLinearPolarizationSensitvity1) 362 Phi_Unpolarized=GeneratePolarizedPhi(aparaSq 362 Phi_Unpolarized=GeneratePolarizedPhi(aparaSquare,aperpSquare,0.); 363 else 363 else 364 Phi_Unpolarized=CLHEP::twopi*G4UniformRand() 364 Phi_Unpolarized=CLHEP::twopi*G4UniformRand(); 365 Direction_Unpolarized.setX(sin(theta)*co 365 Direction_Unpolarized.setX(sin(theta)*cos(Phi_Unpolarized)); 366 Direction_Unpolarized.setY(sin(theta)*si 366 Direction_Unpolarized.setY(sin(theta)*sin(Phi_Unpolarized)); 367 Direction_Unpolarized.setZ(cos(theta)); 367 Direction_Unpolarized.setZ(cos(theta)); 368 Direction_Unpolarized.rotateUz(aDynamicG 368 Direction_Unpolarized.rotateUz(aDynamicGamma->GetMomentumDirection()); 369 Xi1_Prime=(aparaSquare-aperpSquare)/(apa 369 Xi1_Prime=(aparaSquare-aperpSquare)/(aparaSquare+aperpSquare); 370 Polarization_Unpolarized.setX(Xi1_Prime) 370 Polarization_Unpolarized.setX(Xi1_Prime); 371 Polarization_Unpolarized.setY(0.); 371 Polarization_Unpolarized.setY(0.); 372 Polarization_Unpolarized.setZ(0.); 372 Polarization_Unpolarized.setZ(0.); 373 fParticleChange->ProposeMomentumDirectio 373 fParticleChange->ProposeMomentumDirection(Direction_Unpolarized); 374 fParticleChange->ProposePolarization(Pol 374 fParticleChange->ProposePolarization(Polarization_Unpolarized); 375 return; 375 return; 376 } 376 } 377 377 378 //Linear polarization defined by first Stoke 378 //Linear polarization defined by first Stokes parameter 379 G4double InitialAzimuth=aDynamicGamma->GetMo 379 G4double InitialAzimuth=aDynamicGamma->GetMomentumDirection().phi(); 380 if(InitialAzimuth<0) InitialAzimuth=InitialA 380 if(InitialAzimuth<0) InitialAzimuth=InitialAzimuth+CLHEP::twopi; 381 381 382 G4double Phi_Linear1=0.; 382 G4double Phi_Linear1=0.; 383 383 384 Phi_Linear1 = GeneratePolarizedPhi(aparaSqua 384 Phi_Linear1 = GeneratePolarizedPhi(aparaSquare+aperpSquare+Xi1*(aparaSquare-aperpSquare), 385 aparaSquare+aperpSquare-Xi1*(apar 385 aparaSquare+aperpSquare-Xi1*(aparaSquare-aperpSquare),InitialAzimuth); 386 386 387 Xi1_Prime=((aparaSquare-aperpSquare)+Xi1*(ap 387 Xi1_Prime=((aparaSquare-aperpSquare)+Xi1*(aparaSquare+aperpSquare)*cos(2*Phi_Linear1))/ 388 ((aparaSquare+aperpSquare)+Xi1*(aparaSquar 388 ((aparaSquare+aperpSquare)+Xi1*(aparaSquare-aperpSquare)*cos(2*Phi_Linear1)); 389 Xi2_Prime=(-Xi1*apara_aper_Asterisk*sin(2*Ph 389 Xi2_Prime=(-Xi1*apara_aper_Asterisk*sin(2*Phi_Linear1))/ 390 ((aparaSquare+aperpSquare)+Xi1*(aparaSquar 390 ((aparaSquare+aperpSquare)+Xi1*(aparaSquare-aperpSquare)*cos(2*Phi_Linear1)); 391 Xi3_Prime=(-Xi1*img_apara_aper_Asterisk*sin( 391 Xi3_Prime=(-Xi1*img_apara_aper_Asterisk*sin(2*Phi_Linear1))/ 392 ((aparaSquare+aperpSquare)+Xi1*(aparaSquar 392 ((aparaSquare+aperpSquare)+Xi1*(aparaSquare-aperpSquare)*cos(2*Phi_Linear1)); 393 //Store momentum direction and po;arization 393 //Store momentum direction and po;arization 394 Direction_Linear1.setX(sin(theta)*cos(Phi_Li 394 Direction_Linear1.setX(sin(theta)*cos(Phi_Linear1)); 395 Direction_Linear1.setY(sin(theta)*sin(Phi_Li 395 Direction_Linear1.setY(sin(theta)*sin(Phi_Linear1)); 396 Direction_Linear1.setZ(cos(theta)); 396 Direction_Linear1.setZ(cos(theta)); 397 Polarization_Linear1.setX(Xi1_Prime); 397 Polarization_Linear1.setX(Xi1_Prime); 398 Polarization_Linear1.setY(Xi2_Prime); 398 Polarization_Linear1.setY(Xi2_Prime); 399 Polarization_Linear1.setZ(Xi3_Prime); 399 Polarization_Linear1.setZ(Xi3_Prime); 400 400 401 //Set scattered photon polarization sensitiv 401 //Set scattered photon polarization sensitivity 402 Xi1_Prime=Xi1_Prime*fLinearPolarizationSensi 402 Xi1_Prime=Xi1_Prime*fLinearPolarizationSensitvity1; 403 Xi2_Prime=Xi2_Prime*fLinearPolarizationSensi 403 Xi2_Prime=Xi2_Prime*fLinearPolarizationSensitvity2; 404 Xi3_Prime=Xi3_Prime*fCircularPolarizationSen 404 Xi3_Prime=Xi3_Prime*fCircularPolarizationSensitvity; 405 405 406 G4double dsigmaL1=0.0; 406 G4double dsigmaL1=0.0; 407 if(abs(Xi1)>0.0) dsigmaL1=0.25*((aparaSquare 407 if(abs(Xi1)>0.0) dsigmaL1=0.25*((aparaSquare+aperpSquare)* 408 (1+Xi1*Xi1_Prime*cos(2*Phi_Linear1)) 408 (1+Xi1*Xi1_Prime*cos(2*Phi_Linear1))+ 409 (aparaSquare-aperpSquare)*(Xi1*cos(2 409 (aparaSquare-aperpSquare)*(Xi1*cos(2*Phi_Linear1)+Xi1_Prime) 410 -Xi1*Xi2_Prime*apara_aper_Asterisk*s 410 -Xi1*Xi2_Prime*apara_aper_Asterisk*sin(2*Phi_Linear1)- 411 Xi1*Xi3_Prime*img_apara_aper_Asteris 411 Xi1*Xi3_Prime*img_apara_aper_Asterisk*sin(2*Phi_Linear1)); 412 412 413 //Linear polarization defined by second Stok 413 //Linear polarization defined by second Stokes parameter 414 //G4double IntialAzimuth=aDynamicGamma->GetM 414 //G4double IntialAzimuth=aDynamicGamma->GetMomentumDirection().phi(); 415 G4double Phi_Linear2=0.; 415 G4double Phi_Linear2=0.; 416 416 417 InitialAzimuth=InitialAzimuth-CLHEP::pi/4.; 417 InitialAzimuth=InitialAzimuth-CLHEP::pi/4.; 418 if(InitialAzimuth<0) InitialAzimuth=InitialA 418 if(InitialAzimuth<0) InitialAzimuth=InitialAzimuth+CLHEP::twopi; 419 419 420 Phi_Linear2 = GeneratePolarizedPhi(aparaSqua 420 Phi_Linear2 = GeneratePolarizedPhi(aparaSquare+aperpSquare+Xi1*(aparaSquare-aperpSquare) 421 ,aparaSquare+aperpSquare-Xi1*(apa 421 ,aparaSquare+aperpSquare-Xi1*(aparaSquare-aperpSquare),InitialAzimuth); 422 422 423 Xi1_Prime=((aparaSquare-aperpSquare)+Xi2*(ap 423 Xi1_Prime=((aparaSquare-aperpSquare)+Xi2*(aparaSquare+aperpSquare)*sin(2*Phi_Linear2))/ 424 ((aparaSquare+aperpSquare)+Xi2*(aparaSquar 424 ((aparaSquare+aperpSquare)+Xi2*(aparaSquare-aperpSquare)*sin(2*Phi_Linear2)); 425 Xi2_Prime=(Xi2*apara_aper_Asterisk*cos(2*Phi 425 Xi2_Prime=(Xi2*apara_aper_Asterisk*cos(2*Phi_Linear2))/ 426 ((aparaSquare+aperpSquare)+Xi2*(aparaSquar 426 ((aparaSquare+aperpSquare)+Xi2*(aparaSquare-aperpSquare)*sin(2*Phi_Linear2)); 427 Xi3_Prime=(Xi2*img_apara_aper_Asterisk*cos(2 427 Xi3_Prime=(Xi2*img_apara_aper_Asterisk*cos(2*Phi_Linear2))/ 428 ((aparaSquare+aperpSquare)+Xi2*(aparaSquar 428 ((aparaSquare+aperpSquare)+Xi2*(aparaSquare-aperpSquare)*sin(2*Phi_Linear2)); 429 //Store momentum direction and polarization 429 //Store momentum direction and polarization 430 Direction_Linear2.setX(sin(theta)*cos(Phi_Li 430 Direction_Linear2.setX(sin(theta)*cos(Phi_Linear2)); 431 Direction_Linear2.setY(sin(theta)*sin(Phi_Li 431 Direction_Linear2.setY(sin(theta)*sin(Phi_Linear2)); 432 Direction_Linear2.setZ(cos(theta)); 432 Direction_Linear2.setZ(cos(theta)); 433 Polarization_Linear2.setX(Xi1_Prime); 433 Polarization_Linear2.setX(Xi1_Prime); 434 Polarization_Linear2.setY(Xi2_Prime); 434 Polarization_Linear2.setY(Xi2_Prime); 435 Polarization_Linear2.setZ(Xi3_Prime); 435 Polarization_Linear2.setZ(Xi3_Prime); 436 436 437 //Set scattered photon polarization sensitiv 437 //Set scattered photon polarization sensitivity 438 Xi1_Prime=Xi1_Prime*fLinearPolarizationSensi 438 Xi1_Prime=Xi1_Prime*fLinearPolarizationSensitvity1; 439 Xi2_Prime=Xi2_Prime*fLinearPolarizationSensi 439 Xi2_Prime=Xi2_Prime*fLinearPolarizationSensitvity2; 440 Xi3_Prime=Xi3_Prime*fCircularPolarizationSen 440 Xi3_Prime=Xi3_Prime*fCircularPolarizationSensitvity; 441 441 442 G4double dsigmaL2=0.0; 442 G4double dsigmaL2=0.0; 443 if(abs(Xi2)>0.0) 443 if(abs(Xi2)>0.0) 444 dsigmaL2=0.25*((aparaSquare+aperpSquare)*( 444 dsigmaL2=0.25*((aparaSquare+aperpSquare)*(1+Xi2*Xi1_Prime*sin(2*Phi_Linear2))+ 445 (aparaSquare-aperpSquare)*(Xi2*sin(2*Ph 445 (aparaSquare-aperpSquare)*(Xi2*sin(2*Phi_Linear2)+Xi1_Prime) 446 +Xi2*Xi2_Prime*apara_aper_Asterisk*cos( 446 +Xi2*Xi2_Prime*apara_aper_Asterisk*cos(2*Phi_Linear2)- 447 Xi2*Xi3_Prime*img_apara_aper_Asterisk*c 447 Xi2*Xi3_Prime*img_apara_aper_Asterisk*cos(2*Phi_Linear2)); 448 448 449 //Circular polarization 449 //Circular polarization 450 G4double Phi_Circular = CLHEP::twopi*G4Unifo 450 G4double Phi_Circular = CLHEP::twopi*G4UniformRand(); 451 G4double Theta_Circular = 0.; 451 G4double Theta_Circular = 0.; 452 452 453 Xi1_Prime=(aparaSquare-aperpSquare)/(aparaSq 453 Xi1_Prime=(aparaSquare-aperpSquare)/(aparaSquare+aperpSquare); 454 Xi2_Prime=(-Xi3*img_apara_aper_Asterisk)/(ap 454 Xi2_Prime=(-Xi3*img_apara_aper_Asterisk)/(aparaSquare+aperpSquare); 455 Xi3_Prime=(Xi3*apara_aper_Asterisk)/(aparaSq 455 Xi3_Prime=(Xi3*apara_aper_Asterisk)/(aparaSquare+aperpSquare); 456 456 457 Polarization_Circular.setX(Xi1_Prime); 457 Polarization_Circular.setX(Xi1_Prime); 458 Polarization_Circular.setY(Xi2_Prime); 458 Polarization_Circular.setY(Xi2_Prime); 459 Polarization_Circular.setZ(Xi3_Prime); 459 Polarization_Circular.setZ(Xi3_Prime); 460 460 461 //Set scattered photon polarization sensitiv 461 //Set scattered photon polarization sensitivity 462 Xi1_Prime=Xi1_Prime*fLinearPolarizationSensi 462 Xi1_Prime=Xi1_Prime*fLinearPolarizationSensitvity1; 463 Xi2_Prime=Xi2_Prime*fLinearPolarizationSensi 463 Xi2_Prime=Xi2_Prime*fLinearPolarizationSensitvity2; 464 Xi3_Prime=Xi3_Prime*fCircularPolarizationSen 464 Xi3_Prime=Xi3_Prime*fCircularPolarizationSensitvity; 465 465 466 G4double dsigmaC=0.0; 466 G4double dsigmaC=0.0; 467 if(abs(Xi3)>0.0) 467 if(abs(Xi3)>0.0) 468 dsigmaC=0.25*(aparaSquare+aperpSquare+Xi1_ 468 dsigmaC=0.25*(aparaSquare+aperpSquare+Xi1_Prime*(aparaSquare-aperpSquare)- 469 Xi3*Xi2_Prime*img_apara_aper_Asterisk 469 Xi3*Xi2_Prime*img_apara_aper_Asterisk 470 +Xi3*Xi3_Prime*apara_aper_Asterisk); 470 +Xi3*Xi3_Prime*apara_aper_Asterisk); 471 471 472 if (abs(Xi3)==0.0 && abs(Xi1_Prime)==0.0) 472 if (abs(Xi3)==0.0 && abs(Xi1_Prime)==0.0) 473 { 473 { 474 Direction_Circular.setX(sin(theta)*cos(P 474 Direction_Circular.setX(sin(theta)*cos(Phi_Circular)); 475 Direction_Circular.setY(sin(theta)*sin(P 475 Direction_Circular.setY(sin(theta)*sin(Phi_Circular)); 476 Direction_Circular.setZ(cos(theta)); 476 Direction_Circular.setZ(cos(theta)); 477 } 477 } 478 else 478 else 479 { 479 { 480 G4double c1=0, c2=0, c3=0,c4=0; 480 G4double c1=0, c2=0, c3=0,c4=0; 481 for (G4int i=0;i<=180;++i) 481 for (G4int i=0;i<=180;++i) 482 { 482 { 483 c1=Polarized_ES_Data[Z]->at(4*i+300+181*4* 483 c1=Polarized_ES_Data[Z]->at(4*i+300+181*4*(energyindex)); 484 c2=Polarized_ES_Data[Z]->at(4*i+1+300+181* 484 c2=Polarized_ES_Data[Z]->at(4*i+1+300+181*4*(energyindex)); 485 c3=Polarized_ES_Data[Z]->at(4*i+2+300+181* 485 c3=Polarized_ES_Data[Z]->at(4*i+2+300+181*4*(energyindex)); 486 c4=Polarized_ES_Data[Z]->at(4*i+3+300+181* 486 c4=Polarized_ES_Data[Z]->at(4*i+3+300+181*4*(energyindex)); 487 cdistribution[i]=0.25*((c1*c1+c2*c2+c3*c3+ 487 cdistribution[i]=0.25*((c1*c1+c2*c2+c3*c3+c4*c4)+ 488 Xi1_Prime*(c1*c1+c2*c2-c3*c3-c4*c4)- 488 Xi1_Prime*(c1*c1+c2*c2-c3*c3-c4*c4)- 489 Xi3*Xi2_Prime*(2*c1*c4-2*c2*c3) 489 Xi3*Xi2_Prime*(2*c1*c4-2*c2*c3) 490 +Xi3*Xi3_Prime*(2*c1*c4-2*c2*c3)); 490 +Xi3*Xi3_Prime*(2*c1*c4-2*c2*c3)); 491 } 491 } 492 CLHEP::RandGeneral GenTheta_Circ_Dist(cd 492 CLHEP::RandGeneral GenTheta_Circ_Dist(cdistribution,180); 493 Theta_Circular=CLHEP::pi*GenTheta_Circ_D 493 Theta_Circular=CLHEP::pi*GenTheta_Circ_Dist.shoot(); 494 Direction_Circular.setX(sin(Theta_Circul 494 Direction_Circular.setX(sin(Theta_Circular)*cos(Phi_Circular)); 495 Direction_Circular.setY(sin(Theta_Circul 495 Direction_Circular.setY(sin(Theta_Circular)*sin(Phi_Circular)); 496 Direction_Circular.setZ(cos(Theta_Circul 496 Direction_Circular.setZ(cos(Theta_Circular)); 497 } 497 } 498 498 499 // Sampling scattered photon direction based 499 // Sampling scattered photon direction based on asymmetry arising from polarization mixing 500 G4double totalSigma= dsigmaL1+dsigmaL2+dsigm 500 G4double totalSigma= dsigmaL1+dsigmaL2+dsigmaC; 501 G4double prob1=dsigmaL1/totalSigma; 501 G4double prob1=dsigmaL1/totalSigma; 502 G4double prob2=dsigmaL2/totalSigma; 502 G4double prob2=dsigmaL2/totalSigma; 503 G4double probc=1-(prob1+prob2); 503 G4double probc=1-(prob1+prob2); 504 504 505 //Check the Probability of polarization mixi 505 //Check the Probability of polarization mixing 506 if (abs(probc - dsigmaC/totalSigma)>=0.0001) 506 if (abs(probc - dsigmaC/totalSigma)>=0.0001) 507 { 507 { 508 G4Exception("G4JAEAPolarizedElasticScatt 508 G4Exception("G4JAEAPolarizedElasticScatteringModel::SampleSecondaries()","em1007", 509 JustWarning, 509 JustWarning, 510 "WARNING: Polarization mixing might be i 510 "WARNING: Polarization mixing might be incorrect."); 511 } 511 } 512 512 513 // Generate outgoing photon direction 513 // Generate outgoing photon direction 514 G4ThreeVector finaldirection(0.0,0.0,0.0); 514 G4ThreeVector finaldirection(0.0,0.0,0.0); 515 G4ThreeVector outcomingPhotonPolarization(0. 515 G4ThreeVector outcomingPhotonPolarization(0.0,0.0,0.0); 516 516 517 //Polarization mixing 517 //Polarization mixing 518 G4double polmix=G4UniformRand(); 518 G4double polmix=G4UniformRand(); 519 if (polmix<=prob1) 519 if (polmix<=prob1) 520 { 520 { 521 finaldirection.setX(Direction_Linear1.x( 521 finaldirection.setX(Direction_Linear1.x()); 522 finaldirection.setY(Direction_Linear1.y( 522 finaldirection.setY(Direction_Linear1.y()); 523 finaldirection.setZ(Direction_Linear1.z( 523 finaldirection.setZ(Direction_Linear1.z()); 524 outcomingPhotonPolarization.setX(Polariz 524 outcomingPhotonPolarization.setX(Polarization_Linear1.x()); 525 outcomingPhotonPolarization.setY(Polariz 525 outcomingPhotonPolarization.setY(Polarization_Linear1.y()); 526 outcomingPhotonPolarization.setZ(Polariz 526 outcomingPhotonPolarization.setZ(Polarization_Linear1.z()); 527 } 527 } 528 else if ((polmix>prob1) && (polmix<=prob1+pr 528 else if ((polmix>prob1) && (polmix<=prob1+prob2)) 529 { 529 { 530 finaldirection.setX(Direction_Linear2.x( 530 finaldirection.setX(Direction_Linear2.x()); 531 finaldirection.setY(Direction_Linear2.y( 531 finaldirection.setY(Direction_Linear2.y()); 532 finaldirection.setZ(Direction_Linear2.z( 532 finaldirection.setZ(Direction_Linear2.z()); 533 outcomingPhotonPolarization.setX(Polariz 533 outcomingPhotonPolarization.setX(Polarization_Linear2.x()); 534 outcomingPhotonPolarization.setY(Polariz 534 outcomingPhotonPolarization.setY(Polarization_Linear2.y()); 535 outcomingPhotonPolarization.setZ(Polariz 535 outcomingPhotonPolarization.setZ(Polarization_Linear2.z()); 536 } 536 } 537 else if (polmix>prob1+prob2) 537 else if (polmix>prob1+prob2) 538 { 538 { 539 finaldirection.setX(Direction_Circular.x 539 finaldirection.setX(Direction_Circular.x()); 540 finaldirection.setY(Direction_Circular.y 540 finaldirection.setY(Direction_Circular.y()); 541 finaldirection.setZ(Direction_Circular.z 541 finaldirection.setZ(Direction_Circular.z()); 542 outcomingPhotonPolarization.setX(Polariz 542 outcomingPhotonPolarization.setX(Polarization_Circular.x()); 543 outcomingPhotonPolarization.setY(Polariz 543 outcomingPhotonPolarization.setY(Polarization_Circular.y()); 544 outcomingPhotonPolarization.setZ(Polariz 544 outcomingPhotonPolarization.setZ(Polarization_Circular.z()); 545 } 545 } 546 546 547 //Sampling the Final State 547 //Sampling the Final State 548 finaldirection.rotateUz(aDynamicGamma->GetMo 548 finaldirection.rotateUz(aDynamicGamma->GetMomentumDirection()); 549 fParticleChange->ProposeMomentumDirection(fi 549 fParticleChange->ProposeMomentumDirection(finaldirection); 550 fParticleChange->SetProposedKineticEnergy(ph 550 fParticleChange->SetProposedKineticEnergy(photonEnergy0); 551 fParticleChange->ProposePolarization(outcomi 551 fParticleChange->ProposePolarization(outcomingPhotonPolarization); 552 552 553 } 553 } 554 554 555 //....oooOO0OOooo........oooOO0OOooo........oo 555 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 556 556 557 G4double G4JAEAPolarizedElasticScatteringModel 557 G4double G4JAEAPolarizedElasticScatteringModel::GeneratePolarizedPhi(G4double Sigma_para, 558 G4double Sigma_perp, 558 G4double Sigma_perp, 559 G4double initial_Pol_Plan 559 G4double initial_Pol_Plane) 560 { 560 { 561 G4double phi; 561 G4double phi; 562 G4double phiProbability; 562 G4double phiProbability; 563 G4double Probability=Sigma_perp/(Sigma_para+ 563 G4double Probability=Sigma_perp/(Sigma_para+Sigma_perp); 564 if (Probability<=G4UniformRand()) 564 if (Probability<=G4UniformRand()) 565 { 565 { 566 do 566 do 567 { 567 { 568 phi = CLHEP::twopi * G4UniformRand(); 568 phi = CLHEP::twopi * G4UniformRand(); 569 phiProbability = cos(phi+initial_Pol_Plane 569 phiProbability = cos(phi+initial_Pol_Plane)*cos(phi+initial_Pol_Plane); 570 } 570 } 571 while (phiProbability < G4UniformRand()) 571 while (phiProbability < G4UniformRand()); 572 572 573 } 573 } 574 else 574 else 575 { 575 { 576 do 576 do 577 { 577 { 578 phi = CLHEP::twopi * G4UniformRand(); 578 phi = CLHEP::twopi * G4UniformRand(); 579 phiProbability = sin(phi+initial_Pol_Plane 579 phiProbability = sin(phi+initial_Pol_Plane)*sin(phi+initial_Pol_Plane); 580 } 580 } 581 while (phiProbability < G4UniformRand()) 581 while (phiProbability < G4UniformRand()); 582 } 582 } 583 return phi; 583 return phi; 584 584 585 } 585 } 586 //....oooOO0OOooo........oooOO0OOooo........oo 586 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 587 587 588 void 588 void 589 G4JAEAPolarizedElasticScatteringModel::Initial 589 G4JAEAPolarizedElasticScatteringModel::InitialiseForElement(const G4ParticleDefinition*, 590 G4int Z) 590 G4int Z) 591 { 591 { 592 G4AutoLock l(&G4JAEAPolarizedElasticScatteri 592 G4AutoLock l(&G4JAEAPolarizedElasticScatteringModelMutex); 593 // G4cout << "G4JAEAPolarizedElasticScatter 593 // G4cout << "G4JAEAPolarizedElasticScatteringModel::InitialiseForElement Z= " 594 // << Z << G4endl; 594 // << Z << G4endl; 595 if(!dataCS[Z]) { ReadData(Z); } 595 if(!dataCS[Z]) { ReadData(Z); } 596 l.unlock(); 596 l.unlock(); 597 } 597 } 598 598 599 //....oooOO0OOooo........oooOO0OOooo........oo 599 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 600 600