CML2ExpVoxels Class Reference

#include <ML2ExpVoxels.hh>

Collaboration diagram for CML2ExpVoxels:

Public Member Functions
	CML2ExpVoxels (G4bool bHasExperimentalData, G4int saving_in_Selected_Voxels_every_events, G4int seed, G4String FileExperimentalData, G4String FileExperimentalDataOut)
	~CML2ExpVoxels (void)
void	add (G4ThreeVector pos, G4double depEnergy, G4double density)
void	add (const G4Step *aStep)
std::vector< Svoxel >	getVoxels ()
G4int	getMinNumberOfEvents ()
G4int	getMaxNumberOfEvents ()
G4bool	loadData ()
void	setRecycling (int recycling)
void	saveResults (void)
void	resetNEventsInVoxels ()

Private Member Functions
void	saveHeader ()
void	calculateNormalizedEd (std::vector< Svoxel > &voxels)

Private Attributes
std::vector< Svoxel >	voxels
G4int *	nVoxelsgeometry
G4ThreeVector	minZone
G4ThreeVector	maxZone
G4int	nCurves
G4int *	startCurve
G4int *	stopCurve
G4double *	chi2Factor
G4String	headerText1
G4String	headerText2
G4String	fullFileIn
G4String	fullFileOut
G4String	seedName
G4String	loopName
SGeneralData *	generalData
G4int	nParticle
G4int	nTotalEvents
G4int	saving_in_Selected_Voxels_every_events
G4int	nRecycling
G4bool	bHasExperimentalData

Detailed Description

Definition at line 51 of file ML2ExpVoxels.hh.

Constructor & Destructor Documentation

CML2ExpVoxels()

CML2ExpVoxels::CML2ExpVoxels	(	G4bool	bHasExperimentalData,
		G4int	saving_in_Selected_Voxels_every_events,
		G4int	seed,
		G4String	FileExperimentalData,
		G4String	FileExperimentalDataOut
	)

Definition at line 48 of file ML2ExpVoxels.cc.

                                                                          :startCurve(0),
          stopCurve(0),chi2Factor(0)
{
    char a[10];
    sprintf(a,"%d", seed);
    seedName=(G4String)a;
    saving_in_Selected_Voxels_every_events=saveEvents;
    nRecycling=1;

    
    fullFileOut=FileExperimentalDataOut+seedName+".m";
    fullFileIn=FileExperimentalData;
    nParticle=nTotalEvents=0;

// define the extremes of global-volume containing all experimental voxels
    G4double extr=100000000000.;
    minZone.set(extr, extr, extr);
    maxZone.set(-extr, -extr, -extr);
    bHasExperimentalData=bData;
}

Here is the call graph for this function:

~CML2ExpVoxels()

CML2ExpVoxels::~CML2ExpVoxels

(

void

)

Definition at line 70 of file ML2ExpVoxels.cc.

{
    delete [] startCurve;
    delete [] stopCurve;
    delete [] chi2Factor;
        delete [] nVoxelsgeometry;
}

Member Function Documentation

add() [1/2]

void CML2ExpVoxels::add	(	G4ThreeVector	pos,
		G4double	depEnergy,
		G4double	density
	)

Here is the caller graph for this function:

add() [2/2]

void CML2ExpVoxels::add

(

const G4Step *

aStep

)

Definition at line 165 of file ML2ExpVoxels.cc.

