Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAMolecularDecayDisplacer.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4DNAMolecularDecayDisplacer.cc 65022 2012-11-12 16:43:12Z gcosmo $
27 //
28 // Author: Mathieu Karamitros (kara (AT) cenbg . in2p3 . fr)
29 //
30 // WARNING : This class is released as a prototype.
31 // It might strongly evolve or even disapear in the next releases.
32 //
33 // History:
34 // -----------
35 // 10 Oct 2011 M.Karamitros created
36 //
37 // -------------------------------------------------------------------
38 
39 #include "G4DNAMolecularDecayDisplacer.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4H2O.hh"
43 #include "G4H2.hh"
44 #include "G4Hydrogen.hh"
45 #include "G4OH.hh"
46 #include "G4H3O.hh"
47 #include "G4Electron_aq.hh"
48 #include "G4H2O2.hh"
49 #include "Randomize.hh"
50 #include "G4Molecule.hh"
51 
52 using namespace std;
53 
54 const DisplacementType G4DNAMolecularDecayDisplacer::Ionisation_DissociationDecay = G4VMolecularDecayDisplacer::AddDisplacement();
55 const DisplacementType G4DNAMolecularDecayDisplacer::A1B1_DissociationDecay = G4VMolecularDecayDisplacer::AddDisplacement();
56 const DisplacementType G4DNAMolecularDecayDisplacer::B1A1_DissociationDecay = G4VMolecularDecayDisplacer::AddDisplacement();
57 const DisplacementType G4DNAMolecularDecayDisplacer::AutoIonisation = G4VMolecularDecayDisplacer::AddDisplacement();
58 const DisplacementType G4DNAMolecularDecayDisplacer::DissociativeAttachment = G4VMolecularDecayDisplacer::AddDisplacement();
59 
60 G4DNAMolecularDecayDisplacer::G4DNAMolecularDecayDisplacer() :
61  G4VMolecularDecayDisplacer()
62 {;}
63 
64 G4DNAMolecularDecayDisplacer::~G4DNAMolecularDecayDisplacer()
65 {;}
66 
67 G4ThreeVector G4DNAMolecularDecayDisplacer::GetMotherMoleculeDisplacement(const G4MolecularDecayChannel* theDecayChannel) const
68 {
69  G4int decayType = theDecayChannel -> GetDisplacementType();
70 
71  G4double RMSMotherMoleculeDisplacement=0;
72 
73  if(decayType == Ionisation_DissociationDecay)
74  {
75  RMSMotherMoleculeDisplacement = 2.0 * nanometer ;
76  }
77  else if(decayType == A1B1_DissociationDecay)
78  {
79  RMSMotherMoleculeDisplacement = 0. * nanometer ;
80  }
81  else if(decayType == B1A1_DissociationDecay)
82  {
83  RMSMotherMoleculeDisplacement = 0. * nanometer ;
84  }
85  else if(decayType == AutoIonisation)
86  {
87  RMSMotherMoleculeDisplacement = 2.0 * nanometer ;
88  }
89  else if(decayType == DissociativeAttachment)
90  {
91  RMSMotherMoleculeDisplacement = 0. * nanometer ;
92  }
93 
94  if(RMSMotherMoleculeDisplacement==0)
95  {
96  return G4ThreeVector(0,0,0);
97  }
98  G4ThreeVector RandDirection = radialDistributionOfProducts(RMSMotherMoleculeDisplacement);
99 
100  return RandDirection;
101 }
102 
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
104 
105 vector<G4ThreeVector> G4DNAMolecularDecayDisplacer::GetProductsDisplacement(const G4MolecularDecayChannel* theDecayChannel) const
106 {
107  G4int nbProducts = theDecayChannel -> GetNbProducts();
108  vector<G4ThreeVector> theProductDisplacementVector (nbProducts);
109 
110  typedef map<const G4MoleculeDefinition*,G4double> RMSmap ;
111  RMSmap theRMSmap;
112 
113  G4int decayType = theDecayChannel -> GetDisplacementType();
114 
115  if(decayType == Ionisation_DissociationDecay)
116  {
117  if(fVerbose)
118  G4cout<<"Ionisation_DissociationDecay"<<G4endl;
119  G4double RdmValue = G4UniformRand();
120 
121  if(RdmValue< 0.5)
122  {
123  // H3O
124  theRMSmap[G4H3O::Definition()] = 0.* nanometer;
125  // OH
126  theRMSmap[G4OH::Definition()] = 0.8* nanometer;
127  }
128  else
129  {
130  // H3O
131  theRMSmap[G4H3O::Definition()] = 0.8* nanometer;
132  // OH
133  theRMSmap[G4OH::Definition()] = 0.* nanometer;
134  }
135 
136  for(int i = 0; i < nbProducts ; i++)
137  {
138  G4double theRMSDisplacement;
139  const G4Molecule* product = theDecayChannel->GetProduct(i);
140  theRMSDisplacement = theRMSmap[product->GetDefinition()];
141 
142  if(theRMSDisplacement==0)
143  {
144  theProductDisplacementVector[i] = G4ThreeVector();
145  }
146  else
147  {
148  G4ThreeVector RandDirection = radialDistributionOfProducts(theRMSDisplacement);
149  theProductDisplacementVector[i] = RandDirection;
150  }
151  }
152  }
153  else if(decayType == A1B1_DissociationDecay)
154  {
155  if(fVerbose)
156  G4cout<<"A1B1_DissociationDecay"<<G4endl;
157  G4double theRMSDisplacement = 2.4 * nanometer;
158  G4ThreeVector RandDirection = radialDistributionOfProducts(theRMSDisplacement);
159 
160  for(G4int i =0 ; i < nbProducts ; i++)
161  {
162  const G4Molecule* product = theDecayChannel->GetProduct(i);
163  if(product->GetDefinition()== G4OH::Definition())
164  {
165  theProductDisplacementVector[i] = -1./18.*RandDirection;
166  }
167  else if(product->GetDefinition() == G4Hydrogen::Definition())
168  {
169  theProductDisplacementVector[i] = +17./18.*RandDirection;
170  }
171  }
172  }
173  else if(decayType == B1A1_DissociationDecay)
174  {
175  if(fVerbose)
176  G4cout<<"B1A1_DissociationDecay"<<G4endl;
177  G4double theRMSDisplacement = 0.8 * nanometer;
178  G4ThreeVector RandDirection = radialDistributionOfProducts(theRMSDisplacement);
179 
180  G4int NbOfOH = 0;
181  for(G4int i =0 ; i < nbProducts ; i++)
182  {
183  const G4Molecule* product = theDecayChannel->GetProduct(i);
184  if(product->GetDefinition() == G4H2::Definition())
185  {
186  theProductDisplacementVector[i] = -2./18.*RandDirection;
187  }
188  else if(product->GetDefinition() == G4OH::Definition())
189  {
190  G4ThreeVector OxygenDisplacement = +16./18.*RandDirection;
191  G4double OHRMSDisplacement = 1.1 * nanometer;
192 
193  G4ThreeVector OHDisplacement = radialDistributionOfProducts(OHRMSDisplacement) ;
194 
195  if(NbOfOH==0)
196  {
197  OHDisplacement = 1./2.*OHDisplacement;
198  }
199  else
200  {
201  OHDisplacement = -1./2.*OHDisplacement;
202  }
203 
204  theProductDisplacementVector[i] = OHDisplacement + OxygenDisplacement;
205 
206  NbOfOH ++;
207  }
208  }
209  }
210  else if(decayType == AutoIonisation)
211  {
212  if(fVerbose)
213  G4cout<<"AutoIonisation"<<G4endl;
214  G4double RdmValue = G4UniformRand();
215 
216  if(RdmValue< 0.5)
217  {
218  // H3O
219  theRMSmap[G4H3O::Definition()] = 0.* nanometer;
220  // OH
221  theRMSmap[G4OH::Definition()] = 0.8* nanometer;
222  }
223  else
224  {
225  // H3O
226  theRMSmap[G4H3O::Definition()] = 0.8* nanometer;
227  // OH
228  theRMSmap[G4OH::Definition()] = 0.* nanometer;
229  }
230 
231  for(G4int i =0 ; i < nbProducts ; i++)
232  {
233  G4double theRMSDisplacement;
234  const G4Molecule* product = theDecayChannel->GetProduct(i);
235  theRMSDisplacement = theRMSmap[product->GetDefinition()];
236 
237  if(theRMSDisplacement==0)
238  {
239  theProductDisplacementVector[i] = G4ThreeVector();
240  }
241  else
242  {
243  G4ThreeVector RandDirection = radialDistributionOfProducts(theRMSDisplacement);
244  theProductDisplacementVector[i] = RandDirection;
245  }
246  if(product->GetDefinition() == G4Electron_aq::Definition())
247  {
248  theProductDisplacementVector[i]=radialDistributionOfElectron();
249  }
250  }
251  }
252  else if(decayType == DissociativeAttachment)
253  {
254  if(fVerbose)
255  G4cout<<"DissociativeAttachment"<<G4endl;
256  G4double theRMSDisplacement = 0.8 * nanometer;
257  G4ThreeVector RandDirection = radialDistributionOfProducts(theRMSDisplacement);
258 
259  G4int NbOfOH = 0;
260  for(G4int i =0 ; i < nbProducts ; i++)
261  {
262  const G4Molecule* product = theDecayChannel->GetProduct(i);
263  if(product->GetDefinition() == G4H2::Definition())
264  {
265  theProductDisplacementVector[i] = -2./18.*RandDirection;
266  }
267  else if(product->GetDefinition() == G4OH::Definition())
268  {
269  G4ThreeVector OxygenDisplacement = +16./18.*RandDirection;
270  G4double OHRMSDisplacement = 1.1 * nanometer;
271 
272  G4ThreeVector OHDisplacement = radialDistributionOfProducts(OHRMSDisplacement) ;
273 
274  if(NbOfOH==0)
275  {
276  OHDisplacement = 1./2.*OHDisplacement;
277  }
278  else
279  {
280  OHDisplacement = -1./2.*OHDisplacement;
281  }
282 
283  theProductDisplacementVector[i] = OHDisplacement + OxygenDisplacement;
284 
285  NbOfOH ++;
286  }
287  }
288  }
289 
290  return theProductDisplacementVector;
291 }
292 
293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
294 
295 G4ThreeVector G4DNAMolecularDecayDisplacer::radialDistributionOfProducts(G4double Rrms) const
296 {
297  G4double sigma = Rrms/sqrt(3.);
298  G4double expectationValue = 2.*sqrt(2./3.14)*sigma;
299 
300  G4double XValueForfMax = sqrt(2.*sigma*sigma);
301  G4double fMaxValue = sqrt(2./3.14) * 1./(sigma*sigma*sigma) *
302  (XValueForfMax*XValueForfMax)*
303  exp(-1./2. * (XValueForfMax*XValueForfMax)/(sigma*sigma));
304 
305  G4double R(-1.);
306 
307  do
308  {
309  G4double aRandomfValue = fMaxValue*G4UniformRand();
310 
311  G4double sign;
312  if(G4UniformRand() > 0.5)
313  {
314  sign = +1.;
315  }
316  else
317  {
318  sign = -1;
319  }
320 
321  R = expectationValue + sign*3.*sigma* G4UniformRand();
322  G4double f = sqrt(2./3.14) * 1/pow(sigma, 3) * R*R * exp(-1./2. * R*R/(sigma*sigma));
323 
324  if(aRandomfValue < f)
325  {
326  break;
327  }
328  }
329  while(1);
330 
331  G4double costheta = (2.*G4UniformRand()-1.);
332  G4double theta = acos (costheta);
333  G4double phi = 2.*pi*G4UniformRand();
334 
335  G4double xDirection = R*cos(phi)* sin(theta);
336  G4double yDirection = R*sin(theta)*sin(phi);
337  G4double zDirection = R*costheta;
338  G4ThreeVector RandDirection(xDirection, yDirection, zDirection);
339 
340  return RandDirection;
341 }
342 
343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
344 
345 G4ThreeVector G4DNAMolecularDecayDisplacer::radialDistributionOfElectron() const
346 {
347 
348  G4double sigma = 1./2.;
349  G4double expectationValue = 1. ;
350 
351  G4double XValueForfMax = 1./2.;
352  G4double fMaxValue = 4. * XValueForfMax *
353  exp(-2. * XValueForfMax);
354 
355  G4double R(-1.);
356 
357  do
358  {
359  G4double aRandomfValue = fMaxValue*G4UniformRand();
360 
361  G4double sign;
362  if(G4UniformRand() > 0.5)
363  {
364  sign = +1;
365  }
366  else
367  {
368  sign = -1;
369  }
370 
371  R = (expectationValue * G4UniformRand() )+ sign*3*sigma* G4UniformRand();
372  G4double f = 4* R * exp(- 2. * R);
373 
374  if(aRandomfValue < f)
375  {
376  break;
377  }
378  }
379  while(1);
380 
381  G4double Rnano = R *10* nanometer;
382 
383  G4double costheta = (2*G4UniformRand()-1);
384  G4double theta = acos (costheta);
385  G4double phi = 2*pi*G4UniformRand();
386 
387  G4double xDirection = Rnano*cos(phi)* sin(theta);
388  G4double yDirection = Rnano*sin(theta)*sin(phi);
389  G4double zDirection = Rnano*costheta;
390  G4ThreeVector RandDirection(xDirection, yDirection, zDirection);
391 
392  return RandDirection;
393 }
Generated on Sat May 25 2013 14:33:31 for Geant4 by   doxygen 1.8.4