| 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 91584 2015-07-27 13:01:48Z 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" << 44 #include "G4MolecularConfiguration.hh"
46 45
47 #include "G4ITTrackHolder.hh" << 46 using namespace std;
48 47
49 G4DNAMolecularReaction::G4DNAMolecularReaction << 48 G4DNAMolecularReaction::G4DNAMolecularReaction() :
50 : << 49 G4VITReactionProcess(),
51 fMolReactionTable(reference_cast<const G4 << 50 fMolReactionTable(
52 , fpReactionModel(nullptr) << 51 reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
53 { << 52 {
>> 53 //ctor
>> 54 fVerbose = 0;
>> 55 fReactionModel = 0;
>> 56 fReactionRadius = -1;
>> 57 fDistance = -1;
54 } 58 }
55 59
56 G4DNAMolecularReaction::G4DNAMolecularReaction << 60 G4DNAMolecularReaction::~G4DNAMolecularReaction()
57 : G4DNAMolecularReaction() <<
58 { 61 {
59 fpReactionModel = pReactionModel; << 62 //dtor
>> 63 if (fpChanges) delete fpChanges;
60 } 64 }
61 65
62 G4bool G4DNAMolecularReaction::TestReactibilit << 66 G4DNAMolecularReaction::G4DNAMolecularReaction(const G4DNAMolecularReaction& other) :
63 << 67 G4VITReactionProcess(other),
64 << 68 fMolReactionTable(
65 << 69 reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
66 { << 70 {
67 const auto pMoleculeA = GetMolecule(trackA << 71 //copy ctor
68 const auto pMoleculeB = GetMolecule(trackB << 72 fVerbose = other.fVerbose;
69 << 73 fMolReactionTable = other.fMolReactionTable;
70 const G4double reactionRadius = fpReaction << 74 fReactionModel = 0;
71 << 75 fReactionRadius = -1;
72 G4double separationDistance = -1.; << 76 fDistance = -1;
73 <<
74 if (currentStepTime == 0.) <<
75 { <<
76 userStepTimeLimit = false; <<
77 } <<
78 <<
79 G4bool output = fpReactionModel->FindReact <<
80 <<
81 return output; <<
82 } 77 }
83 78
84 std::unique_ptr<G4ITReactionChange> G4DNAMolec << 79 G4DNAMolecularReaction& G4DNAMolecularReaction::operator=(const G4DNAMolecularReaction& rhs)
85 <<
86 { 80 {
87 std::unique_ptr<G4ITReactionChange> pChang << 81 if (this == &rhs) return *this; // handle self assignment
88 pChanges->Initialize(trackA, trackB); <<
89 82
90 const auto pMoleculeA = GetMolecule(trackA << 83 fVerbose = rhs.fVerbose;
91 const auto pMoleculeB = GetMolecule(trackB << 84 fMolReactionTable = rhs.fMolReactionTable;
>> 85 fReactionRadius = -1;
>> 86 fDistance = -1;
92 87
93 const auto pReactionData = fMolReactionTab << 88 //assignment operator
94 << 89 return *this;
95 const G4int nbProducts = pReactionData->Ge <<
96 <<
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 } 90 }
124 91
125 void G4DNAMolecularReaction::SetReactionModel( << 92 G4bool G4DNAMolecularReaction::TestReactibility(const G4Track& trackA,
>> 93 const G4Track& trackB,
>> 94 const double currentStepTime,
>> 95 const double /*previousStepTime*/,
>> 96 bool userStepTimeLimit) /*const*/
126 { 97 {
127 fpReactionModel = pReactionModel; << 98 G4MolecularConfiguration* moleculeA =
>> 99 GetMolecule(trackA)->GetMolecularConfiguration();
>> 100 G4MolecularConfiguration* moleculeB =
>> 101 GetMolecule(trackB)->GetMolecularConfiguration();
>> 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;
128 } 158 }
129 159
130 std::vector<std::unique_ptr<G4ITReactionChange << 160 G4ITReactionChange* G4DNAMolecularReaction::MakeReaction(const G4Track& trackA,
131 G4ITReactionSet* pReactionSet, << 161 const G4Track& trackB)
132 const double currentStepTime, <<
133 const double /*fGlobalTime*/, <<
134 const bool reachedUserStepTimeLimit) <<
135 { 162 {
136 std::vector<std::unique_ptr<G4ITReactionCh << 163 fpChanges = new G4ITReactionChange();
137 fReactionInfo.clear(); << 164 fpChanges->Initialize(trackA, trackB);
138 165
139 if (pReactionSet == nullptr) << 166 G4MolecularConfiguration* moleculeA =
140 { << 167 GetMolecule(trackA)->GetMolecularConfiguration();
141 return fReactionInfo; << 168 G4MolecularConfiguration* moleculeB =
142 } << 169 GetMolecule(trackB)->GetMolecularConfiguration();
>> 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();
143 205
144 G4ITReactionPerTrackMap& reactionPerTrackM << 206 for (G4int j = 0; j < nbProducts; j++)
145 for (auto tracks_i = reactionPerTrackMap.c <<
146 tracks_i != reactionPerTrackMap.cend( <<
147 tracks_i = reactionPerTrackMap.cbegin <<
148 { 207 {
149 G4Track* pTrackA = tracks_i->first; << 208 G4Molecule* product = new G4Molecule(reactionData->GetProduct(j));
150 if (pTrackA->GetTrackStatus() == fStop << 209 G4Track* productTrack = product->BuildTrack(trackA.GetGlobalTime(),
151 { << 210 reactionSite);
152 continue; << 211
153 } << 212 // G4cout << ">> G4DNAMolecularReaction::MakeReaction " << G4endl
154 << 213 // << "\t track A ("<< trackA.GetTrackID() <<"): " << trackA.GetGlobalTime() << G4endl
155 G4ITReactionPerTrackPtr reactionPerTra << 214 // << "\t track B ("<< trackB.GetTrackID() <<"): " << trackB.GetGlobalTime() << G4endl;
156 G4ITReactionList& reactionList = react << 215
157 << 216 productTrack->SetTrackStatus(fAlive);
158 assert(reactionList.cbegin() != reacti << 217
159 << 218 fpChanges->AddSecondary(productTrack);
160 for (auto it = reactionList.cbegin(); << 219 G4MoleculeFinder::Instance()->Push(productTrack);
161 { <<
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 } 220 }
>> 221 }
>> 222
>> 223 fpChanges->KillParents(true);
195 224
196 pReactionSet->CleanAllReaction(); << 225 return fpChanges;
197 return fReactionInfo; <<
198 } 226 }
199 227