MaGeDecay0Interface.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.                                                             *
// *                                                                  *
// ********************************************************************
// History:
// --------
// 25 Oct 2004   L.Pandola    First implementation
// 04 Nov 2004   L.Pandola    Change: not necessary to give immediately the file name

#include "generators/MaGeDecay0Interface.hh"
#include "io/MJLogger.hh"
#include "G4Event.hh"
#include "G4PrimaryVertex.hh"
#include "G4PrimaryParticle.hh"
#include "generators/MaGeDecay0InterfaceMessenger.hh"
#include "G4ThreeVector.hh"
#include "G4ParticleTable.hh"
#include <string.h>

MaGeDecay0Interface::MaGeDecay0Interface(G4String File) :
  fFileName(File)
{
  fTheMessenger = new MaGeDecay0InterfaceMessenger(this);

  G4ThreeVector zero;
  particle_position =  zero;
  particle_time = 0.0;
}

MaGeDecay0Interface::~MaGeDecay0Interface()
{
  delete fTheMessenger;
  if (fInputFile.is_open()) fInputFile.close();
}

void MaGeDecay0Interface::OpenFile()
{
  char header[160];
  fInputFile.open(fFileName);

  if (!(fInputFile.is_open())) {//not open correctly
    MJLog(error) << "File not valid!" << endlog;
    MJLog(fatal) << endlog;
  }
  G4int ifound;
  MJLog(debugging) << "Decay0 file " << fFileName << " found" << endlog;
  //Here one has to read the header
  //Problem: the header has not always the same lenght
  for (G4int j=0;;j++)
    {
      fInputFile.getline(header,160);
      ifound = strncmp (header," First",6);
      if (!ifound) break;
      if (j == 50) {
        MJLog(error) << "Corrupted file? " << endlog;
        MJLog(fatal) << endlog;
      }
    }
  //here one has 2 more lines to skip
  fInputFile.getline(header,160);
  fInputFile.getline(header,160);
  MJLog(trace) << "Decay0 file " << fFileName << " initialized" << endlog;
}

void MaGeDecay0Interface::ChangeFileName(G4String newFile)
{
   if (fFileName != newFile) //check if the new file is equal to the other
    {
      if (fInputFile.is_open()) fInputFile.close(); //close the old file
      fFileName = newFile; 
      OpenFile(); //open the new one
    }
}

void MaGeDecay0Interface::GeneratePrimaryVertex(G4Event* evt)
{
  if (!(fInputFile.is_open())) {
    MJLog(error) << "File must be set!" << endlog;
    MJLog(fatal) << endlog;
  }
  G4int nEntries=0;
  G4int nEvent=0;
  G4double time=0.0;
  G4String particleName;
  G4double px,py,pz,relativeTime;
  G4int particleCode;

  fInputFile >> nEvent >> time >> nEntries;
 
  if (fInputFile.eof())
    {
      fInputFile.close();
      MJLog(error) << "File over: not enough events!" << endlog;
      MJLog(fatal) << endlog;
      return;
    }
  
  
  particle_time = time*s;

  G4PrimaryVertex* vertex = new G4PrimaryVertex (particle_position,
                                                 particle_time);

  G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();

  for (G4int i=0;i<nEntries;i++) //loop of the number of particles per event
    {
      fInputFile >> particleCode >> px >> py >> pz >> relativeTime;
      
      particleName = ConvertCodeToName(particleCode);
     
      //Create primary particle
      G4PrimaryParticle* particle = 
        new G4PrimaryParticle(theParticleTable->FindParticle(particleName),
                              px*MeV,py*MeV,pz*MeV);
     
      particle->SetProperTime(relativeTime*s+particle_time); 

      //Add particle to vertex
      vertex->SetPrimary(particle);
    }

  evt->AddPrimaryVertex(vertex);
}

G4String MaGeDecay0Interface::ConvertCodeToName(G4int code)
{
  if (code == 1)
    return "gamma";
  else if (code == 2)
    return "e+";
  else if (code == 3)
    return "e-";
  else if (code == 47)
    return "alpha";
  else 
    MJLog(error) << "Unknown particle code" << endlog;
  return " ";
}
      
 

---