{
    G4ThreeVector pos;
    G4double depEnergy, density, voxelVolume;
    
    pos=aStep->GetPreStepPoint()->GetPosition();
    depEnergy=aStep->GetTotalEnergyDeposit();
    density=aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial()->GetDensity();

    G4ThreeVector minPos, maxPos;
    G4bool newEvent=false;
    G4double voxelMass, dose;

// check if the event is inside the global-volume 
    if (minZone.getX()<= pos.getX() && pos.getX()<maxZone.getX() && 
        minZone.getY()<= pos.getY() && pos.getY()<maxZone.getY() && 
        minZone.getZ()<= pos.getZ() && pos.getZ()<maxZone.getZ())
    {
// look for the voxel containing the event
        for (int i=0; i<(int)voxels.size(); i++)
        {
            minPos=voxels[i].pos-voxels[i].halfSize;
            maxPos=voxels[i].pos+voxels[i].halfSize;
            if (minPos.getX()<= pos.getX() && pos.getX()<maxPos.getX() && 
                minPos.getY()<= pos.getY() && pos.getY()<maxPos.getY() && 
                minPos.getZ()<= pos.getZ() && pos.getZ()<maxPos.getZ())
            {
                voxelVolume=voxels[i].halfSize.getX()*voxels[i].halfSize.getY()*voxels[i].halfSize.getZ()*8.;
                voxelMass=density*voxelVolume;
// calculate the dose 
                dose=depEnergy/(voxelMass*nRecycling);
                voxels[i].nEvents++;
                                nVoxelsgeometry[i]++;
                voxels[i].depEnergy+=dose;
                voxels[i].depEnergy2+=dose*dose;
                newEvent=true;

                Sparticle *particle=new Sparticle;
                particle->dir=aStep->GetPreStepPoint()->GetMomentumDirection();
                particle->pos=aStep->GetPreStepPoint()->GetPosition();
                particle->kinEnergy=dose; // I use the same kinEnergy name to store the dose
                particle->nPrimaryPart=-1;
                particle->partPDGE=aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
                particle->primaryParticlePDGE=-1;
                particle->volumeId=i; // voxel index where the dose is accumulating
                particle->volumeName="-1";
            }
        }
        if (newEvent)
        {
// save data
            nTotalEvents++;
            if (nTotalEvents%saving_in_Selected_Voxels_every_events==0 && nTotalEvents>0)
            {
                saveResults();
            }
        }
    }
}

Here is the call graph for this function:

calculateNormalizedEd()

void CML2ExpVoxels::calculateNormalizedEd ( std::vector< Svoxel > &  voxels )

private

Definition at line 284 of file ML2ExpVoxels.cc.

{
    int i,j;
    G4double cs, cc;
    int n;
    G4double d2, dd;
        G4double v;
    for (j=0;j<nCurves;j++)
    {
        cs=cc=0.;
        for (i=startCurve[j]-1;i<stopCurve[j];i++)
        {
            cs+=vox[i].depEnergy*vox[i].expDose;
            cc+=vox[i].depEnergy*vox[i].depEnergy;
        }
        if (cc>0.)
        {
            chi2Factor[j]=cs/cc; 
        }
        for (i=startCurve[j]-1;i<stopCurve[j];i++)
        {
            dd=vox[i].depEnergy*vox[i].depEnergy;
            d2=vox[i].depEnergy2;
            n=vox[i].nEvents;
            vox[i].depEnergyNorm=chi2Factor[j]*vox[i].depEnergy;
            v=n*d2-dd;
            if (v<0.){v=0;}
            if (n>1){vox[i].depEnergyNormError=chi2Factor[j]*std::sqrt(v/(n-1));}
            if (n==1){vox[i].depEnergyNormError=vox[i].depEnergyNorm;}
        }
    }
}

Here is the caller graph for this function:

getMaxNumberOfEvents()

G4int CML2ExpVoxels::getMaxNumberOfEvents

(

)

Definition at line 234 of file ML2ExpVoxels.cc.

{
        int n=nVoxelsgeometry[0];
    for (int i=0;i<(int)voxels.size();i++)
    { 
                if (n<nVoxelsgeometry[i]){n = nVoxelsgeometry[i];}
    }
    return n;
}

Here is the caller graph for this function:

getMinNumberOfEvents()

G4int CML2ExpVoxels::getMinNumberOfEvents

(

)

Definition at line 225 of file ML2ExpVoxels.cc.

{
    int n=voxels[0].nEvents;
    for (int i=0;i<(int)voxels.size();i++)
    { 
        if (n>voxels[i].nEvents){n = voxels[i].nEvents;}
    }
    return n;
}

Here is the caller graph for this function:

getVoxels()

std::vector<Svoxel> CML2ExpVoxels::getVoxels ( )

inline

Definition at line 59 of file ML2ExpVoxels.hh.

{return voxels;}

Here is the call graph for this function:

loadData()

G4bool CML2ExpVoxels::loadData

(

void

)

Definition at line 77 of file ML2ExpVoxels.cc.

