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