| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4DNAMolecularReaction.cc 87061 2014-11-24 11:43:34Z gcosmo $
26 // 27 //
27 // Author: Mathieu Karamitros (kara@cenbg.in2p 28 // Author: Mathieu Karamitros (kara@cenbg.in2p3.fr)
28 // 29 //
29 // WARNING : This class is released as a proto 30 // WARNING : This class is released as a prototype.
30 // It might strongly evolve or even disapear i 31 // It might strongly evolve or even disapear in the next releases.
31 // 32 //
32 // History: 33 // History:
33 // ----------- 34 // -----------
34 // 10 Oct 2011 M.Karamitros created 35 // 10 Oct 2011 M.Karamitros created
35 // 36 //
36 // ------------------------------------------- 37 // -------------------------------------------------------------------
37 38
38 #include "G4DNAMolecularReaction.hh" 39 #include "G4DNAMolecularReaction.hh"
39 #include "G4DNAMolecularReactionTable.hh" 40 #include "G4DNAMolecularReactionTable.hh"
40 #include "G4VDNAReactionModel.hh" 41 #include "G4VDNAReactionModel.hh"
41 #include "G4MolecularConfiguration.hh" <<
42 #include "G4Molecule.hh" <<
43 #include "G4MoleculeFinder.hh" 42 #include "G4MoleculeFinder.hh"
44 #include "G4ITReactionChange.hh" << 43 #include "G4UnitsTable.hh"
45 #include "G4ITReaction.hh" <<
46 44
47 #include "G4ITTrackHolder.hh" << 45 using namespace std;
48 46
49 G4DNAMolecularReaction::G4DNAMolecularReaction << 47 G4DNAMolecularReaction::G4DNAMolecularReaction() :
50 : << 48 G4VITReactionProcess(),
51 fMolReactionTable(reference_cast<const G4 << 49 fMolReactionTable(
52 , fpReactionModel(nullptr) << 50 reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
53 { << 51 {
>> 52 //ctor
>> 53 fVerbose = 0;
>> 54 fReactionModel = 0;
>> 55 fReactionRadius = -1;
>> 56 fDistance = -1;
54 } 57 }
55 58
56 G4DNAMolecularReaction::G4DNAMolecularReaction << 59 G4DNAMolecularReaction::~G4DNAMolecularReaction()
57 : G4DNAMolecularReaction() <<
58 { 60 {
59 fpReactionModel = pReactionModel; << 61 //dtor
>> 62 if (fpChanges) delete fpChanges;
60 } 63 }
61 64
62 G4bool G4DNAMolecularReaction::TestReactibilit << 65 G4DNAMolecularReaction::G4DNAMolecularReaction(const G4DNAMolecularReaction& other) :
63 << 66 G4VITReactionProcess(other),
64 << 67 fMolReactionTable(
65 << 68 reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
66 { << 69 {
67 const auto pMoleculeA = GetMolecule(trackA << 70 //copy ctor
68 const auto pMoleculeB = GetMolecule(trackB << 71 fVerbose = other.fVerbose;
69 << 72 fMolReactionTable = other.fMolReactionTable;
70 const G4double reactionRadius = fpReaction << 73 fReactionModel = 0;
71 << 74 fReactionRadius = -1;
72 G4double separationDistance = -1.; << 75 fDistance = -1;
73 <<
74 if (currentStepTime == 0.) <<
75 { <<
76 userStepTimeLimit = false; <<
77 } <<
78 <<
79 G4bool output = fpReactionModel->FindReact <<
80 <<
81 return output; <<
82 } 76 }
83 77
84 std::unique_ptr<G4ITReactionChange> G4DNAMolec << 78 G4DNAMolecularReaction& G4DNAMolecularReaction::operator=(const G4DNAMolecularReaction& rhs)
85 <<
86 { 79 {
87 std::unique_ptr<G4ITReactionChange> pChang << 80 if (this == &rhs) return *this; // handle self assignment
88 pChanges->Initialize(trackA, trackB); <<
89 <<
90 const auto pMoleculeA = GetMolecule(trackA <<
91 const auto pMoleculeB = GetMolecule(trackB <<
92 81
93 const auto pReactionData = fMolReactionTab << 82 fVerbose = rhs.fVerbose;
>> 83 fMolReactionTable = rhs.fMolReactionTable;
>> 84 fReactionRadius = -1;
>> 85 fDistance = -1;
94 86
95 const G4int nbProducts = pReactionData->Ge << 87 //assignment operator
96 << 88 return *this;
97 if (nbProducts != 0) <<
98 { <<
99 const G4double D1 = pMoleculeA->GetDif <<
100 const G4double D2 = pMoleculeB->GetDif <<
101 const G4double sqrD1 = D1 == 0. ? 0. : <<
102 const G4double sqrD2 = D2 == 0. ? 0. : <<
103 const G4double inv_numerator = 1./(sqr <<
104 const G4ThreeVector reactionSite = sqr <<
105 + sqr <<
106 <<
107 for (G4int j = 0; j < nbProducts; ++j) <<
108 { <<
109 auto pProduct = new G4Molecule(pRe <<
110 auto pProductTrack = pProduct->Bui <<
111 <<
112 pProductTrack->SetTrackStatus(fAli <<
113 <<
114 G4ITTrackHolder::Instance()->Push( <<
115 <<
116 pChanges->AddSecondary(pProductTra <<
117 G4MoleculeFinder::Instance()->Push <<
118 } <<
119 } <<
120 <<
121 pChanges->KillParents(true); <<
122 return pChanges; <<
123 } 89 }
124 90
125 void G4DNAMolecularReaction::SetReactionModel( << 91 G4bool G4DNAMolecularReaction::TestReactibility(const G4Track& trackA,
>> 92 const G4Track& trackB,
>> 93 const double currentStepTime,
>> 94 const double /*previousStepTime*/,
>> 95 bool userStepTimeLimit) /*const*/
126 { 96 {
127 fpReactionModel = pReactionModel; << 97 G4Molecule* moleculeA = GetMolecule(trackA);
>> 98 G4Molecule* moleculeB = GetMolecule(trackB);
>> 99
>> 100 if (!fReactionModel)
>> 101 {
>> 102 G4ExceptionDescription exceptionDescription(
>> 103 "You have to give a reaction model to the molecular reaction process");
>> 104 G4Exception("G4DNAMolecularReaction::TestReactibility",
>> 105 "MolecularReaction001", FatalErrorInArgument,
>> 106 exceptionDescription);
>> 107 return false; // makes coverity happy
>> 108 }
>> 109 if (fMolReactionTable == 0)
>> 110 {
>> 111 G4ExceptionDescription exceptionDescription(
>> 112 "You have to give a reaction table to the molecular reaction process");
>> 113 G4Exception("G4DNAMolecularReaction::TestReactibility",
>> 114 "MolecularReaction002", FatalErrorInArgument,
>> 115 exceptionDescription);
>> 116 return false; // makes coverity happy
>> 117 }
>> 118
>> 119 // Retrieve reaction range
>> 120 fReactionRadius = -1; // reaction Range
>> 121 fReactionRadius = fReactionModel->GetReactionRadius(moleculeA, moleculeB);
>> 122
>> 123 fDistance = -1; // separation distance
>> 124
>> 125 if (currentStepTime == 0.)
>> 126 {
>> 127 userStepTimeLimit = false;
>> 128 }
>> 129
>> 130 G4bool output = fReactionModel->FindReaction(trackA, trackB, fReactionRadius,
>> 131 fDistance, userStepTimeLimit);
>> 132
>> 133 #ifdef G4VERBOSE
>> 134 // DEBUG
>> 135 if (fVerbose > 1)
>> 136 {
>> 137 G4cout << "\033[1;39;36m" << "G4MolecularReaction " << G4endl;
>> 138 G4cout << "FindReaction returned : " << G4BestUnit(output,"Length") << G4endl;
>> 139 G4cout << " reaction radius : " << G4BestUnit(fReactionRadius,"Length")
>> 140 << " real distance : " << G4BestUnit((trackA.GetPosition() - trackB.GetPosition()).mag(), "Length")
>> 141 << " calculated distance by model (= -1 if real distance > reaction radius and the user limitation step is not reached) : "
>> 142 << G4BestUnit(fDistance,"Length")
>> 143 << G4endl;
>> 144
>> 145 G4cout << "TrackID A : " << trackA.GetTrackID()
>> 146 << ", TrackID B : " << trackB.GetTrackID()
>> 147 << " | MolA " << moleculeA->GetName()
>> 148 << ", MolB " << moleculeB->GetName()
>> 149 <<"\033[0m\n"
>> 150 << G4endl;
>> 151 G4cout << "--------------------------------------------" << G4endl;
>> 152 }
>> 153 #endif
>> 154 return output;
128 } 155 }
129 156
130 std::vector<std::unique_ptr<G4ITReactionChange << 157 G4ITReactionChange* G4DNAMolecularReaction::MakeReaction(const G4Track& trackA,
131 G4ITReactionSet* pReactionSet, << 158 const G4Track& trackB)
132 const double currentStepTime, <<
133 const double /*fGlobalTime*/, <<
134 const bool reachedUserStepTimeLimit) <<
135 { 159 {
136 std::vector<std::unique_ptr<G4ITReactionCh << 160 fpChanges = new G4ITReactionChange();
137 fReactionInfo.clear(); << 161 fpChanges->Initialize(trackA, trackB);
138 162
139 if (pReactionSet == nullptr) << 163 G4Molecule* moleculeA = GetMolecule(trackA);
140 { << 164 G4Molecule* moleculeB = GetMolecule(trackB);
141 return fReactionInfo; <<
142 } <<
143 165
144 G4ITReactionPerTrackMap& reactionPerTrackM << 166 #ifdef G4VERBOSE
145 for (auto tracks_i = reactionPerTrackMap.c << 167 // DEBUG
146 tracks_i != reactionPerTrackMap.cend( << 168 if (fVerbose)
147 tracks_i = reactionPerTrackMap.cbegin << 169 {
>> 170 G4cout << "G4DNAMolecularReaction::MakeReaction" << G4endl;
>> 171 G4cout<<"TrackA n°" << trackA.GetTrackID()
>> 172 <<"\t | Track B n°" << trackB.GetTrackID() << G4endl;
>> 173
>> 174 G4cout<<"Track A : Position : " << G4BestUnit(trackA.GetPosition(),"Length")
>> 175 <<"\t Global Time : " << G4BestUnit(trackA.GetGlobalTime(), "Time")<< G4endl;
>> 176
>> 177 G4cout<<"Track B : Position : " << G4BestUnit(trackB.GetPosition() ,"Length")
>> 178 <<"\t Global Time : " << G4BestUnit(trackB.GetGlobalTime(), "Time")<< G4endl;
>> 179
>> 180 G4cout<<"Reaction range : " << G4BestUnit(fReactionRadius,"Length")
>> 181 << " \t Separation distance : " << G4BestUnit(fDistance,"Length")<< G4endl;
>> 182 G4cout << "--------------------------------------------" << G4endl;
>> 183 }
>> 184 #endif
>> 185
>> 186 const G4DNAMolecularReactionData* reactionData = fMolReactionTable
>> 187 ->GetReactionData(moleculeA, moleculeB);
>> 188
>> 189 G4int nbProducts = reactionData->GetNbProducts();
>> 190
>> 191 if (nbProducts)
>> 192 {
>> 193 G4double D1 = moleculeA->GetDiffusionCoefficient();
>> 194 G4double D2 = moleculeB->GetDiffusionCoefficient();
>> 195 G4double sqrD1 = sqrt(D1);
>> 196 G4double sqrD2 = sqrt(D2);
>> 197 G4double numerator = sqrD1 + sqrD2;
>> 198 G4ThreeVector reactionSite = sqrD1 / numerator * trackA.GetPosition()
>> 199 + sqrD2 / numerator * trackB.GetPosition();
>> 200
>> 201 for (G4int j = 0; j < nbProducts; j++)
148 { 202 {
149 G4Track* pTrackA = tracks_i->first; << 203 G4Molecule* product = new G4Molecule(*reactionData->GetProduct(j));
150 if (pTrackA->GetTrackStatus() == fStop << 204 G4Track* productTrack = product->BuildTrack(trackA.GetGlobalTime(),
151 { << 205 reactionSite);
152 continue; << 206
153 } << 207 // G4cout << ">> G4DNAMolecularReaction::MakeReaction " << G4endl
154 << 208 // << "\t track A ("<< trackA.GetTrackID() <<"): " << trackA.GetGlobalTime() << G4endl
155 G4ITReactionPerTrackPtr reactionPerTra << 209 // << "\t track B ("<< trackB.GetTrackID() <<"): " << trackB.GetGlobalTime() << G4endl;
156 G4ITReactionList& reactionList = react << 210
157 << 211 productTrack->SetTrackStatus(fAlive);
158 assert(reactionList.cbegin() != reacti << 212
159 << 213 fpChanges->AddSecondary(productTrack);
160 for (auto it = reactionList.cbegin(); << 214 G4MoleculeFinder::Instance()->Push(productTrack);
161 { << 215 // G4ITManager<G4Molecule>::Instance()->Push(productTrack);
162 G4ITReactionPtr reaction(*it); <<
163 G4Track* pTrackB = reaction->GetRe <<
164 if (pTrackB->GetTrackStatus() == f <<
165 { <<
166 continue; <<
167 } <<
168 <<
169 if (pTrackB == pTrackA) <<
170 { <<
171 G4ExceptionDescription excepti <<
172 exceptionDescription <<
173 << "The IT reaction proces <<
174 exceptionDescription << "The p <<
175 G4Exception("G4ITModelProcesso <<
176 "ITModelProcessor0 <<
177 FatalErrorInArgume <<
178 exceptionDescripti <<
179 } <<
180 <<
181 pReactionSet->SelectThisReaction(s <<
182 <<
183 if (TestReactibility(*pTrackA, *pT <<
184 { <<
185 auto pReactionChange = MakeRea <<
186 <<
187 if (pReactionChange) <<
188 { <<
189 fReactionInfo.push_back(st <<
190 break; <<
191 } <<
192 } <<
193 } <<
194 } 216 }
>> 217 }
>> 218
>> 219 fpChanges->KillParents(true);
195 220
196 pReactionSet->CleanAllReaction(); << 221 return fpChanges;
197 return fReactionInfo; <<
198 } 222 }
199 223