microprox/src/TrackerSD.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 TrackerSD.cc
0027 /// \brief Implementation of the TrackerSD 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 publications:
0032 // Med. Phys. 45 (2018) e722-e739
0033 // Phys. Med. 31 (2015) 861-874
0034 // Med. Phys. 37 (2010) 4692-4708
0035 // Int. J. Model. Simul. Sci. Comput. 1 (2010) 157–178
0036 //
0037 // The Geant4-DNA web site is available at http://geant4-dna.org
0038 //
0039
0040 #include "TrackerSD.hh"
0041
0042 #include "G4AnalysisManager.hh"
0043 #include "G4SDManager.hh"
0044 #include "G4SystemOfUnits.hh"
0045 #include "Randomize.hh"
0046
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0048
0049 TrackerSD::TrackerSD(const G4String& name, const G4String& hitsCollectionName)
0050   : G4VSensitiveDetector(name), fHitsCollection(NULL)
0051 {
0052   collectionName.insert(hitsCollectionName);
0053 }
0054
0055 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0056
0057 TrackerSD::~TrackerSD() {}
0058
0059 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0060
0061 void TrackerSD::Initialize(G4HCofThisEvent* hce)
0062 {
0063   // Create hits collection
0064   fHitsCollection = new TrackerHitsCollection(SensitiveDetectorName, collectionName[0]);
0065
0066   // Add this collection in hce
0067   G4int hcID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
0068
0069   hce->AddHitsCollection(hcID, fHitsCollection);
0070 }
0071
0072 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0073
0074 G4bool TrackerSD::ProcessHits(G4Step* aStep, G4TouchableHistory*)
0075 {
0076   // Energy deposit
0077   G4double edep = aStep->GetTotalEnergyDeposit();
0078
0079   if (edep == 0.) return false;
0080
0081   TrackerHit* newHit = new TrackerHit();
0082
0083   newHit->SetTrackID(aStep->GetTrack()->GetTrackID());
0084   newHit->SetEdep(edep);
0085   newHit->SetPos(aStep->GetPostStepPoint()->GetPosition());
0086
0087   if (aStep->GetTrack()->GetTrackID() == 1 && aStep->GetTrack()->GetParentID() == 0)
0088     newHit->SetIncidentEnergy(aStep->GetTrack()->GetVertexKineticEnergy());
0089
0090   fHitsCollection->insert(newHit);
0091
0092   // Print
0093   // newHit->Print();
0094
0095   return true;
0096 }
0097
0098 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0099 void TrackerSD::EndOfEvent(G4HCofThisEvent*)
0100 {
0101   G4int nofHits = fHitsCollection->entries();
0102
0103   G4double Einc = 0;
0104
0105   /*
0106   G4cout << G4endl
0107   << "-------->Hits Collection: in this event they are "
0108   << nofHits
0109   << " hits in the target volume " << G4endl;
0110   */
0111
0112   // PROCESSING OF PROXIMITY FUNCTION t(x)
0113
0114   // *************************************
0115   // Please select herebelow :
0116   // - the minimum value of x radius
0117   // - the maximum value of x radius
0118   // - the number of steps in radius
0119
0120   G4double minRadius = 0.1 * nm;
0121
0122   G4double maxRadius = 10000 * nm;
0123   G4int nRadiusSteps = 101;
0124
0125   //
0126
0127   auto analysisManager = G4AnalysisManager::Instance();
0128
0129   G4double radius(minRadius);
0130   G4double stpRadius(std::pow(maxRadius / radius, 1. / static_cast<G4double>(nRadiusSteps - 1)));
0131   G4int step(nRadiusSteps);
0132   G4int noRadius(0);
0133
0134   // 1) loop on radius
0135
0136   while (step > 0) {
0137     step--;
0138     noRadius = nRadiusSteps - step;
0139
0140     // G4cout << "---radius/nm=" << radius/nm  << G4endl;
0141
0142     // Computation of t(x)
0143
0144     G4double tNum = 0.;
0145     G4double tDenom = 0.;
0146     G4int nbEdep = 0;
0147
0148     // 2) loop on hits
0149
0150     for (G4int k = 0; k < nofHits; k++) {
0151       G4ThreeVector hitPos = (*fHitsCollection)[k]->GetPos();
0152       G4double hitNrj = (*fHitsCollection)[k]->GetEdep();
0153
0154       // G4cout << "======> hit position x/nm =" << hitPos.x()/nm << G4endl;
0155       // G4cout << "======> hit position y/nm =" << hitPos.y()/nm << G4endl;
0156       // G4cout << "======> hit position z/nm =" << hitPos.z()/nm << G4endl;
0157
0158       // 3) loop on all other hits located within shell
0159
0160       G4double localSum = 0.;
0161
0162       for (G4int i = 0; i < nofHits; i++) {
0163         if ((*fHitsCollection)[i]->GetIncidentEnergy() > 0)
0164           Einc = (*fHitsCollection)[i]->GetIncidentEnergy();
0165
0166         G4ThreeVector localPosi = (*fHitsCollection)[i]->GetPos();
0167
0168         if (((localPosi.x() - hitPos.x()) * (localPosi.x() - hitPos.x())
0169                + (localPosi.y() - hitPos.y()) * (localPosi.y() - hitPos.y())
0170                + (localPosi.z() - hitPos.z()) * (localPosi.z() - hitPos.z())
0171              < radius * stpRadius * radius * stpRadius)
0172             && ((localPosi.x() - hitPos.x()) * (localPosi.x() - hitPos.x())
0173                   + (localPosi.y() - hitPos.y()) * (localPosi.y() - hitPos.y())
0174                   + (localPosi.z() - hitPos.z()) * (localPosi.z() - hitPos.z())
0175                 >= radius * radius))
0176
0177         {
0178           localSum = localSum + (*fHitsCollection)[i]->GetEdep();
0179           nbEdep = nbEdep + 1;
0180         }
0181       }
0182
0183       tNum = tNum + localSum * hitNrj;
0184       tDenom = tDenom + hitNrj;
0185
0186     }  // loop on hits
0187
0188     // fill ntuple including weighting
0189     // does not work with ntuple merging...
0190
0191     analysisManager->FillNtupleDColumn(0, 0, radius / nm);
0192     analysisManager->FillNtupleIColumn(0, 1, noRadius);
0193     analysisManager->FillNtupleDColumn(0, 2, nofHits);
0194     analysisManager->FillNtupleDColumn(0, 3, nbEdep);
0195     analysisManager->FillNtupleDColumn(0, 4, (tNum / tDenom) / eV);
0196     analysisManager->FillNtupleDColumn(0, 5, (stpRadius * radius) / nm);
0197     analysisManager->FillNtupleDColumn(0, 6, Einc / eV);
0198     analysisManager->AddNtupleRow();
0199
0200     // G4cout << "---radius/nm=" << radius/nm << G4endl;
0201     // G4cout << "----end of radius--- " << radius/nm << G4endl;
0202
0203     radius *= stpRadius;
0204
0205   }  // loop on radii
0206 }
0207
0208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......