Geant4  9.6.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 64057 2012-10-30 15:04:49Z 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 "G4ITManager.hh"
43 #include "G4UnitsTable.hh"
44 
45 using namespace std;
46 
48  fMolReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
49 {
50  //ctor
51  fVerbose = 0;
52  fReactionModel = 0;
53  fReactionRadius = -1;
54  fDistance = -1;
55 }
56 
58 {
59  //dtor
60  if(fpChanges) delete fpChanges;
61 }
62 
64  fMolReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
65 {
66  //copy ctor
67  fVerbose = other.fVerbose ;
69  fReactionModel = 0;
70  fReactionRadius = -1;
71  fDistance = -1;
72 }
73 
75 {
76  if (this == &rhs) return *this; // handle self assignment
77 
78  fVerbose = rhs.fVerbose ;
80  fReactionRadius = -1;
81  fDistance = -1;
82 
83  //assignment operator
84  return *this;
85 }
86 
88  const G4Track& trackB,
89  const double currentStepTime,
90  const double /*previousStepTime*/,
91  bool userStepTimeLimit) /*const*/
92 {
93  G4Molecule* moleculeA = GetMolecule(trackA);
94  G4Molecule* moleculeB = GetMolecule(trackB);
95 
96  if(!fReactionModel)
97  {
98  G4ExceptionDescription exceptionDescription ("You have to give a reaction model to the molecular reaction process");
99  G4Exception("G4DNAMolecularReaction::TestReactibility","MolecularReaction001",
100  FatalErrorInArgument,exceptionDescription);
101  return false; // makes coverity happy
102  }
103  if(fMolReactionTable==0)
104  {
105  G4ExceptionDescription exceptionDescription ("You have to give a reaction table to the molecular reaction process");
106  G4Exception("G4DNAMolecularReaction::TestReactibility","MolecularReaction002",
107  FatalErrorInArgument,exceptionDescription);
108  return false; // makes coverity happy
109  }
110 
111  // Retrieve reaction range
112  fReactionRadius = -1 ; // reaction Range
113  fReactionRadius = fReactionModel -> GetReactionRadius(moleculeA, moleculeB);
114 
115  fDistance = -1 ; // separation distance
116 
117  if(currentStepTime == 0.)
118  {
119  userStepTimeLimit = false;
120  }
121 
122  G4bool output = fReactionModel->FindReaction(trackA, trackB, fReactionRadius,fDistance, userStepTimeLimit);
123 
124 #ifdef G4VERBOSE
125  // DEBUG
126  if(fVerbose > 1)
127  {
128  G4cout << "\033[1;39;36m" << "G4MolecularReaction "<< G4endl;
129  G4cout << "FindReaction returned : " << G4BestUnit(output,"Length") << G4endl;
130  G4cout << " reaction radius : " << G4BestUnit(fReactionRadius,"Length")
131  << " real distance : " << G4BestUnit((trackA.GetPosition() - trackB.GetPosition()).mag(), "Length")
132  << " calculated distance by model (= -1 if real distance > reaction radius and the user limitation step is not reached) : "
133  << G4BestUnit(fDistance,"Length")
134  << G4endl;
135 
136  G4cout << "TrackID A : " << trackA.GetTrackID()
137  << ", TrackID B : " << trackB.GetTrackID()
138  << " | MolA " << moleculeA->GetName()
139  << ", MolB " << moleculeB->GetName()
140  <<"\033[0m\n"
141  << G4endl;
142  G4cout << "--------------------------------------------" << G4endl;
143  }
144 #endif
145  return output ;
146 }
147 
148 
150 {
152  fpChanges->Initialize(trackA, trackB);
153 
154  G4Molecule* moleculeA = GetMolecule(trackA);
155  G4Molecule* moleculeB = GetMolecule(trackB);
156 
157 #ifdef G4VERBOSE
158  // DEBUG
159  if(fVerbose)
160  {
161  G4cout << "G4DNAMolecularReaction::MakeReaction" << G4endl;
162  G4cout<<"TrackA n°" << trackA.GetTrackID()
163  <<"\t | Track B n°" << trackB.GetTrackID() << G4endl;
164 
165  G4cout<<"Track A : Position : " << G4BestUnit(trackA.GetPosition(),"Length")
166  <<"\t Global Time : " << G4BestUnit(trackA.GetGlobalTime(), "Time")<< G4endl;
167 
168  G4cout<<"Track B : Position : " << G4BestUnit(trackB.GetPosition() ,"Length")
169  <<"\t Global Time : " << G4BestUnit(trackB.GetGlobalTime(), "Time")<< G4endl;
170 
171  G4cout<<"Reaction range : " << G4BestUnit(fReactionRadius,"Length")
172  << " \t Separation distance : " << G4BestUnit(fDistance,"Length")<< G4endl;
173  G4cout << "--------------------------------------------" << G4endl;
174  }
175 #endif
176 
177  const G4DNAMolecularReactionData* reactionData = fMolReactionTable->GetReactionData(moleculeA, moleculeB);
178 
179  G4int nbProducts = reactionData->GetNbProducts();
180 
181  if (nbProducts)
182  {
183  G4double D1 = moleculeA->GetDiffusionCoefficient();
184  G4double D2 = moleculeB->GetDiffusionCoefficient();
185  G4double sqrD1 = sqrt(D1);
186  G4double sqrD2 = sqrt(D2);
187  G4double numerator = sqrD1 + sqrD2;
188  G4ThreeVector reactionSite = sqrD1/numerator * trackA.GetPosition()
189  + sqrD2/numerator * trackB.GetPosition() ;
190 
191  for (G4int j=0 ; j < nbProducts; j++)
192  {
193  G4Molecule* product = new G4Molecule(*reactionData->GetProduct(j));
194  G4Track* productTrack = product->BuildTrack(trackA.GetGlobalTime(),
195  reactionSite);
196 
197  productTrack->SetTrackStatus(fAlive);
198 
199  fpChanges->AddSecondary(productTrack);
200  G4ITManager<G4Molecule>::Instance()->Push(productTrack);
201  }
202  }
203 
204  fpChanges->KillParents(true);
205 
206  return fpChanges;
207 }