Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
G4PAIPhotData Class Reference

#include <G4PAIPhotData.hh>

Public Member Functions

 G4PAIPhotData (G4double tmin, G4double tmax, G4int verbose)
 
 ~G4PAIPhotData ()
 
void Initialise (const G4MaterialCutsCouple *, G4double cut, G4PAIPhotModel *)
 
G4double DEDXPerVolume (G4int coupleIndex, G4double scaledTkin, G4double cut) const
 
G4double CrossSectionPerVolume (G4int coupleIndex, G4double scaledTkin, G4double tcut, G4double tmax) const
 
G4double GetPlasmonRatio (G4int coupleIndex, G4double scaledTkin) const
 
G4double SampleAlongStepTransfer (G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
 
G4double SampleAlongStepPhotonTransfer (G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
 
G4double SampleAlongStepPlasmonTransfer (G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
 
G4double SamplePostStepTransfer (G4int coupleIndex, G4double scaledTkin) const
 
G4double SamplePostStepPhotonTransfer (G4int coupleIndex, G4double scaledTkin) const
 
G4double SamplePostStepPlasmonTransfer (G4int coupleIndex, G4double scaledTkin) const
 

Detailed Description

Definition at line 66 of file G4PAIPhotData.hh.

Constructor & Destructor Documentation

G4PAIPhotData::G4PAIPhotData ( G4double  tmin,
G4double  tmax,
G4int  verbose 
)

Definition at line 59 of file G4PAIPhotData.cc.

References G4cout, G4endl, python.hepunit::GeV, python.hepunit::keV, G4INCL::Math::max(), python.hepunit::MeV, and python.hepunit::TeV.

60 {
61  const G4int nPerDecade = 10;
62  const G4double lowestTkin = 50*keV;
63  const G4double highestTkin = 10*TeV;
64 
65  // fPAIxSection.SetVerbose(ver);
66 
67  fLowestKineticEnergy = std::max(tmin, lowestTkin);
68  fHighestKineticEnergy = tmax;
69 
70  if(tmax < 10*fLowestKineticEnergy)
71  {
72  fHighestKineticEnergy = 10*fLowestKineticEnergy;
73  }
74  else if(tmax > highestTkin)
75  {
76  fHighestKineticEnergy = std::max(highestTkin, 10*fLowestKineticEnergy);
77  }
78  fTotBin = (G4int)(nPerDecade*
79  std::log10(fHighestKineticEnergy/fLowestKineticEnergy));
80 
81  fParticleEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
82  fHighestKineticEnergy,
83  fTotBin);
84  if(0 < ver) {
85  G4cout << "### G4PAIPhotData: Nbins= " << fTotBin
86  << " Tmin(MeV)= " << fLowestKineticEnergy/MeV
87  << " Tmax(GeV)= " << fHighestKineticEnergy/GeV
88  << " tmin(keV)= " << tmin/keV << G4endl;
89  }
90 }
G4int
int G4int
Definition: G4Types.hh:78
G4cout
G4GLOB_DLL std::ostream G4cout
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:116
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4PAIPhotData::~G4PAIPhotData ( )

Definition at line 94 of file G4PAIPhotData.cc.

References n.

95 {
96  //delete fParticleEnergyVector;
97 
98  size_t n = fPAIxscBank.size();
99  if(0 < n)
100  {
101  for(size_t i=0; i<n; ++i)
102  {
103  if(fPAIxscBank[i])
104  {
105  // fPAIxscBank[i]->clearAndDestroy();
106  delete fPAIxscBank[i];
107  }
108  if(fPAIdEdxBank[i])
109  {
110  //fPAIdEdxBank[i]->clearAndDestroy();
111  delete fPAIdEdxBank[i];
112  }
113  //delete fdEdxTable[i];
114  //delete fdNdxCutTable[i];
115  }
116  }
117 }
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144

Member Function Documentation

G4double G4PAIPhotData::CrossSectionPerVolume ( G4int  coupleIndex,
G4double  scaledTkin,
G4double  tcut,
G4double  tmax 
) const

Definition at line 301 of file G4PAIPhotData.cc.

Referenced by G4PAIPhotModel::CrossSectionPerVolume().

304 {
305  G4double cross, xscEl, xscEl2, xscPh, xscPh2;
306 
307  cross=tcut+tmax;
308 
309  // iPlace is in interval from 0 to (N-1)
310 
311  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
312  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
313 
314  G4bool one = true;
315 
316  if( scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
317  else if( scaledTkin > fParticleEnergyVector->Energy(0) ) one = false;
318 
319 
320  xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
321  xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
322 
323  xscPh = xscPh2;
324  xscEl = xscEl2;
325 
326  cross = xscPh + xscEl;
327 
328  if( !one )
329  {
330  xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
331 
332  G4double E1 = fParticleEnergyVector->Energy(iPlace);
333  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
334 
335  G4double W = 1.0/(E2 - E1);
336  G4double W1 = (E2 - scaledTkin)*W;
337  G4double W2 = (scaledTkin - E1)*W;
338 
339  xscEl *= W1;
340  xscEl += W2*xscEl2;
341 
342  xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
343 
344  E1 = fParticleEnergyVector->Energy(iPlace);
345  E2 = fParticleEnergyVector->Energy(iPlace+1);
346 
347  W = 1.0/(E2 - E1);
348  W1 = (E2 - scaledTkin)*W;
349  W2 = (scaledTkin - E1)*W;
350 
351  xscPh *= W1;
352  xscPh += W2*xscPh2;
353 
354  cross = xscEl + xscPh;
355  }
356  if( cross < 0.0) cross = 0.0;
357 
358  return cross;
359 }
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
G4bool
bool G4bool
Definition: G4Types.hh:79
G4PhysicsVector::FindBin
size_t FindBin(G4double energy, size_t idx) const
G4PhysicsVector::Energy
G4double Energy(size_t index) const
G4double
double G4double
Definition: G4Types.hh:76
G4double G4PAIPhotData::DEDXPerVolume ( G4int  coupleIndex,
G4double  scaledTkin,
G4double  cut 
) const

Definition at line 265 of file G4PAIPhotData.cc.

Referenced by G4PAIPhotModel::ComputeDEDXPerVolume().

267 {
268  // VI: iPlace is the low edge index of the bin
269  // iPlace is in interval from 0 to (N-1)
270  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
271  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
272 
273  G4bool one = true;
274  if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
275  else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
276  one = false;
277  }
278 
279  // VI: apply interpolation of the vector
280  G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin);
281  G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
282  if(!one) {
283  G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
284  G4double E1 = fParticleEnergyVector->Energy(iPlace);
285  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
286  G4double W = 1.0/(E2 - E1);
287  G4double W1 = (E2 - scaledTkin)*W;
288  G4double W2 = (scaledTkin - E1)*W;
289  del *= W1;
290  del += W2*del2;
291  }
292  dEdx -= del;
293 
294  if( dEdx < 0.) { dEdx = 0.; }
295  return dEdx;
296 }
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
G4bool
bool G4bool
Definition: G4Types.hh:79
G4PhysicsVector::FindBin
size_t FindBin(G4double energy, size_t idx) const
G4PhysicsVector::Energy
G4double Energy(size_t index) const
G4double
double G4double
Definition: G4Types.hh:76
G4double G4PAIPhotData::GetPlasmonRatio ( G4int  coupleIndex,
G4double  scaledTkin 
) const

Definition at line 364 of file G4PAIPhotData.cc.

Referenced by G4PAIPhotModel::SampleSecondaries().

365 {
366  G4double cross, xscEl, xscEl2, xscPh, xscPh2, plRatio;
367  // iPlace is in interval from 0 to (N-1)
368 
369  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
370  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
371 
372  G4bool one = true;
373 
374  if( scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
375  else if( scaledTkin > fParticleEnergyVector->Energy(0) ) one = false;
376 
377 
378  xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
379  xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
380 
381  xscPh = xscPh2;
382  xscEl = xscEl2;
383 
384  cross = xscPh + xscEl;
385 
386  if( !one )
387  {
388  xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
389 
390  G4double E1 = fParticleEnergyVector->Energy(iPlace);
391  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
392 
393  G4double W = 1.0/(E2 - E1);
394  G4double W1 = (E2 - scaledTkin)*W;
395  G4double W2 = (scaledTkin - E1)*W;
396 
397  xscEl *= W1;
398  xscEl += W2*xscEl2;
399 
400  xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
401 
402  E1 = fParticleEnergyVector->Energy(iPlace);
403  E2 = fParticleEnergyVector->Energy(iPlace+1);
404 
405  W = 1.0/(E2 - E1);
406  W1 = (E2 - scaledTkin)*W;
407  W2 = (scaledTkin - E1)*W;
408 
409  xscPh *= W1;
410  xscPh += W2*xscPh2;
411 
412  cross = xscEl + xscPh;
413  }
414  if( cross <= 0.0)
415  {
416  plRatio = 2.0;
417  }
418  else
419  {
420  plRatio = xscEl/cross;
421 
422  if( plRatio > 1. || plRatio < 0.) plRatio = 2.0;
423  }
424  return plRatio;
425 }
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
G4bool
bool G4bool
Definition: G4Types.hh:79
G4PhysicsVector::FindBin
size_t FindBin(G4double energy, size_t idx) const
G4PhysicsVector::Energy
G4double Energy(size_t index) const
G4double
double G4double
Definition: G4Types.hh:76
void G4PAIPhotData::Initialise ( const G4MaterialCutsCouple *  couple,
G4double  cut,
G4PAIPhotModel *  model 
)

Definition at line 121 of file G4PAIPhotData.cc.

References G4PAIPhotModel::ComputeMaxEnergy(), python.hepunit::eV, G4cout, G4endl, G4ProductionCutsTable::GetEnergyCutsVector(), G4MaterialCutsCouple::GetMaterial(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), idxG4ElectronCut, idxG4GammaCut, G4PhysicsTable::insertAt(), python.hepunit::keV, n, python.hepunit::proton_mass_c2, G4PhysicsFreeVector::PutValue(), G4PhysicsVector::PutValue(), and G4PhysicsVector::Value().

Referenced by G4PAIPhotModel::Initialise().

123 {
124  G4ProductionCutsTable* theCoupleTable=
125  G4ProductionCutsTable::GetProductionCutsTable();
126  size_t numOfCouples = theCoupleTable->GetTableSize();
127  size_t jMatCC;
128 
129  for (jMatCC = 0; jMatCC < numOfCouples; jMatCC++ )
130  {
131  if( couple == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
132  }
133  if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
134 
135  const vector<G4double>* deltaCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
136  const vector<G4double>* photonCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
137  G4double deltaCutInKineticEnergyNow = (*deltaCutInKineticEnergy)[jMatCC];
138  G4double photonCutInKineticEnergyNow = (*photonCutInKineticEnergy)[jMatCC];
139 
140  G4cout<<"G4PAIPhotData::Initialise: "<<"cut = "<<cut/keV<<" keV; cutEl = "
141  <<deltaCutInKineticEnergyNow/keV<<" keV; cutPh = "
142  <<photonCutInKineticEnergyNow/keV<<" keV"<<G4endl;
143 
144  // if( deltaCutInKineticEnergyNow != cut ) deltaCutInKineticEnergyNow = cut; // exception??
145 
146  const G4Material* mat = couple->GetMaterial();
147  fSandia.Initialize(const_cast<G4Material*>(mat));
148 
149  G4PhysicsTable* PAItransferTable = new G4PhysicsTable(fTotBin+1);
150  G4PhysicsTable* PAIphotonTable = new G4PhysicsTable(fTotBin+1);
151  G4PhysicsTable* PAIplasmonTable = new G4PhysicsTable(fTotBin+1);
152 
153  G4PhysicsTable* PAIdEdxTable = new G4PhysicsTable(fTotBin+1);
154  G4PhysicsLogVector* dEdxCutVector =
155  new G4PhysicsLogVector(fLowestKineticEnergy,
156  fHighestKineticEnergy,
157  fTotBin);
158  G4PhysicsLogVector* dEdxMeanVector =
159  new G4PhysicsLogVector(fLowestKineticEnergy,
160  fHighestKineticEnergy,
161  fTotBin);
162 
163  G4PhysicsLogVector* dNdxCutVector =
164  new G4PhysicsLogVector(fLowestKineticEnergy,
165  fHighestKineticEnergy,
166  fTotBin);
167  G4PhysicsLogVector* dNdxCutPhotonVector =
168  new G4PhysicsLogVector(fLowestKineticEnergy,
169  fHighestKineticEnergy,
170  fTotBin);
171  G4PhysicsLogVector* dNdxCutPlasmonVector =
172  new G4PhysicsLogVector(fLowestKineticEnergy,
173  fHighestKineticEnergy,
174  fTotBin);
175 
176  // low energy Sandia interval
177  G4double Tmin = fSandia.GetSandiaMatTablePAI(0,0);
178 
179  // energy safety
180  const G4double deltaLow = 100.*eV;
181 
182  for (G4int i = 0; i <= fTotBin; ++i)
183  {
184  G4double kinEnergy = fParticleEnergyVector->Energy(i);
185  G4double Tmax = model->ComputeMaxEnergy(kinEnergy);
186  G4double tau = kinEnergy/proton_mass_c2;
187  G4double bg2 = tau*( tau + 2. );
188 
189  if ( Tmax < Tmin + deltaLow ) Tmax = Tmin + deltaLow;
190 
191  fPAIxSection.Initialize( mat, Tmax, bg2, &fSandia);
192 
193  //G4cout << i << ". TransferMax(keV)= "<< Tmax/keV << " cut(keV)= "
194  // << cut/keV << " E(MeV)= " << kinEnergy/MeV << G4endl;
195 
196  G4int n = fPAIxSection.GetSplineSize();
197 
198  G4PhysicsFreeVector* transferVector = new G4PhysicsFreeVector(n);
199  G4PhysicsFreeVector* photonVector = new G4PhysicsFreeVector(n);
200  G4PhysicsFreeVector* plasmonVector = new G4PhysicsFreeVector(n);
201 
202  G4PhysicsFreeVector* dEdxVector = new G4PhysicsFreeVector(n);
203 
204  for( G4int k = 0; k < n; k++ )
205  {
206  G4double t = fPAIxSection.GetSplineEnergy(k+1);
207 
208  transferVector->PutValue(k , t,
209  t*fPAIxSection.GetIntegralPAIxSection(k+1));
210  photonVector->PutValue(k , t,
211  t*fPAIxSection.GetIntegralCerenkov(k+1));
212  plasmonVector->PutValue(k , t,
213  t*fPAIxSection.GetIntegralPlasmon(k+1));
214 
215  dEdxVector->PutValue(k, t, fPAIxSection.GetIntegralPAIdEdx(k+1));
216  }
217  // G4cout << *transferVector << G4endl;
218 
219  G4double ionloss = fPAIxSection.GetMeanEnergyLoss();// total <dE/dx>
220 
221  if(ionloss < 0.0) ionloss = 0.0;
222 
223  dEdxMeanVector->PutValue(i,ionloss);
224 
225  G4double dNdxCut = transferVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
226  G4double dNdxCutPhoton = photonVector->Value(photonCutInKineticEnergyNow)/photonCutInKineticEnergyNow;
227  G4double dNdxCutPlasmon = plasmonVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
228 
229  G4double dEdxCut = dEdxVector->Value(cut)/cut;
230  //G4cout << "i= " << i << " x= " << dNdxCut << G4endl;
231 
232  if(dNdxCut < 0.0) { dNdxCut = 0.0; }
233  if(dNdxCutPhoton < 0.0) { dNdxCutPhoton = 0.0; }
234  if(dNdxCutPlasmon < 0.0) { dNdxCutPlasmon = 0.0; }
235 
236  dNdxCutVector->PutValue(i, dNdxCut);
237  dNdxCutPhotonVector->PutValue(i, dNdxCutPhoton);
238  dNdxCutPlasmonVector->PutValue(i, dNdxCutPlasmon);
239 
240  dEdxCutVector->PutValue(i, dEdxCut);
241 
242  PAItransferTable->insertAt(i,transferVector);
243  PAIphotonTable->insertAt(i,photonVector);
244  PAIplasmonTable->insertAt(i,plasmonVector);
245  PAIdEdxTable->insertAt(i,dEdxVector);
246 
247  } // end of Tkin loop
248 
249  fPAIxscBank.push_back(PAItransferTable);
250  fPAIphotonBank.push_back(PAIphotonTable);
251  fPAIplasmonBank.push_back(PAIplasmonTable);
252 
253  fPAIdEdxBank.push_back(PAIdEdxTable);
254  fdEdxTable.push_back(dEdxMeanVector);
255 
256  fdNdxCutTable.push_back(dNdxCutVector);
257  fdNdxCutPhotonTable.push_back(dNdxCutPhotonVector);
258  fdNdxCutPlasmonTable.push_back(dNdxCutPlasmonVector);
259 
260  fdEdxCutTable.push_back(dEdxCutVector);
261 }
G4ProductionCutsTable::GetEnergyCutsVector
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
Definition: G4ProductionCutsTable.hh:250
G4PhysicsFreeVector::PutValue
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
Definition: G4PhysicsFreeVector.cc:94
G4PAIxSection::GetIntegralPlasmon
G4double GetIntegralPlasmon(G4int i) const
Definition: G4PAIxSection.hh:336
G4PAIxSection::GetIntegralPAIxSection
G4double GetIntegralPAIxSection(G4int i) const
Definition: G4PAIxSection.hh:312
idxG4GammaCut
Definition: G4ProductionCuts.hh:53
G4SandiaTable::GetSandiaMatTablePAI
G4double GetSandiaMatTablePAI(G4int, G4int)
Definition: G4SandiaTable.cc:986
G4PAIxSection::GetIntegralPAIdEdx
G4double GetIntegralPAIdEdx(G4int i) const
Definition: G4PAIxSection.hh:318
G4int
int G4int
Definition: G4Types.hh:78
idxG4ElectronCut
Definition: G4ProductionCuts.hh:54
G4cout
G4GLOB_DLL std::ostream G4cout
G4PAIPhotModel::ComputeMaxEnergy
G4double ComputeMaxEnergy(G4double scaledEnergy)
Definition: G4PAIPhotModel.hh:153
G4PAIxSection::GetSplineSize
G4int GetSplineSize() const
Definition: G4PAIxSection.hh:187
G4PAIxSection::GetIntegralCerenkov
G4double GetIntegralCerenkov(G4int i) const
Definition: G4PAIxSection.hh:324
G4PhysicsVector::PutValue
void PutValue(size_t index, G4double theValue)
python.hepunit.proton_mass_c2
float proton_mass_c2
Definition: hepunit.py:275
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4PhysicsVector::Energy
G4double Energy(size_t index) const
G4PhysicsVector::Value
G4double Value(G4double theEnergy, size_t &lastidx) const
Definition: G4PhysicsVector.cc:511
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:63
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition: G4ProductionCutsTable.hh:262
G4PAIxSection::Initialize
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
Definition: G4PAIxSection.cc:606
G4PhysicsTable::insertAt
void insertAt(size_t, G4PhysicsVector *)
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4PAIxSection::GetMeanEnergyLoss
G4double GetMeanEnergyLoss() const
Definition: G4PAIxSection.hh:199
G4SandiaTable::Initialize
void Initialize(G4Material *)
Definition: G4SandiaTable.cc:625
G4PAIxSection::GetSplineEnergy
G4double GetSplineEnergy(G4int i) const
Definition: G4PAIxSection.hh:306
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150
G4double G4PAIPhotData::SampleAlongStepPhotonTransfer ( G4int  coupleIndex,
G4double  kinEnergy,
G4double  scaledTkin,
G4double  stepFactor 
) const

Definition at line 521 of file G4PAIPhotData.cc.

References G4PhysicsVector::Energy(), G4Poisson(), G4UniformRand, and position.

Referenced by G4PAIPhotModel::SampleFluctuations().

525 {
526  G4double loss = 0.0;
527  G4double omega;
528  G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
529 
530  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
531  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
532 
533  G4bool one = true;
534 
535  if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
536  else if(scaledTkin > fParticleEnergyVector->Energy(0)) one = false;
537 
538  G4PhysicsLogVector* vcut = fdNdxCutPhotonTable[coupleIndex];
539  G4PhysicsVector* v1 = (*(fPAIphotonBank[coupleIndex]))(iPlace);
540  G4PhysicsVector* v2 = 0;
541 
542  dNdxCut1 = (*vcut)[iPlace];
543  G4double e1 = v1->Energy(0);
544  G4double e2 = e1;
545 
546  G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
547 
548  meanNumber = meanN1;
549 
550  // G4cout<<"iPlace = "<<iPlace<< " meanN1= " << meanN1
551  // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " dNdxCut1= " << dNdxCut1 << G4endl;
552 
553  dNdxCut2 = dNdxCut1;
554  W1 = 1.0;
555  W2 = 0.0;
556  if(!one)
557  {
558  v2 = (*(fPAIphotonBank[coupleIndex]))(iPlace+1);
559  dNdxCut2 = (*vcut)[iPlace+1];
560  e2 = v2->Energy(0);
561 
562  G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
563 
564  E1 = fParticleEnergyVector->Energy(iPlace);
565  E2 = fParticleEnergyVector->Energy(iPlace+1);
566  W = 1.0/(E2 - E1);
567  W1 = (E2 - scaledTkin)*W;
568  W2 = (scaledTkin - E1)*W;
569  meanNumber = W1*meanN1 + W2*meanN2;
570 
571  //G4cout<<"meanN= " << meanNumber << " meanN2= " << meanN2
572  // << " dNdxCut2= " << dNdxCut2 << G4endl;
573  }
574  if( meanNumber <= 0.0) return 0.0;
575 
576  G4int numOfCollisions = G4Poisson(meanNumber);
577 
578  //G4cout << "N= " << numOfCollisions << G4endl;
579 
580  if( 0 == numOfCollisions) return 0.0;
581 
582  for(G4int i=0; i< numOfCollisions; ++i)
583  {
584  G4double rand = G4UniformRand();
585  position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
586  omega = GetEnergyPhotonTransfer(coupleIndex, iPlace, position);
587 
588  //G4cout << "omega(keV)= " << omega/keV << G4endl;
589 
590  if(!one)
591  {
592  position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
593  G4double omega2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);
594  omega = omega*W1 + omega2*W2;
595  }
596  //G4cout << "omega(keV)= " << omega/keV << G4endl;
597 
598  loss += omega;
599  if( loss > kinEnergy) { break; }
600  }
601 
602  // G4cout<<"PAIPhotData AlongStepLoss = "<<loss/keV<<" keV, on step = "
603  //<<step/mm<<" mm"<<G4endl;
604 
605  if ( loss > kinEnergy) loss = kinEnergy;
606  else if( loss < 0.) loss = 0.;
607 
608  return loss;
609 }
G4Poisson
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4bool
bool G4bool
Definition: G4Types.hh:79
G4PhysicsVector::FindBin
size_t FindBin(G4double energy, size_t idx) const
G4PhysicsVector::Energy
G4double Energy(size_t index) const
position
int position
Definition: filter.cc:7
G4double
double G4double
Definition: G4Types.hh:76
G4double G4PAIPhotData::SampleAlongStepPlasmonTransfer ( G4int  coupleIndex,
G4double  kinEnergy,
G4double  scaledTkin,
G4double  stepFactor 
) const

Definition at line 613 of file G4PAIPhotData.cc.

References G4PhysicsVector::Energy(), G4Poisson(), G4UniformRand, and position.

Referenced by G4PAIPhotModel::SampleFluctuations().

617 {
618  G4double loss = 0.0;
619  G4double omega;
620  G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
621 
622  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
623  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
624 
625  G4bool one = true;
626 
627  if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
628  else if(scaledTkin > fParticleEnergyVector->Energy(0)) one = false;
629 
630  G4PhysicsLogVector* vcut = fdNdxCutPlasmonTable[coupleIndex];
631  G4PhysicsVector* v1 = (*(fPAIplasmonBank[coupleIndex]))(iPlace);
632  G4PhysicsVector* v2 = 0;
633 
634  dNdxCut1 = (*vcut)[iPlace];
635  G4double e1 = v1->Energy(0);
636  G4double e2 = e1;
637 
638  G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
639 
640  meanNumber = meanN1;
641 
642  // G4cout<<"iPlace = "<<iPlace<< " meanN1= " << meanN1
643  // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " dNdxCut1= " << dNdxCut1 << G4endl;
644 
645  dNdxCut2 = dNdxCut1;
646  W1 = 1.0;
647  W2 = 0.0;
648  if(!one)
649  {
650  v2 = (*(fPAIplasmonBank[coupleIndex]))(iPlace+1);
651  dNdxCut2 = (*vcut)[iPlace+1];
652  e2 = v2->Energy(0);
653 
654  G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
655 
656  E1 = fParticleEnergyVector->Energy(iPlace);
657  E2 = fParticleEnergyVector->Energy(iPlace+1);
658  W = 1.0/(E2 - E1);
659  W1 = (E2 - scaledTkin)*W;
660  W2 = (scaledTkin - E1)*W;
661  meanNumber = W1*meanN1 + W2*meanN2;
662 
663  //G4cout<<"meanN= " << meanNumber << " meanN2= " << meanN2
664  // << " dNdxCut2= " << dNdxCut2 << G4endl;
665  }
666  if( meanNumber <= 0.0) return 0.0;
667 
668  G4int numOfCollisions = G4Poisson(meanNumber);
669 
670  //G4cout << "N= " << numOfCollisions << G4endl;
671 
672  if( 0 == numOfCollisions) return 0.0;
673 
674  for(G4int i=0; i< numOfCollisions; ++i)
675  {
676  G4double rand = G4UniformRand();
677  position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
678  omega = GetEnergyPlasmonTransfer(coupleIndex, iPlace, position);
679 
680  //G4cout << "omega(keV)= " << omega/keV << G4endl;
681 
682  if(!one)
683  {
684  position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
685  G4double omega2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1, position);
686  omega = omega*W1 + omega2*W2;
687  }
688  //G4cout << "omega(keV)= " << omega/keV << G4endl;
689 
690  loss += omega;
691  if( loss > kinEnergy) { break; }
692  }
693 
694  // G4cout<<"PAIPhotData AlongStepLoss = "<<loss/keV<<" keV, on step = "
695  //<<step/mm<<" mm"<<G4endl;
696 
697  if ( loss > kinEnergy) loss = kinEnergy;
698  else if( loss < 0.) loss = 0.;
699 
700  return loss;
701 }
G4Poisson
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4bool
bool G4bool
Definition: G4Types.hh:79
G4PhysicsVector::FindBin
size_t FindBin(G4double energy, size_t idx) const
G4PhysicsVector::Energy
G4double Energy(size_t index) const
position
int position
Definition: filter.cc:7
G4double
double G4double
Definition: G4Types.hh:76
G4double G4PAIPhotData::SampleAlongStepTransfer ( G4int  coupleIndex,
G4double  kinEnergy,
G4double  scaledTkin,
G4double  stepFactor 
) const

Definition at line 429 of file G4PAIPhotData.cc.

References G4PhysicsVector::Energy(), G4Poisson(), G4UniformRand, and position.

433 {
434  G4double loss = 0.0;
435  G4double omega;
436  G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
437 
438  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
439  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
440 
441  G4bool one = true;
442 
443  if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
444  else if(scaledTkin > fParticleEnergyVector->Energy(0)) one = false;
445 
446  G4PhysicsLogVector* vcut = fdNdxCutTable[coupleIndex];
447  G4PhysicsVector* v1 = (*(fPAIxscBank[coupleIndex]))(iPlace);
448  G4PhysicsVector* v2 = 0;
449 
450  dNdxCut1 = (*vcut)[iPlace];
451  G4double e1 = v1->Energy(0);
452  G4double e2 = e1;
453 
454  G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
455 
456  meanNumber = meanN1;
457 
458  // G4cout<<"iPlace = "<<iPlace<< " meanN1= " << meanN1
459  // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " dNdxCut1= " << dNdxCut1 << G4endl;
460 
461  dNdxCut2 = dNdxCut1;
462  W1 = 1.0;
463  W2 = 0.0;
464  if(!one)
465  {
466  v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
467  dNdxCut2 = (*vcut)[iPlace+1];
468  e2 = v2->Energy(0);
469 
470  G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
471 
472  E1 = fParticleEnergyVector->Energy(iPlace);
473  E2 = fParticleEnergyVector->Energy(iPlace+1);
474  W = 1.0/(E2 - E1);
475  W1 = (E2 - scaledTkin)*W;
476  W2 = (scaledTkin - E1)*W;
477  meanNumber = W1*meanN1 + W2*meanN2;
478 
479  //G4cout<<"meanN= " << meanNumber << " meanN2= " << meanN2
480  // << " dNdxCut2= " << dNdxCut2 << G4endl;
481  }
482  if( meanNumber <= 0.0) return 0.0;
483 
484  G4int numOfCollisions = G4Poisson(meanNumber);
485 
486  //G4cout << "N= " << numOfCollisions << G4endl;
487 
488  if( 0 == numOfCollisions) return 0.0;
489 
490  for(G4int i=0; i< numOfCollisions; ++i)
491  {
492  G4double rand = G4UniformRand();
493  position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
494  omega = GetEnergyTransfer(coupleIndex, iPlace, position);
495 
496  //G4cout << "omega(keV)= " << omega/keV << G4endl;
497 
498  if(!one)
499  {
500  position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
501  G4double omega2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
502  omega = omega*W1 + omega2*W2;
503  }
504  //G4cout << "omega(keV)= " << omega/keV << G4endl;
505 
506  loss += omega;
507  if( loss > kinEnergy) { break; }
508  }
509 
510  // G4cout<<"PAIPhotData AlongStepLoss = "<<loss/keV<<" keV, on step = "
511  //<<step/mm<<" mm"<<G4endl;
512 
513  if ( loss > kinEnergy) loss = kinEnergy;
514  else if( loss < 0.) loss = 0.;
515 
516  return loss;
517 }
G4Poisson
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4bool
bool G4bool
Definition: G4Types.hh:79
G4PhysicsVector::FindBin
size_t FindBin(G4double energy, size_t idx) const
G4PhysicsVector::Energy
G4double Energy(size_t index) const
position
int position
Definition: filter.cc:7
G4double
double G4double
Definition: G4Types.hh:76
G4double G4PAIPhotData::SamplePostStepPhotonTransfer ( G4int  coupleIndex,
G4double  scaledTkin 
) const

Definition at line 764 of file G4PAIPhotData.cc.

References G4UniformRand, and position.

Referenced by G4PAIPhotModel::SampleSecondaries().

766 {
767  //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
768  G4double transfer = 0.0;
769  G4double rand = G4UniformRand();
770 
771  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
772 
773  // size_t iTransfer, iTr1, iTr2;
774  G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
775 
776  G4PhysicsVector* cutv = fdNdxCutPhotonTable[coupleIndex];
777 
778  // Fermi plato, try from left
779 
780  if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
781  {
782  position = (*cutv)[nPlace]*rand;
783  transfer = GetEnergyPhotonTransfer(coupleIndex, nPlace, position);
784  }
785  else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
786  {
787  position = (*cutv)[0]*rand;
788  transfer = GetEnergyPhotonTransfer(coupleIndex, 0, position);
789  }
790  else
791  {
792  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
793 
794  dNdxCut1 = (*cutv)[iPlace];
795  dNdxCut2 = (*cutv)[iPlace+1];
796 
797  E1 = fParticleEnergyVector->Energy(iPlace);
798  E2 = fParticleEnergyVector->Energy(iPlace+1);
799  W = 1.0/(E2 - E1);
800  W1 = (E2 - scaledTkin)*W;
801  W2 = (scaledTkin - E1)*W;
802 
803  //G4cout<<"iPlace= " << " dNdxCut1 = "<<dNdxCut1
804  // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;
805 
806  position = dNdxCut1*rand;
807 
808  G4double tr1 = GetEnergyPhotonTransfer(coupleIndex, iPlace, position);
809 
810  position = dNdxCut2*rand;
811  G4double tr2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);
812 
813  transfer = tr1*W1 + tr2*W2;
814  }
815  //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl;
816  if(transfer < 0.0 ) { transfer = 0.0; }
817  return transfer;
818 }
G4PhysicsVector::GetMaxEnergy
G4double GetMaxEnergy() const
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4PhysicsVector::FindBin
size_t FindBin(G4double energy, size_t idx) const
G4PhysicsVector::Energy
G4double Energy(size_t index) const
position
int position
Definition: filter.cc:7
G4double
double G4double
Definition: G4Types.hh:76
G4double G4PAIPhotData::SamplePostStepPlasmonTransfer ( G4int  coupleIndex,
G4double  scaledTkin 
) const

