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 ClusteringAlgo.cc 36 /// \file ClusteringAlgo.cc 36 /// \brief Implementation of the ClustreringAl 37 /// \brief Implementation of the ClustreringAlgo class 37 38 38 #include "ClusteringAlgo.hh" 39 #include "ClusteringAlgo.hh" 39 << 40 #include "ClusteringAlgoMessenger.hh" 40 #include "ClusteringAlgoMessenger.hh" 41 41 42 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh" 43 #include "Randomize.hh" 43 #include "Randomize.hh" 44 44 45 #include <map> 45 #include <map> 46 46 47 using namespace std; 47 using namespace std; 48 48 49 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 50 50 51 ClusteringAlgo::ClusteringAlgo(G4double pEps, << 51 ClusteringAlgo::ClusteringAlgo(G4double pEps,G4int pMinPts, 52 G4double pEMinD << 52 G4double pSPointsProb,G4double pEMinDamage, G4double pEMaxDamage) 53 : fEps(pEps), << 53 :fEps(pEps),fMinPts(pMinPts), 54 fMinPts(pMinPts), << 54 fSPointsProb(pSPointsProb),fEMinDamage(pEMinDamage),fEMaxDamage(pEMaxDamage) 55 fSPointsProb(pSPointsProb), << 56 fEMinDamage(pEMinDamage), << 57 fEMaxDamage(pEMaxDamage) << 58 { 55 { 59 fNextSBPointID = 0; 56 fNextSBPointID = 0; 60 fpClustAlgoMessenger = new ClusteringAlgoMes 57 fpClustAlgoMessenger = new ClusteringAlgoMessenger(this); 61 } 58 } 62 59 63 //....oooOO0OOooo........oooOO0OOooo........oo 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 64 61 65 ClusteringAlgo::~ClusteringAlgo() 62 ClusteringAlgo::~ClusteringAlgo() 66 { 63 { 67 delete fpClustAlgoMessenger; 64 delete fpClustAlgoMessenger; 68 Purge(); 65 Purge(); 69 } 66 } 70 67 71 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 72 69 73 // Random sampling in space 70 // Random sampling in space 74 G4bool ClusteringAlgo::IsInSensitiveArea() 71 G4bool ClusteringAlgo::IsInSensitiveArea() 75 { 72 { 76 return fSPointsProb > G4UniformRand(); 73 return fSPointsProb > G4UniformRand(); 77 } 74 } 78 75 79 //....oooOO0OOooo........oooOO0OOooo........oo 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 80 77 81 // Random sampling in energy 78 // Random sampling in energy 82 G4bool ClusteringAlgo::IsEdepSufficient(G4doub 79 G4bool ClusteringAlgo::IsEdepSufficient(G4double pEdep) 83 { 80 { 84 if (pEdep < fEMinDamage) { << 81 if(pEdep<fEMinDamage) >> 82 { 85 return false; 83 return false; 86 } 84 } 87 85 88 else if (pEdep > fEMaxDamage) { << 86 else if(pEdep>fEMaxDamage) >> 87 { 89 return true; 88 return true; 90 } 89 } 91 else { << 90 else 92 G4double proba = (pEdep / eV - fEMinDamage << 91 { 93 return (proba > G4UniformRand()); << 92 G4double proba = (pEdep/eV - fEMinDamage/eV)/ >> 93 (fEMaxDamage/eV-fEMinDamage/eV); >> 94 return (proba>G4UniformRand()); 94 } 95 } >> 96 95 } 97 } 96 98 97 //....oooOO0OOooo........oooOO0OOooo........oo 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 98 100 99 // Add an event interaction to the unregistere 101 // Add an event interaction to the unregistered damage if 100 // good conditions (pos and energy) are met 102 // good conditions (pos and energy) are met 101 // 103 // 102 104 103 void ClusteringAlgo::RegisterDamage(G4ThreeVec << 105 void ClusteringAlgo::RegisterDamage(G4ThreeVector pPos,G4double pEdep) 104 { 106 { 105 if (IsEdepSufficient(pEdep)) { << 107 if(IsEdepSufficient(pEdep)) 106 if (IsInSensitiveArea()) { << 108 { 107 fpSetOfPoints.push_back(new SBPoint(fNex << 109 if(IsInSensitiveArea()) >> 110 { >> 111 fpSetOfPoints.push_back( new SBPoint(fNextSBPointID++, pPos, pEdep)); 108 } 112 } 109 } 113 } 110 } 114 } 111 115 112 //....oooOO0OOooo........oooOO0OOooo........oo 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 113 117 114 map<G4int, G4int> ClusteringAlgo::RunClusterin << 118 map<G4int,G4int> ClusteringAlgo::RunClustering() 115 { 119 { >> 120 116 // quick sort style 121 // quick sort style 117 // create cluster 122 // create cluster 118 std::vector<SBPoint*>::iterator itVisitorPt, 123 std::vector<SBPoint*>::iterator itVisitorPt, itObservedPt; 119 for (itVisitorPt = fpSetOfPoints.begin(); it << 124 for(itVisitorPt = fpSetOfPoints.begin(); >> 125 itVisitorPt != fpSetOfPoints.end(); >> 126 ++itVisitorPt ) >> 127 { 120 itObservedPt = itVisitorPt; 128 itObservedPt = itVisitorPt; 121 itObservedPt++; << 129 itObservedPt ++; 122 while (itObservedPt != fpSetOfPoints.end() << 130 while(itObservedPt != fpSetOfPoints.end() ) >> 131 { 123 // if at least one of the two points has 132 // if at least one of the two points has not a cluster 124 if (!((*itObservedPt)->HasCluster() && ( << 133 if(!((*itObservedPt)->HasCluster() && (*itVisitorPt)->HasCluster())) 125 if (AreOnTheSameCluster((*itObservedPt << 134 { 126 fEps)) << 135 if(AreOnTheSameCluster( (*itObservedPt)->GetPosition(), >> 136 (*itVisitorPt)->GetPosition(),fEps)) 127 { 137 { 128 // if none has a cluster. Create a n 138 // if none has a cluster. Create a new one 129 if (!(*itObservedPt)->HasCluster() & << 139 if(!(*itObservedPt)->HasCluster() && !(*itVisitorPt)->HasCluster()) >> 140 { 130 // create the new cluster 141 // create the new cluster 131 set<SBPoint*> clusterPoints; 142 set<SBPoint*> clusterPoints; 132 clusterPoints.insert((*itObservedP 143 clusterPoints.insert((*itObservedPt)); 133 clusterPoints.insert((*itVisitorPt 144 clusterPoints.insert((*itVisitorPt)); 134 ClusterSBPoints* lCluster = new Cl 145 ClusterSBPoints* lCluster = new ClusterSBPoints(clusterPoints); 135 assert(lCluster); 146 assert(lCluster); 136 fpClusters.push_back(lCluster); 147 fpClusters.push_back(lCluster); 137 assert(lCluster); 148 assert(lCluster); 138 // inform SB point that they are p 149 // inform SB point that they are part of a cluster now 139 assert(lCluster); 150 assert(lCluster); 140 (*itObservedPt)->SetCluster(lClust 151 (*itObservedPt)->SetCluster(lCluster); 141 assert(lCluster); 152 assert(lCluster); 142 (*itVisitorPt)->SetCluster(lCluste 153 (*itVisitorPt)->SetCluster(lCluster); 143 } << 154 }else 144 else { << 155 { 145 // add the point to the existing c 156 // add the point to the existing cluster 146 if ((*itObservedPt)->HasCluster()) << 157 if((*itObservedPt)->HasCluster()) >> 158 { 147 (*itObservedPt)->GetCluster()->A 159 (*itObservedPt)->GetCluster()->AddSBPoint((*itVisitorPt)); 148 (*itVisitorPt)->SetCluster((*itO 160 (*itVisitorPt)->SetCluster((*itObservedPt)->GetCluster()); 149 } 161 } 150 162 151 if ((*itVisitorPt)->HasCluster()) << 163 if((*itVisitorPt)->HasCluster()) >> 164 { 152 (*itVisitorPt)->GetCluster()->Ad 165 (*itVisitorPt)->GetCluster()->AddSBPoint((*itObservedPt)); 153 (*itObservedPt)->SetCluster((*it 166 (*itObservedPt)->SetCluster((*itVisitorPt)->GetCluster()); 154 } 167 } 155 } 168 } 156 } 169 } 157 } 170 } 158 ++itObservedPt; 171 ++itObservedPt; 159 } 172 } 160 } 173 } 161 174 162 // associate isolated points and merge clust 175 // associate isolated points and merge clusters 163 IncludeUnassociatedPoints(); 176 IncludeUnassociatedPoints(); 164 MergeClusters(); 177 MergeClusters(); 165 178 166 // return cluster size distribution 179 // return cluster size distribution 167 return GetClusterSizeDistribution(); 180 return GetClusterSizeDistribution(); 168 } 181 } 169 182 170 //....oooOO0OOooo........oooOO0OOooo........oo 183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 171 184 172 // try to merge cluster between them, based on 185 // try to merge cluster between them, based on the distance between barycenters 173 void ClusteringAlgo::MergeClusters() 186 void ClusteringAlgo::MergeClusters() 174 { 187 { 175 std::vector<ClusterSBPoints*>::iterator itCl 188 std::vector<ClusterSBPoints*>::iterator itCluster1, itCluster2; 176 for (itCluster1 = fpClusters.begin(); itClus << 189 for(itCluster1 = fpClusters.begin(); >> 190 itCluster1 != fpClusters.end(); >> 191 ++itCluster1) >> 192 { 177 G4ThreeVector baryCenterClust1 = (*itClust 193 G4ThreeVector baryCenterClust1 = (*itCluster1)->GetBarycenter(); 178 itCluster2 = itCluster1; 194 itCluster2 = itCluster1; 179 itCluster2++; 195 itCluster2++; 180 while (itCluster2 != fpClusters.end()) { << 196 while(itCluster2 != fpClusters.end()) >> 197 { 181 G4ThreeVector baryCenterClust2 = (*itClu 198 G4ThreeVector baryCenterClust2 = (*itCluster2)->GetBarycenter(); 182 // if we can merge both cluster 199 // if we can merge both cluster 183 if (AreOnTheSameCluster(baryCenterClust1 << 200 if(AreOnTheSameCluster(baryCenterClust1, baryCenterClust2,fEps)) >> 201 { 184 (*itCluster1)->MergeWith(*itCluster2); 202 (*itCluster1)->MergeWith(*itCluster2); 185 delete *itCluster2; 203 delete *itCluster2; 186 fpClusters.erase(itCluster2); 204 fpClusters.erase(itCluster2); 187 return MergeClusters(); 205 return MergeClusters(); 188 } << 206 }else 189 else { << 207 { 190 itCluster2++; 208 itCluster2++; 191 } 209 } 192 } 210 } 193 } 211 } 194 } 212 } 195 213 196 //....oooOO0OOooo........oooOO0OOooo........oo 214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 197 215 198 void ClusteringAlgo::IncludeUnassociatedPoints 216 void ClusteringAlgo::IncludeUnassociatedPoints() 199 { 217 { 200 std::vector<SBPoint*>::iterator itVisitorPt; 218 std::vector<SBPoint*>::iterator itVisitorPt; >> 219 int nbPtSansCluster = 0; 201 // Associate all point not in a cluster if p 220 // Associate all point not in a cluster if possible ( to the first found cluster) 202 for (itVisitorPt = fpSetOfPoints.begin(); it << 221 for(itVisitorPt = fpSetOfPoints.begin(); 203 if (!(*itVisitorPt)->HasCluster()) { << 222 itVisitorPt != fpSetOfPoints.end(); >> 223 ++itVisitorPt) >> 224 { >> 225 if(!(*itVisitorPt)->HasCluster()) >> 226 { >> 227 nbPtSansCluster ++; 204 FindCluster(*itVisitorPt); 228 FindCluster(*itVisitorPt); 205 } 229 } 206 } 230 } 207 } 231 } 208 232 209 //....oooOO0OOooo........oooOO0OOooo........oo 233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 210 234 211 bool ClusteringAlgo::FindCluster(SBPoint* pPt) 235 bool ClusteringAlgo::FindCluster(SBPoint* pPt) 212 { 236 { 213 assert(!pPt->HasCluster()); 237 assert(!pPt->HasCluster()); 214 std::vector<ClusterSBPoints*>::iterator itCl 238 std::vector<ClusterSBPoints*>::iterator itCluster; 215 for (itCluster = fpClusters.begin(); itClust << 239 for(itCluster = fpClusters.begin(); 216 // if((*itCluster)->hasIn(pPt, fEps)) << 240 itCluster != fpClusters.end(); 217 if ((*itCluster)->HasInBarycenter(pPt, fEp << 241 ++itCluster) >> 242 { >> 243 //if((*itCluster)->hasIn(pPt, fEps)) >> 244 if((*itCluster)->HasInBarycenter(pPt, fEps)) >> 245 { 218 (*itCluster)->AddSBPoint(pPt); 246 (*itCluster)->AddSBPoint(pPt); 219 return true; 247 return true; 220 } 248 } 221 } << 249 } 222 return false; 250 return false; 223 } 251 } 224 252 225 //....oooOO0OOooo........oooOO0OOooo........oo 253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 226 254 227 bool ClusteringAlgo::AreOnTheSameCluster(G4Thr << 255 bool ClusteringAlgo::AreOnTheSameCluster(G4ThreeVector pPt1, >> 256 G4ThreeVector pPt2, G4double pMinDist) 228 { 257 { 229 G4double x1 = pPt1.x() / nm; << 258 G4double x1=pPt1.x()/nm; 230 G4double y1 = pPt1.y() / nm; << 259 G4double y1=pPt1.y()/nm; 231 G4double z1 = pPt1.z() / nm; << 260 G4double z1=pPt1.z()/nm; 232 << 261 233 G4double x2 = pPt2.x() / nm; << 262 G4double x2=pPt2.x()/nm; 234 G4double y2 = pPt2.y() / nm; << 263 G4double y2=pPt2.y()/nm; 235 G4double z2 = pPt2.z() / nm; << 264 G4double z2=pPt2.z()/nm; 236 265 237 // if the two points are closed enough 266 // if the two points are closed enough 238 if (((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 << 267 if(((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2))<= 239 <= (pMinDist / nm * pMinDist / nm)) << 268 (pMinDist/nm*pMinDist/nm)) 240 { 269 { 241 return true; 270 return true; 242 } << 271 }else 243 else { << 272 { 244 return false; 273 return false; 245 } 274 } 246 } 275 } 247 276 248 //....oooOO0OOooo........oooOO0OOooo........oo 277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 249 278 250 G4int ClusteringAlgo::GetSSB() const 279 G4int ClusteringAlgo::GetSSB() const 251 { 280 { 252 G4int nbSSB = 0; 281 G4int nbSSB = 0; 253 std::vector<SBPoint*>::const_iterator itSDSP 282 std::vector<SBPoint*>::const_iterator itSDSPt; 254 for (itSDSPt = fpSetOfPoints.begin(); itSDSP << 283 for(itSDSPt = fpSetOfPoints.begin(); 255 if (!(*itSDSPt)->HasCluster()) { << 284 itSDSPt != fpSetOfPoints.end(); >> 285 ++itSDSPt) >> 286 { >> 287 if(!(*itSDSPt)->HasCluster()) >> 288 { 256 nbSSB++; 289 nbSSB++; 257 } 290 } 258 } 291 } 259 return nbSSB; 292 return nbSSB; 260 } 293 } 261 294 262 //....oooOO0OOooo........oooOO0OOooo........oo 295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 263 296 264 G4int ClusteringAlgo::GetComplexSSB() const 297 G4int ClusteringAlgo::GetComplexSSB() const 265 { 298 { 266 G4int nbSSB = 0; 299 G4int nbSSB = 0; 267 std::vector<ClusterSBPoints*>::const_iterato 300 std::vector<ClusterSBPoints*>::const_iterator itCluster; 268 for (itCluster = fpClusters.begin(); itClust << 301 for(itCluster = fpClusters.begin(); 269 if ((*itCluster)->IsSSB()) { << 302 itCluster != fpClusters.end(); 270 nbSSB++; << 303 ++itCluster) >> 304 { >> 305 if((*itCluster)->IsSSB()) >> 306 { >> 307 nbSSB ++; 271 } 308 } 272 } 309 } 273 return nbSSB; 310 return nbSSB; 274 } 311 } 275 312 276 //....oooOO0OOooo........oooOO0OOooo........oo 313 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 277 314 278 G4int ClusteringAlgo::GetDSB() const 315 G4int ClusteringAlgo::GetDSB() const 279 { 316 { 280 G4int nbDSB = 0; 317 G4int nbDSB = 0; 281 std::vector<ClusterSBPoints*>::const_iterato 318 std::vector<ClusterSBPoints*>::const_iterator itCluster; 282 for (itCluster = fpClusters.begin(); itClust << 319 for(itCluster = fpClusters.begin(); 283 if ((*itCluster)->IsDSB()) { << 320 itCluster != fpClusters.end(); 284 nbDSB++; << 321 ++itCluster) >> 322 { >> 323 if((*itCluster)->IsDSB()) >> 324 { >> 325 nbDSB ++; 285 } 326 } 286 } 327 } 287 return nbDSB; 328 return nbDSB; 288 } 329 } 289 330 290 //....oooOO0OOooo........oooOO0OOooo........oo 331 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 291 332 292 map<G4int, G4int> ClusteringAlgo::GetClusterSi << 333 map<G4int,G4int> ClusteringAlgo::GetClusterSizeDistribution() 293 { 334 { 294 std::map<G4int, G4int> sizeDistribution; << 335 std::map<G4int,G4int> sizeDistribution; 295 sizeDistribution[1] = GetSSB(); << 336 sizeDistribution[1]=GetSSB(); 296 std::vector<ClusterSBPoints*>::const_iterato 337 std::vector<ClusterSBPoints*>::const_iterator itCluster; 297 for (itCluster = fpClusters.begin(); itClust << 338 for(itCluster = fpClusters.begin(); >> 339 itCluster != fpClusters.end(); >> 340 itCluster++) >> 341 { 298 sizeDistribution[(*itCluster)->GetSize()]+ 342 sizeDistribution[(*itCluster)->GetSize()]++; 299 } 343 } 300 return sizeDistribution; 344 return sizeDistribution; 301 } 345 } 302 346 303 //....oooOO0OOooo........oooOO0OOooo........oo 347 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 304 348 305 void ClusteringAlgo::Purge() 349 void ClusteringAlgo::Purge() 306 { 350 { 307 fNextSBPointID = 0; 351 fNextSBPointID = 0; 308 std::vector<ClusterSBPoints*>::iterator itCl 352 std::vector<ClusterSBPoints*>::iterator itCluster; 309 for (itCluster = fpClusters.begin(); itClust << 353 for(itCluster = fpClusters.begin(); >> 354 itCluster != fpClusters.end(); >> 355 ++itCluster) >> 356 { 310 delete *itCluster; 357 delete *itCluster; 311 *itCluster = NULL; 358 *itCluster = NULL; 312 } << 359 } 313 fpClusters.clear(); 360 fpClusters.clear(); 314 std::vector<SBPoint*>::iterator itPt; 361 std::vector<SBPoint*>::iterator itPt; 315 for (itPt = fpSetOfPoints.begin(); itPt != f << 362 for(itPt = fpSetOfPoints.begin(); >> 363 itPt != fpSetOfPoints.end(); >> 364 ++itPt) >> 365 { 316 delete *itPt; 366 delete *itPt; 317 *itPt = NULL; 367 *itPt = NULL; 318 } 368 } 319 fpSetOfPoints.clear(); 369 fpSetOfPoints.clear(); 320 } 370 } >> 371 321 372