Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4JAEAPolarizedElasticScatteringModel Class Reference

#include <G4JAEAPolarizedElasticScatteringModel.hh>

+ Inheritance diagram for G4JAEAPolarizedElasticScatteringModel:

Public Member Functions

 G4JAEAPolarizedElasticScatteringModel ()
 
virtual ~G4JAEAPolarizedElasticScatteringModel ()
 
void Initialise (const G4ParticleDefinition *, const G4DataVector &) override
 
void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel) override
 
void InitialiseForElement (const G4ParticleDefinition *, G4int Z) override
 
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX) override
 
void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
 
void SetLowEnergyThreshold (G4double val)
 
void SetPolarizationSensitvity (G4bool, G4bool, G4bool)
 
void SetDebugVerbosity (G4int val)
 
G4JAEAPolarizedElasticScatteringModel & operator= (const G4JAEAPolarizedElasticScatteringModel &right)=delete
 
 G4JAEAPolarizedElasticScatteringModel (const G4JAEAPolarizedElasticScatteringModel &)=delete
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
virtual G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double GetPartialCrossSection (const G4Material *, G4int level, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double ComputeCrossSectionPerShell (const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ChargeSquareRatio (const G4Track &)
 
virtual G4double GetChargeSquareRatio (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual G4double GetParticleCharge (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void StartTracking (G4Track *)
 
virtual void CorrectionsAlongStep (const G4MaterialCutsCouple *, const G4DynamicParticle *, const G4double &length, G4double &eloss)
 
virtual G4double Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
 
virtual G4double MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
 
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void DefineForRegion (const G4Region *)
 
virtual void FillNumberOfSecondaries (G4int &numberOfTriplets, G4int &numberOfRecoil)
 
virtual void ModelDescription (std::ostream &outFile) const
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector< G4EmElementSelector * > * GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector * > *)
 
G4double ComputeDEDX (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
G4double CrossSection (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeMeanFreePath (const G4ParticleDefinition *, G4double kineticEnergy, const G4Material *, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, const G4Element *, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectTargetAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double logKineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * GetCurrentElement (const G4Material *mat=nullptr) const
 
G4int SelectRandomAtomNumber (const G4Material *) const
 
const G4Isotope * GetCurrentIsotope (const G4Element *elm=nullptr) const
 
G4int SelectIsotopeNumber (const G4Element *) const
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementData * GetElementData ()
 
G4PhysicsTable * GetCrossSectionTable ()
 
G4VEmFluctuationModel * GetModelOfFluctuations ()
 
G4VEmAngularDistribution * GetAngularDistribution ()
 
G4VEmModel * GetTripletModel ()
 
void SetTripletModel (G4VEmModel *)
 
void SetAngularDistribution (G4VEmAngularDistribution *)
 
G4double HighEnergyLimit () const
 
G4double LowEnergyLimit () const
 
G4double HighEnergyActivationLimit () const
 
G4double LowEnergyActivationLimit () const
 
G4double PolarAngleLimit () const
 
G4double SecondaryThreshold () const
 
G4bool DeexcitationFlag () const
 
G4bool ForceBuildTableFlag () const
 
G4bool UseAngularGeneratorFlag () const
 
void SetAngularGeneratorFlag (G4bool)
 
void SetHighEnergyLimit (G4double)
 
void SetLowEnergyLimit (G4double)
 
void SetActivationHighEnergyLimit (G4double)
 
void SetActivationLowEnergyLimit (G4double)
 
G4bool IsActive (G4double kinEnergy) const
 
void SetPolarAngleLimit (G4double)
 
void SetSecondaryThreshold (G4double)
 
void SetDeexcitationFlag (G4bool val)
 
void SetForceBuildTable (G4bool val)
 
void SetFluctuationFlag (G4bool val)
 
void SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
void SetUseBaseMaterials (G4bool val)
 
G4bool UseBaseMaterials () const
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4String & GetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 
void SetLPMFlag (G4bool)
 
G4VEmModel & operator= (const G4VEmModel &right)=delete
 
 G4VEmModel (const G4VEmModel &)=delete
 

Additional Inherited Members

- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLoss * GetParticleChangeForLoss ()
 
G4ParticleChangeForGamma * GetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCouple * CurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 
- Protected Attributes inherited from G4VEmModel
G4ElementData * fElementData = nullptr
 
G4VParticleChange * pParticleChange = nullptr
 
G4PhysicsTable * xSectionTable = nullptr
 
const G4Material * pBaseMaterial = nullptr
 
const std::vector< G4double > * theDensityFactor = nullptr
 
const std::vector< G4int > * theDensityIdx = nullptr
 
G4double inveplus
 
G4double pFactor = 1.0
 
std::size_t currentCoupleIndex = 0
 
std::size_t basedCoupleIndex = 0
 
G4bool lossFlucFlag = true
 

Detailed Description

Definition at line 46 of file G4JAEAPolarizedElasticScatteringModel.hh.

Constructor & Destructor Documentation

◆ G4JAEAPolarizedElasticScatteringModel() [1/2]

G4JAEAPolarizedElasticScatteringModel::G4JAEAPolarizedElasticScatteringModel ( )
explicit

Definition at line 49 of file G4JAEAPolarizedElasticScatteringModel.cc.

50 :G4VEmModel("G4JAEAPolarizedElasticScatteringModel"),isInitialised(false)
51{
52 fParticleChange = 0;
53 lowEnergyLimit = 100 * keV; //low energy limit for JAEAElasticScattering cross section data
54 fLinearPolarizationSensitvity1=1;
55 fLinearPolarizationSensitvity2=1;
56 fCircularPolarizationSensitvity=1;
57
58 verboseLevel= 0;
59 // Verbosity scale for debugging purposes:
60 // 0 = nothing
61 // 1 = calculation of cross sections, file openings...
62 // 2 = entering in methods
63
64 if(verboseLevel > 0)
65 {
66 G4cout << "G4JAEAPolarizedElasticScatteringModel is constructed " << G4endl;
67 }
68
69}
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4VEmModel::G4VEmModel
G4VEmModel(const G4String &nam)
Definition G4VEmModel.cc:67

◆ ~G4JAEAPolarizedElasticScatteringModel()

G4JAEAPolarizedElasticScatteringModel::~G4JAEAPolarizedElasticScatteringModel ( )
virtual

Definition at line 73 of file G4JAEAPolarizedElasticScatteringModel.cc.

74{
75 if(IsMaster()) {
76 for(G4int i=0; i<=maxZ; ++i) {
77 if(dataCS[i]) {
78 delete dataCS[i];
79 dataCS[i] = nullptr;
80 }
81 if (Polarized_ES_Data[i]){
82 delete Polarized_ES_Data[i];
83 Polarized_ES_Data[i] = nullptr;
84 }
85 }
86 }
87}
G4int
int G4int
Definition G4Types.hh:85
G4VEmModel::IsMaster
G4bool IsMaster() const
Definition G4VEmModel.hh:717

◆ G4JAEAPolarizedElasticScatteringModel() [2/2]

G4JAEAPolarizedElasticScatteringModel::G4JAEAPolarizedElasticScatteringModel ( const G4JAEAPolarizedElasticScatteringModel & )
delete

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

G4double G4JAEAPolarizedElasticScatteringModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition * ,
G4double kinEnergy,
G4double Z,
G4double A = 0,
G4double cut = 0,
G4double emax = DBL_MAX )
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 212 of file G4JAEAPolarizedElasticScatteringModel.cc.

217{
218 //Select the energy-grid point closest to the photon energy
219 // G4double *whichenergy = lower_bound(ESdata[0],ESdata[0]+300,GammaEnergy);
220 // int energyindex = max(0,(int)(whichenergy-ESdata[0]-1));
221
222 if (verboseLevel > 1)
223 {
224 G4cout << "G4JAEAPolarizedElasticScatteringModel::ComputeCrossSectionPerAtom()"
225 << G4endl;
226 }
227
228 if(GammaEnergy < lowEnergyLimit) { return 0.0; }
229
230 G4double xs = 0.0;
231
232 G4int intZ = G4lrint(Z);
233
234 if(intZ < 1 || intZ > maxZ) { return xs; }
235
236 G4PhysicsFreeVector* pv = dataCS[intZ];
237
238 // if element was not initialised
239 // do initialisation safely for MT mode
240 if(!pv) {
241 InitialiseForElement(0, intZ);
242 pv = dataCS[intZ];
243 if(!pv) { return xs; }
244 }
245
246 std::size_t n = pv->GetVectorLength() - 1;
247
248 G4double e = GammaEnergy;
249 if(e >= pv->Energy(n)) {
250 xs = (*pv)[n];
251 } else if(e >= pv->Energy(0)) {
252 xs = pv->Value(e);
253 }
254
255 if(verboseLevel > 0)
256 {
257 G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)="
258 << e << G4endl;
259 G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
260 G4cout << " -> first E*E*cs value in CS data file (iu) =" << (*pv)[0]
261 << G4endl;
262 G4cout << " -> last E*E*cs value in CS data file (iu) =" << (*pv)[n]
263 << G4endl;
264 G4cout << "*********************************************************"
265 << G4endl;
266 }
267
268 return (xs);
269}
G4double
double G4double
Definition G4Types.hh:83
G4JAEAPolarizedElasticScatteringModel::InitialiseForElement
void InitialiseForElement(const G4ParticleDefinition *, G4int Z) override
Definition G4JAEAPolarizedElasticScatteringModel.cc:589
G4PhysicsVector::Energy
G4double Energy(const std::size_t index) const
G4PhysicsVector::Value
G4double Value(const G4double energy, std::size_t &lastidx) const
G4PhysicsVector::GetVectorLength
std::size_t GetVectorLength() const
G4lrint
int G4lrint(double ad)
Definition templates.hh:134

◆ Initialise()

void G4JAEAPolarizedElasticScatteringModel::Initialise ( const G4ParticleDefinition * particle,
const G4DataVector & cuts )
overridevirtual

Implements G4VEmModel.

Definition at line 91 of file G4JAEAPolarizedElasticScatteringModel.cc.

93{
94 if (verboseLevel > 1)
95 {
96 G4cout << "Calling Initialise() of G4JAEAPolarizedElasticScatteringModel." << G4endl
97 << "Energy range: "
98 << LowEnergyLimit() / eV << " eV - "
99 << HighEnergyLimit() / GeV << " GeV"
100 << G4endl;
101 }
102
103 if(IsMaster()) {
104
105 // Initialise element selector
106 InitialiseElementSelectors(particle, cuts);
107
108 // Access to elements
109 const char* path = G4FindDataDir("G4LEDATA");
110 G4ProductionCutsTable* theCoupleTable =
111 G4ProductionCutsTable::GetProductionCutsTable();
112 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
113
114 for(G4int i=0; i<numOfCouples; ++i)
115 {
116 const G4MaterialCutsCouple* couple =
117 theCoupleTable->GetMaterialCutsCouple(i);
118 const G4Material* material = couple->GetMaterial();
119 const G4ElementVector* theElementVector = material->GetElementVector();
120 std::size_t nelm = material->GetNumberOfElements();
121
122 for (std::size_t j=0; j<nelm; ++j)
123 {
124 G4int Z = G4lrint((*theElementVector)[j]->GetZ());
125 if(Z < 1) { Z = 1; }
126 else if(Z > maxZ) { Z = maxZ; }
127 if( (!dataCS[Z]) ) { ReadData(Z, path); }
128 }
129 }
130 }
131
132 if(isInitialised) { return; }
133 fParticleChange = GetParticleChangeForGamma();
134 isInitialised = true;
135
136}
G4ElementVector
std::vector< const G4Element * > G4ElementVector
Definition G4ElementVector.hh:33
G4FindDataDir
const char * G4FindDataDir(const char *)
Definition G4FindDataDir.cc:94
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition G4MaterialCutsCouple.hh:166
G4Material::GetElementVector
const G4ElementVector * GetElementVector() const
Definition G4Material.hh:183
G4Material::GetNumberOfElements
std::size_t GetNumberOfElements() const
Definition G4Material.hh:180
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition G4ProductionCutsTable.hh:256
G4ProductionCutsTable::GetTableSize
std::size_t GetTableSize() const
Definition G4ProductionCutsTable.hh:249
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition G4ProductionCutsTable.cc:58
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition G4VEmModel.cc:120
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition G4VEmModel.hh:640
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition G4VEmModel.hh:633
G4VEmModel::InitialiseElementSelectors
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition G4VEmModel.cc:135

◆ InitialiseForElement()

void G4JAEAPolarizedElasticScatteringModel::InitialiseForElement ( const G4ParticleDefinition * ,
G4int Z )
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 589 of file G4JAEAPolarizedElasticScatteringModel.cc.

591{
592 G4AutoLock l(&G4JAEAPolarizedElasticScatteringModelMutex);
593 // G4cout << "G4JAEAPolarizedElasticScatteringModel::InitialiseForElement Z= "
594 // << Z << G4endl;
595 if(!dataCS[Z]) { ReadData(Z); }
596 l.unlock();
597}

Referenced by ComputeCrossSectionPerAtom().

◆ InitialiseLocal()

void G4JAEAPolarizedElasticScatteringModel::InitialiseLocal ( const G4ParticleDefinition * ,
G4VEmModel * masterModel )
overridevirtual

Reimplemented from G4VEmModel.

Definition at line 140 of file G4JAEAPolarizedElasticScatteringModel.cc.

142{
143 SetElementSelectors(masterModel->GetElementSelectors());
144}
G4VEmModel::SetElementSelectors
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition G4VEmModel.hh:816
G4VEmModel::GetElementSelectors
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition G4VEmModel.hh:808

◆ operator=()

G4JAEAPolarizedElasticScatteringModel & G4JAEAPolarizedElasticScatteringModel::operator= ( const G4JAEAPolarizedElasticScatteringModel & right)
delete

◆ SampleSecondaries()

void G4JAEAPolarizedElasticScatteringModel::SampleSecondaries ( std::vector< G4DynamicParticle * > * ,
const G4MaterialCutsCouple * couple,
const G4DynamicParticle * aDynamicGamma,
G4double tmin,
G4double maxEnergy )
overridevirtual

Implements G4VEmModel.

Definition at line 273 of file G4JAEAPolarizedElasticScatteringModel.cc.

278{
279 if (verboseLevel > 1) {
280
281 G4cout << "Calling SampleSecondaries() of G4JAEAPolarizedElasticScatteringModel."
282 << G4endl;
283 }
284 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
285
286 // absorption of low-energy gamma
287 if (photonEnergy0 <= lowEnergyLimit)
288 {
289 fParticleChange->ProposeTrackStatus(fStopAndKill);
290 fParticleChange->SetProposedKineticEnergy(0.);
291 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
292 return ;
293 }
294
295 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
296 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
297 G4int Z = G4lrint(elm->GetZ());
298
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*4*(energyindex));
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+181*4*(energyindex));
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;
309 }
310
311 CLHEP::RandGeneral GenThetaDist(distribution,180);
312 //Intial sampling of the scattering angle. To be updated for the circular polarization
313 G4double theta = CLHEP::pi*GenThetaDist.shoot();
314 //G4double theta =45.*CLHEP::pi/180.;
315 //Theta is in degree to call scattering amplitudes
316 G4int theta_in_degree =round(theta*180./CLHEP::pi);
317
318 //theta_in_degree=45;
319
320 G4double am1=0,am2=0,am3=0,am4=0,aparaSquare=0,aperpSquare=0,
321 apara_aper_Asterisk=0,img_apara_aper_Asterisk=0;
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_degree+1+300+181*4*(energyindex));
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_degree+3+300+181*4*(energyindex));
326 aparaSquare=am1*am1+am2*am2;
327 aperpSquare=am3*am3+am4*am4;
328 apara_aper_Asterisk=2*a1*a3+2*a2*a4;
329 img_apara_aper_Asterisk=2*a1*a4-2*a2*a3;
330
331 G4ThreeVector Direction_Unpolarized(0.,0.,0.);
332 G4ThreeVector Direction_Linear1(0.,0.,0.);
333 G4ThreeVector Direction_Linear2(0.,0.,0.);
334 G4ThreeVector Direction_Circular(0.,0.,0.);
335 G4ThreeVector Polarization_Unpolarized(0.,0.,0.);
336 G4ThreeVector Polarization_Linear1(0.,0.,0.);
337 G4ThreeVector Polarization_Linear2(0.,0.,0.);
338 G4ThreeVector Polarization_Circular(0.,0.,0.);
339
340 //Stokes parameters for the incoming and outgoing photon
341 G4double Xi1=0, Xi2=0, Xi3=0, Xi1_Prime=0,Xi2_Prime=0,Xi3_Prime=0;
342
343 //Getting the Stokes parameters for the incoming photon
344 G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
345 Xi1=gammaPolarization0.x();
346 Xi2=gammaPolarization0.y();
347 Xi3=gammaPolarization0.z();
348
349 //Polarization vector must be unit vector (5% tolerance)
350 if ((gammaPolarization0.mag())>1.05 || (Xi1*Xi1>1.05) || (Xi2*Xi2>1.05) || (Xi3*Xi3>1.05))
351 {
352 G4Exception("G4JAEAPolarizedElasticScatteringModel::SampleSecondaries()","em1006",
353 JustWarning,
354 "WARNING: G4JAEAPolarizedElasticScatteringModel is only compatible with a unit polarization vector.");
355 return;
356 }
357 //Unpolarized gamma rays
358 if (Xi1==0 && Xi2==0 && Xi3==0)
359 {
360 G4double Phi_Unpolarized=0;
361 if (fLinearPolarizationSensitvity1)
362 Phi_Unpolarized=GeneratePolarizedPhi(aparaSquare,aperpSquare,0.);
363 else
364 Phi_Unpolarized=CLHEP::twopi*G4UniformRand();
365 Direction_Unpolarized.setX(sin(theta)*cos(Phi_Unpolarized));
366 Direction_Unpolarized.setY(sin(theta)*sin(Phi_Unpolarized));
367 Direction_Unpolarized.setZ(cos(theta));
368 Direction_Unpolarized.rotateUz(aDynamicGamma->GetMomentumDirection());
369 Xi1_Prime=(aparaSquare-aperpSquare)/(aparaSquare+aperpSquare);
370 Polarization_Unpolarized.setX(Xi1_Prime);
371 Polarization_Unpolarized.setY(0.);
372 Polarization_Unpolarized.setZ(0.);
373 fParticleChange->ProposeMomentumDirection(Direction_Unpolarized);
374 fParticleChange->ProposePolarization(Polarization_Unpolarized);
375 return;
376 }
377
378 //Linear polarization defined by first Stokes parameter
379 G4double InitialAzimuth=aDynamicGamma->GetMomentumDirection().phi();
380 if(InitialAzimuth<0) InitialAzimuth=InitialAzimuth+CLHEP::twopi;
381
382 G4double Phi_Linear1=0.;
383
384 Phi_Linear1 = GeneratePolarizedPhi(aparaSquare+aperpSquare+Xi1*(aparaSquare-aperpSquare),
385 aparaSquare+aperpSquare-Xi1*(aparaSquare-aperpSquare),InitialAzimuth);
386
387 Xi1_Prime=((aparaSquare-aperpSquare)+Xi1*(aparaSquare+aperpSquare)*cos(2*Phi_Linear1))/
388 ((aparaSquare+aperpSquare)+Xi1*(aparaSquare-aperpSquare)*cos(2*Phi_Linear1));
389 Xi2_Prime=(-Xi1*apara_aper_Asterisk*sin(2*Phi_Linear1))/
390 ((aparaSquare+aperpSquare)+Xi1*(aparaSquare-aperpSquare)*cos(2*Phi_Linear1));
391 Xi3_Prime=(-Xi1*img_apara_aper_Asterisk*sin(2*Phi_Linear1))/
392 ((aparaSquare+aperpSquare)+Xi1*(aparaSquare-aperpSquare)*cos(2*Phi_Linear1));
393 //Store momentum direction and po;arization
394 Direction_Linear1.setX(sin(theta)*cos(Phi_Linear1));
395 Direction_Linear1.setY(sin(theta)*sin(Phi_Linear1));
396 Direction_Linear1.setZ(cos(theta));
397 Polarization_Linear1.setX(Xi1_Prime);
398 Polarization_Linear1.setY(Xi2_Prime);
399 Polarization_Linear1.setZ(Xi3_Prime);
400
401 //Set scattered photon polarization sensitivity
402 Xi1_Prime=Xi1_Prime*fLinearPolarizationSensitvity1;
403 Xi2_Prime=Xi2_Prime*fLinearPolarizationSensitvity2;
404 Xi3_Prime=Xi3_Prime*fCircularPolarizationSensitvity;
405
406 G4double dsigmaL1=0.0;
407 if(abs(Xi1)>0.0) dsigmaL1=0.25*((aparaSquare+aperpSquare)*
408 (1+Xi1*Xi1_Prime*cos(2*Phi_Linear1))+
409 (aparaSquare-aperpSquare)*(Xi1*cos(2*Phi_Linear1)+Xi1_Prime)
410 -Xi1*Xi2_Prime*apara_aper_Asterisk*sin(2*Phi_Linear1)-
411 Xi1*Xi3_Prime*img_apara_aper_Asterisk*sin(2*Phi_Linear1));
412
413 //Linear polarization defined by second Stokes parameter
414 //G4double IntialAzimuth=aDynamicGamma->GetMomentumDirection().phi();
415 G4double Phi_Linear2=0.;
416
417 InitialAzimuth=InitialAzimuth-CLHEP::pi/4.;
418 if(InitialAzimuth<0) InitialAzimuth=InitialAzimuth+CLHEP::twopi;
419
420 Phi_Linear2 = GeneratePolarizedPhi(aparaSquare+aperpSquare+Xi1*(aparaSquare-aperpSquare)
421 ,aparaSquare+aperpSquare-Xi1*(aparaSquare-aperpSquare),InitialAzimuth);
422
423 Xi1_Prime=((aparaSquare-aperpSquare)+Xi2*(aparaSquare+aperpSquare)*sin(2*Phi_Linear2))/
424 ((aparaSquare+aperpSquare)+Xi2*(aparaSquare-aperpSquare)*sin(2*Phi_Linear2));
425 Xi2_Prime=(Xi2*apara_aper_Asterisk*cos(2*Phi_Linear2))/
426 ((aparaSquare+aperpSquare)+Xi2*(aparaSquare-aperpSquare)*sin(2*Phi_Linear2));
427 Xi3_Prime=(Xi2*img_apara_aper_Asterisk*cos(2*Phi_Linear2))/
428 ((aparaSquare+aperpSquare)+Xi2*(aparaSquare-aperpSquare)*sin(2*Phi_Linear2));
429 //Store momentum direction and polarization
430 Direction_Linear2.setX(sin(theta)*cos(Phi_Linear2));
431 Direction_Linear2.setY(sin(theta)*sin(Phi_Linear2));
432 Direction_Linear2.setZ(cos(theta));
433 Polarization_Linear2.setX(Xi1_Prime);
434 Polarization_Linear2.setY(Xi2_Prime);
435 Polarization_Linear2.setZ(Xi3_Prime);
436
437 //Set scattered photon polarization sensitivity
438 Xi1_Prime=Xi1_Prime*fLinearPolarizationSensitvity1;
439 Xi2_Prime=Xi2_Prime*fLinearPolarizationSensitvity2;
440 Xi3_Prime=Xi3_Prime*fCircularPolarizationSensitvity;
441
442 G4double dsigmaL2=0.0;
443 if(abs(Xi2)>0.0)
444 dsigmaL2=0.25*((aparaSquare+aperpSquare)*(1+Xi2*Xi1_Prime*sin(2*Phi_Linear2))+
445 (aparaSquare-aperpSquare)*(Xi2*sin(2*Phi_Linear2)+Xi1_Prime)
446 +Xi2*Xi2_Prime*apara_aper_Asterisk*cos(2*Phi_Linear2)-
447 Xi2*Xi3_Prime*img_apara_aper_Asterisk*cos(2*Phi_Linear2));
448
449 //Circular polarization
450 G4double Phi_Circular = CLHEP::twopi*G4UniformRand();
451 G4double Theta_Circular = 0.;
452
453 Xi1_Prime=(aparaSquare-aperpSquare)/(aparaSquare+aperpSquare);
454 Xi2_Prime=(-Xi3*img_apara_aper_Asterisk)/(aparaSquare+aperpSquare);
455 Xi3_Prime=(Xi3*apara_aper_Asterisk)/(aparaSquare+aperpSquare);
456
457 Polarization_Circular.setX(Xi1_Prime);
458 Polarization_Circular.setY(Xi2_Prime);
459 Polarization_Circular.setZ(Xi3_Prime);
460
461 //Set scattered photon polarization sensitivity
462 Xi1_Prime=Xi1_Prime*fLinearPolarizationSensitvity1;
463 Xi2_Prime=Xi2_Prime*fLinearPolarizationSensitvity2;
464 Xi3_Prime=Xi3_Prime*fCircularPolarizationSensitvity;
465
466 G4double dsigmaC=0.0;
467 if(abs(Xi3)>0.0)
468 dsigmaC=0.25*(aparaSquare+aperpSquare+Xi1_Prime*(aparaSquare-aperpSquare)-
469 Xi3*Xi2_Prime*img_apara_aper_Asterisk
470 +Xi3*Xi3_Prime*apara_aper_Asterisk);
471
472 if (abs(Xi3)==0.0 && abs(Xi1_Prime)==0.0)
473 {
474 Direction_Circular.setX(sin(theta)*cos(Phi_Circular));
475 Direction_Circular.setY(sin(theta)*sin(Phi_Circular));
476 Direction_Circular.setZ(cos(theta));
477 }
478 else
479 {
480 G4double c1=0, c2=0, c3=0,c4=0;
481 for (G4int i=0;i<=180;++i)
482 {
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*4*(energyindex));
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*4*(energyindex));
487 cdistribution[i]=0.25*((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)
490 +Xi3*Xi3_Prime*(2*c1*c4-2*c2*c3));
491 }
492 CLHEP::RandGeneral GenTheta_Circ_Dist(cdistribution,180);
493 Theta_Circular=CLHEP::pi*GenTheta_Circ_Dist.shoot();
494 Direction_Circular.setX(sin(Theta_Circular)*cos(Phi_Circular));
495 Direction_Circular.setY(sin(Theta_Circular)*sin(Phi_Circular));
496 Direction_Circular.setZ(cos(Theta_Circular));
497 }
498
499 // Sampling scattered photon direction based on asymmetry arising from polarization mixing
500 G4double totalSigma= dsigmaL1+dsigmaL2+dsigmaC;
501 G4double prob1=dsigmaL1/totalSigma;
502 G4double prob2=dsigmaL2/totalSigma;
503 G4double probc=1-(prob1+prob2);
504
505 //Check the Probability of polarization mixing
506 if (abs(probc - dsigmaC/totalSigma)>=0.0001)
507 {
508 G4Exception("G4JAEAPolarizedElasticScatteringModel::SampleSecondaries()","em1007",
509 JustWarning,
510 "WARNING: Polarization mixing might be incorrect.");
511 }
512
513 // Generate outgoing photon direction
514 G4ThreeVector finaldirection(0.0,0.0,0.0);
515 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(Polarization_Linear1.x());
525 outcomingPhotonPolarization.setY(Polarization_Linear1.y());
526 outcomingPhotonPolarization.setZ(Polarization_Linear1.z());
527 }
528 else if ((polmix>prob1) && (polmix<=prob1+prob2))
529 {
530 finaldirection.setX(Direction_Linear2.x());
531 finaldirection.setY(Direction_Linear2.y());
532 finaldirection.setZ(Direction_Linear2.z());
533 outcomingPhotonPolarization.setX(Polarization_Linear2.x());
534 outcomingPhotonPolarization.setY(Polarization_Linear2.y());
535 outcomingPhotonPolarization.setZ(Polarization_Linear2.z());
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(Polarization_Circular.x());
543 outcomingPhotonPolarization.setY(Polarization_Circular.y());
544 outcomingPhotonPolarization.setZ(Polarization_Circular.z());
545 }
546
547 //Sampling the Final State
548 finaldirection.rotateUz(aDynamicGamma->GetMomentumDirection());
549 fParticleChange->ProposeMomentumDirection(finaldirection);
550 fParticleChange->SetProposedKineticEnergy(photonEnergy0);
551 fParticleChange->ProposePolarization(outcomingPhotonPolarization);
552
553}
JustWarning
@ JustWarning
Definition G4ExceptionSeverity.hh:71
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
fStopAndKill
@ fStopAndKill
Definition G4TrackStatus.hh:46
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::phi
double phi() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::mag
double mag() const
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4DynamicParticle::GetPolarization
const G4ThreeVector & GetPolarization() const
G4Element::GetZ
G4double GetZ() const
Definition G4Element.hh:119
G4ParticleChangeForGamma::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition G4ParticleChangeForGamma.hh:106
G4ParticleChangeForGamma::ProposePolarization
void ProposePolarization(const G4ThreeVector &dir)
Definition G4ParticleChangeForGamma.hh:132
G4ParticleChangeForGamma::ProposeMomentumDirection
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
Definition G4ParticleChangeForGamma.hh:120
G4VEmModel::SelectRandomAtom
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition G4VEmModel.hh:560
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)

◆ SetDebugVerbosity()

void G4JAEAPolarizedElasticScatteringModel::SetDebugVerbosity ( G4int val)
inline

Definition at line 76 of file G4JAEAPolarizedElasticScatteringModel.hh.

76{verboseLevel = val;};

◆ SetLowEnergyThreshold()

void G4JAEAPolarizedElasticScatteringModel::SetLowEnergyThreshold ( G4double val)
inline

Definition at line 73 of file G4JAEAPolarizedElasticScatteringModel.hh.

73{lowEnergyLimit = val;};

◆ SetPolarizationSensitvity()

void G4JAEAPolarizedElasticScatteringModel::SetPolarizationSensitvity ( G4bool linear1,
G4bool linear2,
G4bool circular )
inline

Definition at line 104 of file G4JAEAPolarizedElasticScatteringModel.hh.

107{
108 fLinearPolarizationSensitvity1=linear1;
109 fLinearPolarizationSensitvity2=linear2;
110 fCircularPolarizationSensitvity=circular;
111}
The documentation for this class was generated from the following files: