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

geant4.9.4 beta rc0

File:
1 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 
Note: See TracChangeset for help on using the changeset viewer.