ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNASecondOrderReaction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNASecondOrderReaction.cc
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 //
28 
29 #include <G4VScheduler.hh>
30 #include "G4SystemOfUnits.hh"
31 #include "G4Molecule.hh"
34 #include "G4DNADamage.hh"
35 #include "G4UnitsTable.hh"
36 #include "G4TrackingInformation.hh"
37 
38 #ifndef State
39 #define State(theXInfo) (GetState<SecondOrderReactionState>()->theXInfo)
40 #endif
41 
43 {
45  enableAtRestDoIt = false;
46  enableAlongStepDoIt = false;
47  enablePostStepDoIt = true;
48 
50 
52  // meaning G4DNASecondOrderReaction contains a class inheriting from G4ProcessState
53 
54  fIsInitialized = false;
56  fpMaterial = 0;
57  fReactionRate = -1.;
58  fConcentration = -1.;
60  fProposesTimeStep = true;
63 
64  verboseLevel = 0;
65 }
66 
68  G4VITProcess(aName,type)
69 {
70  Create();
71 }
72 
74  G4VITProcess(rhs)
75 {
76  Create();
77 }
78 
80 {
81  ;
82 }
84 {
85  if (this == &rhs) return *this; // handle self assignment
86 
87  //assignment operator
88  return *this;
89 }
90 
92 {
94  fIsInGoodMaterial = false;
95 }
96 
98 {
101  fIsInitialized = true;
102 }
103 
104 void
106 {
110 }
111 
112 void
114  const G4Material* mat, double reactionRate)
115 {
116  if(fIsInitialized)
117  {
118  G4ExceptionDescription exceptionDescription ;
119  exceptionDescription << "G4DNASecondOrderReaction was already initialised. ";
120  exceptionDescription << "You cannot set a reaction after initialisation.";
121  G4Exception("G4DNASecondOrderReaction::SetReaction","G4DNASecondOrderReaction001",
122  FatalErrorInArgument,exceptionDescription);
123  }
124  fpMolecularConfiguration = molConf;
125  fpMaterial = mat;
126  fReactionRate = reactionRate;
127 }
128 
130  G4double /*previousStepSize*/,
131  G4ForceCondition* pForceCond)
132 {
133  // G4cout << "G4DNASecondOrderReaction::PostStepGetPhysicalInteractionLength" << G4endl;
134  // G4cout << "For reaction : " << fpMaterial->GetName() << " + " << fpMolecularConfiguration->GetName() << G4endl;
135 
136  //_______________________________________________________________________
137  // Check whether the track is in the good material (maybe composite material)
138  const G4Material* material = track.GetMaterial();
139 
140  G4Molecule* mol = GetMolecule(track);
141  if(!mol) return DBL_MAX;
143  {
144  // G4cout <<"mol->GetMolecularConfiguration() != fpMolecularConfiguration" << G4endl;
145  return DBL_MAX;
146  }
147 
148  G4double molDensity = (*fpMoleculeDensity)[material->GetIndex()];
149 
150  if(molDensity == 0.0) // ie : not found
151  {
152  if(State(fIsInGoodMaterial))
153  {
155  //State(fPreviousTimeAtPreStepPoint) = -1;
156  State(fIsInGoodMaterial) = false;
157  }
158 
159  // G4cout << " Material " << fpMaterial->GetName() << " not found "
160  // <<" | name of current material : " << material->GetName()
161  // << G4endl;
162 
163  return DBL_MAX; // Becareful return here !!
164  }
165 
166  // G4cout << " Va calculer le temps d'interaction " << G4endl;
167 
168  State(fIsInGoodMaterial) = true;
169 
170  // fConcentration = molDensity/fMolarMassOfMaterial;
171  fConcentration = molDensity/CLHEP::Avogadro;
172  // G4cout << "Concentration : " << fConcentration / (g/mole)<< G4endl;
173 
174  //_______________________________________________________________________
175  // Either initialize the lapse of time left
176  // meaning
177  // => the track enters for the first time in the material
178  // or substract the previous time step to the previously calculated lapse
179  // of time left
180  // meaning
181  // => the track has not left this material since the previous call
182  G4double previousTimeStep(-1.);
183 
184  if(State(fPreviousTimeAtPreStepPoint) != -1)
185  {
186  previousTimeStep = track.GetGlobalTime() -
187  State(fPreviousTimeAtPreStepPoint) ;
188  }
189 
190  State(fPreviousTimeAtPreStepPoint) = track.GetGlobalTime();
191 
192  // condition is set to "Not Forced"
193  *pForceCond = NotForced;
194 
195  if (
196  (previousTimeStep < 0.0) ||
197  (fpState->theNumberOfInteractionLengthLeft<=0.0)) {
198  // beggining of tracking (or just after DoIt of this process)
200  } else if ( previousTimeStep > 0.0) {
201  // get mean free path
202  // subtract NumberOfInteractionLengthLeft
203  SubtractNumberOfInteractionLengthLeft( previousTimeStep );
204  } else {
205  // zero time step
206  // Force trigerring the process
207  //*pForceCond = Forced;
208  }
209 
210  fpState->currentInteractionLength = 1/(fReactionRate*fConcentration);
211 
212  // G4cout << "fpState->currentInteractionLength = "
213  // << fpState->currentInteractionLength << G4endl;
214 
215  G4double value;
216  if (fpState->currentInteractionLength <DBL_MAX) {
217  value = fpState->theNumberOfInteractionLengthLeft
218  * (fpState->currentInteractionLength);
219  } else {
220  value = DBL_MAX;
221  }
222 #ifdef G4VERBOSE
223  if (verboseLevel>2) {
224  G4cout << "G4VITRestDiscreteProcess::PostStepGetPhysicalInteractionLength ";
225  G4cout << "[ " << GetProcessName() << "]" <<G4endl;
226  track.GetDynamicParticle()->DumpInfo();
227  G4cout << " in Material " << track.GetMaterial()->GetName() <<G4endl;
228  G4cout << "InteractionLength= " << value/cm <<"[cm] " <<G4endl;
229  }
230 #endif
231 
232 // G4cout << "currentInteractionLength : " << fpState->currentInteractionLength << G4endl;
233 // G4cout << "Material : " << fpMaterial->GetName()
234 // << "ID: " << track.GetTrackID()
235 // << " Returned time : " << G4BestUnit(value,"Time") << G4endl;
236 
237  if(value < fReturnedValue)
239 
240  return value*-1;
241  // multiple by -1 to indicate to the tracking system that we are returning a time
242 }
243 
245 {
246  G4Molecule* molecule = GetMolecule(track);
247 #ifdef G4VERBOSE
248  if(verboseLevel > 1)
249  {
250  G4cout << "___________" << G4endl;
251  G4cout << ">>> Beginning of G4DNASecondOrderReaction verbose" << G4endl;
252  G4cout << ">>> Returned value : " << G4BestUnit(fReturnedValue,"Time") << G4endl;
253  G4cout << ">>> Time Step : " << G4BestUnit(G4VScheduler::Instance()->GetTimeStep(),"Time") << G4endl;
254  G4cout << ">>> Reaction : " << molecule->GetName() << " + " << fpMaterial->GetName() << G4endl;
255  G4cout << ">>> End of G4DNASecondOrderReaction verbose <<<" << G4endl;
256  }
257 #endif
262  State(fPreviousTimeAtPreStepPoint) = -1;
263  return &fParticleChange;
264 }
265