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 <iomanip>
00040
00041 #include "G4DNAMolecularReactionTable.hh"
00042 #include "G4PhysicalConstants.hh"
00043 #include "G4SystemOfUnits.hh"
00044 #include "G4UIcommand.hh"
00045 #include "G4VDNAReactionModel.hh"
00046 #include "G4MoleculeHandleManager.hh"
00047
00048 using namespace std;
00049
00050 G4DNAMolecularReactionTable* G4DNAMolecularReactionTable::fInstance(0);
00051
00052 G4DNAMolecularReactionData::G4DNAMolecularReactionData():
00053 fReactive1(),fReactive2(),
00054 fReactionRate(0.),fReducedReactionRadius(0.),
00055 fProducts(0)
00056 {;}
00057
00058 G4DNAMolecularReactionData::G4DNAMolecularReactionData(G4double reactionRate,
00059 const G4Molecule* reactive1,
00060 const G4Molecule* reactive2):fProducts(0)
00061 {
00062 fReactionRate = reactionRate;
00063 SetReactive1(reactive1);
00064 SetReactive2(reactive2);
00065
00066 G4double sumDiffCoeff(0.);
00067
00068 if(*reactive1 == *reactive2)
00069 {
00070 sumDiffCoeff = reactive1->GetDiffusionCoefficient();
00071 fReducedReactionRadius = fReactionRate/(4*pi* reactive1->GetDiffusionCoefficient() * Avogadro);
00072 }
00073 else
00074 {
00075 sumDiffCoeff = reactive1->GetDiffusionCoefficient() + reactive2->GetDiffusionCoefficient();
00076 fReducedReactionRadius = fReactionRate/(4*pi* sumDiffCoeff * Avogadro);
00077 }
00078 }
00079
00080 G4DNAMolecularReactionData::~G4DNAMolecularReactionData()
00081 {
00082 if(fProducts)
00083 {
00084 fProducts->clear() ;
00085 delete fProducts;
00086 fProducts = 0;
00087 }
00088 }
00089
00090 void G4DNAMolecularReactionData::SetReactive1(const G4Molecule* reactive)
00091 {
00092 fReactive1 = G4MoleculeHandleManager::Instance()->GetMoleculeHandle(reactive);
00093 }
00094 void G4DNAMolecularReactionData::SetReactive2(const G4Molecule* reactive)
00095 {
00096 fReactive2 = G4MoleculeHandleManager::Instance()->GetMoleculeHandle(reactive);
00097 }
00098 void G4DNAMolecularReactionData::SetReactive(const G4Molecule* reactive1,
00099 const G4Molecule* reactive2)
00100 {
00101 fReactive1 = G4MoleculeHandleManager::Instance()->GetMoleculeHandle(reactive1);
00102 fReactive2 = G4MoleculeHandleManager::Instance()->GetMoleculeHandle(reactive2);
00103 }
00104
00105 void G4DNAMolecularReactionData::AddProduct(const G4Molecule* molecule)
00106 {
00107 if(!fProducts) fProducts = new std::vector<G4MoleculeHandle>();
00108 fProducts->push_back(G4MoleculeHandleManager::Instance()->GetMoleculeHandle(molecule));
00109 }
00110
00111 G4DNAMolecularReactionTable* G4DNAMolecularReactionTable::GetReactionTable()
00112 {
00113 if(!fInstance)
00114 {
00115 fInstance = new G4DNAMolecularReactionTable();
00116 }
00117 return fInstance;
00118 }
00119
00120 void G4DNAMolecularReactionTable::DeleteInstance()
00121 {
00122
00123
00124 if(fInstance)
00125 delete fInstance;
00126 fInstance = 0;
00127 }
00128
00129 G4DNAMolecularReactionTable::G4DNAMolecularReactionTable() : G4ITReactionTable(),
00130 fMoleculeHandleManager(G4MoleculeHandleManager::Instance())
00131 {
00132
00133 fVerbose = false;
00134 return;
00135 }
00136
00137 G4DNAMolecularReactionTable::~G4DNAMolecularReactionTable()
00138 {
00139
00140
00141 ReactionDataMap::iterator it1 = fReactionData.begin();
00142
00143 std::map<const G4Molecule*,
00144 const G4DNAMolecularReactionData*,
00145 compMoleculeP>::iterator it2;
00146
00147 for(;it1!=fReactionData.end();it1++)
00148 {
00149 for(it2 = it1->second.begin();it2 != it1->second.end();it2++)
00150 {
00151 const G4DNAMolecularReactionData* reactionData = it2->second;
00152 if(reactionData)
00153 {
00154 const G4Molecule* reactive1 = reactionData->GetReactive1();
00155 const G4Molecule* reactive2 = reactionData->GetReactive2();
00156
00157 fReactionData[reactive1][reactive2] = 0;
00158 fReactionData[reactive2][reactive1] = 0;
00159
00160 delete reactionData;
00161 }
00162 }
00163 }
00164
00165 fReactionDataMV.clear();
00166 fReactionData.clear();
00167 fReactivesMV.clear();
00168 }
00169
00170 void G4DNAMolecularReactionTable::SetReaction(G4DNAMolecularReactionData* reactionData)
00171 {
00172 const G4Molecule* reactive1 = reactionData->GetReactive1() ;
00173 const G4Molecule* reactive2 = reactionData->GetReactive2() ;
00174
00175 fReactionData[reactive1][reactive2] = reactionData;
00176 fReactivesMV[reactive1].push_back(reactive2);
00177 fReactionDataMV[reactive1].push_back(reactionData);
00178
00179 if(reactive1 != reactive2)
00180 {
00181 fReactionData[reactive2][reactive1] = reactionData;
00182 fReactivesMV[reactive2].push_back(reactive1);
00183 fReactionDataMV[reactive2].push_back(reactionData);
00184 }
00185 }
00186
00187 void G4DNAMolecularReactionTable::SetReaction(G4double reactionRate,
00188 const G4Molecule* reactive1,
00189 const G4Molecule* reactive2)
00190 {
00191 G4DNAMolecularReactionData* reactionData = new G4DNAMolecularReactionData(reactionRate, reactive1, reactive2);
00192 SetReaction(reactionData);
00193 }
00194
00195 void G4DNAMolecularReactionTable::PrintTable(G4VDNAReactionModel* pReactionModel)
00196 {
00197
00198
00199 if(pReactionModel)
00200 {
00201 if(!(pReactionModel->GetReactionTable()))
00202 pReactionModel -> SetReactionTable(this);
00203 }
00204
00205 ReactivesMV::iterator itReactives;
00206
00207 map<G4Molecule*,map<G4Molecule*, G4bool> > alreadyPrint;
00208
00209 G4cout<<"Nombre particules intervenants dans les reactions = "<< fReactivesMV.size() <<G4endl;
00210
00211 G4int nbPrintable = fReactivesMV.size()*fReactivesMV.size();
00212
00213 G4String *outputReaction = new G4String[nbPrintable];
00214 G4String *outputReactionRate = new G4String[nbPrintable];
00215 G4String *outputRange = new G4String[nbPrintable];
00216 G4int n = 0;
00217
00218 for(itReactives = fReactivesMV.begin() ; itReactives != fReactivesMV.end() ; itReactives++)
00219 {
00220 G4Molecule* moleculeA = (G4Molecule*) itReactives->first;
00221 const vector<const G4Molecule*>* reactivesVector = CanReactWith(moleculeA);
00222
00223 if(pReactionModel)
00224 pReactionModel -> InitialiseToPrint(moleculeA);
00225
00226 G4int nbReactants = fReactivesMV[itReactives->first].size();
00227
00228 for(G4int iReact = 0 ; iReact < nbReactants ; iReact++)
00229 {
00230
00231 G4Molecule* moleculeB = (G4Molecule*) (*reactivesVector)[iReact];
00232
00233 const G4DNAMolecularReactionData* reactionData = fReactionData[moleculeA][moleculeB];
00234
00235
00236
00237 if(!alreadyPrint[moleculeA][moleculeB])
00238 {
00239 outputReaction[n]=
00240 moleculeA->GetName()
00241 +" + " +
00242 moleculeB->GetName();
00243
00244 G4int nbProducts = reactionData->GetNbProducts();
00245
00246 if(nbProducts)
00247 {
00248 outputReaction[n] += " -> "+ reactionData->GetProduct(0)->GetName();
00249
00250 for(G4int j = 1 ; j < nbProducts ; j++)
00251 {
00252 outputReaction[n]+=" + "+reactionData->GetProduct(j)->GetName();
00253 }
00254 }
00255 else
00256 {
00257 outputReaction[n]+=" -> No product";
00258 }
00259
00260
00261
00262 outputReactionRate[n] = G4UIcommand::ConvertToString(reactionData->GetReactionRate()/(1e-3*m3/(mole*s)));
00263
00264
00265
00266 G4double interactionRange = -1;
00267 if(pReactionModel)
00268 interactionRange = pReactionModel->GetReactionRadius(iReact);
00269
00270 if(interactionRange!=-1)
00271 {
00272 outputRange[n] = G4UIcommand::ConvertToString(interactionRange/nanometer);
00273 }
00274 else
00275 {
00276 outputRange[n] = "";
00277 }
00278
00279 alreadyPrint[moleculeB][moleculeA] = TRUE;
00280 n++;
00281 }
00282 }
00283 }
00284 G4cout<<"Number of possible reactions: "<< n << G4endl;
00285
00287
00288
00290
00291 G4int maxlengthOutputReaction = -1;
00292 G4int maxlengthOutputReactionRate = -1;
00293
00294 for(G4int i = 0 ; i < n ; i++)
00295 {
00296 if(maxlengthOutputReaction < (G4int) outputReaction[i].length())
00297 {
00298 maxlengthOutputReaction = outputReaction[i].length();
00299 }
00300 if(maxlengthOutputReactionRate < (G4int)outputReactionRate[i].length())
00301 {
00302 maxlengthOutputReactionRate = outputReactionRate[i].length();
00303 }
00304 }
00305
00306 maxlengthOutputReaction+=2;
00307 maxlengthOutputReactionRate+=2;
00308
00309 if(maxlengthOutputReaction<10) maxlengthOutputReaction = 10;
00310 if(maxlengthOutputReactionRate<30) maxlengthOutputReactionRate = 30;
00311
00312 G4String title[3];
00313
00314 title[0] = "Reaction";
00315 title[1] = "Reaction Rate [dm3/(mol*s)]";
00316 title[2] = "Interaction Range for chosen reaction model";
00317
00318 G4cout<< setfill(' ')
00319 << setw(maxlengthOutputReaction) << left << title[0]
00320 << setw(maxlengthOutputReactionRate) << left << title[1]
00321 << setw(2) << left << title[2]
00322 << G4endl;
00323
00324 G4cout.fill('-');
00325 G4cout.width(maxlengthOutputReaction+2+maxlengthOutputReactionRate+2+(G4int)title[2].length());
00326 G4cout<<"-"<<G4endl;
00327 G4cout.fill(' ');
00328
00329 for(G4int i = 0 ; i < n ; i ++)
00330 {
00331 G4cout<< setw(maxlengthOutputReaction)<< left << outputReaction[i]
00332 << setw(maxlengthOutputReactionRate) << left << outputReactionRate[i]
00333 << setw(2) << left <<outputRange[i]
00334 <<G4endl;
00335
00336 G4cout.fill('-');
00337 G4cout.width(maxlengthOutputReaction+2+maxlengthOutputReactionRate+2+(G4int)title[2].length());
00338 G4cout<<"-"<<G4endl;
00339 G4cout.fill(' ');
00340 }
00341
00342 delete [] outputReaction;
00343 delete [] outputReactionRate;
00344 delete [] outputRange;
00345 }
00346
00347
00348
00349 const G4DNAMolecularReactionData*
00350 G4DNAMolecularReactionTable::GetReactionData(const G4Molecule* reactive1,
00351 const G4Molecule* reactive2) const
00352 {
00353 if(fReactionData.empty())
00354 {
00355 G4String errMsg = "No reaction table was implemented";
00356 G4Exception("G4MolecularInteractionTable::CanInteractWith","",FatalErrorInArgument, errMsg);
00357 return 0;
00358 }
00359
00360 ReactionDataMap::const_iterator it1 = fReactionData.find(reactive1);
00361
00362 if(it1 == fReactionData.end())
00363 {
00364 G4cout<<"Nom : " << reactive1->GetName()<<G4endl;
00365 G4String errMsg = "No reaction table was implemented for this molecule Definition : "
00366 + reactive1 -> GetName();
00367 G4Exception("G4MolecularInteractionTable::CanReactWith","",FatalErrorInArgument, errMsg);
00368 }
00369
00370 std::map<const G4Molecule*,
00371 const G4DNAMolecularReactionData*,
00372 compMoleculeP>::const_iterator it2 = it1->second.find(reactive2);
00373
00374 if(it2 == it1->second.end())
00375 {
00376 G4cout<<"Nom : " << reactive2->GetName()<<G4endl;
00377 G4String errMsg = "No reaction table was implemented for this molecule Definition : "
00378 + reactive2 -> GetName();
00379 G4Exception("G4MolecularInteractionTable::CanReactWith","",FatalErrorInArgument, errMsg);
00380 }
00381
00382 return (it2->second);
00383 }
00384
00385 const std::vector<const G4Molecule*>*
00386 G4DNAMolecularReactionTable::CanReactWith(const G4Molecule * aMolecule) const
00387 {
00388 if(fReactivesMV.empty())
00389 {
00390 G4String errMsg = "No reaction table was implemented";
00391 G4Exception("G4MolecularInteractionTable::CanReactWith","",FatalErrorInArgument, errMsg);
00392 return 0;
00393 }
00394
00395 ReactivesMV::const_iterator itReactivesMap = fReactivesMV.find(aMolecule) ;
00396
00397 if(itReactivesMap == fReactivesMV.end())
00398 {
00399 G4cout<<"Nom : " << aMolecule->GetName()<<G4endl;
00400 G4String errMsg = "No reaction table was implemented for this molecule Definition : "
00401 + aMolecule -> GetName();
00402 G4Exception("G4MolecularInteractionTable::CanReactWith","",FatalErrorInArgument, errMsg);
00403 return 0;
00404 }
00405 else
00406 {
00407 if(fVerbose)
00408 {
00409 G4cout<< " G4MolecularInteractionTable::CanReactWith :"<<G4endl;
00410 G4cout<<"You are checking reactants for : " << aMolecule->GetName()<<G4endl;
00411 G4cout<<" the number of reactants is : " << itReactivesMap->second.size()<<G4endl;
00412
00413 std::vector<const G4Molecule*>::const_iterator itProductsVector =
00414 itReactivesMap->second.begin();
00415
00416 for( ; itProductsVector != itReactivesMap->second.end(); itProductsVector++)
00417 {
00418 G4cout<<(*itProductsVector)->GetName()<<G4endl;
00419 }
00420 }
00421 return &(itReactivesMap->second);
00422 }
00423 return 0;
00424 }
00425
00426
00427 const std::map<const G4Molecule*, const G4DNAMolecularReactionData*, compMoleculeP>*
00428 G4DNAMolecularReactionTable::GetReativesNData(const G4Molecule* molecule) const
00429 {
00430
00431 if(fReactionData.empty())
00432 {
00433 G4String errMsg = "No reaction table was implemented";
00434 G4Exception("G4MolecularInteractionTable::CanInteractWith","",FatalErrorInArgument, errMsg);
00435 return 0;
00436 }
00437
00438 ReactionDataMap::const_iterator itReactivesMap = fReactionData.find(molecule) ;
00439
00440 if(itReactivesMap == fReactionData.end())
00441 {
00442 G4cout<<"Nom : " << molecule->GetName()<<G4endl;
00443 G4String errMsg = "No reaction table was implemented for this molecule Definition : "
00444 + molecule -> GetName();
00445 G4Exception("G4MolecularInteractionTable::CanReactWith","",FatalErrorInArgument, errMsg);
00446 }
00447 else
00448 {
00449 if(fVerbose)
00450 {
00451 G4cout<< " G4MolecularInteractionTable::CanReactWith :"<<G4endl;
00452 G4cout<<"You are checking reactants for : " << molecule->GetName()<<G4endl;
00453 G4cout<<" the number of reactants is : " << itReactivesMap->second.size()<<G4endl;
00454
00455 std::map<const G4Molecule*,
00456 const G4DNAMolecularReactionData*,
00457 compMoleculeP>::const_iterator itProductsVector =
00458 itReactivesMap->second.begin();
00459
00460 for( ; itProductsVector != itReactivesMap->second.end(); itProductsVector++)
00461 {
00462 G4cout<<itProductsVector->first->GetName()<<G4endl;
00463 }
00464 }
00465 return &(itReactivesMap->second);
00466 }
00467
00468 return 0;
00469 }
00470
00471 const std::vector<const G4DNAMolecularReactionData*>*
00472 G4DNAMolecularReactionTable::GetReactionData(const G4Molecule* molecule) const
00473 {
00474 if(fReactionDataMV.empty())
00475 {
00476 G4String errMsg = "No reaction table was implemented";
00477 G4Exception("G4MolecularInteractionTable::CanInteractWith","",FatalErrorInArgument, errMsg);
00478 return 0 ;
00479 }
00480 ReactionDataMV::const_iterator it = fReactionDataMV.find(molecule) ;
00481
00482 if(it == fReactionDataMV.end())
00483 {
00484 G4cout<<"Nom : " << molecule->GetName()<<G4endl;
00485 G4String errMsg = "No reaction table was implemented for this molecule Definition : "
00486 + molecule -> GetName();
00487 G4Exception("G4MolecularInteractionTable::GetReactionData","",FatalErrorInArgument, errMsg);
00488 return 0;
00489 }
00490
00491 return &(it->second);
00492 }
