MaGeGeneratorCosmicRayMuons.cc

Go to the documentation of this file.

//
// ********************************************************************
// * DISCLAIMER                                                       *
// *                                                                  *
// * 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.                                                             *
// *                                                                  *
// ********************************************************************
//
// This is supposed to be a generator for cosmic ray muons (with the
// correct spectrum) in the tunnel. It is work in progress
//
// History:
// --------
// 28 Oct 2004   L.Pandola    First implementation (not working)

#include "generators/MaGeGeneratorCosmicRayMuons.hh"
#include "io/MJLogger.hh"

#include "G4PrimaryParticle.hh"
#include "G4Event.hh"
#include "Randomize.hh"
//#include <math.h>
//#include <vector.h>
#include "G4TransportationManager.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4MuonMinus.hh"
#include "G4MuonPlus.hh"
#include "G4DataVector.hh"

MaGeGeneratorCosmicRayMuons::MaGeGeneratorCosmicRayMuons(G4String fileE,
                                                         G4String fileA) : 
  fFileEnergyName(fileE),fFileAngularName(fileA)
{
  //particle initialization
  NumberOfParticlesToBeGenerated = 1;
  particle_definition = G4MuonMinus::MuonMinusDefinition();
  G4ThreeVector zero(0., 0., 0.);
  particle_momentum_direction = G4ParticleMomentum(0., 0., 1.);
  particle_energy = 500*GeV;
  particle_position = zero;
  particle_time = 0.0;
  particle_polarization = zero;
  particle_charge = particle_definition->GetPDGCharge();

  //geometry initialization
  halfz = 10.0*m;
  Radius = 5.0*m;
  verbosityLevel = 0;

  fEnergy = new G4DataVector();
  fEnergySpectrum = new G4DataVector();

  //open data files
  G4String path = "generators/data/";
  OpenEnergyFile(path+fFileEnergyName);
  OpenAngularFile(path+fFileAngularName);
  

  theMessenger = new MaGeGeneratorCosmicRayMuonsMessenger(this);
  gNavigator = G4TransportationManager::GetTransportationManager()
    ->GetNavigatorForTracking();
}

MaGeGeneratorCosmicRayMuons::~MaGeGeneratorCosmicRayMuons()
{
  delete fEnergy;
  delete fEnergySpectrum;
  delete theMessenger;
}

void MaGeGeneratorCosmicRayMuons::SetHalfZ(G4double zhalf)
{
  halfz = zhalf;
}

void MaGeGeneratorCosmicRayMuons::SetRadius(G4double rad)
{
  Radius = rad;
}

void MaGeGeneratorCosmicRayMuons::SampleInitialPosition()
{
  G4ThreeVector RandPos;
  G4double x=0., y=0., z=0.;
  
  //Bisogna generare le posizioni da un cerchio

  RandPos.setX(x);
  RandPos.setY(y);
  RandPos.setZ(z);
  G4ThreeVector CentreCoords(0.,0.,halfz);
  particle_position = CentreCoords + RandPos;
}

void MaGeGeneratorCosmicRayMuons::SetVerbosity(int vL)
{
  verbosityLevel = vL;
  G4cout << "Verbosity Set to: " << verbosityLevel << G4endl;
}

void MaGeGeneratorCosmicRayMuons::GenerateAngularSpectrum()
{;
//bisogna leggere da un file 
}

void MaGeGeneratorCosmicRayMuons::GenerateEnergySpectrum()
{
  G4double totIntegral = (*fEnergySpectrum)[fEnergySpectrum->size()-1]; 
  //final element = total area
  G4double random = G4UniformRand()*totIntegral;
  G4int j=0;
  for (size_t i=0;i<fEnergySpectrum->size()-2;i++)
    {
      //finds the right place
      if (random < (*fEnergySpectrum)[i+1] && random >= (*fEnergySpectrum)[i])
        j=i;
    }
 
  particle_energy=G4UniformRand()
    *((*fEnergy)[j+1]-(*fEnergy)[j])+(*fEnergy)[j];
  return;
}

