G4PAIPhotData Class Reference

#include <G4PAIPhotData.hh>

Collaboration diagram for G4PAIPhotData:

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

Private Member Functions
G4double	GetEnergyTransfer (G4int coupleIndex, size_t iPlace, G4double position) const
G4double	GetEnergyPhotonTransfer (G4int coupleIndex, size_t iPlace, G4double position) const
G4double	GetEnergyPlasmonTransfer (G4int coupleIndex, size_t iPlace, G4double position) const
G4PAIPhotData &	operator= (const G4PAIPhotData &right)
	G4PAIPhotData (const G4PAIPhotData &)

Private Attributes
G4int	fTotBin
G4double	fLowestKineticEnergy
G4double	fHighestKineticEnergy
G4PhysicsLogVector *	fParticleEnergyVector
G4PAIxSection	fPAIxSection
G4SandiaTable	fSandia
std::vector< G4PhysicsTable * >	fPAIxscBank
std::vector< G4PhysicsTable * >	fPAIphotonBank
std::vector< G4PhysicsTable * >	fPAIplasmonBank
std::vector< G4PhysicsTable * >	fPAIdEdxBank
std::vector< G4PhysicsLogVector * >	fdEdxTable
std::vector< G4PhysicsLogVector * >	fdNdxCutTable
std::vector< G4PhysicsLogVector * >	fdNdxCutPhotonTable
std::vector< G4PhysicsLogVector * >	fdNdxCutPlasmonTable
std::vector< G4PhysicsLogVector * >	fdEdxCutTable

Detailed Description

Definition at line 66 of file G4PAIPhotData.hh.

Constructor & Destructor Documentation

G4PAIPhotData() [1/2]

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

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;
  }
}

~G4PAIPhotData()

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;  
}

G4PAIPhotData() [2/2]

G4PAIPhotData::G4PAIPhotData ( const G4PAIPhotData &  )

private

Member Function Documentation

CrossSectionPerVolume()

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:

DEDXPerVolume()

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:

GetEnergyPhotonTransfer()

G4double G4PAIPhotData::GetEnergyPhotonTransfer ( G4int  coupleIndex, size_t  iPlace, G4double  position  ) const

private

Definition at line 941 of file G4PAIPhotData.cc.

{ 
  G4PhysicsVector* v = (*(fPAIphotonBank[coupleIndex]))(iPlace); 
  if(position*v->Energy(0) >= (*v)[0])  return v->Energy(0); 

  size_t iTransferMax = v->GetVectorLength() - 1;

  size_t iTransfer;
  G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;

  for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) 
  {
    x2 = v->Energy(iTransfer);
    y2 = (*v)[iTransfer]/x2;
    if(position >= y2)  break; 
  }
  x1 = v->Energy(iTransfer-1);
  y1 = (*v)[iTransfer-1]/x1;

  //G4cout << "i= " << iTransfer << " imax= " << iTransferMax
  //     << " x1= " << x1 << " x2= " << x2 << G4endl;

  energyTransfer = x1;

  if ( x1 != x2 ) 
  {
    if ( y1 == y2  ) 
    {
      energyTransfer += (x2 - x1)*G4UniformRand();
    } 
    else 
    {
      if( x1*1.1 < x2 ) 
      {
    const G4int nbins = 5;
        G4double del = (x2 - x1)/G4int(nbins);
        x2  = x1;

        for(G4int i=1; i<=nbins; ++i) 
        {
          x2 += del;
          y2 = v->Value(x2)/x2;
          if(position >= y2) { break; }
          x1 = x2;
          y1 = y2;
    }
      }
      energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
    }
  }
  //  G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
  //     << " y1= " << y1 << " y2= " << y2 << " pos= " << position
  //     << " E(keV)= " << energyTransfer/keV << G4endl;
 
  return energyTransfer;
}

Here is the call graph for this function:

GetEnergyPlasmonTransfer()

G4double G4PAIPhotData::GetEnergyPlasmonTransfer ( G4int  coupleIndex, size_t  iPlace, G4double  position  ) const

private

Definition at line 1002 of file G4PAIPhotData.cc.

{ 
  G4PhysicsVector* v = (*(fPAIplasmonBank[coupleIndex]))(iPlace); 

  if( position*v->Energy(0) >= (*v)[0] )  return v->Energy(0); 

  size_t iTransferMax = v->GetVectorLength() - 1;

  size_t iTransfer;
  G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;

  for(iTransfer = 1; iTransfer <= iTransferMax; ++iTransfer) 
  {
    x2 = v->Energy(iTransfer);
    y2 = (*v)[iTransfer]/x2;
    if(position >= y2)  break; 
  }
  x1 = v->Energy(iTransfer-1);
  y1 = (*v)[iTransfer-1]/x1;

  //G4cout << "i= " << iTransfer << " imax= " << iTransferMax
  //     << " x1= " << x1 << " x2= " << x2 << G4endl;

  energyTransfer = x1;

  if ( x1 != x2 ) 
  {
    if ( y1 == y2  ) 
    {
      energyTransfer += (x2 - x1)*G4UniformRand();
    } 
    else 
    {
      if(x1*1.1 < x2) 
      {
    const G4int nbins = 5;
        G4double del = (x2 - x1)/G4int(nbins);
        x2  = x1;

        for( G4int i = 1; i <= nbins; ++i ) 
        {
          x2 += del;
          y2 = v->Value(x2)/x2;

          if(position >= y2)  break; 

          x1 = x2;
          y1 = y2;
    }
      }
      energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
    }
  }
  //  G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
  //     << " y1= " << y1 << " y2= " << y2 << " pos= " << position
  //     << " E(keV)= " << energyTransfer/keV << G4endl; 

  return energyTransfer;
}

Here is the call graph for this function:

GetEnergyTransfer()

G4double G4PAIPhotData::GetEnergyTransfer ( G4int  coupleIndex, size_t  iPlace, G4double  position  ) const

private

Definition at line 890 of file G4PAIPhotData.cc.

{ 
  G4PhysicsVector* v = (*(fPAIxscBank[coupleIndex]))(iPlace); 
  if(position*v->Energy(0) >= (*v)[0]) { return v->Energy(0); }

  size_t iTransferMax = v->GetVectorLength() - 1;

  size_t iTransfer;
  G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;

  for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
    x2 = v->Energy(iTransfer);
    y2 = (*v)[iTransfer]/x2;
    if(position >= y2) { break; }
  }

  x1 = v->Energy(iTransfer-1);
  y1 = (*v)[iTransfer-1]/x1;
  //G4cout << "i= " << iTransfer << " imax= " << iTransferMax
  //     << " x1= " << x1 << " x2= " << x2 << G4endl;

  energyTransfer = x1;
  if ( x1 != x2 ) {
    if ( y1 == y2  ) {
      energyTransfer += (x2 - x1)*G4UniformRand();
    } else {
      if(x1*1.1 < x2) {
    const G4int nbins = 5;
        G4double del = (x2 - x1)/G4int(nbins);
        x2  = x1;
        for(G4int i=1; i<=nbins; ++i) {
          x2 += del;
          y2 = v->Value(x2)/x2;
          if(position >= y2) { break; }
          x1 = x2;
          y1 = y2;
    }
      }
      energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
    }
  }
  //  G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
  //     << " y1= " << y1 << " y2= " << y2 << " pos= " << position
  //     << " E(keV)= " << energyTransfer/keV << G4endl; 
  return energyTransfer;
}

Here is the call graph for this function:

GetPlasmonRatio()

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:

Initialise()

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:

operator=()

G4PAIPhotData& G4PAIPhotData::operator= ( const G4PAIPhotData &  right )

private

SampleAlongStepPhotonTransfer()

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:

SampleAlongStepPlasmonTransfer()

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:

SampleAlongStepTransfer()

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:

SamplePostStepPhotonTransfer()

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:

SamplePostStepPlasmonTransfer()

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:

SamplePostStepTransfer()

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;
}

Member Data Documentation

fdEdxCutTable

std::vector<G4PhysicsLogVector*> G4PAIPhotData::fdEdxCutTable

private

Definition at line 135 of file G4PAIPhotData.hh.

fdEdxTable

std::vector<G4PhysicsLogVector*> G4PAIPhotData::fdEdxTable

private

Definition at line 129 of file G4PAIPhotData.hh.

fdNdxCutPhotonTable

std::vector<G4PhysicsLogVector*> G4PAIPhotData::fdNdxCutPhotonTable

private

Definition at line 132 of file G4PAIPhotData.hh.

fdNdxCutPlasmonTable

std::vector<G4PhysicsLogVector*> G4PAIPhotData::fdNdxCutPlasmonTable

private

Definition at line 133 of file G4PAIPhotData.hh.

fdNdxCutTable

std::vector<G4PhysicsLogVector*> G4PAIPhotData::fdNdxCutTable

private

Definition at line 131 of file G4PAIPhotData.hh.

fHighestKineticEnergy

G4double G4PAIPhotData::fHighestKineticEnergy

private

Definition at line 117 of file G4PAIPhotData.hh.

fLowestKineticEnergy

G4double G4PAIPhotData::fLowestKineticEnergy

private

Definition at line 116 of file G4PAIPhotData.hh.

fPAIdEdxBank

std::vector<G4PhysicsTable*> G4PAIPhotData::fPAIdEdxBank

private

Definition at line 128 of file G4PAIPhotData.hh.

fPAIphotonBank

std::vector<G4PhysicsTable*> G4PAIPhotData::fPAIphotonBank

private

Definition at line 125 of file G4PAIPhotData.hh.

fPAIplasmonBank

std::vector<G4PhysicsTable*> G4PAIPhotData::fPAIplasmonBank

private

Definition at line 126 of file G4PAIPhotData.hh.

fPAIxscBank

std::vector<G4PhysicsTable*> G4PAIPhotData::fPAIxscBank

private

Definition at line 124 of file G4PAIPhotData.hh.

fPAIxSection

G4PAIxSection G4PAIPhotData::fPAIxSection

private

Definition at line 121 of file G4PAIPhotData.hh.

fParticleEnergyVector

G4PhysicsLogVector* G4PAIPhotData::fParticleEnergyVector

private

Definition at line 119 of file G4PAIPhotData.hh.

fSandia

G4SandiaTable G4PAIPhotData::fSandia

private

Definition at line 122 of file G4PAIPhotData.hh.

fTotBin

G4int G4PAIPhotData::fTotBin

private

Definition at line 115 of file G4PAIPhotData.hh.

---

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

• G4PAIPhotData.hh
• G4PAIPhotData.cc