microtrack/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 #include "TrackerSD.hh"
0030
0031 #include "G4AnalysisManager.hh"
0032 #include "G4Exception.hh"
0033 #include "G4ParticleDefinition.hh"
0034 #include "G4ProcessType.hh"
0035 #include "G4RunManager.hh"
0036 #include "G4SDManager.hh"
0037 #include "G4Step.hh"
0038 #include "G4StepPoint.hh"
0039 #include "G4SystemOfUnits.hh"
0040 #include "G4VProcess.hh"
0041 #include "Randomize.hh"
0042 #include "globals.hh"
0043
0044 #include "DetectorConstruction.hh"
0045
0046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0047 TrackerSD::TrackerSD(const G4String& sdName, const G4String& hitsCollectionName,
0048                      const int /*depthIndex*/)
0049   : G4VSensitiveDetector(sdName)
0050 {
0051   collectionName.insert(hitsCollectionName);
0052 }
0053
0054 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0055
0056 TrackerSD::~TrackerSD() = default;
0057
0058 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0059
0060 void TrackerSD::Initialize(G4HCofThisEvent* hce)
0061 {
0062   // Create hits collection
0063   fHitsCollection =
0064     new TrackerHitColl(SensitiveDetectorName, collectionName[0]);
0065
0066   // Add this collection in hce
0067   G4int collID =
0068     G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
0069
0070   hce->AddHitsCollection(collID, fHitsCollection);
0071 }
0072
0073 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0074
0075 G4bool TrackerSD::ProcessHits(G4Step* aStep, G4TouchableHistory*)
0076 {
0077   // Get the track
0078   G4Track* aTrack = aStep->GetTrack();
0079
0080   // Get Pre and Post step points
0081   G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0082   G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
0083
0084   // Get ParentID
0085   G4int parentID = aTrack->GetParentID();
0086
0087   // Only for primary particles
0088   if (parentID == 0) {
0089     // Check if the particle is entering or exiting
0090     if (preStepPoint->GetStepStatus() == fGeomBoundary) {
0091       // Get the pre-step energy
0092       G4double preEnergy = preStepPoint->GetKineticEnergy();
0093
0094       // Fill the histogram
0095       G4AnalysisManager::Instance()->FillH1(12, preEnergy);
0096     }
0097
0098     if (postStepPoint->GetStepStatus() == fGeomBoundary) {
0099       // Get the post-step energy
0100       G4double postEnergy = postStepPoint->GetKineticEnergy();
0101
0102       // Fill the histogram
0103       G4AnalysisManager::Instance()->FillH1(13, postEnergy);
0104     }
0105   }
0106
0107   // Get the energy deposit
0108   G4double edep = aStep->GetTotalEnergyDeposit();
0109   if (edep == 0.) return false;
0110
0111   // Declare a new hit
0112   auto newHit = new TrackerHit();
0113
0114   newHit->SetEdep(edep);
0115
0116   // Get the position of the hit
0117   G4ThreeVector Pos = postStepPoint->GetPosition();
0118   newHit->SetPosition(Pos);
0119
0120   // Add the hit to the collection
0121   fHitsCollection->insert(newHit);
0122
0123   return true;
0124 }
0125
0126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0127
0128 void TrackerSD::EndOfEvent(G4HCofThisEvent*)
0129 {
0130   const auto nofHits = fHitsCollection->entries();
0131
0132   if (verboseLevel > 0) {
0133     G4cout << "  Hits Collection: in this event, we have " << nofHits
0134            << " hits in the tracker." << G4endl;
0135   }
0136
0137   // Get the DetectorConstruction instance
0138   const DetectorConstruction* detConstruction =
0139     static_cast<const DetectorConstruction*>(
0140       G4RunManager::GetRunManager()->GetUserDetectorConstruction());
0141
0142   // Declare variables for the dimensions of the detector
0143   G4double HitSelRegZ = 0.;
0144   G4double HitSelRegXY = 0.;
0145   G4double site_radius = 0.;
0146   G4double DetDensity = 0.;
0147   G4ThreeVector hitPos;
0148   G4ThreeVector randCenterPos;
0149
0150   if (detConstruction) {
0151     site_radius = detConstruction->GetSiteRadius();
0152     HitSelRegZ = detConstruction->GetHitSelRegZ();
0153     HitSelRegXY = detConstruction->GetHitSelRegXY();
0154     auto mat = G4Material::GetMaterial(detConstruction->GetMaterial());
0155     if (mat)
0156       DetDensity = mat->GetDensity();
0157     else {
0158       G4ExceptionDescription msg;
0159       msg << "Material not found."
0160           << "Something unexpected has occurred." << G4endl;
0161       G4Exception("TrackerSD::EndOfEvent()", "TrackerSD002", JustWarning, msg);
0162     }
0163   }
0164   else {
0165     G4ExceptionDescription msg;
0166     msg << "Detector construction not found."
0167         << "Something unexpected has occurred." << G4endl;
0168     G4Exception("TrackerSD::EndOfEvent()", "TrackerSD001", JustWarning, msg);
0169   }
0170   G4int nHsel = 0;  // number of hits in the selection region
0171   G4int nHsite = 0;  // number of hits in the site
0172   G4int nHint = 0;  // number of hits both in site and selection region
0173   G4double evtEdep = 0.;  // energy deposit in the event
0174
0175   if (nofHits > 0) {
0176     const G4int maxTries = 1000;
0177     G4int tries = 0;
0178     G4bool found = false;
0179
0180     while (tries < maxTries && !found) {
0181       std::size_t randHit = static_cast<std::size_t>(G4UniformRand() * nofHits);
0182       auto hit = (*fHitsCollection)[randHit];
0183       G4ThreeVector hitPosition = hit->GetPosition();
0184
0185       if (hitPosition.x() > -HitSelRegXY / 2
0186           && hitPosition.x() < HitSelRegXY / 2
0187           && hitPosition.y() > -HitSelRegXY / 2
0188           && hitPosition.y() < HitSelRegXY / 2
0189           && hitPosition.z() > -HitSelRegZ / 2
0190           && hitPosition.z() < HitSelRegZ / 2)
0191       {
0192         hitPos = hitPosition;  // store valid hit position
0193
0194         found = true;
0195       }
0196       ++tries;
0197     }
0198
0199     if (found) {
0200       G4double site_radius2 = site_radius * site_radius;
0201
0202       // Random placement of a sphere center containing the hit
0203       G4double xRand, yRand, zRand, randRad2;
0204       do {
0205         xRand = (2 * G4UniformRand() - 1) * site_radius;
0206         yRand = (2 * G4UniformRand() - 1) * site_radius;
0207         zRand = (2 * G4UniformRand() - 1) * site_radius;
0208         randRad2 = xRand * xRand + yRand * yRand + zRand * zRand;
0209       } while (randRad2 > site_radius2);
0210
0211       randCenterPos = G4ThreeVector(xRand + hitPos.x(), yRand + hitPos.y(),
0212                                     zRand + hitPos.z());
0213     }
0214
0215     else {
0216       G4cout
0217         << "In this event, no hits were found in the hit selection "
0218         << "region after trying " << tries << " times." << G4endl
0219         << "Please consider increasing the size of the hit selection region."
0220         << G4endl;
0221       return;  // Skip analysis if no valid hit was found
0222     }
0223   }
0224
0225   // Loop over hits
0226   for (std::size_t jj = 0; jj < nofHits; jj++) {
0227     auto hit2 = (*fHitsCollection)[jj];
0228     G4ThreeVector hit2Position = hit2->GetPosition();
0229
0230     G4bool inSlab =
0231       (hit2Position.x() > -HitSelRegXY / 2 && hit2Position.x() < HitSelRegXY / 2
0232        && hit2Position.y() > -HitSelRegXY / 2
0233        && hit2Position.y() < HitSelRegXY / 2
0234        && hit2Position.z() > -HitSelRegZ / 2
0235        && hit2Position.z() < HitSelRegZ / 2);
0236
0237     // Check if the hit is within the site
0238     G4double dist = (hit2Position - randCenterPos).mag();
0239     if (dist < site_radius) {
0240       nHsite++;
0241       evtEdep += hit2->GetEdep();
0242     }
0243
0244     // Check if the hit is within the selection region
0245     if (inSlab) nHsel++;
0246
0247     // Check if the hit is within the selection region and the site
0248     if (dist < site_radius && inSlab) nHint++;
0249   }
0250
0251   // Access the analysis manager
0252   auto analysisManager = G4AnalysisManager::Instance();
0253
0254   // Calculate lineal energy
0255   G4double y = (evtEdep) / ((4. / 3.) * site_radius);  // in keV/um
0256
0257   // Calculate specific energy
0258   G4double mass = ((4. / 3.) * CLHEP::pi * site_radius * site_radius
0259                    * site_radius * DetDensity);
0260   G4double z = (evtEdep / mass);
0261
0262   // Fill histograms
0263   if (nHint > 0) {
0264     // Define the weight
0265     G4double wght = G4double(nHsel) / G4double(nHint);
0266
0267     // Histogram 0: Single-event energy imparted in keV
0268     analysisManager->FillH1(0, evtEdep / keV, wght);
0269
0270     // Histogram 1: Weighted single-event energy imparted in keV
0271     analysisManager->FillH1(1, evtEdep / keV, (evtEdep * wght / keV));
0272
0273     // Histogram 2: Squared weighted energy imparted per event in keV^2
0274     analysisManager->FillH1(2, evtEdep / keV,
0275                             ((evtEdep * evtEdep * wght) / (keV * keV)));
0276
0277     // Histogram 3: Lineal energy in keV/um
0278     analysisManager->FillH1(3, y / (keV / um), wght);
0279
0280     // Histogram 4: Dose-weighted lineal energy in keV/um
0281     analysisManager->FillH1(4, y / (keV / um), (y * wght / (keV / um)));
0282
0283     // Histogram 5: Squared weighted lineal energy in (keV/um)^2
0284     analysisManager->FillH1(5, y / (keV / um),
0285                             (y * y * wght / ((keV / um) * (keV / um))));
0286
0287     // Histogram 6: Specific energy in keV/g
0288     analysisManager->FillH1(6, z / gray, wght);
0289
0290     // Histogram 7: Weighted specific energy in keV/g
0291     analysisManager->FillH1(7, z / gray, (z * wght / gray));
0292
0293     // Histogram 8: Squared-weighted specific energy in (keV/g)^2
0294     analysisManager->FillH1(8, z / gray, (z * z * wght / (gray * gray)));
0295
0296     // Histogram 9: Number of hits in the selection region
0297     analysisManager->FillH1(9, nHsel, 1.);
0298
0299     // Histogram 10: Number of hits in the site
0300     analysisManager->FillH1(10, nHsite, 1.);
0301
0302     // Histogram 11: Number of hits both in site and selection region
0303     analysisManager->FillH1(11, nHint, 1.);
0304
0305     // Histogram 14: Number of hits in site vs Energy imparted per event
0306     analysisManager->FillH2(0, evtEdep / keV, nHsite, 1.);
0307   }
0308   else {
0309     G4ExceptionDescription msg;
0310     msg << "In this event, we had nHint == 0. "
0311         << "Something unexpected has occurred." << G4endl;
0312     G4Exception("TrackerSD::EndOfEvent()", "TrackerSD003", JustWarning, msg);
0313   }
0314 }