Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAMolecularReaction.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: G4DNAMolecularReaction.cc 91584 2015-07-27 13:01:48Z gcosmo $
27 //
28 // Author: Mathieu Karamitros (kara@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 #include "G4VDNAReactionModel.hh"
42 #include "G4MoleculeFinder.hh"
43 #include "G4UnitsTable.hh"
45 
46 using namespace std;
47 
50  fMolReactionTable(
52 {
53  //ctor
54  fVerbose = 0;
55  fReactionModel = 0;
56  fReactionRadius = -1;
57  fDistance = -1;
58 }
59 
61 {
62  //dtor
63  if (fpChanges) delete fpChanges;
64 }
65 
67  G4VITReactionProcess(other),
68  fMolReactionTable(
70 {
71  //copy ctor
72  fVerbose = other.fVerbose;
74  fReactionModel = 0;
75  fReactionRadius = -1;
76  fDistance = -1;
77 }
78 
80 {
81  if (this == &rhs) return *this; // handle self assignment
82 
83  fVerbose = rhs.fVerbose;
85  fReactionRadius = -1;
86  fDistance = -1;
87 
88  //assignment operator
89  return *this;
90 }
91 
93  const G4Track& trackB,
94  const double currentStepTime,
95  const double /*previousStepTime*/,
96  bool userStepTimeLimit) /*const*/
97 {
98  G4MolecularConfiguration* moleculeA =
100  G4MolecularConfiguration* moleculeB =
102 
103  if (!fReactionModel)
104  {
105  G4ExceptionDescription exceptionDescription(
106  "You have to give a reaction model to the molecular reaction process");
107  G4Exception("G4DNAMolecularReaction::TestReactibility",
108  "MolecularReaction001", FatalErrorInArgument,
109  exceptionDescription);
110  return false; // makes coverity happy
111  }
112  if (fMolReactionTable == 0)
113  {
114  G4ExceptionDescription exceptionDescription(
115  "You have to give a reaction table to the molecular reaction process");
116  G4Exception("G4DNAMolecularReaction::TestReactibility",
117  "MolecularReaction002", FatalErrorInArgument,
118  exceptionDescription);
119  return false; // makes coverity happy
120  }
121 
122  // Retrieve reaction range
123  fReactionRadius = -1; // reaction Range
124  fReactionRadius = fReactionModel->GetReactionRadius(moleculeA, moleculeB);
125 
126  fDistance = -1; // separation distance
127 
128  if (currentStepTime == 0.)
129  {
130  userStepTimeLimit = false;
131  }
132 
133  G4bool output = fReactionModel->FindReaction(trackA, trackB, fReactionRadius,
134  fDistance, userStepTimeLimit);
135 
136 #ifdef G4VERBOSE
137  // DEBUG
138  if (fVerbose > 1)
139  {
140  G4cout << "\033[1;39;36m" << "G4MolecularReaction " << G4endl;
141  G4cout << "FindReaction returned : " << G4BestUnit(output,"Length") << G4endl;
142  G4cout << " reaction radius : " << G4BestUnit(fReactionRadius,"Length")
143  << " real distance : " << G4BestUnit((trackA.GetPosition() - trackB.GetPosition()).mag(), "Length")
144  << " calculated distance by model (= -1 if real distance > reaction radius and the user limitation step is not reached) : "
145  << G4BestUnit(fDistance,"Length")
146  << G4endl;
147 
148  G4cout << "TrackID A : " << trackA.GetTrackID()
149  << ", TrackID B : " << trackB.GetTrackID()
150  << " | MolA " << moleculeA->GetName()
151  << ", MolB " << moleculeB->GetName()
152  <<"\033[0m\n"
153  << G4endl;
154  G4cout << "--------------------------------------------" << G4endl;
155  }
156 #endif
157  return output;
158 }
159 
161  const G4Track& trackB)
162 {
164  fpChanges->Initialize(trackA, trackB);
165 
166  G4MolecularConfiguration* moleculeA =
168  G4MolecularConfiguration* moleculeB =
170 
171 #ifdef G4VERBOSE
172  // DEBUG
173  if (fVerbose)
174  {
175  G4cout << "G4DNAMolecularReaction::MakeReaction" << G4endl;
176  G4cout<<"TrackA n°" << trackA.GetTrackID()
177  <<"\t | Track B n°" << trackB.GetTrackID() << G4endl;
178 
179  G4cout<<"Track A : Position : " << G4BestUnit(trackA.GetPosition(),"Length")
180  <<"\t Global Time : " << G4BestUnit(trackA.GetGlobalTime(), "Time")<< G4endl;
181 
182  G4cout<<"Track B : Position : " << G4BestUnit(trackB.GetPosition() ,"Length")
183  <<"\t Global Time : " << G4BestUnit(trackB.GetGlobalTime(), "Time")<< G4endl;
184 
185  G4cout<<"Reaction range : " << G4BestUnit(fReactionRadius,"Length")
186  << " \t Separation distance : " << G4BestUnit(fDistance,"Length")<< G4endl;
187  G4cout << "--------------------------------------------" << G4endl;
188  }
189 #endif
190 
191  const G4DNAMolecularReactionData* reactionData = fMolReactionTable
192  ->GetReactionData(moleculeA, moleculeB);
193 
194  G4int nbProducts = reactionData->GetNbProducts();
195 
196  if (nbProducts)
197  {
198  G4double D1 = moleculeA->GetDiffusionCoefficient();
199  G4double D2 = moleculeB->GetDiffusionCoefficient();
200  G4double sqrD1 = sqrt(D1);
201  G4double sqrD2 = sqrt(D2);
202  G4double numerator = sqrD1 + sqrD2;
203  G4ThreeVector reactionSite = sqrD1 / numerator * trackA.GetPosition()
204  + sqrD2 / numerator * trackB.GetPosition();
205 
206  for (G4int j = 0; j < nbProducts; j++)
207  {
208  G4Molecule* product = new G4Molecule(reactionData->GetProduct(j));
209  G4Track* productTrack = product->BuildTrack(trackA.GetGlobalTime(),
210  reactionSite);
211 
212 // G4cout << ">> G4DNAMolecularReaction::MakeReaction " << G4endl
213 // << "\t track A ("<< trackA.GetTrackID() <<"): " << trackA.GetGlobalTime() << G4endl
214 // << "\t track B ("<< trackB.GetTrackID() <<"): " << trackB.GetGlobalTime() << G4endl;
215 
216  productTrack->SetTrackStatus(fAlive);
217 
218  fpChanges->AddSecondary(productTrack);
219  G4MoleculeFinder::Instance()->Push(productTrack);
220  }
221  }
222 
223  fpChanges->KillParents(true);
224 
225  return fpChanges;
226 }
void SetTrackStatus(const G4TrackStatus aTrackStatus)
G4DNAMolecularReaction & operator=(const G4DNAMolecularReaction &other)
virtual G4bool TestReactibility(const G4Track &, const G4Track &, const double currentStepTime, const double previousStepTime, bool userStepTimeLimit)
const G4String & GetName() const
virtual void Push(G4Track *track)
const G4DNAMolecularReactionData * GetReactionData(G4MolecularConfiguration *, G4MolecularConfiguration *) const
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
const G4ThreeVector & GetPosition() const
ReturnType & reference_cast(OriginalType &source)
static G4ITFinder * Instance()
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
G4MolecularConfiguration * GetMolecularConfiguration() const
Definition: G4Molecule.cc:576
const G4ThreeVector const G4double const
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
void Initialize(const G4Track &, const G4Track &, G4VParticleChange *particleChangeA=0, G4VParticleChange *particleChangeB=0)
bool G4bool
Definition: G4Types.hh:79
virtual G4double GetReactionRadius(G4MolecularConfiguration *, G4MolecularConfiguration *)=0
G4int GetTrackID() const
G4double GetGlobalTime() const
G4Molecule * GetMolecule(const G4Track &track)
Definition: G4Molecule.cc:76
const G4DNAMolecularReactionTable *& fMolReactionTable
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4Track * BuildTrack(G4double globalTime, const G4ThreeVector &Position)
Definition: G4Molecule.cc:391
virtual G4bool FindReaction(const G4Track &, const G4Track &, const G4double, G4double &, const G4bool)=0
G4VDNAReactionModel * fReactionModel
G4ITReactionChange * fpChanges
#define G4endl
Definition: G4ios.hh:61
virtual G4ITReactionChange * MakeReaction(const G4Track &, const G4Track &)
double G4double
Definition: G4Types.hh:76
G4MolecularConfiguration * GetProduct(G4int i) const
void AddSecondary(G4Track *aSecondary)