G4DNAIRT_geometries Class Reference

#include <G4DNAIRT_geometries.hh>

Inheritance diagram for G4DNAIRT_geometries:

Public Member Functions
G4int	FindBin (G4int, G4double, G4double, G4double)
std::vector< std::unique_ptr< G4ITReactionChange > >	FindReaction (G4ITReactionSet *, const G4double, const G4double, const G4bool) override
	G4DNAIRT_geometries ()
	G4DNAIRT_geometries (const G4DNAIRT_geometries &other)=delete
	G4DNAIRT_geometries (G4VDNAReactionModel *)
G4double	GetIndependentReactionTime (const G4MolecularConfiguration *, const G4MolecularConfiguration *, G4double)
void	Initialize () override
void	IRTSampling ()
virtual G4bool	IsApplicable (const G4ITType &, const G4ITType &) const
std::unique_ptr< G4ITReactionChange >	MakeReaction (const G4Track &, const G4Track &) override
G4DNAIRT_geometries &	operator= (const G4DNAIRT_geometries &other)=delete
G4double	SamplePDC (G4double, G4double)
void	Sampling (G4Track *)
void	SetReactionModel (G4VDNAReactionModel *)
virtual void	SetReactionTable (const G4ITReactionTable *)
void	SpaceBinning ()
G4bool	TestReactibility (const G4Track &, const G4Track &, G4double, G4bool) override
	~G4DNAIRT_geometries () override

Protected Attributes
const G4DNAMolecularReactionTable *&	fMolReactionTable
G4String	fName
G4VDNAReactionModel *	fpReactionModel
const G4ITReactionTable *	fpReactionTable = nullptr

Private Attributes
G4ErrorFunction *	erfc
G4VDNAMolecularGeometry *	fGeometry
G4int	fNx
G4int	fNy
G4int	fNz
G4double	fRCutOff
G4ITReactionSet *	fReactionSet
G4ITTrackHolder *	fTrackHolder
G4double	fXMax
G4double	fXMin
G4double	fYMax
G4double	fYMin
G4double	fZMax
G4double	fZMin
std::vector< std::pair< G4ThreeVector, G4Track * > >	positionMap
std::map< G4int, std::map< G4int, std::map< G4int, std::vector< G4Track * > > > >	spaceBinned
G4double	timeMax
G4double	timeMin
G4int	xendIndex
G4int	xiniIndex
G4int	yendIndex
G4int	yiniIndex
G4int	zendIndex
G4int	ziniIndex

Detailed Description

Definition at line 65 of file G4DNAIRT_geometries.hh.

Constructor & Destructor Documentation

G4DNAIRT_geometries() [1/3]

G4DNAIRT_geometries::G4DNAIRT_geometries

(

)

Definition at line 54 of file G4DNAIRT_geometries.cc.

                                          :
G4VITReactionProcess(),
fMolReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable)),
fpReactionModel(nullptr),
fTrackHolder(G4ITTrackHolder::Instance()),
fReactionSet(nullptr),
fGeometry(nullptr)
{
  timeMin = G4Scheduler::Instance()->GetStartTime();
  timeMax = G4Scheduler::Instance()->GetEndTime();
 
  fXMin = 1e9*nm;
  fYMin = 1e9*nm;
  fZMin = 1e9*nm;
 
  fXMax = 0e0*nm;
  fYMax = 0e0*nm;
  fZMax = 0e0*nm;
 
  fNx = 0;
  fNy = 0;
  fNz = 0;
 
  xiniIndex = 0, yiniIndex = 0, ziniIndex = 0;
  xendIndex = 0, yendIndex = 0, zendIndex = 0;
 
  fRCutOff =
    1.45 * nm + 2 * std::sqrt(8*9.46e9*nm*nm/s * timeMax); // 95% confidence level
 
  erfc = new G4ErrorFunction();
}

References erfc, fNx, fNy, fNz, fRCutOff, fXMax, fXMin, fYMax, fYMin, fZMax, fZMin, G4Scheduler::GetEndTime(), G4Scheduler::GetStartTime(), G4Scheduler::Instance(), nm, s, timeMax, timeMin, xendIndex, xiniIndex, yendIndex, yiniIndex, zendIndex, and ziniIndex.

G4DNAIRT_geometries() [2/3]

G4DNAIRT_geometries::G4DNAIRT_geometries ( G4VDNAReactionModel *  pReactionModel )

explicit

Definition at line 87 of file G4DNAIRT_geometries.cc.

  : G4DNAIRT_geometries()
{
  fpReactionModel = pReactionModel;
}

References fpReactionModel.

~G4DNAIRT_geometries()

G4DNAIRT_geometries::~G4DNAIRT_geometries ( )

override

Definition at line 93 of file G4DNAIRT_geometries.cc.

{
  delete erfc;
}

References erfc.

G4DNAIRT_geometries() [3/3]

G4DNAIRT_geometries::G4DNAIRT_geometries ( const G4DNAIRT_geometries &  other )

delete

Member Function Documentation

FindBin()

G4int G4DNAIRT_geometries::FindBin	(	G4int	n,
		G4double	xmin,
		G4double	xmax,
		G4double	value
	)

(xmax < value) ) //value >= xmax )

