BesPrimaryGeneratorAction.cc

Go to the documentation of this file.

//---------------------------------------------------------------------------//
//      BOOST --- BESIII Object_Oriented Simulation Tool                     //
//---------------------------------------------------------------------------//
// Description: Mimic the BES generator, TESTER
// Author: Liuhm
// Created: June, 2003
// Modified: Liuhm, 7 April, 2004, re-randomize phi distribution
// Comment: Use Geant4 ParticleGun
//---------------------------------------------------------------------------//
 
#include "GenSim/BesPrimaryGeneratorAction.hh"
#include "G4Event.hh"
#include "G4HEPEvtInterface.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleGun.hh"
#include "G4ParticleMomentum.hh"
#include "G4ParticleTable.hh"
#include "GenSim/BesPrimaryGeneratorMessenger.hh"
#include "Randomize.hh"
#include "TFile.h"
#include "TH1F.h"
#include "globals.hh"
 
using namespace CLHEP;
 
BesPrimaryGeneratorAction::BesPrimaryGeneratorAction() {
  nParticle    = 1;
  particleName = "pi-";
  minCos       = -0.8;
  maxCos       = 0.8;
  phiStart     = 0.;
  phiEnd       = 360.;
  pMomentum    = 1.;
  deltaP       = 0.;
  posX         = 0;
  posY         = 0;
  posZ         = 0;
 
  particleGun   = new G4ParticleGun( 1 );
  isGenbes      = true;
  HEPEvt        = 0;
  generatorName = "tester";
  // TESTER parameters passed by this messenger
  messenger = new BesPrimaryGeneratorMessenger( this );
 
  if ( generatorName == "cosmic" )
  {
    std::string path = getenv( "GENSIMROOT" );
    G4cout << "path: " << path << G4endl;
 
    path += "/root/";
    std::string pFile     = path + "ppdc.root";
    std::string thetaFile = path + "theta.root";
    std::string phiFile   = path + "phi.root";
 
    TFile* f1 = new TFile( pFile.c_str() );
    h1        = (TH1F*)f1->Get( "htemp" );
 
    TFile* f2 = new TFile( thetaFile.c_str() );
    h2        = (TH1F*)f2->Get( "htemp" );
 
    TFile* f3 = new TFile( phiFile.c_str() );
    h3        = (TH1F*)f3->Get( "htemp" );
 
    // ftest = new TFile("ftest.root","recreate");
    // tuple = new TNtuple("test","test","p:theta:phi");
    // counter = 0;
  }
}
 
BesPrimaryGeneratorAction::~BesPrimaryGeneratorAction() {
  delete particleGun;
  if ( messenger ) delete messenger;
  if ( generatorName == "genbes" ) delete HEPEvt;
}
 
void BesPrimaryGeneratorAction::GeneratePrimaries( G4Event* anEvent ) {
  if ( generatorName == "tester" )
  {
    G4ParticleTable*      particleTable = G4ParticleTable::GetParticleTable();
    G4ParticleDefinition* particle      = particleTable->FindParticle( particleName );
    particleGun->SetParticleDefinition( particle );
    particleGun->SetParticlePosition( G4ThreeVector( posX, posY, posZ ) );
    particleGun->SetParticleTime( 648 * ns );
 
    // set TESTER parameters
    for ( G4int i = 0; i < nParticle; i++ )
    {
      G4double pMag = pMomentum;
      // randomize p
      if ( deltaP > 0. ) pMag = pMomentum - deltaP * ( 1.0 - 2.0 * G4UniformRand() );
      pMag = pMag * GeV;
      // randomize cos(theta)
      G4double costheta = minCos + ( maxCos - minCos ) * G4UniformRand();
      // randomize phi
      G4double phi      = phiStart + ( phiEnd - phiStart ) * G4UniformRand();
      phi               = phi * degree;
      G4double sintheta = sqrt( 1. - costheta * costheta );
      // computer 3-vector momentum
      G4ParticleMomentum aMomentum;
      aMomentum[0] = pMag * sintheta * cos( phi );
      aMomentum[1] = pMag * sintheta * sin( phi );
      aMomentum[2] = pMag * costheta;
      // use ParticleGun to generate event
      particleGun->SetParticleMomentum( aMomentum );
      particleGun->GeneratePrimaryVertex( anEvent );
    }
  }
  else if ( generatorName == "cosmic" )
  {
    G4cout << "generatorName: " << generatorName << G4endl;
    G4ParticleTable*      particleTable = G4ParticleTable::GetParticleTable();
    G4ParticleDefinition* particle      = particleTable->FindParticle( particleName );
    G4cout << "particleName: " << particleName << G4endl;
 
    particleGun->SetParticleDefinition( particle );
    particleGun->SetParticlePosition( G4ThreeVector( posX, posY, posZ ) );
 
    /*std::string path = getenv("GENSIMROOT");
    G4cout<<"path: "<<path<<G4endl;
 
    path += "/root/";
    std::string pFile = path + "ppdc.root";
    std::string thetaFile = path + "theta.root";
    std::string phiFile = path +"phi.root";*/
 
    // TFile* f1 = new TFile(pFile.c_str());
    // H1F* h1 = (TH1F*)f1->Get("htemp");
    G4double pMag = h1->GetRandom() * GeV;
    G4cout << "pMag: " << pMag << G4endl;
    // f1->Close();
 
    // TFile* f2 = new TFile(thetaFile.c_str());
    // TH1F* h2 = (TH1F*)f2->Get("htemp");
    // randomize cos(theta)
    G4double theta = (Double_t)h2->GetRandom();
    G4cout << "theta: " << theta << G4endl;
    G4double costheta = cos( theta );
    // f2->Close();
 
    // TFile* f3 = new TFile(phiFile.c_str());
    // TH1F* h3 = (TH1F*)f3->Get("htemp");
    // randomize phi
    G4double phi = (Double_t)h3->GetRandom();
    G4cout << "phi: " << phi << G4endl;
    // f3->Close();
 
    // f1->Close();
    // f2->Close();
    // f3->Close();
    // ftest->ReOpen("update");
    // tuple->Fill(pMag,theta,phi);
    // counter++;
    // if(counter==2000)
    //   tuple->Write();
 
    G4double sintheta = sqrt( 1. - costheta * costheta );
    // computer 3-vector momentum
    G4ParticleMomentum aMomentum;
    aMomentum[0] = pMag * sintheta * cos( phi );
    aMomentum[1] = pMag * sintheta * sin( phi );
    aMomentum[2] = pMag * costheta;
    // use ParticleGun to generate event
    particleGun->SetParticleMomentum( aMomentum );
    particleGun->GeneratePrimaryVertex( anEvent );
  }
 
  else if ( generatorName == "genbes" )
  {
    G4cout << "genbes called" << G4endl;
    if ( isGenbes )
    {
      isGenbes = false;
      HEPEvt   = new G4HEPEvtInterface( genbesName );
    }
    HEPEvt->SetParticleTime( 648 * ns );
    HEPEvt->GeneratePrimaryVertex( anEvent );
  }
}