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