In This Package:

Li9He8Decayerator.cc

Go to the documentation of this file.

#include "Li9He8Decayerator.h"
#include "TVector3.h"
#include "GaudiKernel/IRndmGenSvc.h"
#include "GaudiKernel/ParticleProperty.h"

#include "HepMC/GenVertex.h"
#include "HepMC/GenEvent.h"
#include "CLHEP/Units/SystemOfUnits.h"
#include <CLHEP/Units/PhysicalConstants.h>

#include <sstream>

using namespace std;
using namespace CLHEP;

Li9He8Decayerator::Li9He8Decayerator(const std::string& type,
                             const std::string& name, 
                             const IInterface* parent)
    : GaudiTool(type,name,parent)
{
  declareInterface<IHepMCEventMutator>(this);
  
  declareProperty("Li9fraction", m_Li9fraction=0.90, 
                  "Li9 fraction Li9/(Li9+He8)");
  declareProperty("CompleteDecay",m_completedecay=0,
                  "Complete Decay include decay channels without neutron");
  
}

Li9He8Decayerator::~Li9He8Decayerator()
{
}

HepMC::GenParticle* Li9He8Decayerator::make_particle(TVector3& vect, 
                                                     double mass, int pid)
{
  double energy = sqrt(vect.Mag2()+mass*mass);
  HepMC::FourVector fourmom(vect.x()/MeV,vect.y()/MeV,vect.z()/MeV,
                            energy/MeV);
  //the new GenParticle will be deleted by GenVertex destructor
  return new HepMC::GenParticle(fourmom,pid,1/*=status*/);
}

StatusCode Li9He8Decayerator::mutate(HepMC::GenEvent& event)
{
    TVector3 pElectron(0,0,0);
    TVector3 pNeutron(0,0,0); 
    TVector3 pGamma(0,0,0); 
    TVector3 pAlpha1(0,0,0); 
    TVector3 pAlpha2(0,0,0);

    // PDG code for Ions (Geant4 version)
    // Nuclear codes are given as 10-digit numbers +-100ZZZAAAI.
    //For a nucleus consisting of np protons and nn neutrons
    // A = np + nn and Z = np.
    // I gives the isomer level, with I = 0 corresponding
    // to the ground state and I >0 to excitations

    int mo_pid = 0;
    Li9He8 decay(rgs);
    if( m_uni() < m_Li9fraction ){//Li9 decay
      mo_pid = 1000030090;
      decay.Li9Decay(pElectron, pNeutron, pAlpha1, pAlpha2, 
                     alpha_mass, m_completedecay);
    }
    else{
      mo_pid = 1000020080;
      decay.He8Decay(pElectron, pNeutron, pGamma, m_completedecay);
    }

    //the new GenVertex will be deleted by GenEvent destructor
    HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(0,0,0,0));
    event.add_vertex(vertex);
    event.set_signal_process_vertex(vertex);

    int pid = 0;
    if( pElectron.Mag2()>0 ){
      pid = 11;
      HepMC::GenParticle* particle = make_particle(pElectron, 
                                                   electron_mass_c2, pid);
      vertex->add_particle_out(particle);
    }
    if( pNeutron.Mag2()>0 ){
      pid = 2112;
      HepMC::GenParticle* particle = make_particle(pNeutron, 
                                                   neutron_mass_c2, pid);
      vertex->add_particle_out(particle);
    }
    if( pGamma.Mag2()>0 ){
      pid = 22;
      HepMC::GenParticle* particle = make_particle(pGamma, 0, pid);
      vertex->add_particle_out(particle);
    }
    if( pAlpha1.Mag2()>0 ){
      pid = 1000020040;
      HepMC::GenParticle* particle = make_particle(pAlpha1, alpha_mass, pid);
      vertex->add_particle_out(particle);
    }
    if( pAlpha2.Mag2()>0 ){
      pid = 1000020040;
      HepMC::GenParticle* particle = make_particle(pAlpha2, alpha_mass, pid);
      vertex->add_particle_out(particle);
    }

    HepMC::GenParticle* particle_in = 
      new HepMC::GenParticle(HepMC::FourVector(0,0,0,0),mo_pid,2);
    // Save mother particle as info particle (default status=2)
    vertex->add_particle_in(particle_in);

    return StatusCode::SUCCESS;
}


StatusCode Li9He8Decayerator::initialize()
{
  //random number
  if (service("RndmGenSvc",rgs,true).isFailure()) {
    fatal() << "Failed to get random service" << endreq;
    return StatusCode::FAILURE;        
  }
  StatusCode sc = m_uni.initialize(rgs, Rndm::Flat(0,1));
  if (sc.isFailure()) {
    throw GaudiException("Failed to initialize uniform random numbers",
                         "GtDecayerator::DecayRates",StatusCode::FAILURE);
  }
  //particle property
  m_pps = svc<IParticlePropertySvc>( "ParticlePropertySvc", true ) ;
  ParticleProperty*  ppp = m_pps->findByStdHepID( 1000020040 );//alpha
  if ( ppp ) {
    alpha_mass = ppp->mass();
  } else {
    alpha_mass = 3727.379;
  }

  return StatusCode::SUCCESS;
}

StatusCode Li9He8Decayerator::finalize()
{
  return StatusCode::SUCCESS;
}

---