| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // This example is provided by the Geant4-DNA
27 // Any report or published results obtained us
28 // shall cite the following Geant4-DNA collabo
29 // Med. Phys. 37 (2010) 4692-4708
30 // The Geant4-DNA web site is available at htt
31 //
32 // Authors: Henri Payno and Yann Perrot
33 //
34 //
35 /// \file ClusteringAlgo.cc
36 /// \brief Implementation of the ClustreringAl
37
38 #include "ClusteringAlgo.hh"
39
40 #include "ClusteringAlgoMessenger.hh"
41
42 #include "G4SystemOfUnits.hh"
43 #include "Randomize.hh"
44
45 #include <map>
46
47 using namespace std;
48
49 //....oooOO0OOooo........oooOO0OOooo........oo
50
51 ClusteringAlgo::ClusteringAlgo(G4double pEps,
52 G4double pEMinD
53 : fEps(pEps),
54 fMinPts(pMinPts),
55 fSPointsProb(pSPointsProb),
56 fEMinDamage(pEMinDamage),
57 fEMaxDamage(pEMaxDamage)
58 {
59 fNextSBPointID = 0;
60 fpClustAlgoMessenger = new ClusteringAlgoMes
61 }
62
63 //....oooOO0OOooo........oooOO0OOooo........oo
64
65 ClusteringAlgo::~ClusteringAlgo()
66 {
67 delete fpClustAlgoMessenger;
68 Purge();
69 }
70
71 //....oooOO0OOooo........oooOO0OOooo........oo
72
73 // Random sampling in space
74 G4bool ClusteringAlgo::IsInSensitiveArea()
75 {
76 return fSPointsProb > G4UniformRand();
77 }
78
79 //....oooOO0OOooo........oooOO0OOooo........oo
80
81 // Random sampling in energy
82 G4bool ClusteringAlgo::IsEdepSufficient(G4doub
83 {
84 if (pEdep < fEMinDamage) {
85 return false;
86 }
87
88 else if (pEdep > fEMaxDamage) {
89 return true;
90 }
91 else {
92 G4double proba = (pEdep / eV - fEMinDamage
93 return (proba > G4UniformRand());
94 }
95 }
96
97 //....oooOO0OOooo........oooOO0OOooo........oo
98
99 // Add an event interaction to the unregistere
100 // good conditions (pos and energy) are met
101 //
102
103 void ClusteringAlgo::RegisterDamage(G4ThreeVec
104 {
105 if (IsEdepSufficient(pEdep)) {
106 if (IsInSensitiveArea()) {
107 fpSetOfPoints.push_back(new SBPoint(fNex
108 }
109 }
110 }
111
112 //....oooOO0OOooo........oooOO0OOooo........oo
113
114 map<G4int, G4int> ClusteringAlgo::RunClusterin
115 {
116 // quick sort style
117 // create cluster
118 std::vector<SBPoint*>::iterator itVisitorPt,
119 for (itVisitorPt = fpSetOfPoints.begin(); it
120 itObservedPt = itVisitorPt;
121 itObservedPt++;
122 while (itObservedPt != fpSetOfPoints.end()
123 // if at least one of the two points has
124 if (!((*itObservedPt)->HasCluster() && (
125 if (AreOnTheSameCluster((*itObservedPt
126 fEps))
127 {
128 // if none has a cluster. Create a n
129 if (!(*itObservedPt)->HasCluster() &
130 // create the new cluster
131 set<SBPoint*> clusterPoints;
132 clusterPoints.insert((*itObservedP
133 clusterPoints.insert((*itVisitorPt
134 ClusterSBPoints* lCluster = new Cl
135 assert(lCluster);
136 fpClusters.push_back(lCluster);
137 assert(lCluster);
138 // inform SB point that they are p
139 assert(lCluster);
140 (*itObservedPt)->SetCluster(lClust
141 assert(lCluster);
142 (*itVisitorPt)->SetCluster(lCluste
143 }
144 else {
145 // add the point to the existing c
146 if ((*itObservedPt)->HasCluster())
147 (*itObservedPt)->GetCluster()->A
148 (*itVisitorPt)->SetCluster((*itO
149 }
150
151 if ((*itVisitorPt)->HasCluster())
152 (*itVisitorPt)->GetCluster()->Ad
153 (*itObservedPt)->SetCluster((*it
154 }
155 }
156 }
157 }
158 ++itObservedPt;
159 }
160 }
161
162 // associate isolated points and merge clust
163 IncludeUnassociatedPoints();
164 MergeClusters();
165
166 // return cluster size distribution
167 return GetClusterSizeDistribution();
168 }
169
170 //....oooOO0OOooo........oooOO0OOooo........oo
171
172 // try to merge cluster between them, based on
173 void ClusteringAlgo::MergeClusters()
174 {
175 std::vector<ClusterSBPoints*>::iterator itCl
176 for (itCluster1 = fpClusters.begin(); itClus
177 G4ThreeVector baryCenterClust1 = (*itClust
178 itCluster2 = itCluster1;
179 itCluster2++;
180 while (itCluster2 != fpClusters.end()) {
181 G4ThreeVector baryCenterClust2 = (*itClu
182 // if we can merge both cluster
183 if (AreOnTheSameCluster(baryCenterClust1
184 (*itCluster1)->MergeWith(*itCluster2);
185 delete *itCluster2;
186 fpClusters.erase(itCluster2);
187 return MergeClusters();
188 }
189 else {
190 itCluster2++;
191 }
192 }
193 }
194 }
195
196 //....oooOO0OOooo........oooOO0OOooo........oo
197
198 void ClusteringAlgo::IncludeUnassociatedPoints
199 {
200 std::vector<SBPoint*>::iterator itVisitorPt;
201 // Associate all point not in a cluster if p
202 for (itVisitorPt = fpSetOfPoints.begin(); it
203 if (!(*itVisitorPt)->HasCluster()) {
204 FindCluster(*itVisitorPt);
205 }
206 }
207 }
208
209 //....oooOO0OOooo........oooOO0OOooo........oo
210
211 bool ClusteringAlgo::FindCluster(SBPoint* pPt)
212 {
213 assert(!pPt->HasCluster());
214 std::vector<ClusterSBPoints*>::iterator itCl
215 for (itCluster = fpClusters.begin(); itClust
216 // if((*itCluster)->hasIn(pPt, fEps))
217 if ((*itCluster)->HasInBarycenter(pPt, fEp
218 (*itCluster)->AddSBPoint(pPt);
219 return true;
220 }
221 }
222 return false;
223 }
224
225 //....oooOO0OOooo........oooOO0OOooo........oo
226
227 bool ClusteringAlgo::AreOnTheSameCluster(G4Thr
228 {
229 G4double x1 = pPt1.x() / nm;
230 G4double y1 = pPt1.y() / nm;
231 G4double z1 = pPt1.z() / nm;
232
233 G4double x2 = pPt2.x() / nm;
234 G4double y2 = pPt2.y() / nm;
235 G4double z2 = pPt2.z() / nm;
236
237 // if the two points are closed enough
238 if (((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1
239 <= (pMinDist / nm * pMinDist / nm))
240 {
241 return true;
242 }
243 else {
244 return false;
245 }
246 }
247
248 //....oooOO0OOooo........oooOO0OOooo........oo
249
250 G4int ClusteringAlgo::GetSSB() const
251 {
252 G4int nbSSB = 0;
253 std::vector<SBPoint*>::const_iterator itSDSP
254 for (itSDSPt = fpSetOfPoints.begin(); itSDSP
255 if (!(*itSDSPt)->HasCluster()) {
256 nbSSB++;
257 }
258 }
259 return nbSSB;
260 }
261
262 //....oooOO0OOooo........oooOO0OOooo........oo
263
264 G4int ClusteringAlgo::GetComplexSSB() const
265 {
266 G4int nbSSB = 0;
267 std::vector<ClusterSBPoints*>::const_iterato
268 for (itCluster = fpClusters.begin(); itClust
269 if ((*itCluster)->IsSSB()) {
270 nbSSB++;
271 }
272 }
273 return nbSSB;
274 }
275
276 //....oooOO0OOooo........oooOO0OOooo........oo
277
278 G4int ClusteringAlgo::GetDSB() const
279 {
280 G4int nbDSB = 0;
281 std::vector<ClusterSBPoints*>::const_iterato
282 for (itCluster = fpClusters.begin(); itClust
283 if ((*itCluster)->IsDSB()) {
284 nbDSB++;
285 }
286 }
287 return nbDSB;
288 }
289
290 //....oooOO0OOooo........oooOO0OOooo........oo
291
292 map<G4int, G4int> ClusteringAlgo::GetClusterSi
293 {
294 std::map<G4int, G4int> sizeDistribution;
295 sizeDistribution[1] = GetSSB();
296 std::vector<ClusterSBPoints*>::const_iterato
297 for (itCluster = fpClusters.begin(); itClust
298 sizeDistribution[(*itCluster)->GetSize()]+
299 }
300 return sizeDistribution;
301 }
302
303 //....oooOO0OOooo........oooOO0OOooo........oo
304
305 void ClusteringAlgo::Purge()
306 {
307 fNextSBPointID = 0;
308 std::vector<ClusterSBPoints*>::iterator itCl
309 for (itCluster = fpClusters.begin(); itClust
310 delete *itCluster;
311 *itCluster = NULL;
312 }
313 fpClusters.clear();
314 std::vector<SBPoint*>::iterator itPt;
315 for (itPt = fpSetOfPoints.begin(); itPt != f
316 delete *itPt;
317 *itPt = NULL;
318 }
319 fpSetOfPoints.clear();
320 }
321