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