cellularPhantom/src/RunAction.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 //       MONTE CARLO SIMULATION OF REALISTIC GEOMETRY FROM MICROSCOPES IMAGES
0028 //
0029 // Authors and contributors:
0030 // P. Barberet (a), S. Incerti (a), N. H. Tran (a), L. Morelli (a,b)
0031 //
0032 // a) University of Bordeaux, CNRS, LP2i, UMR5797, Gradignan, France
0033 // b) Politecnico di Milano, Italy
0034 //
0035 // If you use this code, please cite the following publication:
0036 // P. Barberet et al.,
0037 // "Monte-Carlo dosimetry on a realistic cell monolayer
0038 // geometry exposed to alpha particles."
0039 // Ph. Barberet et al 2012 Phys. Med. Biol. 57 2189
0040 // doi: 110.1088/0031-9155/57/8/2189
0041 // -----------------------------------------------------------------------------
0042
0043 #include "RunAction.hh"
0044
0045 #include "G4UnitsTable.hh"
0046
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0048
0049 RunAction::RunAction()
0050 :G4UserRunAction()
0051 {
0052   if (isMaster)
0053   {
0054     // Declare ntuples
0055     auto man = G4AnalysisManager::Instance();
0056     man->SetDefaultFileType("root");
0057     man->SetFirstNtupleId(1);
0058
0059     // Create 1st ntuple (id = 1)
0060     man->CreateNtuple("ntuple1", "RED");
0061     man->CreateNtupleDColumn("x");
0062     man->CreateNtupleDColumn("y");
0063     man->CreateNtupleDColumn("z");
0064     man->CreateNtupleDColumn("energy");
0065     man->CreateNtupleDColumn("dose");
0066     man->CreateNtupleIColumn("voxelID");
0067     man->FinishNtuple();
0068
0069     // Create 2nd ntuple (id = 2)
0070     man->CreateNtuple("ntuple2", "GREEN");
0071     man->CreateNtupleDColumn("x");
0072     man->CreateNtupleDColumn("y");
0073     man->CreateNtupleDColumn("z");
0074     man->CreateNtupleDColumn("energy");
0075     man->CreateNtupleDColumn("dose");
0076     man->CreateNtupleIColumn("voxelID");
0077     man->FinishNtuple();
0078
0079     // Create 3rd ntuple (id = 3)
0080     man->CreateNtuple("ntuple3", "BLUE");
0081     man->CreateNtupleDColumn("x");
0082     man->CreateNtupleDColumn("y");
0083     man->CreateNtupleDColumn("z");
0084     man->CreateNtupleDColumn("energy");
0085     man->CreateNtupleDColumn("dose");
0086     man->CreateNtupleIColumn("voxelID");
0087     man->FinishNtuple();
0088   }
0089 }
0090
0091 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0092
0093 RunAction::~RunAction()
0094 {}
0095
0096 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0097
0098 G4Run* RunAction::GenerateRun()
0099 {
0100   fRun = new Run();
0101   return fRun;
0102 }
0103
0104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0105
0106 void RunAction::BeginOfRunAction(const G4Run *)
0107 {
0108   if (isMaster)
0109   {
0110     // Analysis manager
0111     auto man = G4AnalysisManager::Instance();
0112     man->OpenFile("phantom");
0113
0114     // Access phantom singleton
0115     fMyPhantomParam = CellParameterisation::Instance();
0116   }
0117 }
0118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0119
0120 void RunAction::EndOfRunAction(const G4Run * /*aRun*/)
0121 {
0122   if (isMaster)
0123   {
0124     // Display results from merged local runs
0125     fRun->EndOfRun();
0126
0127     // Fill ntuples
0128     auto man = G4AnalysisManager::Instance();
0129
0130     G4double X, Y, Z;
0131
0132     // Total mass of voxels
0133     G4double redMassTot=0.;
0134     G4double greenMassTot=0.;
0135     G4double blueMassTot=0.;
0136
0137     redMassTot = fMyPhantomParam->GetRedMass();
0138     greenMassTot = fMyPhantomParam->GetGreenMass();
0139     blueMassTot = fMyPhantomParam->GetBlueMass();
0140
0141     // (Optional) Numbers of voxel
0142     //G4double redVox=0;
0143     //G4double greenVox=0;
0144     //G4double blueVox=0;
0145     //redVox = fMyPhantomParam->GetRedTotalPixels();
0146     //greenVox = fMyPhantomParam->GetGreenTotalPixels();
0147     //blueVox = fMyPhantomParam->GetBlueTotalPixels();
0148
0149     // (Optional) Single voxel mass
0150     //G4double redMass=0.;
0151     //G4double greenMass=0.;
0152     //G4double blueMass=0.;
0153     //redMass = redMassTot/redVox;
0154     //greenMass = greenMassTot/greenVox;
0155     //blueMass = blueMassTot/blueVox;
0156
0157     // Save x, y, z and energy for every voxel having absorbed an energy above 0.
0158     // Energy is in keV
0159     // Dose is in Gy
0160
0161     G4double cumulatedDeposit = 0;
0162
0163     // Loop on voxels and collect energy and dose from merged local runs
0164     for (G4int i = 0; i < fMyPhantomParam->GetPhantomTotalPixels(); i++)
0165     {
0166       cumulatedDeposit = fRun->GetVoxelEdeposit(i);
0167
0168       if (cumulatedDeposit > 0.)
0169       {
0170         X = (fMyPhantomParam->GetVoxelThreeVectorOriginal(i).x()) / um;
0171         Y = (fMyPhantomParam->GetVoxelThreeVectorOriginal(i).y()) / um;
0172         Z = (fMyPhantomParam->GetVoxelThreeVectorOriginal(i).z()) / um;
0173
0174         if (fMyPhantomParam->GetMaterial(i) == 1)
0175         {
0176           man->FillNtupleDColumn(1,0,X);
0177           man->FillNtupleDColumn(1,1,Y);
0178           man->FillNtupleDColumn(1,2,Z);
0179           man->FillNtupleDColumn(1,3,cumulatedDeposit/keV);
0180           man->FillNtupleDColumn(1,4,((cumulatedDeposit/joule)/(redMassTot/kg)));
0181           man->FillNtupleIColumn(1,5,i);
0182           man->AddNtupleRow(1);
0183         }
0184         else if (fMyPhantomParam->GetMaterial(i) == 2)
0185         {
0186           man->FillNtupleDColumn(2,0,X);
0187           man->FillNtupleDColumn(2,1,Y);
0188           man->FillNtupleDColumn(2,2,Z);
0189           man->FillNtupleDColumn(2,3,cumulatedDeposit/keV);
0190           man->FillNtupleDColumn(2,4,((cumulatedDeposit/joule)/(greenMassTot/kg)));
0191           man->FillNtupleIColumn(2,5,i);
0192           man->AddNtupleRow(2);
0193         }
0194         else if (fMyPhantomParam->GetMaterial(i) == 3)
0195         {
0196           man->FillNtupleDColumn(3,0,X);
0197           man->FillNtupleDColumn(3,1,Y);
0198           man->FillNtupleDColumn(3,2,Z);
0199           man->FillNtupleDColumn(3,3,cumulatedDeposit/keV);
0200           man->FillNtupleDColumn(3,4,((cumulatedDeposit/joule)/(blueMassTot/kg)));
0201           man->FillNtupleIColumn(3,5,i);
0202           man->AddNtupleRow(3);
0203         }
0204       }
0205     }
0206
0207     // Save histograms
0208     man->Write();
0209     man->CloseFile();
0210     man->Clear();
0211   }
0212 }