Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | List of all members
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 
)
explicit

Definition at line 59 of file G4PAIPhotData.cc.

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
TeV
static constexpr double TeV
Definition: G4SIunits.hh:218
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:112
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double
double G4double
Definition: G4Types.hh:76
keV
static constexpr double keV
Definition: G4SIunits.hh:216

Here is the call graph for this function:

G4PAIPhotData::~G4PAIPhotData ( )

Definition at line 94 of file G4PAIPhotData.cc.

95 {
96  //G4cout << "G4PAIPhotData::~G4PAIPhotData() " << this << G4endl;
97  size_t n = fPAIxscBank.size();
98  if(0 < n)
99  {
100  for(size_t i=0; i<n; ++i)
101  {
102  if(fPAIxscBank[i])
103  {
104  fPAIxscBank[i]->clearAndDestroy();
105  delete fPAIxscBank[i];
106  fPAIxscBank[i] = 0;
107  }
108  if(fPAIdEdxBank[i])
109  {
110  fPAIdEdxBank[i]->clearAndDestroy();
111  delete fPAIdEdxBank[i];
112  fPAIdEdxBank[i]= 0;
113  }
114  delete fdEdxTable[i];
115  delete fdNdxCutTable[i];
116  fdEdxTable[i] = 0;
117  fdNdxCutTable[i] = 0;
118  }
119  }
120  delete fParticleEnergyVector;
121  fParticleEnergyVector = 0;
122  //G4cout << "G4PAIPhotData::~G4PAIPhotData() done for " << this << G4endl;
123 }

Member Function Documentation

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

Definition at line 308 of file G4PAIPhotData.cc.

311 {
312  G4double cross, xscEl, xscEl2, xscPh, xscPh2;
313 
314  cross=tcut+tmax;
315 
316  // iPlace is in interval from 0 to (N-1)
317 
318  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
319  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
320 
321  G4bool one = true;
322 
323  if( scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
324  else if( scaledTkin > fParticleEnergyVector->Energy(0) ) one = false;
325 
326 
327  xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
328  xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
329 
330  xscPh = xscPh2;
331  xscEl = xscEl2;
332 
333  cross = xscPh + xscEl;
334 
335  if( !one )
336  {
337  xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
338 
339  G4double E1 = fParticleEnergyVector->Energy(iPlace);
340  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
341 
342  G4double W = 1.0/(E2 - E1);
343  G4double W1 = (E2 - scaledTkin)*W;
344  G4double W2 = (scaledTkin - E1)*W;
345 
346  xscEl *= W1;
347  xscEl += W2*xscEl2;
348 
349  xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
350 
351  E1 = fParticleEnergyVector->Energy(iPlace);
352  E2 = fParticleEnergyVector->Energy(iPlace+1);
353 
354  W = 1.0/(E2 - E1);
355  W1 = (E2 - scaledTkin)*W;
356  W2 = (scaledTkin - E1)*W;
357 
358  xscPh *= W1;
359  xscPh += W2*xscPh2;
360 
361  cross = xscEl + xscPh;
362  }
363  if( cross < 0.0) cross = 0.0;
364 
365  return cross;
366 }
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

Here is the caller graph for this function:

G4double G4PAIPhotData::DEDXPerVolume ( G4int  coupleIndex,
G4double  scaledTkin,
G4double  cut 
) const

Definition at line 272 of file G4PAIPhotData.cc.

274 {
275  // VI: iPlace is the low edge index of the bin
276  // iPlace is in interval from 0 to (N-1)
277  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
278  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
279 
280  G4bool one = true;
281  if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
282  else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
283  one = false;
284  }
285 
286  // VI: apply interpolation of the vector
287  G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin);
288  G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
289  if(!one) {
290  G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
291  G4double E1 = fParticleEnergyVector->Energy(iPlace);
292  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
293  G4double W = 1.0/(E2 - E1);
294  G4double W1 = (E2 - scaledTkin)*W;
295  G4double W2 = (scaledTkin - E1)*W;
296  del *= W1;
297  del += W2*del2;
298  }
299  dEdx -= del;
300 
301  if( dEdx < 0.) { dEdx = 0.; }
302  return dEdx;
303 }
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

Here is the caller graph for this function:

G4double G4PAIPhotData::GetPlasmonRatio ( G4int  coupleIndex,
G4double  scaledTkin 
) const

Definition at line 371 of file G4PAIPhotData.cc.

