clustering/src/ClusterSBPoints.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 // This example is provided by the Geant4-DNA collaboration
0027 // Any report or published results obtained using the Geant4-DNA software
0028 // shall cite the following Geant4-DNA collaboration publication:
0029 // Med. Phys. 37 (2010) 4692-4708
0030 // The Geant4-DNA web site is available at http://geant4-dna.org
0031 //
0032 // Authors: Henri Payno and Yann Perrot
0033 //
0034 //
0035 /// \file ClusterSBPoints.cc
0036 /// \brief Implementation of the ClusterSBPoints class
0037
0038 #include "ClusterSBPoints.hh"
0039
0040 #include "G4SystemOfUnits.hh"
0041
0042 #include <iostream>
0043
0044 using namespace std;
0045
0046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0047
0048 ClusterSBPoints::ClusterSBPoints(std::set<SBPoint*> pSBPoints) : fpRegisteredSBPoints()
0049 {
0050   UpdateDoubleStrand();
0051   std::set<SBPoint*>::iterator itPt;
0052   for (itPt = pSBPoints.begin(); itPt != pSBPoints.end(); ++itPt) {
0053     AddSBPoint(*itPt);
0054   }
0055 }
0056
0057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0058
0059 ClusterSBPoints::~ClusterSBPoints()
0060 {
0061   Clear();
0062 }
0063
0064 void ClusterSBPoints::Clear()
0065 {
0066   std::set<SBPoint*>::iterator itPt;
0067   for (itPt = fpRegisteredSBPoints.begin(); itPt != fpRegisteredSBPoints.end(); ++itPt) {
0068     (*itPt)->CleanCluster();
0069   }
0070   fpRegisteredSBPoints.clear();
0071 }
0072
0073 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0074
0075 void ClusterSBPoints::AddSBPoint(SBPoint* pSBPoint)
0076 {
0077   assert(pSBPoint);
0078   fpRegisteredSBPoints.insert(pSBPoint);
0079   pSBPoint->SetCluster(this);
0080
0081   UpdateDoubleStrand();
0082 }
0083
0084 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0085
0086 G4ThreeVector ClusterSBPoints::GetBarycenter() const
0087 {
0088   G4double x = 0;
0089   G4double y = 0;
0090   G4double z = 0;
0091
0092   std::set<SBPoint*>::iterator itSDSPt;
0093   for (itSDSPt = fpRegisteredSBPoints.begin(); itSDSPt != fpRegisteredSBPoints.end(); ++itSDSPt) {
0094     x += (*itSDSPt)->GetPosition().x();
0095     y += (*itSDSPt)->GetPosition().y();
0096     z += (*itSDSPt)->GetPosition().z();
0097   }
0098
0099   return G4ThreeVector(x / fpRegisteredSBPoints.size(), y / fpRegisteredSBPoints.size(),
0100                        z / fpRegisteredSBPoints.size());
0101 }
0102
0103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0104
0105 G4double ClusterSBPoints::GetEdep() const
0106 {
0107   G4double res = 0;
0108   std::set<SBPoint*>::iterator itSDSPt;
0109   for (itSDSPt = fpRegisteredSBPoints.begin(); itSDSPt != fpRegisteredSBPoints.end(); ++itSDSPt) {
0110     res += (*itSDSPt)->GetEdep();
0111   }
0112   return res;
0113 }
0114
0115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0116
0117 void ClusterSBPoints::UpdateDoubleStrand()
0118 {
0119   fIsDoubleSB = false;
0120   bool firstStrandTouch = false;
0121   bool secondStrandTouch = false;
0122
0123   std::set<SBPoint*>::iterator itSDSPt;
0124   for (itSDSPt = fpRegisteredSBPoints.begin(); itSDSPt != fpRegisteredSBPoints.end(); ++itSDSPt) {
0125     // if the SDSPoint is localized on the first strand
0126     if (((*itSDSPt)->GetTouchedStrand() == 0) && !firstStrandTouch) {
0127       firstStrandTouch = true;
0128       if (secondStrandTouch) {
0129         fIsDoubleSB = true;
0130         return;
0131       }
0132     }
0133     // if the SDSPoint is localized on the second strand
0134     if (((*itSDSPt)->GetTouchedStrand() == 1) && !secondStrandTouch) {
0135       secondStrandTouch = true;
0136       if (firstStrandTouch) {
0137         fIsDoubleSB = true;
0138         return;
0139       }
0140     }
0141   }
0142 }
0143
0144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0145
0146 bool AreOnTheSameCluster(const SBPoint* pPt1, const SBPoint* pPt2, G4double pMinDist)
0147 {
0148   assert(pPt1);
0149   assert(pPt2);
0150
0151   G4double x1 = pPt1->GetPosition().x() / nm;
0152   G4double y1 = pPt1->GetPosition().y() / nm;
0153   G4double z1 = pPt1->GetPosition().z() / nm;
0154
0155   G4double x2 = pPt2->GetPosition().x() / nm;
0156   G4double y2 = pPt2->GetPosition().y() / nm;
0157   G4double z2 = pPt2->GetPosition().z() / nm;
0158
0159   // if the two points are closed enough
0160   if (((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2))
0161       <= (pMinDist / nm * pMinDist / nm))
0162   {
0163     return true;
0164   }
0165   else {
0166     return false;
0167   }
0168 }
0169
0170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0171
0172 void ClusterSBPoints::FindAllPointsPossible(std::vector<SBPoint*>* pPtsToCheck, G4int pMinPts,
0173                                             G4double pMinDist)
0174 {
0175   assert((unsigned int)pMinPts > this->GetSize());
0176   std::vector<SBPoint*>::iterator itPt = pPtsToCheck->begin();
0177   while (itPt != pPtsToCheck->end()) {
0178     // If 1- each SBpoint is part of only one cluster
0179     //    2- the point isn't already in the cluster
0180     //    3- the point is close enough of the barycenter
0181     if ((!(*itPt)->HasCluster()) && (fpRegisteredSBPoints.find(*itPt) == fpRegisteredSBPoints.end())
0182         && HasInBarycenter(*itPt, pMinDist))  // first version used HasIn method
0183     {
0184       // the point is added
0185       this->AddSBPoint(*itPt);
0186       if (this->GetSize() >= (unsigned int)pMinPts) {
0187         return;
0188       }
0189       // restart from scratch
0190       itPt = pPtsToCheck->begin();
0191     }
0192     else {
0193       ++itPt;
0194     }
0195   }
0196 }
0197
0198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0199
0200 bool ClusterSBPoints::HasIn(const SBPoint* pPtToCheck, G4double pMinDist)
0201 {
0202   // check if the given point is near one of the cluster's point
0203   std::set<SBPoint*>::iterator itClusPt;
0204   for (itClusPt = fpRegisteredSBPoints.begin(); itClusPt != fpRegisteredSBPoints.end(); ++itClusPt)
0205   {
0206     // if are two different pts
0207     if ((*pPtToCheck != *(*itClusPt))) {
0208       // if close enought of an include point of the cluster
0209       if (AreOnTheSameCluster(pPtToCheck, *itClusPt, pMinDist)) {
0210         return true;
0211       }
0212     }
0213   }
0214   return false;
0215 }
0216
0217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0218
0219 bool ClusterSBPoints::HasInBarycenter(const SBPoint* pPtToCheck, G4double pMinDist)
0220 {
0221   G4double x1 = pPtToCheck->GetPosition().x() / nm;
0222   G4double y1 = pPtToCheck->GetPosition().y() / nm;
0223   G4double z1 = pPtToCheck->GetPosition().z() / nm;
0224
0225   G4double x2 = this->GetBarycenter().x() / nm;
0226   G4double y2 = this->GetBarycenter().y() / nm;
0227   G4double z2 = this->GetBarycenter().z() / nm;
0228
0229   // if the two points are closed enough
0230   if (((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2))
0231       <= (pMinDist / nm * pMinDist / nm))
0232   {
0233     return true;
0234   }
0235   else {
0236     return false;
0237   }
0238 }
0239
0240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0241
0242 /// this will insert all registredSBPoint
0243 /// from the given cluster to this cluster.
0244 void ClusterSBPoints::MergeWith(ClusterSBPoints* pCluster)
0245 {
0246   std::set<SBPoint*> points = pCluster->GetRegistredSBPoints();
0247   pCluster->Clear();
0248   std::set<SBPoint*>::iterator itPt;
0249   for (itPt = points.begin(); itPt != points.end(); ++itPt) {
0250     this->AddSBPoint(*itPt);
0251   }
0252 }