G4DNASmoluchowskiReactionModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id: G4DNASmoluchowskiReactionModel.cc 65695 2012-11-27 11:39:12Z gunter $
00027 //
00028 #include "G4DNASmoluchowskiReactionModel.hh"
00029 #include "Randomize.hh"
00030 #include "G4Track.hh"
00031 #include "G4DNAMolecularReactionTable.hh"
00032 
00033 G4DNASmoluchowskiReactionModel::G4DNASmoluchowskiReactionModel() : G4VDNAReactionModel()
00034 {
00035     fReactionData = 0 ;
00036 }
00037 
00038 G4DNASmoluchowskiReactionModel::G4DNASmoluchowskiReactionModel(const G4DNASmoluchowskiReactionModel& __right) :
00039     G4VDNAReactionModel(__right)
00040 {
00041     fReactionData = 0 ;
00042 }
00043 
00044 G4DNASmoluchowskiReactionModel& G4DNASmoluchowskiReactionModel::operator=(const G4DNASmoluchowskiReactionModel& right)
00045 {
00046     if(this == &right) return *this;
00047     fReactionData = 0;
00048     return *this;
00049 }
00050 
00051 G4DNASmoluchowskiReactionModel::~G4DNASmoluchowskiReactionModel()
00052 {
00053     fReactionData = 0 ;
00054 }
00055 
00056 void G4DNASmoluchowskiReactionModel::Initialise(const G4Molecule* __molecule, const G4Track&)
00057 {
00058     fReactionData = fReactionTable->GetReactionData(__molecule);
00059 }
00060 
00061 void G4DNASmoluchowskiReactionModel::InitialiseToPrint(const G4Molecule* __molecule)
00062 {
00063     fReactionData = fReactionTable->GetReactionData(__molecule);
00064 }
00065 
00066 G4double G4DNASmoluchowskiReactionModel::GetReactionRadius(const G4Molecule* mol1,
00067                                                                const G4Molecule* mol2)
00068 {
00069     G4double __output = fReactionTable -> GetReactionData(mol1,mol2)->GetReducedReactionRadius();
00070     return  __output ;
00071 }
00072 
00073 G4double G4DNASmoluchowskiReactionModel::GetReactionRadius(const G4int __i)
00074 {
00075     G4double __output = (*fReactionData)[__i] -> GetReducedReactionRadius();
00076     return  __output ;
00077 }
00078 
00079 G4bool G4DNASmoluchowskiReactionModel::FindReaction(const G4Track& __trackA,
00080         const G4Track& __trackB,
00081         const G4double __R,
00082         G4double& __r,
00083         const G4bool __alongStepReaction)
00084 {
00085     G4double __postStepSeparation = 0;
00086     bool __do_break = false ;
00087     G4double __R2 = __R*__R ;
00088     int k = 0 ;
00089 
00090     for(; k < 3 ; k++)
00091     {
00092         __postStepSeparation += std::pow(__trackA.GetPosition()[k] - __trackB.GetPosition()[k],2);
00093 
00094         if(__postStepSeparation > __R2)
00095         {
00096             __do_break = true  ;
00097             break ;
00098         }
00099     }
00100 
00101     if(__do_break == false)
00102     {
00103          // The loop was not break
00104          // => __r^2 < __R^2
00105         __r = std::sqrt(__postStepSeparation);
00106         return true;
00107     }
00108     else if(__alongStepReaction == true)
00109     {
00110         //G4cout << "alongStepReaction==true" << G4endl;
00111         //Along step cheack and
00112         // the loop has break
00113 
00114         // Continue loop
00115         for(; k < 3 ; k++)
00116         {
00117             __postStepSeparation += std::pow(__trackA.GetPosition()[k] - __trackB.GetPosition()[k],2);
00118         }
00119         // Use Green approach : the Brownian bridge
00120         __r = (__postStepSeparation = std::sqrt(__postStepSeparation) );
00121 
00122         G4Molecule* __moleculeA = GetMolecule(__trackA);
00123         G4Molecule* __moleculeB = GetMolecule(__trackB);
00124 
00125         G4double __D = __moleculeA->GetDiffusionCoefficient() + __moleculeB->GetDiffusionCoefficient();
00126 
00127         G4ThreeVector __preStepPositionA = __trackA.GetStep()->GetPreStepPoint() ->GetPosition();
00128         G4ThreeVector __preStepPositionB = __trackB.GetStep()->GetPreStepPoint() ->GetPosition();
00129 
00130         if(__preStepPositionA == __trackA.GetPosition())
00131         {
00132             G4ExceptionDescription exceptionDescription ;
00133             exceptionDescription << "The molecule : " <<  __moleculeA->GetName();
00134             exceptionDescription << " did not move since the previous step ";
00135             G4Exception("G4DNASmoluchowskiReactionModel::FindReaction","G4DNASmoluchowskiReactionModel",
00136                         FatalErrorInArgument,exceptionDescription);
00137         }
00138 
00139         G4double __preStepSeparation = (__preStepPositionA - __preStepPositionB).mag();
00140 
00141         G4double __probabiltyOfEncounter =  std::exp(-(__preStepSeparation - __R)*(__postStepSeparation - __R)
00142                                         / (__D* (__trackB.GetStep()->GetDeltaTime())));
00143         G4double __selectedPOE = G4UniformRand();
00144 
00145         if(__selectedPOE<=__probabiltyOfEncounter)  return true;
00146     }
00147 
00148     return false ;
00149 }

Generated on Mon May 27 17:48:04 2013 for Geant4 by  doxygen 1.4.7