ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAOneStepThermalizationModel.hpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNAOneStepThermalizationModel.hpp
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 //
27 // Author: Mathieu Karamitros
28 //
29 // WARNING : This class is released as a prototype.
30 // It might strongly evolve or even disappear in the next releases.
31 //
32 // History:
33 // -----------
34 // 13 Nov 2016 M.Karamitros created
35 //
36 // -------------------------------------------------------------------
37 
38 #include "G4PhysicalConstants.hh"
39 #include "G4SystemOfUnits.hh"
42 #include "G4NistManager.hh"
43 #include "G4DNAChemistryManager.hh"
46 #include "G4ITNavigator.hh"
47 #include "G4Navigator.hh"
48 
49 //#define MODEL_VERBOSE
50 
51 //------------------------------------------------------------------------------
52 
53 template<typename MODEL>
56  const G4String& nam) :
57 G4VEmModel(nam), fIsInitialised(false)
58 {
59  fVerboseLevel = 0;
62  SetHighEnergyLimit(exStructure.ExcitationEnergy(0));
64  fpWaterDensity = 0;
65 }
66 
67 //------------------------------------------------------------------------------
68 
69 template<typename MODEL>
71 {
72  // if(fpNavigator && fpNavigator->GetNavigatorState())
73  // delete fpNavigator->GetNavigatorState();
74 }
75 
76 //------------------------------------------------------------------------------
77 template<typename MODEL>
79 Initialise(const G4ParticleDefinition* particleDefinition,
80  const G4DataVector&)
81 {
82 #ifdef MODEL_VERBOSE
83  if(fVerboseLevel)
84  G4cout << "Calling G4DNAOneStepThermalizationModel::Initialise()"
85  << G4endl;
86 #endif
87  if (particleDefinition->GetParticleName() != "e-")
88  {
90  errMsg << "G4DNAOneStepThermalizationModel can only be applied "
91  "to electrons";
92  G4Exception("G4DNAOneStepThermalizationModel::CrossSectionPerVolume",
93  "G4DNAOneStepThermalizationModel001",
94  FatalErrorInArgument,errMsg);
95  return;
96  }
97 
98  if(!fIsInitialised)
99  {
100  fIsInitialised = true;
101  fpParticleChangeForGamma = GetParticleChangeForGamma();
102  }
103 
106  GetNavigatorForTracking();
107 
108  fpNavigator.reset(new G4Navigator());
109 
110  if(navigator){ // add these checks for testing mode
111  auto world=navigator->GetWorldVolume();
112  if(world){
113  fpNavigator->SetWorldVolume(world);
114  //fNavigator->NewNavigatorState();
115  }
116  }
117 
118  fpWaterDensity =
120  GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
121 }
122 
123 //------------------------------------------------------------------------------
124 template<typename MODEL>
127  const G4ParticleDefinition*,
128  G4double ekin,
129  G4double,
130  G4double)
131 {
132 #ifdef MODEL_VERBOSE
133  if(fVerboseLevel > 1)
134  G4cout << "Calling CrossSectionPerVolume() of G4DNAOneStepThermalizationModel"
135  << G4endl;
136 #endif
137 
138  if(ekin > HighEnergyLimit()){
139  return 0.0;
140  }
141 
142  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
143 
144  if(waterDensity!= 0.0){
145  return DBL_MAX;
146  }
147  return 0.;
148 }
149 
150 //------------------------------------------------------------------------------
151 template<typename MODEL>
153  return MODEL::GetRmean(k);
154 }
155 
156 
157 //------------------------------------------------------------------------------
158 
159 template<typename MODEL>
162 {
163  return MODEL::GetPenetration(k, displacement);
164 }
165 
166 //------------------------------------------------------------------------------
167 template<typename MODEL>
169 SampleSecondaries(std::vector<G4DynamicParticle*>*,
170  const G4MaterialCutsCouple*,
172  G4double,
173  G4double)
174 {
175 #ifdef MODEL_VERBOSE
176  if(fVerboseLevel)
177  G4cout << "Calling SampleSecondaries() of G4DNAOneStepThermalizationModel"
178  << G4endl;
179 #endif
180 
181  G4double k = particle->GetKineticEnergy();
182 
183  if (k <= HighEnergyLimit())
184  {
185  fpParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
186  fpParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
187 
189  {
190  G4ThreeVector displacement(0,0,0);
191  GetPenetration(k, displacement);
192 
193  //______________________________________________________________
194  const G4Track * theIncomingTrack =
195  fpParticleChangeForGamma->GetCurrentTrack();
196  G4ThreeVector finalPosition(theIncomingTrack->GetPosition()+displacement);
197 
198  fpNavigator->SetWorldVolume(theIncomingTrack->GetTouchable()->
199  GetVolume(theIncomingTrack->GetTouchable()->
200  GetHistoryDepth()));
201 
202  double displacementMag = displacement.mag();
203  double safety = DBL_MAX;
204  G4ThreeVector direction = displacement/displacementMag;
205 
206  //--
207  // 6/09/16 - recupere de molecular dissocation
208  double mag_displacement = displacement.mag();
209  G4ThreeVector displacement_direction = displacement/mag_displacement;
210 
211  // double step = DBL_MAX;
212  // step = fNavigator->CheckNextStep(theIncomingTrack->GetPosition(),
213  // displacement_direction,
214  // mag_displacement,
215  // safety);
216  //
217  //
218  // if(safety < mag_displacement)
219  // {
221  // finalPosition = theIncomingTrack->GetPosition()
222  // + (displacement/displacementMag)*safety*0.80;
223  // }
224  //--
225 
226  fpNavigator->ResetHierarchyAndLocate(theIncomingTrack->GetPosition(),
227  direction,
228  *((G4TouchableHistory*)
229  theIncomingTrack->GetTouchable()));
230 
231  fpNavigator->ComputeStep(theIncomingTrack->GetPosition(),
232  displacement/displacementMag,
233  displacementMag,
234  safety);
235 
236  if(safety <= displacementMag)
237  {
238  finalPosition = theIncomingTrack->GetPosition()
239  + (displacement/displacementMag)*safety*0.80;
240  }
241 
243  &finalPosition);
244 
245  fpParticleChangeForGamma->SetProposedKineticEnergy(25.e-3*eV);
246  }
247  }
248 }