void MaGeGeneratorCosmicRayMuons::SetParticleDefinition
  (G4ParticleDefinition* aParticleDefinition)
{
  particle_definition = aParticleDefinition;
  particle_charge = particle_definition->GetPDGCharge();
}


void MaGeGeneratorCosmicRayMuons::GeneratePrimaryVertex(G4Event *evt)
{

  if(particle_definition==NULL) {
    G4cout << "No particle has been defined!" << G4endl;
    return;
  }
  
  SampleInitialPosition();

  // Angular stuff
  GenerateAngularSpectrum();

  //Energy stuff
  GenerateEnergySpectrum();
  
  // create a new vertex
  G4PrimaryVertex* vertex = 
    new G4PrimaryVertex(particle_position,particle_time);

  if(verbosityLevel >= 2)
    G4cout << "Creating primaries and assigning to vertex" << G4endl;

  // create new primaries and set them to the vertex
  G4double mass =  particle_definition->GetPDGMass();
  G4double energy = particle_energy + mass;
  G4double pmom = sqrt(energy*energy-mass*mass);
  G4double px = pmom*particle_momentum_direction.x();
  G4double py = pmom*particle_momentum_direction.y();
  G4double pz = pmom*particle_momentum_direction.z();

  if(verbosityLevel >= 1){
    G4cout << "Particle name: " 
           << particle_definition->GetParticleName() << G4endl; 
    G4cout << "       Energy: "<<particle_energy << G4endl;
    G4cout << "     Position: "<<particle_position<< G4endl; 
    G4cout << "    Direction: "<<particle_momentum_direction << G4endl;
    G4cout << " NumberOfParticlesToBeGenerated: "
           << NumberOfParticlesToBeGenerated << G4endl;
  }

  if (G4UniformRand() > 0.5) 
    particle_definition = G4MuonMinus::MuonMinusDefinition();
  else
    particle_definition = G4MuonPlus::MuonPlusDefinition();

  for( G4int i=0; i<NumberOfParticlesToBeGenerated; i++ ) {
    G4PrimaryParticle* particle =
      new G4PrimaryParticle(particle_definition,px,py,pz);
    particle->SetMass( mass );
    particle->SetCharge( particle_charge );
    particle->SetPolarization(particle_polarization.x(),
                              particle_polarization.y(),
                              particle_polarization.z());
    vertex->SetPrimary( particle );
  }

  evt->AddPrimaryVertex( vertex );
  if(verbosityLevel > 1)
    G4cout << " Primary Vetex generated "<< G4endl;   
}

void MaGeGeneratorCosmicRayMuons::OpenEnergyFile(G4String filename)
{
  G4DataVector* temporary = new G4DataVector();

  fEnergyInputFile.open(filename);
  if (!(fEnergyInputFile.is_open())) //not open correctly
    {
      MJLog(fatal) << "File could not be open" << endlog;
    }
  // The data format must be:
  // energy (GeV), spectrum
  // with evenly-spaced bins
  G4int nbin,i;
  G4double dat1, dat2;
  for (i=0; ;i++)
    {
      if (fEnergyInputFile.eof()) break;
      fEnergyInputFile >> dat1 >> dat2;
      fEnergy->push_back(dat1*GeV);
      temporary->push_back(dat2);
    }
  
  if (fEnergy->size() < 3) 
    MJLog(error) << "Invalid data file" << endlog;
  
  fEnergyInputFile.close();
  
  G4double binWidth = (*fEnergy)[1]-(*fEnergy)[0];
  fEnergy->push_back((*fEnergy)[fEnergy->size()-1]+0.5*binWidth); //
  //right edge of the last energy bin

  //Calculate cumulative distribution
  G4double cumul =0.0;
  for (i=0;i<temporary->size();i++) {
    fEnergySpectrum->push_back(cumul);
    cumul += (*temporary)[i];
  }
  fEnergySpectrum->push_back(cumul); //last one
  delete temporary;
}

void MaGeGeneratorCosmicRayMuons::OpenAngularFile(G4String filename)
{;}

---