Definition at line 822 of file G4PAIPhotData.cc.

References G4UniformRand, and position.

Referenced by G4PAIPhotModel::SampleSecondaries().

824 {
825  //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
826  G4double transfer = 0.0;
827  G4double rand = G4UniformRand();
828 
829  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
830 
831  // size_t iTransfer, iTr1, iTr2;
832  G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
833 
834  G4PhysicsVector* cutv = fdNdxCutPlasmonTable[coupleIndex];
835 
836  // Fermi plato, try from left
837  if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
838  {
839  position = (*cutv)[nPlace]*rand;
840  transfer = GetEnergyPlasmonTransfer(coupleIndex, nPlace, position);
841  }
842  else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
843  {
844  position = (*cutv)[0]*rand;
845  transfer = GetEnergyPlasmonTransfer(coupleIndex, 0, position);
846  }
847  else
848  {
849  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
850 
851  dNdxCut1 = (*cutv)[iPlace];
852  dNdxCut2 = (*cutv)[iPlace+1];
853 
854  E1 = fParticleEnergyVector->Energy(iPlace);
855  E2 = fParticleEnergyVector->Energy(iPlace+1);
856  W = 1.0/(E2 - E1);
857  W1 = (E2 - scaledTkin)*W;
858  W2 = (scaledTkin - E1)*W;
859 
860  //G4cout<<"iPlace= " << " dNdxCut1 = "<<dNdxCut1
861  // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;
862 
863  position = dNdxCut1*rand;
864  G4double tr1 = GetEnergyPlasmonTransfer(coupleIndex, iPlace, position);
865 
866  position = dNdxCut2*rand;
867  G4double tr2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1, position);
868 
869  transfer = tr1*W1 + tr2*W2;
870  }
871  //G4cout<<"PAImodel PostStepPlasmonTransfer = "<<transfer/keV<<" keV"<<G4endl;
872 
873  if(transfer < 0.0 ) transfer = 0.0;
874 
875  return transfer;
876 }
G4PhysicsVector::GetMaxEnergy
G4double GetMaxEnergy() const
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4PhysicsVector::FindBin
size_t FindBin(G4double energy, size_t idx) const
G4PhysicsVector::Energy
G4double Energy(size_t index) const
position
int position
Definition: filter.cc:7
G4double
double G4double
Definition: G4Types.hh:76
G4double G4PAIPhotData::SamplePostStepTransfer ( G4int  coupleIndex,
G4double  scaledTkin 
) const

