Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ClusterSBPoints.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 "ClusterSBPoints.hh"
40 #include "G4SystemOfUnits.hh"
41 #include <iostream>
42 
43 using namespace std;
44 
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46 
47 ClusterSBPoints::ClusterSBPoints(std::set<SBPoint*> pSBPoints):
48  fpRegisteredSBPoints()
49 {
50  UpdateDoubleStrand();
51  std::set<SBPoint*>::iterator itPt;
52  for(itPt = pSBPoints.begin(); itPt != pSBPoints.end(); ++itPt)
53  {
54  AddSBPoint(*itPt);
55  }
56 }
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59 
61 {
62  Clear();
63 }
64 
66 {
67  std::set<SBPoint*>::iterator itPt;
68  for(itPt = fpRegisteredSBPoints.begin();
69  itPt != fpRegisteredSBPoints.end();
70  ++itPt)
71  {
72  (*itPt)->CleanCluster();
73  }
74  fpRegisteredSBPoints.clear();
75 }
76 
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 
80 {
81  assert(pSBPoint);
82  fpRegisteredSBPoints.insert(pSBPoint);
83  pSBPoint->SetCluster(this);
84 
85  UpdateDoubleStrand();
86 }
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
89 
91 {
92  G4double x = 0;
93  G4double y = 0;
94  G4double z = 0;
95 
96  std::set<SBPoint* >::iterator itSDSPt;
97  for(itSDSPt = fpRegisteredSBPoints.begin();
98  itSDSPt != fpRegisteredSBPoints.end();
99  ++itSDSPt)
100  {
101  x+= (*itSDSPt)->GetPosition().x();
102  y+= (*itSDSPt)->GetPosition().y();
103  z+= (*itSDSPt)->GetPosition().z();
104  }
105 
106  return G4ThreeVector(
107  x/fpRegisteredSBPoints.size(),
108  y/fpRegisteredSBPoints.size(),
109  z/fpRegisteredSBPoints.size());
110 }
111 
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
113 
115 {
116  G4double res = 0;
117  std::set<SBPoint* >::iterator itSDSPt;
118  for(itSDSPt = fpRegisteredSBPoints.begin();
119  itSDSPt != fpRegisteredSBPoints.end();
120  ++itSDSPt)
121  {
122  res += (*itSDSPt)->GetEdep();
123  }
124  return res;
125 }
126 
127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
128 
129 void ClusterSBPoints::UpdateDoubleStrand()
130 {
131  fIsDoubleSB = false;
132  bool firstStrandTouch = false;
133  bool secondStrandTouch = false;
134 
135  std::set<SBPoint* >::iterator itSDSPt;
136  for(itSDSPt = fpRegisteredSBPoints.begin();
137  itSDSPt != fpRegisteredSBPoints.end();
138  ++itSDSPt)
139  {
140  // if the SDSPoint is localized on the first strand
141  if( ((*itSDSPt)->GetTouchedStrand() == 0 ) && !firstStrandTouch )
142  {
143  firstStrandTouch = true;
144  if(secondStrandTouch)
145  {
146  fIsDoubleSB = true;
147  return;
148  }
149  }
150  // if the SDSPoint is localized on the second strand
151  if( ((*itSDSPt)->GetTouchedStrand() == 1 ) && !secondStrandTouch )
152  {
153  secondStrandTouch = true;
154  if(firstStrandTouch)
155  {
156  fIsDoubleSB = true;
157  return;
158  }
159  }
160  }
161 }
162 
163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
164 
165 bool AreOnTheSameCluster(const SBPoint* pPt1, const SBPoint* pPt2,
166  G4double pMinDist)
167 {
168  assert(pPt1);
169  assert(pPt2);
170 
171  G4double x1=pPt1->GetPosition().x()/nm;
172  G4double y1=pPt1->GetPosition().y()/nm;
173  G4double z1=pPt1->GetPosition().z()/nm;
174 
175  G4double x2=pPt2->GetPosition().x()/nm;
176  G4double y2=pPt2->GetPosition().y()/nm;
177  G4double z2=pPt2->GetPosition().z()/nm;
178 
179  // if the two points are closed enough
180  if(((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2))<=
181  (pMinDist/nm*pMinDist/nm))
182  {
183  return true;
184  }else
185  {
186  return false;
187  }
188 }
189 
190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191 
192 void ClusterSBPoints::FindAllPointsPossible(std::vector<SBPoint*>* pPtsToCheck,
193  G4int pMinPts, G4double pMinDist)
194 {
195  assert((unsigned int)pMinPts > this->GetSize());
196  std::vector<SBPoint*>::iterator itPt = pPtsToCheck->begin();
197  while(itPt != pPtsToCheck->end() )
198  {
199 
200  // If 1- each SBpoint is part of only one cluster
201  // 2- the point isn't already in the cluster
202  // 3- the point is close enough of the barycenter
203  if( (!(*itPt)->HasCluster())
204  && (fpRegisteredSBPoints.find(*itPt) == fpRegisteredSBPoints.end())
205  && HasInBarycenter(*itPt, pMinDist)) // first version used HasIn method
206  {
207  // the point is added
208  this->AddSBPoint(*itPt);
209  if(this->GetSize() >= (unsigned int)pMinPts)
210  {
211  return;
212  }
213  // restart from scratch
214  itPt = pPtsToCheck->begin();
215  }else
216  {
217  ++itPt;
218  }
219  }
220 }
221 
222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
223 
224 bool ClusterSBPoints::HasIn(const SBPoint* pPtToCheck, G4double pMinDist)
225 {
226  // check if the given point is near one of the cluster's point
227  std::set<SBPoint*>::iterator itClusPt;
228  for(itClusPt = fpRegisteredSBPoints.begin();
229  itClusPt != fpRegisteredSBPoints.end();
230  ++itClusPt)
231  {
232  // if are two different pts
233  if((*pPtToCheck != *(*itClusPt)))
234  {
235  // if close enought of an include point of the cluster
236  if( AreOnTheSameCluster(pPtToCheck, *itClusPt, pMinDist) )
237  {
238  return true;
239  }
240  }
241  }
242  return false;
243 }
244 
245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
246 
248  G4double pMinDist)
249 {
250 
251  G4double x1=pPtToCheck->GetPosition().x()/nm;
252  G4double y1=pPtToCheck->GetPosition().y()/nm;
253  G4double z1=pPtToCheck->GetPosition().z()/nm;
254 
255  G4double x2=this->GetBarycenter().x()/nm;
256  G4double y2=this->GetBarycenter().y()/nm;
257  G4double z2=this->GetBarycenter().z()/nm;
258 
259  // if the two points are closed enough
260  if(((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2))<=
261  (pMinDist/nm*pMinDist/nm))
262  {
263  return true;
264  }else
265  {
266  return false;
267  }
268 }
269 
270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
271 
275 {
276  std::set<SBPoint*> points = pCluster->GetRegistredSBPoints();
277  pCluster->Clear();
278  std::set<SBPoint*>::iterator itPt;
279  for(itPt = points.begin(); itPt != points.end(); ++itPt)
280  {
281  this->AddSBPoint(*itPt);
282  }
283 }
284 
bool HasIn(const SBPoint *, G4double)
void AddSBPoint(SBPoint *pSBPoint)
Definition of the ClusterSBPoints class.
std::set< SBPoint * > GetRegistredSBPoints() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetEdep() const
define a cluster of SB Points
ClusterSBPoints(std::set< SBPoint * > pPoints)
double x() const
unsigned int GetSize() const
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
double z() const
void SetCluster(ClusterSBPoints *pCluster)
Definition: SBPoint.hh:74
void FindAllPointsPossible(std::vector< SBPoint * > *ptToCheck, G4int minPts, G4double minDist)
void MergeWith(ClusterSBPoints *)
bool HasInBarycenter(const SBPoint *, G4double)
static constexpr double nm
Definition: G4SIunits.hh:112
double y() const
G4ThreeVector GetBarycenter() const
tuple z
Definition: test.py:28
double G4double
Definition: G4Types.hh:76
G4ThreeVector GetPosition() const
Definition: SBPoint.hh:60
defines a point of energy deposition which defines a damage to the DNA.
Definition: SBPoint.hh:47
bool AreOnTheSameCluster(const SBPoint *pPt1, const SBPoint *pPt2, G4double pMinDist)