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