Definition at line 708 of file G4PAIPhotData.cc.

References G4UniformRand, and position.

710 {
711  //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
712  G4double transfer = 0.0;
713  G4double rand = G4UniformRand();
714 
715  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
716 
717  // size_t iTransfer, iTr1, iTr2;
718  G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
719 
720  G4PhysicsVector* cutv = fdNdxCutTable[coupleIndex];
721 
722  // Fermi plato, try from left
723  if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
724  {
725  position = (*cutv)[nPlace]*rand;
726  transfer = GetEnergyTransfer(coupleIndex, nPlace, position);
727  }
728  else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
729  {
730  position = (*cutv)[0]*rand;
731  transfer = GetEnergyTransfer(coupleIndex, 0, position);
732  }
733  else
734  {
735  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
736 
737  dNdxCut1 = (*cutv)[iPlace];
738  dNdxCut2 = (*cutv)[iPlace+1];
739 
740  E1 = fParticleEnergyVector->Energy(iPlace);
741  E2 = fParticleEnergyVector->Energy(iPlace+1);
742  W = 1.0/(E2 - E1);
743  W1 = (E2 - scaledTkin)*W;
744  W2 = (scaledTkin - E1)*W;
745 
746  //G4cout<<"iPlace= " << " dNdxCut1 = "<<dNdxCut1
747  // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;
748 
749  position = dNdxCut1*rand;
750  G4double tr1 = GetEnergyTransfer(coupleIndex, iPlace, position);
751 
752  position = dNdxCut2*rand;
753  G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
754 
755  transfer = tr1*W1 + tr2*W2;
756  }
757  //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl;
758  if(transfer < 0.0 ) { transfer = 0.0; }
759  return transfer;
760 }
G4PhysicsVector::GetMaxEnergy
G4double GetMaxEnergy() const
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4PhysicsVector::FindBin
size_t FindBin(G4double energy, size_t idx) const
G4PhysicsVector::Energy
G4double Energy(size_t index) const
position
int position
Definition: filter.cc:7
G4double
double G4double
Definition: G4Types.hh:76
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:56:51 for Geant4.10 by   doxygen 1.8.7