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.

{ 
  const G4int nPerDecade     = 10; 
  const G4double lowestTkin  = 50*keV;
  const G4double highestTkin = 10*TeV;

  // fPAIxSection.SetVerbose(ver);

  fLowestKineticEnergy  = std::max(tmin, lowestTkin);
  fHighestKineticEnergy = tmax;

  if(tmax < 10*fLowestKineticEnergy) 
  { 
    fHighestKineticEnergy = 10*fLowestKineticEnergy;
  } 
  else if(tmax > highestTkin) 
  {
    fHighestKineticEnergy = std::max(highestTkin, 10*fLowestKineticEnergy);
  }
  fTotBin = (G4int)(nPerDecade*
            std::log10(fHighestKineticEnergy/fLowestKineticEnergy));

  fParticleEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
                         fHighestKineticEnergy,
                         fTotBin);
  if(0 < ver) {
    G4cout << "### G4PAIPhotData: Nbins= " << fTotBin
       << " Tmin(MeV)= " << fLowestKineticEnergy/MeV
       << " Tmax(GeV)= " << fHighestKineticEnergy/GeV 
       << "  tmin(keV)= " << tmin/keV << G4endl;
  }
}

Here is the call graph for this function:

G4PAIPhotData::~G4PAIPhotData

(

)

Definition at line 94 of file G4PAIPhotData.cc.

{
  //G4cout << "G4PAIPhotData::~G4PAIPhotData() " << this << G4endl;
  size_t n = fPAIxscBank.size();
  if(0 < n) 
  {
    for(size_t i=0; i<n; ++i) 
    {
      if(fPAIxscBank[i]) 
      {
    fPAIxscBank[i]->clearAndDestroy();
    delete fPAIxscBank[i];
    fPAIxscBank[i] = 0;
      }
      if(fPAIdEdxBank[i]) 
      {
    fPAIdEdxBank[i]->clearAndDestroy();
    delete fPAIdEdxBank[i];
    fPAIdEdxBank[i]= 0;
      }
      delete fdEdxTable[i];
      delete fdNdxCutTable[i];
      fdEdxTable[i] = 0;
      fdNdxCutTable[i] = 0;
    }
  }
  delete fParticleEnergyVector;
  fParticleEnergyVector = 0;
  //G4cout << "G4PAIPhotData::~G4PAIPhotData() done for " << this << G4endl;  
}

Member Function Documentation

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

Definition at line 308 of file G4PAIPhotData.cc.

{
  G4double cross, xscEl, xscEl2, xscPh, xscPh2;

  cross=tcut+tmax;

  // iPlace is in interval from 0 to (N-1)

  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;

  G4bool one = true;

  if(      scaledTkin >= fParticleEnergyVector->Energy(nPlace))  iPlace = nPlace; 
  else if( scaledTkin > fParticleEnergyVector->Energy(0)      )   one   = false; 
  

  xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
  xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);

  xscPh = xscPh2;
  xscEl = xscEl2;

  cross  = xscPh + xscEl;
 
  if( !one ) 
  {
    xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);

    G4double E1 = fParticleEnergyVector->Energy(iPlace); 
    G4double E2 = fParticleEnergyVector->Energy(iPlace+1);

    G4double W  = 1.0/(E2 - E1);
    G4double W1 = (E2 - scaledTkin)*W;
    G4double W2 = (scaledTkin - E1)*W;

    xscEl *= W1;
    xscEl += W2*xscEl2;

    xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);

    E1 = fParticleEnergyVector->Energy(iPlace); 
    E2 = fParticleEnergyVector->Energy(iPlace+1);

    W  = 1.0/(E2 - E1);
    W1 = (E2 - scaledTkin)*W;
    W2 = (scaledTkin - E1)*W;

    xscPh *= W1;
    xscPh += W2*xscPh2;

    cross = xscEl + xscPh;
  }
  if( cross < 0.0)  cross = 0.0; 

  return cross;
}

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.

{
  // VI: iPlace is the low edge index of the bin
  // iPlace is in interval from 0 to (N-1)
  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;

  G4bool one = true;
  if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
  else if(scaledTkin > fParticleEnergyVector->Energy(0)) { 
    one = false; 
  }

  // VI: apply interpolation of the vector
  G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin);
  G4double del  = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
  if(!one) {
    G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
    G4double E1 = fParticleEnergyVector->Energy(iPlace); 
    G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
    G4double W  = 1.0/(E2 - E1);
    G4double W1 = (E2 - scaledTkin)*W;
    G4double W2 = (scaledTkin - E1)*W;
    del *= W1;
    del += W2*del2;
  }
  dEdx -= del;

  if( dEdx < 0.) { dEdx = 0.; }
  return dEdx;
}

Here is the caller graph for this function:

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

Definition at line 371 of file G4PAIPhotData.cc.

{
  G4double cross, xscEl, xscEl2, xscPh, xscPh2, plRatio;
  // iPlace is in interval from 0 to (N-1)

  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;

  G4bool one = true;

  if(      scaledTkin >= fParticleEnergyVector->Energy(nPlace))  iPlace = nPlace; 
  else if( scaledTkin > fParticleEnergyVector->Energy(0)      )   one   = false; 
  

  xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
  xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);

  xscPh = xscPh2;
  xscEl = xscEl2;

  cross  = xscPh + xscEl;
 
  if( !one ) 
  {
    xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);

    G4double E1 = fParticleEnergyVector->Energy(iPlace); 
    G4double E2 = fParticleEnergyVector->Energy(iPlace+1);

    G4double W  = 1.0/(E2 - E1);
    G4double W1 = (E2 - scaledTkin)*W;
    G4double W2 = (scaledTkin - E1)*W;

    xscEl *= W1;
    xscEl += W2*xscEl2;

    xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);

    E1 = fParticleEnergyVector->Energy(iPlace); 
    E2 = fParticleEnergyVector->Energy(iPlace+1);

    W  = 1.0/(E2 - E1);
    W1 = (E2 - scaledTkin)*W;
    W2 = (scaledTkin - E1)*W;

    xscPh *= W1;
    xscPh += W2*xscPh2;

    cross = xscEl + xscPh;
  }
  if( cross <= 0.0)  
  {
    plRatio = 2.0; 
  }
  else
  {
    plRatio = xscEl/cross;

    if( plRatio > 1. || plRatio < 0.) plRatio = 2.0;
  }
  return plRatio;
}

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.

