| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collaboration publication:
29 // Med. Phys. 37 (2010) 4692-4708
30 // The Geant4-DNA web site is available at http://geant4-dna.org
31 //
32 // Authors: Henri Payno and Yann Perrot
33 //
34 //
35 /// \file ClusteringAlgo.cc
36 /// \brief Implementation of the ClustreringAlgo class
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........oooOO0OOooo........oooOO0OOooo......
50
51 ClusteringAlgo::ClusteringAlgo(G4double pEps, G4int pMinPts, G4double pSPointsProb,
52 G4double pEMinDamage, G4double pEMaxDamage)
53 : fEps(pEps),
54 fMinPts(pMinPts),
55 fSPointsProb(pSPointsProb),
56 fEMinDamage(pEMinDamage),
57 fEMaxDamage(pEMaxDamage)
58 {
59 fNextSBPointID = 0;
60 fpClustAlgoMessenger = new ClusteringAlgoMessenger(this);
61 }
62
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64
65 ClusteringAlgo::~ClusteringAlgo()
66 {
67 delete fpClustAlgoMessenger;
68 Purge();
69 }
70
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72
73 // Random sampling in space
74 G4bool ClusteringAlgo::IsInSensitiveArea()
75 {
76 return fSPointsProb > G4UniformRand();
77 }
78
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80
81 // Random sampling in energy
82 G4bool ClusteringAlgo::IsEdepSufficient(G4double pEdep)
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 / eV) / (fEMaxDamage / eV - fEMinDamage / eV);
93 return (proba > G4UniformRand());
94 }
95 }
96
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98
99 // Add an event interaction to the unregistered damage if
100 // good conditions (pos and energy) are met
101 //
102
103 void ClusteringAlgo::RegisterDamage(G4ThreeVector pPos, G4double pEdep)
104 {
105 if (IsEdepSufficient(pEdep)) {
106 if (IsInSensitiveArea()) {
107 fpSetOfPoints.push_back(new SBPoint(fNextSBPointID++, pPos, pEdep));
108 }
109 }
110 }
111
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
113
114 map<G4int, G4int> ClusteringAlgo::RunClustering()
115 {
116 // quick sort style
117 // create cluster
118 std::vector<SBPoint*>::iterator itVisitorPt, itObservedPt;
119 for (itVisitorPt = fpSetOfPoints.begin(); itVisitorPt != fpSetOfPoints.end(); ++itVisitorPt) {
120 itObservedPt = itVisitorPt;
121 itObservedPt++;
122 while (itObservedPt != fpSetOfPoints.end()) {
123 // if at least one of the two points has not a cluster
124 if (!((*itObservedPt)->HasCluster() && (*itVisitorPt)->HasCluster())) {
125 if (AreOnTheSameCluster((*itObservedPt)->GetPosition(), (*itVisitorPt)->GetPosition(),
126 fEps))
127 {
128 // if none has a cluster. Create a new one
129 if (!(*itObservedPt)->HasCluster() && !(*itVisitorPt)->HasCluster()) {
130 // create the new cluster
131 set<SBPoint*> clusterPoints;
132 clusterPoints.insert((*itObservedPt));
133 clusterPoints.insert((*itVisitorPt));
134 ClusterSBPoints* lCluster = new ClusterSBPoints(clusterPoints);
135 assert(lCluster);
136 fpClusters.push_back(lCluster);
137 assert(lCluster);
138 // inform SB point that they are part of a cluster now
139 assert(lCluster);
140 (*itObservedPt)->SetCluster(lCluster);
141 assert(lCluster);
142 (*itVisitorPt)->SetCluster(lCluster);
143 }
144 else {
145 // add the point to the existing cluster
146 if ((*itObservedPt)->HasCluster()) {
147 (*itObservedPt)->GetCluster()->AddSBPoint((*itVisitorPt));
148 (*itVisitorPt)->SetCluster((*itObservedPt)->GetCluster());
149 }
150
151 if ((*itVisitorPt)->HasCluster()) {
152 (*itVisitorPt)->GetCluster()->AddSBPoint((*itObservedPt));
153 (*itObservedPt)->SetCluster((*itVisitorPt)->GetCluster());
154 }
155 }
156 }
157 }
158 ++itObservedPt;
159 }
160 }
161
162 // associate isolated points and merge clusters
163 IncludeUnassociatedPoints();
164 MergeClusters();
165
166 // return cluster size distribution
167 return GetClusterSizeDistribution();
168 }
169
170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
171
172 // try to merge cluster between them, based on the distance between barycenters
173 void ClusteringAlgo::MergeClusters()
174 {
175 std::vector<ClusterSBPoints*>::iterator itCluster1, itCluster2;
176 for (itCluster1 = fpClusters.begin(); itCluster1 != fpClusters.end(); ++itCluster1) {
177 G4ThreeVector baryCenterClust1 = (*itCluster1)->GetBarycenter();
178 itCluster2 = itCluster1;
179 itCluster2++;
180 while (itCluster2 != fpClusters.end()) {
181 G4ThreeVector baryCenterClust2 = (*itCluster2)->GetBarycenter();
182 // if we can merge both cluster
183 if (AreOnTheSameCluster(baryCenterClust1, baryCenterClust2, fEps)) {
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........oooOO0OOooo........oooOO0OOooo......
197
198 void ClusteringAlgo::IncludeUnassociatedPoints()
199 {
200 std::vector<SBPoint*>::iterator itVisitorPt;
201 // Associate all point not in a cluster if possible ( to the first found cluster)
202 for (itVisitorPt = fpSetOfPoints.begin(); itVisitorPt != fpSetOfPoints.end(); ++itVisitorPt) {
203 if (!(*itVisitorPt)->HasCluster()) {
204 FindCluster(*itVisitorPt);
205 }
206 }
207 }
208
209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
210
211 bool ClusteringAlgo::FindCluster(SBPoint* pPt)
212 {
213 assert(!pPt->HasCluster());
214 std::vector<ClusterSBPoints*>::iterator itCluster;
215 for (itCluster = fpClusters.begin(); itCluster != fpClusters.end(); ++itCluster) {
216 // if((*itCluster)->hasIn(pPt, fEps))
217 if ((*itCluster)->HasInBarycenter(pPt, fEps)) {
218 (*itCluster)->AddSBPoint(pPt);
219 return true;
220 }
221 }
222 return false;
223 }
224
225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
226
227 bool ClusteringAlgo::AreOnTheSameCluster(G4ThreeVector pPt1, G4ThreeVector pPt2, G4double pMinDist)
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 - y2) + (z1 - z2) * (z1 - z2))
239 <= (pMinDist / nm * pMinDist / nm))
240 {
241 return true;
242 }
243 else {
244 return false;
245 }
246 }
247
248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
249
250 G4int ClusteringAlgo::GetSSB() const
251 {
252 G4int nbSSB = 0;
253 std::vector<SBPoint*>::const_iterator itSDSPt;
254 for (itSDSPt = fpSetOfPoints.begin(); itSDSPt != fpSetOfPoints.end(); ++itSDSPt) {
255 if (!(*itSDSPt)->HasCluster()) {
256 nbSSB++;
257 }
258 }
259 return nbSSB;
260 }
261
262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
263
264 G4int ClusteringAlgo::GetComplexSSB() const
265 {
266 G4int nbSSB = 0;
267 std::vector<ClusterSBPoints*>::const_iterator itCluster;
268 for (itCluster = fpClusters.begin(); itCluster != fpClusters.end(); ++itCluster) {
269 if ((*itCluster)->IsSSB()) {
270 nbSSB++;
271 }
272 }
273 return nbSSB;
274 }
275
276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
277
278 G4int ClusteringAlgo::GetDSB() const
279 {
280 G4int nbDSB = 0;
281 std::vector<ClusterSBPoints*>::const_iterator itCluster;
282 for (itCluster = fpClusters.begin(); itCluster != fpClusters.end(); ++itCluster) {
283 if ((*itCluster)->IsDSB()) {
284 nbDSB++;
285 }
286 }
287 return nbDSB;
288 }
289
290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
291
292 map<G4int, G4int> ClusteringAlgo::GetClusterSizeDistribution()
293 {
294 std::map<G4int, G4int> sizeDistribution;
295 sizeDistribution[1] = GetSSB();
296 std::vector<ClusterSBPoints*>::const_iterator itCluster;
297 for (itCluster = fpClusters.begin(); itCluster != fpClusters.end(); itCluster++) {
298 sizeDistribution[(*itCluster)->GetSize()]++;
299 }
300 return sizeDistribution;
301 }
302
303 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
304
305 void ClusteringAlgo::Purge()
306 {
307 fNextSBPointID = 0;
308 std::vector<ClusterSBPoints*>::iterator itCluster;
309 for (itCluster = fpClusters.begin(); itCluster != fpClusters.end(); ++itCluster) {
310 delete *itCluster;
311 *itCluster = NULL;
312 }
313 fpClusters.clear();
314 std::vector<SBPoint*>::iterator itPt;
315 for (itPt = fpSetOfPoints.begin(); itPt != fpSetOfPoints.end(); ++itPt) {
316 delete *itPt;
317 *itPt = NULL;
318 }
319 fpSetOfPoints.clear();
320 }
321