radiobiology/src/LETAccumulable.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 LETAccumulable.cc
0027 /// \brief Implementation of the RadioBio::LETAccumulable class
0028
0029 #include "LETAccumulable.hh"
0030
0031 #include "G4EmCalculator.hh"
0032 #include "G4LogicalVolume.hh"
0033 #include "G4ParticleDefinition.hh"
0034
0035 #include "Hit.hh"
0036 #include "Manager.hh"
0037 #include "VoxelizedSensitiveDetector.hh"
0038 #include <G4SystemOfUnits.hh>
0039
0040 #include <tuple>
0041
0042 namespace RadioBio
0043 {
0044
0045 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0046
0047 LETAccumulable::LETAccumulable() : VRadiobiologicalAccumulable("LET") {}
0048
0049 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0050
0051 void LETAccumulable::Merge(const G4VAccumulable& rhs)
0052 {
0053   if (GetVerboseLevel() > 0) {
0054     G4cout << "LETAccumulable::Merge()" << G4endl;
0055   }
0056   const LETAccumulable& other = dynamic_cast<const LETAccumulable&>(rhs);
0057
0058   // Merges standard counters
0059   fTotalLETT += other.GetTotalLETT();
0060   fDTotalLETT += other.GetDTotalLETT();
0061   fDTotalLETD += other.GetDTotalLETD();
0062   fTotalLETD += other.GetTotalLETD();
0063
0064   // Merges ion counters
0065   auto rhsIonLetStore = other.GetIonLetStore();
0066
0067   // Loop over rhs ions
0068   for (unsigned int l = 0; l < rhsIonLetStore.size(); l++) {
0069     G4int PDGencoding = rhsIonLetStore[l].GetPDGencoding();
0070     size_t q;
0071     // Loop over lhs ions to find the one
0072     for (q = 0; q < fIonLetStore.size(); q++) {
0073       // Check if primary status is the same
0074       if (fIonLetStore[q].GetPDGencoding() == PDGencoding) {
0075         if (rhsIonLetStore[l].IsPrimary() == fIonLetStore[q].IsPrimary()) break;
0076       }
0077     }
0078
0079     if (q == fIonLetStore.size())  // Just copy the IonLet element in lhs store
0080       fIonLetStore.push_back(rhsIonLetStore[l]);
0081     else  // Merges rhs data with lhs ones
0082       fIonLetStore[q].Merge(&(rhsIonLetStore[l]));
0083   }
0084 }
0085
0086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0087
0088 void LETAccumulable::Reset()
0089 {
0090   if (GetVerboseLevel() > 0) {
0091     G4cout << "LETAccumulable::Reset()" << G4endl;
0092   }
0093   if (fInitialized) {
0094     fTotalLETT = 0.0;
0095     fDTotalLETT = 0.0;
0096     fDTotalLETD = 0.0;
0097     fTotalLETD = 0.0;
0098     fIonLetStore.clear();
0099   }
0100   else {
0101     Initialize();
0102   }
0103 }
0104
0105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0106
0107 // To accumulate given the hit
0108 void LETAccumulable::Accumulate(Hit* hit)
0109 {
0110   if (GetVerboseLevel() > 1) {
0111     G4cout << "LETAccumulable::Accumulate()" << G4endl;
0112   }
0113
0114   // No need to accumulate if no energy deposit or track length is zero.
0115   if (hit->GetDeltaE() == 0. || hit->GetTrackLength() == 0.) return;
0116
0117   // Atomic number
0118   G4int Z = hit->GetPartType()->GetAtomicNumber();
0119   if (Z < 1) return;  // Calculate only protons and ions
0120
0121   G4int PDGencoding = hit->GetPartType()->GetPDGEncoding();
0122
0123   if (PDGencoding == 22 || PDGencoding == 11) return;  // do not accumulate for electrons or gammas
0124
0125   // Get simple Particle data Group ID without excitation level
0126   PDGencoding -= PDGencoding % 10;
0127
0128   // Get voxel number
0129   G4int xIndex = hit->GetXindex();
0130   G4int yIndex = hit->GetYindex();
0131   ;
0132   G4int zIndex = hit->GetZindex();
0133   ;
0134
0135   G4int voxel =
0136     VoxelizedSensitiveDetector::GetInstance()->GetThisVoxelNumber(xIndex, yIndex, zIndex);
0137
0138   // Get mean kinetic energy
0139   G4double ekinMean = hit->GetEkinMean();
0140
0141   // Get particle definition
0142   const G4ParticleDefinition* particleDef = hit->GetPartType();
0143
0144   // Get material
0145   G4Material* mat = hit->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial();
0146
0147   // ICRU stopping power calculation
0148   G4EmCalculator emCal;
0149   // Use the mean kinetic energy of ions in a step to calculate ICRU stopping power
0150   G4double Lsn = emCal.ComputeElectronicDEDX(ekinMean, particleDef, mat);
0151   // G4cout << "ERASE ME. Lsn = " << Lsn << G4endl;
0152   // G4cout << "ERASE ME. ekinMean = " << ekinMean << G4endl;
0153
0154   // Get deposited energy
0155   G4double DE = hit->GetDeltaE();
0156
0157   // Get deposited energy considering secondary electrons
0158   G4double DEElectrons = hit->GetElectronEnergy();
0159
0160   // Get track length
0161   G4double DX = hit->GetTrackLength();
0162
0163   // Get track ID
0164   G4int trackID = hit->GetTrackID();
0165
0166   // Total LET calculation...
0167   // Total dose LET Numerator, including secondary electrons energy deposit
0168   fTotalLETD[voxel] += (DE + DEElectrons) * Lsn;
0169   // Total dose LET Denominator, including secondary electrons energy deposit
0170   fDTotalLETD[voxel] += DE + DEElectrons;
0171   // Total track LET Numerator
0172   fTotalLETT[voxel] += DX * Lsn;
0173   // Total track LET Denominator
0174   fDTotalLETT[voxel] += DX;
0175
0176   size_t l;
0177   for (l = 0; l < fIonLetStore.size(); l++) {
0178     // Judge species of the current particle and store it
0179     if (fIonLetStore[l].GetPDGencoding() == PDGencoding)
0180       if (((trackID == 1) && (fIonLetStore[l].IsPrimary()))
0181           || ((trackID != 1) && (!fIonLetStore[l].IsPrimary())))
0182         break;
0183   }
0184
0185   // Just another type of ion/particle for our store...
0186   if (l == fIonLetStore.size()) {
0187     // Mass number
0188     G4int A = particleDef->GetAtomicMass();
0189
0190     // Particle name
0191     G4String fullName = particleDef->GetParticleName();
0192     G4String name = fullName.substr(0, fullName.find("["));  // Cut excitation energy [x.y]
0193
0194     IonLet ion(trackID, PDGencoding, fullName, name, Z, A,
0195                VoxelizedSensitiveDetector::GetInstance()->GetTotalVoxelNumber());
0196     fIonLetStore.push_back(ion);
0197   }
0198
0199   // Calculate ions LET and store them
0200   fIonLetStore[l].Update(voxel, DE, DEElectrons, Lsn, DX);
0201 }
0202
0203 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0204
0205 G4int LETAccumulable::GetVerboseLevel() const
0206 {
0207   // return same level of LET class
0208   return Manager::GetInstance()->GetQuantity("LET")->GetVerboseLevel();
0209 }
0210
0211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0212
0213 void LETAccumulable::Initialize()
0214 {
0215   if (GetVerboseLevel() > 0) {
0216     G4cout << "LETAccumulable::Initialize(): " << G4endl;
0217   }
0218
0219   G4int voxNumber = VoxelizedSensitiveDetector::GetInstance()->GetTotalVoxelNumber();
0220
0221   fTotalLETT = array_type(0.0, voxNumber);
0222   fTotalLETD = array_type(0.0, voxNumber);
0223   fDTotalLETT = array_type(0.0, voxNumber);
0224   fDTotalLETD = array_type(0.0, voxNumber);
0225   fInitialized = true;
0226 }
0227
0228 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0229
0230 }  // namespace RadioBio