Definition at line 430 of file G4DNAIRT_geometries.cc.

                                                                                        {
 
  G4int bin = -1;
  if ( value <= xmin )
    bin = 0; //1;
  else if ( value >= xmax)  
    bin = n-1; //n;
  else
    bin = G4int( n * ( value - xmin )/( xmax - xmin ) ); //bin = 1 + G4int( n * ( value - xmin )/( xmax - xmin ) );
 
  if ( bin < 0 ) bin = 0;
  if ( bin >= n ) bin = n-1;
 
  return bin;
}

References CLHEP::detail::n.

Referenced by IRTSampling(), MakeReaction(), and Sampling().

FindReaction()

std::vector< std::unique_ptr< G4ITReactionChange > > G4DNAIRT_geometries::FindReaction ( G4ITReactionSet *  pReactionSet, const  G4double, const G4double  fGlobalTime, const  G4bool  )

overridevirtual

Implements G4VITReactionProcess.

Definition at line 590 of file G4DNAIRT_geometries.cc.

{
  std::vector<std::unique_ptr<G4ITReactionChange>> fReactionInfo;
  fReactionInfo.clear();
 
  if (pReactionSet == nullptr)
  {
    return fReactionInfo;
  }
 
  auto fReactionsetInTime = pReactionSet->GetReactionsPerTime();
  assert(fReactionsetInTime.begin() != fReactionsetInTime.end());
 
  auto it_begin = fReactionsetInTime.begin();
  while(it_begin != fReactionsetInTime.end())
  {
    G4double irt = it_begin->get()->GetTime();
 
    if(fGlobalTime < irt) break;
 
    pReactionSet->SelectThisReaction(*it_begin);
 
    G4Track* pTrackA = it_begin->get()->GetReactants().first;
    G4Track* pTrackB = it_begin->get()->GetReactants().second;
    auto pReactionChange = MakeReaction(*pTrackA, *pTrackB);
 
    if(pReactionChange){
      fReactionInfo.push_back(std::move(pReactionChange));
    }
 
    fReactionsetInTime = pReactionSet->GetReactionsPerTime();
    it_begin = fReactionsetInTime.begin();
  }
 
  return fReactionInfo;
}

References G4ITReactionSet::GetReactionsPerTime(), MakeReaction(), and G4ITReactionSet::SelectThisReaction().

GetIndependentReactionTime()

G4double G4DNAIRT_geometries::GetIndependentReactionTime	(	const G4MolecularConfiguration *	molA,
		const G4MolecularConfiguration *	molB,
		G4double	distance
	)

Definition at line 371 of file G4DNAIRT_geometries.cc.

                                                                                                                                                      {
  const auto pMoleculeA = molA;
  const auto pMoleculeB = molB;
  auto fReactionData = fMolReactionTable->GetReactionData(pMoleculeA, pMoleculeB);
  G4int reactionType = fReactionData->GetReactionType();
  G4double r0 = distance;
  if(r0 == 0) r0 += 1e-3*nm;
  G4double irt = -1 * ps;
  G4double D = molA->GetDiffusionCoefficient() +
               molB->GetDiffusionCoefficient();
  if(D == 0) D += 1e-20*(m2/s);
  G4double rc = fReactionData->GetOnsagerRadius();
 
  if ( reactionType == 0){
    G4double sigma = fReactionData->GetEffectiveReactionRadius();
 
    if(sigma > r0) return 0; // contact reaction
    if( rc != 0) r0 = -rc / (1-std::exp(rc/r0));
 
    G4double Winf = sigma/r0;
    G4double W = G4UniformRand();
 
    if ( W > 0 && W < Winf ) irt = (0.25/D) * std::pow( (r0-sigma)/erfc->erfcInv(r0*W/sigma), 2 );
 
    return irt;
  }
  else if ( reactionType == 1 ){
    G4double sigma = fReactionData->GetReactionRadius();
    G4double kact = fReactionData->GetActivationRateConstant();
    G4double kdif = fReactionData->GetDiffusionRateConstant();
    G4double kobs = fReactionData->GetObservedReactionRateConstant();
 
    G4double a, b, Winf;
 
    if ( rc == 0 ) {
      a = 1/sigma * kact / kobs;
      b = (r0 - sigma) / 2;
    } else {
      G4double v = kact/Avogadro/(4*CLHEP::pi*pow(sigma,2) * exp(-rc / sigma));
      G4double alpha = v+rc*D/(pow(sigma,2)*(1-exp(-rc/sigma)));
      a = 4*pow(sigma,2)*alpha/(D*pow(rc,2))*pow(sinh(rc/(2*sigma)),2);
      b = rc/4*(cosh(rc/(2*r0))/sinh(rc/(2*r0))-cosh(rc/(2*sigma))/sinh(rc/(2*sigma)));
      r0 = -rc/(1-std::exp(rc/r0));
      sigma = fReactionData->GetEffectiveReactionRadius();
    }
 
    if(sigma > r0){
      if(fReactionData->GetProbability() > G4UniformRand()) return 0;
      else return irt;
    }
    Winf = sigma / r0 * kobs / kdif;
 
    if(Winf > G4UniformRand()) irt = SamplePDC(a,b)/D;
    return irt;
  }
 
  return -1 * ps;
}

References alpha, source.hepunit::Avogadro, D(), erfc, G4ErrorFunction::erfcInv(), fMolReactionTable, G4UniformRand, G4MolecularConfiguration::GetDiffusionCoefficient(), G4DNAMolecularReactionTable::GetReactionData(), G4DNAMolecularReactionData::GetReactionType(), m2, nm, CLHEP::pi, ps, s, and SamplePDC().

Referenced by Sampling().

Initialize()

void G4DNAIRT_geometries::Initialize ( )

overridevirtual

First initialization (done once for all at the begin of the run) eg. check if the reaction table is given ...

Reimplemented from G4VITReactionProcess.

Definition at line 98 of file G4DNAIRT_geometries.cc.

                                    {
 
  timeMin = G4Scheduler::Instance()->GetGlobalTime();
  timeMax = std::min(timeMin + G4Scheduler::Instance()->GetLimitingTimeStep(),
                     G4Scheduler::Instance()->GetEndTime());
  if(timeMin == 0) return;
 
  fTrackHolder = G4ITTrackHolder::Instance();
  if(fTrackHolder->GetMainList()->size() == 0) return;
 
  fReactionSet = G4ITReactionSet::Instance();
  fReactionSet->CleanAllReaction();
  fReactionSet->SortByTime();
 
  spaceBinned.clear();
  positionMap.clear();
 
  fRCutOff =
    1.45 * nm + 2 * std::sqrt(8*9.46e9*nm*nm/s * (timeMax - timeMin));
 
  xiniIndex = 0;
  yiniIndex = 0;
  ziniIndex = 0;
  xendIndex = 0;
  yendIndex = 0;
  zendIndex = 0;
 
  fXMin = 1e9*nm;
  fYMin = 1e9*nm;
  fZMin = 1e9*nm;
 
  fXMax = 0e0*nm;
  fYMax = 0e0*nm;
  fZMax = 0e0*nm;
 
  fNx = 0;
  fNy = 0;
  fNz = 0;
 
  fGeometry = fMolReactionTable->GetGeometry();
 
  SpaceBinning();       // 1. binning the space
  IRTSampling();        // 2. Sampling of the IRT
 
}

References G4ITReactionSet::CleanAllReaction(), fGeometry, fMolReactionTable, fNx, fNy, fNz, fRCutOff, fReactionSet, fTrackHolder, fXMax, fXMin, fYMax, fYMin, fZMax, fZMin, G4DNAMolecularReactionTable::GetGeometry(), G4Scheduler::GetGlobalTime(), G4ITTrackHolder::GetMainList(), G4ITReactionSet::Instance(), G4ITTrackHolder::Instance(), G4Scheduler::Instance(), IRTSampling(), G4INCL::Math::min(), nm, positionMap, s, G4FastList< OBJECT >::size(), G4ITReactionSet::SortByTime(), spaceBinned, SpaceBinning(), timeMax, timeMin, xendIndex, xiniIndex, yendIndex, yiniIndex, zendIndex, and ziniIndex.

IRTSampling()

void G4DNAIRT_geometries::IRTSampling

(

)

Definition at line 180 of file G4DNAIRT_geometries.cc.

                                     {
 
  auto it_begin = fTrackHolder->GetMainList()->begin();
  while(it_begin != fTrackHolder->GetMainList()->end()){
    G4int I = FindBin(fNx, fXMin, fXMax, it_begin->GetPosition().x());
    G4int J = FindBin(fNy, fYMin, fYMax, it_begin->GetPosition().y());
    G4int K = FindBin(fNz, fZMin, fZMax, it_begin->GetPosition().z());
 
    spaceBinned[I][J][K].push_back(*it_begin);
 
    Sampling(*it_begin);
    ++it_begin;
  }
}

References G4FastList< OBJECT >::begin(), G4FastList< OBJECT >::end(), FindBin(), fNx, fNy, fNz, fTrackHolder, fXMax, fXMin, fYMax, fYMin, fZMax, fZMin, G4ITTrackHolder::GetMainList(), Sampling(), and spaceBinned.

Referenced by Initialize().

IsApplicable()

G4bool G4VITReactionProcess::IsApplicable ( const G4ITType &  , const G4ITType &    ) const

virtualinherited

Definition at line 43 of file G4VITReactionProcess.cc.

{
    return true;
}

MakeReaction()

std::unique_ptr< G4ITReactionChange > G4DNAIRT_geometries::MakeReaction ( const G4Track &  trackA, const G4Track &  trackB  )

overridevirtual

Implements G4VITReactionProcess.

Definition at line 479 of file G4DNAIRT_geometries.cc.

{
 
  std::unique_ptr<G4ITReactionChange> pChanges(new G4ITReactionChange());
  pChanges->Initialize(trackA, trackB);
 
  const auto pMoleculeA = GetMolecule(trackA)->GetMolecularConfiguration();
  const auto pMoleculeB = GetMolecule(trackB)->GetMolecularConfiguration();
  const auto pReactionData = fMolReactionTable->GetReactionData(pMoleculeA, pMoleculeB);
 
  G4double globalTime = G4Scheduler::Instance()->GetGlobalTime();
  G4double effectiveReactionRadius = pReactionData->GetEffectiveReactionRadius();
 
  const G4double D1 = pMoleculeA->GetDiffusionCoefficient();
  const G4double D2 = pMoleculeB->GetDiffusionCoefficient();
 
  G4ThreeVector r1 = trackA.GetPosition();
  G4ThreeVector r2 = trackB.GetPosition();
 
  if(r1 == r2) r2 += G4ThreeVector(0,0,1e-3*nm);
 
  G4ThreeVector S1 = r1 - r2;
 
  G4double r0 = S1.mag();
 
  S1.setMag(effectiveReactionRadius);
 
  G4double dt = globalTime - trackA.GetGlobalTime();
 
  if(dt != 0 && (D1 + D2) != 0 && r0 != 0){
    G4double s12 = 2.0 * D1 * dt;
    G4double s22 = 2.0 * D2 * dt;
    if(s12 == 0) r2 = r1;
    else if(s22 == 0) r1 = r2;
    else{
      G4double alpha = effectiveReactionRadius * r0 / (2*(D1 + D2)*dt);
      G4ThreeVector S2 = (r1 + (s12 / s22)*r2) + G4ThreeVector(G4RandGauss::shoot(0, s12 + s22 * s22 / s12),
                                                               G4RandGauss::shoot(0, s12 + s22 * s22 / s12),
                                                               G4RandGauss::shoot(0, s12 + s22 * s22 / s12));
 
      S1.setPhi(rad * G4UniformRand() * 2.0 * CLHEP::pi);
      S1.setTheta(rad * std::acos(1.0 + 1./alpha * std::log(1.0 - G4UniformRand() * (1 - std::exp(-2.0 * alpha)))));
 
      r1 = (D1 * S1 + D2 * S2) / (D1 + D2);
      r2 = D2 * (S2 - S1) / (D1 + D2);
    }
  }
 
  auto pTrackA = const_cast<G4Track*>(pChanges->GetTrackA());
  auto pTrackB = const_cast<G4Track*>(pChanges->GetTrackB());
 
  pTrackA->SetPosition(r1);
  pTrackB->SetPosition(r2);
 
  pTrackA->SetGlobalTime(globalTime);
  pTrackB->SetGlobalTime(globalTime);
 
  pTrackA->SetTrackStatus(fStopButAlive);
  pTrackB->SetTrackStatus(fStopButAlive);
 
  const G4int nbProducts = pReactionData->GetNbProducts();
 
  if(nbProducts){
 
    const G4double sqrD1 = D1 == 0. ? 0. : std::sqrt(D1);
    const G4double sqrD2 = D2 == 0. ? 0. : std::sqrt(D2);
    const G4double inv_numerator = 1./(sqrD1 + sqrD2);
    const G4ThreeVector reactionSite = sqrD2 * inv_numerator * trackA.GetPosition()
                                     + sqrD1 * inv_numerator * trackB.GetPosition();
 
    std::vector<G4ThreeVector> position;
 
    if(nbProducts == 1){
      position.push_back(reactionSite);
    }else if(nbProducts == 2){
      position.push_back(trackA.GetPosition());
      position.push_back(trackB.GetPosition());
    }else if (nbProducts == 3){
      position.push_back(reactionSite);
      position.push_back(trackA.GetPosition());
      position.push_back(trackB.GetPosition());
    }
 
    for(G4int u = 0; u < nbProducts; ++u){
 
      auto product = new G4Molecule(pReactionData->GetProduct(u));
      auto productTrack = product->BuildTrack(globalTime,
                                              position[u]);
 
      productTrack->SetTrackStatus(fAlive);
      fTrackHolder->Push(productTrack);
 
      pChanges->AddSecondary(productTrack);
 
      G4int I = FindBin(fNx, fXMin, fXMax, position[u].x());
      G4int J = FindBin(fNy, fYMin, fYMax, position[u].y());
      G4int K = FindBin(fNz, fZMin, fZMax, position[u].z());
 
      spaceBinned[I][J][K].push_back(productTrack);
 
      Sampling(productTrack);
    }
  }
 
  fTrackHolder->MergeSecondariesWithMainList();
  pChanges->KillParents(true);
  return pChanges;
}

References alpha, fAlive, FindBin(), fMolReactionTable, fNx, fNy, fNz, fStopButAlive, fTrackHolder, fXMax, fXMin, fYMax, fYMin, fZMax, fZMin, G4UniformRand, G4Scheduler::GetGlobalTime(), G4Track::GetGlobalTime(), G4Molecule::GetMolecularConfiguration(), GetMolecule(), G4Track::GetPosition(), G4DNAMolecularReactionTable::GetReactionData(), G4Scheduler::Instance(), CLHEP::Hep3Vector::mag(), G4ITTrackHolder::MergeSecondariesWithMainList(), nm, CLHEP::pi, position, G4ITTrackHolder::Push(), rad, Sampling(), CLHEP::Hep3Vector::setMag(), CLHEP::Hep3Vector::setPhi(), G4Track::SetPosition(), CLHEP::Hep3Vector::setTheta(), G4INCL::DeJongSpin::shoot(), and spaceBinned.

Referenced by FindReaction().

operator=()

G4DNAIRT_geometries & G4DNAIRT_geometries::operator= ( const G4DNAIRT_geometries &  other )

delete

SamplePDC()

G4double G4DNAIRT_geometries::SamplePDC	(	G4double	a,
		G4double	b
	)

Definition at line 446 of file G4DNAIRT_geometries.cc.

                                                              {
 
  G4double p = 2.0 * std::sqrt(2.0*b/a);
  G4double q = 2.0 / std::sqrt(2.0*b/a);
  G4double M = max(1.0/(a*a),3.0*b/a);
 
  G4double X, U, lambdax;
 
  G4int ntrials = 0;
  while(1) {
 
    // Generate X
    U = G4UniformRand();
    if ( U < p/(p + q * M) ) X = pow(U * (p + q * M) / 2, 2);
    else X = pow(2/((1-U)*(p+q*M)/M),2);
 
    U = G4UniformRand();
 
    lambdax = std::exp(-b*b/X) * ( 1.0 - a * std::sqrt(CLHEP::pi * X) * erfc->erfcx(b/std::sqrt(X) + a*std::sqrt(X)));
 
    if ((X <= 2.0*b/a && U <= lambdax) ||
        (X >= 2.0*b/a && U*M/X <= lambdax)) break;
 
    ntrials++;
 
    if ( ntrials > 10000 ){
      G4cout<<"Totally rejected"<<'\n';
      return -1.0;
    }
  }
  return X;
}

References erfc, G4ErrorFunction::erfcx(), G4cout, G4UniformRand, M, G4INCL::Math::max(), and CLHEP::pi.

Referenced by GetIndependentReactionTime().

Sampling()

void G4DNAIRT_geometries::Sampling

(

G4Track *

track

)

Definition at line 195 of file G4DNAIRT_geometries.cc.

                                                {
  G4Molecule* molA = G4Molecule::GetMolecule(track);
  const G4MolecularConfiguration* molConfA = molA->GetMolecularConfiguration();
  if(molConfA->GetDiffusionCoefficient() == 0) return;
 
  const vector<const G4MolecularConfiguration*>* reactivesVector =
    fMolReactionTable->CanReactWith(molConfA);
 
  if(reactivesVector == nullptr) return;
 
  G4double globalTime = G4Scheduler::Instance()->GetGlobalTime();
  G4double minTime = timeMax;
 
  xiniIndex = FindBin(fNx, fXMin, fXMax, track->GetPosition().x()-fRCutOff);
  xendIndex = FindBin(fNx, fXMin, fXMax, track->GetPosition().x()+fRCutOff);
  yiniIndex = FindBin(fNy, fYMin, fYMax, track->GetPosition().y()-fRCutOff);
  yendIndex = FindBin(fNy, fYMin, fYMax, track->GetPosition().y()+fRCutOff);
  ziniIndex = FindBin(fNz, fZMin, fZMax, track->GetPosition().z()-fRCutOff);
  zendIndex = FindBin(fNz, fZMin, fZMax, track->GetPosition().z()+fRCutOff);
 
  for ( G4int ii = xiniIndex; ii <= xendIndex; ++ii ) {
    for ( G4int jj = yiniIndex; jj <= yendIndex; ++jj ) {
      for ( G4int kk = ziniIndex; kk <= zendIndex; ++kk ) {
 
        std::vector<G4Track*> spaceBin = spaceBinned[ii][jj][kk];
        for ( G4int n = 0; n < (G4int)spaceBinned[ii][jj][kk].size(); ++n ) {
          if(!spaceBin[n] || track == spaceBin[n]) continue;
          if(spaceBin[n]->GetTrackStatus() == fStopButAlive) continue;
 
          G4Molecule* molB = G4Molecule::GetMolecule(spaceBin[n]);
          if(!molB) continue;
 
          const G4MolecularConfiguration* molConfB = molB->GetMolecularConfiguration();
          if(molConfB->GetDiffusionCoefficient() == 0) continue;
 
          auto it = std::find(reactivesVector->begin(), reactivesVector->end(), molConfB);
          if(it == reactivesVector->end()) continue;
 
          G4ThreeVector orgPosB = spaceBin[n]->GetPosition();
          G4double dt = track->GetGlobalTime() - spaceBin[n]->GetGlobalTime();
          G4ThreeVector newPosB = orgPosB;
 
          if(dt > 0){
            G4double sigma, x, y, z;
            G4double diffusionCoefficient = G4Molecule::GetMolecule(spaceBin[n])->GetDiffusionCoefficient();
 
            sigma = std::sqrt(2.0 * diffusionCoefficient * dt);
 
            x = G4RandGauss::shoot(0., 1.0)*sigma;
            y = G4RandGauss::shoot(0., 1.0)*sigma;
            z = G4RandGauss::shoot(0., 1.0)*sigma;
 
            newPosB = orgPosB + G4ThreeVector(x,y,z);
          }else if(dt < 0) continue;
 
          G4double r0 = (newPosB - track->GetPosition()).mag();
          G4double irt = GetIndependentReactionTime(molConfA,
                                                    molConfB,
                                                    r0);
          if(irt>=0 && irt<timeMax - globalTime)
          {
            irt += globalTime;
            if(irt < minTime) minTime = irt;
#ifdef DEBUG
            G4cout<<irt<<'\t'<<molConfA->GetName()<<" "<<track->GetTrackID()<<'\t'<<molConfB->GetName()<<" "<<spaceBin[n]->GetTrackID()<<'\n';
#endif
            fReactionSet->AddReaction(irt,track,spaceBin[n]);
          }
        }
        spaceBin.clear();
      }
    }
  }
 
// Scavenging & first order reactions
 
  auto fReactionDatas = fMolReactionTable->GetReactionData(molConfA);
  G4int index = -1;
 
  for(size_t u=0; u<fReactionDatas->size();++u){
    auto molB = (*fReactionDatas)[u]->GetReactant2();
    if(molB == G4MoleculeTable::Instance()->GetConfiguration("H2O(B)") ||
       molB == G4MoleculeTable::Instance()->GetConfiguration("H3Op(B)") ||
       molB == G4MoleculeTable::Instance()->GetConfiguration("OHm(B)")){
      G4double kObs = (*fReactionDatas)[u]->GetObservedReactionRateConstant();
      if(kObs == 0) continue;
      G4double time = -(std::log(1.0 - G4UniformRand())/kObs) + globalTime;
      if( time < minTime && time >= globalTime && time < timeMax){
        minTime = time;
        index = (G4int)u;
      }
    }
  }
 
  if(index != -1){
#ifdef DEBUG
    G4cout<<"scavenged: "<<minTime<<'\t'<<molConfA->GetName()<<it_begin->GetTrackID()<<'\n';
#endif
    G4Molecule* fakeMol = new G4Molecule((*fReactionDatas)[index]->GetReactant2());
    G4Track* fakeTrack = fakeMol->BuildTrack(globalTime,track->GetPosition());
    fTrackHolder->Push(fakeTrack);
    fReactionSet->AddReaction(minTime, track, fakeTrack);
  }
 
 
  // DNA reactions
  if(fGeometry == nullptr) return;
 
  const G4VTouchable* touchable = track->GetTouchable();
  if(touchable == nullptr) return;
 
  const G4LogicalVolume* logicalVolume = touchable->GetVolume()->GetLogicalVolume();
 
  const G4ThreeVector& globalPos = track->GetPosition();
  const G4ThreeVector& localPos = touchable->GetHistory()->GetTopTransform().TransformPoint(globalPos);
 
  G4double D = GetMolecule(track)->GetDiffusionCoefficient();
  G4double time_step = abs(timeMax - timeMin);
  G4double search_range = 2*sqrt(2*D*time_step);
 
  std::vector<G4VPhysicalVolume*> result_pv;
  result_pv.clear();
  fGeometry->FindNearbyMolecules(logicalVolume,
                                 localPos,
                                 result_pv,
                                 search_range);
 
  if(result_pv.empty()) return;
 
  for(auto physicalVolume : result_pv){
    const G4Material* material = physicalVolume->GetLogicalVolume()->GetMaterial();
    G4MolecularConfiguration* dna_molConf =
      G4DNAMolecularMaterial::Instance()->GetMolecularConfiguration(material);
 
    auto it = std::find(reactivesVector->begin(), reactivesVector->end(), dna_molConf);
    if(it == reactivesVector->end()) continue;
    G4ThreeVector pos = physicalVolume->GetTranslation();
 
    G4ThreeVector globalPos_DNA = touchable->GetHistory()->GetTopTransform().Inverse().TransformPoint(pos);
 
    G4double r0 = (pos - localPos).mag();
    G4double irt = GetIndependentReactionTime(molConfA,dna_molConf,r0);
 
    if(irt>=0 && irt<timeMax - globalTime)
    {
 
      index = -1;
      for(size_t i=0;i<positionMap.size();++i){
        if(globalPos_DNA == positionMap[i].first){
          index = (G4int)i;
          break;
        }
      }
 
      G4Track* DNATrack;
      if(index == -1){
        auto DNAMol = new G4Molecule(dna_molConf);
        DNATrack = DNAMol->BuildTrack(globalTime,globalPos_DNA);
        DNATrack->SetTrackStatus(fAlive);
        fTrackHolder->Push(DNATrack);
        positionMap.push_back(std::make_pair(globalPos_DNA,DNATrack));
      }else{
        DNATrack = positionMap[index].second;
      }
 
      irt += globalTime;
      if(irt < minTime) minTime = irt;
#ifdef DEBUG
      G4cout<<irt<<'\t'<<globalPos_DNA<<'\t'<<molConfA->GetName()<<" "<<track->GetTrackID()<<'\t'<<dna_molConf->GetName()<<" "<<DNATrack->GetTrackID()<<'\n';
#endif
      fReactionSet->AddReaction(irt,track,DNATrack);
    }
  }
}

References G4ITReactionSet::AddReaction(), G4Molecule::BuildTrack(), G4DNAMolecularReactionTable::CanReactWith(), D(), fAlive, fGeometry, FindBin(), G4VDNAMolecularGeometry::FindNearbyMolecules(), fMolReactionTable, fNx, fNy, fNz, fRCutOff, fReactionSet, fStopButAlive, fTrackHolder, fXMax, fXMin, fYMax, fYMin, fZMax, fZMin, G4cout, G4UniformRand, G4MolecularConfiguration::GetDiffusionCoefficient(), G4Molecule::GetDiffusionCoefficient(), G4Scheduler::GetGlobalTime(), G4Track::GetGlobalTime(), G4VTouchable::GetHistory(), GetIndependentReactionTime(), G4VPhysicalVolume::GetLogicalVolume(), G4Molecule::GetMolecularConfiguration(), G4DNAMolecularMaterial::GetMolecularConfiguration(), GetMolecule(), G4Molecule::GetMolecule(), G4MolecularConfiguration::GetName(), G4Track::GetPosition(), G4DNAMolecularReactionData::GetReactant2(), G4DNAMolecularReactionTable::GetReactionData(), G4NavigationHistory::GetTopTransform(), G4Track::GetTouchable(), G4Track::GetTrackID(), G4VTouchable::GetVolume(), G4Scheduler::Instance(), G4MoleculeTable::Instance(), G4DNAMolecularMaterial::Instance(), G4AffineTransform::Inverse(), eplot::material, CLHEP::detail::n, pos, positionMap, G4ITTrackHolder::Push(), G4Track::SetTrackStatus(), G4INCL::DeJongSpin::shoot(), spaceBinned, timeMax, timeMin, G4AffineTransform::TransformPoint(), CLHEP::Hep3Vector::x(), xendIndex, xiniIndex, CLHEP::Hep3Vector::y(), yendIndex, yiniIndex, CLHEP::Hep3Vector::z(), zendIndex, and ziniIndex.

Referenced by IRTSampling(), and MakeReaction().

SetReactionModel()

void G4DNAIRT_geometries::SetReactionModel

(

G4VDNAReactionModel *

model

)

Definition at line 639 of file G4DNAIRT_geometries.cc.

{
  fpReactionModel = model;
}

References fpReactionModel.

SetReactionTable()

void G4VITReactionProcess::SetReactionTable ( const G4ITReactionTable *  pReactionTable )

virtualinherited

Definition at line 38 of file G4VITReactionProcess.cc.

{
    fpReactionTable = pReactionTable;
}

References G4VITReactionProcess::fpReactionTable.

SpaceBinning()

void G4DNAIRT_geometries::SpaceBinning

(

)

Definition at line 144 of file G4DNAIRT_geometries.cc.

                                      {
  G4double time_step = G4Scheduler::Instance()->GetPreviousTimeStep()*ps;
 
  auto it_begin = fTrackHolder->GetMainList()->begin();
  while(it_begin != fTrackHolder->GetMainList()->end()){
 
    // for diffusion
    G4double D = GetMolecule(*it_begin)->GetDiffusionCoefficient();
    G4double sqrt_2Dt = sqrt(2 * D * time_step);
 
    G4double x = G4RandGauss::shoot(0,sqrt_2Dt);
    G4double y = G4RandGauss::shoot(0,sqrt_2Dt);
    G4double z = G4RandGauss::shoot(0,sqrt_2Dt);
 
    G4ThreeVector position_ori = it_begin->GetPosition();
    G4ThreeVector position = position_ori + G4ThreeVector(x,y,z);
    it_begin->SetPosition(position);
    it_begin->SetGlobalTime(timeMax);
 
    if ( fXMin > position.x() ) fXMin = position.x(); 
    if ( fYMin > position.y() ) fYMin = position.y();
    if ( fZMin > position.z() ) fZMin = position.z();
 
    if ( fXMax < position.x() ) fXMax = position.x();
    if ( fYMax < position.y() ) fYMax = position.y();
    if ( fZMax < position.z() ) fZMax = position.z();
 
    ++it_begin;
  }
 
  fNx = G4int((fXMax-fXMin)/fRCutOff) == 0 ? 1 : G4int((fXMax-fXMin)/fRCutOff);
  fNy = G4int((fYMax-fYMin)/fRCutOff) == 0 ? 1 : G4int((fYMax-fYMin)/fRCutOff);
  fNz = G4int((fZMax-fZMin)/fRCutOff) == 0 ? 1 : G4int((fZMax-fZMin)/fRCutOff);
 
}

References G4FastList< OBJECT >::begin(), D(), G4FastList< OBJECT >::end(), fNx, fNy, fNz, fRCutOff, fTrackHolder, fXMax, fXMin, fYMax, fYMin, fZMax, fZMin, G4Molecule::GetDiffusionCoefficient(), G4ITTrackHolder::GetMainList(), GetMolecule(), G4Scheduler::GetPreviousTimeStep(), G4Scheduler::Instance(), ps, G4INCL::DeJongSpin::shoot(), and timeMax.

Referenced by Initialize().

TestReactibility()

G4bool G4DNAIRT_geometries::TestReactibility ( const G4Track &  , const G4Track &  , G4double  , G4bool    )

overridevirtual

Implements G4VITReactionProcess.

Definition at line 631 of file G4DNAIRT_geometries.cc.

{
  return true;
}

Field Documentation

erfc

G4ErrorFunction* G4DNAIRT_geometries::erfc

private

Definition at line 100 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), GetIndependentReactionTime(), SamplePDC(), and ~G4DNAIRT_geometries().

fGeometry

G4VDNAMolecularGeometry* G4DNAIRT_geometries::fGeometry

private

Definition at line 115 of file G4DNAIRT_geometries.hh.

Referenced by Initialize(), and Sampling().

fMolReactionTable

const G4DNAMolecularReactionTable*& G4DNAIRT_geometries::fMolReactionTable

protected

Definition at line 94 of file G4DNAIRT_geometries.hh.

Referenced by GetIndependentReactionTime(), Initialize(), MakeReaction(), and Sampling().

fName

G4String G4VITReactionProcess::fName

protectedinherited

Definition at line 91 of file G4VITReactionProcess.hh.

fNx

G4int G4DNAIRT_geometries::fNx

private

Definition at line 111 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), IRTSampling(), MakeReaction(), Sampling(), and SpaceBinning().

fNy

G4int G4DNAIRT_geometries::fNy

private

Definition at line 111 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), IRTSampling(), MakeReaction(), Sampling(), and SpaceBinning().

fNz

G4int G4DNAIRT_geometries::fNz

private

Definition at line 111 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), IRTSampling(), MakeReaction(), Sampling(), and SpaceBinning().

fpReactionModel

G4VDNAReactionModel* G4DNAIRT_geometries::fpReactionModel

protected

Definition at line 95 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), and SetReactionModel().

fpReactionTable

const G4ITReactionTable* G4VITReactionProcess::fpReactionTable = nullptr

protectedinherited

Definition at line 90 of file G4VITReactionProcess.hh.

Referenced by G4VITReactionProcess::SetReactionTable().

fRCutOff

G4double G4DNAIRT_geometries::fRCutOff

private

Definition at line 105 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), Sampling(), and SpaceBinning().

fReactionSet

G4ITReactionSet* G4DNAIRT_geometries::fReactionSet

private

Definition at line 99 of file G4DNAIRT_geometries.hh.

Referenced by Initialize(), and Sampling().

fTrackHolder

G4ITTrackHolder* G4DNAIRT_geometries::fTrackHolder

private

Definition at line 98 of file G4DNAIRT_geometries.hh.

Referenced by Initialize(), IRTSampling(), MakeReaction(), Sampling(), and SpaceBinning().

fXMax

G4double G4DNAIRT_geometries::fXMax

private

Definition at line 110 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), IRTSampling(), MakeReaction(), Sampling(), and SpaceBinning().

fXMin

G4double G4DNAIRT_geometries::fXMin

private

Definition at line 109 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), IRTSampling(), MakeReaction(), Sampling(), and SpaceBinning().

fYMax

G4double G4DNAIRT_geometries::fYMax

private

Definition at line 110 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), IRTSampling(), MakeReaction(), Sampling(), and SpaceBinning().

fYMin

G4double G4DNAIRT_geometries::fYMin

private

Definition at line 109 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), IRTSampling(), MakeReaction(), Sampling(), and SpaceBinning().

fZMax

G4double G4DNAIRT_geometries::fZMax

private

Definition at line 110 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), IRTSampling(), MakeReaction(), Sampling(), and SpaceBinning().

fZMin

G4double G4DNAIRT_geometries::fZMin

private

Definition at line 109 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), IRTSampling(), MakeReaction(), Sampling(), and SpaceBinning().

positionMap

std::vector<std::pair<G4ThreeVector,G4Track*> > G4DNAIRT_geometries::positionMap

private

Definition at line 103 of file G4DNAIRT_geometries.hh.

Referenced by Initialize(), and Sampling().

spaceBinned

std::map<G4int,std::map<G4int,std::map<G4int,std::vector<G4Track*> > > > G4DNAIRT_geometries::spaceBinned

private

Definition at line 102 of file G4DNAIRT_geometries.hh.

Referenced by Initialize(), IRTSampling(), MakeReaction(), and Sampling().

timeMax

G4double G4DNAIRT_geometries::timeMax

private

Definition at line 107 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), Sampling(), and SpaceBinning().

timeMin

G4double G4DNAIRT_geometries::timeMin

private

Definition at line 106 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), and Sampling().

xendIndex

G4int G4DNAIRT_geometries::xendIndex

private

Definition at line 113 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), and Sampling().

xiniIndex

G4int G4DNAIRT_geometries::xiniIndex

private

Definition at line 112 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), and Sampling().

yendIndex

G4int G4DNAIRT_geometries::yendIndex

private

Definition at line 113 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), and Sampling().

yiniIndex

G4int G4DNAIRT_geometries::yiniIndex

private

Definition at line 112 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), and Sampling().

zendIndex

G4int G4DNAIRT_geometries::zendIndex

private

Definition at line 113 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), and Sampling().

ziniIndex

G4int G4DNAIRT_geometries::ziniIndex

private

Definition at line 112 of file G4DNAIRT_geometries.hh.

Referenced by G4DNAIRT_geometries(), Initialize(), and Sampling().

---

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

• source/processes/electromagnetic/dna/models/include/G4DNAIRT_geometries.hh
• source/processes/electromagnetic/dna/models/src/G4DNAIRT_geometries.cc