Geant4_10
G4DNABrownianTransportation.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: G4DNABrownianTransportation.cc 74551 2013-10-14 12:59:14Z 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 
41 
42 #include <CLHEP/Random/Stat.h>
43 
45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh"
47 #include "Randomize.hh"
48 #include "G4Molecule.hh"
49 #include "G4RandomDirection.hh"
50 #include "G4ParticleTable.hh"
51 #include "G4SafetyHelper.hh"
53 #include "G4UnitsTable.hh"
54 #include "G4NistManager.hh"
56 
57 using namespace std;
58 
59 #ifndef State
60 #define State(theXInfo) (fpBrownianState->theXInfo)
61 #endif
62 
63 //#ifndef State
64 //#define State(theXInfo) (fTransportationState->theXInfo)
65 //#endif
66 
67 //COLOR FOR DEBUGING
68 //#define RED_ON_WHITE "\033[0;31m"
69 //#define GREEN "\033[32;40m"
70 #define GREEN_ON_BLUE "\033[1;32;44m"
71 #define RESET "\033[0m"
72 
73 static double InvErf(double x)
74 {
76 }
77 
78 static double InvErfc(double x)
79 {
80  return CLHEP::HepStat::inverseErf(1.-x);
81 }
82 
84  G4ITTransportation(aName, verbosity),
85  InitProcessState(fpBrownianState, fTransportationState)
86 {
87  //ctor
89  verboseLevel = 1;
92  fpWaterDensity = 0;
93 }
94 
96 {;}
97 
99  G4ITTransportation(right),
100  InitProcessState(fpBrownianState, fTransportationState)
101 {
102  //copy ctor
103  SetProcessSubType(61);
104  verboseLevel = right.verboseLevel;
106  fNistWater = right.fNistWater;
108 }
109 
111 {
112  if (this == &rhs) return *this; // handle self assignment
113  //assignment operator
114  return *this;
115 }
116 
118 {
119  fPathLengthWasCorrected = false;
120 }
121 
123 {
124  fpState = new G4ITBrownianState();
127 }
128 
130 {
131  if(verboseLevel > 0)
132  {
133  G4cout << G4endl << GetProcessName() << ": for "
134  << setw(24) << particle.GetParticleName()
135  << "\tSubType= " << GetProcessSubType() << G4endl;
136  }
137  // Initialize water density pointer
139 }
140 
142  const G4Step& step,
143  const double timeStep,
144  double& spaceStep)
145 {
146  // G4cout << "G4ITBrownianTransportation::ComputeStep" << G4endl;
147 
148  /* If this method is called, this step
149  * cannot be geometry limited.
150  * In case the step is limited by the geometry,
151  * this method should not be called.
152  * The fTransportEndPosition calculated in
153  * the method AlongStepIL should be taken
154  * into account.
155  * In order to do so, the flag IsLeadingStep
156  * is on. Meaning : this track has the minimum
157  * interaction length over all others.
158  */
159  if(GetIT(track)->GetTrackingInfo()->IsLeadingStep())
160  {
161  const G4VITProcess* ITProc = ((const G4VITProcess*) step.GetPostStepPoint()->GetProcessDefinedStep());
162  bool makeException = true;
163 
164  if(ITProc && ITProc->ProposesTimeStep()) makeException = false;
165 
166  if(makeException)
167  {
168 
169  G4ExceptionDescription exceptionDescription ;
170  exceptionDescription << "ComputeStep is called while the track has the minimum interaction time";
171  exceptionDescription << " so it should not recompute a timeStep ";
172  G4Exception("G4DNABrownianTransportation::ComputeStep","G4DNABrownianTransportation001",
173  FatalErrorInArgument,exceptionDescription);
174  }
175  }
176 
177  State(fGeometryLimitedStep) = false;
178  // TODO : generalize this process to all kind of brownian objects
179  // G4ITBrownianObject* ITBrown = GetITBrownianObject(track) ;
180  // G4double diffCoeff = ITBrown->GetDiffusionCoefficient(track.GetMaterial());
181  G4Molecule* molecule = GetMolecule(track) ;
182  G4double diffCoeff = molecule->GetDiffusionCoefficient();
183 
184  if(timeStep > 0)
185  {
186  spaceStep = DBL_MAX;
187 
188  while(spaceStep > State(endpointDistance))
189  // Probably inefficient when the track is close close to boundaries
190  // it goes with fUserMaximumTimeBeforeReachingBoundary == false
191  // fUserMaximumTimeBeforeReachingBoundary == true, it should never loop
192  {
193  G4double x = G4RandGauss::shoot(0,sqrt(2*diffCoeff*timeStep));
194  G4double y = G4RandGauss::shoot(0,sqrt(2*diffCoeff*timeStep));
195  G4double z = G4RandGauss::shoot(0,sqrt(2*diffCoeff*timeStep));
196 
197  spaceStep = sqrt(x*x + y*y + z*z);
198  }
199  // State(fTransportEndPosition).set(x + track.GetPosition().x(),
200  // y + track.GetPosition().y(),
201  // z + track.GetPosition().z());
202 
203  State(fTransportEndPosition)= spaceStep*step.GetPostStepPoint()->GetMomentumDirection() + track.GetPosition();
204  }
205  else
206  {
207  spaceStep = 0. ;
208  State(fTransportEndPosition) = track.GetPosition() ;
209  }
210 
211  State(fCandidateEndGlobalTime) = step.GetPreStepPoint()->GetGlobalTime() + timeStep ;
212  State(fEndGlobalTimeComputed) = true ;
213 
214 #ifdef G4VERBOSE
215  // DEBUG
216  if(fVerboseLevel>1)
217  {
219  << "G4ITBrownianTransportation::ComputeStep() : "
220  << " trackID : " << track.GetTrackID()
221  << " : Molecule name: " << molecule-> GetName()
222  << G4endl
223  << "Diffusion length : " << G4BestUnit(spaceStep, "Length")
224  << " within time step : " << G4BestUnit(timeStep,"Time")
225  << RESET
226  << G4endl<< G4endl;
227  }
228 #endif
229 }
230 
232 {
234 
235 #ifdef G4VERBOSE
236  // DEBUG
237  if(fVerboseLevel>1)
238  {
240  << "G4ITBrownianTransportation::PostStepDoIt() :"
241  << " trackID : " << track.GetTrackID()
242  << " Molecule name: " << GetMolecule(track)-> GetName() << G4endl;
243  G4cout<< "Diffusion length : "<< G4BestUnit(step.GetStepLength(),"Length") <<" within time step : "
244  << G4BestUnit(step.GetDeltaTime(),"Time") << "\t"
245  << " Current global time : " << G4BestUnit(track.GetGlobalTime(),"Time")
246  << RESET
247  << G4endl<< G4endl;
248  }
249 #endif
250  return &fParticleChange ;
251 }
252 
254  const G4Track& track)
255 {
256 
257 #ifdef G4VERBOSE
258  // DEBUG
259  if (fVerboseLevel>1)
260  {
262  << setw(18)<< "G4DNABrownianTransportation::Diffusion :"
263  << setw(8) << GetIT(track)->GetName()
264  << "\t trackID:" << track.GetTrackID() <<"\t"
265  << " Global Time = " << G4BestUnit(track.GetGlobalTime(),"Time")
266  << RESET
267  << G4endl<< G4endl;
268  }
269 #endif
270 
271  G4Material* material = track.GetMaterial();
272 // if (material != fNistWater && material->GetBaseMaterial() != fNistWater)
273 
274  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
275 
276  if(waterDensity == 0.0)
277 // if (material == nistwater || material->GetBaseMaterial() == nistwater)
278  {
279  G4cout << "A track is outside water material : trackID"<< track.GetTrackID() << " (" << GetMolecule(track)->GetName() <<")" << G4endl;
280  G4cout << "Local Time : "<< (track.GetLocalTime()) /s<<G4endl;
281  G4cout<< "Step Number :" << track.GetCurrentStepNumber() <<G4endl;
282 
285 
286  // Got pb with :
287  // fParticleChange.ProposeTrackStatus(fStopAndKill);
288  // It makes the tracks going straight without killing them
289 
290  return ; // &fParticleChange is the final returned object
291  }
292 
293  G4double costheta = (2*G4UniformRand()-1);
294  G4double theta = acos (costheta);
295  G4double phi = 2*pi*G4UniformRand();
296 
297  G4double xMomentum = cos(phi)* sin(theta);
298  G4double yMomentum = sin(theta)*sin(phi);
299  G4double zMomentum = costheta;
300 
301  fParticleChange.ProposeMomentumDirection(xMomentum, yMomentum, zMomentum);
302  State(fMomentumChanged) = true;
304 
305  // G4cout << "BROWN : Propose new direction :" << G4ThreeVector(xMomentum, yMomentum, zMomentum) << G4endl;
306 
307  // Alternative
308  //fParticleChange.ProposeMomentumDirection(G4RandomDirection());
309 
310  return; // &fParticleChange is the final returned object
311 }
312 
313 
315  const G4Track& track,
316  G4double previousStepSize,
317  G4double currentMinimumStep,
318  G4double& currentSafety,
319  G4GPILSelection* selection)
320 {
322  previousStepSize,
323  currentMinimumStep,
324  currentSafety,
325  selection);
326 
327  G4double diffusionCoefficient = GetMolecule(track)->GetDiffusionCoefficient();
328  // G4double diffusionCoefficient = GetITBrownianObject(track)->GetDiffusionCoefficient(track.GetMaterial());
329 
330  if(State(fGeometryLimitedStep))
331  {
332  // 99 % of the space step distribution is lower than
333  // d_99 = 8 * sqrt(D*t)
334  // where t is the corresponding time step
335  // so by inversion :
337  {
338  State(theInteractionTimeLeft) = (geometryStepLength*geometryStepLength)/(64 * diffusionCoefficient);
339  }
340  else // Will use a random time
341  {
342  State(theInteractionTimeLeft) = 1/(4*diffusionCoefficient) * pow(geometryStepLength/InvErfc(G4UniformRand()),2);
343  }
344 
345  State(fCandidateEndGlobalTime) = track.GetGlobalTime() + State(theInteractionTimeLeft);
346  State(fPathLengthWasCorrected) = false;
347  }
348  else
349  {
350  geometryStepLength = 2*sqrt(diffusionCoefficient*State(theInteractionTimeLeft) ) *InvErf(G4UniformRand());
351  State(fPathLengthWasCorrected) = true;
352  State(endpointDistance) = geometryStepLength;
353  }
354 
355  return geometryStepLength ;
356 }
357 
359 //
360 // Initialize ParticleChange (by setting all its members equal
361 // to corresponding members in G4Track)
363  const G4Step& step )
364 {
366  Diffusion(track);
367  return &fParticleChange;
368 }
ThreeVector shoot(const G4int Ap, const G4int Af)
const std::vector< double > * GetDensityTableFor(const G4Material *) const
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
G4double GetLocalTime() const
G4ParticleChangeForTransport fParticleChange
G4DNABrownianTransportation & operator=(const G4DNABrownianTransportation &other)
G4int verboseLevel
Definition: G4VProcess.hh:368
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetStepLength() const
const XML_Char * s
Definition: expat.h:262
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:578
virtual G4VParticleChange * AlongStepDoIt(const G4Track &track, const G4Step &)
size_t GetIndex() const
Definition: G4Material.hh:260
static double inverseErf(double t)
const G4ThreeVector & GetPosition() const
virtual G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &)
G4double GetDiffusionCoefficient() const
Definition: G4Molecule.cc:389
virtual void StartTracking(G4Track *aTrack)
virtual G4VParticleChange * AlongStepDoIt(const G4Track &track, const G4Step &stepData)
G4ProcessState * fpState
tuple x
Definition: test.py:50
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
virtual const G4String & GetName() const =0
virtual G4double AlongStepGetPhysicalInteractionLength(const G4Track &track, G4double, G4double currentMinimumStep, G4double &currentSafety, G4GPILSelection *selection)
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
const G4String & GetParticleName() const
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
Double_t y
Definition: plot.C:279
string material
Definition: eplot.py:19
G4StepPoint * GetPreStepPoint() const
#define State(theXInfo)
const G4String & GetName() const
Definition: G4Molecule.cc:243
virtual void StartTracking(G4Track *aTrack)
const G4ThreeVector & GetMomentumDirection() const
G4IT * GetIT(const G4Track *track)
Definition: G4IT.cc:48
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
G4int GetCurrentStepNumber() const
void SetInstantiateProcessState(G4bool flag)
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
Definition: G4Step.hh:76
G4double GetDeltaTime() const
G4int GetTrackID() const
G4double GetGlobalTime() const
G4Molecule * GetMolecule(const G4Track &track)
Definition: G4Molecule.cc:67
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
G4Material * GetMaterial() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4DNAMolecularMaterial * Instance()
const G4VProcess * GetProcessDefinedStep() const
#define InitProcessState(destination, source)
Definition: G4VITProcess.hh:53
G4StepPoint * GetPostStepPoint() const
virtual void ComputeStep(const G4Track &, const G4Step &, const double, double &)
void ProposeEnergy(G4double finalEnergy)
void Diffusion(const G4Track &track)
tuple z
Definition: test.py:28
G4DNABrownianTransportation(const G4String &aName="DNABrownianTransportation", G4int verbosityLevel=1)
#define G4endl
Definition: G4ios.hh:61
G4double GetGlobalTime() const
#define GREEN_ON_BLUE
virtual G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4bool ProposesTimeStep() const
void SetMomentumChanged(G4bool b)
#define DBL_MAX
Definition: templates.hh:83
G4int GetProcessSubType() const
Definition: G4VProcess.hh:426
const std::vector< G4double > * fpWaterDensity
virtual G4double AlongStepGetPhysicalInteractionLength(const G4Track &, G4double, G4double, G4double &, G4GPILSelection *)
G4GPILSelection