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