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: G4DNAMolecularDecayDisplacer.cc 65022 2012-11-12 16:43:12Z gcosmo $ 00027 // 00028 // Author: Mathieu Karamitros (kara (AT) cenbg . in2p3 . fr) 00029 // 00030 // WARNING : This class is released as a prototype. 00031 // It might strongly evolve or even disapear in the next releases. 00032 // 00033 // History: 00034 // ----------- 00035 // 10 Oct 2011 M.Karamitros created 00036 // 00037 // ------------------------------------------------------------------- 00038 00039 #include "G4DNAMolecularDecayDisplacer.hh" 00040 #include "G4PhysicalConstants.hh" 00041 #include "G4SystemOfUnits.hh" 00042 #include "G4H2O.hh" 00043 #include "G4H2.hh" 00044 #include "G4Hydrogen.hh" 00045 #include "G4OH.hh" 00046 #include "G4H3O.hh" 00047 #include "G4Electron_aq.hh" 00048 #include "G4H2O2.hh" 00049 #include "Randomize.hh" 00050 #include "G4Molecule.hh" 00051 00052 using namespace std; 00053 00054 const DisplacementType G4DNAMolecularDecayDisplacer::Ionisation_DissociationDecay = G4VMolecularDecayDisplacer::AddDisplacement(); 00055 const DisplacementType G4DNAMolecularDecayDisplacer::A1B1_DissociationDecay = G4VMolecularDecayDisplacer::AddDisplacement(); 00056 const DisplacementType G4DNAMolecularDecayDisplacer::B1A1_DissociationDecay = G4VMolecularDecayDisplacer::AddDisplacement(); 00057 const DisplacementType G4DNAMolecularDecayDisplacer::AutoIonisation = G4VMolecularDecayDisplacer::AddDisplacement(); 00058 const DisplacementType G4DNAMolecularDecayDisplacer::DissociativeAttachment = G4VMolecularDecayDisplacer::AddDisplacement(); 00059 00060 G4DNAMolecularDecayDisplacer::G4DNAMolecularDecayDisplacer() : 00061 G4VMolecularDecayDisplacer() 00062 {;} 00063 00064 G4DNAMolecularDecayDisplacer::~G4DNAMolecularDecayDisplacer() 00065 {;} 00066 00067 G4ThreeVector G4DNAMolecularDecayDisplacer::GetMotherMoleculeDisplacement(const G4MolecularDecayChannel* theDecayChannel) const 00068 { 00069 G4int decayType = theDecayChannel -> GetDisplacementType(); 00070 00071 G4double RMSMotherMoleculeDisplacement=0; 00072 00073 if(decayType == Ionisation_DissociationDecay) 00074 { 00075 RMSMotherMoleculeDisplacement = 2.0 * nanometer ; 00076 } 00077 else if(decayType == A1B1_DissociationDecay) 00078 { 00079 RMSMotherMoleculeDisplacement = 0. * nanometer ; 00080 } 00081 else if(decayType == B1A1_DissociationDecay) 00082 { 00083 RMSMotherMoleculeDisplacement = 0. * nanometer ; 00084 } 00085 else if(decayType == AutoIonisation) 00086 { 00087 RMSMotherMoleculeDisplacement = 2.0 * nanometer ; 00088 } 00089 else if(decayType == DissociativeAttachment) 00090 { 00091 RMSMotherMoleculeDisplacement = 0. * nanometer ; 00092 } 00093 00094 if(RMSMotherMoleculeDisplacement==0) 00095 { 00096 return G4ThreeVector(0,0,0); 00097 } 00098 G4ThreeVector RandDirection = radialDistributionOfProducts(RMSMotherMoleculeDisplacement); 00099 00100 return RandDirection; 00101 } 00102 00103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 00104 00105 vector<G4ThreeVector> G4DNAMolecularDecayDisplacer::GetProductsDisplacement(const G4MolecularDecayChannel* theDecayChannel) const 00106 { 00107 G4int nbProducts = theDecayChannel -> GetNbProducts(); 00108 vector<G4ThreeVector> theProductDisplacementVector (nbProducts); 00109 00110 typedef map<const G4MoleculeDefinition*,G4double> RMSmap ; 00111 RMSmap theRMSmap; 00112 00113 G4int decayType = theDecayChannel -> GetDisplacementType(); 00114 00115 if(decayType == Ionisation_DissociationDecay) 00116 { 00117 if(fVerbose) 00118 G4cout<<"Ionisation_DissociationDecay"<<G4endl; 00119 G4double RdmValue = G4UniformRand(); 00120 00121 if(RdmValue< 0.5) 00122 { 00123 // H3O 00124 theRMSmap[G4H3O::Definition()] = 0.* nanometer; 00125 // OH 00126 theRMSmap[G4OH::Definition()] = 0.8* nanometer; 00127 } 00128 else 00129 { 00130 // H3O 00131 theRMSmap[G4H3O::Definition()] = 0.8* nanometer; 00132 // OH 00133 theRMSmap[G4OH::Definition()] = 0.* nanometer; 00134 } 00135 00136 for(int i = 0; i < nbProducts ; i++) 00137 { 00138 G4double theRMSDisplacement; 00139 const G4Molecule* product = theDecayChannel->GetProduct(i); 00140 theRMSDisplacement = theRMSmap[product->GetDefinition()]; 00141 00142 if(theRMSDisplacement==0) 00143 { 00144 theProductDisplacementVector[i] = G4ThreeVector(); 00145 } 00146 else 00147 { 00148 G4ThreeVector RandDirection = radialDistributionOfProducts(theRMSDisplacement); 00149 theProductDisplacementVector[i] = RandDirection; 00150 } 00151 } 00152 } 00153 else if(decayType == A1B1_DissociationDecay) 00154 { 00155 if(fVerbose) 00156 G4cout<<"A1B1_DissociationDecay"<<G4endl; 00157 G4double theRMSDisplacement = 2.4 * nanometer; 00158 G4ThreeVector RandDirection = radialDistributionOfProducts(theRMSDisplacement); 00159 00160 for(G4int i =0 ; i < nbProducts ; i++) 00161 { 00162 const G4Molecule* product = theDecayChannel->GetProduct(i); 00163 if(product->GetDefinition()== G4OH::Definition()) 00164 { 00165 theProductDisplacementVector[i] = -1./18.*RandDirection; 00166 } 00167 else if(product->GetDefinition() == G4Hydrogen::Definition()) 00168 { 00169 theProductDisplacementVector[i] = +17./18.*RandDirection; 00170 } 00171 } 00172 } 00173 else if(decayType == B1A1_DissociationDecay) 00174 { 00175 if(fVerbose) 00176 G4cout<<"B1A1_DissociationDecay"<<G4endl; 00177 G4double theRMSDisplacement = 0.8 * nanometer; 00178 G4ThreeVector RandDirection = radialDistributionOfProducts(theRMSDisplacement); 00179 00180 G4int NbOfOH = 0; 00181 for(G4int i =0 ; i < nbProducts ; i++) 00182 { 00183 const G4Molecule* product = theDecayChannel->GetProduct(i); 00184 if(product->GetDefinition() == G4H2::Definition()) 00185 { 00186 theProductDisplacementVector[i] = -2./18.*RandDirection; 00187 } 00188 else if(product->GetDefinition() == G4OH::Definition()) 00189 { 00190 G4ThreeVector OxygenDisplacement = +16./18.*RandDirection; 00191 G4double OHRMSDisplacement = 1.1 * nanometer; 00192 00193 G4ThreeVector OHDisplacement = radialDistributionOfProducts(OHRMSDisplacement) ; 00194 00195 if(NbOfOH==0) 00196 { 00197 OHDisplacement = 1./2.*OHDisplacement; 00198 } 00199 else 00200 { 00201 OHDisplacement = -1./2.*OHDisplacement; 00202 } 00203 00204 theProductDisplacementVector[i] = OHDisplacement + OxygenDisplacement; 00205 00206 NbOfOH ++; 00207 } 00208 } 00209 } 00210 else if(decayType == AutoIonisation) 00211 { 00212 if(fVerbose) 00213 G4cout<<"AutoIonisation"<<G4endl; 00214 G4double RdmValue = G4UniformRand(); 00215 00216 if(RdmValue< 0.5) 00217 { 00218 // H3O 00219 theRMSmap[G4H3O::Definition()] = 0.* nanometer; 00220 // OH 00221 theRMSmap[G4OH::Definition()] = 0.8* nanometer; 00222 } 00223 else 00224 { 00225 // H3O 00226 theRMSmap[G4H3O::Definition()] = 0.8* nanometer; 00227 // OH 00228 theRMSmap[G4OH::Definition()] = 0.* nanometer; 00229 } 00230 00231 for(G4int i =0 ; i < nbProducts ; i++) 00232 { 00233 G4double theRMSDisplacement; 00234 const G4Molecule* product = theDecayChannel->GetProduct(i); 00235 theRMSDisplacement = theRMSmap[product->GetDefinition()]; 00236 00237 if(theRMSDisplacement==0) 00238 { 00239 theProductDisplacementVector[i] = G4ThreeVector(); 00240 } 00241 else 00242 { 00243 G4ThreeVector RandDirection = radialDistributionOfProducts(theRMSDisplacement); 00244 theProductDisplacementVector[i] = RandDirection; 00245 } 00246 if(product->GetDefinition() == G4Electron_aq::Definition()) 00247 { 00248 theProductDisplacementVector[i]=radialDistributionOfElectron(); 00249 } 00250 } 00251 } 00252 else if(decayType == DissociativeAttachment) 00253 { 00254 if(fVerbose) 00255 G4cout<<"DissociativeAttachment"<<G4endl; 00256 G4double theRMSDisplacement = 0.8 * nanometer; 00257 G4ThreeVector RandDirection = radialDistributionOfProducts(theRMSDisplacement); 00258 00259 G4int NbOfOH = 0; 00260 for(G4int i =0 ; i < nbProducts ; i++) 00261 { 00262 const G4Molecule* product = theDecayChannel->GetProduct(i); 00263 if(product->GetDefinition() == G4H2::Definition()) 00264 { 00265 theProductDisplacementVector[i] = -2./18.*RandDirection; 00266 } 00267 else if(product->GetDefinition() == G4OH::Definition()) 00268 { 00269 G4ThreeVector OxygenDisplacement = +16./18.*RandDirection; 00270 G4double OHRMSDisplacement = 1.1 * nanometer; 00271 00272 G4ThreeVector OHDisplacement = radialDistributionOfProducts(OHRMSDisplacement) ; 00273 00274 if(NbOfOH==0) 00275 { 00276 OHDisplacement = 1./2.*OHDisplacement; 00277 } 00278 else 00279 { 00280 OHDisplacement = -1./2.*OHDisplacement; 00281 } 00282 00283 theProductDisplacementVector[i] = OHDisplacement + OxygenDisplacement; 00284 00285 NbOfOH ++; 00286 } 00287 } 00288 } 00289 00290 return theProductDisplacementVector; 00291 } 00292 00293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 00294 00295 G4ThreeVector G4DNAMolecularDecayDisplacer::radialDistributionOfProducts(G4double Rrms) const 00296 { 00297 G4double sigma = Rrms/sqrt(3.); 00298 G4double expectationValue = 2.*sqrt(2./3.14)*sigma; 00299 00300 G4double XValueForfMax = sqrt(2.*sigma*sigma); 00301 G4double fMaxValue = sqrt(2./3.14) * 1./(sigma*sigma*sigma) * 00302 (XValueForfMax*XValueForfMax)* 00303 exp(-1./2. * (XValueForfMax*XValueForfMax)/(sigma*sigma)); 00304 00305 G4double R(-1.); 00306 00307 do 00308 { 00309 G4double aRandomfValue = fMaxValue*G4UniformRand(); 00310 00311 G4double sign; 00312 if(G4UniformRand() > 0.5) 00313 { 00314 sign = +1.; 00315 } 00316 else 00317 { 00318 sign = -1; 00319 } 00320 00321 R = expectationValue + sign*3.*sigma* G4UniformRand(); 00322 G4double f = sqrt(2./3.14) * 1/pow(sigma, 3) * R*R * exp(-1./2. * R*R/(sigma*sigma)); 00323 00324 if(aRandomfValue < f) 00325 { 00326 break; 00327 } 00328 } 00329 while(1); 00330 00331 G4double costheta = (2.*G4UniformRand()-1.); 00332 G4double theta = acos (costheta); 00333 G4double phi = 2.*pi*G4UniformRand(); 00334 00335 G4double xDirection = R*cos(phi)* sin(theta); 00336 G4double yDirection = R*sin(theta)*sin(phi); 00337 G4double zDirection = R*costheta; 00338 G4ThreeVector RandDirection(xDirection, yDirection, zDirection); 00339 00340 return RandDirection; 00341 } 00342 00343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 00344 00345 G4ThreeVector G4DNAMolecularDecayDisplacer::radialDistributionOfElectron() const 00346 { 00347 00348 G4double sigma = 1./2.; 00349 G4double expectationValue = 1. ; 00350 00351 G4double XValueForfMax = 1./2.; 00352 G4double fMaxValue = 4. * XValueForfMax * 00353 exp(-2. * XValueForfMax); 00354 00355 G4double R(-1.); 00356 00357 do 00358 { 00359 G4double aRandomfValue = fMaxValue*G4UniformRand(); 00360 00361 G4double sign; 00362 if(G4UniformRand() > 0.5) 00363 { 00364 sign = +1; 00365 } 00366 else 00367 { 00368 sign = -1; 00369 } 00370 00371 R = (expectationValue * G4UniformRand() )+ sign*3*sigma* G4UniformRand(); 00372 G4double f = 4* R * exp(- 2. * R); 00373 00374 if(aRandomfValue < f) 00375 { 00376 break; 00377 } 00378 } 00379 while(1); 00380 00381 G4double Rnano = R *10* nanometer; 00382 00383 G4double costheta = (2*G4UniformRand()-1); 00384 G4double theta = acos (costheta); 00385 G4double phi = 2*pi*G4UniformRand(); 00386 00387 G4double xDirection = Rnano*cos(phi)* sin(theta); 00388 G4double yDirection = Rnano*sin(theta)*sin(phi); 00389 G4double zDirection = Rnano*costheta; 00390 G4ThreeVector RandDirection(xDirection, yDirection, zDirection); 00391 00392 return RandDirection; 00393 }