{
  G4ProductionCutsTable* theCoupleTable=
        G4ProductionCutsTable::GetProductionCutsTable();
  size_t numOfCouples = theCoupleTable->GetTableSize();
  size_t jMatCC;

  for (jMatCC = 0; jMatCC < numOfCouples; jMatCC++ )
  {
    if( couple == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
  }
  if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;

  const vector<G4double>*  deltaCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
  const vector<G4double>*  photonCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
  G4double deltaCutInKineticEnergyNow  = (*deltaCutInKineticEnergy)[jMatCC];
  G4double photonCutInKineticEnergyNow = (*photonCutInKineticEnergy)[jMatCC];

  G4cout<<"G4PAIPhotData::Initialise: "<<"cut = "<<cut/keV<<" keV; cutEl = "
        <<deltaCutInKineticEnergyNow/keV<<" keV; cutPh = "
    <<photonCutInKineticEnergyNow/keV<<" keV"<<G4endl;

  // if( deltaCutInKineticEnergyNow != cut ) deltaCutInKineticEnergyNow = cut; // exception??

  G4PhysicsLogVector* dEdxCutVector =
    new G4PhysicsLogVector(fLowestKineticEnergy,
               fHighestKineticEnergy,
               fTotBin);

  G4PhysicsLogVector* dNdxCutVector = 
    new G4PhysicsLogVector(fLowestKineticEnergy,
               fHighestKineticEnergy,
               fTotBin);
  G4PhysicsLogVector* dNdxCutPhotonVector = 
    new G4PhysicsLogVector(fLowestKineticEnergy,
               fHighestKineticEnergy,
               fTotBin);
  G4PhysicsLogVector* dNdxCutPlasmonVector = 
    new G4PhysicsLogVector(fLowestKineticEnergy,
               fHighestKineticEnergy,
               fTotBin);

  const G4Material* mat = couple->GetMaterial();     
  fSandia.Initialize(const_cast<G4Material*>(mat));

  G4PhysicsTable* PAItransferTable = new G4PhysicsTable(fTotBin+1);
  G4PhysicsTable* PAIphotonTable = new G4PhysicsTable(fTotBin+1);
  G4PhysicsTable* PAIplasmonTable = new G4PhysicsTable(fTotBin+1);

  G4PhysicsTable* PAIdEdxTable = new G4PhysicsTable(fTotBin+1);
  G4PhysicsLogVector* dEdxMeanVector =
    new G4PhysicsLogVector(fLowestKineticEnergy,
               fHighestKineticEnergy,
               fTotBin);

  // low energy Sandia interval
  G4double Tmin = fSandia.GetSandiaMatTablePAI(0,0); 

  // energy safety
  const G4double deltaLow = 100.*eV; 

  for (G4int i = 0; i <= fTotBin; ++i) 
  {
    G4double kinEnergy = fParticleEnergyVector->Energy(i);
    G4double Tmax = model->ComputeMaxEnergy(kinEnergy);
    G4double tau = kinEnergy/proton_mass_c2;
    G4double bg2 = tau*( tau + 2. );

    if ( Tmax < Tmin + deltaLow ) Tmax = Tmin + deltaLow; 

    fPAIxSection.Initialize( mat, Tmax, bg2, &fSandia);

    //G4cout << i << ". TransferMax(keV)= "<< Tmax/keV << "  cut(keV)= " 
    //     << cut/keV << "  E(MeV)= " << kinEnergy/MeV << G4endl;

    G4int n = fPAIxSection.GetSplineSize();

    G4PhysicsFreeVector* transferVector = new G4PhysicsFreeVector(n);
    G4PhysicsFreeVector* photonVector   = new G4PhysicsFreeVector(n);
    G4PhysicsFreeVector* plasmonVector  = new G4PhysicsFreeVector(n);

    G4PhysicsFreeVector* dEdxVector     = new G4PhysicsFreeVector(n);

    for( G4int k = 0; k < n; k++ )
    {
      G4double t = fPAIxSection.GetSplineEnergy(k+1);

      transferVector->PutValue(k , t, 
                               t*fPAIxSection.GetIntegralPAIxSection(k+1));
      photonVector->PutValue(k , t, 
                               t*fPAIxSection.GetIntegralCerenkov(k+1));
      plasmonVector->PutValue(k , t, 
                               t*fPAIxSection.GetIntegralPlasmon(k+1));

      dEdxVector->PutValue(k, t, fPAIxSection.GetIntegralPAIdEdx(k+1));
    }
    // G4cout << *transferVector << G4endl;

    G4double ionloss = fPAIxSection.GetMeanEnergyLoss();//  total <dE/dx>

    if(ionloss < 0.0) ionloss = 0.0; 

    dEdxMeanVector->PutValue(i,ionloss);

    G4double dNdxCut = transferVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
    G4double dNdxCutPhoton = photonVector->Value(photonCutInKineticEnergyNow)/photonCutInKineticEnergyNow;
    G4double dNdxCutPlasmon = plasmonVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;

    G4double dEdxCut = dEdxVector->Value(cut)/cut;
    //G4cout << "i= " << i << " x= " << dNdxCut << G4endl;

    if(dNdxCut < 0.0) { dNdxCut = 0.0; }
    if(dNdxCutPhoton < 0.0) { dNdxCutPhoton = 0.0; }
    if(dNdxCutPlasmon < 0.0) { dNdxCutPlasmon = 0.0; }

    dNdxCutVector->PutValue(i, dNdxCut);
    dNdxCutPhotonVector->PutValue(i, dNdxCutPhoton);
    dNdxCutPlasmonVector->PutValue(i, dNdxCutPlasmon);

    dEdxCutVector->PutValue(i, dEdxCut);

    PAItransferTable->insertAt(i,transferVector);
    PAIphotonTable->insertAt(i,photonVector);
    PAIplasmonTable->insertAt(i,plasmonVector);
    PAIdEdxTable->insertAt(i,dEdxVector);

  } // end of Tkin loop

  fPAIxscBank.push_back(PAItransferTable);
  fPAIphotonBank.push_back(PAIphotonTable);
  fPAIplasmonBank.push_back(PAIplasmonTable);

  fPAIdEdxBank.push_back(PAIdEdxTable);
  fdEdxTable.push_back(dEdxMeanVector);

  fdNdxCutTable.push_back(dNdxCutVector);
  fdNdxCutPhotonTable.push_back(dNdxCutPhotonVector);
  fdNdxCutPlasmonTable.push_back(dNdxCutPlasmonVector);

  fdEdxCutTable.push_back(dEdxCutVector);
}

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.

{
  G4double loss = 0.0;
  G4double omega; 
  G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;

  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
 
  G4bool one = true;

  if     (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace; 
  else if(scaledTkin > fParticleEnergyVector->Energy(0))          one = false; 

  G4PhysicsLogVector* vcut = fdNdxCutPhotonTable[coupleIndex];
  G4PhysicsVector*      v1 = (*(fPAIphotonBank[coupleIndex]))(iPlace);
  G4PhysicsVector*      v2 = 0;

  dNdxCut1    = (*vcut)[iPlace];
  G4double e1 = v1->Energy(0);
  G4double e2 = e1;

  G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;

  meanNumber = meanN1;

  // G4cout<<"iPlace = "<<iPlace<< " meanN1= " << meanN1 
  //    <<"    (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " dNdxCut1= " << dNdxCut1 << G4endl;

  dNdxCut2 = dNdxCut1;
  W1 = 1.0;
  W2 = 0.0;
  if(!one) 
  {
    v2 = (*(fPAIphotonBank[coupleIndex]))(iPlace+1);
    dNdxCut2 = (*vcut)[iPlace+1];
    e2 = v2->Energy(0);

    G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;

    E1 = fParticleEnergyVector->Energy(iPlace); 
    E2 = fParticleEnergyVector->Energy(iPlace+1);
    W = 1.0/(E2 - E1);
    W1 = (E2 - scaledTkin)*W;
    W2 = (scaledTkin - E1)*W;
    meanNumber = W1*meanN1 + W2*meanN2;

    //G4cout<<"meanN= " <<  meanNumber << " meanN2= " << meanN2 
    //    << " dNdxCut2= " << dNdxCut2 << G4endl;
  }
  if( meanNumber <= 0.0) return 0.0; 

  G4int numOfCollisions = G4Poisson(meanNumber);

  //G4cout << "N= " << numOfCollisions << G4endl;

  if( 0 == numOfCollisions) return 0.0; 

  for(G4int i=0; i< numOfCollisions; ++i) 
  {
    G4double rand = G4UniformRand();
    position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
    omega = GetEnergyPhotonTransfer(coupleIndex, iPlace, position);

    //G4cout << "omega(keV)= " << omega/keV << G4endl;

    if(!one) 
    {
      position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
      G4double omega2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);
      omega = omega*W1 + omega2*W2;
    }
    //G4cout << "omega(keV)= " << omega/keV << G4endl;

    loss += omega;
    if( loss > kinEnergy) { break; }
  }
  
  // G4cout<<"PAIPhotData AlongStepLoss = "<<loss/keV<<" keV, on step = "
  //<<step/mm<<" mm"<<G4endl; 

  if     ( loss > kinEnergy) loss = kinEnergy; 
  else if( loss < 0.)        loss = 0.;
 
  return loss;
}

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.

{
  G4double loss = 0.0;
  G4double omega; 
  G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;

  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
 
  G4bool one = true;

  if     (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace; 
  else if(scaledTkin > fParticleEnergyVector->Energy(0))          one = false; 

  G4PhysicsLogVector* vcut = fdNdxCutPlasmonTable[coupleIndex];
  G4PhysicsVector*      v1 = (*(fPAIplasmonBank[coupleIndex]))(iPlace);
  G4PhysicsVector*      v2 = 0;

  dNdxCut1    = (*vcut)[iPlace];
  G4double e1 = v1->Energy(0);
  G4double e2 = e1;

  G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;

  meanNumber = meanN1;

  // G4cout<<"iPlace = "<<iPlace<< " meanN1= " << meanN1 
  //    <<"    (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " dNdxCut1= " << dNdxCut1 << G4endl;

  dNdxCut2 = dNdxCut1;
  W1 = 1.0;
  W2 = 0.0;
  if(!one) 
  {
    v2 = (*(fPAIplasmonBank[coupleIndex]))(iPlace+1);
    dNdxCut2 = (*vcut)[iPlace+1];
    e2 = v2->Energy(0);

    G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;

    E1 = fParticleEnergyVector->Energy(iPlace); 
    E2 = fParticleEnergyVector->Energy(iPlace+1);
    W = 1.0/(E2 - E1);
    W1 = (E2 - scaledTkin)*W;
    W2 = (scaledTkin - E1)*W;
    meanNumber = W1*meanN1 + W2*meanN2;

    //G4cout<<"meanN= " <<  meanNumber << " meanN2= " << meanN2 
    //    << " dNdxCut2= " << dNdxCut2 << G4endl;
  }
  if( meanNumber <= 0.0) return 0.0; 

  G4int numOfCollisions = G4Poisson(meanNumber);

  //G4cout << "N= " << numOfCollisions << G4endl;

  if( 0 == numOfCollisions) return 0.0; 

  for(G4int i=0; i< numOfCollisions; ++i) 
  {
    G4double rand = G4UniformRand();
    position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
    omega = GetEnergyPlasmonTransfer(coupleIndex, iPlace, position);

    //G4cout << "omega(keV)= " << omega/keV << G4endl;

    if(!one) 
    {
      position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
      G4double omega2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1, position);
      omega = omega*W1 + omega2*W2;
    }
    //G4cout << "omega(keV)= " << omega/keV << G4endl;

    loss += omega;
    if( loss > kinEnergy) { break; }
  }
  
  // G4cout<<"PAIPhotData AlongStepLoss = "<<loss/keV<<" keV, on step = "
  //<<step/mm<<" mm"<<G4endl; 

  if     ( loss > kinEnergy) loss = kinEnergy; 
  else if( loss < 0.)        loss = 0.;
 
  return loss;
}

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.

{
  G4double loss = 0.0;
  G4double omega; 
  G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;

  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
 
  G4bool one = true;

  if     (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace; 
  else if(scaledTkin > fParticleEnergyVector->Energy(0))          one = false; 

  G4PhysicsLogVector* vcut = fdNdxCutTable[coupleIndex];
  G4PhysicsVector*      v1 = (*(fPAIxscBank[coupleIndex]))(iPlace);
  G4PhysicsVector*      v2 = 0;

  dNdxCut1    = (*vcut)[iPlace];
  G4double e1 = v1->Energy(0);
  G4double e2 = e1;

  G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;

  meanNumber = meanN1;

  // G4cout<<"iPlace = "<<iPlace<< " meanN1= " << meanN1 
  //    <<"    (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " dNdxCut1= " << dNdxCut1 << G4endl;

  dNdxCut2 = dNdxCut1;
  W1 = 1.0;
  W2 = 0.0;
  if(!one) 
  {
    v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
    dNdxCut2 = (*vcut)[iPlace+1];
    e2 = v2->Energy(0);

    G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;

    E1 = fParticleEnergyVector->Energy(iPlace); 
    E2 = fParticleEnergyVector->Energy(iPlace+1);
    W = 1.0/(E2 - E1);
    W1 = (E2 - scaledTkin)*W;
    W2 = (scaledTkin - E1)*W;
    meanNumber = W1*meanN1 + W2*meanN2;

    //G4cout<<"meanN= " <<  meanNumber << " meanN2= " << meanN2 
    //    << " dNdxCut2= " << dNdxCut2 << G4endl;
  }
  if( meanNumber <= 0.0) return 0.0; 

  G4int numOfCollisions = G4Poisson(meanNumber);

  //G4cout << "N= " << numOfCollisions << G4endl;

  if( 0 == numOfCollisions) return 0.0; 

  for(G4int i=0; i< numOfCollisions; ++i) 
  {
    G4double rand = G4UniformRand();
    position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
    omega = GetEnergyTransfer(coupleIndex, iPlace, position);

    //G4cout << "omega(keV)= " << omega/keV << G4endl;

    if(!one) 
    {
      position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
      G4double omega2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
      omega = omega*W1 + omega2*W2;
    }
    //G4cout << "omega(keV)= " << omega/keV << G4endl;

    loss += omega;
    if( loss > kinEnergy) { break; }
  }
  
  // G4cout<<"PAIPhotData AlongStepLoss = "<<loss/keV<<" keV, on step = "
  //<<step/mm<<" mm"<<G4endl; 

  if     ( loss > kinEnergy) loss = kinEnergy; 
  else if( loss < 0.)        loss = 0.;
 
  return loss;
}

Here is the call graph for this function:

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

Definition at line 771 of file G4PAIPhotData.cc.

{  
  //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
  G4double transfer = 0.0;
  G4double rand = G4UniformRand();

  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;

  //  size_t iTransfer, iTr1, iTr2;
  G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;

  G4PhysicsVector* cutv = fdNdxCutPhotonTable[coupleIndex];

  // Fermi plato, try from left

  if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy()) 
  {
    position = (*cutv)[nPlace]*rand;
    transfer = GetEnergyPhotonTransfer(coupleIndex, nPlace, position);
  }
  else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
  {
    position = (*cutv)[0]*rand;
    transfer = GetEnergyPhotonTransfer(coupleIndex, 0, position);
  }
  else 
  {  
    size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);

    dNdxCut1 = (*cutv)[iPlace];  
    dNdxCut2 = (*cutv)[iPlace+1];  

    E1 = fParticleEnergyVector->Energy(iPlace); 
    E2 = fParticleEnergyVector->Energy(iPlace+1);
    W  = 1.0/(E2 - E1);
    W1 = (E2 - scaledTkin)*W;
    W2 = (scaledTkin - E1)*W;

    //G4cout<<"iPlace= " << "  dNdxCut1 = "<<dNdxCut1 
    //    <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;

    position = dNdxCut1*rand;

    G4double tr1 = GetEnergyPhotonTransfer(coupleIndex, iPlace, position);

    position = dNdxCut2*rand;
    G4double tr2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);

    transfer = tr1*W1 + tr2*W2;
  }
  //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl; 
  if(transfer < 0.0 ) { transfer = 0.0; }
  return transfer;
}

