source: trunk/examples/advanced/hadrontherapy/src/HadrontherapyAnalysisManager.cc@ 1344

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// $Id: HadrontherapyAnalisysManager.cc;
// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy

#include "HadrontherapyAnalysisManager.hh"
#include "HadrontherapyMatrix.hh"
#include "HadrontherapyAnalysisFileMessenger.hh"
#include <time.h>

HadrontherapyAnalysisManager* HadrontherapyAnalysisManager::instance = 0;

HadrontherapyAnalysisManager::HadrontherapyAnalysisManager(): 
#ifdef G4ANALYSIS_USE_ROOT 
analysisFileName("DoseDistribution.root"),theTFile(0), histo1(0), histo2(0), histo3(0),
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),
kinFragNtuple(0),
kineticEnergyPrimaryNtuple(0),
doseFragNtuple(0),
fluenceFragNtuple(0),
letFragNtuple(0),
theROOTNtuple(0),
theROOTIonTuple(0),
fragmentNtuple(0),
metaData(0),
eventCounter(0)
{
    fMess = new HadrontherapyAnalysisFileMessenger(this);
}
#endif
/////////////////////////////////////////////////////////////////////////////

HadrontherapyAnalysisManager::~HadrontherapyAnalysisManager()
{
    delete(fMess); 
#ifdef G4ANALYSIS_USE_ROOT
    Clear();    
#endif
}

HadrontherapyAnalysisManager* HadrontherapyAnalysisManager::GetInstance()
{
        if (instance == 0) instance = new HadrontherapyAnalysisManager;
        return instance;
}
#ifdef G4ANALYSIS_USE_ROOT
void HadrontherapyAnalysisManager::Clear()
{
    // XXX delete ALL variables!!
    // but this seems to be automatically done by 
    // theTFile->Close() command!
    // so we get a *double free or corruption* 

    delete metaData;
    metaData = 0;

    delete fragmentNtuple;
    fragmentNtuple = 0;

    delete theROOTIonTuple;
    theROOTIonTuple = 0;

    delete theROOTNtuple;
    theROOTNtuple = 0;

    delete histo16;
    histo16 = 0;

    delete histo15;
    histo15 = 0;

    delete histo14;
    histo14 = 0;

    delete histo13;
    histo13 = 0;

    delete histo12;
    histo12 = 0;

    delete histo11;
    histo11 = 0;

    delete histo10;
    histo10 = 0;

    delete histo9;
    histo9 = 0;

    delete histo8;
    histo8 = 0;

    delete histo7;
    histo7 = 0;

    delete histo6;
    histo6 = 0;

    delete histo5;
    histo5 = 0;

    delete histo4;
    histo4 = 0;

    delete histo3;
    histo3 = 0;

    delete histo2;
    histo2 = 0;

    delete histo1;
    histo1 = 0;
}
/////////////////////////////////////////////////////////////////////////////

