radiobiology/src/LET.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 //
0027 /// \file radiobiology/src/LET.cc
0028 /// \brief Implementation of the RadioBio::LET class
0029
0030 #include "LET.hh"
0031
0032 #include "G4EmCalculator.hh"
0033 #include "G4RunManager.hh"
0034 #include "G4SystemOfUnits.hh"
0035
0036 #include "DetectorConstruction.hh"
0037 #include "Hit.hh"
0038 #include "LETAccumulable.hh"
0039 #include "LETMessenger.hh"
0040 #include "VoxelizedSensitiveDetector.hh"
0041
0042 #include <cmath>
0043
0044 namespace RadioBio
0045 {
0046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0047
0048 LET::LET() : VRadiobiologicalQuantity()
0049 {
0050   fPath = "RadioBio_LET.out";  // Default output filename
0051
0052   fMessenger = new LETMessenger(this);
0053
0054   Initialize();
0055 }
0056
0057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0058
0059 LET::~LET()
0060 {
0061   delete fMessenger;
0062 }
0063
0064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0065
0066 void LET::Initialize()
0067 {
0068   G4int VoxelNumber = VoxelizedSensitiveDetector::GetInstance()->GetTotalVoxelNumber();
0069
0070   // Instances for Total LET
0071   fNTotalLETT.resize(VoxelNumber);
0072   fNTotalLETD.resize(VoxelNumber);
0073   fDTotalLETT.resize(VoxelNumber);
0074   fDTotalLETD.resize(VoxelNumber);
0075
0076   fTotalLETD.resize(VoxelNumber);
0077   fTotalLETT.resize(VoxelNumber);
0078
0079   fInitialized = true;
0080 }
0081
0082 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0083
0084 void LET::Compute()
0085 {
0086   // Skip LET computation if calculation not enabled.
0087   if (!fCalculationEnabled) {
0088     if (fVerboseLevel > 0) {
0089       G4cout << "LET::Compute() called but skipped as calculation not enabled" << G4endl;
0090     }
0091     return;
0092   }
0093   if (fVerboseLevel > 0) G4cout << "LET::Compute()" << G4endl;
0094
0095   if (VoxelizedSensitiveDetector::GetInstance() == nullptr)
0096     G4Exception("LET::ComputeLET", "PointerNotAvailable", FatalException,
0097                 "Computing LET without voxelized geometry pointer!");
0098
0099   G4int VoxelNumber = VoxelizedSensitiveDetector::GetInstance()->GetTotalVoxelNumber();
0100
0101   for (G4int v = 0; v < VoxelNumber; v++) {
0102     if (fVerboseLevel > 1) G4cout << "COMPUTING LET of voxel " << v << G4endl;
0103
0104     // Compute total LET
0105     if (fDTotalLETD[v] > 0.) fTotalLETD[v] = fNTotalLETD[v] / fDTotalLETD[v];
0106     // G4cout << "ERASE ME. DEBUG: voxel = " << v << " fNTotalLETD[v] = " << fNTotalLETD[v]
0107     //        << " fDTotalLETD[v] = " << fDTotalLETD[v] << G4endl;
0108     if (fDTotalLETT[v] > 0.) fTotalLETT[v] = fNTotalLETT[v] / fDTotalLETT[v];
0109   }
0110
0111   // Sort ions by A and then by Z ...
0112   std::sort(fIonLetStore.begin(), fIonLetStore.end());
0113
0114   // Compute LET Track and LET Dose for ions
0115   for (size_t ion = 0; ion < fIonLetStore.size(); ion++) {
0116     fIonLetStore[ion].Calculate();
0117   }
0118
0119   fCalculated = true;
0120 }
0121
0122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0123
0124 // save LET
0125 void LET::Store()
0126 {
0127   // Skip LET storing if calculation not enabled.
0128   if (!fCalculationEnabled) {
0129     if (fVerboseLevel > 0) {
0130       G4cout << "LET::Store() called but skipped as calculation not enabled" << G4endl;
0131     }
0132     return;
0133   }
0134
0135   if (fSaved == true)
0136     G4Exception("LET::StoreLET", "NtuplesAlreadySaved", FatalException,
0137                 "Overwriting LET file. For multiple runs, change filename.");
0138
0139   Compute();
0140 #define width 15L
0141   if (fVerboseLevel > 0) G4cout << "LET::StoreLET" << G4endl;
0142
0143   // If there is at least one ion
0144   if (fIonLetStore.size()) {
0145     std::ofstream ofs(fPath);
0146     // ofs.open(fPath, std::ios::out);
0147     if (ofs.is_open()) {
0148       // Write the voxels index and total LETs and the list of particles/ions
0149       ofs << std::setprecision(6) << std::left << "i\tj\tk\t";
0150       ofs << std::setw(width) << "LDT";
0151       ofs << std::setw(width) << "LTT";
0152
0153       for (size_t l = 0; l < fIonLetStore.size(); l++)  // Write ions name
0154       {
0155         G4String a = (fIonLetStore[l].IsPrimary()) ? "_1_D" : "_D";
0156         ofs << std::setw(width) << fIonLetStore[l].GetName() + a;
0157         G4String b = (fIonLetStore[l].IsPrimary()) ? "_1_T" : "_T";
0158         ofs << std::setw(width) << fIonLetStore[l].GetName() + b;
0159       }
0160       ofs << G4endl;
0161
0162       // Write data
0163       G4AnalysisManager* LETFragmentTuple = G4AnalysisManager::Instance();
0164
0165       LETFragmentTuple->SetVerboseLevel(1);
0166       LETFragmentTuple->SetFirstHistoId(1);
0167       LETFragmentTuple->SetFirstNtupleId(1);
0168
0169       LETFragmentTuple->SetDefaultFileType("xml");
0170       LETFragmentTuple->OpenFile("LET");
0171
0172       LETFragmentTuple->CreateNtuple("coordinate", "LET");
0173
0174       LETFragmentTuple->CreateNtupleIColumn("i");  // 0
0175       LETFragmentTuple->CreateNtupleIColumn("j");  // 1
0176       LETFragmentTuple->CreateNtupleIColumn("k");  // 2
0177       LETFragmentTuple->CreateNtupleDColumn("TotalLETD");  // 3
0178       LETFragmentTuple->CreateNtupleDColumn("TotalLETT");  // 4
0179       LETFragmentTuple->CreateNtupleIColumn("A");  // 5
0180       LETFragmentTuple->CreateNtupleIColumn("Z");  // 6
0181       LETFragmentTuple->CreateNtupleDColumn("IonLetD");  // 7
0182       LETFragmentTuple->CreateNtupleDColumn("IonLetT");  // 8
0183       LETFragmentTuple->FinishNtuple();
0184
0185       auto voxSensDet = VoxelizedSensitiveDetector::GetInstance();
0186
0187       for (G4int i = 0; i < voxSensDet->GetVoxelNumberAlongX(); i++)
0188         for (G4int j = 0; j < voxSensDet->GetVoxelNumberAlongY(); j++)
0189           for (G4int k = 0; k < voxSensDet->GetVoxelNumberAlongZ(); k++) {
0190             LETFragmentTuple->FillNtupleIColumn(1, 0, i);
0191             LETFragmentTuple->FillNtupleIColumn(1, 1, j);
0192             LETFragmentTuple->FillNtupleIColumn(1, 2, k);
0193
0194             G4int v = voxSensDet->GetThisVoxelNumber(i, j, k);
0195
0196             // Write total LETs and voxels index
0197             ofs << G4endl;
0198             ofs << i << '\t' << j << '\t' << k << '\t';
0199             ofs << std::setw(width) << fTotalLETD[v] / (keV / um);
0200             ofs << std::setw(width) << fTotalLETT[v] / (keV / um);
0201
0202             // Write ions LETs
0203             for (size_t l = 0; l < fIonLetStore.size(); l++) {
0204               // Write ions LETs
0205               ofs << std::setw(width) << fIonLetStore[l].GetLETD()[v] / (keV / um);
0206               ofs << std::setw(width) << fIonLetStore[l].GetLETT()[v] / (keV / um);
0207             }
0208
0209             LETFragmentTuple->FillNtupleDColumn(1, 3, fTotalLETD[v] / (keV / um));
0210             LETFragmentTuple->FillNtupleDColumn(1, 4, fTotalLETT[v] / (keV / um));
0211
0212             for (size_t ll = 0; ll < fIonLetStore.size(); ll++) {
0213               LETFragmentTuple->FillNtupleIColumn(1, 5, fIonLetStore[ll].GetA());
0214               LETFragmentTuple->FillNtupleIColumn(1, 6, fIonLetStore[ll].GetZ());
0215
0216               LETFragmentTuple->FillNtupleDColumn(1, 7, fIonLetStore[ll].GetLETD()[v] / (keV / um));
0217               LETFragmentTuple->FillNtupleDColumn(1, 8, fIonLetStore[ll].GetLETT()[v] / (keV / um));
0218               LETFragmentTuple->AddNtupleRow(1);
0219             }
0220           }
0221       ofs.close();
0222
0223       // LETFragmentTuple->Write();
0224       LETFragmentTuple->CloseFile();
0225     }
0226   }
0227
0228   fSaved = true;
0229 }
0230
0231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0232
0233 void LET::Reset()
0234 {
0235   if (fVerboseLevel > 1) {
0236     G4cout << "LET::Reset(): ";
0237   }
0238   fNTotalLETT = 0.0;
0239   fNTotalLETD = 0.0;
0240   fDTotalLETT = 0.0;
0241   fDTotalLETD = 0.0;
0242
0243   fTotalLETD = 0.0;
0244   fTotalLETT = 0.0;
0245
0246   fIonLetStore.clear();
0247   fCalculated = false;
0248 }
0249
0250 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0251
0252 // Add data taken from accumulables
0253 void LET::AddFromAccumulable(G4VAccumulable* GenAcc)
0254 {
0255   LETAccumulable* acc = (LETAccumulable*)GenAcc;
0256
0257   AddNTotalLETT(acc->GetTotalLETT());
0258   AddDTotalLETT(acc->GetDTotalLETT());
0259   AddNTotalLETD(acc->GetTotalLETD());
0260   AddDTotalLETD(acc->GetDTotalLETD());
0261
0262   // Merges ion counters
0263   // Loop over rhs ions
0264   for (unsigned int l = 0; l < acc->GetIonLetStore().size(); l++) {
0265     G4int PDGencoding = acc->GetIonLetStore()[l].GetPDGencoding();
0266     size_t q;
0267     // Loop over lhs ions to find the one
0268     for (q = 0; q < fIonLetStore.size(); q++) {
0269       // If he found it, sums values
0270       if (fIonLetStore[q].GetPDGencoding() == PDGencoding) {
0271         if (acc->GetIonLetStore()[l].IsPrimary() == fIonLetStore[q].IsPrimary()) break;
0272       }
0273     }
0274     // If ion is missing, copy it
0275     if (q == fIonLetStore.size())
0276       fIonLetStore.push_back(acc->GetIonLetStore()[l]);
0277     else  // Merge rhs data with lhs ones
0278       fIonLetStore[q].Merge(&(acc->GetIonLetStore()[l]));
0279   }
0280   fCalculated = false;
0281 }
0282
0283 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0284
0285 void LET::PrintParameters()
0286 {
0287   G4cout << "*******************************************" << G4endl
0288          << "******* Parameters of the class LET *******" << G4endl
0289          << "*******************************************" << G4endl;
0290   PrintVirtualParameters();
0291 }
0292
0293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0294
0295 }  // namespace RadioBio