Here is the caller graph for this function:

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

Definition at line 829 of file G4PAIPhotData.cc.

{  
  //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
  G4double transfer = 0.0;
  G4double rand = G4UniformRand();

  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;

  //  size_t iTransfer, iTr1, iTr2;
  G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;

  G4PhysicsVector* cutv = fdNdxCutPlasmonTable[coupleIndex];

  // Fermi plato, try from left
  if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy()) 
  {
    position = (*cutv)[nPlace]*rand;
    transfer = GetEnergyPlasmonTransfer(coupleIndex, nPlace, position);
  }
  else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
  {
    position = (*cutv)[0]*rand;
    transfer = GetEnergyPlasmonTransfer(coupleIndex, 0, position);
  }
  else 
  {  
    size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);

    dNdxCut1 = (*cutv)[iPlace];  
    dNdxCut2 = (*cutv)[iPlace+1];  

    E1 = fParticleEnergyVector->Energy(iPlace); 
    E2 = fParticleEnergyVector->Energy(iPlace+1);
    W  = 1.0/(E2 - E1);
    W1 = (E2 - scaledTkin)*W;
    W2 = (scaledTkin - E1)*W;

    //G4cout<<"iPlace= " << "  dNdxCut1 = "<<dNdxCut1 
    //    <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;

    position = dNdxCut1*rand;
    G4double tr1 = GetEnergyPlasmonTransfer(coupleIndex, iPlace, position);

    position = dNdxCut2*rand;
    G4double tr2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1, position);

    transfer = tr1*W1 + tr2*W2;
  }
  //G4cout<<"PAImodel PostStepPlasmonTransfer = "<<transfer/keV<<" keV"<<G4endl; 

  if(transfer < 0.0 )  transfer = 0.0;
 
  return transfer;
}