void HadrontherapyAnalysisManager::SetAnalysisFileName(G4String aFileName)
{
        this->analysisFileName = aFileName;
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::book()
{
        delete theTFile; // this is similar to theTFile->Close() => delete all associated variables created via new, moreover it delete itself.  

        theTFile = new TFile(analysisFileName, "RECREATE");

        // Create the histograms with the energy deposit along the X axis
        histo1 = createHistogram1D("braggPeak","slice, energy", 400, 0., 80); //<different waterthicknesses are accoutned for in ROOT-analysis stage
        histo2 = createHistogram1D("h20","Secondary protons - slice, energy", 400, 0., 400.);
        histo3 = createHistogram1D("h30","Secondary neutrons - slice, energy", 400, 0., 400.);
        histo4 = createHistogram1D("h40","Secondary alpha - slice, energy", 400, 0., 400.);
        histo5 = createHistogram1D("h50","Secondary gamma - slice, energy", 400, 0., 400.);
        histo6 = createHistogram1D("h60","Secondary electron - slice, energy", 400, 0., 400.);
        histo7 = createHistogram1D("h70","Secondary triton - slice, energy", 400, 0., 400.);
        histo8 = createHistogram1D("h80","Secondary deuteron - slice, energy", 400, 0., 400.);
        histo9 = createHistogram1D("h90","Secondary pion - slice, energy", 400, 0., 400.);
        histo10 = createHistogram1D("h100","Energy distribution of secondary electrons", 70, 0., 70.);
        histo11 = createHistogram1D("h110","Energy distribution of secondary photons", 70, 0., 70.);
        histo12 = createHistogram1D("h120","Energy distribution of secondary deuterons", 70, 0., 70.);
        histo13 = createHistogram1D("h130","Energy distribution of secondary tritons", 70, 0., 70.);
        histo14 = createHistogram1D("h140","Energy distribution of secondary alpha particles", 70, 0., 70.);
        histo15 = createHistogram1D("heliumEnergyAfterPhantom","Energy distribution of secondary helium fragments after the phantom",
                70, 0., 500.);
        histo16 = createHistogram1D("hydrogenEnergyAfterPhantom","Energy distribution of secondary helium fragments after the phantom",
                70, 0., 500.);

        kinFragNtuple  = new TNtuple("kinFragNtuple", 
                "Kinetic energy by voxel & fragment", 
                "i:j:k:A:Z:kineticEnergy");
        kineticEnergyPrimaryNtuple= new TNtuple("kineticEnergyPrimaryNtuple", 
                "Kinetic energy by voxel of primary", 
                "i:j:k:kineticEnergy");
        doseFragNtuple = new TNtuple("doseFragNtuple",
                "Energy deposit by voxel & fragment",
                "i:j:k:A:Z:energy");

        fluenceFragNtuple = new TNtuple("fluenceFragNtuple", 
                "Fluence by voxel & fragment",
                "i:j:k:A:Z:fluence");

        letFragNtuple = new TNtuple("letFragNtuple", 
                "Let by voxel & fragment",
                "i:j:k:A:Z:letT:letD");

        theROOTNtuple =   new TNtuple("theROOTNtuple", 
                "Energy deposit by slice",
                "i:j:k:energy");

        theROOTIonTuple = new TNtuple("theROOTIonTuple",
                "Generic ion information",
                "a:z:occupancy:energy");

        fragmentNtuple =  new TNtuple("fragmentNtuple",
                "Fragments",
                "A:Z:energy:posX:posY:posZ");

        metaData =        new TNtuple("metaData",
                "Metadata",
                "events:detectorDistance:waterThickness:beamEnergy:energyError:phantomCenterDistance");
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::FillEnergyDeposit(G4int i,
                                                                                                         G4int j,
                                                                                                         G4int k,
                                                                                                         G4double energy)
{
        if (theROOTNtuple) {
                theROOTNtuple->Fill(i, j, k, energy);
        }
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::BraggPeak(G4int slice, G4double energy)
{
        histo1->SetBinContent(slice, energy); //This uses setbincontent instead of fill to get labels correct
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::SecondaryProtonEnergyDeposit(G4int slice, G4double energy)
{
        histo2->Fill(slice, energy);
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::SecondaryNeutronEnergyDeposit(G4int slice, G4double energy)
{
        histo3->Fill(slice, energy);
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::SecondaryAlphaEnergyDeposit(G4int slice, G4double energy)
{
        histo4->Fill(slice, energy);
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::SecondaryGammaEnergyDeposit(G4int slice, G4double energy)
{
        histo5->Fill(slice, energy);
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::SecondaryElectronEnergyDeposit(G4int slice, G4double energy)
{
        histo6->Fill(slice, energy);
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::SecondaryTritonEnergyDeposit(G4int slice, G4double energy)
{
        histo7->Fill(slice, energy);
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::SecondaryDeuteronEnergyDeposit(G4int slice, G4double energy)
{
        histo8->Fill(slice, energy);
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::SecondaryPionEnergyDeposit(G4int slice, G4double energy)
{
        histo9->Fill(slice, energy);
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::electronEnergyDistribution(G4double energy)
{
        histo10->Fill(energy);
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::gammaEnergyDistribution(G4double energy)
{
        histo11->Fill(energy);
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::deuteronEnergyDistribution(G4double energy)
{
        histo12->Fill(energy);
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::tritonEnergyDistribution(G4double energy)
{
        histo13->Fill(energy);
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::alphaEnergyDistribution(G4double energy)
{
        histo14->Fill(energy);
}
        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::heliumEnergy(G4double secondaryParticleKineticEnergy)
{
        histo15->Fill(secondaryParticleKineticEnergy);
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::hydrogenEnergy(G4double secondaryParticleKineticEnergy)
{
        histo16->Fill(secondaryParticleKineticEnergy);
}

        /////////////////////////////////////////////////////////////////////////////
        // FillKineticFragmentTuple create an ntuple where the voxel indexs, the atomic number and mass and the kinetic
        // energy of all the particles interacting with the phantom, are stored
void HadrontherapyAnalysisManager::FillKineticFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double kinEnergy)
{
    kinFragNtuple -> Fill(i, j, k, A, Z, kinEnergy);
}

        /////////////////////////////////////////////////////////////////////////////
        // FillKineticEnergyPrimaryNTuple creates a ntuple where the voxel indexs and the kinetic
        // energies of ONLY primary particles interacting with the phantom, are stored
void HadrontherapyAnalysisManager::FillKineticEnergyPrimaryNTuple(G4int i, G4int j, G4int k, G4double kinEnergy)
{
    kineticEnergyPrimaryNtuple -> Fill(i, j, k, kinEnergy);
}

        /////////////////////////////////////////////////////////////////////////////
        // This function is called only if ROOT is activated.
        // It is called by the HadrontherapyMatric.cc class file and it is used to create two ntuples containing 
        // the total energy deposited and the fluence values, in each voxel and per any particle (primary beam
        // and secondaries ) 
void HadrontherapyAnalysisManager::FillVoxelFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double energy, G4double fluence)
{
                // Fill the ntuple containing the voxel, mass and atomic number and the energy deposited
    doseFragNtuple ->    Fill( i, j, k, A, Z, energy );
        
                // Fill the ntuple containing the voxel, mass and atomic number and the fluence
        if (i==1 && Z==1) {
                fluenceFragNtuple -> Fill( i, j, k, A, Z, fluence );

        }
}

void HadrontherapyAnalysisManager::FillLetFragmentTuple(G4int i, G4int j, G4int k, G4int A, G4double Z, G4double letT, G4double letD)
{
        letFragNtuple -> Fill( i, j, k, A, Z, letT, letD);

}
        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::FillFragmentTuple(G4int A, G4double Z, G4double energy, G4double posX, G4double posY, G4double posZ)
{
        fragmentNtuple->Fill(A, Z, energy, posX, posY, posZ);
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::genericIonInformation(G4int a,
                                                                                                                 G4double z,
                                                                                                                 G4int electronOccupancy,
                                                                                                                 G4double energy)
{
        if (theROOTIonTuple) {
                theROOTIonTuple->Fill(a, z, electronOccupancy, energy);
        }
}

        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::startNewEvent()
{
        eventCounter++;
}
        /////////////////////////////////////////////////////////////////////////////
void HadrontherapyAnalysisManager::setGeometryMetaData(G4double endDetectorPosition, G4double waterThickness, G4double phantomCenter)
{
        this->detectorDistance = endDetectorPosition;
        this->phantomDepth = waterThickness;
        this->phantomCenterDistance = phantomCenter;
}
void HadrontherapyAnalysisManager::setBeamMetaData(G4double meanKineticEnergy,G4double sigmaEnergy)
{
        this->beamEnergy = meanKineticEnergy;
        this->energyError = sigmaEnergy;
}
/////////////////////////////////////////////////////////////////////////////
// Flush data & close the file
void HadrontherapyAnalysisManager::flush()
{
    if (theTFile)
    {
        theTFile -> Write(); 
        theTFile -> Close();
    }
    eventCounter = 0;
}

#endif