{
    bHasExperimentalData=true;
    std::ifstream in;

    Svoxel voxel; voxel.volumeId=0;
    G4ThreeVector pos, halfSize;
    G4double expDose;

    in.open(fullFileIn, std::ios::in);
    if (in)
    {
        G4String appo;
        char a[1000];
        in.getline(a,1000,'\n'); headerText1=(G4String)a;
        in.getline(a,1000,'\n');
        in >> nCurves;
        startCurve=new G4int[nCurves];
        stopCurve=new G4int[nCurves]; 
        chi2Factor=new G4double[nCurves];
        for (int i=0; i< nCurves; i++)
        {
            chi2Factor[i]=0.;
            in >> startCurve[i]; 
            in >> stopCurve[i]; 
            in >> chi2Factor[i];
        }
        in.getline(a,1000,'\n');
        in.getline(a,1000,'\n'); headerText2=(G4String)a;

        while (!in.eof())
        {
            in >> pos;
            in >> halfSize; 
            if (bHasExperimentalData)
            {
                in >> expDose;
                voxel.expDose=expDose/100.*(joule/kg); // input data in cGy
            }
            else
            {
                voxel.expDose=0.;
            }
            voxel.pos=pos;
            voxel.halfSize=halfSize;
            voxel.depEnergy=0.;
            voxel.depEnergy2=0.;
            voxel.nEvents=0;
            voxel.depEnergyNorm=0.;
            voxel.depEnergyNormError=0.;
            voxels.push_back(voxel);

// calculate the actual extremes of the global-volume containing all the experimental data
            if (minZone.getX()>pos.getX()-halfSize.getX())
            {minZone.setX(pos.getX()-halfSize.getX());}
            if (maxZone.getX()<pos.getX()+halfSize.getX())
            {maxZone.setX(pos.getX()+halfSize.getX());}

            if (minZone.getY()>pos.getY()-halfSize.getY())
            {minZone.setY(pos.getY()-halfSize.getY());}
            if (maxZone.getY()<pos.getY()+halfSize.getY())
            {maxZone.setY(pos.getY()+halfSize.getY());}

            if (minZone.getZ()>pos.getZ()-halfSize.getZ())
            {minZone.setZ(pos.getZ()-halfSize.getZ());}
            if (maxZone.getZ()<pos.getZ()+halfSize.getZ())
            {maxZone.setZ(pos.getZ()+halfSize.getZ());}

        }
    }
    else
    {
        std::cout << "ERROR I can't find the experimental data file" << G4endl;
        return false;
    }
    in.close();

        nVoxelsgeometry=new G4int[(G4int) voxels.size()];
        resetNEventsInVoxels();

    return true;
}

Here is the call graph for this function:

Here is the caller graph for this function:

resetNEventsInVoxels()

void CML2ExpVoxels::resetNEventsInVoxels

(

)

Definition at line 159 of file ML2ExpVoxels.cc.

{
    for (int i=0; i<(int) voxels.size(); i++ )
    {nVoxelsgeometry[i]=0;}
}

Here is the caller graph for this function:

saveHeader()

void CML2ExpVoxels::saveHeader ( )

private

Definition at line 243 of file ML2ExpVoxels.cc.

{
    std::ofstream out;
    out.open(fullFileOut, std::ios::out);
    out <<"% "<< headerText1 << G4endl;
    out <<"n"<< seedName<<"="<< nCurves<<";" << G4endl;
    out <<"fh"<< seedName<<"=["<< G4endl;
    for (int i=0; i< nCurves; i++)
    {
        out << startCurve[i] << '\t';
        out << stopCurve[i] << '\t';
        out << chi2Factor[i]<< G4endl;
    }
    out << "];"<<G4endl;
    out <<"% x [mm], y [mm], z [mm], Dx [mm], Dy [mm], Dz [mm], expDose [Gy], Calculated dose [Gy], Calculated dose2 [Gy^2], nEvents, normDose [Gy], normDoseError [Gy]";
    out << G4endl;
    out.close();    
}

Here is the caller graph for this function:

saveResults()

void CML2ExpVoxels::saveResults

(

void

)

Definition at line 262 of file ML2ExpVoxels.cc.

