clustering/src/ClusteringAlgo.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file ClusteringAlgo.cc
0027 /// \brief Implementation of the ClusteringAlgo class
0028
0029 // This example is provided by the Geant4-DNA collaboration
0030 // Any report or published results obtained using the Geant4-DNA software
0031 // shall cite the following Geant4-DNA collaboration publication:
0032 // Med. Phys. 37 (2010) 4692-4708
0033 // The Geant4-DNA web site is available at http://geant4-dna.org
0034 //
0035 // Authors: Henri Payno and Yann Perrot
0036 //
0037 //
0038
0039 #include "ClusteringAlgo.hh"
0040
0041 #include "ClusteringAlgoMessenger.hh"
0042
0043 #include "G4SystemOfUnits.hh"
0044 #include "Randomize.hh"
0045
0046 #include <map>
0047
0048 using namespace std;
0049
0050 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0051
0052 ClusteringAlgo::ClusteringAlgo(G4double pEps, G4int pMinPts, G4double pSPointsProb,
0053                                G4double pEMinDamage, G4double pEMaxDamage)
0054   : fEps(pEps),
0055     fMinPts(pMinPts),
0056     fSPointsProb(pSPointsProb),
0057     fEMinDamage(pEMinDamage),
0058     fEMaxDamage(pEMaxDamage)
0059 {
0060   fNextSBPointID = 0;
0061   fpClustAlgoMessenger = new ClusteringAlgoMessenger(this);
0062 }
0063
0064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0065
0066 ClusteringAlgo::~ClusteringAlgo()
0067 {
0068   delete fpClustAlgoMessenger;
0069   Purge();
0070 }
0071
0072 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0073
0074 // Random sampling in space
0075 G4bool ClusteringAlgo::IsInSensitiveArea()
0076 {
0077   return fSPointsProb > G4UniformRand();
0078 }
0079
0080 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0081
0082 // Random sampling in energy
0083 G4bool ClusteringAlgo::IsEdepSufficient(G4double pEdep)
0084 {
0085   if (pEdep < fEMinDamage) {
0086     return false;
0087   }
0088
0089   else if (pEdep > fEMaxDamage) {
0090     return true;
0091   }
0092   else {
0093     G4double proba = (pEdep / eV - fEMinDamage / eV) / (fEMaxDamage / eV - fEMinDamage / eV);
0094     return (proba > G4UniformRand());
0095   }
0096 }
0097
0098 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0099
0100 // Add an event interaction to the unregistered damage if
0101 // good conditions (pos and energy) are met
0102 //
0103
0104 void ClusteringAlgo::RegisterDamage(G4ThreeVector pPos, G4double pEdep)
0105 {
0106   if (IsEdepSufficient(pEdep)) {
0107     if (IsInSensitiveArea()) {
0108       fpSetOfPoints.push_back(new SBPoint(fNextSBPointID++, pPos, pEdep));
0109     }
0110   }
0111 }
0112
0113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0114
0115 map<G4int, G4int> ClusteringAlgo::RunClustering()
0116 {
0117   // quick sort style
0118   // create cluster
0119   std::vector<SBPoint*>::iterator itVisitorPt, itObservedPt;
0120   for (itVisitorPt = fpSetOfPoints.begin(); itVisitorPt != fpSetOfPoints.end(); ++itVisitorPt) {
0121     itObservedPt = itVisitorPt;
0122     itObservedPt++;
0123     while (itObservedPt != fpSetOfPoints.end()) {
0124       // if at least one of the two points has not a cluster
0125       if (!((*itObservedPt)->HasCluster() && (*itVisitorPt)->HasCluster())) {
0126         if (AreOnTheSameCluster((*itObservedPt)->GetPosition(), (*itVisitorPt)->GetPosition(),
0127                                 fEps))
0128         {
0129           // if none has a cluster. Create a new one
0130           if (!(*itObservedPt)->HasCluster() && !(*itVisitorPt)->HasCluster()) {
0131             // create the new cluster
0132             set<SBPoint*> clusterPoints;
0133             clusterPoints.insert((*itObservedPt));
0134             clusterPoints.insert((*itVisitorPt));
0135             ClusterSBPoints* lCluster = new ClusterSBPoints(clusterPoints);
0136             assert(lCluster);
0137             fpClusters.push_back(lCluster);
0138             assert(lCluster);
0139             // inform SB point that they are part of a cluster now
0140             assert(lCluster);
0141             (*itObservedPt)->SetCluster(lCluster);
0142             assert(lCluster);
0143             (*itVisitorPt)->SetCluster(lCluster);
0144           }
0145           else {
0146             // add the point to the existing cluster
0147             if ((*itObservedPt)->HasCluster()) {
0148               (*itObservedPt)->GetCluster()->AddSBPoint((*itVisitorPt));
0149               (*itVisitorPt)->SetCluster((*itObservedPt)->GetCluster());
0150             }
0151
0152             if ((*itVisitorPt)->HasCluster()) {
0153               (*itVisitorPt)->GetCluster()->AddSBPoint((*itObservedPt));
0154               (*itObservedPt)->SetCluster((*itVisitorPt)->GetCluster());
0155             }
0156           }
0157         }
0158       }
0159       ++itObservedPt;
0160     }
0161   }
0162
0163   // associate isolated points and merge clusters
0164   IncludeUnassociatedPoints();
0165   MergeClusters();
0166
0167   // return cluster size distribution
0168   return GetClusterSizeDistribution();
0169 }
0170
0171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0172
0173 // try to merge cluster between them, based on the distance between barycenters
0174 void ClusteringAlgo::MergeClusters()
0175 {
0176   std::vector<ClusterSBPoints*>::iterator itCluster1, itCluster2;
0177   for (itCluster1 = fpClusters.begin(); itCluster1 != fpClusters.end(); ++itCluster1) {
0178     G4ThreeVector baryCenterClust1 = (*itCluster1)->GetBarycenter();
0179     itCluster2 = itCluster1;
0180     itCluster2++;
0181     while (itCluster2 != fpClusters.end()) {
0182       G4ThreeVector baryCenterClust2 = (*itCluster2)->GetBarycenter();
0183       // if we can merge both cluster
0184       if (AreOnTheSameCluster(baryCenterClust1, baryCenterClust2, fEps)) {
0185         (*itCluster1)->MergeWith(*itCluster2);
0186         delete *itCluster2;
0187         fpClusters.erase(itCluster2);
0188         return MergeClusters();
0189       }
0190       else {
0191         itCluster2++;
0192       }
0193     }
0194   }
0195 }
0196
0197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0198
0199 void ClusteringAlgo::IncludeUnassociatedPoints()
0200 {
0201   std::vector<SBPoint*>::iterator itVisitorPt;
0202   // Associate all point not in a cluster if possible ( to the first found cluster)
0203   for (itVisitorPt = fpSetOfPoints.begin(); itVisitorPt != fpSetOfPoints.end(); ++itVisitorPt) {
0204     if (!(*itVisitorPt)->HasCluster()) {
0205       FindCluster(*itVisitorPt);
0206     }
0207   }
0208 }
0209
0210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0211
0212 bool ClusteringAlgo::FindCluster(SBPoint* pPt)
0213 {
0214   assert(!pPt->HasCluster());
0215   std::vector<ClusterSBPoints*>::iterator itCluster;
0216   for (itCluster = fpClusters.begin(); itCluster != fpClusters.end(); ++itCluster) {
0217     // if((*itCluster)->hasIn(pPt, fEps))
0218     if ((*itCluster)->HasInBarycenter(pPt, fEps)) {
0219       (*itCluster)->AddSBPoint(pPt);
0220       return true;
0221     }
0222   }
0223   return false;
0224 }
0225
0226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0227
0228 bool ClusteringAlgo::AreOnTheSameCluster(G4ThreeVector pPt1, G4ThreeVector pPt2, G4double pMinDist)
0229 {
0230   G4double x1 = pPt1.x() / nm;
0231   G4double y1 = pPt1.y() / nm;
0232   G4double z1 = pPt1.z() / nm;
0233
0234   G4double x2 = pPt2.x() / nm;
0235   G4double y2 = pPt2.y() / nm;
0236   G4double z2 = pPt2.z() / nm;
0237
0238   // if the two points are closed enough
0239   if (((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2))
0240       <= (pMinDist / nm * pMinDist / nm))
0241   {
0242     return true;
0243   }
0244   else {
0245     return false;
0246   }
0247 }
0248
0249 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0250
0251 G4int ClusteringAlgo::GetSSB() const
0252 {
0253   G4int nbSSB = 0;
0254   std::vector<SBPoint*>::const_iterator itSDSPt;
0255   for (itSDSPt = fpSetOfPoints.begin(); itSDSPt != fpSetOfPoints.end(); ++itSDSPt) {
0256     if (!(*itSDSPt)->HasCluster()) {
0257       nbSSB++;
0258     }
0259   }
0260   return nbSSB;
0261 }
0262
0263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0264
0265 G4int ClusteringAlgo::GetComplexSSB() const
0266 {
0267   G4int nbSSB = 0;
0268   std::vector<ClusterSBPoints*>::const_iterator itCluster;
0269   for (itCluster = fpClusters.begin(); itCluster != fpClusters.end(); ++itCluster) {
0270     if ((*itCluster)->IsSSB()) {
0271       nbSSB++;
0272     }
0273   }
0274   return nbSSB;
0275 }
0276
0277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0278
0279 G4int ClusteringAlgo::GetDSB() const
0280 {
0281   G4int nbDSB = 0;
0282   std::vector<ClusterSBPoints*>::const_iterator itCluster;
0283   for (itCluster = fpClusters.begin(); itCluster != fpClusters.end(); ++itCluster) {
0284     if ((*itCluster)->IsDSB()) {
0285       nbDSB++;
0286     }
0287   }
0288   return nbDSB;
0289 }
0290
0291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0292
0293 map<G4int, G4int> ClusteringAlgo::GetClusterSizeDistribution()
0294 {
0295   std::map<G4int, G4int> sizeDistribution;
0296   sizeDistribution[1] = GetSSB();
0297   std::vector<ClusterSBPoints*>::const_iterator itCluster;
0298   for (itCluster = fpClusters.begin(); itCluster != fpClusters.end(); itCluster++) {
0299     sizeDistribution[(*itCluster)->GetSize()]++;
0300   }
0301   return sizeDistribution;
0302 }
0303
0304 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0305
0306 void ClusteringAlgo::Purge()
0307 {
0308   fNextSBPointID = 0;
0309   std::vector<ClusterSBPoints*>::iterator itCluster;
0310   for (itCluster = fpClusters.begin(); itCluster != fpClusters.end(); ++itCluster) {
0311     delete *itCluster;
0312     *itCluster = NULL;
0313   }
0314   fpClusters.clear();
0315   std::vector<SBPoint*>::iterator itPt;
0316   for (itPt = fpSetOfPoints.begin(); itPt != fpSetOfPoints.end(); ++itPt) {
0317     delete *itPt;
0318     *itPt = NULL;
0319   }
0320   fpSetOfPoints.clear();
0321 }