Here is the caller graph for this function:

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

Definition at line 715 of file G4PAIPhotData.cc.

{  
  //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
  G4double transfer = 0.0;
  G4double rand = G4UniformRand();

  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;

  //  size_t iTransfer, iTr1, iTr2;
  G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;

  G4PhysicsVector* cutv = fdNdxCutTable[coupleIndex];

  // Fermi plato, try from left
  if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy()) 
  {
    position = (*cutv)[nPlace]*rand;
    transfer = GetEnergyTransfer(coupleIndex, nPlace, position);
  }
  else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
  {
    position = (*cutv)[0]*rand;
    transfer = GetEnergyTransfer(coupleIndex, 0, position);
  }
  else 
  {  
    size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);

    dNdxCut1 = (*cutv)[iPlace];  
    dNdxCut2 = (*cutv)[iPlace+1];  

    E1 = fParticleEnergyVector->Energy(iPlace); 
    E2 = fParticleEnergyVector->Energy(iPlace+1);
    W  = 1.0/(E2 - E1);
    W1 = (E2 - scaledTkin)*W;
    W2 = (scaledTkin - E1)*W;

    //G4cout<<"iPlace= " << "  dNdxCut1 = "<<dNdxCut1 
    //    <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;

    position = dNdxCut1*rand;
    G4double tr1 = GetEnergyTransfer(coupleIndex, iPlace, position);

    position = dNdxCut2*rand;
    G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);

    transfer = tr1*W1 + tr2*W2;
  }
  //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl; 
  if(transfer < 0.0 ) { transfer = 0.0; }
  return transfer;
}

---

The documentation for this class was generated from the following files:

• source/geant4.10.03.p03/source/processes/electromagnetic/standard/include/G4PAIPhotData.hh
• source/geant4.10.03.p03/source/processes/electromagnetic/standard/src/G4PAIPhotData.cc