{
    if (nTotalEvents>0)
    {
        calculateNormalizedEd(voxels);
        saveHeader();
        std::ofstream out;
        out.open(fullFileOut, std::ios::app);
        out <<"d"<< seedName<<"=["<< G4endl;
        for (int i=0; i<(int)voxels.size(); i++)
        {
            out <<voxels[i].pos.getX()/mm<<'\t'<<voxels[i].pos.getY()/mm<<'\t'<<voxels[i].pos.getZ()/mm<<'\t';
            out <<voxels[i].halfSize.getX()/mm<<'\t'<<voxels[i].halfSize.getY()/mm<<'\t'<<voxels[i].halfSize.getZ()/mm<<'\t';
            out <<voxels[i].expDose/(joule/kg)<<'\t'<<voxels[i].depEnergy/(joule/kg)<<'\t'<<voxels[i].depEnergy2/((joule/kg)*(joule/kg))<<'\t'<<voxels[i].nEvents<< '\t';
            out <<voxels[i].depEnergyNorm/(joule/kg)<<'\t'<<voxels[i].depEnergyNormError/(joule/kg);
            out << G4endl;
        }
        out << "];"<<G4endl;
        out.close();
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

setRecycling()

void CML2ExpVoxels::setRecycling ( int  recycling )

inline

Definition at line 66 of file ML2ExpVoxels.hh.

{nRecycling=recycling;}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

bHasExperimentalData

G4bool CML2ExpVoxels::bHasExperimentalData

private

Definition at line 83 of file ML2ExpVoxels.hh.

chi2Factor

G4double* CML2ExpVoxels::chi2Factor

private

Definition at line 77 of file ML2ExpVoxels.hh.

fullFileIn

G4String CML2ExpVoxels::fullFileIn

private

Definition at line 78 of file ML2ExpVoxels.hh.

fullFileOut

G4String CML2ExpVoxels::fullFileOut

private

Definition at line 78 of file ML2ExpVoxels.hh.

generalData

SGeneralData* CML2ExpVoxels::generalData

private

Definition at line 80 of file ML2ExpVoxels.hh.

headerText1

G4String CML2ExpVoxels::headerText1

private

Definition at line 78 of file ML2ExpVoxels.hh.

headerText2

G4String CML2ExpVoxels::headerText2

private

Definition at line 78 of file ML2ExpVoxels.hh.

loopName

G4String CML2ExpVoxels::loopName

private

Definition at line 79 of file ML2ExpVoxels.hh.

maxZone

G4ThreeVector CML2ExpVoxels::maxZone

private

Definition at line 74 of file ML2ExpVoxels.hh.

minZone

G4ThreeVector CML2ExpVoxels::minZone

private

Definition at line 74 of file ML2ExpVoxels.hh.

nCurves

G4int CML2ExpVoxels::nCurves

private

Definition at line 75 of file ML2ExpVoxels.hh.

nParticle

G4int CML2ExpVoxels::nParticle

private

Definition at line 81 of file ML2ExpVoxels.hh.

nRecycling

G4int CML2ExpVoxels::nRecycling

private

Definition at line 82 of file ML2ExpVoxels.hh.

nTotalEvents

G4int CML2ExpVoxels::nTotalEvents

private

Definition at line 82 of file ML2ExpVoxels.hh.

nVoxelsgeometry

G4int* CML2ExpVoxels::nVoxelsgeometry

private

Definition at line 73 of file ML2ExpVoxels.hh.

saving_in_Selected_Voxels_every_events

G4int CML2ExpVoxels::saving_in_Selected_Voxels_every_events

private

Definition at line 82 of file ML2ExpVoxels.hh.

seedName

G4String CML2ExpVoxels::seedName

private

Definition at line 79 of file ML2ExpVoxels.hh.

startCurve

G4int* CML2ExpVoxels::startCurve

private

Definition at line 76 of file ML2ExpVoxels.hh.

stopCurve

G4int * CML2ExpVoxels::stopCurve

private

Definition at line 76 of file ML2ExpVoxels.hh.

voxels

std::vector<Svoxel> CML2ExpVoxels::voxels

private

Definition at line 72 of file ML2ExpVoxels.hh.

---

The documentation for this class was generated from the following files:

• ML2ExpVoxels.hh
• ML2ExpVoxels.cc