| 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 64057 2012-10-30 15:04:49Z 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 "G4ITManager.hh"
42 #include "G4Molecule.hh" << 43 #include "G4UnitsTable.hh"
43 #include "G4MoleculeFinder.hh" << 44
44 #include "G4ITReactionChange.hh" << 45 using namespace std;
45 #include "G4ITReaction.hh" << 46
46 << 47 G4DNAMolecularReaction::G4DNAMolecularReaction():G4VITReactionProcess(),
47 #include "G4ITTrackHolder.hh" << 48 fMolReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
48 <<
49 G4DNAMolecularReaction::G4DNAMolecularReaction <<
50 : <<
51 fMolReactionTable(reference_cast<const G4 <<
52 , fpReactionModel(nullptr) <<
53 { 49 {
>> 50 //ctor
>> 51 fVerbose = 0;
>> 52 fReactionModel = 0;
>> 53 fReactionRadius = -1;
>> 54 fDistance = -1;
54 } 55 }
55 56
56 G4DNAMolecularReaction::G4DNAMolecularReaction << 57 G4DNAMolecularReaction::~G4DNAMolecularReaction()
57 : G4DNAMolecularReaction() <<
58 { 58 {
59 fpReactionModel = pReactionModel; << 59 //dtor
>> 60 if(fpChanges) delete fpChanges;
60 } 61 }
61 62
62 G4bool G4DNAMolecularReaction::TestReactibilit << 63 G4DNAMolecularReaction::G4DNAMolecularReaction(const G4DNAMolecularReaction& other):G4VITReactionProcess(other),
63 << 64 fMolReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable))
64 <<
65 <<
66 { 65 {
67 const auto pMoleculeA = GetMolecule(trackA << 66 //copy ctor
68 const auto pMoleculeB = GetMolecule(trackB << 67 fVerbose = other.fVerbose ;
69 << 68 fMolReactionTable = other.fMolReactionTable;
70 const G4double reactionRadius = fpReaction << 69 fReactionModel = 0;
>> 70 fReactionRadius = -1;
>> 71 fDistance = -1;
>> 72 }
71 73
72 G4double separationDistance = -1.; << 74 G4DNAMolecularReaction& G4DNAMolecularReaction::operator=(const G4DNAMolecularReaction& rhs)
>> 75 {
>> 76 if (this == &rhs) return *this; // handle self assignment
73 77
74 if (currentStepTime == 0.) << 78 fVerbose = rhs.fVerbose ;
75 { << 79 fMolReactionTable = rhs.fMolReactionTable;
76 userStepTimeLimit = false; << 80 fReactionRadius = -1;
77 } << 81 fDistance = -1;
78 82
79 G4bool output = fpReactionModel->FindReact << 83 //assignment operator
80 << 84 return *this;
81 return output; <<
82 } 85 }
83 86
84 std::unique_ptr<G4ITReactionChange> G4DNAMolec << 87 G4bool G4DNAMolecularReaction::TestReactibility(const G4Track& trackA,
85 << 88 const G4Track& trackB,
>> 89 const double currentStepTime,
>> 90 const double /*previousStepTime*/,
>> 91 bool userStepTimeLimit) /*const*/
86 { 92 {
87 std::unique_ptr<G4ITReactionChange> pChang << 93 G4Molecule* moleculeA = GetMolecule(trackA);
88 pChanges->Initialize(trackA, trackB); << 94 G4Molecule* moleculeB = GetMolecule(trackB);
89 95
90 const auto pMoleculeA = GetMolecule(trackA << 96 if(!fReactionModel)
91 const auto pMoleculeB = GetMolecule(trackB << 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 }
92 110
93 const auto pReactionData = fMolReactionTab << 111 // Retrieve reaction range
>> 112 fReactionRadius = -1 ; // reaction Range
>> 113 fReactionRadius = fReactionModel -> GetReactionRadius(moleculeA, moleculeB);
94 114
95 const G4int nbProducts = pReactionData->Ge << 115 fDistance = -1 ; // separation distance
96 116
97 if (nbProducts != 0) << 117 if(currentStepTime == 0.)
98 { 118 {
99 const G4double D1 = pMoleculeA->GetDif << 119 userStepTimeLimit = false;
100 const G4double D2 = pMoleculeB->GetDif << 120 }
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 121
114 G4ITTrackHolder::Instance()->Push( << 122 G4bool output = fReactionModel->FindReaction(trackA, trackB, fReactionRadius,fDistance, userStepTimeLimit);
115 123
116 pChanges->AddSecondary(pProductTra << 124 #ifdef G4VERBOSE
117 G4MoleculeFinder::Instance()->Push << 125 // DEBUG
118 } << 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;
119 } 143 }
120 << 144 #endif
121 pChanges->KillParents(true); << 145 return output ;
122 return pChanges; <<
123 } 146 }
124 147
125 void G4DNAMolecularReaction::SetReactionModel( <<
126 { <<
127 fpReactionModel = pReactionModel; <<
128 } <<
129 148
130 std::vector<std::unique_ptr<G4ITReactionChange << 149 G4ITReactionChange* G4DNAMolecularReaction::MakeReaction(const G4Track& trackA, const G4Track& trackB)
131 G4ITReactionSet* pReactionSet, <<
132 const double currentStepTime, <<
133 const double /*fGlobalTime*/, <<
134 const bool reachedUserStepTimeLimit) <<
135 { 150 {
136 std::vector<std::unique_ptr<G4ITReactionCh << 151 fpChanges = new G4ITReactionChange();
137 fReactionInfo.clear(); << 152 fpChanges->Initialize(trackA, trackB);
138 153
139 if (pReactionSet == nullptr) << 154 G4Molecule* moleculeA = GetMolecule(trackA);
>> 155 G4Molecule* moleculeB = GetMolecule(trackB);
>> 156
>> 157 #ifdef G4VERBOSE
>> 158 // DEBUG
>> 159 if(fVerbose)
140 { 160 {
141 return fReactionInfo; << 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;
142 } 174 }
>> 175 #endif
143 176
144 G4ITReactionPerTrackMap& reactionPerTrackM << 177 const G4DNAMolecularReactionData* reactionData = fMolReactionTable->GetReactionData(moleculeA, moleculeB);
145 for (auto tracks_i = reactionPerTrackMap.c <<
146 tracks_i != reactionPerTrackMap.cend( <<
147 tracks_i = reactionPerTrackMap.cbegin <<
148 { <<
149 G4Track* pTrackA = tracks_i->first; <<
150 if (pTrackA->GetTrackStatus() == fStop <<
151 { <<
152 continue; <<
153 } <<
154 178
155 G4ITReactionPerTrackPtr reactionPerTra << 179 G4int nbProducts = reactionData->GetNbProducts();
156 G4ITReactionList& reactionList = react <<
157 180
158 assert(reactionList.cbegin() != reacti << 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() ;
159 190
160 for (auto it = reactionList.cbegin(); << 191 for (G4int j=0 ; j < nbProducts; j++)
161 { 192 {
162 G4ITReactionPtr reaction(*it); << 193 G4Molecule* product = new G4Molecule(*reactionData->GetProduct(j));
163 G4Track* pTrackB = reaction->GetRe << 194 G4Track* productTrack = product->BuildTrack(trackA.GetGlobalTime(),
164 if (pTrackB->GetTrackStatus() == f << 195 reactionSite);
165 { << 196
166 continue; << 197 productTrack->SetTrackStatus(fAlive);
167 } << 198
168 << 199 fpChanges->AddSecondary(productTrack);
169 if (pTrackB == pTrackA) << 200 G4ITManager<G4Molecule>::Instance()->Push(productTrack);
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 } 201 }
194 } 202 }
195 203
196 pReactionSet->CleanAllReaction(); << 204 fpChanges->KillParents(true);
197 return fReactionInfo; << 205
>> 206 return fpChanges;
198 } 207 }
199 208