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