ML2PrimaryGenerationAction.cc

Go to the documentation of this file.

//
// ********************************************************************
// * 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.          *
// ********************************************************************
//
// The code was written by :
//      ^Claudio Andenna  claudio.andenna@ispesl.it, claudio.andenna@iss.infn.it
//      *Barbara Caccia barbara.caccia@iss.it
//      with the support of Pablo Cirrone (LNS, INFN Catania Italy)
//      with the contribute of Alessandro Occhigrossi*
//
// ^INAIL DIPIA - ex ISPESL and INFN Roma, gruppo collegato Sanità, Italy
// *Istituto Superiore di Sanità and INFN Roma, gruppo collegato Sanità, Italy
//  Viale Regina Elena 299, 00161 Roma (Italy)
//  tel (39) 06 49902246
//  fax (39) 06 49387075
//
// more information:
// http://g4advancedexamples.lngs.infn.it/Examples/medical-linac
//
//*******************************************************//

#include "ML2PrimaryGenerationAction.hh"

using namespace CLHEP;

CML2PrimaryGenerationAction::CML2PrimaryGenerationAction(void)
:particleGun(0),gamma(0),electron(0),positron(0),primaryParticleData(0),particles(0)
{
}
CML2PrimaryGenerationAction* CML2PrimaryGenerationAction::instance = 0;

CML2PrimaryGenerationAction* CML2PrimaryGenerationAction::GetInstance(void)
{
  if (instance == 0)
    {
      instance = new CML2PrimaryGenerationAction();
    }
  return instance;
}
void CML2PrimaryGenerationAction::inizialize(SPrimaryParticle *pData)
{
        rm=new G4RotationMatrix();
        PrimaryGenerationActionMessenger=new CML2PrimaryGenerationActionMessenger(this);
        particle=new Sparticle;
        nParticle=nPhSpParticles=nRandomParticles=0;

        G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();

        gamma=particleTable->FindParticle("gamma");
        electron=particleTable->FindParticle("e-");
        positron=particleTable->FindParticle("e+");
        particleGun=new G4ParticleGun();

        primaryParticleData=pData;
        primaryParticleData->nPrimaryParticle=0;
        primaryParticleData->partPDGE=0;
}

void CML2PrimaryGenerationAction::design(G4double aTZ)
{
        accTargetZPosition=aTZ;
        switch (idParticleSource)
        {
        case id_randomTarget:
                setGunRandom();
                break;
        case id_phaseSpace:
                setGunCalculatedPhaseSpace();
                break;
        }
}

void CML2PrimaryGenerationAction::setGunRandom()
{
        particleGun->SetParticleDefinition(electron);
        particleGun->SetNumberOfParticles(1);
        idCurrentParticleSource=idParticleSource;
}

void CML2PrimaryGenerationAction::setGunCalculatedPhaseSpace()
{
        particles=new Sparticle[nMaxParticlesInRamPhaseSpace];
        particleGun->SetNumberOfParticles(1);
        idCurrentParticleSource=idParticleSource;
}

