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 "G4DNAChemistryManager.hh"
00040 #include "G4SystemOfUnits.hh"
00041 #include "G4Molecule.hh"
00042 #include "G4ITTrackHolder.hh"
00043 #include "G4H2O.hh"
00044 #include "G4DNAMolecularReactionTable.hh"
00045 #include "G4DNAWaterExcitationStructure.hh"
00046 #include "G4DNAWaterIonisationStructure.hh"
00047 #include "G4Electron_aq.hh"
00048 #include "G4ITManager.hh"
00049 #include "G4MolecularConfiguration.hh"
00050 #include "G4MoleculeCounter.hh"
00051
00052 using namespace std;
00053
00054 auto_ptr<G4DNAChemistryManager> G4DNAChemistryManager::fInstance(0);
00055
00056 G4DNAChemistryManager::G4DNAChemistryManager() :
00057 fActiveChemistry(false)
00058 {
00059 fExcitationLevel = 0;
00060 fIonisationLevel = 0;
00061 fWriteFile = false;
00062 }
00063
00064 G4DNAChemistryManager* G4DNAChemistryManager::Instance()
00065 {
00066 if(!fInstance.get()) fInstance = auto_ptr<G4DNAChemistryManager>(new G4DNAChemistryManager());
00067 return fInstance.get();
00068 }
00069
00070 G4DNAChemistryManager::~G4DNAChemistryManager()
00071 {
00072 if (fOutput.is_open())
00073 {
00074 fOutput.close();
00075 }
00076 if(fIonisationLevel) delete fIonisationLevel;
00077 if(fExcitationLevel) delete fExcitationLevel;
00078 G4DNAMolecularReactionTable::DeleteInstance();
00079 G4MoleculeHandleManager::DeleteInstance();
00080 G4MolecularConfiguration::DeleteManager();
00081 fInstance.release();
00082 G4MoleculeCounter::DeleteInstance();
00083 }
00084
00085 void G4DNAChemistryManager::DeleteInstance()
00086 {
00087 if(fInstance.get())
00088 fInstance.reset();
00089 }
00090
00091 void G4DNAChemistryManager::WriteInto(const G4String& output,
00092 ios_base::openmode mode)
00093 {
00094 fOutput.open(output.data(), mode);
00095 fOutput << std::setprecision(6) << std::scientific;
00096
00097 fOutput << setw(11) << left << "#Parent ID"
00098 << setw(10) << "Molecule"
00099 << setw(14) << "Elec Modif"
00100 << setw(13) << "Energy (eV)"
00101 << setw(22) << "X pos of parent [nm]"
00102 << setw(22) << "Y pos of parent [nm]"
00103 << setw(22) << "Z pos of parent [nm]"
00104 << setw(14) << "X pos [nm]"
00105 << setw(14) << "Y pos [nm]"
00106 << setw(14) << "Z pos [nm]"
00107 << G4endl
00108 << setw(21) << "#"
00109 << setw(13) << "1)io/ex=0/1"
00110 << G4endl
00111 << setw(21) << "#"
00112 << setw(13) << "2)level=0...5"
00113 << G4endl;
00114 fWriteFile = true;
00115 }
00116
00117 void G4DNAChemistryManager::CloseFile()
00118 {
00119 if (fOutput.is_open())
00120 {
00121 fOutput.close();
00122 }
00123 fWriteFile = false;
00124 }
00125
00126 G4DNAWaterExcitationStructure* G4DNAChemistryManager::GetExcitationLevel()
00127 {
00128 if(!fExcitationLevel)
00129 {
00130 fExcitationLevel = new G4DNAWaterExcitationStructure;
00131 }
00132 return fExcitationLevel;
00133 }
00134
00135 G4DNAWaterIonisationStructure* G4DNAChemistryManager::GetIonisationLevel()
00136 {
00137 if(!fIonisationLevel)
00138 {
00139 fIonisationLevel = new G4DNAWaterIonisationStructure;
00140 }
00141 return fIonisationLevel;
00142 }
00143
00144 void G4DNAChemistryManager::CreateWaterMolecule(ElectronicModification modification,
00145 G4int electronicLevel,
00146 const G4Track* theIncomingTrack)
00147 {
00148 if(fWriteFile)
00149 {
00150 G4double energy = -1.;
00151
00152 switch (modification)
00153 {
00154 case eDissociativeAttachment:
00155 energy = -1;
00156 break;
00157 case eExcitedMolecule :
00158 energy = GetExcitationLevel()->ExcitationEnergy(electronicLevel);
00159 break;
00160 case eIonizedMolecule :
00161 energy = GetIonisationLevel()->IonisationEnergy(electronicLevel);
00162 break;
00163 }
00164
00165 fOutput << setw(11) << left << theIncomingTrack->GetTrackID()
00166 << setw(10) << "H2O"
00167 << left << modification
00168 << internal <<":"
00169 << right <<electronicLevel
00170 << left
00171 << setw(11) << ""
00172 << std::setprecision(2) << std::fixed
00173 << setw(13) << energy/eV
00174 << std::setprecision(6) << std::scientific
00175 << setw(22) << (theIncomingTrack->GetPosition().x())/nanometer
00176 << setw(22) << (theIncomingTrack->GetPosition().y())/nanometer
00177 << setw(22) << (theIncomingTrack->GetPosition().z())/nanometer
00178 << G4endl;
00179 }
00180
00181 if(fActiveChemistry)
00182 {
00183 G4Molecule * H2O = new G4Molecule (G4H2O::Definition());
00184
00185 switch (modification)
00186 {
00187 case eDissociativeAttachment:
00188 H2O -> AddElectron(5,1);
00189 break;
00190 case eExcitedMolecule :
00191 H2O -> ExciteMolecule(electronicLevel);
00192 break;
00193 case eIonizedMolecule :
00194 H2O -> IonizeMolecule(electronicLevel);
00195 break;
00196 }
00197
00198 G4Track * H2OTrack = H2O->BuildTrack(1*picosecond,
00199 theIncomingTrack->GetPosition());
00200
00201 H2OTrack -> SetParentID(theIncomingTrack->GetTrackID());
00202 H2OTrack -> SetTrackStatus(fStopButAlive);
00203 H2OTrack -> SetKineticEnergy(0.);
00204
00205 G4ITTrackHolder::Instance()->PushTrack(H2OTrack);
00206 }
00207 }
00208
00209 void G4DNAChemistryManager::CreateSolvatedElectron(const G4Track* theIncomingTrack,
00210 G4ThreeVector* finalPosition)
00211
00212 {
00213 if(fWriteFile)
00214 {
00215 fOutput << setw(11)<< theIncomingTrack->GetTrackID()
00216 << setw(10)<< "e_aq"
00217 << setw(14)<< -1
00218 << std::setprecision(2) << std::fixed
00219 << setw(13)<< theIncomingTrack->GetKineticEnergy()/eV
00220 << std::setprecision(6) << std::scientific
00221 << setw(22)<< (theIncomingTrack->GetPosition().x())/nanometer
00222 << setw(22)<< (theIncomingTrack->GetPosition().y())/nanometer
00223 << setw(22)<< (theIncomingTrack->GetPosition().z())/nanometer ;
00224
00225 if(finalPosition != 0)
00226 {
00227 fOutput<< setw(14)<< (finalPosition->x())/nanometer
00228 << setw(14)<< (finalPosition->y())/nanometer
00229 << setw(14)<< (finalPosition->z())/nanometer ;
00230 }
00231
00232 fOutput << G4endl;
00233 }
00234
00235 if(fActiveChemistry)
00236 {
00237 G4Molecule* e_aq = new G4Molecule(G4Electron_aq::Definition());
00238 G4Track * e_aqTrack(0);
00239 if(finalPosition)
00240 {
00241 e_aqTrack = e_aq->BuildTrack(picosecond,*finalPosition);
00242 }
00243 else
00244 {
00245 e_aqTrack = e_aq->BuildTrack(picosecond,theIncomingTrack->GetPosition());
00246 }
00247 e_aqTrack -> SetTrackStatus(fAlive);
00248 e_aqTrack -> SetParentID(theIncomingTrack->GetTrackID());
00249 G4ITTrackHolder::Instance()->PushTrack(e_aqTrack);
00250 G4ITManager<G4Molecule>::Instance()->Push(e_aqTrack);
00251 }
00252 }
00253
00254
00255 void G4DNAChemistryManager::PushMolecule(G4Molecule*& molecule, double time,
00256 const G4ThreeVector& position, int parentID)
00257 {
00258 if(fWriteFile)
00259 {
00260 fOutput << setw(11)<< parentID
00261 << setw(10)<< molecule->GetName()
00262 << setw(14)<< -1
00263 << std::setprecision(2) << std::fixed
00264 << setw(13)<< -1
00265 << std::setprecision(6) << std::scientific
00266 << setw(22)<< (position.x())/nanometer
00267 << setw(22)<< (position.y())/nanometer
00268 << setw(22)<< (position.z())/nanometer;
00269 fOutput << G4endl;
00270 }
00271
00272 if(fActiveChemistry)
00273 {
00274 G4Track* track = molecule->BuildTrack(time,position);
00275 track -> SetTrackStatus(fAlive);
00276 track -> SetParentID(parentID);
00277 G4ITTrackHolder::Instance()->PushTrack(track);
00278 G4ITManager<G4Molecule>::Instance()->Push(track);
00279 }
00280 else
00281 {
00282 delete molecule;
00283 molecule = 0;
00284 }
00285 }
00286
00287 void G4DNAChemistryManager::PushMoleculeAtParentTimeAndPlace(G4Molecule*& molecule,
00288 const G4Track* theIncomingTrack)
00289 {
00290 if(fWriteFile)
00291 {
00292 fOutput << setw(11)<< theIncomingTrack->GetTrackID()
00293 << setw(10)<< molecule->GetName()
00294 << setw(14)<< -1
00295 << std::setprecision(2) << std::fixed
00296 << setw(13)<< theIncomingTrack->GetKineticEnergy()/eV
00297 << std::setprecision(6) << std::scientific
00298 << setw(22)<< (theIncomingTrack->GetPosition().x())/nanometer
00299 << setw(22)<< (theIncomingTrack->GetPosition().y())/nanometer
00300 << setw(22)<< (theIncomingTrack->GetPosition().z())/nanometer ;
00301 fOutput << G4endl;
00302 }
00303
00304 if(fActiveChemistry)
00305 {
00306 G4Track* track = molecule->BuildTrack(theIncomingTrack->GetGlobalTime(),theIncomingTrack->GetPosition());
00307 track -> SetTrackStatus(fAlive);
00308 track -> SetParentID(theIncomingTrack->GetTrackID());
00309 G4ITTrackHolder::Instance()->PushTrack(track);
00310 G4ITManager<G4Molecule>::Instance()->Push(track);
00311 }
00312 else
00313 {
00314 delete molecule;
00315 molecule = 0;
00316 }
00317 }