Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/models/src/G4DiffusionControlledReactionModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/dna/models/src/G4DiffusionControlledReactionModel.cc (Version 11.3.0) and /processes/electromagnetic/dna/models/src/G4DiffusionControlledReactionModel.cc (Version 9.5.p1)


  1 // *******************************************      1 
  2 // * License and Disclaimer                       
  3 // *                                              
  4 // * The  Geant4 software  is  copyright of th    
  5 // * the Geant4 Collaboration.  It is provided    
  6 // * conditions of the Geant4 Software License    
  7 // * LICENSE and available at  http://cern.ch/    
  8 // * include a list of copyright holders.         
  9 // *                                              
 10 // * Neither the authors of this software syst    
 11 // * institutes,nor the agencies providing fin    
 12 // * work  make  any representation or  warran    
 13 // * regarding  this  software system or assum    
 14 // * use.  Please see the license in the file     
 15 // * for the full disclaimer and the limitatio    
 16 // *                                              
 17 // * This  code  implementation is the result     
 18 // * technical work of the GEANT4 collaboratio    
 19 // * By using,  copying,  modifying or  distri    
 20 // * any work based  on the software)  you  ag    
 21 // * use  in  resulting  scientific  publicati    
 22 // * acceptance of all terms of the Geant4 Sof    
 23 // *******************************************    
 24 //                                                
 25 //                                                
 26 // Author: Hoang TRAN                             
 27                                                   
 28 #include "G4DiffusionControlledReactionModel.h    
 29 #include "G4Track.hh"                             
 30 #include "G4DNAMolecularReactionTable.hh"         
 31 #include "G4PhysicalConstants.hh"                 
 32 #include "G4Exp.hh"                               
 33 #include "G4IRTUtils.hh"                          
 34 #include "G4SystemOfUnits.hh"                     
 35 #include "G4Electron_aq.hh"                       
 36 #include "Randomize.hh"                           
 37 #include "G4Molecule.hh"                          
 38 #include "G4ErrorFunction.hh"                     
 39                                                   
 40 G4DiffusionControlledReactionModel::G4Diffusio    
 41                                                   
 42 G4DiffusionControlledReactionModel::~G4Diffusi    
 43                                                   
 44 void G4DiffusionControlledReactionModel::Initi    
 45   const G4MolecularConfiguration* pMolecule, c    
 46 {                                                 
 47   fpReactionData = fpReactionTable->GetReactio    
 48 }                                                 
 49                                                   
 50 void G4DiffusionControlledReactionModel::Initi    
 51   const G4MolecularConfiguration* pMolecule)      
 52 {                                                 
 53   fpReactionData = fpReactionTable->GetReactio    
 54 }                                                 
 55                                                   
 56 G4double G4DiffusionControlledReactionModel::G    
 57   const G4MolecularConfiguration* pMol1, const    
 58 {                                                 
 59   auto reactionData = fpReactionTable->GetReac    
 60   if(reactionData == nullptr)                     
 61   {                                               
 62     G4ExceptionDescription exceptionDescriptio    
 63     exceptionDescription << "No reactionData"     
 64                          << " for : " << pMol1    
 65                          << pMol2->GetName();     
 66     G4Exception("G4DiffusionControlledReaction    
 67                 "::GetReactionRadius()",          
 68                 "G4DiffusionControlledReaction    
 69                 exceptionDescription);            
 70     return 0.;                                    
 71   }                                               
 72   return reactionData->GetEffectiveReactionRad    
 73 }                                                 
 74                                                   
 75 G4double G4DiffusionControlledReactionModel::G    
 76 {                                                 
 77   auto pMol1 = (*fpReactionData)[i]->GetReacta    
 78   auto pMol2 = (*fpReactionData)[i]->GetReacta    
 79   return GetReactionRadius(pMol1, pMol2);         
 80 }                                                 
 81                                                   
 82 G4double G4DiffusionControlledReactionModel::G    
 83   const G4Track& trackA, const G4Track& trackB    
 84 {                                                 
 85   auto pMolConfA = GetMolecule(trackA)->GetMol    
 86   auto pMolConfB = GetMolecule(trackB)->GetMol    
 87                                                   
 88   G4double D =                                    
 89     pMolConfA->GetDiffusionCoefficient() + pMo    
 90                                                   
 91   if(D == 0)                                      
 92   {                                               
 93     G4ExceptionDescription exceptionDescriptio    
 94     exceptionDescription << "The total diffusi    
 95                          << pMolConfA->GetName    
 96                          << pMolConfB->GetName    
 97     G4Exception("G4DiffusionControlledReaction    
 98                 "::GetTimeToEncounter()",         
 99                 "G4DiffusionControlledReaction    
100                 exceptionDescription);            
101   }                                               
102                                                   
103   auto reactionData = G4DNAMolecularReactionTa    
104     pMolConfA, pMolConfB);                        
105   G4double kobs     = reactionData->GetObserve    
106   G4double distance = (trackA.GetPosition() -     
107   G4double SmoluchowskiRadius = reactionData->    
108                                                   
109   if(distance == 0 || distance < SmoluchowskiR    
110   {                                               
111     G4ExceptionDescription exceptionDescriptio    
112     exceptionDescription << "distance = " << d    
113                          << " Reff = " << Smol    
114                          << " for : " << pMolC    
115                          << pMolConfB->GetName    
116     G4Exception("G4DiffusionControlledReaction    
117                 "::GetTimeToEncounter()",         
118                 "G4DiffusionControlledReaction    
119                 exceptionDescription);            
120   }                                               
121   else                                            
122   {                                               
123     G4double Winf  = SmoluchowskiRadius / dist    
124     G4double U     = G4UniformRand();             
125     G4double X     = 0;                           
126     G4double irt_1 = -1.0 * ps;                   
127     if(Winf > 0 && U < Winf)                      
128     {                                             
129       G4double erfcIn = G4ErrorFunction::erfcI    
130       if(erfcIn != 0)                             
131       {                                           
132         G4double d =                              
133           (distance - SmoluchowskiRadius) / er    
134         irt_1 = (1.0 / (4 * D)) * d * d;          
135       }                                           
136     }                                             
137                                                   
138     if(reactionData->GetReactionType() == 0)      
139     {                                             
140       return irt_1;                               
141     }                                             
142                                                   
143     if(irt_1 < 0)                                 
144     {                                             
145       return irt_1;                               
146     }                                             
147                                                   
148     G4double kdif = 4 * CLHEP::pi * D * Smoluc    
149                                                   
150     if(pMolConfA == pMolConfB)                    
151     {                                             
152       kdif /= 2;                                  
153     }                                             
154     G4double kact   = G4IRTUtils::GetKact(kobs    
155     G4double sumOfk = kact + kdif;                
156     if(sumOfk != 0)                               
157     {                                             
158       G4double rateFactor = kact / sumOfk;        
159       if(G4UniformRand() > rateFactor)            
160       {                                           
161         return -1.0 * ps;                         
162       }                                           
163       G4double Y = std::abs(G4RandGauss::shoot    
164                                                   
165       if(Y > 0)                                   
166       {                                           
167         X = -(G4Log(G4UniformRand())) / Y;        
168       }                                           
169       G4double f     = X * SmoluchowskiRadius     
170       G4double irt_2 = (f * f) / D;               
171       return irt_1 + irt_2;                       
172     }                                             
173   }                                               
174   return -1.0 * ps;                               
175 }                                                 
176