| 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 // This example is provided by the Geant4-DNA 26 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained us 27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collabo 28 // shall cite the following Geant4-DNA collaboration publication:
29 // Med. Phys. 37 (2010) 4692-4708 29 // Med. Phys. 37 (2010) 4692-4708
30 // The Geant4-DNA web site is available at htt 30 // The Geant4-DNA web site is available at http://geant4-dna.org
31 // 31 //
32 // Authors: Henri Payno and Yann Perrot 32 // Authors: Henri Payno and Yann Perrot
33 // 33 //
>> 34 // $Id$
34 // 35 //
35 /// \file ClusteringAlgo.cc 36 /// \file ClusteringAlgo.cc
36 /// \brief Implementation of the ClustreringAl 37 /// \brief Implementation of the ClustreringAlgo class
37 38
38 #include "ClusteringAlgo.hh" 39 #include "ClusteringAlgo.hh"
39 <<
40 #include "ClusteringAlgoMessenger.hh" 40 #include "ClusteringAlgoMessenger.hh"
41 41
42 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh"
43 #include "Randomize.hh" 43 #include "Randomize.hh"
44 44
45 #include <map> 45 #include <map>
46 46
47 using namespace std; 47 using namespace std;
48 48
49 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 50
51 ClusteringAlgo::ClusteringAlgo(G4double pEps, << 51 ClusteringAlgo::ClusteringAlgo(G4double pEps,G4int pMinPts,
52 G4double pEMinD << 52 G4double pSPointsProb,G4double pEMinDamage, G4double pEMaxDamage)
53 : fEps(pEps), << 53 :fEps(pEps),fMinPts(pMinPts),
54 fMinPts(pMinPts), << 54 fSPointsProb(pSPointsProb),fEMinDamage(pEMinDamage),fEMaxDamage(pEMaxDamage)
55 fSPointsProb(pSPointsProb), <<
56 fEMinDamage(pEMinDamage), <<
57 fEMaxDamage(pEMaxDamage) <<
58 { 55 {
59 fNextSBPointID = 0; 56 fNextSBPointID = 0;
60 fpClustAlgoMessenger = new ClusteringAlgoMes 57 fpClustAlgoMessenger = new ClusteringAlgoMessenger(this);
61 } 58 }
62 59
63 //....oooOO0OOooo........oooOO0OOooo........oo 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 61
65 ClusteringAlgo::~ClusteringAlgo() 62 ClusteringAlgo::~ClusteringAlgo()
66 { 63 {
67 delete fpClustAlgoMessenger; 64 delete fpClustAlgoMessenger;
68 Purge(); 65 Purge();
69 } 66 }
70 67
71 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72 69
73 // Random sampling in space 70 // Random sampling in space
74 G4bool ClusteringAlgo::IsInSensitiveArea() 71 G4bool ClusteringAlgo::IsInSensitiveArea()
75 { 72 {
76 return fSPointsProb > G4UniformRand(); 73 return fSPointsProb > G4UniformRand();
77 } 74 }
78 75
79 //....oooOO0OOooo........oooOO0OOooo........oo 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80 77
81 // Random sampling in energy 78 // Random sampling in energy
82 G4bool ClusteringAlgo::IsEdepSufficient(G4doub 79 G4bool ClusteringAlgo::IsEdepSufficient(G4double pEdep)
83 { 80 {
84 if (pEdep < fEMinDamage) { << 81 if(pEdep<fEMinDamage)
>> 82 {
85 return false; 83 return false;
86 } 84 }
87 85
88 else if (pEdep > fEMaxDamage) { << 86 else if(pEdep>fEMaxDamage)
>> 87 {
89 return true; 88 return true;
90 } 89 }
91 else { << 90 else
92 G4double proba = (pEdep / eV - fEMinDamage << 91 {
93 return (proba > G4UniformRand()); << 92 G4double proba = (pEdep/eV - fEMinDamage/eV)/
>> 93 (fEMaxDamage/eV-fEMinDamage/eV);
>> 94 return (proba>G4UniformRand());
94 } 95 }
>> 96
95 } 97 }
96 98
97 //....oooOO0OOooo........oooOO0OOooo........oo 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98 100
99 // Add an event interaction to the unregistere 101 // Add an event interaction to the unregistered damage if
100 // good conditions (pos and energy) are met 102 // good conditions (pos and energy) are met
101 // 103 //
102 104
103 void ClusteringAlgo::RegisterDamage(G4ThreeVec << 105 void ClusteringAlgo::RegisterDamage(G4ThreeVector pPos,G4double pEdep)
104 { 106 {
105 if (IsEdepSufficient(pEdep)) { << 107 if(IsEdepSufficient(pEdep))
106 if (IsInSensitiveArea()) { << 108 {
107 fpSetOfPoints.push_back(new SBPoint(fNex << 109 if(IsInSensitiveArea())
>> 110 {
>> 111 fpSetOfPoints.push_back( new SBPoint(fNextSBPointID++, pPos, pEdep));
108 } 112 }
109 } 113 }
110 } 114 }
111 115
112 //....oooOO0OOooo........oooOO0OOooo........oo 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
113 117
114 map<G4int, G4int> ClusteringAlgo::RunClusterin << 118 map<G4int,G4int> ClusteringAlgo::RunClustering()
115 { 119 {
>> 120
116 // quick sort style 121 // quick sort style
117 // create cluster 122 // create cluster
118 std::vector<SBPoint*>::iterator itVisitorPt, 123 std::vector<SBPoint*>::iterator itVisitorPt, itObservedPt;
119 for (itVisitorPt = fpSetOfPoints.begin(); it << 124 for(itVisitorPt = fpSetOfPoints.begin();
>> 125 itVisitorPt != fpSetOfPoints.end();
>> 126 ++itVisitorPt )
>> 127 {
120 itObservedPt = itVisitorPt; 128 itObservedPt = itVisitorPt;
121 itObservedPt++; << 129 itObservedPt ++;
122 while (itObservedPt != fpSetOfPoints.end() << 130 while(itObservedPt != fpSetOfPoints.end() )
>> 131 {
123 // if at least one of the two points has 132 // if at least one of the two points has not a cluster
124 if (!((*itObservedPt)->HasCluster() && ( << 133 if(!((*itObservedPt)->HasCluster() && (*itVisitorPt)->HasCluster()))
125 if (AreOnTheSameCluster((*itObservedPt << 134 {
126 fEps)) << 135 if(AreOnTheSameCluster( (*itObservedPt)->GetPosition(),
>> 136 (*itVisitorPt)->GetPosition(),fEps))
127 { 137 {
128 // if none has a cluster. Create a n 138 // if none has a cluster. Create a new one
129 if (!(*itObservedPt)->HasCluster() & << 139 if(!(*itObservedPt)->HasCluster() && !(*itVisitorPt)->HasCluster())
>> 140 {
130 // create the new cluster 141 // create the new cluster
131 set<SBPoint*> clusterPoints; 142 set<SBPoint*> clusterPoints;
132 clusterPoints.insert((*itObservedP 143 clusterPoints.insert((*itObservedPt));
133 clusterPoints.insert((*itVisitorPt 144 clusterPoints.insert((*itVisitorPt));
134 ClusterSBPoints* lCluster = new Cl 145 ClusterSBPoints* lCluster = new ClusterSBPoints(clusterPoints);
135 assert(lCluster); 146 assert(lCluster);
136 fpClusters.push_back(lCluster); 147 fpClusters.push_back(lCluster);
137 assert(lCluster); 148 assert(lCluster);
138 // inform SB point that they are p 149 // inform SB point that they are part of a cluster now
139 assert(lCluster); 150 assert(lCluster);
140 (*itObservedPt)->SetCluster(lClust 151 (*itObservedPt)->SetCluster(lCluster);
141 assert(lCluster); 152 assert(lCluster);
142 (*itVisitorPt)->SetCluster(lCluste 153 (*itVisitorPt)->SetCluster(lCluster);
143 } << 154 }else
144 else { << 155 {
145 // add the point to the existing c 156 // add the point to the existing cluster
146 if ((*itObservedPt)->HasCluster()) << 157 if((*itObservedPt)->HasCluster())
>> 158 {
147 (*itObservedPt)->GetCluster()->A 159 (*itObservedPt)->GetCluster()->AddSBPoint((*itVisitorPt));
148 (*itVisitorPt)->SetCluster((*itO 160 (*itVisitorPt)->SetCluster((*itObservedPt)->GetCluster());
149 } 161 }
150 162
151 if ((*itVisitorPt)->HasCluster()) << 163 if((*itVisitorPt)->HasCluster())
>> 164 {
152 (*itVisitorPt)->GetCluster()->Ad 165 (*itVisitorPt)->GetCluster()->AddSBPoint((*itObservedPt));
153 (*itObservedPt)->SetCluster((*it 166 (*itObservedPt)->SetCluster((*itVisitorPt)->GetCluster());
154 } 167 }
155 } 168 }
156 } 169 }
157 } 170 }
158 ++itObservedPt; 171 ++itObservedPt;
159 } 172 }
160 } 173 }
161 174
162 // associate isolated points and merge clust 175 // associate isolated points and merge clusters
163 IncludeUnassociatedPoints(); 176 IncludeUnassociatedPoints();
164 MergeClusters(); 177 MergeClusters();
165 178
166 // return cluster size distribution 179 // return cluster size distribution
167 return GetClusterSizeDistribution(); 180 return GetClusterSizeDistribution();
168 } 181 }
169 182
170 //....oooOO0OOooo........oooOO0OOooo........oo 183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
171 184
172 // try to merge cluster between them, based on 185 // try to merge cluster between them, based on the distance between barycenters
173 void ClusteringAlgo::MergeClusters() 186 void ClusteringAlgo::MergeClusters()
174 { 187 {
175 std::vector<ClusterSBPoints*>::iterator itCl 188 std::vector<ClusterSBPoints*>::iterator itCluster1, itCluster2;
176 for (itCluster1 = fpClusters.begin(); itClus << 189 for(itCluster1 = fpClusters.begin();
>> 190 itCluster1 != fpClusters.end();
>> 191 ++itCluster1)
>> 192 {
177 G4ThreeVector baryCenterClust1 = (*itClust 193 G4ThreeVector baryCenterClust1 = (*itCluster1)->GetBarycenter();
178 itCluster2 = itCluster1; 194 itCluster2 = itCluster1;
179 itCluster2++; 195 itCluster2++;
180 while (itCluster2 != fpClusters.end()) { << 196 while(itCluster2 != fpClusters.end())
>> 197 {
181 G4ThreeVector baryCenterClust2 = (*itClu 198 G4ThreeVector baryCenterClust2 = (*itCluster2)->GetBarycenter();
182 // if we can merge both cluster 199 // if we can merge both cluster
183 if (AreOnTheSameCluster(baryCenterClust1 << 200 if(AreOnTheSameCluster(baryCenterClust1, baryCenterClust2,fEps))
>> 201 {
184 (*itCluster1)->MergeWith(*itCluster2); 202 (*itCluster1)->MergeWith(*itCluster2);
185 delete *itCluster2; 203 delete *itCluster2;
186 fpClusters.erase(itCluster2); 204 fpClusters.erase(itCluster2);
187 return MergeClusters(); 205 return MergeClusters();
188 } << 206 }else
189 else { << 207 {
190 itCluster2++; 208 itCluster2++;
191 } 209 }
192 } 210 }
193 } 211 }
194 } 212 }
195 213
196 //....oooOO0OOooo........oooOO0OOooo........oo 214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
197 215
198 void ClusteringAlgo::IncludeUnassociatedPoints 216 void ClusteringAlgo::IncludeUnassociatedPoints()
199 { 217 {
200 std::vector<SBPoint*>::iterator itVisitorPt; 218 std::vector<SBPoint*>::iterator itVisitorPt;
>> 219 int nbPtSansCluster = 0;
201 // Associate all point not in a cluster if p 220 // Associate all point not in a cluster if possible ( to the first found cluster)
202 for (itVisitorPt = fpSetOfPoints.begin(); it << 221 for(itVisitorPt = fpSetOfPoints.begin();
203 if (!(*itVisitorPt)->HasCluster()) { << 222 itVisitorPt != fpSetOfPoints.end();
>> 223 ++itVisitorPt)
>> 224 {
>> 225 if(!(*itVisitorPt)->HasCluster())
>> 226 {
>> 227 nbPtSansCluster ++;
204 FindCluster(*itVisitorPt); 228 FindCluster(*itVisitorPt);
205 } 229 }
206 } 230 }
207 } 231 }
208 232
209 //....oooOO0OOooo........oooOO0OOooo........oo 233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
210 234
211 bool ClusteringAlgo::FindCluster(SBPoint* pPt) 235 bool ClusteringAlgo::FindCluster(SBPoint* pPt)
212 { 236 {
213 assert(!pPt->HasCluster()); 237 assert(!pPt->HasCluster());
214 std::vector<ClusterSBPoints*>::iterator itCl 238 std::vector<ClusterSBPoints*>::iterator itCluster;
215 for (itCluster = fpClusters.begin(); itClust << 239 for(itCluster = fpClusters.begin();
216 // if((*itCluster)->hasIn(pPt, fEps)) << 240 itCluster != fpClusters.end();
217 if ((*itCluster)->HasInBarycenter(pPt, fEp << 241 ++itCluster)
>> 242 {
>> 243 //if((*itCluster)->hasIn(pPt, fEps))
>> 244 if((*itCluster)->HasInBarycenter(pPt, fEps))
>> 245 {
218 (*itCluster)->AddSBPoint(pPt); 246 (*itCluster)->AddSBPoint(pPt);
219 return true; 247 return true;
220 } 248 }
221 } << 249 }
222 return false; 250 return false;
223 } 251 }
224 252
225 //....oooOO0OOooo........oooOO0OOooo........oo 253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
226 254
227 bool ClusteringAlgo::AreOnTheSameCluster(G4Thr << 255 bool ClusteringAlgo::AreOnTheSameCluster(G4ThreeVector pPt1,
>> 256 G4ThreeVector pPt2, G4double pMinDist)
228 { 257 {
229 G4double x1 = pPt1.x() / nm; << 258 G4double x1=pPt1.x()/nm;
230 G4double y1 = pPt1.y() / nm; << 259 G4double y1=pPt1.y()/nm;
231 G4double z1 = pPt1.z() / nm; << 260 G4double z1=pPt1.z()/nm;
232 << 261
233 G4double x2 = pPt2.x() / nm; << 262 G4double x2=pPt2.x()/nm;
234 G4double y2 = pPt2.y() / nm; << 263 G4double y2=pPt2.y()/nm;
235 G4double z2 = pPt2.z() / nm; << 264 G4double z2=pPt2.z()/nm;
236 265
237 // if the two points are closed enough 266 // if the two points are closed enough
238 if (((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 << 267 if(((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2))<=
239 <= (pMinDist / nm * pMinDist / nm)) << 268 (pMinDist/nm*pMinDist/nm))
240 { 269 {
241 return true; 270 return true;
242 } << 271 }else
243 else { << 272 {
244 return false; 273 return false;
245 } 274 }
246 } 275 }
247 276
248 //....oooOO0OOooo........oooOO0OOooo........oo 277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
249 278
250 G4int ClusteringAlgo::GetSSB() const 279 G4int ClusteringAlgo::GetSSB() const
251 { 280 {
252 G4int nbSSB = 0; 281 G4int nbSSB = 0;
253 std::vector<SBPoint*>::const_iterator itSDSP 282 std::vector<SBPoint*>::const_iterator itSDSPt;
254 for (itSDSPt = fpSetOfPoints.begin(); itSDSP << 283 for(itSDSPt = fpSetOfPoints.begin();
255 if (!(*itSDSPt)->HasCluster()) { << 284 itSDSPt != fpSetOfPoints.end();
>> 285 ++itSDSPt)
>> 286 {
>> 287 if(!(*itSDSPt)->HasCluster())
>> 288 {
256 nbSSB++; 289 nbSSB++;
257 } 290 }
258 } 291 }
259 return nbSSB; 292 return nbSSB;
260 } 293 }
261 294
262 //....oooOO0OOooo........oooOO0OOooo........oo 295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
263 296
264 G4int ClusteringAlgo::GetComplexSSB() const 297 G4int ClusteringAlgo::GetComplexSSB() const
265 { 298 {
266 G4int nbSSB = 0; 299 G4int nbSSB = 0;
267 std::vector<ClusterSBPoints*>::const_iterato 300 std::vector<ClusterSBPoints*>::const_iterator itCluster;
268 for (itCluster = fpClusters.begin(); itClust << 301 for(itCluster = fpClusters.begin();
269 if ((*itCluster)->IsSSB()) { << 302 itCluster != fpClusters.end();
270 nbSSB++; << 303 ++itCluster)
>> 304 {
>> 305 if((*itCluster)->IsSSB())
>> 306 {
>> 307 nbSSB ++;
271 } 308 }
272 } 309 }
273 return nbSSB; 310 return nbSSB;
274 } 311 }
275 312
276 //....oooOO0OOooo........oooOO0OOooo........oo 313 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
277 314
278 G4int ClusteringAlgo::GetDSB() const 315 G4int ClusteringAlgo::GetDSB() const
279 { 316 {
280 G4int nbDSB = 0; 317 G4int nbDSB = 0;
281 std::vector<ClusterSBPoints*>::const_iterato 318 std::vector<ClusterSBPoints*>::const_iterator itCluster;
282 for (itCluster = fpClusters.begin(); itClust << 319 for(itCluster = fpClusters.begin();
283 if ((*itCluster)->IsDSB()) { << 320 itCluster != fpClusters.end();
284 nbDSB++; << 321 ++itCluster)
>> 322 {
>> 323 if((*itCluster)->IsDSB())
>> 324 {
>> 325 nbDSB ++;
285 } 326 }
286 } 327 }
287 return nbDSB; 328 return nbDSB;
288 } 329 }
289 330
290 //....oooOO0OOooo........oooOO0OOooo........oo 331 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
291 332
292 map<G4int, G4int> ClusteringAlgo::GetClusterSi << 333 map<G4int,G4int> ClusteringAlgo::GetClusterSizeDistribution()
293 { 334 {
294 std::map<G4int, G4int> sizeDistribution; << 335 std::map<G4int,G4int> sizeDistribution;
295 sizeDistribution[1] = GetSSB(); << 336 sizeDistribution[1]=GetSSB();
296 std::vector<ClusterSBPoints*>::const_iterato 337 std::vector<ClusterSBPoints*>::const_iterator itCluster;
297 for (itCluster = fpClusters.begin(); itClust << 338 for(itCluster = fpClusters.begin();
>> 339 itCluster != fpClusters.end();
>> 340 itCluster++)
>> 341 {
298 sizeDistribution[(*itCluster)->GetSize()]+ 342 sizeDistribution[(*itCluster)->GetSize()]++;
299 } 343 }
300 return sizeDistribution; 344 return sizeDistribution;
301 } 345 }
302 346
303 //....oooOO0OOooo........oooOO0OOooo........oo 347 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
304 348
305 void ClusteringAlgo::Purge() 349 void ClusteringAlgo::Purge()
306 { 350 {
307 fNextSBPointID = 0; 351 fNextSBPointID = 0;
308 std::vector<ClusterSBPoints*>::iterator itCl 352 std::vector<ClusterSBPoints*>::iterator itCluster;
309 for (itCluster = fpClusters.begin(); itClust << 353 for(itCluster = fpClusters.begin();
>> 354 itCluster != fpClusters.end();
>> 355 ++itCluster)
>> 356 {
310 delete *itCluster; 357 delete *itCluster;
311 *itCluster = NULL; 358 *itCluster = NULL;
312 } << 359 }
313 fpClusters.clear(); 360 fpClusters.clear();
314 std::vector<SBPoint*>::iterator itPt; 361 std::vector<SBPoint*>::iterator itPt;
315 for (itPt = fpSetOfPoints.begin(); itPt != f << 362 for(itPt = fpSetOfPoints.begin();
>> 363 itPt != fpSetOfPoints.end();
>> 364 ++itPt)
>> 365 {
316 delete *itPt; 366 delete *itPt;
317 *itPt = NULL; 367 *itPt = NULL;
318 } 368 }
319 fpSetOfPoints.clear(); 369 fpSetOfPoints.clear();
320 } 370 }
>> 371
321 372