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