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 "G4DNAMoleculeEncounterStepper.hh"
00040 #include "G4VDNAReactionModel.hh"
00041 #include "G4DNAMolecularReactionTable.hh"
00042 #include "G4H2O.hh"
00043 #include "G4UnitsTable.hh"
00044
00045 using namespace std;
00046
00047 G4DNAMoleculeEncounterStepper::G4DNAMoleculeEncounterStepper() :
00048 G4VITTimeStepper(),
00049 fMolecularReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable)),
00050 fReactionModel(0)
00051 {
00052 fVerbose = 0;
00053 fHasAlreadyReachedNullTime = false;
00054 }
00055
00056 G4DNAMoleculeEncounterStepper& G4DNAMoleculeEncounterStepper::operator=(const G4DNAMoleculeEncounterStepper& rhs)
00057 {
00058 if(this == &rhs) return *this;
00059 fReactionModel = 0;
00060 fVerbose = rhs.fVerbose;
00061 fMolecularReactionTable = rhs.fMolecularReactionTable;
00062 fHasAlreadyReachedNullTime = false;
00063 return *this;
00064 }
00065
00066 G4DNAMoleculeEncounterStepper::~G4DNAMoleculeEncounterStepper()
00067 {}
00068
00069 G4DNAMoleculeEncounterStepper::G4DNAMoleculeEncounterStepper(const G4DNAMoleculeEncounterStepper& right) :
00070 G4VITTimeStepper(right),
00071 fMolecularReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
00072 {
00073 fVerbose = right.fVerbose ;
00074 fMolecularReactionTable = right.fMolecularReactionTable;
00075 fReactionModel = 0;
00076 fHasAlreadyReachedNullTime = false;
00077 }
00078
00079 void G4DNAMoleculeEncounterStepper::Prepare()
00080 {
00081
00082
00083 G4VITTimeStepper::Prepare();
00084 G4ITManager<G4Molecule>::Instance()->UpdatePositionMap();
00085 }
00086
00087 G4double G4DNAMoleculeEncounterStepper::CalculateStep(const G4Track& trackA, const G4double& userMinTimeStep)
00088 {
00089
00090
00091
00092 G4Molecule* moleculeA = GetMolecule(trackA);
00093
00094
00095 #ifdef G4VERBOSE
00096 if(fVerbose)
00097 {
00098 G4cout << "_______________________________________________________________________" << G4endl;
00099 G4cout << "G4DNAMoleculeEncounterStepper::CalculateStep" << G4endl;
00100 G4cout << "Incident molecule : " << moleculeA->GetName()
00101 << " (" << trackA.GetTrackID() << ") "
00102 << G4endl;
00103 }
00104 #endif
00105
00106
00107
00108 const vector<const G4Molecule*>* reactivesVector =
00109 fMolecularReactionTable -> CanReactWith(moleculeA);
00110
00111 if(!reactivesVector)
00112 {
00113 #ifdef G4VERBOSE
00114
00115 if(fVerbose > 1)
00116 {
00117 G4cout << "!!!!!!!!!!!!!!!!!!!!"<<G4endl;
00118 G4cout << "!!! WARNING" << G4endl;
00119 G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity for the reaction because the molecule "
00120 << moleculeA->GetName()
00121 << " does not have any reactants given in the reaction table."
00122 << G4endl;
00123 G4cout << "!!!!!!!!!!!!!!!!!!!!"<<G4endl;
00124 }
00125 #endif
00126 return DBL_MAX;
00127 }
00128
00129 G4int nbReactives = reactivesVector->size();
00130
00131 if(nbReactives == 0)
00132 {
00133 #ifdef G4VERBOSE
00134
00135 if(fVerbose)
00136 {
00137 G4cout << "!!!!!!!!!!!!!!!!!!!!"<<G4endl;
00138 G4cout << "!!! WARNING" << G4endl;
00139 G4cout << "G4MoleculeEncounterStepper::CalculateStep will return infinity for the reaction because the molecule "
00140 << moleculeA->GetName()
00141 << " does not have any reactants given in the reaction table."
00142 << "This message can also result from a wrong implementation of the reaction table."
00143 << G4endl;
00144 G4cout << "!!!!!!!!!!!!!!!!!!!!"<<G4endl;
00145 }
00146 #endif
00147 return DBL_MAX;
00148 }
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159 fUserMinTimeStep = userMinTimeStep ;
00160 if(fReactants) fReactants = 0 ;
00161 fReactants = new vector<G4Track*>();
00162
00163 fSampledMinTimeStep = DBL_MAX;
00164 fHasAlreadyReachedNullTime = false;
00165
00166 fReactionModel -> Initialise(moleculeA, trackA) ;
00167
00168
00169
00170 for (G4int i=0 ; i<nbReactives ; i++)
00171 {
00172 const G4Molecule* moleculeB = (*reactivesVector)[i];
00173
00174
00175
00176 G4double R = -1 ;
00177 R = fReactionModel -> GetReactionRadius(i);
00178
00179
00180
00181 G4KDTreeResultHandle results (G4ITManager<G4Molecule>::Instance()
00182 -> FindNearest(moleculeA, moleculeB));
00183
00184 RetrieveResults(trackA,moleculeA,moleculeB,R,results, true);
00185 }
00186
00187 #ifdef G4VERBOSE
00188
00189 if(fVerbose)
00190 {
00191 G4cout << "G4MoleculeEncounterStepper::CalculateStep will finally return :"
00192 << G4BestUnit(fSampledMinTimeStep, "Time")<< G4endl;
00193
00194 if(fVerbose > 1)
00195 {
00196 G4cout << "TrackA: " << moleculeA->GetName() << " (" << trackA.GetTrackID() << ") can react with: ";
00197
00198 vector<G4Track*>::iterator it;
00199 for(it = fReactants->begin() ; it != fReactants->end() ; it++)
00200 {
00201 G4Track* trackB = *it;
00202 G4cout << GetMolecule(trackB)->GetName() << " ("
00203 << trackB->GetTrackID() << ") \t ";
00204 }
00205 G4cout << G4endl;
00206 }
00207 }
00208 #endif
00209 return fSampledMinTimeStep ;
00210 }
00211
00212
00213 void G4DNAMoleculeEncounterStepper::RetrieveResults(const G4Track& trackA, const G4Molecule* moleculeA,
00214 const G4Molecule* moleculeB, const G4double R,
00215 G4KDTreeResultHandle& results, G4bool iterate)
00216 {
00217 if(!results)
00218 {
00219 #ifdef G4VERBOSE
00220
00221 if(fVerbose > 1)
00222 {
00223 G4cout << "No molecule " << moleculeB->GetName()
00224 << " found to react with "
00225 << moleculeA->GetName()
00226 << G4endl;
00227 }
00228 #endif
00229 return ;
00230 }
00231
00232 G4double DA = moleculeA->GetDiffusionCoefficient() ;
00233 G4double DB = moleculeB->GetDiffusionCoefficient() ;
00234
00235 for(results->Rewind();
00236 !results->End();
00237 results->Next())
00238 {
00239
00240 G4IT* reactiveB = (G4IT*) results->GetItemData() ;
00241
00242 if (reactiveB==0)
00243 {
00244
00245
00246 continue ;
00247 }
00248
00249 G4Track *trackB = reactiveB->GetTrack();
00250
00251 if(trackB == 0)
00252 {
00253 G4ExceptionDescription exceptionDescription ;
00254 exceptionDescription << "The reactant B found using the ITManager does not have a valid track "
00255 << " attached to it. If this is done on purpose, please do not record this "
00256 << " molecule in the ITManager."
00257 << G4endl;
00258 G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults","MoleculeEncounterStepper001",
00259 FatalErrorInArgument,exceptionDescription);
00260 continue ;
00261 }
00262
00263 if (trackB->GetTrackStatus() != fAlive)
00264 {
00265 G4ExceptionDescription exceptionDescription ;
00266 exceptionDescription << "The track status of one of the nearby reactants is not fAlive" << G4endl;
00267 exceptionDescription << "The incomming trackID "
00268 << "(trackA entering in G4DNAMoleculeEncounterStepper and "
00269 << "for which you are looking reactant for) is : "
00270 << trackA.GetTrackID() <<"("<< GetMolecule(trackA)->GetName()<<")" << G4endl;
00271 exceptionDescription << "And the trackID of the reactant (trackB) is: "
00272 << trackB->GetTrackID() <<"("<< GetMolecule(trackB)->GetName()<<")" << G4endl;
00273 G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults","MoleculeEncounterStepper002",
00274 FatalErrorInArgument,exceptionDescription);
00275 continue ;
00276 }
00277
00278 if(trackB == &trackA)
00279 {
00280
00281 G4ExceptionDescription exceptionDescription ;
00282 exceptionDescription << "A track is reacting with itself (which is impossible) ie trackA == trackB"
00283 << G4endl ;
00284 exceptionDescription << "Molecule A (and B) is of type : "
00285 << moleculeA->GetName()
00286 << " with trackID : "
00287 << trackA.GetTrackID() << G4endl;
00288
00289 G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults","MoleculeEncounterStepper003",
00290 FatalErrorInArgument,exceptionDescription);
00291
00292 }
00293
00294 if(fabs(trackB->GetGlobalTime() - trackA.GetGlobalTime()) > trackA.GetGlobalTime()*(1-1/100) )
00295 {
00296
00297 G4ExceptionDescription exceptionDescription;
00298 exceptionDescription << "The interacting tracks are not synchronized in time"<< G4endl;
00299 exceptionDescription<< "trackB->GetGlobalTime() != trackA.GetGlobalTime()"
00300 << G4endl;
00301
00302 exceptionDescription << "trackA : trackID : " << trackA.GetTrackID()
00303 << "\t Name :" << moleculeA->GetName()
00304 <<"\t trackA->GetGlobalTime() = "
00305 << G4BestUnit(trackA.GetGlobalTime(), "Time") << G4endl;
00306
00307 exceptionDescription << "trackB : trackID : " << trackB->GetTrackID()
00308 << "\t Name :" << moleculeB->GetName()
00309 << "\t trackB->GetGlobalTime() = "
00310 << G4BestUnit(trackB->GetGlobalTime(), "Time")<< G4endl;
00311
00312 G4Exception("G4DNAMoleculeEncounterStepper::RetrieveResults","MoleculeEncounterStepper004",
00313 FatalErrorInArgument,exceptionDescription);
00314 }
00315
00316 G4double r2 = results->GetDistanceSqr() ;
00317 #ifdef G4VERBOSE
00318 if(fVerbose > 1)
00319 {
00320 G4cout << "\t ************************************************** " << G4endl;
00321 G4cout <<"\t Reaction between "
00322 << moleculeA->GetName() << " (" << trackA.GetTrackID() << ") "
00323 << " & " << moleculeB->GetName() << " (" << trackB->GetTrackID() << "), "
00324 << "Interaction Range = "
00325 << G4BestUnit(R, "Length")<<G4endl;
00326 G4cout <<"\t Real distance between reactants = "
00327 << G4BestUnit((trackA.GetPosition() - trackB->GetPosition()).mag(), "Length")<<G4endl;
00328 G4cout <<"\t Distance between reactants calculated by nearest neighbor algorithm = "
00329 << G4BestUnit(sqrt(r2), "Length")<<G4endl;
00330
00331 }
00332 #endif
00333
00334 if(r2 <= R*R)
00335 {
00336
00337
00338
00339
00340
00341 if(fHasAlreadyReachedNullTime == false)
00342 {
00343 fReactants->clear();
00344 fHasAlreadyReachedNullTime = true;
00345 }
00346
00347 if(iterate)
00348 {
00349 fSampledMinTimeStep = 0.;
00350 G4KDTreeResultHandle results2 (G4ITManager<G4Molecule>::Instance()
00351 -> FindNearestInRange(moleculeA, moleculeB,R));
00352 RetrieveResults(trackA, moleculeA, moleculeB, R, results2, false);
00353 }
00354 else
00355 {
00356 fReactants->push_back(trackB);
00357 }
00358
00359 continue;
00360 }
00361 else
00362 {
00363 G4double r = sqrt(r2);
00364 G4double tempMinET = pow(r - R,2)
00365 /(16 * (DA + DB + 2*sqrt(DA*DB)));
00366
00367 if(tempMinET <= fSampledMinTimeStep)
00368 {
00369 if(tempMinET <= fUserMinTimeStep)
00370 {
00371 if(fSampledMinTimeStep > fUserMinTimeStep)
00372 fReactants->clear();
00373 fSampledMinTimeStep = fUserMinTimeStep;
00374 fReactants->push_back(trackB);
00375
00376 }
00377 else
00378 {
00379 fSampledMinTimeStep = tempMinET;
00380 if(tempMinET < fSampledMinTimeStep)
00381 fReactants->clear();
00382 fReactants->push_back(trackB);
00383 }
00384 }
00385 }
00386 }
00387 }