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