Ignore:
Timestamp:
Jun 14, 2010, 3:54:58 PM (15 years ago)
Author:
garnier
Message:

geant4.9.4 beta rc0

Location:
trunk/examples/advanced/hadrontherapy/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/examples/advanced/hadrontherapy/src/HadrontherapyAnalysisManager.cc

    r1230 r1313  
    1 //
    2 // ********************************************************************
    3 // * License and Disclaimer                                           *
    4 // *                                                                  *
    5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
    6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
    7 // * conditions of the Geant4 Software License,  included in the file *
    8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
    9 // * include a list of copyright holders.                             *
    10 // *                                                                  *
    11 // * Neither the authors of this software system, nor their employing *
    12 // * institutes,nor the agencies providing financial support for this *
    13 // * work  make  any representation or  warranty, express or implied, *
    14 // * regarding  this  software system or assume any liability for its *
    15 // * use.  Please see the license in the file  LICENSE  and URL above *
    16 // * for the full disclaimer and the limitation of liability.         *
    17 // *                                                                  *
    18 // * This  code  implementation is the result of  the  scientific and *
    19 // * technical work of the GEANT4 collaboration.                      *
    20 // * By using,  copying,  modifying or  distributing the software (or *
    21 // * any work based  on the software)  you  agree  to acknowledge its *
    22 // * use  in  resulting  scientific  publications,  and indicate your *
    23 // * acceptance of all terms of the Geant4 Software license.          *
    24 // ********************************************************************
    25 //
    26 // $Id: HadrontherapyAnalisysManager.cc;
    27 // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
     1        //
     2        // ********************************************************************
     3        // * License and Disclaimer                                           *
     4        // *                                                                  *
     5        // * The  Geant4 software  is  copyright of the Copyright Holders  of *
     6        // * the Geant4 Collaboration.  It is provided  under  the terms  and *
     7        // * conditions of the Geant4 Software License,  included in the file *
     8        // * LICENSE and available at  http://cern.ch/geant4/license .  These *
     9        // * include a list of copyright holders.                             *
     10        // *                                                                  *
     11        // * Neither the authors of this software system, nor their employing *
     12        // * institutes,nor the agencies providing financial support for this *
     13        // * work  make  any representation or  warranty, express or implied, *
     14        // * regarding  this  software system or assume any liability for its *
     15        // * use.  Please see the license in the file  LICENSE  and URL above *
     16        // * for the full disclaimer and the limitation of liability.         *
     17        // *                                                                  *
     18        // * This  code  implementation is the result of  the  scientific and *
     19        // * technical work of the GEANT4 collaboration.                      *
     20        // * By using,  copying,  modifying or  distributing the software (or *
     21        // * any work based  on the software)  you  agree  to acknowledge its *
     22        // * use  in  resulting  scientific  publications,  and indicate your *
     23        // * acceptance of all terms of the Geant4 Software license.          *
     24        // ********************************************************************
     25        //
     26        // $Id: HadrontherapyAnalisysManager.cc;
     27        // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
    2828
    2929#include "HadrontherapyAnalysisManager.hh"
    … …  
    3131#include "HadrontherapyAnalysisFileMessenger.hh"
    3232#include <time.h>
    33 #ifdef ANALYSIS_USE
     33
    3434HadrontherapyAnalysisManager* HadrontherapyAnalysisManager::instance = 0;
    3535
    36 #ifdef G4ANALYSIS_USE_ROOT
    37 #undef G4ANALYSIS_USE
    38 #endif
    39 
    40 /////////////////////////////////////////////////////////////////////////////
    41 
    42 #ifdef G4ANALYSIS_USE
    43 HadrontherapyAnalysisManager::HadrontherapyAnalysisManager() :
    44   analysisFileName("DoseDistribution.root"), aFact(0), theTree(0), histFact(0), tupFact(0), h1(0), h2(0), h3(0),
    45   h4(0), h5(0), h6(0), h7(0), h8(0), h9(0), h10(0), h11(0), h12(0), h13(0), h14(0), h15(0), h16(0), ntuple(0),
    46   ionTuple(0),
    47   fragmentTuple(0),
    48   eventCounter(0)
    49 {
    50         fMess = new HadrontherapyAnalysisFileMessenger(this);
    51 }
    52 #endif
    53 #ifdef G4ANALYSIS_USE_ROOT
    54 HadrontherapyAnalysisManager::HadrontherapyAnalysisManager() :
    55   analysisFileName("DoseDistribution.root"),theTFile(0), histo1(0), histo2(0), histo3(0),
    56   histo4(0), histo5(0), histo6(0), histo7(0), histo8(0), histo9(0), histo10(0), histo11(0), histo12(0), histo13(0), histo14(0), histo15(0), histo16(0),
    57   theROOTNtuple(0),
    58   theROOTIonTuple(0),
    59   fragmentNtuple(0),
    60   metaData(0),
    61   eventCounter(0)
    62 {
    63   fMess = new HadrontherapyAnalysisFileMessenger(this);
     36HadrontherapyAnalysisManager::HadrontherapyAnalysisManager():
     37#ifdef G4ANALYSIS_USE_ROOT
     38analysisFileName("DoseDistribution.root"),theTFile(0), histo1(0), histo2(0), histo3(0),
     39histo4(0), histo5(0), histo6(0), histo7(0), histo8(0), histo9(0), histo10(0), histo11(0), histo12(0), histo13(0), histo14(0), histo15(0), histo16(0),
     40kinFragNtuple(0),
     41kineticEnergyPrimaryNtuple(0),
     42doseFragNtuple(0),
     43fluenceFragNtuple(0),
     44letFragNtuple(0),
     45theROOTNtuple(0),
     46theROOTIonTuple(0),
     47fragmentNtuple(0),
     48metaData(0),
     49eventCounter(0)
     50{
     51    fMess = new HadrontherapyAnalysisFileMessenger(this);
    6452}
    6553#endif
    6654/////////////////////////////////////////////////////////////////////////////
     55
    6756HadrontherapyAnalysisManager::~HadrontherapyAnalysisManager()
    6857{
    69 delete(fMess); //kill the messenger
    70 #ifdef G4ANALYSIS_USE
    71   delete fragmentTuple;
    72   fragmentTuple = 0;
    73 
    74   delete ionTuple;
    75   ionTuple = 0;
    76 
    77   delete ntuple;
    78   ntuple = 0;
    79 
    80   delete h16;
    81   h16 = 0;
    82 
    83   delete h15;
    84   h15 = 0;
    85 
    86   delete h14;
    87   h14 = 0;
    88 
    89   delete h13;
    90   h13 = 0;
    91 
    92   delete h12;
    93   h12 = 0;
    94 
    95   delete h11;
    96   h11 = 0;
    97 
    98   delete h10;
    99   h10 = 0;
    100 
    101   delete h9;
    102   h9 = 0;
    103 
    104   delete h8;
    105   h8 = 0;
    106 
    107   delete h7;
    108   h7 = 0;
    109 
    110   delete h6;
    111   h6 = 0;
    112 
    113   delete h5;
    114   h5 = 0;
    115 
    116   delete h4;
    117   h4 = 0;
    118 
    119   delete h3;
    120   h3 = 0;
    121 
    122   delete h2;
    123   h2 = 0;
    124 
    125   delete h1;
    126   h1 = 0;
    127 
    128   delete tupFact;
    129   tupFact = 0;
    130 
    131   delete histFact;
    132   histFact = 0;
    133 
    134   delete theTree;
    135   theTree = 0;
    136 
    137   delete aFact;
    138   aFact = 0;
     58    delete(fMess);
     59#ifdef G4ANALYSIS_USE_ROOT
     60    Clear();   
    13961#endif
     62}
     63
     64HadrontherapyAnalysisManager* HadrontherapyAnalysisManager::GetInstance()
     65{
     66        if (instance == 0) instance = new HadrontherapyAnalysisManager;
     67        return instance;
     68}
    14069#ifdef G4ANALYSIS_USE_ROOT
    141   delete metaData;
    142   metaData = 0;
    143 
    144   delete fragmentNtuple;
    145   fragmentNtuple = 0;
    146 
    147   delete theROOTIonTuple;
    148   theROOTIonTuple = 0;
    149 
    150   delete theROOTNtuple;
    151   theROOTNtuple = 0;
    152  
    153   delete histo16;
    154   histo14 = 0;
    155 
    156   delete histo15;
    157   histo14 = 0;
    158 
    159   delete histo14;
    160   histo14 = 0;
    161 
    162   delete histo13;
    163   histo13 = 0;
    164 
    165   delete histo12;
    166   histo12 = 0;
    167 
    168   delete histo11;
    169   histo11 = 0;
    170 
    171   delete histo10;
    172   histo10 = 0;
    173 
    174   delete histo9;
    175   histo9 = 0;
    176 
    177   delete histo8;
    178   histo8 = 0;
    179 
    180   delete histo7;
    181   histo7 = 0;
    182 
    183   delete histo6;
    184   histo6 = 0;
    185 
    186   delete histo5;
    187   histo5 = 0;
    188 
    189   delete histo4;
    190   histo4 = 0;
    191 
    192   delete histo3;
    193   histo3 = 0;
    194 
    195   delete histo2;
    196   histo2 = 0;
    197 
    198   delete histo1;
    199   histo1 = 0;
     70void HadrontherapyAnalysisManager::Clear()
     71{
     72    // XXX delete ALL variables!!
     73    // but this seems to be automatically done by
     74    // theTFile->Close() command!
     75    // so we get a *double free or corruption*
     76
     77    delete metaData;
     78    metaData = 0;
     79
     80    delete fragmentNtuple;
     81    fragmentNtuple = 0;
     82
     83    delete theROOTIonTuple;
     84    theROOTIonTuple = 0;
     85
     86    delete theROOTNtuple;
     87    theROOTNtuple = 0;
     88
     89    delete histo16;
     90    histo16 = 0;
     91
     92    delete histo15;
     93    histo15 = 0;
     94
     95    delete histo14;
     96    histo14 = 0;
     97
     98    delete histo13;
     99    histo13 = 0;
     100
     101    delete histo12;
     102    histo12 = 0;
     103
     104    delete histo11;
     105    histo11 = 0;
     106
     107    delete histo10;
     108    histo10 = 0;
     109
     110    delete histo9;
     111    histo9 = 0;
     112
     113    delete histo8;
     114    histo8 = 0;
     115
     116    delete histo7;
     117    histo7 = 0;
     118
     119    delete histo6;
     120    histo6 = 0;
     121
     122    delete histo5;
     123    histo5 = 0;
     124
     125    delete histo4;
     126    histo4 = 0;
     127
     128    delete histo3;
     129    histo3 = 0;
     130
     131    delete histo2;
     132    histo2 = 0;
     133
     134    delete histo1;
     135    histo1 = 0;
     136}
     137/////////////////////////////////////////////////////////////////////////////
     138
     139void HadrontherapyAnalysisManager::SetAnalysisFileName(G4String aFileName)
     140{
     141        this->analysisFileName = aFileName;
     142}
     143
     144        /////////////////////////////////////////////////////////////////////////////
     145void HadrontherapyAnalysisManager::book()
     146{
     147        delete theTFile; // this is similar to theTFile->Close() => delete all associated variables created via new, moreover it delete itself. 
     148
     149        theTFile = new TFile(analysisFileName, "RECREATE");
     150
     151        // Create the histograms with the energy deposit along the X axis
     152        histo1 = createHistogram1D("braggPeak","slice, energy", 400, 0., 80); //<different waterthicknesses are accoutned for in ROOT-analysis stage
     153        histo2 = createHistogram1D("h20","Secondary protons - slice, energy", 400, 0., 400.);
     154        histo3 = createHistogram1D("h30","Secondary neutrons - slice, energy", 400, 0., 400.);
     155        histo4 = createHistogram1D("h40","Secondary alpha - slice, energy", 400, 0., 400.);
     156        histo5 = createHistogram1D("h50","Secondary gamma - slice, energy", 400, 0., 400.);
     157        histo6 = createHistogram1D("h60","Secondary electron - slice, energy", 400, 0., 400.);
     158        histo7 = createHistogram1D("h70","Secondary triton - slice, energy", 400, 0., 400.);
     159        histo8 = createHistogram1D("h80","Secondary deuteron - slice, energy", 400, 0., 400.);
     160        histo9 = createHistogram1D("h90","Secondary pion - slice, energy", 400, 0., 400.);
     161        histo10 = createHistogram1D("h100","Energy distribution of secondary electrons", 70, 0., 70.);
     162        histo11 = createHistogram1D("h110","Energy distribution of secondary photons", 70, 0., 70.);
     163        histo12 = createHistogram1D("h120","Energy distribution of secondary deuterons", 70, 0., 70.);
     164        histo13 = createHistogram1D("h130","Energy distribution of secondary tritons", 70, 0., 70.);
     165        histo14 = createHistogram1D("h140","Energy distribution of secondary alpha particles", 70, 0., 70.);
     166        histo15 = createHistogram1D("heliumEnergyAfterPhantom","Energy distribution of secondary helium fragments after the phantom",
     167                70, 0., 500.);
     168        histo16 = createHistogram1D("hydrogenEnergyAfterPhantom","Energy distribution of secondary helium fragments after the phantom",
     169                70, 0., 500.);
     170
     171        kinFragNtuple  = new TNtuple("kinFragNtuple",
     172                "Kinetic energy by voxel & fragment",
     173                "i:j:k:A:Z:kineticEnergy");
     174        kineticEnergyPrimaryNtuple= new TNtuple("kineticEnergyPrimaryNtuple",
     175                "Kinetic energy by voxel of primary",
     176                "i:j:k:kineticEnergy");
     177        doseFragNtuple = new TNtuple("doseFragNtuple",
     178                "Energy deposit by voxel & fragment",
     179                "i:j:k:A:Z:energy");
     180
     181        fluenceFragNtuple = new TNtuple("fluenceFragNtuple",
     182                "Fluence by voxel & fragment",
     183                "i:j:k:A:Z:fluence");
     184
     185        letFragNtuple = new TNtuple("letFragNtuple",
     186                "Let by voxel & fragment",
     187                "i:j:k:A:Z:letT:letD");
     188
     189        theROOTNtuple =   new TNtuple("theROOTNtuple",
     190                "Energy deposit by slice",
     191                "i:j:k:energy");
     192
     193        theROOTIonTuple = new TNtuple("theROOTIonTuple",
     194                "Generic ion information",
     195                "a:z:occupancy:energy");
     196
     197        fragmentNtuple =  new TNtuple("fragmentNtuple",
     198                "Fragments",
     199                "A:Z:energy:posX:posY:posZ");
     200
     201        metaData =        new TNtuple("metaData",
     202                "Metadata",
     203                "events:detectorDistance:waterThickness:beamEnergy:energyError:phantomCenterDistance");
     204}
     205
     206        /////////////////////////////////////////////////////////////////////////////
     207void HadrontherapyAnalysisManager::FillEnergyDeposit(G4int i,
     208                                                                                                         G4int j,
     209                                                                                                         G4int k,
     210                                                                                                         G4double energy)
     211{
     212        if (theROOTNtuple) {
     213                theROOTNtuple->Fill(i, j, k, energy);
     214        }
     215}
     216
     217        /////////////////////////////////////////////////////////////////////////////
     218void HadrontherapyAnalysisManager::BraggPeak(G4int slice, G4double energy)
     219{
     220        histo1->SetBinContent(slice, energy); //This uses setbincontent instead of fill to get labels correct
     221}
     222
     223        /////////////////////////////////////////////////////////////////////////////
     224void HadrontherapyAnalysisManager::SecondaryProtonEnergyDeposit(G4int slice, G4double energy)
     225{
     226        histo2->Fill(slice, energy);
     227}
     228
     229        /////////////////////////////////////////////////////////////////////////////
     230void HadrontherapyAnalysisManager::SecondaryNeutronEnergyDeposit(G4int slice, G4double energy)
     231{
     232        histo3->Fill(slice, energy);
     233}
     234
     235        /////////////////////////////////////////////////////////////////////////////
     236void HadrontherapyAnalysisManager::SecondaryAlphaEnergyDeposit(G4int slice, G4double energy)
     237{
     238        histo4->Fill(slice, energy);
     239}
     240
     241        /////////////////////////////////////////////////////////////////////////////
     242void HadrontherapyAnalysisManager::SecondaryGammaEnergyDeposit(G4int slice, G4double energy)
     243{
     244        histo5->Fill(slice, energy);
     245}
     246
     247        /////////////////////////////////////////////////////////////////////////////
     248void HadrontherapyAnalysisManager::SecondaryElectronEnergyDeposit(G4int slice, G4double energy)
     249{
     250        histo6->Fill(slice, energy);
     251}
     252
     253        /////////////////////////////////////////////////////////////////////////////
     254void HadrontherapyAnalysisManager::SecondaryTritonEnergyDeposit(G4int slice, G4double energy)
     255{
     256        histo7->Fill(slice, energy);
     257}
     258
     259        /////////////////////////////////////////////////////////////////////////////
     260void HadrontherapyAnalysisManager::SecondaryDeuteronEnergyDeposit(G4int slice, G4double energy)
     261{
     262        histo8->Fill(slice, energy);
     263}
     264
     265        /////////////////////////////////////////////////////////////////////////////
     266void HadrontherapyAnalysisManager::SecondaryPionEnergyDeposit(G4int slice, G4double energy)
     267{
     268        histo9->Fill(slice, energy);
     269}
     270
     271        /////////////////////////////////////////////////////////////////////////////
     272void HadrontherapyAnalysisManager::electronEnergyDistribution(G4double energy)
     273{
     274        histo10->Fill(energy);
     275}
     276
     277        /////////////////////////////////////////////////////////////////////////////
     278void HadrontherapyAnalysisManager::gammaEnergyDistribution(G4double energy)
     279{
     280        histo11->Fill(energy);
     281}
     282
     283        /////////////////////////////////////////////////////////////////////////////
     284void HadrontherapyAnalysisManager::deuteronEnergyDistribution(G4double energy)
     285{
     286        histo12->Fill(energy);
     287}
     288
     289        /////////////////////////////////////////////////////////////////////////////
     290void HadrontherapyAnalysisManager::tritonEnergyDistribution(G4double energy)
     291{
     292        histo13->Fill(energy);
     293}
     294
     295        /////////////////////////////////////////////////////////////////////////////
     296void HadrontherapyAnalysisManager::alphaEnergyDistribution(G4double energy)
     297{
     298        histo14->Fill(energy);
     299}
     300        /////////////////////////////////////////////////////////////////////////////
     301void HadrontherapyAnalysisManager::heliumEnergy(G4double secondaryParticleKineticEnergy)
     302{
     303        histo15->Fill(secondaryParticleKineticEnergy);
     304}
     305
     306        /////////////////////////////////////////////////////////////////////////////
     307void HadrontherapyAnalysisManager::hydrogenEnergy(G4double secondaryParticleKineticEnergy)
     308{
     309        histo16->Fill(secondaryParticleKineticEnergy);
     310}
     311
     312        /////////////////////////////////////////////////////////////////////////////
     313        // FillKineticFragmentTuple create an ntuple where the voxel indexs, the atomic number and mass and the kinetic
     314        // energy of all the particles interacting with the phantom, are stored
     315void HadrontherapyAnalysisManager::FillKineticFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double kinEnergy)
     316{
     317    kinFragNtuple -> Fill(i, j, k, A, Z, kinEnergy);
     318}
     319
     320        /////////////////////////////////////////////////////////////////////////////
     321        // FillKineticEnergyPrimaryNTuple creates a ntuple where the voxel indexs and the kinetic
     322        // energies of ONLY primary particles interacting with the phantom, are stored
     323void HadrontherapyAnalysisManager::FillKineticEnergyPrimaryNTuple(G4int i, G4int j, G4int k, G4double kinEnergy)
     324{
     325    kineticEnergyPrimaryNtuple -> Fill(i, j, k, kinEnergy);
     326}
     327
     328        /////////////////////////////////////////////////////////////////////////////
     329        // This function is called only if ROOT is activated.
     330        // It is called by the HadrontherapyMatric.cc class file and it is used to create two ntuples containing
     331        // the total energy deposited and the fluence values, in each voxel and per any particle (primary beam
     332        // and secondaries )
     333void HadrontherapyAnalysisManager::FillVoxelFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double energy, G4double fluence)
     334{
     335                // Fill the ntuple containing the voxel, mass and atomic number and the energy deposited
     336    doseFragNtuple ->    Fill( i, j, k, A, Z, energy );
     337       
     338                // Fill the ntuple containing the voxel, mass and atomic number and the fluence
     339        if (i==1 && Z==1) {
     340                fluenceFragNtuple -> Fill( i, j, k, A, Z, fluence );
     341
     342        }
     343}
     344
     345void HadrontherapyAnalysisManager::FillLetFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double letT, G4double letD)
     346{
     347        letFragNtuple -> Fill( i, j, k, A, Z, letT, letD);
     348
     349}
     350        /////////////////////////////////////////////////////////////////////////////
     351void HadrontherapyAnalysisManager::FillFragmentTuple(G4int A, G4double Z, G4double energy, G4double posX, G4double posY, G4double posZ)
     352{
     353        fragmentNtuple->Fill(A, Z, energy, posX, posY, posZ);
     354}
     355
     356        /////////////////////////////////////////////////////////////////////////////
     357void HadrontherapyAnalysisManager::genericIonInformation(G4int a,
     358                                                                                                                 G4double z,
     359                                                                                                                 G4int electronOccupancy,
     360                                                                                                                 G4double energy)
     361{
     362        if (theROOTIonTuple) {
     363                theROOTIonTuple->Fill(a, z, electronOccupancy, energy);
     364        }
     365}
     366
     367        /////////////////////////////////////////////////////////////////////////////
     368void HadrontherapyAnalysisManager::startNewEvent()
     369{
     370        eventCounter++;
     371}
     372        /////////////////////////////////////////////////////////////////////////////
     373void HadrontherapyAnalysisManager::setGeometryMetaData(G4double endDetectorPosition, G4double waterThickness, G4double phantomCenter)
     374{
     375        this->detectorDistance = endDetectorPosition;
     376        this->phantomDepth = waterThickness;
     377        this->phantomCenterDistance = phantomCenter;
     378}
     379void HadrontherapyAnalysisManager::setBeamMetaData(G4double meanKineticEnergy,G4double sigmaEnergy)
     380{
     381        this->beamEnergy = meanKineticEnergy;
     382        this->energyError = sigmaEnergy;
     383}
     384/////////////////////////////////////////////////////////////////////////////
     385// Flush data & close the file
     386void HadrontherapyAnalysisManager::flush()
     387{
     388    if (theTFile)
     389    {
     390        theTFile -> Write();
     391        theTFile -> Close();
     392    }
     393    eventCounter = 0;
     394}
     395
    200396#endif
    201 }
    202 /////////////////////////////////////////////////////////////////////////////
    203 HadrontherapyAnalysisManager* HadrontherapyAnalysisManager::getInstance()
    204 {
    205   if (instance == 0) instance = new HadrontherapyAnalysisManager;
    206   return instance;
    207 }
    208 
    209 /////////////////////////////////////////////////////////////////////////////
    210 void HadrontherapyAnalysisManager::SetAnalysisFileName(G4String aFileName)
    211 {
    212   this->analysisFileName = aFileName;
    213 }
    214 
    215 /////////////////////////////////////////////////////////////////////////////
    216 void HadrontherapyAnalysisManager::book()
    217 {
    218 #ifdef G4ANALYSIS_USE
    219   // Build up  the  analysis factory
    220   aFact = AIDA_createAnalysisFactory();
    221   AIDA::ITreeFactory* treeFact = aFact -> createTreeFactory();
    222 
    223   // Create the .hbk or the .root file
    224   G4String fileName = "DoseDistribution.hbk";
    225 
    226   std::string opts = "export=root";
    227 
    228   theTree = treeFact -> create(fileName,"hbook",false,true);
    229   theTree = treeFact -> create(analysisFileName,"ROOT",false,true,opts);
    230 
    231   // Factories are not "managed" by an AIDA analysis system.
    232   // They must be deleted by the AIDA user code.
    233   delete treeFact;
    234 
    235   // Create the histogram and the ntuple factory
    236   histFact = aFact -> createHistogramFactory(*theTree);
    237   tupFact = aFact -> createTupleFactory(*theTree);
    238 
    239   // Create the histograms with the enrgy deposit along the X axis
    240   h1 = histFact -> createHistogram1D("10","slice, energy", 400, 0., 400. );
    241 
    242   h2 = histFact -> createHistogram1D("20","Secondary protons - slice, energy", 400, 0., 400. );
    243 
    244   h3 = histFact -> createHistogram1D("30","Secondary neutrons - slice, energy", 400, 0., 400. );
    245 
    246   h4 = histFact -> createHistogram1D("40","Secondary alpha - slice, energy", 400, 0., 400. );
    247 
    248   h5 = histFact -> createHistogram1D("50","Secondary gamma - slice, energy", 400, 0., 400. );
    249 
    250   h6 = histFact -> createHistogram1D("60","Secondary electron - slice, energy", 400, 0., 400. );
    251 
    252   h7 = histFact -> createHistogram1D("70","Secondary triton - slice, energy", 400, 0., 400. );
    253 
    254   h8 = histFact -> createHistogram1D("80","Secondary deuteron - slice, energy", 400, 0., 400. );
    255 
    256   h9 = histFact -> createHistogram1D("90","Secondary pion - slice, energy", 400, 0., 400. );
    257 
    258   h10 = histFact -> createHistogram1D("100","Energy distribution of secondary electrons", 70, 0., 70. );
    259 
    260   h11 = histFact -> createHistogram1D("110","Energy distribution of secondary photons", 70, 0., 70. );
    261 
    262   h12 = histFact -> createHistogram1D("120","Energy distribution of secondary deuterons", 70, 0., 70. );
    263 
    264   h13 = histFact -> createHistogram1D("130","Energy distribution of secondary tritons", 70, 0., 70. );
    265 
    266   h14 = histFact -> createHistogram1D("140","Energy distribution of secondary alpha particles", 70, 0., 70. );
    267 
    268   h15 = histFact -> createHistogram1D("150","Energy distribution of helium fragments after the phantom", 70, 0., 500.);
    269 
    270   h16 = histFact -> createHistogram1D("160","Energy distribution of hydrogen fragments after the phantom", 70, 0., 500.);
    271 
    272   // Create the ntuple
    273   G4String columnNames = "int i; int j; int k; double energy;";
    274   G4String options = "";
    275   if (tupFact) ntuple = tupFact -> create("1","1",columnNames, options);
    276 
    277   // Create the ntuple
    278   G4String columnNames2 = "int a; double z;  int occupancy; double energy;";
    279   G4String options2 = "";
    280   if (tupFact) ionTuple = tupFact -> create("2","2", columnNames2, options2);
    281 
    282   // Create the fragment ntuple
    283   G4String columnNames3 = "int a; double z; double energy; double posX; double posY; double posZ;";
    284   G4String options3 = "";
    285   if (tupFact) fragmentTuple = tupFact -> create("3","3", columnNames3, options3);
    286 #endif
    287 #ifdef G4ANALYSIS_USE_ROOT
    288   // Use ROOT
    289   theTFile = new TFile(analysisFileName, "RECREATE");
    290 
    291   // Create the histograms with the energy deposit along the X axis
    292   histo1 = createHistogram1D("braggPeak","slice, energy", 400, 0., 27.9); //<different waterthicknesses are accoutned for in ROOT-analysis stage
    293   histo2 = createHistogram1D("h20","Secondary protons - slice, energy", 400, 0., 400.);
    294   histo3 = createHistogram1D("h30","Secondary neutrons - slice, energy", 400, 0., 400.);
    295   histo4 = createHistogram1D("h40","Secondary alpha - slice, energy", 400, 0., 400.);
    296   histo5 = createHistogram1D("h50","Secondary gamma - slice, energy", 400, 0., 400.);
    297   histo6 = createHistogram1D("h60","Secondary electron - slice, energy", 400, 0., 400.);
    298   histo7 = createHistogram1D("h70","Secondary triton - slice, energy", 400, 0., 400.);
    299   histo8 = createHistogram1D("h80","Secondary deuteron - slice, energy", 400, 0., 400.);
    300   histo9 = createHistogram1D("h90","Secondary pion - slice, energy", 400, 0., 400.);
    301   histo10 = createHistogram1D("h100","Energy distribution of secondary electrons", 70, 0., 70.);
    302   histo11 = createHistogram1D("h110","Energy distribution of secondary photons", 70, 0., 70.);
    303   histo12 = createHistogram1D("h120","Energy distribution of secondary deuterons", 70, 0., 70.);
    304   histo13 = createHistogram1D("h130","Energy distribution of secondary tritons", 70, 0., 70.);
    305   histo14 = createHistogram1D("h140","Energy distribution of secondary alpha particles", 70, 0., 70.);
    306   histo15 = createHistogram1D("heliumEnergyAfterPhantom","Energy distribution of secondary helium fragments after the phantom",
    307                            70, 0., 500.);
    308   histo16 = createHistogram1D("hydrogenEnergyAfterPhantom","Energy distribution of secondary helium fragments after the phantom",
    309                            70, 0., 500.);
    310 
    311   theROOTNtuple = new TNtuple("theROOTNtuple", "Energy deposit by slice", "i:j:k:energy");
    312   theROOTIonTuple = new TNtuple("theROOTIonTuple", "Generic ion information", "a:z:occupancy:energy");
    313   fragmentNtuple = new TNtuple("fragmentNtuple", "Fragments", "A:Z:energy:posX:posY:posZ");
    314   metaData = new TNtuple("metaData", "Metadata", "events:detectorDistance:waterThickness:beamEnergy:energyError:phantomCenterDistance");
    315 #endif
    316 }
    317 
    318 /////////////////////////////////////////////////////////////////////////////
    319 void HadrontherapyAnalysisManager::FillEnergyDeposit(G4int i,
    320                                                      G4int j,
    321                                                      G4int k,
    322                                                      G4double energy)
    323 {
    324 #ifdef G4ANALYSIS_USE
    325  if (ntuple)     {
    326        G4int iSlice = ntuple -> findColumn("i");
    327        G4int jSlice = ntuple -> findColumn("j");
    328        G4int kSlice = ntuple -> findColumn("k");
    329        G4int iEnergy = ntuple -> findColumn("energy");
    330 
    331        ntuple -> fill(iSlice,i);
    332        ntuple -> fill(jSlice,j);
    333        ntuple -> fill(kSlice,k);
    334        ntuple -> fill(iEnergy, energy);     }
    335 
    336   ntuple -> addRow();
    337 #endif
    338 #ifdef G4ANALYSIS_USE_ROOT
    339   if (theROOTNtuple) {
    340     theROOTNtuple->Fill(i, j, k, energy);
    341   }
    342 #endif
    343 }
    344 
    345 /////////////////////////////////////////////////////////////////////////////
    346 void HadrontherapyAnalysisManager::BraggPeak(G4int slice, G4double energy)
    347 {
    348 #ifdef G4ANALYSIS_USE
    349   h1 -> fill(slice,energy);
    350 #endif
    351 #ifdef G4ANALYSIS_USE_ROOT
    352   histo1->SetBinContent(slice, energy); //This uses setbincontent instead of fill to get labels correct
    353 #endif
    354 }
    355 
    356 /////////////////////////////////////////////////////////////////////////////
    357 void HadrontherapyAnalysisManager::SecondaryProtonEnergyDeposit(G4int slice, G4double energy)
    358 {
    359 #ifdef G4ANALYSIS_USE
    360   h2 -> fill(slice,energy);
    361 #endif
    362 #ifdef G4ANALYSIS_USE_ROOT
    363   histo2->Fill(slice, energy);
    364 #endif
    365 }
    366 
    367 /////////////////////////////////////////////////////////////////////////////
    368 void HadrontherapyAnalysisManager::SecondaryNeutronEnergyDeposit(G4int slice, G4double energy)
    369 {
    370 #ifdef G4ANALYSIS_USE
    371   h3 -> fill(slice,energy);
    372 #endif
    373 #ifdef G4ANALYSIS_USE_ROOT
    374   histo3->Fill(slice, energy);
    375 #endif
    376 }
    377 
    378 /////////////////////////////////////////////////////////////////////////////
    379 void HadrontherapyAnalysisManager::SecondaryAlphaEnergyDeposit(G4int slice, G4double energy)
    380 {
    381 #ifdef G4ANALYSIS_USE
    382   h4 -> fill(slice,energy);
    383 #endif
    384 #ifdef G4ANALYSIS_USE_ROOT
    385   histo4->Fill(slice, energy);
    386 #endif
    387 }
    388 
    389 /////////////////////////////////////////////////////////////////////////////
    390 void HadrontherapyAnalysisManager::SecondaryGammaEnergyDeposit(G4int slice, G4double energy)
    391 {
    392 #ifdef G4ANALYSIS_USE
    393   h5 -> fill(slice,energy);
    394 #endif
    395 #ifdef G4ANALYSIS_USE_ROOT
    396   histo5->Fill(slice, energy);
    397 #endif
    398 }
    399 
    400 /////////////////////////////////////////////////////////////////////////////
    401 void HadrontherapyAnalysisManager::SecondaryElectronEnergyDeposit(G4int slice, G4double energy)
    402 {
    403 #ifdef G4ANALYSIS_USE
    404   h6 -> fill(slice,energy);
    405 #endif
    406 #ifdef G4ANALYSIS_USE_ROOT
    407   histo6->Fill(slice, energy);
    408 #endif
    409 }
    410 
    411 /////////////////////////////////////////////////////////////////////////////
    412 void HadrontherapyAnalysisManager::SecondaryTritonEnergyDeposit(G4int slice, G4double energy)
    413 {
    414 #ifdef G4ANALYSIS_USE
    415   h7 -> fill(slice,energy);
    416 #endif
    417 #ifdef G4ANALYSIS_USE_ROOT
    418   histo7->Fill(slice, energy);
    419 #endif
    420 }
    421 
    422 /////////////////////////////////////////////////////////////////////////////
    423 void HadrontherapyAnalysisManager::SecondaryDeuteronEnergyDeposit(G4int slice, G4double energy)
    424 {
    425 #ifdef G4ANALYSIS_USE
    426   h8 -> fill(slice,energy);
    427 #endif
    428 #ifdef G4ANALYSIS_USE_ROOT
    429   histo8->Fill(slice, energy);
    430 #endif
    431 }
    432 
    433 /////////////////////////////////////////////////////////////////////////////
    434 void HadrontherapyAnalysisManager::SecondaryPionEnergyDeposit(G4int slice, G4double energy)
    435 {
    436 #ifdef G4ANALYSIS_USE
    437   h9 -> fill(slice,energy);
    438 #endif
    439 #ifdef G4ANALYSIS_USE_ROOT
    440   histo9->Fill(slice, energy);
    441 #endif
    442 }
    443 
    444 /////////////////////////////////////////////////////////////////////////////
    445 void HadrontherapyAnalysisManager::electronEnergyDistribution(G4double energy)
    446 {
    447 #ifdef G4ANALYSIS_USE
    448   h10 -> fill(energy);
    449 #endif
    450 #ifdef G4ANALYSIS_USE_ROOT
    451   histo10->Fill(energy);
    452 #endif
    453 }
    454 
    455 /////////////////////////////////////////////////////////////////////////////
    456 void HadrontherapyAnalysisManager::gammaEnergyDistribution(G4double energy)
    457 {
    458 #ifdef G4ANALYSIS_USE
    459   h11 -> fill(energy);
    460 #endif
    461 #ifdef G4ANALYSIS_USE_ROOT
    462   histo11->Fill(energy);
    463 #endif
    464 }
    465 
    466 /////////////////////////////////////////////////////////////////////////////
    467 void HadrontherapyAnalysisManager::deuteronEnergyDistribution(G4double energy)
    468 {
    469 #ifdef G4ANALYSIS_USE
    470   h12 -> fill(energy);
    471 #endif
    472 #ifdef G4ANALYSIS_USE_ROOT
    473   histo12->Fill(energy);
    474 #endif
    475 }
    476 
    477 /////////////////////////////////////////////////////////////////////////////
    478 void HadrontherapyAnalysisManager::tritonEnergyDistribution(G4double energy)
    479 {
    480 #ifdef G4ANALYSIS_USE
    481   h13 -> fill(energy);
    482 #endif
    483 #ifdef G4ANALYSIS_USE_ROOT
    484   histo13->Fill(energy);
    485 #endif
    486 }
    487 
    488 /////////////////////////////////////////////////////////////////////////////
    489 void HadrontherapyAnalysisManager::alphaEnergyDistribution(G4double energy)
    490 {
    491 #ifdef G4ANALYSIS_USE
    492   h14 -> fill(energy);
    493 #endif
    494 #ifdef G4ANALYSIS_USE_ROOT
    495   histo14->Fill(energy);
    496 #endif
    497 }
    498 /////////////////////////////////////////////////////////////////////////////
    499 void HadrontherapyAnalysisManager::heliumEnergy(G4double secondaryParticleKineticEnergy)
    500 {
    501 #ifdef G4ANALYSIS_USE
    502   h15->fill(secondaryParticleKineticEnergy);
    503 #endif
    504 #ifdef G4ANALYSIS_USE_ROOT
    505   histo15->Fill(secondaryParticleKineticEnergy);
    506 #endif
    507 }
    508 
    509 /////////////////////////////////////////////////////////////////////////////
    510 void HadrontherapyAnalysisManager::hydrogenEnergy(G4double secondaryParticleKineticEnergy)
    511 {
    512 #ifdef G4ANALYSIS_USE
    513   h16->fill(secondaryParticleKineticEnergy);
    514 #endif
    515 #ifdef G4ANALYSIS_USE_ROOT
    516   histo16->Fill(secondaryParticleKineticEnergy);
    517 #endif
    518 }
    519 
    520 
    521 
    522 void HadrontherapyAnalysisManager::fillFragmentTuple(G4int A, G4double Z, G4double energy, G4double posX, G4double posY, G4double posZ)
    523 {
    524 #ifdef G4ANALYSIS_USE
    525   if (fragmentTuple)    {
    526        G4int aIndex = fragmentTuple -> findColumn("a");
    527        G4int zIndex = fragmentTuple -> findColumn("z");
    528        G4int energyIndex = fragmentTuple -> findColumn("energy");
    529        G4int posXIndex = fragmentTuple -> findColumn("posX");
    530        G4int posYIndex = fragmentTuple -> findColumn("posY");
    531        G4int posZIndex = fragmentTuple -> findColumn("posZ");
    532 
    533        fragmentTuple -> fill(aIndex,A);
    534        fragmentTuple -> fill(zIndex,Z);
    535        fragmentTuple -> fill(energyIndex, energy);
    536        fragmentTuple -> fill(posXIndex, posX);
    537        fragmentTuple -> fill(posYIndex, posY);
    538        fragmentTuple -> fill(posZIndex, posZ);
    539        fragmentTuple -> addRow();
    540   }
    541 #endif
    542 #ifdef G4ANALYSIS_USE_ROOT
    543   //G4cout <<" A = " << A << "  Z = " << Z << " energy = " << energy << G4endl;
    544   fragmentNtuple->Fill(A, Z, energy, posX, posY, posZ);
    545 #endif
    546 }
    547 
    548 /////////////////////////////////////////////////////////////////////////////
    549 void HadrontherapyAnalysisManager::genericIonInformation(G4int a,
    550                                                          G4double z,
    551                                                          G4int electronOccupancy,
    552                                                          G4double energy)
    553 {
    554 #ifdef G4ANALYSIS_USE
    555   if (ionTuple)    {
    556        G4int aIndex = ionTuple -> findColumn("a");
    557        G4int zIndex = ionTuple -> findColumn("z");
    558        G4int electronIndex = ionTuple -> findColumn("occupancy");
    559        G4int energyIndex = ionTuple -> findColumn("energy");
    560 
    561        ionTuple -> fill(aIndex,a);
    562       ionTuple -> fill(zIndex,z);
    563       ionTuple -> fill(electronIndex, electronOccupancy);
    564        ionTuple -> fill(energyIndex, energy);
    565      }
    566    ionTuple -> addRow();
    567 #endif
    568 #ifdef G4ANALYSIS_USE_ROOT
    569    if (theROOTIonTuple) {
    570      theROOTIonTuple->Fill(a, z, electronOccupancy, energy);
    571    }
    572 #endif
    573 }
    574 
    575 /////////////////////////////////////////////////////////////////////////////
    576 void HadrontherapyAnalysisManager::startNewEvent()
    577 {
    578   eventCounter++;
    579 }
    580 /////////////////////////////////////////////////////////////////////////////
    581 void HadrontherapyAnalysisManager::setGeometryMetaData(G4double endDetectorPosition, G4double waterThickness, G4double phantomCenter)
    582 {
    583   this->detectorDistance = endDetectorPosition;
    584   this->phantomDepth = waterThickness;
    585   this->phantomCenterDistance = phantomCenter;
    586 }
    587 void HadrontherapyAnalysisManager::setBeamMetaData(G4double meanKineticEnergy,G4double sigmaEnergy)
    588 {
    589   this->beamEnergy = meanKineticEnergy;
    590   this->energyError = sigmaEnergy;
    591 }
    592 /////////////////////////////////////////////////////////////////////////////
    593 void HadrontherapyAnalysisManager::flush()
    594 {
    595   HadrontherapyMatrix* matrix = HadrontherapyMatrix::getInstance();
    596 
    597   matrix->TotalEnergyDeposit();
    598 #ifdef G4ANALYSIS_USE
    599   theTree -> commit();
    600   theTree ->close();
    601 #endif
    602 #ifdef G4ANALYSIS_USE_ROOT
    603   metaData->Fill((Float_t) eventCounter,(Float_t) detectorDistance, (Float_t) phantomDepth, (Float_t) beamEnergy,(Float_t) energyError, (Float_t) phantomCenterDistance);
    604   metaData->Write();
    605   theROOTNtuple->Write();
    606   theROOTIonTuple->Write();
    607   fragmentNtuple->Write();
    608   theTFile->Write();
    609   //  theTFile->Clear();
    610   theTFile->Close();
    611 #endif
    612   eventCounter = 0;
    613   matrix->flush();
    614 }
    615 /////////////////////////////////////////////////////////////////////////////
    616 void HadrontherapyAnalysisManager::finish()
    617 {
    618 #ifdef G4ANALYSIS_USE
    619   // Write all histograms to file
    620   theTree -> commit();
    621   // Close (will again commit)
    622   theTree ->close();
    623 #endif
    624 #ifdef G4ANALYSIS_USE_ROOT
    625   metaData->Fill((Float_t) eventCounter,(Float_t) detectorDistance, (Float_t) phantomDepth, (Float_t) beamEnergy,(Float_t) energyError, (Float_t) phantomCenterDistance);
    626   metaData->Write();
    627   theROOTNtuple->Write();
    628   theROOTIonTuple->Write();
    629   fragmentNtuple->Write();
    630   theTFile->Write();
    631   theTFile->Close();
    632 #endif
    633   eventCounter = 0;
    634 }
    635 
    636 #endif
    637 
    638 
    639 
    640 
    641 
    642 
    643 
    644 
    645 

  • trunk/examples/advanced/hadrontherapy/src/HadrontherapyDetectorConstruction.cc

    r1230 r1313  
    2828// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
    2929
     30
    3031#include "G4SDManager.hh"
    3132#include "G4RunManager.hh"
     33#include "G4GeometryManager.hh"
     34#include "G4SolidStore.hh"
     35#include "G4PhysicalVolumeStore.hh"
     36#include "G4LogicalVolumeStore.hh"
    3237#include "G4Box.hh"
    3338#include "G4LogicalVolume.hh"
    … …  
    4247#include "G4VisAttributes.hh"
    4348#include "G4NistManager.hh"
     49
     50#include "HadrontherapyDetectorConstruction.hh"
    4451#include "HadrontherapyDetectorROGeometry.hh"
    4552#include "HadrontherapyDetectorMessenger.hh"
    4653#include "HadrontherapyDetectorSD.hh"
    47 #include "HadrontherapyDetectorConstruction.hh"
    4854#include "HadrontherapyMatrix.hh"
     55#include "HadrontherapyAnalysisManager.hh"
     56
     57#include <cmath>
    4958
    5059/////////////////////////////////////////////////////////////////////////////
    5160HadrontherapyDetectorConstruction::HadrontherapyDetectorConstruction(G4VPhysicalVolume* physicalTreatmentRoom)
    52   : motherPhys(physicalTreatmentRoom),
     61  : motherPhys(physicalTreatmentRoom), // pointer to WORLD volume
    5362    detectorSD(0), detectorROGeometry(0), matrix(0),
    54     phantomPhysicalVolume(0),
    55     detectorLogicalVolume(0), detectorPhysicalVolume(0),
    56     phantomSizeX(20.*cm), phantomSizeY(20.*cm), phantomSizeZ(20.*cm), // Default half dimensions
    57     detectorSizeX(2.*cm), detectorSizeY(2.*cm), detectorSizeZ(2.*cm),
    58     phantomPosition(20.*cm, 0.*cm, 0.*cm),
    59     detectorToPhantomPosition(0.*cm,18.*cm,18.*cm)// Default displacement of the detector respect to the phantom
    60 {
     63    phantom(0), detector(0),
     64    phantomLogicalVolume(0), detectorLogicalVolume(0),
     65    phantomPhysicalVolume(0), detectorPhysicalVolume(0),
     66    aRegion(0)
     67{
     68  HadrontherapyAnalysisManager::GetInstance();
     69
    6170  // NOTE! that the HadrontherapyDetectorConstruction class
    6271  // does NOT inherit from G4VUserDetectorConstruction G4 class
    … …  
    6978  // Default detector voxels size
    7079  // 200 slabs along the beam direction (X)
    71   sizeOfVoxelAlongX = 200 *um; 
    72   sizeOfVoxelAlongY = 2 * detectorSizeY;
    73   sizeOfVoxelAlongZ = 2 * detectorSizeZ;
    74 
    75   // Calculate (and eventually set) detector position by displacement, phantom size and detector size
    76   SetDetectorPosition();
    77 
    78   // Build phantom and associated detector
    79   ConstructPhantom();
    80   ConstructDetector();
    81   // Set number of the detector voxels along X Y and Z directions. 
    82   // This will construct also the sensitive detector, the ROGeometry
    83   // and the matrix where the energy deposited is collected!
    84   SetNumberOfVoxelBySize(sizeOfVoxelAlongX, sizeOfVoxelAlongY, sizeOfVoxelAlongZ);
     80  sizeOfVoxelAlongX = 200 *um;
     81  sizeOfVoxelAlongY = 4 *cm;
     82  sizeOfVoxelAlongZ = 4 *cm;
     83
     84  // Define here the material of the water phantom and of the detector
     85  SetPhantomMaterial("G4_WATER");
     86  // Construct geometry (messenger commands)
     87  SetDetectorSize(4.*cm, 4.*cm, 4.*cm);
     88  SetPhantomSize(40. *cm, 40. *cm, 40. *cm);
     89  SetPhantomPosition(G4ThreeVector(20. *cm, 0. *cm, 0. *cm));
     90  SetDetectorToPhantomPosition(G4ThreeVector(0. *cm, 18. *cm, 18. *cm));
     91
     92
     93  // Write virtual parameters to the real ones and check for consistency     
     94  UpdateGeometry();
    8595}
    8696
    … …  
    8898HadrontherapyDetectorConstruction::~HadrontherapyDetectorConstruction()
    8999{
    90     delete detectorROGeometry;// This should be safe in C++ even if the argument is a NULL pointer 
     100    delete detectorROGeometry; 
    91101    delete matrix; 
    92102    delete detectorMessenger;
    93103}
    94104
     105/////////////////////////////////////////////////////////////////////////////
     106// ConstructPhantom() is the method that reconstuct a water box (called phantom
     107// (or water phantom) in the usual Medical physicists slang).
     108// A water phantom can be considered a good
     109// approximation of a an human body.
    95110void HadrontherapyDetectorConstruction::ConstructPhantom()
    96111{
    97   //----------------------------------------
    98   // Phantom:
    99   // A box used to approximate tissues
    100   //----------------------------------------
    101 
    102     G4bool isotopes =  false;
    103     G4Material* waterNist = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER", isotopes);
    104     phantom = new G4Box("Phantom",phantomSizeX, phantomSizeY, phantomSizeZ);
     112    // Definition of the solid volume of the Phantom
     113    phantom = new G4Box("Phantom",
     114                        phantomSizeX/2,
     115                        phantomSizeY/2,
     116                        phantomSizeZ/2);
     117   
     118// Definition of the logical volume of the Phantom
    105119    phantomLogicalVolume = new G4LogicalVolume(phantom,
    106                                              waterNist,
     120                                             phantomMaterial,
    107121                                             "phantomLog", 0, 0, 0);
    108122 
     123    // Definition of the physics volume of the Phantom
    109124    phantomPhysicalVolume = new G4PVPlacement(0,
    110125                                            phantomPosition,
    … …  
    119134    red -> SetVisibility(true);
    120135    red -> SetForceSolid(true);
    121 //red -> SetForceWireframe(true);
     136    //red -> SetForceWireframe(true);
    122137    phantomLogicalVolume -> SetVisAttributes(red);
    123138}
    124139
    125140/////////////////////////////////////////////////////////////////////////////
     141// ConstructDetector() it the method the reconstruct a detector region
     142// inside the water phantom. It is a volume, located inside the water phantom
     143// and with two coincident faces:
     144//
     145//           **************************
     146//           *   water phantom        *
     147//           *                        *
     148//           *                        *
     149//           *---------------         *
     150//  Beam     *              -         *
     151//  ----->   * detector     -         *
     152//           *              -         *
     153//           *---------------         *
     154//           *                        *
     155//           *                        *
     156//           *                        *
     157//           **************************
     158//
     159// The detector is the volume that can be dived in slices or voxelized
     160// and in it we can collect a number of usefull information:
     161// dose distribution, fluence distribution, LET and so on
    126162void HadrontherapyDetectorConstruction::ConstructDetector()
    127163{
    128   //-----------
    129   // Detector
    130   //-----------
    131     G4bool isotopes =  false;
    132     G4Material* waterNist = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER", isotopes);
    133     detector = new G4Box("Detector",detectorSizeX,detectorSizeY,detectorSizeZ);
     164
     165    // Definition of the solid volume of the Detector
     166    detector = new G4Box("Detector",
     167                         detectorSizeX/2,
     168                         detectorSizeY/2,
     169                         detectorSizeZ/2);
     170   
     171    // Definition of the logic volume of the Phantom
    134172    detectorLogicalVolume = new G4LogicalVolume(detector,
    135                                                 waterNist,
     173                                                detectorMaterial,
    136174                                                "DetectorLog",
    137175                                                0,0,0);
    138 // Detector is attached by default to the phantom face directly exposed to the beam
    139     detectorPhysicalVolume = new G4PVPlacement(0,
    140                                              detectorPosition, // Setted by displacement
    141                                             "DetectorPhys",
    142                                              detectorLogicalVolume,
    143                                              phantomPhysicalVolume,
    144                                              false,0);
    145  
     176// Definition of the physical volume of the Phantom
     177    detectorPhysicalVolume = new G4PVPlacement(0, 
     178                                               detectorPosition, // Setted by displacement
     179                                               "DetectorPhys",
     180                                               detectorLogicalVolume,
     181                                               phantomPhysicalVolume,
     182                                               false,0);
     183
    146184// Visualisation attributes of the detector
    147185    skyBlue = new G4VisAttributes( G4Colour(135/255. , 206/255. ,  235/255. ));
    148186    skyBlue -> SetVisibility(true);
    149187    skyBlue -> SetForceSolid(true);
    150 //skyBlue -> SetForceWireframe(true);
     188    //skyBlue -> SetForceWireframe(true);
    151189    detectorLogicalVolume -> SetVisAttributes(skyBlue);
     190
     191  // **************
     192  // Cut per Region
     193  // **************
    152194 
    153 }
    154 /////////////////////////////////////////////////////////////////////////////
     195  // A smaller cut is fixed in the phantom to calculate the energy deposit with the
     196  // required accuracy
     197    if (!aRegion)
     198    {
     199        aRegion = new G4Region("DetectorLog");
     200        detectorLogicalVolume -> SetRegion(aRegion);
     201        aRegion -> AddRootLogicalVolume(detectorLogicalVolume);
     202    }
     203
     204}
     205
     206/////////////////////////////////////////////////////////////////////////////
     207
    155208void  HadrontherapyDetectorConstruction::ConstructSensitiveDetector(G4ThreeVector detectorToWorldPosition)
    156209{ 
    157210    // Install new Sensitive Detector and ROGeometry
    158211    delete detectorROGeometry; // this should be safe in C++ also if we have a NULL pointer
    159 
     212    //if (detectorSD) detectorSD->PrintAll();
     213    //delete detectorSD;
    160214    // Sensitive Detector and ReadOut geometry definition
    161215    G4SDManager* sensitiveDetectorManager = G4SDManager::GetSDMpointer();
    162216
    163     G4String sensitiveDetectorName = "Detector";
     217    static G4String sensitiveDetectorName = "Detector";
    164218    if (!detectorSD)
    165219        {
    … …  
    168222        }
    169223    // The Read Out Geometry is instantiated
    170     G4String ROGeometryName = "DetectorROGeometry";
     224    static G4String ROGeometryName = "DetectorROGeometry";
    171225    detectorROGeometry = new HadrontherapyDetectorROGeometry(ROGeometryName,
    172226                                                            detectorToWorldPosition,
    173                                                             detectorSizeX,
    174                                                             detectorSizeY,
    175                                                             detectorSizeZ,
     227                                                            detectorSizeX/2,
     228                                                            detectorSizeY/2,
     229                                                            detectorSizeZ/2,
    176230                                                            numberOfVoxelsAlongX,
    177231                                                            numberOfVoxelsAlongY,
    … …  
    193247        }
    194248}
    195 
     249void  HadrontherapyDetectorConstruction::ParametersCheck()
     250{
     251    // Check phantom/detector sizes & relative position
     252    if (!IsInside(detectorSizeX,
     253                detectorSizeY,
     254                detectorSizeZ,
     255                phantomSizeX,
     256                phantomSizeY,
     257                phantomSizeZ,
     258                detectorToPhantomPosition
     259                ))
     260        G4Exception("Error at HadrontherapyDetectorConstruction::ParametersCheck(). Detector is not fully inside Phantom!");
     261
     262    // Check Detector sizes respect to the voxel ones
     263
     264    if ( detectorSizeX < sizeOfVoxelAlongX) {
     265        G4Exception("Error at HadrontherapyDetectorConstruction::ParametersCheck(). Detector X size must be bigger or equal than that of Voxel X");
     266    }
     267    if ( detectorSizeY < sizeOfVoxelAlongY) {
     268        G4Exception("Error at HadrontherapyDetectorConstruction::ParametersCheck(). Detector Y size must be bigger or equal than that of Voxel Y");
     269    }
     270    if ( detectorSizeZ < sizeOfVoxelAlongZ) {
     271        G4Exception("Error at HadrontherapyDetectorConstruction::ParametersCheck(). Detector Z size must be bigger or equal than that of Voxel Z");
     272    }
     273
     274}
    196275/////////////////
    197276// MESSENGERS //
    198277////////////////
    199 G4bool HadrontherapyDetectorConstruction::SetNumberOfVoxelBySize(G4double sizeX, G4double sizeY, G4double sizeZ)
    200 {
    201     // Only change positive dimensions
    202     // XXX numberOfVoxels must be an integer, warn the user
    203 
    204     if (sizeX > 0)
     278
     279G4bool HadrontherapyDetectorConstruction::SetPhantomMaterial(G4String material)
     280{
     281
     282    if (G4Material* pMat = G4NistManager::Instance()->FindOrBuildMaterial(material, false) )
    205283    {
    206         if (sizeX > 2*detectorSizeX)
     284        phantomMaterial  = pMat;
     285        detectorMaterial = pMat;
     286        if (detectorLogicalVolume && phantomLogicalVolume)
    207287        {
    208             G4cout << "WARNING: Voxel X size must be smaller or equal than that of detector X" << G4endl;
    209             return false;
     288            detectorLogicalVolume -> SetMaterial(pMat);
     289            phantomLogicalVolume ->  SetMaterial(pMat);
     290
     291            G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
     292            G4RunManager::GetRunManager() -> GeometryHasBeenModified();
     293            G4cout << "The material of Phantom/Detector has been changed to " << material << G4endl;
    210294        }
    211         // Round to the nearest integer
    212         numberOfVoxelsAlongX = lrint(2 * detectorSizeX / sizeX);
    213         sizeOfVoxelAlongX = (2 * detectorSizeX / numberOfVoxelsAlongX );
    214         if(sizeOfVoxelAlongX!=sizeX) G4cout << "Rounding " <<
    215                                                 G4BestUnit(sizeX, "Length") << " to " <<
    216                                                 G4BestUnit(sizeOfVoxelAlongX, "Length") << G4endl;
    217     }
    218 
    219     if (sizeY > 0)
     295    }
     296    else
    220297    {
    221         if (sizeY > 2*detectorSizeY)
    222         {
    223             G4cout << "WARNING: Voxel Y size must be smaller or equal than that of detector Y" << G4endl;
    224             return false;
    225         }
    226         numberOfVoxelsAlongY = lrint(2 * detectorSizeY / sizeY);
    227         sizeOfVoxelAlongY = (2 * detectorSizeY / numberOfVoxelsAlongY );
    228         if(sizeOfVoxelAlongY!=sizeY) G4cout << "Rounding " <<
    229                                                 G4BestUnit(sizeY, "Length") << " to " <<
    230                                                 G4BestUnit(sizeOfVoxelAlongY, "Length") << G4endl;
    231     }
    232     if (sizeZ > 0)
     298        G4cout << "WARNING: material \"" << material << "\" doesn't exist in NIST elements/materials"
     299            " table [located in $G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc]" << G4endl;
     300        G4cout << "Use command \"/parameter/nist\" to see full materials list!" << G4endl;
     301        return false;
     302    }
     303
     304    return true;
     305}
     306/////////////////////////////////////////////////////////////////////////////
     307void HadrontherapyDetectorConstruction::SetPhantomSize(G4double sizeX, G4double sizeY, G4double sizeZ)
     308{
     309    if (sizeX > 0.) phantomSizeX = sizeX;
     310    if (sizeY > 0.) phantomSizeY = sizeY;
     311    if (sizeZ > 0.) phantomSizeZ = sizeZ;
     312}
     313/////////////////////////////////////////////////////////////////////////////
     314/////////////////////////////////////////////////////////////////////////////
     315void HadrontherapyDetectorConstruction::SetDetectorSize(G4double sizeX, G4double sizeY, G4double sizeZ)
     316{
     317    if (sizeX > 0.) {detectorSizeX = sizeX;}
     318    if (sizeY > 0.) {detectorSizeY = sizeY;}
     319    if (sizeZ > 0.) {detectorSizeZ = sizeZ;}
     320    SetVoxelSize(sizeOfVoxelAlongX, sizeOfVoxelAlongY, sizeOfVoxelAlongZ);
     321}
     322/////////////////////////////////////////////////////////////////////////////
     323
     324void HadrontherapyDetectorConstruction::SetVoxelSize(G4double sizeX, G4double sizeY, G4double sizeZ)
     325{
     326    if (sizeX > 0.) {sizeOfVoxelAlongX = sizeX;}
     327    if (sizeY > 0.) {sizeOfVoxelAlongY = sizeY;}
     328    if (sizeZ > 0.) {sizeOfVoxelAlongZ = sizeZ;}
     329}
     330void HadrontherapyDetectorConstruction::SetPhantomPosition(G4ThreeVector pos)
     331{
     332    phantomPosition = pos;
     333}
     334
     335/////////////////////////////////////////////////////////////////////////////
     336void HadrontherapyDetectorConstruction::SetDetectorToPhantomPosition(G4ThreeVector displ)
     337{
     338    detectorToPhantomPosition = displ;
     339}
     340/////////////////////////////////////////////////////////////////////////////
     341void HadrontherapyDetectorConstruction::UpdateGeometry()
     342{
     343    ParametersCheck();
     344
     345    //G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
     346    G4GeometryManager::GetInstance() -> OpenGeometry();
     347    if (phantom)
    233348    {
    234         if (sizeZ > 2*detectorSizeZ)
    235         {
    236             G4cout << "WARNING: Voxel Z size must be smaller or equal than that of detector Z" << G4endl;
    237             return false;
    238         }
    239         numberOfVoxelsAlongZ = lrint(2 * detectorSizeZ / sizeZ);
    240         sizeOfVoxelAlongZ = (2 * detectorSizeZ / numberOfVoxelsAlongZ );
    241         if(sizeOfVoxelAlongZ!=sizeZ) G4cout << "Rounding " <<
    242                                                 G4BestUnit(sizeZ, "Length") << " to " <<
    243                                                 G4BestUnit(sizeOfVoxelAlongZ, "Length") << G4endl;
    244     }
    245 
    246     G4cout << "The (X, Y, Z) sizes of the Voxels are: (" <<
    247                 G4BestUnit(sizeOfVoxelAlongX, "Length")  << ", " <<
    248                 G4BestUnit(sizeOfVoxelAlongY, "Length")  << ", " <<
    249                 G4BestUnit(sizeOfVoxelAlongZ, "Length") << ')' << G4endl;
     349            phantom -> SetXHalfLength(phantomSizeX/2);
     350            phantom -> SetYHalfLength(phantomSizeY/2);
     351            phantom -> SetZHalfLength(phantomSizeZ/2);
     352            phantomPhysicalVolume -> SetTranslation(phantomPosition);
     353    }
     354    else   ConstructPhantom();
     355
     356    // Get the center of the detector
     357    SetDetectorPosition();
     358    if (detector)
     359    {
     360                        detector -> SetXHalfLength(detectorSizeX/2);
     361                        detector -> SetYHalfLength(detectorSizeY/2);
     362                        detector -> SetZHalfLength(detectorSizeZ/2);
     363                        detectorPhysicalVolume -> SetTranslation(detectorPosition);
     364    }
     365    else    ConstructDetector();
     366   
     367// Round to nearest integer number of voxel
     368        numberOfVoxelsAlongX = lrint(detectorSizeX / sizeOfVoxelAlongX);
     369        sizeOfVoxelAlongX = ( detectorSizeX / numberOfVoxelsAlongX );
     370
     371        numberOfVoxelsAlongY = lrint(detectorSizeY / sizeOfVoxelAlongY);
     372        sizeOfVoxelAlongY = ( detectorSizeY / numberOfVoxelsAlongY );
     373
     374        numberOfVoxelsAlongZ = lrint(detectorSizeZ / sizeOfVoxelAlongZ);
     375        sizeOfVoxelAlongZ = ( detectorSizeZ / numberOfVoxelsAlongZ );
     376
     377    //G4cout << "*************** DetectorToWorldPosition " << GetDetectorToWorldPosition()/cm << "\n";
     378    ConstructSensitiveDetector(GetDetectorToWorldPosition());
     379
     380    volumeOfVoxel = sizeOfVoxelAlongX * sizeOfVoxelAlongY * sizeOfVoxelAlongZ;
     381    massOfVoxel = detectorMaterial -> GetDensity() * volumeOfVoxel;
     382    //  This will clear the existing matrix (together with all data inside it)!
     383    matrix = HadrontherapyMatrix::GetInstance(numberOfVoxelsAlongX,
     384                                              numberOfVoxelsAlongY,
     385                                              numberOfVoxelsAlongZ,
     386                                              massOfVoxel);
     387
     388    // Inform the kernel about the new geometry
     389    G4RunManager::GetRunManager() -> GeometryHasBeenModified();
     390    G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
     391
     392    PrintParameters();
     393}
     394
     395void HadrontherapyDetectorConstruction::PrintParameters()
     396{
     397
     398    G4cout << "The (X,Y,Z) dimensions of the phantom are : (" <<
     399        G4BestUnit( phantom -> GetXHalfLength()*2., "Length") << ',' <<
     400        G4BestUnit( phantom -> GetYHalfLength()*2., "Length") << ',' <<
     401        G4BestUnit( phantom -> GetZHalfLength()*2., "Length") << ')' << G4endl;
     402   
     403    G4cout << "The (X,Y,Z) dimensions of the detector are : (" <<
     404        G4BestUnit( detector -> GetXHalfLength()*2., "Length") << ',' <<
     405        G4BestUnit( detector -> GetYHalfLength()*2., "Length") << ',' <<
     406        G4BestUnit( detector -> GetZHalfLength()*2., "Length") << ')' << G4endl;
     407
     408    G4cout << "Displacement between Phantom and World is: ";
     409    G4cout << "DX= "<< G4BestUnit(phantomPosition.getX(),"Length") <<
     410        "DY= "<< G4BestUnit(phantomPosition.getY(),"Length") <<
     411        "DZ= "<< G4BestUnit(phantomPosition.getZ(),"Length") << G4endl;
     412
     413    G4cout << "The (X,Y,Z) sizes of the Voxels are: (" <<
     414        G4BestUnit(sizeOfVoxelAlongX, "Length")  << ',' <<
     415        G4BestUnit(sizeOfVoxelAlongY, "Length")  << ',' <<
     416        G4BestUnit(sizeOfVoxelAlongZ, "Length") << ')' << G4endl;
    250417
    251418    G4cout << "The number of Voxels along (X,Y,Z) is: (" <<
    252                 numberOfVoxelsAlongX  << ", " <<
    253                 numberOfVoxelsAlongY  << ", "  <<
    254                 numberOfVoxelsAlongZ  << ')' << G4endl;
    255 
    256     //  This will clear the existing matrix (together with data inside it)!
    257     matrix = HadrontherapyMatrix::getInstance(numberOfVoxelsAlongX,
    258                                               numberOfVoxelsAlongY,
    259                                               numberOfVoxelsAlongZ);
    260 
    261     // Here construct the Sensitive Detector and Read Out Geometry
    262     ConstructSensitiveDetector(GetDetectorToWorldPosition());
    263     G4RunManager::GetRunManager() -> GeometryHasBeenModified();
    264     return true;
    265 }
    266 /////////////////////////////////////////////////////////////////////////////
    267 G4bool HadrontherapyDetectorConstruction::SetDetectorSize(G4double sizeX, G4double sizeY, G4double sizeZ)
    268 {
    269 // Check that the detector stay inside the phantom
    270         if (sizeX > 0 && sizeX < sizeOfVoxelAlongX) {G4cout << "WARNING: Detector X size must be bigger than that of Voxel X" << G4endl; return false;}
    271         if (sizeY > 0 && sizeY < sizeOfVoxelAlongY) {G4cout << "WARNING: Detector Y size must be bigger than that of Voxel Y" << G4endl; return false;}
    272         if (sizeZ > 0 && sizeZ < sizeOfVoxelAlongZ) {G4cout << "WARNING: Detector Z size must be bigger than that of Voxel Z" << G4endl; return false;}
    273 
    274         if (!IsInside(sizeX/2,
    275                       sizeY/2,
    276                       sizeZ/2,
    277                       phantomSizeX,
    278                       phantomSizeY,
    279                       phantomSizeZ,
    280                       detectorToPhantomPosition))
    281         {return false;}
    282 // Negative or null values mean don't change it!
    283     if (sizeX > 0) {
    284                         detectorSizeX = sizeX/2;
    285                         detector -> SetXHalfLength(detectorSizeX);
    286                    }
    287 
    288     if (sizeY > 0) {
    289                         detectorSizeY = sizeY/2;
    290                         detector -> SetYHalfLength(detectorSizeY);
    291                    }
    292 
    293     if (sizeZ > 0) {
    294                         detectorSizeZ = sizeZ/2;
    295                         detector -> SetZHalfLength(detectorSizeZ);
    296                     }
    297 
    298 
    299     G4cout << "The (X, Y, Z) dimensions of the detector are : (" <<
    300                   G4BestUnit( detector -> GetXHalfLength()*2., "Length") << ", " <<
    301                   G4BestUnit( detector -> GetYHalfLength()*2., "Length") << ", " <<
    302                   G4BestUnit( detector -> GetZHalfLength()*2., "Length") << ')' << G4endl;
    303 // Adjust detector position
    304     SetDetectorPosition();
    305 // Adjust voxels number accordingly to new detector geometry
    306 // Matrix will be re-instantiated!
    307 // Voxels and ROGeometry must follow the detector!
    308     SetNumberOfVoxelBySize(sizeOfVoxelAlongX, sizeOfVoxelAlongY, sizeOfVoxelAlongZ);
    309     G4RunManager::GetRunManager() -> GeometryHasBeenModified();
    310     return true;
    311 }
    312 
    313 /////////////////////////////////////////////////////////////////////////////
    314 G4bool HadrontherapyDetectorConstruction::SetPhantomSize(G4double sizeX, G4double sizeY, G4double sizeZ)
    315 {
    316 
    317         if (!IsInside(detectorSizeX,
    318                                   detectorSizeY,
    319                                   detectorSizeZ,
    320                                   sizeX/2,//method parameters
    321                                   sizeY/2,
    322                                   sizeZ/2,
    323                                   detectorToPhantomPosition
    324                                   ))
    325         return false;
    326 
    327 // Only change positive dimensions
    328     if (sizeX > 0) {
    329                      phantomSizeX = sizeX/2;
    330                      phantom -> SetXHalfLength(phantomSizeX);
    331                    }
    332     if (sizeY > 0) {
    333                      phantomSizeY = sizeY/2;
    334                      phantom -> SetYHalfLength(phantomSizeY);
    335                    }
    336 
    337     if (sizeZ > 0) {
    338                      phantomSizeZ = sizeZ/2;
    339                      phantom -> SetZHalfLength(phantomSizeZ);
    340                    }
    341  
    342 
    343     G4cout << "The (X, Y, Z) dimensions of the phantom are : (" <<
    344                   G4BestUnit( phantom -> GetXHalfLength()*2., "Length") << ", " <<
    345                   G4BestUnit( phantom -> GetYHalfLength()*2., "Length") << ", " <<
    346                   G4BestUnit( phantom -> GetZHalfLength()*2., "Length") << ')' << G4endl;
    347 //G4cout << '\n' << "Coordinate volume: " << phantomPhysicalVolume -> GetTranslation() << G4endl;
    348 // Adjust detector position inside phantom
    349     SetDetectorPosition();
    350 
    351     ConstructSensitiveDetector(GetDetectorToWorldPosition());
    352     G4RunManager::GetRunManager() -> GeometryHasBeenModified();
    353     return true;
    354 }
    355 /////////////////////////////////////////////////////////////////////////////
    356 G4bool HadrontherapyDetectorConstruction::SetPhantomPosition(G4ThreeVector displacement)
    357 {
    358 // Set Phantom position respect to the World
    359 // TODO check for overlap!
    360     phantomPosition = displacement;
    361     if (phantomPhysicalVolume)
    362         {
    363             phantomPhysicalVolume -> SetTranslation(phantomPosition);
    364             G4cout << "Displacement between Phantom and World is: ";
    365             G4cout << "DX= "<< G4BestUnit(phantomPosition.getX(),"Length") << ", " <<
    366                       "DY= "<< G4BestUnit(phantomPosition.getY(),"Length") << ", " <<
    367                       "DZ= "<< G4BestUnit(phantomPosition.getZ(),"Length") << G4endl;
    368 
    369 // Redraw ROGeometry!
    370             ConstructSensitiveDetector(GetDetectorToWorldPosition());
    371             G4RunManager::GetRunManager() -> GeometryHasBeenModified();
    372         }
    373     return true;
    374 }
    375 
    376 /////////////////////////////////////////////////////////////////////////////
    377 G4bool HadrontherapyDetectorConstruction::SetDetectorToPhantomPosition(G4ThreeVector displacement)
    378 {
    379 // Ignore negative values
    380     if (displacement[0] < 0.) displacement[0] = detectorToPhantomPosition[0];
    381     if (displacement[1] < 0.) displacement[1] = detectorToPhantomPosition[1];
    382     if (displacement[2] < 0.) displacement[2] = detectorToPhantomPosition[2];
    383 
    384     if (!IsInside(detectorSizeX,
    385                   detectorSizeY,
    386                   detectorSizeZ,
    387                   phantomSizeX,
    388                   phantomSizeY,
    389                   phantomSizeZ,
    390                   displacement // method parameter!
    391                   ))
    392     {return false;}
    393     detectorToPhantomPosition = displacement;
    394 
    395 // Adjust detector position inside phantom
    396     SetDetectorPosition();
    397 
    398     ConstructSensitiveDetector(GetDetectorToWorldPosition());
    399     G4RunManager::GetRunManager() -> GeometryHasBeenModified();
    400     return true;
    401 }
     419        numberOfVoxelsAlongX  << ',' <<
     420        numberOfVoxelsAlongY  <<','  <<
     421        numberOfVoxelsAlongZ  << ')' << G4endl;
     422
     423}

  • trunk/examples/advanced/hadrontherapy/src/HadrontherapyDetectorMessenger.cc

    r1230 r1313  
    3131#include "HadrontherapyDetectorConstruction.hh"
    3232#include "G4UIdirectory.hh"
    33 #include "G4UIcmdWithADoubleAndUnit.hh"
     33#include "G4UIcmdWith3VectorAndUnit.hh"
     34#include "G4UIcmdWithoutParameter.hh"
    3435#include "G4UIcmdWithAString.hh"
    35 #include "G4UIcmdWith3VectorAndUnit.hh"
    3636
    3737/////////////////////////////////////////////////////////////////////////////
    … …  
    4949                                                "PhantomSizeAlongZ", false);
    5050    changeThePhantomSizeCmd -> SetDefaultUnit("mm");
    51     changeThePhantomSizeCmd -> SetUnitCandidates("um mm cm");
     51    changeThePhantomSizeCmd -> SetUnitCandidates("nm um mm cm");
    5252    changeThePhantomSizeCmd -> AvailableForStates(G4State_Idle);
    5353
     54
     55    // Change Phantom material
     56    changeThePhantomMaterialCmd = new G4UIcmdWithAString("/changePhantom/material", this);
     57    changeThePhantomMaterialCmd -> SetGuidance("Change the Phantom and the detector material");
     58    changeThePhantomMaterialCmd -> SetParameterName("PhantomMaterial", false);
     59    changeThePhantomMaterialCmd -> SetDefaultValue("G4_WATER");
     60    changeThePhantomMaterialCmd -> AvailableForStates(G4State_Idle);
    5461
    5562    // Change Phantom position
    … …  
    6168                                                    "PositionAlongZ", false);
    6269    changeThePhantomPositionCmd -> SetDefaultUnit("mm");
    63     changeThePhantomPositionCmd -> SetUnitCandidates("mm cm m");
     70    changeThePhantomPositionCmd -> SetUnitCandidates("um mm cm m");
    6471    changeThePhantomPositionCmd -> AvailableForStates(G4State_Idle);
    6572
     73
     74    updateCmd = new G4UIcmdWithoutParameter("/changePhantom/update",this);
     75    updateCmd->SetGuidance("Update Phantom/Detector geometry.");
     76    updateCmd->SetGuidance("This command MUST be applied before \"beamOn\" ");
     77    updateCmd->SetGuidance("if you changed geometrical value(s).");
     78    updateCmd->AvailableForStates(G4State_Idle);
    6679
    6780    //  Change detector size
    … …  
    7487    changeTheDetectorSizeCmd -> SetParameterName("DetectorSizeAlongX", "DetectorSizeAlongY", "DetectorSizeAlongZ", false);
    7588    changeTheDetectorSizeCmd -> SetDefaultUnit("mm");
    76     changeTheDetectorSizeCmd -> SetUnitCandidates("um mm cm");
     89    changeTheDetectorSizeCmd -> SetUnitCandidates("nm um mm cm");
    7790    changeTheDetectorSizeCmd -> AvailableForStates(G4State_Idle);
    7891
    … …  
    8598                                                              "DisplacementAlongZ", false);
    8699    changeTheDetectorToPhantomPositionCmd -> SetDefaultUnit("mm");
    87     changeTheDetectorToPhantomPositionCmd -> SetUnitCandidates("um mm cm");
     100    changeTheDetectorToPhantomPositionCmd -> SetUnitCandidates("nm um mm cm");
    88101    changeTheDetectorToPhantomPositionCmd -> AvailableForStates(G4State_Idle);
    89102   
    … …  
    94107    changeTheDetectorVoxelCmd -> SetParameterName("VoxelSizeAlongX", "VoxelSizeAlongY", "VoxelSizeAlongZ", false);
    95108    changeTheDetectorVoxelCmd -> SetDefaultUnit("mm");
    96     changeTheDetectorVoxelCmd -> SetUnitCandidates("um mm cm");
     109    changeTheDetectorVoxelCmd -> SetUnitCandidates("nm um mm cm");
    97110    changeTheDetectorVoxelCmd -> AvailableForStates(G4State_Idle);
    98111   }
    … …  
    104117    delete changeThePhantomSizeCmd;
    105118    delete changeThePhantomPositionCmd;
     119    delete changeThePhantomMaterialCmd;
     120    delete updateCmd;
    106121    delete changeTheDetectorDir;
    107122    delete changeTheDetectorSizeCmd;
    … …  
    124139         hadrontherapyDetector -> SetPhantomPosition(size);
    125140  }
     141  else if (command == changeThePhantomMaterialCmd)
     142  {
     143      hadrontherapyDetector -> SetPhantomMaterial(newValue);
     144  }
    126145  else if (command == changeTheDetectorSizeCmd)
    127146  {
    … …  
    137156  {
    138157        G4ThreeVector size = changeTheDetectorVoxelCmd  -> GetNew3VectorValue(newValue);
    139         hadrontherapyDetector -> SetNumberOfVoxelBySize(size.getX(),size.getY(),size.getZ());
     158        hadrontherapyDetector -> SetVoxelSize(size.getX(),size.getY(),size.getZ());
     159  }
     160  else if (command == updateCmd)
     161  {
     162      hadrontherapyDetector -> UpdateGeometry();
    140163  }
    141164}

  • trunk/examples/advanced/hadrontherapy/src/HadrontherapyEventAction.cc

    r1230 r1313  
    2424// ********************************************************************
    2525//
    26 // $Id: HadrontherapyEventAction.cc;
    27 // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
    28 
     26// $Id: HadrontherapyEventAction.cc;
     27//
     28// See more at: http://workgroup.lngs.infn.it/geant4lns
     29//
     30// ----------------------------------------------------------------------------
     31//                 GEANT 4 - Hadrontherapy example
     32// ----------------------------------------------------------------------------
     33// Code developed by:
     34//
     35// G.A.P. Cirrone(a)*, G. Candiano, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
     36//
     37// (a) Laboratori Nazionali del Sud
     38//     of the National Institute for Nuclear Physics, Catania, Italy
     39// (b) National Institute for Nuclear Physics Section of Genova, genova, Italy
     40//
     41// * cirrone@lns.infn.it
     42// --------------------------------------------------------------
    2943#include "G4Event.hh"
    3044#include "G4EventManager.hh"
    3145#include "G4HCofThisEvent.hh"
    3246#include "G4VHitsCollection.hh"
    33 #include "G4TrajectoryContainer.hh"
    34 #include "G4Trajectory.hh"
    35 #include "G4VVisManager.hh"
    3647#include "G4SDManager.hh"
    3748#include "G4VVisManager.hh"
     49
    3850#include "HadrontherapyEventAction.hh"
    3951#include "HadrontherapyDetectorHit.hh"
    … …  
    6476  //printing survey
    6577  if (evtNb%printModulo == 0)
    66     G4cout << "\n---> Begin of Event: " << evtNb << G4endl;
     78     G4cout << "\n---> Begin of Event: " << evtNb << G4endl;
    6779   
    6880  G4SDManager* pSDManager = G4SDManager::GetSDMpointer();
    6981  if(hitsCollectionID == -1)
    7082    hitsCollectionID = pSDManager -> GetCollectionID("HadrontherapyDetectorHitsCollection");
     83 
    7184}
    7285
    7386/////////////////////////////////////////////////////////////////////////////
    7487void HadrontherapyEventAction::EndOfEventAction(const G4Event* evt)
    75 {  
     88{
    7689  if(hitsCollectionID < 0)
    7790  return;
    7891  G4HCofThisEvent* HCE = evt -> GetHCofThisEvent();
    79  
     92
     93  // Clear voxels hit list
     94  HadrontherapyMatrix* matrix = HadrontherapyMatrix::GetInstance();
     95  if (matrix) matrix -> ClearHitTrack();
     96
    8097  if(HCE)
    8198  {
    … …  
    83100    if(CHC)
    84101     {
    85            matrix = HadrontherapyMatrix::getInstance(); 
    86102       if(matrix)
    87103          {
    … …  
    101117    }
    102118  }
    103   // Extract the trajectories and draw them in the visualisation
    104 
    105   if (G4VVisManager::GetConcreteInstance())
    106     {
    107       G4TrajectoryContainer * trajectoryContainer = evt -> GetTrajectoryContainer();
    108       G4int n_trajectories = 0;
    109       if (trajectoryContainer) n_trajectories = trajectoryContainer -> entries();
    110      
    111       for (G4int i = 0; i < n_trajectories; i++)
    112         {
    113           G4Trajectory* trj = (G4Trajectory*)
    114             ((*(evt -> GetTrajectoryContainer()))[i]);
    115           if(drawFlag == "all") trj -> DrawTrajectory(50);
    116           else if((drawFlag == "charged")&&(trj -> GetCharge() != 0.))
    117             trj -> DrawTrajectory(50);
    118           else if ((drawFlag == "neutral")&&(trj -> GetCharge() == 0.))
    119             trj -> DrawTrajectory(50);               
    120         }
    121     }
    122119}
    123120

  • trunk/examples/advanced/hadrontherapy/src/HadrontherapyMatrix.cc

    r1230 r1313  
    1 //
    2 // ********************************************************************
    3 // * License and Disclaimer                                           *
    4 // *                                                                  *
    5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
    6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
    7 // * conditions of the Geant4 Software License,  included in the file *
    8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
    9 // * include a list of copyright holders.                             *
    10 // *                                                                  *
    11 // * Neither the authors of this software system, nor their employing *
    12 // * institutes,nor the agencies providing financial support for this *
    13 // * work  make  any representation or  warranty, express or implied, *
    14 // * regarding  this  software system or assume any liability for its *
    15 // * use.  Please see the license in the file  LICENSE  and URL above *
    16 // * for the full disclaimer and the limitation of liability.         *
    17 // *                                                                  *
    18 // * This  code  implementation is the result of  the  scientific and *
    19 // * technical work of the GEANT4 collaboration.                      *
    20 // * By using,  copying,  modifying or  distributing the software (or *
    21 // * any work based  on the software)  you  agree  to acknowledge its *
    22 // * use  in  resulting  scientific  publications,  and indicate your *
    23 // * acceptance of all terms of the Geant4 Software license.          *
    24 // ********************************************************************
    25 //
    26 // $Id: HadrontherapyMatrix.cc;
    27 // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy//
     1        //
     2        // ********************************************************************
     3        // * License and Disclaimer                                           *
     4        // *                                                                  *
     5        // * The  Geant4 software  is  copyright of the Copyright Holders  of *
     6        // * the Geant4 Collaboration.  It is provided  under  the terms  and *
     7        // * conditions of the Geant4 Software License,  included in the file*
     8        // * LICENSE and available at  http://cern.ch/geant4/license .  These *
     9        // * include a list of copyright holders.                             *
     10        // *                                                                  *
     11        // * Neither the authors of this software system, nor their employing *
     12        // * institutes,nor the agencies providing financial support for this *
     13        // * work  make  any representation or  warranty, express or implied, *
     14        // * regarding  this  software system or assume any liability for its *
     15        // * use.  Please see the license in the file  LICENSE  and URL above *
     16        // * for the full disclaimer and the limitation of liability.         *
     17        // *                                                                  *
     18        // * This  code  implementation is the result of  the  scientific and *
     19        // * technical work of the GEANT4 collaboration.                      *
     20        // * By using,  copying,  modifying or  distributing the software (or *
     21        // * any work based  on the software)  you  agree  to acknowledge its *
     22        // * use  in  resulting  scientific  publications,  and indicate your *
     23        // * acceptance of all terms of the Geant4 Software license.          *
     24        // ********************************************************************
     25        //
     26        // $Id: HadrontherapyMatrix.cc;
     27        // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy//
    2828
    2929#include "HadrontherapyMatrix.hh"
    3030#include "HadrontherapyAnalysisManager.hh"
     31#include "G4RunManager.hh"
     32#include "HadrontherapyPrimaryGeneratorAction.hh"
     33#include "G4ParticleGun.hh"
     34
     35#include <fstream>
     36#include <unistd.h>
     37#include <iostream>
     38#include <sstream>
     39#include <iomanip>
    3140#include "globals.hh"
    32 #include <fstream>
    33 
     41
     42// Units definition: CLHEP/Units/SystemOfUnits.h
     43//
    3444HadrontherapyMatrix* HadrontherapyMatrix::instance = NULL;
    35 // Static method that only return a pointer to the matrix object
    36 HadrontherapyMatrix* HadrontherapyMatrix::getInstance()
    37 {
    38   return instance;
    39 }
    40 // This STATIC method delete (!) the old matrix and rewrite a new object returning a pointer to it
    41 HadrontherapyMatrix* HadrontherapyMatrix::getInstance(G4int voxelX, G4int voxelY, G4int voxelZ) 
    42 {
    43     if (instance) delete instance;
    44     instance = new HadrontherapyMatrix(voxelX, voxelY, voxelZ);
    45     instance -> Initialize();   
    46     return instance;
    47 }
    48 
    49 HadrontherapyMatrix::HadrontherapyMatrix(G4int voxelX, G4int voxelY, G4int voxelZ):
    50     matrix(0)
     45G4bool HadrontherapyMatrix::secondaries = false;
     46        // Only return a pointer to matrix
     47HadrontherapyMatrix* HadrontherapyMatrix::GetInstance()
     48{
     49        return instance;
     50}
     51        // This STATIC method delete (!) the old matrix and rewrite a new object returning a pointer to it
     52        // TODO A check on the parameters is required!
     53HadrontherapyMatrix* HadrontherapyMatrix::GetInstance(G4int voxelX, G4int voxelY, G4int voxelZ, G4double mass) 
     54{
     55        if (instance) delete instance;
     56        instance = new HadrontherapyMatrix(voxelX, voxelY, voxelZ, mass);
     57        instance -> Initialize();// XXX
     58        return instance;
     59}
     60HadrontherapyMatrix::HadrontherapyMatrix(G4int voxelX, G4int voxelY, G4int voxelZ, G4double mass):
     61    filename("Dose.out"),
     62    doseUnit(MeV/g)
    5163{ 
    52 // Number of the voxels of the phantom
    53 // For Y = Z = 1 the phantom is divided in slices (and not in voxels)
    54 // orthogonal to the beam axis
    55   numberVoxelX = voxelX;
    56   numberVoxelY = voxelY;
    57   numberVoxelZ = voxelZ;
    58   // Create the dose matrix
    59   matrix = new G4double[numberVoxelX*numberVoxelY*numberVoxelZ];
    60   if (matrix) G4cout << "Matrix: Memory space to store physical dose into " << 
    61                         numberVoxelX*numberVoxelY*numberVoxelZ <<
    62                         " voxels has been allocated " << G4endl;
    63   else G4Exception("Can't allocate memory to store physical dose!");
    64 }
    65 
     64        // Number of the voxels of the phantom
     65        // For Y = Z = 1 the phantom is divided in slices (and not in voxels)
     66        // orthogonal to the beam axis
     67        numberOfVoxelAlongX = voxelX;
     68        numberOfVoxelAlongY = voxelY;
     69        numberOfVoxelAlongZ = voxelZ;
     70        massOfVoxel = mass;
     71        // Create the dose matrix
     72        matrix = new G4double[numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ];
     73        if (matrix)
     74        {
     75                G4cout << "HadrontherapyMatrix: Memory space to store physical dose into " << 
     76                numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ <<
     77                " voxels has been allocated " << G4endl;
     78        }
     79        else G4Exception(" HadrontherapyMatrix::HadrontherapyMatrix. Can't allocate memory to store physical dose!");
     80                // Hit voxel (TrackID) marker
     81                // This array mark the status of voxel, if a hit occur, with the trackID of the particle
     82                // Must be initialized
     83        hitTrack = new G4int[numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ];
     84        ClearHitTrack();
     85}
     86
     87/////////////////////////////////////////////////////////////////////////////
    6688HadrontherapyMatrix::~HadrontherapyMatrix()
    6789{
    6890    delete[] matrix;
    69 }
    70 
    71 void HadrontherapyMatrix::flush()
    72 {
    73         if(matrix)
    74         for(int i=0; i<numberVoxelX*numberVoxelY*numberVoxelZ; i++)
    75         {
    76           matrix[i] = 0;
    77         }
    78 }
     91    delete[] hitTrack;
     92                // clear fluences/dose data
     93    Clear();
     94}
     95
     96/////////////////////////////////////////////////////////////////////////////
     97void HadrontherapyMatrix::Clear()
     98{
     99    for (size_t i=0; i<ionStore.size(); i++)
     100    {
     101        delete[] ionStore[i].dose;
     102        delete[] ionStore[i].fluence;
     103    }
     104    ionStore.clear();
     105}
     106
     107/////////////////////////////////////////////////////////////////////////////
     108// Initialise the elements of the matrix to zero
    79109void HadrontherapyMatrix::Initialize()
    80110{
    81 // Initialise the elements of the matrix to zero
    82   for(G4int i = 0; i < numberVoxelX; i++)
    83     {
    84       for(G4int j = 0; j < numberVoxelY; j++)
    85            {
    86               for(G4int k = 0; k < numberVoxelZ; k++)
    87 
    88               matrix[Index(i,j,k)] = 0.;
    89            }
    90     }
    91 }
     111    Clear();
     112    for(int i=0;i<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ;i++)
     113    {
     114                matrix[i] = 0;
     115    }
     116}
     117        /////////////////////////////////////////////////////////////////////////////
     118        /////////////////////////////////////////////////////////////////////////////
     119        // Print generated nuclides list
     120void HadrontherapyMatrix::PrintNuclides()
     121{
     122    for (size_t i=0; i<ionStore.size(); i++)
     123    {
     124                G4cout << ionStore[i].name << G4endl;
     125    }
     126}
     127        /////////////////////////////////////////////////////////////////////////////
     128        // Clear Hit voxel (TrackID) markers
     129void HadrontherapyMatrix::ClearHitTrack()
     130{
     131        for(G4int i=0; i<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ; i++) hitTrack[i] = 0;
     132}
     133        // Return Hit status
     134G4int* HadrontherapyMatrix::GetHitTrack(G4int i, G4int j, G4int k)
     135{
     136        return &(hitTrack[Index(i,j,k)]);
     137}
     138/////////////////////////////////////////////////////////////////////////////
     139// Dose methods...
     140// Fill DOSE/fluence matrix for particle/ion:
     141// If fluence parameter is true (default value is FALSE) then fluence at voxel (i, j, k) is increased.
     142// The energyDeposit parameter fill the dose matrix for voxel (i,j,k)
     143/////////////////////////////////////////////////////////////////////////////
     144
     145G4bool HadrontherapyMatrix::Fill(G4int trackID,
     146                                 G4ParticleDefinition* particleDef,
     147                                 G4int i, G4int j, G4int k,
     148                                 G4double energyDeposit,
     149                                 G4bool fluence)
     150{
     151    if (energyDeposit <=0. && !fluence || !secondaries) return false;
     152    // Get Particle Data Group particle ID
     153    G4int PDGencoding = particleDef -> GetPDGEncoding();
     154    PDGencoding -= PDGencoding%10;
     155       
     156    // Search for already allocated data...
     157    for (size_t l=0; l < ionStore.size(); l++)
     158    {
     159                if (ionStore[l].PDGencoding == PDGencoding )
     160                {   // Is it a primary or a secondary particle?
     161                        if ( trackID ==1 && ionStore[l].isPrimary || trackID !=1 && !ionStore[l].isPrimary)
     162                        {
     163                                if (energyDeposit > 0.) ionStore[l].dose[Index(i, j, k)] += energyDeposit/massOfVoxel;
     164                               
     165                                        // Fill a matrix per each ion with the fluence
     166                                if (fluence) ionStore[l].fluence[Index(i, j, k)]++;
     167                                return true;
     168                        }
     169                }
     170    }
     171       
     172    G4int Z = particleDef-> GetAtomicNumber();
     173    G4int A = particleDef-> GetAtomicMass();
     174
     175    G4String fullName = particleDef -> GetParticleName();
     176    G4String name = fullName.substr (0, fullName.find("[") ); // cut excitation energy 
     177  // Let's put a new particle in our store...
     178    ion newIon =
     179    {
     180                (trackID == 1) ? true:false,
     181                PDGencoding,
     182                name,
     183                name.length(),
     184                Z,
     185                A,
     186                new G4double[numberOfVoxelAlongX * numberOfVoxelAlongY * numberOfVoxelAlongZ],
     187                new unsigned int[numberOfVoxelAlongX * numberOfVoxelAlongY * numberOfVoxelAlongZ]
     188    };
     189                // Initialize data
     190    if (newIon.dose && newIon.fluence)
     191    {
     192                for(G4int m=0; m<numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ; m++)
     193                {
     194                        newIon.dose[m] = 0.;
     195                        newIon.fluence[m] = 0;
     196                }
     197                if (energyDeposit > 0.) newIon.dose[Index(i, j, k)] += energyDeposit/massOfVoxel;
     198                if (fluence) newIon.fluence[Index(i, j, k)]++;
     199               
     200                ionStore.push_back(newIon);
     201               
     202                // TODO Put some verbosity check
     203                /*
     204                 G4cout << "Memory space to store the DOSE/FLUENCE into " << 
     205                 numberOfVoxelAlongX*numberOfVoxelAlongY*numberOfVoxelAlongZ <<
     206                 " voxels has been allocated for the nuclide " << newIon.name <<
     207                 " (Z = " << Z << ", A = " << A << ")" << G4endl ;
     208                 */
     209                return true;
     210    }
     211    else // XXX Out of memory! XXX
     212    {
     213                return false;
     214    }
     215       
     216}
     217       
     218        /////////////////////////////////////////////////////////////////////////////
     219        /////////////////////////////////////////////////////////////////////////////
     220        // Methods to store data to filenames...
     221        ////////////////////////////////////////////////////////////////////////////
     222        ////////////////////////////////////////////////////////////////////////////
     223        //
     224        // General method to store matrix data to filename
     225void HadrontherapyMatrix::StoreMatrix(G4String file, void* data, size_t psize)
     226{
     227    if (data)
     228    {
     229                ofs.open(file, std::ios::out);
     230                if (ofs.is_open())
     231                {
     232                        for(G4int i = 0; i < numberOfVoxelAlongX; i++)
     233                                for(G4int j = 0; j < numberOfVoxelAlongY; j++)
     234                                        for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
     235                                        {
     236                                                G4int n = Index(i, j, k);
     237                                                        // Check for data type: u_int, G4double, XXX
     238                                                if (psize == sizeof(unsigned int))
     239                                                {
     240                                                        unsigned int* pdata = (unsigned int*)data;
     241                                                        if (pdata[n]) ofs << i << '\t' << j << '\t' <<
     242                                                                k << '\t' << pdata[n] << G4endl;
     243                                                }
     244                                                else if (psize == sizeof(G4double))
     245                                                {
     246                                                        G4double* pdata = (G4double*)data;
     247                                                        if (pdata[n]) ofs << i << '\t' << j << '\t' <<
     248                                                                k << '\t' << pdata[n] << G4endl;
     249                                                }
     250                                        }
     251                        ofs.close();
     252                }
     253    }
     254}
     255
     256        // Store fluence per single ion in multiple files
     257void HadrontherapyMatrix::StoreFluenceData()
     258{
     259    for (size_t i=0; i < ionStore.size(); i++){
     260                StoreMatrix(ionStore[i].name + "_Fluence.out", ionStore[i].fluence, sizeof(unsigned int));
     261    }
     262}
     263        // Store dose per single ion in multiple files
     264void HadrontherapyMatrix::StoreDoseData()
     265{
     266       
     267    for (size_t i=0; i < ionStore.size(); i++){
     268                StoreMatrix(ionStore[i].name + "_Dose.out", ionStore[i].dose, sizeof(G4double));
     269    }
     270}
     271        /////////////////////////////////////////////////////////////////////////
     272        // Store dose for all ions into a single file and into ntuples.
     273        // Please note that this function is called via messenger commands
     274        // defined in the HadrontherapyAnalysisFileMessenger.cc class file
     275void HadrontherapyMatrix::StoreDoseFluenceAscii()
     276{
     277    // Sort like periodic table
     278    std::sort(ionStore.begin(), ionStore.end());
     279#define width 15L
     280    ofs.open(filename, std::ios::out);
     281    if (ofs.is_open())
     282    {
     283
     284        //
     285        // Write the voxels index and the list of particles/ions
     286        ofs << std::setprecision(6) << std::left <<
     287            "i\tj\tk\t";
     288/*         
     289            G4RunManager *runManager = G4RunManager::GetRunManager();
     290            HadrontherapyPrimaryGeneratorAction *pPGA = (HadrontherapyPrimaryGeneratorAction*)runManager -> GetUserPrimaryGeneratorAction();
     291            G4String name = pPGA -> GetParticleGun() -> GetParticleDefinition() -> GetParticleName();
     292*/         
     293        // Total dose
     294        ofs << std::setw(width) << "Dose";
     295
     296        for (size_t l=0; l < ionStore.size(); l++)
     297        {
     298            G4String a = (ionStore[l].isPrimary) ? "_1":""; // is it a primary?
     299            ofs << std::setw(width) << ionStore[l].name + a <<
     300                std::setw(width) << ionStore[l].name  + a;
     301        }
     302        ofs << G4endl;
     303
     304/*
     305        // PDGencoding
     306        ofs << std::setprecision(6) << std::left <<
     307        "0\t0\t0\t";
     308
     309        // Total dose
     310        ofs << std::setw(width) << '0';
     311        for (size_t l=0; l < ionStore.size(); l++)
     312        {
     313        ofs << std::setw(width) << ionStore[l].PDGencoding  <<
     314        std::setw(width) << ionStore[l].PDGencoding;
     315        }
     316        ofs << G4endl;
     317*/
     318        // Write data
     319        for(G4int i = 0; i < numberOfVoxelAlongX; i++)
     320            for(G4int j = 0; j < numberOfVoxelAlongY; j++)
     321                for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
     322                {
     323                    G4int n = Index(i, j, k);
     324                    // Write only not identically null data lines
     325                    if (matrix[n])
     326                    {
     327                        ofs << G4endl;
     328                        ofs << i << '\t' << j << '\t' << k << '\t';
     329                        // Total dose
     330                        ofs << std::setw(width) << matrix[n]/(doseUnit);
     331                        {
     332                            for (size_t l=0; l < ionStore.size(); l++)
     333                            {
     334                                // Fill ASCII file rows
     335                                ofs << std::setw(width) << ionStore[l].dose[n]/(doseUnit) <<
     336                                    std::setw(width) << ionStore[l].fluence[n];
     337                            }
     338                        }
     339                    }
     340                }
     341        ofs.close();
     342    }
     343}
     344/////////////////////////////////////////////////////////////////////////////
     345
     346#ifdef G4ANALYSIS_USE_ROOT
     347void HadrontherapyMatrix::StoreDoseFluenceRoot()
     348{
     349    HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance();
     350    for(G4int i = 0; i < numberOfVoxelAlongX; i++)
     351        for(G4int j = 0; j < numberOfVoxelAlongY; j++)
     352            for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
     353            {
     354                G4int n = Index(i, j, k);
     355                for (size_t l=0; l < ionStore.size(); l++)
     356
     357                {
     358                    // Do the same work for .root file: fill dose/fluence ntuple 
     359                    analysis -> FillVoxelFragmentTuple( i, j, k,
     360                                                        ionStore[l].A,
     361                                                        ionStore[l].Z,
     362                                                        ionStore[l].dose[n]/(doseUnit),
     363                                                        ionStore[l].fluence[n] );
     364
     365
     366                }
     367            }
     368}
     369#endif
    92370
    93371void HadrontherapyMatrix::Fill(G4int i, G4int j, G4int k,
    94372                               G4double energyDeposit)
    95373{
    96   if (matrix)
    97     matrix[Index(i,j,k)] += energyDeposit;
    98  
    99   // Store the energy deposit in the matrix element corresponding
    100   // to the phantom voxel  i, j, k
    101 }
    102 
     374    if (matrix)
     375                matrix[Index(i,j,k)] += energyDeposit;
     376       
     377    // Store the energy deposit in the matrix element corresponding
     378    // to the phantom voxel 
     379}
    103380void HadrontherapyMatrix::TotalEnergyDeposit()
    104381{
     382  // Convert energy deposited to dose.
    105383  // Store the information of the matrix in a ntuple and in
    106384  // a 1D Histogram
    107385
    108   G4int k;
    109   G4int j;
    110   G4int i;
    111  
    112   if (matrix)
     386#ifdef G4ANALYSIS_USE_ROOT
     387    HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance();
     388#endif
     389    if (matrix)
    113390    { 
    114       std::ofstream ofs;       
    115       ofs.open("DoseDistribution.out");
    116      
    117       for(G4int l = 0; l < numberVoxelZ; l++)
    118         {
    119           k = l;
    120          
    121           for(G4int m = 0; m < numberVoxelY; m++)
    122             {
    123               j = m * numberVoxelZ + k;
    124                
    125                 for(G4int n = 0; n <  numberVoxelX; n++)
    126                   {
    127                     i =  n* numberVoxelZ * numberVoxelY + j;
    128                     if(matrix[i] != 0)
    129                       {
    130                         ofs << n << '\t' << m << '\t' <<
    131                           k << '\t' << matrix[i] << G4endl;
    132                        
    133 #ifdef ANALYSIS_USE
    134                         HadrontherapyAnalysisManager* analysis =
    135                         HadrontherapyAnalysisManager::getInstance();
    136                         analysis -> FillEnergyDeposit(n, m, k, matrix[i]);
    137                         analysis -> BraggPeak(n, matrix[i]);
     391        for(G4int i = 0; i < numberOfVoxelAlongX; i++)
     392            for(G4int j = 0; j < numberOfVoxelAlongY; j++)
     393                for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
     394                {
     395                    G4int n = Index(i,j,k);
     396                    matrix[n] = matrix[n] / massOfVoxel;
     397#ifdef G4ANALYSIS_USE_ROOT
     398                    analysis -> FillEnergyDeposit(i, j, k, matrix[n]/doseUnit);
     399                    analysis -> BraggPeak(i, matrix[n]/doseUnit);
    138400#endif
    139                       }
    140 }       
    141               }
    142           }
    143     }
    144 }
     401                }
     402    }
     403}
     404

  • trunk/examples/advanced/hadrontherapy/src/HadrontherapyPhysicsList.cc

    r1230 r1313  
    7676//                        and local physic for ion-ion enelastic processes)
    7777
     78#include "G4RunManager.hh"
     79#include "G4Region.hh"
     80#include "G4RegionStore.hh"
    7881#include "HadrontherapyPhysicsList.hh"
    7982#include "HadrontherapyPhysicsListMessenger.hh"
    … …  
    8689#include "LocalINCLIonIonInelasticPhysic.hh"         // Physic dedicated to the ion-ion inelastic processes using INCL/ABLA
    8790
    88 // Physic lists (contained inside the Geant4 distribution)
     91// Physic lists (contained inside the Geant4 source code, in the 'physicslists folder')
    8992#include "G4EmStandardPhysics_option3.hh"
    9093#include "G4EmLivermorePhysics.hh"
    … …  
    173176{
    174177  // transportation
    175   //
    176178  AddTransportation();
    177179
    178180  // electromagnetic physics list
    179   //
    180181  emPhysicsList->ConstructProcess();
    181182  em_config.AddModels();
    182183
    183184  // decay physics list
    184   //
    185185  decPhysicsList->ConstructProcess();
    186186
    … …  
    207207  //   ELECTROMAGNETIC MODELS
    208208  /////////////////////////////////////////////////////////////////////////////
    209 
    210   if (name == "standard_opt3") {
     209    if (name == "standard_opt3") {
    211210    emName = name;
    212211    delete emPhysicsList;
    213212    emPhysicsList = new G4EmStandardPhysics_option3();
     213    G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
    214214    G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmStandardPhysics_option3" << G4endl;
    215215
    … …  
    218218    delete emPhysicsList;
    219219    emPhysicsList = new G4EmLivermorePhysics();
     220    G4RunManager::GetRunManager()-> PhysicsHasBeenModified();
    220221    G4cout << "THE FOLLOWING ELECTROMAGNETIC PHYSICS LIST HAS BEEN ACTIVATED: G4EmLivermorePhysics" << G4endl;
    221222
    … …  
    308309  SetCutValue(cutForPositron, "e+");
    309310
     311  // Set cuts for detector
     312  SetDetectorCut(defaultCutValue);
    310313  if (verboseLevel>0) DumpCutValuesTable();
    311314}
    … …  
    331334  SetParticleCuts(cutForPositron, G4Positron::Positron());
    332335}
     336
     337void HadrontherapyPhysicsList::SetDetectorCut(G4double cut)
     338{
     339
     340  G4String regionName = "DetectorLog";
     341  G4Region* region = G4RegionStore::GetInstance()->GetRegion(regionName);
     342
     343  G4ProductionCuts* cuts = new G4ProductionCuts ;
     344  cuts -> SetProductionCut(cut,G4ProductionCuts::GetIndex("gamma"));
     345  cuts -> SetProductionCut(cut,G4ProductionCuts::GetIndex("e-"));
     346  cuts -> SetProductionCut(cut,G4ProductionCuts::GetIndex("e+"));
     347  region -> SetProductionCuts(cuts);
     348}
     349

  • trunk/examples/advanced/hadrontherapy/src/HadrontherapyPhysicsListMessenger.cc

    r1230 r1313  
    6969  allCutCmd->AvailableForStates(G4State_PreInit,G4State_Idle); 
    7070
     71  allDetectorCmd = new G4UIcmdWithADoubleAndUnit("/physic/setDetectorCuts",this); 
     72  allDetectorCmd->SetGuidance("Set cut for all. into Detector");
     73  allDetectorCmd->SetParameterName("cut",false);
     74  allDetectorCmd->SetUnitCategory("Length");
     75  allDetectorCmd->SetRange("cut>0.0");
     76  allDetectorCmd->AvailableForStates(G4State_PreInit,G4State_Idle); 
     77
    7178  pListCmd = new G4UIcmdWithAString("/physic/addPhysics",this); 
    7279  pListCmd->SetGuidance("Add physics list.");
    7380  pListCmd->SetParameterName("PList",false);
    74   pListCmd->AvailableForStates(G4State_PreInit); 
     81  pListCmd->AvailableForStates(G4State_PreInit, G4State_Idle); 
    7582
    7683  packageListCmd = new G4UIcmdWithAString("/physic/addPackage",this);
    7784  packageListCmd->SetGuidance("Add physics package.");
    7885  packageListCmd->SetParameterName("package",false);
    79   packageListCmd->AvailableForStates(G4State_PreInit);
     86  packageListCmd->AvailableForStates(G4State_PreInit, G4State_Idle);
    8087}
    8188
    … …  
    8794  delete protoCutCmd;
    8895  delete allCutCmd;
     96  delete allDetectorCmd;
    8997  delete pListCmd;
    9098  delete physDir;   
    … …  
    99107   { pPhysicsList->SetCutForGamma(gammaCutCmd->GetNewDoubleValue(newValue));}
    100108     
    101   if( command == electCutCmd )
     109  else if( command == electCutCmd )
    102110   { pPhysicsList->SetCutForElectron(electCutCmd->GetNewDoubleValue(newValue));}
    103111     
    104   if( command == protoCutCmd )
     112  else if( command == protoCutCmd )
    105113   { pPhysicsList->SetCutForPositron(protoCutCmd->GetNewDoubleValue(newValue));}
    106114
    107   if( command == allCutCmd )
     115  else if( command == allCutCmd )
    108116    {
    109117      G4double cut = allCutCmd->GetNewDoubleValue(newValue);
    … …  
    112120      pPhysicsList->SetCutForPositron(cut);
    113121    }
    114 
    115   if( command == pListCmd )
     122  else if( command == allDetectorCmd)
     123  {
     124      G4double cut = allDetectorCmd -> GetNewDoubleValue(newValue);
     125      pPhysicsList -> SetDetectorCut(cut);
     126  }
     127  else if( command == pListCmd )
    116128   { pPhysicsList->AddPhysicsList(newValue);}
    117129
    118130
    119   if( command == packageListCmd )
     131  else if( command == packageListCmd )
    120132   { pPhysicsList->AddPackage(newValue);}
    121133

  • trunk/examples/advanced/hadrontherapy/src/HadrontherapyPrimaryGeneratorAction.cc

    r1230 r1313  
    8686  sigmaEnergy = defaultsigmaEnergy;
    8787 
    88 #ifdef ANALYSIS_USE
     88#ifdef G4ANALYSIS_USE_ROOT
    8989  // Write these values into the analysis if needed. Have to be written separately on change.
    90   HadrontherapyAnalysisManager::getInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
     90  HadrontherapyAnalysisManager::GetInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
    9191#endif
    9292
    … …  
    120120void HadrontherapyPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
    121121{
    122 #ifdef ANALYSIS_USE
     122#ifdef G4ANALYSIS_USE_ROOT
    123123  // Increment the event counter
    124   HadrontherapyAnalysisManager::getInstance()->startNewEvent();
     124  HadrontherapyAnalysisManager::GetInstance()->startNewEvent();
    125125#endif
    126126
    … …  
    178178{
    179179        meanKineticEnergy = val;
    180 #ifdef ANALYSIS_USE
     180#ifdef G4ANALYSIS_USE_ROOT
    181181  // Update the beam-data in the analysis manager
    182   HadrontherapyAnalysisManager::getInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
     182  HadrontherapyAnalysisManager::GetInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
    183183#endif
    184184
    … …  
    188188{
    189189        sigmaEnergy = val;
    190 #ifdef ANALYSIS_USE
     190#ifdef G4ANALYSIS_USE_ROOT
    191191  // Update the sigmaenergy in the metadata.
    192   HadrontherapyAnalysisManager::getInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
     192  HadrontherapyAnalysisManager::GetInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
    193193#endif
    194194}
    … …  
    217217G4double HadrontherapyPrimaryGeneratorAction::GetmeanKineticEnergy(void)
    218218{ return meanKineticEnergy;}
     219

  • trunk/examples/advanced/hadrontherapy/src/HadrontherapyRunAction.cc

    r1230 r1313  
    2727// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
    2828
     29
    2930#include "HadrontherapyRunAction.hh"
    3031#include "HadrontherapyEventAction.hh"
    … …  
    3839#include "HadrontherapyRunAction.hh"
    3940
     41
     42#include "HadrontherapyAnalysisManager.hh"
     43#include "HadrontherapyMatrix.hh"
     44
    4045HadrontherapyRunAction::HadrontherapyRunAction()
    4146{
    … …  
    4853void HadrontherapyRunAction::BeginOfRunAction(const G4Run* aRun)
    4954{       
    50    G4RunManager::GetRunManager()->SetRandomNumberStore(true);
    51    G4cout << "Run " << aRun -> GetRunID() << " starts ..." << G4endl;
     55    G4RunManager::GetRunManager()-> SetRandomNumberStore(true);
     56    G4cout << "Run " << aRun -> GetRunID() << " starts ..." << G4endl;
    5257
    53    electromagnetic = 0;
    54    hadronic = 0;
     58    // Warning! any beamOn will reset all data (dose, fluence, histograms, etc)!
     59    //
     60
     61    // Initialize matrix with energy of primaries clearing data inside
     62    if (HadrontherapyMatrix::GetInstance())
     63    {
     64        HadrontherapyMatrix::GetInstance() -> Initialize();
     65    }
     66
     67#ifdef G4ANALYSIS_USE_ROOT
     68    HadrontherapyAnalysisManager::GetInstance() -> flush();     // Finalize the root file
     69    // Initialize root analysis ----> book
     70    HadrontherapyAnalysisManager::GetInstance() -> book();
     71#endif
     72    electromagnetic = 0;
     73    hadronic = 0;
    5574}
    5675
    5776void HadrontherapyRunAction::EndOfRunAction(const G4Run*)
    5877{
    59   //   G4cout << " Summary of Run " << aRun -> GetRunID() <<" :"<< G4endl;
    60   //G4cout << "Number of electromagnetic processes of primary particles in the phantom:"
    61   //     << electromagnetic << G4endl;
    62   //G4cout << "Number of hadronic processes of primary particles in the phantom:"
    63   //     << hadronic << G4endl;
     78    // Store dose & fluence data to ASCII & ROOT files
     79    if ( HadrontherapyMatrix::GetInstance() )
     80    {
     81        HadrontherapyMatrix::GetInstance() -> TotalEnergyDeposit();
     82        HadrontherapyMatrix::GetInstance() -> StoreDoseFluenceAscii();
     83
     84#ifdef G4ANALYSIS_USE_ROOT
     85        if (HadrontherapyAnalysisManager::GetInstance())
     86        {
     87            HadrontherapyMatrix::GetInstance() -> StoreDoseFluenceRoot();
     88        }
     89#endif
     90    }
     91
     92    //G4cout << " Summary of Run " << aRun -> GetRunID() <<" :"<< G4endl;
     93    //G4cout << "Number of electromagnetic processes of primary particles in the phantom:"
     94    //     << electromagnetic << G4endl;
     95    //G4cout << "Number of hadronic processes of primary particles in the phantom:"
     96    //     << hadronic << G4endl;
    6497}
    6598void HadrontherapyRunAction::AddEMProcess()
    … …  
    72105}
    73106
    74 
    75 

  • trunk/examples/advanced/hadrontherapy/src/HadrontherapySteppingAction.cc

    r1230 r1313  
    3939#include "G4ParticleDefinition.hh"
    4040#include "G4ParticleTypes.hh"
     41#include "G4UserEventAction.hh"
    4142
    4243#include "HadrontherapyAnalysisManager.hh"
    … …  
    4748HadrontherapySteppingAction::HadrontherapySteppingAction( HadrontherapyRunAction *run)
    4849{
    49   runAction = run;
     50    runAction = run;
    5051}
    5152
    … …  
    5859void HadrontherapySteppingAction::UserSteppingAction(const G4Step* aStep)
    5960{
     61    /*
     62    // USEFULL METHODS TO RETRIEVE INFORMATION DURING THE STEPS
     63    if( (aStep->GetTrack()->GetVolume()->GetName() == "DetectorPhys")
     64    && aStep->GetTrack()->GetDefinition()->GetParticleName() == "proton")
     65    //G4int evtNb = G4RunManager::GetRunManager()->GetCurrentEvent() -> GetEventID();
     66    {
     67    G4cout << "ENERGIA:   " << aStep->GetTrack()->GetKineticEnergy()
     68    << "   VOLUME  " << aStep->GetTrack()->GetVolume()->GetName()
     69    << "   MATERIALE    " <<  aStep -> GetTrack() -> GetMaterial() -> GetName()
     70    << "   EVENTO     " << G4RunManager::GetRunManager()->GetCurrentEvent() -> GetEventID()
     71    << "   POS     " << aStep->GetTrack()->GetPosition().x()
     72    << G4endl;
     73    }
     74    */
     75
    6076    if( aStep->GetTrack()->GetVolume()->GetName() == "NewDetectorPhys"){
    61 #ifdef ANALYSIS_USE
    62       G4ParticleDefinition *def = aStep->GetTrack()->GetDefinition();
    63       G4double secondaryParticleKineticEnergy =  aStep->GetTrack()->GetKineticEnergy();     
    64       G4String particleType = def->GetParticleType(); // particle type = nucleus for d, t, He3, alpha, and heavier nuclei
    65       G4String particleName = def->GetParticleName(); // e.g. for alpha: the name = "alpha" and type = "nucleus"
    66       if(particleType == "nucleus") {
    67         G4int A = def->GetBaryonNumber();
    68         G4double Z = def->GetPDGCharge();
    69         G4double posX = aStep->GetTrack()->GetPosition().x() / cm;
    70         G4double posY = aStep->GetTrack()->GetPosition().y() / cm;
    71         G4double posZ = aStep->GetTrack()->GetPosition().z() / cm;
    72         G4double energy = secondaryParticleKineticEnergy / A / MeV;
    73 
    74         HadrontherapyAnalysisManager* analysisMgr =  HadrontherapyAnalysisManager::getInstance();   
    75         analysisMgr->fillFragmentTuple(A, Z, energy, posX, posY, posZ);
    76       } else if(particleName == "proton") {   // proton (hydrogen-1) is a special case
    77         G4double posX = aStep->GetTrack()->GetPosition().x() / cm ;
    78         G4double posY = aStep->GetTrack()->GetPosition().y() / cm ;
    79         G4double posZ = aStep->GetTrack()->GetPosition().z() / cm ;
    80         G4double energy = secondaryParticleKineticEnergy * MeV;    // Hydrogen-1: A = 1, Z = 1
    81         HadrontherapyAnalysisManager::getInstance()->fillFragmentTuple(1, 1.0, energy, posX, posY, posZ);
    82       }
    83 
    84       G4String secondaryParticleName =  def -> GetParticleName(); 
    85       //G4cout <<"Particle: " << secondaryParticleName << G4endl;
    86       //G4cout <<"Energy: " << secondaryParticleKineticEnergy << G4endl;
    87         HadrontherapyAnalysisManager* analysis =  HadrontherapyAnalysisManager::getInstance();   
     77#ifdef G4ANALYSIS_USE_ROOT
     78        G4ParticleDefinition *def = aStep->GetTrack()->GetDefinition();
     79        G4double secondaryParticleKineticEnergy =  aStep->GetTrack()->GetKineticEnergy();     
     80        G4String particleType = def->GetParticleType(); // particle type = nucleus for d, t, He3, alpha, and heavier nuclei
     81        G4String particleName = def->GetParticleName(); // e.g. for alpha: the name = "alpha" and type = "nucleus"
     82        if(particleType == "nucleus") {
     83            G4int A = def->GetBaryonNumber();
     84            G4double Z = def->GetPDGCharge();
     85            G4double posX = aStep->GetTrack()->GetPosition().x() / cm;
     86            G4double posY = aStep->GetTrack()->GetPosition().y() / cm;
     87            G4double posZ = aStep->GetTrack()->GetPosition().z() / cm;
     88            G4double energy = secondaryParticleKineticEnergy / A / MeV;
     89
     90            HadrontherapyAnalysisManager* analysisMgr =  HadrontherapyAnalysisManager::GetInstance();   
     91            analysisMgr->FillFragmentTuple(A, Z, energy, posX, posY, posZ);
     92        } else if(particleName == "proton") {   // proton (hydrogen-1) is a special case
     93            G4double posX = aStep->GetTrack()->GetPosition().x() / cm ;
     94            G4double posY = aStep->GetTrack()->GetPosition().y() / cm ;
     95            G4double posZ = aStep->GetTrack()->GetPosition().z() / cm ;
     96            G4double energy = secondaryParticleKineticEnergy * MeV;    // Hydrogen-1: A = 1, Z = 1
     97            HadrontherapyAnalysisManager::GetInstance()->FillFragmentTuple(1, 1.0, energy, posX, posY, posZ);
     98        }
     99
     100        G4String secondaryParticleName =  def -> GetParticleName(); 
     101        //G4cout <<"Particle: " << secondaryParticleName << G4endl;
     102        //G4cout <<"Energy: " << secondaryParticleKineticEnergy << G4endl;
     103        HadrontherapyAnalysisManager* analysis =  HadrontherapyAnalysisManager::GetInstance();   
    88104        //There is a bunch of stuff recorded with the energy 0, something should perhaps be done about this.
    89105        if(secondaryParticleName == "proton") {
    90           analysis->hydrogenEnergy(secondaryParticleKineticEnergy / MeV);
     106            analysis->hydrogenEnergy(secondaryParticleKineticEnergy / MeV);
    91107        }
    92108        if(secondaryParticleName == "deuteron") {
    93           analysis->hydrogenEnergy((secondaryParticleKineticEnergy/2) / MeV);
     109            analysis->hydrogenEnergy((secondaryParticleKineticEnergy/2) / MeV);
    94110        }
    95111        if(secondaryParticleName == "triton") {
    96           analysis->hydrogenEnergy((secondaryParticleKineticEnergy/3) / MeV);
     112            analysis->hydrogenEnergy((secondaryParticleKineticEnergy/3) / MeV);
    97113        }
    98114        if(secondaryParticleName == "alpha") {
    99           analysis->heliumEnergy((secondaryParticleKineticEnergy/4) / MeV);
     115            analysis->heliumEnergy((secondaryParticleKineticEnergy/4) / MeV);
    100116        }
    101117        if(secondaryParticleName == "He3"){
    102           analysis->heliumEnergy((secondaryParticleKineticEnergy/3) / MeV);             
     118            analysis->heliumEnergy((secondaryParticleKineticEnergy/3) / MeV);           
    103119        }
    104120#endif
    … …  
    107123    }
    108124
    109   // Electromagnetic and hadronic processes of primary particles in the phantom
    110   //setting phantomPhys correctly will break something here fixme
    111   if ((aStep -> GetTrack() -> GetTrackID() == 1) &&
    112     (aStep -> GetTrack() -> GetVolume() -> GetName() == "PhantomPhys") &&
    113     (aStep -> GetPostStepPoint() -> GetProcessDefinedStep() != NULL))
     125    // Electromagnetic and hadronic processes of primary particles in the phantom
     126    //setting phantomPhys correctly will break something here fixme
     127    if ((aStep -> GetTrack() -> GetTrackID() == 1) &&
     128            (aStep -> GetTrack() -> GetVolume() -> GetName() == "PhantomPhys") &&
     129            (aStep -> GetPostStepPoint() -> GetProcessDefinedStep() != NULL))
     130    {
     131        G4String process = aStep -> GetPostStepPoint() ->
     132            GetProcessDefinedStep() -> GetProcessName();
     133
     134        if ((process == "Transportation") || (process == "StepLimiter")) {;}
     135        else
     136        {
     137            if ((process == "msc") || (process == "hLowEIoni") || (process == "hIoni"))
     138            {
     139                runAction -> AddEMProcess();
     140            }
     141            else 
    114142            {
    115               G4String process = aStep -> GetPostStepPoint() ->
    116                 GetProcessDefinedStep() -> GetProcessName();
    117  
    118               if ((process == "Transportation") || (process == "StepLimiter")) {;}
    119               else {
    120                 if ((process == "msc") || (process == "hLowEIoni") || (process == "hIoni"))
    121                   {
    122                     runAction -> AddEMProcess();
    123                   }
    124                 else 
    125                  {
    126                    runAction -> AddHadronicProcess();
    127 
    128                    if ( (process != "LElastic") && (process != "ProtonInelastic") && (process != "hElastic") )
    129                      G4cout << "Warning! Unknown proton process: "<< process << G4endl;
    130                  }
    131               }         
     143                runAction -> AddHadronicProcess();
     144
     145                if ( (process != "LElastic") && (process != "ProtonInelastic") && (process != "hElastic") )
     146                    G4cout << "Warning! Unknown proton process: "<< process << G4endl;
    132147            }
    133 
    134  // Retrieve information about the secondaries originated in the phantom
    135 
    136  G4SteppingManager*  steppingManager = fpSteppingManager;
    137  
    138   // check if it is alive
    139   //if(theTrack-> GetTrackStatus() == fAlive) { return; }
    140 
    141   // Retrieve the secondary particles
    142   G4TrackVector* fSecondary = steppingManager -> GetfSecondary();
    143      
    144   for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++)
     148        }         
     149    }
     150
     151    // Retrieve information about the secondaries originated in the phantom
     152
     153    G4SteppingManager*  steppingManager = fpSteppingManager;
     154
     155    // check if it is alive
     156    //if(theTrack-> GetTrackStatus() == fAlive) { return; }
     157
     158    // Retrieve the secondary particles
     159    G4TrackVector* fSecondary = steppingManager -> GetfSecondary();
     160
     161    for(size_t lp1=0;lp1<(*fSecondary).size(); lp1++)
    145162    {
    146       G4String volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName();
    147  
    148       if (volumeName == "phantomPhys")
     163        G4String volumeName = (*fSecondary)[lp1] -> GetVolume() -> GetName();
     164
     165        if (volumeName == "phantomPhys")
    149166        {
    150 #ifdef ANALYSIS_USE   
    151           G4String secondaryParticleName =  (*fSecondary)[lp1]->GetDefinition() -> GetParticleName(); 
    152           G4double secondaryParticleKineticEnergy =  (*fSecondary)[lp1] -> GetKineticEnergy();     
    153 
    154           HadrontherapyAnalysisManager* analysis =  HadrontherapyAnalysisManager::getInstance();   
    155        
    156           if (secondaryParticleName == "e-")
    157             analysis -> electronEnergyDistribution(secondaryParticleKineticEnergy/MeV);
    158                    
    159           if (secondaryParticleName == "gamma")
    160             analysis -> gammaEnergyDistribution(secondaryParticleKineticEnergy/MeV);
    161 
    162           if (secondaryParticleName == "deuteron")
    163             analysis -> deuteronEnergyDistribution(secondaryParticleKineticEnergy/MeV);
    164        
    165           if (secondaryParticleName == "triton")
    166             analysis -> tritonEnergyDistribution(secondaryParticleKineticEnergy/MeV);
    167            
    168           if (secondaryParticleName == "alpha")
    169             analysis -> alphaEnergyDistribution(secondaryParticleKineticEnergy/MeV);
    170        
    171           G4double z = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetPDGCharge();
    172           if (z > 0.)
    173             {     
    174               G4int a = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetBaryonNumber();
    175               G4int electronOccupancy = (*fSecondary)[lp1] ->  GetDynamicParticle() -> GetTotalOccupancy();
    176              
    177               // If a generic ion is originated in the detector, its baryonic number, PDG charge,
    178               // total number of electrons in the orbitals are stored in a ntuple
    179               analysis -> genericIonInformation(a, z, electronOccupancy, secondaryParticleKineticEnergy/MeV);
     167#ifdef G4ANALYSIS_USE_ROOT   
     168            G4String secondaryParticleName =  (*fSecondary)[lp1]->GetDefinition() -> GetParticleName(); 
     169            G4double secondaryParticleKineticEnergy =  (*fSecondary)[lp1] -> GetKineticEnergy();     
     170
     171            HadrontherapyAnalysisManager* analysis =  HadrontherapyAnalysisManager::GetInstance();   
     172
     173            if (secondaryParticleName == "e-")
     174                analysis -> electronEnergyDistribution(secondaryParticleKineticEnergy/MeV);
     175
     176            if (secondaryParticleName == "gamma")
     177                analysis -> gammaEnergyDistribution(secondaryParticleKineticEnergy/MeV);
     178
     179            if (secondaryParticleName == "deuteron")
     180                analysis -> deuteronEnergyDistribution(secondaryParticleKineticEnergy/MeV);
     181
     182            if (secondaryParticleName == "triton")
     183                analysis -> tritonEnergyDistribution(secondaryParticleKineticEnergy/MeV);
     184
     185            if (secondaryParticleName == "alpha")
     186                analysis -> alphaEnergyDistribution(secondaryParticleKineticEnergy/MeV);
     187
     188            G4double z = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetPDGCharge();
     189            if (z > 0.)
     190            {     
     191                G4int a = (*fSecondary)[lp1]-> GetDynamicParticle() -> GetDefinition() -> GetBaryonNumber();
     192                G4int electronOccupancy = (*fSecondary)[lp1] ->  GetDynamicParticle() -> GetTotalOccupancy();
     193
     194                // If a generic ion is originated in the detector, its baryonic number, PDG charge,
     195                // total number of electrons in the orbitals are stored in a ntuple
     196                analysis -> genericIonInformation(a, z, electronOccupancy, secondaryParticleKineticEnergy/MeV);
    180197            }
    181198#endif
Note: See TracChangeset for help on using the changeset viewer.