Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 94218 2015-11-09 08:24:48Z 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 
46 #include <G4Scheduler.hh>
47 #include "G4PhysicalConstants.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "Randomize.hh"
50 #include "G4Molecule.hh"
51 #include "G4RandomDirection.hh"
52 #include "G4ParticleTable.hh"
53 #include "G4SafetyHelper.hh"
55 #include "G4UnitsTable.hh"
56 #include "G4NistManager.hh"
58 #include "G4ITNavigator.hh"
59 #include "G4ITSafetyHelper.hh" // Not used yet
60 #include "G4TrackingInformation.hh"
61 
62 using namespace std;
63 
64 #ifndef State
65 #define State(theXInfo) (GetState<G4ITBrownianState>()->theXInfo)
66 #endif
67 
68 //#ifndef State
69 //#define State(theXInfo) (fTransportationState->theXInfo)
70 //#endif
71 
72 //#define USE_COLOR 1
73 
74 #ifdef USE_COLOR
75 #define RED "\033[0;31m"
76 #define LIGHT_RED "\33[1;31m"
77 #define GREEN "\033[32;40m"
78 #define GREEN_ON_BLUE "\033[1;32;44m"
79 #define RESET_COLOR "\033[0m"
80 #else
81 #define RED ""
82 #define LIGHT_RED ""
83 #define GREEN ""
84 #define GREEN_ON_BLUE ""
85 #define RESET_COLOR ""
86 #endif
87 
88 //#define DEBUG_MEM 1
89 
90 #ifdef DEBUG_MEM
91 #include "G4MemStat.hh"
92 using namespace G4MemStat;
93 using G4MemStat::MemStat;
94 #endif
95 
96 static double InvErf(double x)
97 {
99 }
100 
101 static double InvErfc(double x)
102 {
103  return CLHEP::HepStat::inverseErf(1. - x);
104 }
105 
106 //static double Erf(double x)
107 //{
108 // return CLHEP::HepStat::erf(x);
109 //}
110 
111 static double Erfc(double x)
112 {
113  return 1 - CLHEP::HepStat::erf(1. - x);
114 }
115 
116 #ifndef State
117 #define State(theXInfo) (GetState<G4ITTransportationState>()->theXInfo)
118 #endif
119 
121  G4int verbosity) :
122  G4ITTransportation(aName, verbosity)
123 {
124 
125  fVerboseLevel = 0;
126 
127  fpState.reset(new G4ITBrownianState());
128 
129  //ctor
130  SetProcessSubType(61);
131 
133  fpWaterDensity = 0;
134 
137  fSpeedMeUp = true;
138 
140  fpBrownianAction = 0;
141 }
142 
144 {
146 }
147 
149  G4ITTransportation(right)
150 {
151  //copy ctor
152  SetProcessSubType(61);
156  fNistWater = right.fNistWater;
159  fSpeedMeUp = right.fSpeedMeUp;
161 }
162 
164 {
165  if(this == &rhs) return *this; // handle self assignment
166  //assignment operator
167  return *this;
168 }
169 
172 {
173  fPathLengthWasCorrected = false;
174  fTimeStepReachedLimit = false;
175  fComputeLastPosition = false;
176  fRandomNumber = -1;
177 }
178 
180 {
181  fpState.reset(new G4ITBrownianState());
182 // G4cout << "G4DNABrownianTransportation::StartTracking : "
183 // "Initialised track State" << G4endl;
186 }
187 
189 {
190  if(verboseLevel > 0)
191  {
192  G4cout << G4endl<< GetProcessName() << ": for "
193  << setw(24) << particle.GetParticleName()
194  << "\tSubType= " << GetProcessSubType() << G4endl;
195  }
196  // Initialize water density pointer
198  GetDensityTableFor(G4Material::GetMaterial("G4_WATER"));
199 
202 }
203 
205  const G4Step& step,
206  const double timeStep,
207  double& spaceStep)
208 {
209  // G4cout << "G4ITBrownianTransportation::ComputeStep" << G4endl;
210 
211  /* If this method is called, this step
212  * cannot be geometry limited.
213  * In case the step is limited by the geometry,
214  * this method should not be called.
215  * The fTransportEndPosition calculated in
216  * the method AlongStepIL should be taken
217  * into account.
218  * In order to do so, the flag IsLeadingStep
219  * is on. Meaning : this track has the minimum
220  * interaction length over all others.
221  */
222  if(GetIT(track)->GetTrackingInfo()->IsLeadingStep())
223  {
224  const G4VITProcess* ITProc = ((const G4VITProcess*) step.GetPostStepPoint()
226  bool makeException = true;
227 
228  if(ITProc && ITProc->ProposesTimeStep()) makeException = false;
229 
230  if(makeException)
231  {
232 
233  G4ExceptionDescription exceptionDescription;
234  exceptionDescription << "ComputeStep is called while the track has"
235  "the minimum interaction time";
236  exceptionDescription << " so it should not recompute a timeStep ";
237  G4Exception("G4DNABrownianTransportation::ComputeStep",
238  "G4DNABrownianTransportation001", FatalErrorInArgument,
239  exceptionDescription);
240  }
241  }
242 
243  State(fGeometryLimitedStep) = false;
244 
245  G4Molecule* molecule = GetMolecule(track);
246 
247  if(timeStep > 0)
248  {
249  spaceStep = DBL_MAX;
250 
251  // TODO : generalize this process to all kind of Brownian objects
252  G4double diffCoeff = molecule->GetDiffusionCoefficient(track.GetMaterial(),
253  track.GetMaterial()->GetTemperature());
254 
255  static double sqrt_2 = sqrt(2.);
256  double sqrt_Dt = sqrt(diffCoeff*timeStep);
257  double sqrt_2Dt = sqrt_2*sqrt_Dt;
258 
259  if(State(fTimeStepReachedLimit)== true)
260  {
261  //========================================================================
262  State(fGeometryLimitedStep) = true;// important
263  //========================================================================
264  spaceStep = State(fEndPointDistance);
265  // G4cout << "State(fTimeStepReachedLimit)== true" << G4endl;
266  }
267  else
268  {
269  G4double x = G4RandGauss::shoot(0,sqrt_2Dt);
270  G4double y = G4RandGauss::shoot(0,sqrt_2Dt);
271  G4double z = G4RandGauss::shoot(0,sqrt_2Dt);
272 
273  spaceStep = sqrt(x*x + y*y + z*z);
274 
275  if(spaceStep >= State(fEndPointDistance))
276  {
277  //G4cout << "spaceStep >= State(fEndPointDistance)" << G4endl;
278  //======================================================================
279  State(fGeometryLimitedStep) = true;// important
280  //======================================================================
281 /*
282  if(fSpeedMeUp)
283  {
284  G4cout << "SpeedMeUp" << G4endl;
285  }
286  else
287 */
288  if(fUseSchedulerMinTimeSteps == false)// jump over barrier NOT used
289  {
290 #ifdef G4VERBOSE
291  if (fVerboseLevel > 1)
292  {
294  << "G4ITBrownianTransportation::ComputeStep() : "
295  << "Step was limited to boundary"
296  << RESET_COLOR
297  << G4endl;
298  }
299 #endif
300  //TODO
301  if(State(fRandomNumber)>=0) // CDF is used
302  {
303  /*
304  //=================================================================
305  State(fGeometryLimitedStep) = true;// important
306  //=================================================================
307  spaceStep = State(fEndPointDistance);
308  */
309 
310  //==================================================================
311  // BE AWARE THAT THE TECHNIQUE USED BELOW IS A 1D APPROXIMATION
312  // Cumulative density function for the 3D case is not yet
313  // implemented
314  //==================================================================
315 // G4cout << GREEN_ON_BLUE
316 // << "G4ITBrownianTransportation::ComputeStep() : "
317 // << "A random number was selected"
318 // << RESET_COLOR
319 // << G4endl;
320  double value = State(fRandomNumber)+(1-State(fRandomNumber))*G4UniformRand();
321  double invErfc = InvErfc(value);
322  spaceStep = invErfc*2*sqrt_Dt;
323 
324  if(State(fTimeStepReachedLimit)== false)
325  {
326  //================================================================
327  State(fGeometryLimitedStep) = false;// important
328  //================================================================
329  }
330  //==================================================================
331  // DEBUG
332 // if(spaceStep > State(fEndPointDistance))
333 // {
334 // G4cout << "value = " << value << G4endl;
335 // G4cout << "invErfc = " << invErfc << G4endl;
336 // G4cout << "spaceStep = " << G4BestUnit(spaceStep, "Length")
337 // << G4endl;
338 // G4cout << "end point distance= " << G4BestUnit(State(fEndPointDistance), "Length")
339 // << G4endl;
340 // }
341 //
342 // assert(spaceStep <= State(fEndPointDistance));
343  //==================================================================
344 
345  }
346  else if(fUseMaximumTimeBeforeReachingBoundary == false) // CDF is used
347  {
348  double min_randomNumber = Erfc(State(fEndPointDistance)/2*sqrt_Dt);
349  double value = min_randomNumber+(1-min_randomNumber)*G4UniformRand();
350  double invErfc = InvErfc(value);
351  spaceStep = invErfc*2*sqrt_Dt;
352  if(spaceStep >= State(fEndPointDistance))
353  {
354  //================================================================
355  State(fGeometryLimitedStep) = true;// important
356  //================================================================
357  }
358  else if(State(fTimeStepReachedLimit)== false)
359  {
360  //================================================================
361  State(fGeometryLimitedStep) = false;// important
362  //================================================================
363  }
364  }
365  else // CDF is NOT used
366  {
367  //==================================================================
368  State(fGeometryLimitedStep) = true;// important
369  //==================================================================
370  spaceStep = State(fEndPointDistance);
371  //TODO
372 
373  /*
374  //==================================================================
375  // 1D approximation to place the brownian between its starting point
376  // and the geometry boundary
377  //==================================================================
378  double min_randomNumber = Erfc(State(fEndPointDistance)/2*sqrt_Dt);
379  double value = State(fRandomNumber)+(1-State(fRandomNumber))*G4UniformRand();
380  double invErfc = InvErfc(value*G4UniformRand());
381  spaceStep = invErfc*2*sqrt_Dt;
382  State(fGeometryLimitedStep) = false;
383  */
384  }
385  }
386 
387  State(fTransportEndPosition)= spaceStep*
388 // step.GetPostStepPoint()->GetMomentumDirection()
389  track.GetMomentumDirection()
390  + track.GetPosition();
391  }
392  else
393  {
394  //======================================================================
395  State(fGeometryLimitedStep) = false;// important
396  //======================================================================
397  State(fTransportEndPosition)= spaceStep*step.GetPostStepPoint()->
398  GetMomentumDirection() + track.GetPosition();
399  }
400  }
401  }
402  else
403  {
404  spaceStep = 0.;
405  State(fTransportEndPosition) = track.GetPosition();
406  State(fGeometryLimitedStep) = false;
407  }
408 
409  State(fCandidateEndGlobalTime) = step.GetPreStepPoint()->GetGlobalTime()
410  + timeStep;
411  State(fEndGlobalTimeComputed) = true;
412 
413 #ifdef G4VERBOSE
414  // DEBUG
415  if (fVerboseLevel > 1)
416  {
418  << "G4ITBrownianTransportation::ComputeStep() : "
419  << " trackID : " << track.GetTrackID() << " : Molecule name: "
420  << molecule->GetName() << G4endl
421  << "Initial position:" << G4BestUnit(track.GetPosition(), "Length")
422  << G4endl
423  << "Initial direction:" << track.GetMomentumDirection() << G4endl
424  << "Final position:" << G4BestUnit(State(fTransportEndPosition), "Length")
425  << G4endl
426  << "Initial magnitude:" << G4BestUnit(track.GetPosition().mag(), "Length")
427  << G4endl
428  << "Final magnitude:" << G4BestUnit(State(fTransportEndPosition).mag(), "Length")
429  << G4endl
430  << "Diffusion length : "
431  << G4BestUnit(spaceStep, "Length")
432  << " within time step : " << G4BestUnit(timeStep,"Time")
433  << G4endl
434  << "State(fTimeStepReachedLimit)= " << State(fTimeStepReachedLimit) << G4endl
435  << "State(fGeometryLimitedStep)=" << State(fGeometryLimitedStep) << G4endl
436  << "End point distance was: " << G4BestUnit(State(fEndPointDistance), "Length")
437  << G4endl
438  << RESET_COLOR
439  << G4endl<< G4endl;
440  }
441 #endif
442 
443 //==============================================================================
444 // DEBUG
445 //assert(spaceStep < State(fEndPointDistance)
446 // || (spaceStep >= State(fEndPointDistance) && State(fGeometryLimitedStep)));
447 //assert(track.GetMomentumDirection() == State(fTransportEndMomentumDir));
448 //==============================================================================
449 }
450 
452  const G4Step& step)
453 {
455 
456 #ifdef G4VERBOSE
457  // DEBUG
458  if (fVerboseLevel > 1)
459  {
460  G4cout << GREEN_ON_BLUE << "G4ITBrownianTransportation::PostStepDoIt() :"
461  << " trackID : " << track.GetTrackID() << " Molecule name: "
462  << GetMolecule(track)->GetName() << G4endl;
463  G4cout << "Diffusion length : "
464  << G4BestUnit(step.GetStepLength(), "Length")
465  <<" within time step : " << G4BestUnit(step.GetDeltaTime(),"Time")
466  << "\t Current global time : "
467  << G4BestUnit(track.GetGlobalTime(),"Time")
468  << RESET_COLOR
469  << G4endl<< G4endl;
470  }
471 #endif
472  return &fParticleChange;
473 }
474 
476 {
477 
478 #ifdef DEBUG_MEM
479  MemStat mem_first = MemoryUsage();
480 #endif
481 
482 #ifdef G4VERBOSE
483  // DEBUG
484  if (fVerboseLevel > 1)
485  {
486  G4cout << GREEN_ON_BLUE << setw(18)
487  << "G4DNABrownianTransportation::Diffusion :" << setw(8)
488  << GetIT(track)->GetName() << "\t trackID:" << track.GetTrackID()
489  << "\t" << " Global Time = "
490  << G4BestUnit(track.GetGlobalTime(), "Time")
491  << RESET_COLOR
492  << G4endl
493  << G4endl;
494  }
495 #endif
496 
497 /*
498  fParticleChange.ProposePosition(State(fTransportEndPosition));
499  //fParticleChange.ProposeEnergy(State(fTransportEndKineticEnergy));
500  fParticleChange.SetMomentumChanged(State(fMomentumChanged));
501 
502  fParticleChange.ProposeGlobalTime(State(fCandidateEndGlobalTime));
503  fParticleChange.ProposeLocalTime(State(fCandidateEndGlobalTime));
504  fParticleChange.ProposeTrueStepLength(track.GetStepLength());
505 */
506  G4Material* material = track.GetMaterial();
507 
508  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
509 
510  if (waterDensity == 0.0)
511  {
512  if(fpBrownianAction)
513  {
514  // Let the user Brownian action class decide what to do
517  return;
518  }
519  else
520  {
521 #ifdef G4VERBOSE
522  if(fVerboseLevel)
523  {
524  G4cout << "A track is outside water material : trackID = "
525  << track.GetTrackID() << " (" << GetMolecule(track)->GetName() <<")"
526  << G4endl;
527  G4cout << "Local Time : " << G4BestUnit(track.GetGlobalTime(), "Time")
528  << G4endl;
529  G4cout << "Step Number :" << track.GetCurrentStepNumber() << G4endl;
530  }
531 #endif
534  return;// &fParticleChange is the final returned object
535  }
536  }
537 
538 
539  #ifdef DEBUG_MEM
540  MemStat mem_intermediaire = MemoryUsage();
541  mem_diff = mem_intermediaire-mem_first;
542  G4cout << "\t\t\t >> || MEM || In G4DNABrownianTransportation::Diffusion "
543  "after dealing with waterDensity for "<< track.GetTrackID()
544  << ", diff is : " << mem_diff << G4endl;
545  #endif
546 
548  State(fMomentumChanged) = true;
550  //
551  #ifdef DEBUG_MEM
552  mem_intermediaire = MemoryUsage();
553  mem_diff = mem_intermediaire-mem_first;
554  G4cout << "\t\t\t >> || MEM || In G4DNABrownianTransportation::"
555  "After proposing new direction to fParticleChange for "
556  << track.GetTrackID() << ", diff is : " << mem_diff << G4endl;
557  #endif
558 
559  return;// &fParticleChange is the final returned object
560 }
561 
562 // NOT USED
564  G4double& presafety,
565  G4double limit)
566 {
567  G4double res = DBL_MAX;
568  if(track.GetVolume() != fpSafetyHelper->GetWorldVolume())
569  {
570  G4TrackStateManager& trackStateMan = GetIT(track)->GetTrackingInfo()
572  fpSafetyHelper->LoadTrackState(trackStateMan);
574  track.GetStep()->GetPreStepPoint()->GetPosition(),
575  track.GetMomentumDirection(),
576  limit, presafety);
578  }
579  return res;
580 }
581 
583  G4double previousStepSize,
584  G4double currentMinimumStep,
585  G4double& currentSafety,
586  G4GPILSelection* selection)
587 {
588 #ifdef G4VERBOSE
589  if(fVerboseLevel)
590  {
591  G4cout << " G4DNABrownianTransportation::AlongStepGetPhysicalInteractionLength - track ID: "
592  << track.GetTrackID() << G4endl;
593  G4cout << "In volume : " << track.GetVolume()->GetName()
594  << " position : " << G4BestUnit(track.GetPosition() , "Length") << G4endl;
595  }
596 #endif
597 
598  G4double geometryStepLength =
600  track, previousStepSize, currentMinimumStep, currentSafety,
601  selection);
602 
603  if(geometryStepLength==0)
604  {
605 // G4cout << "geometryStepLength==0" << G4endl;
606  if(State(fGeometryLimitedStep))
607  {
608 // G4cout << "if(State(fGeometryLimitedStep))" << G4endl;
609  G4TouchableHandle newTouchable = new G4TouchableHistory;
610 
611  newTouchable->UpdateYourself(State(fCurrentTouchableHandle)->GetVolume(),
612  State(fCurrentTouchableHandle)->GetHistory());
613 
614  fLinearNavigator->SetGeometricallyLimitedStep();
615  fLinearNavigator->LocateGlobalPointAndUpdateTouchableHandle(
616  track.GetPosition(), track.GetMomentumDirection(),
617  newTouchable, true);
618 
619  if(newTouchable->GetVolume() == 0)
620  {
621  return 0;
622  }
623 
624  State(fCurrentTouchableHandle) = newTouchable;
625 
626  //=======================================================================
627  // TODO: speed up navigation update
628 // geometryStepLength = fLinearNavigator->ComputeStep(track.GetPosition(),
629 // track.GetMomentumDirection(),
630 // currentMinimumStep,
631 // currentSafety);
632  //=======================================================================
633 
634 
635  //=======================================
636  // Longer but safer ...
637  geometryStepLength =
639  track, previousStepSize, currentMinimumStep, currentSafety,
640  selection);
641 
642  }
643  }
644 
645  //============================================================================
646  // DEBUG
647  // G4cout << "geometryStepLength: " << G4BestUnit(geometryStepLength, "Length")
648  // << " | trackID: " << track.GetTrackID()
649  // << G4endl;
650  //============================================================================
651 
652  G4double diffusionCoefficient = 0;
653 
654  /*
655  G4Material* material = track.GetMaterial();
656  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
657 
658  if (waterDensity == 0.0)
659  {
660  if(fpBrownianAction)
661  {
662  diffusionCoefficient = fpBrownianAction->GetDiffusionCoefficient(material,
663  GetMolecule(track));
664  }
665 
666  if(diffusionCoefficient <= 0)
667  {
668  State(fGeometryLimitedStep) = false;
669  State(theInteractionTimeLeft) = 0;
670  State(fTransportEndPosition) = track.GetPosition();
671  return 0;
672  }
673 
674  }
675  else
676  */
677  //{
678  diffusionCoefficient = GetMolecule(track)->GetDiffusionCoefficient();
679  //}
680 
681  State(fComputeLastPosition) = false;
682  State(fTimeStepReachedLimit) = false;
683 
684  if (State(fGeometryLimitedStep))
685  {
686  // 95 % of the space step distribution is lower than
687  // d_95 = 2 * sqrt(2*D*t)
688  // where t is the corresponding time step
689  // so by inversion :
691  {
692  if(fSpeedMeUp)
693  {
694  State(theInteractionTimeLeft) = (geometryStepLength * geometryStepLength)
695  / (diffusionCoefficient); // d_50 - use straight line
696  }
697  else
698  {
699  State(theInteractionTimeLeft) = (currentSafety * currentSafety)
700  / (diffusionCoefficient); // d_50 - use safety
701 
702  //=====================================================================
703  // State(theInteractionTimeLeft) = (currentSafety * currentSafety)
704  // / (8 * diffusionCoefficient); // d_95
705  //=====================================================================
706  }
707  State(fComputeLastPosition) = true;
708  }
709  else
710  // Will use a random time - this is precise but long to compute in certain
711  // circumstances (many particles - small volumes)
712  {
713  State(fRandomNumber) = G4UniformRand();
714  State(theInteractionTimeLeft) = 1 / (4 * diffusionCoefficient)
715  * pow(geometryStepLength / InvErfc(State(fRandomNumber)),2);
716 
717  State(fTransportEndPosition) = geometryStepLength*
718  track.GetMomentumDirection() + track.GetPosition();
719  }
720 
722  {
723  double minTimeStepAllowed = G4VScheduler::Instance()->GetLimitingTimeStep();
724  //========================================================================
725  // TODO
726 // double currentMinTimeStep = G4VScheduler::Instance()->GetTimeStep();
727  //========================================================================
728 
729  if (State(theInteractionTimeLeft) < minTimeStepAllowed)
730  {
731  State(theInteractionTimeLeft) = minTimeStepAllowed;
732  State(fTimeStepReachedLimit) = true;
733  State(fComputeLastPosition) = true;
734  }
735  }
736  else if(State(theInteractionTimeLeft) < fInternalMinTimeStep)
737  // TODO: find a better way when fForceLimitOnMinTimeSteps is not used
738  {
739  State(fTimeStepReachedLimit) = true;
740  State(theInteractionTimeLeft) = fInternalMinTimeStep;
742  {
743  State(fComputeLastPosition) = true;
744  }
745  }
746 
747  State(fCandidateEndGlobalTime) =
748  track.GetGlobalTime() + State(theInteractionTimeLeft);
749 
750  State(fEndGlobalTimeComputed) = true; // MK: ADDED ON 20/11/2014
751 
752  State(fPathLengthWasCorrected) = false;
753  }
754  else
755  {
756  // Transform geometrical step
757  geometryStepLength = 2
758  * sqrt(diffusionCoefficient * State(theInteractionTimeLeft))
759  * InvErf(G4UniformRand());
760  State(fPathLengthWasCorrected) = true;
761  //State(fEndPointDistance) = geometryStepLength;
762  State(fTransportEndPosition) = geometryStepLength*
763  track.GetMomentumDirection() + track.GetPosition();
764  }
765 
766 #ifdef G4VERBOSE
767  // DEBUG
768  if (fVerboseLevel > 1)
769  {
771  << "G4DNABrownianTransportation::AlongStepGetPhysicalInteractionLength = "
772  << G4BestUnit(geometryStepLength, "Length")
773  << " | trackID = "
774  << track.GetTrackID()
775  << RESET_COLOR
776  << G4endl;
777  }
778 #endif
779 
780 // assert(geometryStepLength < State(fEndPointDistance)
781 // || (geometryStepLength >= State(fEndPointDistance) && State(fGeometryLimitedStep)));
782 
783  return geometryStepLength;
784 }
785 
787 //
788 // Initialize ParticleChange (by setting all its members equal
789 // to corresponding members in G4Track)
792  const G4Step& step)
793 {
794 #ifdef DEBUG_MEM
795  MemStat mem_first, mem_second, mem_diff;
796 #endif
797 
798 #ifdef DEBUG_MEM
799  mem_first = MemoryUsage();
800 #endif
801 
802  if (GetIT(track)->GetTrackingInfo()->IsLeadingStep()
803  && State(fComputeLastPosition))
804  {
805  //==========================================================================
806  // DEBUG
807  //
808 // assert(fabs(State(theInteractionTimeLeft)-
809 // G4VScheduler::Instance()->GetTimeStep()) < DBL_EPSILON);
810  //==========================================================================
811 
812  double spaceStep = DBL_MAX;
813 
814  if(State(theInteractionTimeLeft) <= fInternalMinTimeStep)
815  {
816  spaceStep = State(fEndPointDistance);
817  State(fGeometryLimitedStep) = true;
818  }
819  else
820  {
821  G4double diffusionCoefficient = GetMolecule(track)->
822  GetDiffusionCoefficient();
823 
824  double sqrt_2Dt= sqrt(2 * diffusionCoefficient * State(theInteractionTimeLeft));
825  G4double x = G4RandGauss::shoot(0, sqrt_2Dt);
826  G4double y = G4RandGauss::shoot(0, sqrt_2Dt);
827  G4double z = G4RandGauss::shoot(0, sqrt_2Dt);
828 
829  spaceStep = sqrt(x * x + y * y + z * z);
830 
831  if(spaceStep >= State(fEndPointDistance))
832  {
833  State(fGeometryLimitedStep) = true;
834  if (
835  //fSpeedMeUp == false&&
837  && spaceStep >= State(fEndPointDistance))
838  {
839  spaceStep = State(fEndPointDistance);
840  }
841  }
842  else
843  {
844  State(fGeometryLimitedStep) = false;
845  }
846  }
847 
848 // assert( (spaceStep < State(fEndPointDistance) && State(fGeometryLimitedStep) == false)
849 //|| (spaceStep >= State(fEndPointDistance) && State(fGeometryLimitedStep)));
850 
851  // Calculate final position
852  //
853  State(fTransportEndPosition) = track.GetPosition()
854  + spaceStep * track.GetMomentumDirection();
855  }
856 
857  if(fVerboseLevel)
858  {
860  << "G4DNABrownianTransportation::AlongStepDoIt: "
861  "GeometryLimitedStep = "
862  << State(fGeometryLimitedStep)
863  << RESET_COLOR
864  << G4endl;
865  }
866 
867 // static G4ThreadLocal G4int noCalls = 0;
868 // noCalls++;
869 
870 // fParticleChange.Initialize(track);
871 
873 
874 #ifdef DEBUG_MEM
875  MemStat mem_intermediaire = MemoryUsage();
876  mem_diff = mem_intermediaire-mem_first;
877  G4cout << "\t\t\t >> || MEM || After calling G4ITTransportation::"
878  "AlongStepDoIt for "<< track.GetTrackID() << ", diff is : "
879  << mem_diff << G4endl;
880 #endif
881 
882  if(track.GetStepLength() != 0 // && State(fGeometryLimitedStep)
883  //========================================================================
884  // TODO: avoid changing direction after too small time steps
885 // && (G4VScheduler::Instance()->GetTimeStep() > fInternalMinTimeStep
886 // || fSpeedMeUp == false)
887  //========================================================================
888  )
889  {
890  Diffusion(track);
891  }
892  //else
893  //{
894  // fParticleChange.ProposeMomentumDirection(State(fTransportEndMomentumDir));
895  //}
896 /*
897  if (State(fParticleIsLooping))
898  {
899  if ((State(fNoLooperTrials)>= fThresholdTrials))
900  {
901  fParticleChange.ProposeTrackStatus(fStopAndKill);
902  State(fNoLooperTrials) = 0;
903 #ifdef G4VERBOSE
904  if ((fVerboseLevel > 1))
905  {
906  G4cout
907  << " G4DNABrownianTransportation is killing track that is looping or stuck "
908  << G4endl;
909  G4cout << " Number of trials = " << State(fNoLooperTrials)
910  << " No of calls to AlongStepDoIt = " << noCalls
911  << G4endl;
912  }
913 #endif
914  }
915  else
916  {
917  State(fNoLooperTrials)++;
918  }
919  }
920  else
921  {
922  State(fNoLooperTrials)=0;
923  }
924 */
925 #ifdef DEBUG_MEM
926  mem_intermediaire = MemoryUsage();
927  mem_diff = mem_intermediaire-mem_first;
928  G4cout << "\t\t\t >> || MEM || After calling G4DNABrownianTransportation::"
929  "Diffusion for "<< track.GetTrackID() << ", diff is : "
930  << mem_diff << G4endl;
931 #endif
932 
933  return &fParticleChange;
934 }
935 
virtual void UpdateYourself(G4VPhysicalVolume *pPhysVol, const G4NavigationHistory *history=0)
Definition: G4VTouchable.cc:72
virtual void Transport(const G4Track &, G4ParticleChangeForTransport &)=0
ThreeVector shoot(const G4int Ap, const G4int Af)
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
G4ParticleChangeForTransport fParticleChange
G4int verboseLevel
Definition: G4VProcess.hh:368
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
static G4VScheduler * Instance()
Definition: G4VScheduler.cc:48
G4double GetStepLength() const
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:602
virtual G4VParticleChange * AlongStepDoIt(const G4Track &track, const G4Step &)
G4double CheckNextStep(const G4ThreeVector &position, const G4ThreeVector &direction, const G4double currentMaxStep, G4double &newSafety)
size_t GetIndex() const
Definition: G4Material.hh:262
#define RESET_COLOR
static double inverseErf(double t)
const G4ThreeVector & GetPosition() const
virtual G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &)
G4double GetDiffusionCoefficient() const
Definition: G4Molecule.cc:560
virtual void StartTracking(G4Track *aTrack)
G4ThreeVector G4RandomDirection()
virtual G4VParticleChange * AlongStepDoIt(const G4Track &track, const G4Step &stepData)
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)
const G4Step * GetStep() const
int G4int
Definition: G4Types.hh:78
MemStat MemoryUsage()
Definition: G4MemStat.cc:55
static double InvErfc(double x)
static G4NistManager * Instance()
const G4String & GetParticleName() const
static constexpr double picosecond
Definition: G4SIunits.hh:161
G4ITNavigator * fLinearNavigator
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
string material
Definition: eplot.py:19
G4StepPoint * GetPreStepPoint() const
#define State(theXInfo)
const G4String & GetName() const
Definition: G4Molecule.cc:356
static double Erfc(double x)
virtual void StartTracking(G4Track *aTrack)
G4double ComputeGeomLimit(const G4Track &track, G4double &presafety, G4double limit)
G4IT * GetIT(const G4Track *track)
Definition: G4IT.cc:49
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4int GetCurrentStepNumber() const
const XML_Char int const XML_Char * value
Definition: expat.h:331
const G4String & GetName() const
const G4ThreeVector & GetPosition() const
virtual void ResetTrackState()
G4DNABrownianTransportation & operator=(const G4DNABrownianTransportation &)
void SetInstantiateProcessState(G4bool flag)
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
G4shared_ptr< G4ProcessState > fpState
Definition: G4Step.hh:76
G4double GetDeltaTime() const
G4int GetTrackID() const
G4double GetGlobalTime() const
G4Molecule * GetMolecule(const G4Track &track)
Definition: G4Molecule.cc:76
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
G4ITSafetyHelper * fpSafetyHelper
G4Material * GetMaterial() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4TrackingInformation * GetTrackingInfo()
Definition: G4IT.hh:144
static G4DNAMolecularMaterial * Instance()
const G4VProcess * GetProcessDefinedStep() const
const G4ThreeVector & GetMomentumDirection() const
virtual G4VPhysicalVolume * GetVolume(G4int depth=0) const
Definition: G4VTouchable.cc:44
virtual double GetLimitingTimeStep() const
Definition: G4VScheduler.hh:84
G4StepPoint * GetPostStepPoint() const
virtual void ComputeStep(const G4Track &, const G4Step &, const double, double &)
void ProposeEnergy(G4double finalEnergy)
G4VPhysicalVolume * GetVolume() const
void Diffusion(const G4Track &track)
tuple z
Definition: test.py:28
G4double GetTemperature() const
Definition: G4Material.hh:182
#define G4endl
Definition: G4ios.hh:61
G4double GetGlobalTime() const
#define GREEN_ON_BLUE
static double erf(double x)
virtual G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &)
static double InvErf(double x)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4bool ProposesTimeStep() const
G4TrackStateManager & GetTrackStateManager()
void SetMomentumChanged(G4bool b)
double mag() const
#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 *)
G4double GetStepLength() const
virtual void LoadTrackState(G4TrackStateManager &manager)
G4GPILSelection
G4VPhysicalVolume * GetWorldVolume()
G4DNABrownianTransportation(const G4String &aName="DNABrownianTransportation", G4int verbosityLevel=0)