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