372 {
373  G4double cross, xscEl, xscEl2, xscPh, xscPh2, plRatio;
374  // iPlace is in interval from 0 to (N-1)
375 
376  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
377  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
378 
379  G4bool one = true;
380 
381  if( scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
382  else if( scaledTkin > fParticleEnergyVector->Energy(0) ) one = false;
383 
384 
385  xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
386  xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
387 
388  xscPh = xscPh2;
389  xscEl = xscEl2;
390 
391  cross = xscPh + xscEl;
392 
393  if( !one )
394  {
395  xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
396 
397  G4double E1 = fParticleEnergyVector->Energy(iPlace);
398  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
399 
400  G4double W = 1.0/(E2 - E1);
401  G4double W1 = (E2 - scaledTkin)*W;
402  G4double W2 = (scaledTkin - E1)*W;
403 
404  xscEl *= W1;
405  xscEl += W2*xscEl2;
406 
407  xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
408 
409  E1 = fParticleEnergyVector->Energy(iPlace);
410  E2 = fParticleEnergyVector->Energy(iPlace+1);
411 
412  W = 1.0/(E2 - E1);
413  W1 = (E2 - scaledTkin)*W;
414  W2 = (scaledTkin - E1)*W;
415 
416  xscPh *= W1;
417  xscPh += W2*xscPh2;
418 
419  cross = xscEl + xscPh;
420  }
421  if( cross <= 0.0)
422  {
423  plRatio = 2.0;
424  }
425  else
426  {
427  plRatio = xscEl/cross;
428 
429  if( plRatio > 1. || plRatio < 0.) plRatio = 2.0;
430  }
431  return plRatio;
432 }
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

Here is the caller graph for this function:

void G4PAIPhotData::Initialise ( const G4MaterialCutsCouple *  couple,
G4double  cut,
G4PAIPhotModel *  model 
)

Definition at line 127 of file G4PAIPhotData.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4PAIPhotData::SampleAlongStepPhotonTransfer ( G4int  coupleIndex,
G4double  kinEnergy,
G4double  scaledTkin,
G4double  stepFactor 
) const

Definition at line 528 of file G4PAIPhotData.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4PAIPhotData::SampleAlongStepPlasmonTransfer ( G4int  coupleIndex,
G4double  kinEnergy,
G4double  scaledTkin,
G4double  stepFactor 
) const

Definition at line 620 of file G4PAIPhotData.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4PAIPhotData::SampleAlongStepTransfer ( G4int  coupleIndex,
G4double  kinEnergy,
G4double  scaledTkin,
G4double  stepFactor 
) const

Definition at line 436 of file G4PAIPhotData.cc.

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

Here is the call graph for this function:

G4double G4PAIPhotData::SamplePostStepPhotonTransfer ( G4int  coupleIndex,
G4double  scaledTkin 
) const

Definition at line 771 of file G4PAIPhotData.cc.

773 {
774  //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
775  G4double transfer = 0.0;
776  G4double rand = G4UniformRand();
777 
778  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
779 
780  // size_t iTransfer, iTr1, iTr2;
781  G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
782 
783  G4PhysicsVector* cutv = fdNdxCutPhotonTable[coupleIndex];
784 
785  // Fermi plato, try from left
786 
787  if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
788  {
789  position = (*cutv)[nPlace]*rand;
790  transfer = GetEnergyPhotonTransfer(coupleIndex, nPlace, position);
791  }
792  else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
793  {
794  position = (*cutv)[0]*rand;
795  transfer = GetEnergyPhotonTransfer(coupleIndex, 0, position);
796  }
797  else
798  {
799  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
800 
801  dNdxCut1 = (*cutv)[iPlace];
802  dNdxCut2 = (*cutv)[iPlace+1];
803 
804  E1 = fParticleEnergyVector->Energy(iPlace);
805  E2 = fParticleEnergyVector->Energy(iPlace+1);
806  W = 1.0/(E2 - E1);
807  W1 = (E2 - scaledTkin)*W;
808  W2 = (scaledTkin - E1)*W;
809 
810  //G4cout<<"iPlace= " << " dNdxCut1 = "<<dNdxCut1
811  // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;
812 
813  position = dNdxCut1*rand;
814 
815  G4double tr1 = GetEnergyPhotonTransfer(coupleIndex, iPlace, position);
816 
817  position = dNdxCut2*rand;
818  G4double tr2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);
819 
820  transfer = tr1*W1 + tr2*W2;
821  }
822  //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl;
823  if(transfer < 0.0 ) { transfer = 0.0; }
824  return transfer;
825 }
G4PhysicsVector::GetMaxEnergy
G4double GetMaxEnergy() const
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
position
#define position
Definition: xmlparse.cc:622
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
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