CML2PrimaryGenerationAction::~CML2PrimaryGenerationAction(void)
{
        delete particleGun;
        delete [] particles;
        delete particles;
}
void CML2PrimaryGenerationAction::GeneratePrimaries(G4Event *anEvent)
{
        static int currentRecycle=nRecycling;
        static G4ThreeVector pos0, dir0;
        if (currentRecycle==nRecycling)
        {
                currentRecycle=0;
                switch (idCurrentParticleSource)
                {
                case id_randomTarget:
                                GenerateFromRandom();
                        break;
                case id_phaseSpace:
                                GenerateFromCalculatedPhaseSpace();
                        break;
                }
                pos0=pos;
                dir0=dir;
        }
        currentRecycle++;
        pos=pos0;
        dir=dir0;
        applySourceRotation(); // to follow the accelerator rotation

        primaryParticleData->partPDGE=particleGun->GetParticleDefinition()->GetPDGEncoding();
        primaryParticleData->nPrimaryParticle++;

        particleGun->SetParticleEnergy(ek*MeV);
        particleGun->SetParticlePosition(pos*mm);
        particleGun->SetParticleMomentumDirection((G4ParticleMomentum)dir);
        particleGun->GeneratePrimaryVertex(anEvent);
}
void CML2PrimaryGenerationAction::GenerateFromRandom()
{
        sinTheta=RandGauss::shoot(0., 0.003);
        cosTheta=std::sqrt(1 - sinTheta*sinTheta);
        phi=twopi*G4UniformRand();
        dir.set(sinTheta*std::cos(phi), sinTheta*std::sin(phi), cosTheta);

        ro=G4UniformRand()*GunRadious; 
        alfa=G4UniformRand()*twopi;

        pos.setX(ro*std::sin(alfa));
        pos.setY(ro*std::cos(alfa));
        pos.setZ(-(accTargetZPosition +5.)*mm); // the primary electrons are generated 5 mm before the target
        ek=RandGauss::shoot(GunMeanEnegy, GunStdEnegy);
        nRandomParticles++;
}
void CML2PrimaryGenerationAction::GenerateFromCalculatedPhaseSpace()
{
        static bool bFirstTime=true;
        if (bFirstTime)
        {bFirstTime=false;fillParticlesContainer();}
        if (nParticle==nMaxParticlesInRamPhaseSpace) // once all the particles stored in RAM hae been processed a new set is loaded
        {
                fillParticlesContainer();
                nParticle=0;
        }

        pos=particles[nParticle].pos;
        dir=particles[nParticle].dir;
        ek=particles[nParticle].kinEnergy;
        switch (particles[nParticle].partPDGE)
        {
        case -11:
                particleGun->SetParticleDefinition(positron);
                break;
        case 11:
                particleGun->SetParticleDefinition(electron);
                break;
        case 22:
                particleGun->SetParticleDefinition(gamma);
                break;
        }
        nPhSpParticles++;
        nParticle++;
}
void CML2PrimaryGenerationAction::applySourceRotation()
{
        pos=*rm*pos;
        dir=*rm*dir;
}
void CML2PrimaryGenerationAction::fillParticlesContainer()
{
        static int currentFilePosition=0;
        static int currentFileSize=0;
        int startDataFilePosition;
        std::ifstream in;
        in.open(calculatedPhaseSpaceFileIN, std::ios::in);
        if (in)
        {
                std::cout <<"ERROR phase space file: "<< calculatedPhaseSpaceFileIN << " NOT found. Run abort "<< G4endl;
                G4RunManager::GetRunManager()->AbortRun(true);
        }

        static bool bFirstTime=true;
        if (bFirstTime)
        {
                in.seekg(-1,std::ios::end);
                currentFileSize=in.tellg();
                in.seekg(0,std::ios::beg);
                bFirstTime=false;
        }

        char a[1000];
        in.getline(a,1000,'\n');
        in.getline(a,1000,'\n');
        startDataFilePosition=in.tellg();
        if (currentFilePosition>0)
        {in.seekg(currentFilePosition, std::ios::beg);}
        int i;
        G4double x,y,z;
        G4int d;

        static bool checkFileRewind=false;
        static bool bRewindTheFile=false;
        static int nPhSpFileRewind=0;

        for (i=0;i<nMaxParticlesInRamPhaseSpace;i++)
        {
                if (bRewindTheFile) // to read the phase space file again to fill the container
                {
                        in.close();
                        in.open(calculatedPhaseSpaceFileIN, std::ios::in);
                        in.seekg(startDataFilePosition, std::ios::beg);
                        checkFileRewind=true;
                        bRewindTheFile=false;
                        std::cout<<"\n################\nI have reached the end of the phase space file "<<++nPhSpFileRewind <<" times, I rewind the file\n" << G4endl;
                        std::cout <<"loaded " <<i <<"/"<< nMaxParticlesInRamPhaseSpace<<" particles" << G4endl;
                }
                in >> d; 
                in >> x; in >>y; in >> z; 
/*                      std::cout <<"x:" <<x << G4endl;
                        std::cout <<"y:" <<y << G4endl;
                        std::cout <<"z:" <<z << G4endl;*/
                particles[i].pos.set(x,y,z-accTargetZPosition);
                in >> x; in >>y; in >> z; 
                particles[i].dir.set(x,y,z);
                in >> x; 
                particles[i].kinEnergy=x;
                in >> d; 
                particles[i].partPDGE=d;
                in >> d; in >> d; 
                if (in.eof())   {bRewindTheFile=true;}
                if (checkFileRewind)    {checkFileRewind=false;}
        }
        std::cout <<"loaded " <<i <<"/"<< nMaxParticlesInRamPhaseSpace<<" particles" << G4endl;
        currentFilePosition=in.tellg(); // to remind the actual position in the phase space file
        if (currentFilePosition>=currentFileSize) // to read the phase space file again
        {currentFilePosition=startDataFilePosition;} 
        in.close();
}