hadrontherapy/src/HadrontherapyLet.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 // Hadrontherapy advanced example for Geant4
0027 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy
0028
0029 #include "HadrontherapyDetectorConstruction.hh"
0030 #include "HadrontherapyLet.hh"
0031
0032 #include "HadrontherapyMatrix.hh"
0033 #include "HadrontherapyInteractionParameters.hh"
0034 #include "HadrontherapyPrimaryGeneratorAction.hh"
0035 #include "HadrontherapyMatrix.hh"
0036 #include "G4AnalysisManager.hh"
0037 #include "G4RunManager.hh"
0038 #include "G4SystemOfUnits.hh"
0039 #include <cmath>
0040
0041 HadrontherapyLet* HadrontherapyLet::instance = NULL;
0042 G4bool HadrontherapyLet::doCalculation = false;
0043
0044 HadrontherapyLet* HadrontherapyLet::GetInstance(HadrontherapyDetectorConstruction *pDet)
0045 {
0046     if (instance) delete instance;
0047     instance = new HadrontherapyLet(pDet);
0048     return instance;
0049 }
0050
0051 HadrontherapyLet* HadrontherapyLet::GetInstance()
0052 {
0053     return instance;
0054 }
0055
0056 HadrontherapyLet::HadrontherapyLet(HadrontherapyDetectorConstruction* pDet)
0057 :filename("Let.out"),matrix(0) // Default output filename
0058 {
0059
0060     matrix = HadrontherapyMatrix::GetInstance();
0061
0062     if (!matrix)
0063         G4Exception("HadrontherapyLet::HadrontherapyLet",
0064                     "Hadrontherapy0005", FatalException,
0065                     "HadrontherapyMatrix not found. Firstly create an instance of it.");
0066
0067     nVoxels = matrix -> GetNvoxel();
0068
0069     numberOfVoxelAlongX = matrix -> GetNumberOfVoxelAlongX();
0070     numberOfVoxelAlongY = matrix -> GetNumberOfVoxelAlongY();
0071     numberOfVoxelAlongZ = matrix -> GetNumberOfVoxelAlongZ();
0072
0073     G4RunManager *runManager = G4RunManager::GetRunManager();
0074     pPGA = (HadrontherapyPrimaryGeneratorAction*)runManager -> GetUserPrimaryGeneratorAction();
0075     // Pointer to the detector material
0076     detectorMat = pDet -> GetDetectorLogicalVolume() -> GetMaterial();
0077     density = detectorMat -> GetDensity();
0078     // Instances for Total LET
0079     totalLetD =      new G4double[nVoxels];
0080     DtotalLetD =     new G4double[nVoxels];
0081     totalLetT =      new G4double[nVoxels];
0082     DtotalLetT =     new G4double[nVoxels];
0083
0084 }
0085
0086 HadrontherapyLet::~HadrontherapyLet()
0087 {
0088     Clear();
0089     delete [] totalLetD;
0090     delete [] DtotalLetD;
0091     delete [] totalLetT;
0092     delete [] DtotalLetT;
0093 }
0094
0095 // Fill energy spectrum for every voxel (local energy spectrum)
0096 void HadrontherapyLet::Initialize()
0097 {
0098     for(G4int v=0; v < nVoxels; v++) totalLetD[v] = DtotalLetD[v] = totalLetT[v] = DtotalLetT[v] = 0.;
0099     Clear();
0100 }
0101 /**
0102  * Clear all stored data
0103  */
0104 void HadrontherapyLet::Clear()
0105 {
0106     for (size_t i=0; i < ionLetStore.size(); i++)
0107     {
0108         delete [] ionLetStore[i].letDN; // Let Dose Numerator
0109         delete [] ionLetStore[i].letDD; // Let Dose Denominator
0110         delete [] ionLetStore[i].letTN; // Let Track Numerator
0111         delete [] ionLetStore[i].letTD; // Let Track Denominator
0112     }
0113     ionLetStore.clear();
0114 }
0115 void  HadrontherapyLet::FillEnergySpectrum(G4int trackID,
0116                                            G4ParticleDefinition* particleDef,
0117                                            G4double ekinMean,
0118                                            const G4Material* mat,
0119                                            G4double DE,
0120                                            G4double DEEletrons,
0121                                            G4double DX,
0122                                            G4int i, G4int j, G4int k)
0123 {
0124     if (DE <= 0. || DX <=0. ) return; // calculate only energy deposit
0125     if (!doCalculation) return;
0126
0127     // atomic number
0128     G4int Z = particleDef -> GetAtomicNumber();
0129     if (Z<1) return; // calculate only protons and ions
0130
0131     G4int PDGencoding = particleDef -> GetPDGEncoding();
0132     PDGencoding -= PDGencoding%10; // simple Particle data group id  without excitation level
0133
0134     G4int voxel = matrix -> Index(i,j,k);
0135
0136     // ICRU stopping power calculation
0137     G4EmCalculator emCal;
0138     // use the mean kinetic energy of ions in a step to calculate ICRU stopping power
0139     G4double Lsn = emCal.ComputeElectronicDEDX(ekinMean, particleDef, mat);
0140
0141
0142     // Total LET calculation...
0143     totalLetD[voxel]  += (DE + DEEletrons) * Lsn; // total dose Let Numerator, including secondary electrons energy deposit
0144     DtotalLetD[voxel] += DE + DEEletrons;         // total dose Let Denominator, including secondary electrons energy deposit
0145     totalLetT[voxel]  += DX * Lsn;                // total track Let Numerator
0146     DtotalLetT[voxel] += DX;                      // total track Let Denominator
0147
0148     // store primary ions and secondary ions
0149     size_t l;
0150     for (l=0; l < ionLetStore.size(); l++)
0151     {
0152         // judge species of the current particle and store it
0153         if (ionLetStore[l].PDGencoding == PDGencoding)
0154             if ( ((trackID ==1) && (ionLetStore[l].isPrimary)) || ((trackID !=1) && (!ionLetStore[l].isPrimary)))
0155                 break;
0156     }
0157
0158     if (l == ionLetStore.size()) // Just another type of ion/particle for our store...
0159     {
0160         // mass number
0161         G4int A = particleDef -> GetAtomicMass();
0162
0163         // particle name
0164         G4String fullName = particleDef -> GetParticleName();
0165         G4String name = fullName.substr (0, fullName.find("[") ); // cut excitation energy [x.y]
0166
0167         ionLet ion =
0168         {
0169             (trackID == 1) ? true:false, // is it the primary particle?
0170             PDGencoding,
0171             fullName,
0172             name,
0173             Z,
0174             A,
0175             new G4double[nVoxels], // Let Dose Numerator
0176             new G4double[nVoxels],  // Let Dose Denominator
0177             new G4double[nVoxels], // Let Track Numerator
0178             new G4double[nVoxels],  // Let Track Denominator
0179         };
0180
0181         // Initialize let and other parameters
0182         for(G4int v=0; v < nVoxels; v++)
0183         {
0184             ion.letDN[v] = ion.letDD[v] = ion.letTN[v] = ion.letTD[v] = 0.;
0185         }
0186
0187
0188         ionLetStore.push_back(ion);
0189     }
0190
0191     // calculate ions let and store them
0192     ionLetStore[l].letDN[voxel] += (DE + DEEletrons)* Lsn; // ions dose Let Numerator, including secondary electrons energy deposit
0193     ionLetStore[l].letDD[voxel] += DE + DEEletrons;        // ions dose Let Denominator, including secondary electrons energy deposit
0194     ionLetStore[l].letTN[voxel] += DX* Lsn;                // ions track Let Numerator
0195     ionLetStore[l].letTD[voxel] += DX;                     // ions track Let Denominator
0196
0197 }
0198
0199
0200
0201
0202 void HadrontherapyLet::LetOutput()
0203 {
0204     for(G4int v=0; v < nVoxels; v++)
0205     {
0206         // compute total let
0207         if (DtotalLetD[v]>0.) totalLetD[v] = totalLetD[v]/DtotalLetD[v];
0208         if (DtotalLetT[v]>0.) totalLetT[v] = totalLetT[v]/DtotalLetT[v];
0209     }
0210
0211     // Sort ions by A and then by Z ...
0212     std::sort(ionLetStore.begin(), ionLetStore.end());
0213
0214
0215     // Compute Let Track and Let Dose for ions
0216
0217     for(G4int v=0; v < nVoxels; v++)
0218     {
0219
0220         for (size_t ion=0; ion < ionLetStore.size(); ion++)
0221         {
0222             // compute ions let
0223             if (ionLetStore[ion].letDD[v] >0.) ionLetStore[ion].letDN[v] = ionLetStore[ion].letDN[v] / ionLetStore[ion].letDD[v];
0224             if (ionLetStore[ion].letTD[v] >0.) ionLetStore[ion].letTN[v] = ionLetStore[ion].letTN[v] / ionLetStore[ion].letTD[v];
0225         } // end loop over ions
0226     }
0227 } // end loop over voxels
0228
0229
0230
0231 void HadrontherapyLet::StoreLetAscii()
0232 {
0233 #define width 15L
0234
0235     if(ionLetStore.size())
0236     {
0237         ofs.open(filename, std::ios::out);
0238         if (ofs.is_open())
0239         {
0240
0241             // Write the voxels index and total Lets and the list of particles/ions
0242             ofs << "i" << '\t' << "j" << '\t' << "k";
0243
0244             ofs <<  '\t' << "LDT";
0245             ofs <<  '\t' << "LTT";
0246
0247             for (size_t l=0; l < ionLetStore.size(); l++) // write ions name
0248             {
0249                 G4String a = (ionLetStore[l].isPrimary) ? "_1_D":"_D";
0250                 ofs << '\t' << ionLetStore[l].name  + a ;
0251                 G4String b = (ionLetStore[l].isPrimary) ? "_1_T":"_T";
0252                 ofs << '\t' << ionLetStore[l].name  + b ;
0253             }
0254
0255
0256             // Write data
0257
0258             G4AnalysisManager*  LetFragmentTuple = G4AnalysisManager::Instance();
0259
0260             LetFragmentTuple->SetVerboseLevel(1);
0261             LetFragmentTuple->SetFirstHistoId(1);
0262             LetFragmentTuple->SetFirstNtupleId(1);
0263             LetFragmentTuple ->OpenFile("Let.csv");
0264
0265
0266             LetFragmentTuple ->CreateNtuple("coordinate", "Let");
0267
0268
0269             LetFragmentTuple ->CreateNtupleIColumn("i");//0
0270             LetFragmentTuple ->CreateNtupleIColumn("j");//1
0271             LetFragmentTuple ->CreateNtupleIColumn("k");//2
0272             LetFragmentTuple ->CreateNtupleDColumn("TotalLetD");//3
0273             LetFragmentTuple ->CreateNtupleDColumn("TotalLetT");//4
0274             LetFragmentTuple ->CreateNtupleIColumn("A");//5
0275             LetFragmentTuple ->CreateNtupleIColumn("Z");//6
0276             LetFragmentTuple ->CreateNtupleDColumn("IonLETD");//7
0277             LetFragmentTuple ->CreateNtupleDColumn("IonLETT");//8
0278             LetFragmentTuple ->FinishNtuple();
0279
0280
0281             for(G4int i = 0; i < numberOfVoxelAlongX; i++)
0282                 for(G4int j = 0; j < numberOfVoxelAlongY; j++)
0283                     for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
0284                     {
0285                         LetFragmentTuple->FillNtupleIColumn(1,0, i);
0286                         LetFragmentTuple->FillNtupleIColumn(1,1, j);
0287                         LetFragmentTuple->FillNtupleIColumn(1,2, k);
0288
0289                         G4int v = matrix -> Index(i, j, k);
0290
0291                         // write total Lets and voxels index
0292                         ofs << G4endl;
0293                         ofs << i << '\t' << j << '\t' << k ;
0294                         ofs << '\t' << totalLetD[v]/(keV/um);
0295                         ofs << '\t' << totalLetT[v]/(keV/um);
0296
0297
0298                         // write ions Lets
0299                         for (size_t l=0; l < ionLetStore.size(); l++)
0300                         {
0301
0302                             // Write only not identically null data line
0303                             if(ionLetStore[l].letDN)
0304                             {
0305                                 // write ions Lets
0306                                 ofs << '\t' << ionLetStore[l].letDN[v]/(keV/um) ;
0307                                 ofs << '\t' << ionLetStore[l].letTN[v]/(keV/um);
0308                             }
0309                         }
0310
0311                         LetFragmentTuple->FillNtupleDColumn(1,3, totalLetD[v]/(keV/um));
0312                         LetFragmentTuple->FillNtupleDColumn(1,4, totalLetT[v]/(keV/um));
0313
0314
0315                         for (size_t ll=0; ll < ionLetStore.size(); ll++)
0316                         {
0317
0318
0319                             LetFragmentTuple->FillNtupleIColumn(1,5, ionLetStore[ll].A);
0320                             LetFragmentTuple->FillNtupleIColumn(1,6, ionLetStore[ll].Z);
0321
0322
0323                             LetFragmentTuple->FillNtupleDColumn(1,7, ionLetStore[ll].letDN[v]/(keV/um));
0324                             LetFragmentTuple->FillNtupleDColumn(1,8, ionLetStore[ll].letTN[v]/(keV/um));
0325                             LetFragmentTuple->AddNtupleRow(1);
0326                         }
0327                     }
0328             ofs.close();
0329
0330             LetFragmentTuple->Write();
0331             LetFragmentTuple->CloseFile();
0332         }
0333
0334     }
0335
0336 }
0337