Here is the caller graph for this function:

G4double G4PAIPhotData::SamplePostStepPlasmonTransfer ( G4int  coupleIndex,
G4double  scaledTkin 
) const

Definition at line 829 of file G4PAIPhotData.cc.

831 {
832  //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
833  G4double transfer = 0.0;
834  G4double rand = G4UniformRand();
835 
836  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
837 
838  // size_t iTransfer, iTr1, iTr2;
839  G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
840 
841  G4PhysicsVector* cutv = fdNdxCutPlasmonTable[coupleIndex];
842 
843  // Fermi plato, try from left
844  if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
845  {
846  position = (*cutv)[nPlace]*rand;
847  transfer = GetEnergyPlasmonTransfer(coupleIndex, nPlace, position);
848  }
849  else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
850  {
851  position = (*cutv)[0]*rand;
852  transfer = GetEnergyPlasmonTransfer(coupleIndex, 0, position);
853  }
854  else
855  {
856  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
857 
858  dNdxCut1 = (*cutv)[iPlace];
859  dNdxCut2 = (*cutv)[iPlace+1];
860 
861  E1 = fParticleEnergyVector->Energy(iPlace);
862  E2 = fParticleEnergyVector->Energy(iPlace+1);
863  W = 1.0/(E2 - E1);
864  W1 = (E2 - scaledTkin)*W;
865  W2 = (scaledTkin - E1)*W;
866 
867  //G4cout<<"iPlace= " << " dNdxCut1 = "<<dNdxCut1
868  // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;
869 
870  position = dNdxCut1*rand;
871  G4double tr1 = GetEnergyPlasmonTransfer(coupleIndex, iPlace, position);
872 
873  position = dNdxCut2*rand;
874  G4double tr2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1, position);
875 
876  transfer = tr1*W1 + tr2*W2;
877  }
878  //G4cout<<"PAImodel PostStepPlasmonTransfer = "<<transfer/keV<<" keV"<<G4endl;
879 
880  if(transfer < 0.0 ) transfer = 0.0;
881 
882  return transfer;
883 }
G4PhysicsVector::GetMaxEnergy
G4double GetMaxEnergy() const
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
position
#define position
Definition: xmlparse.cc:622
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
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

Here is the caller graph for this function:

G4double G4PAIPhotData::SamplePostStepTransfer ( G4int  coupleIndex,
G4double  scaledTkin 
) const

Definition at line 715 of file G4PAIPhotData.cc.

717 {
718  //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
719  G4double transfer = 0.0;
720  G4double rand = G4UniformRand();
721 
722  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
723 
724  // size_t iTransfer, iTr1, iTr2;
725  G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
726 
727  G4PhysicsVector* cutv = fdNdxCutTable[coupleIndex];
728 
729  // Fermi plato, try from left
730  if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
731  {
732  position = (*cutv)[nPlace]*rand;
733  transfer = GetEnergyTransfer(coupleIndex, nPlace, position);
734  }
735  else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
736  {
737  position = (*cutv)[0]*rand;
738  transfer = GetEnergyTransfer(coupleIndex, 0, position);
739  }
740  else
741  {
742  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
743 
744  dNdxCut1 = (*cutv)[iPlace];
745  dNdxCut2 = (*cutv)[iPlace+1];
746 
747  E1 = fParticleEnergyVector->Energy(iPlace);
748  E2 = fParticleEnergyVector->Energy(iPlace+1);
749  W = 1.0/(E2 - E1);
750  W1 = (E2 - scaledTkin)*W;
751  W2 = (scaledTkin - E1)*W;
752 
753  //G4cout<<"iPlace= " << " dNdxCut1 = "<<dNdxCut1
754  // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;
755 
756  position = dNdxCut1*rand;
757  G4double tr1 = GetEnergyTransfer(coupleIndex, iPlace, position);
758 
759  position = dNdxCut2*rand;
760  G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
761 
762  transfer = tr1*W1 + tr2*W2;
763  }
764  //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl;
765  if(transfer < 0.0 ) { transfer = 0.0; }
766  return transfer;
767 }
G4PhysicsVector::GetMaxEnergy
G4double GetMaxEnergy() const
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
position
#define position
Definition: xmlparse.cc:622
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
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
The documentation for this class was generated from the following files: