53fpReactionModel(nullptr),
173 const vector<const G4MolecularConfiguration*>* reactivesVector =
176 if(reactivesVector ==
nullptr)
return;
192 std::vector<G4Track*> spaceBin =
spaceBinned[ii][jj][kk];
194 if(!spaceBin[
n] || track == spaceBin[
n])
continue;
203 auto it = std::find(reactivesVector->begin(), reactivesVector->end(), molConfB);
204 if(it == reactivesVector->end())
continue;
214 sigma = std::sqrt(2.0 * diffusionCoefficient * dt);
221 }
else if(dt < 0)
continue;
227 if(irt>=0 && irt<
timeMax - globalTime)
230 if(irt < minTime) minTime = irt;
247 for(
size_t u=0; u<fReactionDatas->size();u++){
248 if((*fReactionDatas)[u]->GetReactant2()->GetDiffusionCoefficient() == 0){
250 if(kObs == 0)
continue;
252 if( time < minTime && time >= globalTime && time <
timeMax){
261 G4cout<<
"scavenged: "<<minTime<<
'\t'<<molConfA->
GetName()<<it_begin->GetTrackID()<<
'\n';
272 const auto pMoleculeA = molA;
273 const auto pMoleculeB = molB;
277 if(r0 == 0) r0 += 1e-3*
nm;
281 if(
D == 0)
D += 1e-20*(
m2/
s);
282 G4double rc = fReactionData->GetOnsagerRadius();
284 if ( reactionType == 0){
285 G4double sigma = fReactionData->GetEffectiveReactionRadius();
287 if(sigma > r0)
return 0;
288 if( rc != 0) r0 = -rc / (1-std::exp(rc/r0));
293 if ( W > 0 && W < Winf ) irt = (0.25/
D) * std::pow( (r0-sigma)/
erfc->
erfcInv(r0*W/sigma), 2 );
297 else if ( reactionType == 1 ){
298 G4double sigma = fReactionData->GetReactionRadius();
299 G4double kact = fReactionData->GetActivationRateConstant();
300 G4double kdif = fReactionData->GetDiffusionRateConstant();
301 G4double kobs = fReactionData->GetObservedReactionRateConstant();
306 a = 1/sigma * kact / kobs;
307 b = (r0 - sigma) / 2;
311 a = 4*pow(sigma,2)*
alpha/(
D*pow(rc,2))*pow(sinh(rc/(2*sigma)),2);
312 b = rc/4*(cosh(rc/(2*r0))/sinh(rc/(2*r0))-cosh(rc/(2*sigma))/sinh(rc/(2*sigma)));
313 r0 = -rc/(1-std::exp(rc/r0));
314 sigma = fReactionData->GetEffectiveReactionRadius();
318 if(fReactionData->GetProbability() >
G4UniformRand())
return 0;
321 Winf = sigma / r0 * kobs / kdif;
335 else if ( value >= xmax)
338 bin =
G4int(
n * ( value - xmin )/( xmax - xmin ) );
340 if ( bin < 0 ) bin = 0;
341 if ( bin >=
n ) bin =
n-1;
348 G4double p = 2.0 * std::sqrt(2.0*b/a);
349 G4double q = 2.0 / std::sqrt(2.0*b/a);
359 if ( U < p/(p + q *
M) ) X = pow(U * (p + q *
M) / 2, 2);
360 else X = pow(2/((1-U)*(p+q*
M)/
M),2);
364 lambdax = std::exp(-b*b/X) * ( 1.0 - a * std::sqrt(
CLHEP::pi * X) *
erfc->
erfcx(b/std::sqrt(X) + a*std::sqrt(X)));
366 if ((X <= 2.0*b/a && U <= lambdax) ||
367 (X >= 2.0*b/a && U*
M/X <= lambdax))
break;
371 if ( ntrials > 10000 ){
372 G4cout<<
"Totally rejected"<<
'\n';
384 pChanges->Initialize(trackA, trackB);
391 G4double effectiveReactionRadius = pReactionData->GetEffectiveReactionRadius();
393 const G4double D1 = pMoleculeA->GetDiffusionCoefficient();
394 const G4double D2 = pMoleculeB->GetDiffusionCoefficient();
405 S1.
setMag(effectiveReactionRadius);
409 if(dt != 0 && (D1 + D2) != 0 && r0 != 0){
412 if(s12 == 0) r2 = r1;
413 else if(s22 == 0) r1 = r2;
415 G4double alpha = effectiveReactionRadius * r0 / (2*(D1 + D2)*dt);
426 r1 = (D1 * S1 + D2 * S2) / (D1 + D2);
427 r2 = D2 * (S2 - S1) / (D1 + D2);
431 auto pTrackA =
const_cast<G4Track*
>(pChanges->GetTrackA());
432 auto pTrackB =
const_cast<G4Track*
>(pChanges->GetTrackB());
435 pTrackB->SetPosition(r2);
437 pTrackA->SetGlobalTime(globalTime);
438 pTrackB->SetGlobalTime(globalTime);
443 const G4int nbProducts = pReactionData->GetNbProducts();
447 const G4double sqrD1 = D1 == 0. ? 0. : std::sqrt(D1);
448 const G4double sqrD2 = D2 == 0. ? 0. : std::sqrt(D2);
449 if((sqrD1 + sqrD2) == 0){
452 const G4double inv_numerator = 1./(sqrD1 + sqrD2);
456 std::vector<G4ThreeVector>
position;
460 }
else if(nbProducts == 2){
463 }
else if (nbProducts == 3){
469 for(
G4int u = 0; u < nbProducts; u++){
471 auto product =
new G4Molecule(pReactionData->GetProduct(u));
472 auto productTrack = product->BuildTrack(globalTime,
475 productTrack->SetTrackStatus(
fAlive);
478 pChanges->AddSecondary(productTrack);
491 pChanges->KillParents(
true);
502 std::vector<std::unique_ptr<G4ITReactionChange>> fReactionInfo;
503 fReactionInfo.clear();
505 if (pReactionSet ==
nullptr)
507 return fReactionInfo;
511 assert(fReactionsetInTime.begin() != fReactionsetInTime.end());
513 auto it_begin = fReactionsetInTime.begin();
514 while(it_begin != fReactionsetInTime.end())
516 G4double irt = it_begin->get()->GetTime();
518 if(fGlobalTime < irt)
break;
522 G4Track* pTrackA = it_begin->get()->GetReactants().first;
523 G4Track* pTrackB = it_begin->get()->GetReactants().second;
524 auto pReactionChange =
MakeReaction(*pTrackA, *pTrackB);
527 fReactionInfo.push_back(std::move(pReactionChange));
531 it_begin = fReactionsetInTime.begin();
534 return fReactionInfo;
G4double D(G4double temp)
static const G4double alpha
G4Molecule * GetMolecule(const G4Track &track)
ReturnType & reference_cast(OriginalType &source)
static constexpr double nm
static constexpr double rad
static constexpr double s
static constexpr double ps
static constexpr double m2
CLHEP::Hep3Vector G4ThreeVector
G4GLOB_DLL std::ostream G4cout
G4ITTrackHolder * fTrackHolder
G4ITReactionSet * fReactionSet
void SetReactionModel(G4VDNAReactionModel *)
G4bool TestReactibility(const G4Track &, const G4Track &, G4double, G4bool) override
G4double SamplePDC(G4double, G4double)
const G4DNAMolecularReactionTable *& fMolReactionTable
void Initialize() override
std::map< G4int, std::map< G4int, std::map< G4int, std::vector< G4Track * > > > > spaceBinned
G4int FindBin(G4int, G4double, G4double, G4double)
G4double GetIndependentReactionTime(const G4MolecularConfiguration *, const G4MolecularConfiguration *, G4double)
std::vector< std::unique_ptr< G4ITReactionChange > > FindReaction(G4ITReactionSet *, const G4double, const G4double, const G4bool) override
std::unique_ptr< G4ITReactionChange > MakeReaction(const G4Track &, const G4Track &) override
G4VDNAReactionModel * fpReactionModel
G4int GetReactionType() const
G4double GetObservedReactionRateConstant() const
Data * GetReactionData(Reactant *, Reactant *) const
const ReactantList * CanReactWith(Reactant *) const
static G4double erfcx(G4double x)
static G4double erfcInv(G4double x)
void AddReaction(G4double time, G4Track *trackA, G4Track *trackB)
static G4ITReactionSet * Instance()
void SelectThisReaction(G4ITReactionPtr reaction)
G4ITReactionPerTime & GetReactionsPerTime()
G4TrackList * GetMainList(Key)
virtual void Push(G4Track *)
void MergeSecondariesWithMainList()
static G4ITTrackHolder * Instance()
const G4String & GetName() const
G4double GetDiffusionCoefficient() const
static G4Molecule * GetMolecule(const G4Track *)
const G4MolecularConfiguration * GetMolecularConfiguration() const
G4Track * BuildTrack(G4double globalTime, const G4ThreeVector &Position)
G4double GetDiffusionCoefficient() const
static G4Scheduler * Instance()
G4double GetEndTime() const
G4double GetGlobalTime() const
G4double GetStartTime() const
void SetPosition(const G4ThreeVector &aValue)
const G4ThreeVector & GetPosition() const
G4double GetGlobalTime() const
static constexpr double pi
ThreeVector shoot(const G4int Ap, const G4int Af)
T max(const T t1, const T t2)
brief Return the largest of the two arguments