Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-02-26 09:19:38

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 <fstream>
0030 #include <iostream>
0031 #include <sstream>
0032 #include <iomanip>
0033 
0034 #include "HadrontherapyMatrix.hh"
0035 #include "HadrontherapyPrimaryGeneratorAction.hh"
0036 #include "globals.hh"
0037 #include "G4SystemOfUnits.hh"
0038 #include "G4RunManager.hh"
0039 #include "G4ParticleGun.hh"
0040 #include "HadrontherapySteppingAction.hh"
0041 #include "HadrontherapyAnalysisFileMessenger.hh"
0042 #include "G4SystemOfUnits.hh"
0043 #include <time.h>
0044 
0045 HadrontherapyAnalysis* HadrontherapyAnalysis::instance = 0;
0046 /////////////////////////////////////////////////////////////////////////////
0047 
0048 HadrontherapyAnalysis::HadrontherapyAnalysis()
0049 {
0050     fMess = new HadrontherapyAnalysisFileMessenger(this);
0051 }
0052 
0053 /////////////////////////////////////////////////////////////////////////////
0054 HadrontherapyAnalysis::~HadrontherapyAnalysis()
0055 {
0056     delete fMess;
0057 }
0058 
0059 /////////////////////////////////////////////////////////////////////////////
0060 HadrontherapyAnalysis* HadrontherapyAnalysis::GetInstance(){
0061 
0062     if (instance == 0) instance = new HadrontherapyAnalysis;
0063     return instance;
0064 }
0065 
0066 HadrontherapyMatrix* HadrontherapyMatrix::instance = NULL;
0067 G4bool HadrontherapyMatrix::secondary = false;
0068 
0069 
0070 // Only return a pointer to matrix
0071 HadrontherapyMatrix* HadrontherapyMatrix::GetInstance()
0072 {
0073     return instance;
0074 }
0075 
0076 /////////////////////////////////////////////////////////////////////////////
0077 // This STATIC method delete (!) the old matrix and rewrite a new object returning a pointer to it
0078 // TODO A check on the parameters is required!
0079 HadrontherapyMatrix* HadrontherapyMatrix::GetInstance(G4int voxelX, G4int voxelY, G4int voxelZ, G4double mass)
0080 {
0081     if (instance) delete instance;
0082     instance = new HadrontherapyMatrix(voxelX, voxelY, voxelZ, mass);
0083     instance -> Initialize();
0084     return instance;
0085 }
0086 
0087 /////////////////////////////////////////////////////////////////////////////
0088 HadrontherapyMatrix::HadrontherapyMatrix(G4int voxelX, G4int voxelY, G4int voxelZ, G4double mass):
0089     stdFile("Dose.out"),
0090     doseUnit(gray)
0091 {
0092     // Number of the voxels of the phantom
0093     // For Y = Z = 1 the phantom is divided in slices (and not in voxels)
0094     // orthogonal to the beam axis
0095     numberOfVoxelAlongX = voxelX;
0096     numberOfVoxelAlongY = voxelY;
0097     numberOfVoxelAlongZ = voxelZ;
0098     massOfVoxel = mass;
0099 
0100 
0101     // Create the dose matrix
0102     matrix = new G4double[numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ];
0103     if (matrix)
0104     {
0105         G4cout << "HadrontherapyMatrix: Memory space to store physical dose into " <<
0106                   numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ <<
0107                   " voxels has been allocated " << G4endl;
0108     }
0109 
0110     else G4Exception("HadrontherapyMatrix::HadrontherapyMatrix()", "Hadrontherapy0005", FatalException, "Can't allocate memory to store physical dose!");
0111 
0112 
0113     // Hit voxel (TrackID) marker
0114     // This array mark the status of voxel, if a hit occur, with the trackID of the particle
0115     // Must be initialized
0116 
0117     hitTrack = new G4int[numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ];
0118     ClearHitTrack();
0119 }
0120 
0121 /////////////////////////////////////////////////////////////////////////////
0122 HadrontherapyMatrix::~HadrontherapyMatrix()
0123 {
0124     delete[] matrix;
0125     delete[] hitTrack;
0126     Clear();
0127 }
0128 
0129 /////////////////////////////////////////////////////////////////////////////
0130 void HadrontherapyMatrix::Clear()
0131 {
0132     for (size_t i=0; i<ionStore.size(); i++)
0133     {
0134         delete[] ionStore[i].dose;
0135         delete[] ionStore[i].fluence;
0136     }
0137     ionStore.clear();
0138 }
0139 
0140 /////////////////////////////////////////////////////////////////////////////
0141 // Initialise the elements of the matrix to zero
0142 
0143 void HadrontherapyMatrix::Initialize()
0144 {
0145     // Clear ions store
0146     Clear();
0147     // Clear dose
0148     for(int i=0;i<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ;i++)
0149     {
0150         matrix[i] = 0;
0151     }
0152 }
0153 
0154 /////////////////////////////////////////////////////////////////////////////
0155 // Print generated nuclides list
0156 
0157 /////////////////////////////////////////////////////////////////////////////
0158 void HadrontherapyMatrix::PrintNuclides()
0159 {
0160     for (size_t i=0; i<ionStore.size(); i++)
0161     {
0162         G4cout << ionStore[i].name << G4endl;
0163     }
0164 }
0165 
0166 /////////////////////////////////////////////////////////////////////////////
0167 // Clear Hit voxel (TrackID) markers
0168 
0169 void HadrontherapyMatrix::ClearHitTrack()
0170 {
0171     for(G4int i=0; i<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ; i++) hitTrack[i] = 0;
0172 }
0173 
0174 // Return Hit status
0175 G4int* HadrontherapyMatrix::GetHitTrack(G4int i, G4int j, G4int k)
0176 {
0177     return &(hitTrack[Index(i,j,k)]);
0178 }
0179 
0180 /////////////////////////////////////////////////////////////////////////////
0181 // Dose methods...
0182 // Fill DOSE/fluence matrix for secondary particles:
0183 // If fluence parameter is true (default value is FALSE) then fluence at voxel (i, j, k) is increased.
0184 // The energyDeposit parameter fill the dose matrix for voxel (i,j,k)
0185 /////////////////////////////////////////////////////////////////////////////
0186 G4bool HadrontherapyMatrix::Fill(G4int trackID,
0187                                  G4ParticleDefinition* particleDef,
0188                                  G4int i, G4int j, G4int k,
0189                                  G4double energyDeposit,
0190                                  G4bool fluence)
0191 {
0192 
0193     if ( (energyDeposit <=0. && !fluence) || !secondary) return false;
0194 
0195     // Get Particle Data Group particle ID
0196     G4int PDGencoding = particleDef -> GetPDGEncoding();
0197     PDGencoding -= PDGencoding%10;
0198 
0199     // Search for already allocated data...
0200     for (size_t l=0; l < ionStore.size(); l++)
0201     {
0202         if (ionStore[l].PDGencoding == PDGencoding )
0203         {   // Is it a primary or a secondary particle?
0204 
0205             if ( (trackID ==1 && ionStore[l].isPrimary) || (trackID !=1 && !ionStore[l].isPrimary))
0206             {
0207                 if (energyDeposit > 0.)
0208 
0209                     ionStore[l].dose[Index(i, j, k)] += energyDeposit;
0210 
0211                 // Fill a matrix per each ion with the fluence
0212 
0213                 if (fluence) ionStore[l].fluence[Index(i, j, k)]++;
0214                 return true;
0215             }
0216         }
0217     }
0218     G4int Z = particleDef-> GetAtomicNumber();
0219     G4int A = particleDef-> GetAtomicMass();
0220     G4String fullName = particleDef -> GetParticleName();
0221     G4String name = fullName.substr (0, fullName.find("[") ); // cut excitation energy
0222 
0223     // Let's put a new particle in our store...
0224     ion newIon =
0225     {
0226         (trackID == 1) ? true:false,
0227         PDGencoding,
0228         name,
0229         name.length(),
0230         Z,
0231         A,
0232         new G4double[numberOfVoxelAlongX * numberOfVoxelAlongY * numberOfVoxelAlongZ],
0233         new unsigned int[numberOfVoxelAlongX * numberOfVoxelAlongY * numberOfVoxelAlongZ]
0234     };
0235 
0236 
0237     // Initialize data
0238     if (newIon.dose && newIon.fluence)
0239     {
0240         for(G4int q=0; q<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ; q++)
0241         {
0242             newIon.dose[q] = 0.;
0243             newIon.fluence[q] = 0;
0244         }
0245 
0246         if (energyDeposit > 0.) newIon.dose[Index(i, j, k)] += energyDeposit;
0247         if (fluence) newIon.fluence[Index(i, j, k)]++;
0248 
0249         ionStore.push_back(newIon);
0250         return true;
0251     }
0252 
0253     else // XXX Out of memory! XXX
0254 
0255     {
0256         return false;
0257     }
0258 }
0259 
0260 /////////////////////////////////////////////////////////////////////////////
0261 /////////////////////////////////////////////////////////////////////////////
0262 // Methods to store data to filenames...
0263 ////////////////////////////////////////////////////////////////////////////
0264 ////////////////////////////////////////////////////////////////////////////
0265 //
0266 // General method to store matrix data to filename
0267 void HadrontherapyMatrix::StoreMatrix(G4String file, void* data, size_t psize)
0268 {
0269     if (data)
0270     {
0271         ofs.open(file, std::ios::out);
0272         if (ofs.is_open())
0273         {
0274             for(G4int i = 0; i < numberOfVoxelAlongX; i++)
0275                 for(G4int j = 0; j < numberOfVoxelAlongY; j++)
0276                     for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
0277                     {
0278                         G4int n = Index(i, j, k);
0279 
0280                         if (psize == sizeof(unsigned int))
0281                         {
0282                             unsigned int* pdata = (unsigned int*)data;
0283 
0284                             if (pdata[n])
0285 
0286                                 ofs << i << '\t' << j << '\t' << k << '\t' << pdata[n] << G4endl;
0287                         }
0288 
0289                         else if (psize == sizeof(G4double))
0290 
0291                         {
0292                             G4double* pdata = (G4double*)data;
0293                             if (pdata[n]) ofs << i << '\t' << j << '\t' << k << '\t' << pdata[n] << G4endl;
0294                         }
0295                     }
0296             ofs.close();
0297         }
0298     }
0299 }
0300 
0301 /////////////////////////////////////////////////////////////////////////////
0302 // Store fluence per single ion in multiple files
0303 void HadrontherapyMatrix::StoreFluenceData()
0304 {
0305     for (size_t i=0; i < ionStore.size(); i++){
0306         StoreMatrix(ionStore[i].name + "_Fluence.out", ionStore[i].fluence, sizeof(unsigned int));
0307     }
0308 }
0309 
0310 /////////////////////////////////////////////////////////////////////////////
0311 // Store dose per single ion in multiple files
0312 void HadrontherapyMatrix::StoreDoseData()
0313 {
0314 
0315     for (size_t i=0; i < ionStore.size(); i++){
0316         StoreMatrix(ionStore[i].name + "_Dose.out", ionStore[i].dose, sizeof(G4double));
0317     }
0318 }
0319 
0320 ////////////////////////////////////////////////////////////////////////
0321 // Store dose into a single file
0322 // or in histograms. Please, note that this function is called via
0323 // messenger commands
0324 // defined in the HadrontherapyAnalysisFileMessenger.cc class file
0325 void HadrontherapyMatrix::StoreDoseFluenceAscii(G4String file)
0326 {
0327 #define width 15L
0328     filename = (file=="") ? stdFile:file;
0329 
0330     // Sort like periodic table
0331 
0332     std::sort(ionStore.begin(), ionStore.end());
0333     G4cout << "Dose is being written to " << filename << G4endl;
0334     ofs.open(filename, std::ios::out);
0335 
0336     if (ofs.is_open())
0337     {
0338         // Write the voxels index and the list of particles/ions
0339         //ofs << std::setprecision(6) << std::left << "i\tj\tk\t";
0340         ofs << "i" << '\t' << "j" << '\t' << "k";
0341         G4cout << "i" << '\t' << "j" << '\t' << "k";
0342 
0343         // Total dose
0344         ofs <<'\t' <<"Dose(Gy)";
0345         //ofs << std::setw(width) << "Dose(Gy)";
0346         G4cout << '\t' << "Dose(Gy)";
0347 
0348         G4String fluence = "_f";
0349         if (secondary)
0350         {
0351             for (size_t l=0; l < ionStore.size(); l++)
0352             {
0353                 G4String a = (ionStore[l].isPrimary) ? "_1":"";     // is it a primary?
0354 
0355                // ofs << std::setw(width) << ionStore[l].name + a <<
0356                //        std::setw(width) << ionStore[l].name  + a + fluence;
0357 
0358                 ofs << '\t' << ionStore[l].name + a <<
0359                           '\t' << ionStore[l].name  + a + fluence;
0360 
0361                 G4cout << '\t' << ionStore[l].name + a <<
0362                           '\t' << ionStore[l].name  + a + fluence;
0363 
0364             }
0365             //ofs << G4endl;
0366         }
0367 
0368         // Write data
0369         for(G4int i = 0; i < numberOfVoxelAlongX; i++)
0370             for(G4int j = 0; j < numberOfVoxelAlongY; j++)
0371                 for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
0372                 {
0373                     G4int n = Index(i, j, k);
0374 
0375                     if (matrix[n])
0376                     {
0377                         ofs << G4endl;
0378                         ofs << i << '\t' << j << '\t' << k << '\t';
0379 
0380                         // Total dose
0381                         //ofs << std::setw(width) << (matrix[n]/massOfVoxel)/doseUnit;
0382                         ofs << (matrix[n]/massOfVoxel)/doseUnit;
0383 
0384                         if (secondary)
0385                         {
0386                             for (size_t l=0; l < ionStore.size(); l++)
0387                             {
0388                                 // Fill ASCII file rows
0389                                 //ofs << std::setw(width) << ionStore[l].dose[n]/massOfVoxel/doseUnit <<
0390                                 //       std::setw(width) << ionStore[l].fluence[n];
0391 
0392                                  ofs << '\t' << ionStore[l].dose[n]/massOfVoxel/doseUnit <<
0393                                      '\t' << ionStore[l].fluence[n];
0394                             }
0395                         }
0396                     }
0397                 }
0398         ofs.close();
0399     }
0400 }
0401 //////////////////////////////////////////////////////////////////////////////
0402 void HadrontherapyMatrix::Fill(G4int i, G4int j, G4int k,
0403                                G4double energyDeposit)
0404 {
0405     if (matrix)
0406         matrix[Index(i,j,k)] += energyDeposit;
0407 
0408     // Store the energy deposit in the matrix element corresponding
0409     // to the phantom voxel
0410 }
0411 
0412 
0413 
0414