00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #include "G4DNAMolecularReaction.hh"
00040 #include "G4DNAMolecularReactionTable.hh"
00041 #include "G4VDNAReactionModel.hh"
00042 #include "G4ITManager.hh"
00043 #include "G4UnitsTable.hh"
00044
00045 using namespace std;
00046
00047 G4DNAMolecularReaction::G4DNAMolecularReaction():G4VITReactionProcess(),
00048 fMolReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
00049 {
00050
00051 fVerbose = 0;
00052 fReactionModel = 0;
00053 fReactionRadius = -1;
00054 fDistance = -1;
00055 }
00056
00057 G4DNAMolecularReaction::~G4DNAMolecularReaction()
00058 {
00059
00060 if(fpChanges) delete fpChanges;
00061 }
00062
00063 G4DNAMolecularReaction::G4DNAMolecularReaction(const G4DNAMolecularReaction& other):G4VITReactionProcess(other),
00064 fMolReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
00065 {
00066
00067 fVerbose = other.fVerbose ;
00068 fMolReactionTable = other.fMolReactionTable;
00069 fReactionModel = 0;
00070 fReactionRadius = -1;
00071 fDistance = -1;
00072 }
00073
00074 G4DNAMolecularReaction& G4DNAMolecularReaction::operator=(const G4DNAMolecularReaction& rhs)
00075 {
00076 if (this == &rhs) return *this;
00077
00078 fVerbose = rhs.fVerbose ;
00079 fMolReactionTable = rhs.fMolReactionTable;
00080 fReactionRadius = -1;
00081 fDistance = -1;
00082
00083
00084 return *this;
00085 }
00086
00087 G4bool G4DNAMolecularReaction::TestReactibility(const G4Track& trackA,
00088 const G4Track& trackB,
00089 const double currentStepTime,
00090 const double ,
00091 bool userStepTimeLimit)
00092 {
00093 G4Molecule* moleculeA = GetMolecule(trackA);
00094 G4Molecule* moleculeB = GetMolecule(trackB);
00095
00096 if(!fReactionModel)
00097 {
00098 G4ExceptionDescription exceptionDescription ("You have to give a reaction model to the molecular reaction process");
00099 G4Exception("G4DNAMolecularReaction::TestReactibility","MolecularReaction001",
00100 FatalErrorInArgument,exceptionDescription);
00101 return false;
00102 }
00103 if(fMolReactionTable==0)
00104 {
00105 G4ExceptionDescription exceptionDescription ("You have to give a reaction table to the molecular reaction process");
00106 G4Exception("G4DNAMolecularReaction::TestReactibility","MolecularReaction002",
00107 FatalErrorInArgument,exceptionDescription);
00108 return false;
00109 }
00110
00111
00112 fReactionRadius = -1 ;
00113 fReactionRadius = fReactionModel -> GetReactionRadius(moleculeA, moleculeB);
00114
00115 fDistance = -1 ;
00116
00117 if(currentStepTime == 0.)
00118 {
00119 userStepTimeLimit = false;
00120 }
00121
00122 G4bool output = fReactionModel->FindReaction(trackA, trackB, fReactionRadius,fDistance, userStepTimeLimit);
00123
00124 #ifdef G4VERBOSE
00125
00126 if(fVerbose > 1)
00127 {
00128 G4cout << "\033[1;39;36m" << "G4MolecularReaction "<< G4endl;
00129 G4cout << "FindReaction returned : " << G4BestUnit(output,"Length") << G4endl;
00130 G4cout << " reaction radius : " << G4BestUnit(fReactionRadius,"Length")
00131 << " real distance : " << G4BestUnit((trackA.GetPosition() - trackB.GetPosition()).mag(), "Length")
00132 << " calculated distance by model (= -1 if real distance > reaction radius and the user limitation step is not reached) : "
00133 << G4BestUnit(fDistance,"Length")
00134 << G4endl;
00135
00136 G4cout << "TrackID A : " << trackA.GetTrackID()
00137 << ", TrackID B : " << trackB.GetTrackID()
00138 << " | MolA " << moleculeA->GetName()
00139 << ", MolB " << moleculeB->GetName()
00140 <<"\033[0m\n"
00141 << G4endl;
00142 G4cout << "--------------------------------------------" << G4endl;
00143 }
00144 #endif
00145 return output ;
00146 }
00147
00148
00149 G4ITReactionChange* G4DNAMolecularReaction::MakeReaction(const G4Track& trackA, const G4Track& trackB)
00150 {
00151 fpChanges = new G4ITReactionChange();
00152 fpChanges->Initialize(trackA, trackB);
00153
00154 G4Molecule* moleculeA = GetMolecule(trackA);
00155 G4Molecule* moleculeB = GetMolecule(trackB);
00156
00157 #ifdef G4VERBOSE
00158
00159 if(fVerbose)
00160 {
00161 G4cout << "G4DNAMolecularReaction::MakeReaction" << G4endl;
00162 G4cout<<"TrackA n°" << trackA.GetTrackID()
00163 <<"\t | Track B n°" << trackB.GetTrackID() << G4endl;
00164
00165 G4cout<<"Track A : Position : " << G4BestUnit(trackA.GetPosition(),"Length")
00166 <<"\t Global Time : " << G4BestUnit(trackA.GetGlobalTime(), "Time")<< G4endl;
00167
00168 G4cout<<"Track B : Position : " << G4BestUnit(trackB.GetPosition() ,"Length")
00169 <<"\t Global Time : " << G4BestUnit(trackB.GetGlobalTime(), "Time")<< G4endl;
00170
00171 G4cout<<"Reaction range : " << G4BestUnit(fReactionRadius,"Length")
00172 << " \t Separation distance : " << G4BestUnit(fDistance,"Length")<< G4endl;
00173 G4cout << "--------------------------------------------" << G4endl;
00174 }
00175 #endif
00176
00177 const G4DNAMolecularReactionData* reactionData = fMolReactionTable->GetReactionData(moleculeA, moleculeB);
00178
00179 G4int nbProducts = reactionData->GetNbProducts();
00180
00181 if (nbProducts)
00182 {
00183 G4double D1 = moleculeA->GetDiffusionCoefficient();
00184 G4double D2 = moleculeB->GetDiffusionCoefficient();
00185 G4double sqrD1 = sqrt(D1);
00186 G4double sqrD2 = sqrt(D2);
00187 G4double numerator = sqrD1 + sqrD2;
00188 G4ThreeVector reactionSite = sqrD1/numerator * trackA.GetPosition()
00189 + sqrD2/numerator * trackB.GetPosition() ;
00190
00191 for (G4int j=0 ; j < nbProducts; j++)
00192 {
00193 G4Molecule* product = new G4Molecule(*reactionData->GetProduct(j));
00194 G4Track* productTrack = product->BuildTrack(trackA.GetGlobalTime(),
00195 reactionSite);
00196
00197 productTrack->SetTrackStatus(fAlive);
00198
00199 fpChanges->AddSecondary(productTrack);
00200 G4ITManager<G4Molecule>::Instance()->Push(productTrack);
00201 }
00202 }
00203
00204 fpChanges->KillParents(true);
00205
00206 return fpChanges;
00207 }