Mcgpj.cxx

Go to the documentation of this file.

//*****************************************************************************
//
// Mcgpj.cxx
//
// General 2-body generator for e^+e^-, mu^+mu^-, pi^+pi^-, tau^+tau^-,
// K_S K_L, K^+K^-, gamma gamma with precision better 0.2 %
//
// Mar 2009 Original BES3 code by Alexei Sibidanov
//
//*****************************************************************************
 
// Header for this module:-
#include "Mcgpj.h"
 
// Framework Related Headers:-
#include "HepMC/GenEvent.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
 
#include "BesRndmGenSvc/IBesRndmGenSvc.h"
#include "CLHEP/Vector/LorentzVector.h"
 
#include "GaudiKernel/SmartDataPtr.h"
 
#include "GeneratorObject/McGenEvent.h"
 
#include "TEPCrossPart.h"
#include "TGGCrossPart.h"
#include "TKcCrossPart.h"
#include "TKinemCut.h"
#include "TKnCrossPart.h"
#include "TMuCrossPart.h"
#include "TPiCrossPart.h"
#include "TRadCor.h"
#include "TRadGlobal.h"
#include "TRandom3.h"
 
#include "T3piCrossPart.h"
#include "T4piCrossPart.h"
#include "T5piCrossPart.h"
#include "TDFun_o.h"
#include "TPhoton_o.h"
 
using namespace CLHEP;
 
TRadCor*     gRad;
TRadGlobal*  gGlobal;
TKinemCut*   gCut;
TConstants*  gConst;
TVCrossPart* MatrixElements;
TCrossPart*  CS;
//-------------------------
 
DECLARE_COMPONENT( Mcgpj )
Mcgpj::Mcgpj( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ) {
  declareProperty( "Data", m_datapath = "${MCGPJROOT}/data" ); // Path to data files
  declareProperty( "VpolFname",
                   m_vpolfname = "vpol.dat" ); // Vacuum polarization data file name
  declareProperty( "CMEnergy", cmE = 3.097 );  // 2*Ebeam [GeV]
  declareProperty( "Process", proc = 0 );      // process
  declareProperty( "NRad", NRad = 20000 );
  declareProperty( "HardPhoton", IsHardPhoton = 1 );
  declareProperty( "NoVacPol", IsNoVacPol = 0 );
  declareProperty( "FSR", IsFSR = 1 );
  declareProperty( "AcolAngle", pc = 0. );
  declareProperty( "DeltaE", de = -1. );                  // [GeV]
  declareProperty( "NTheta0", nt0 = -1. );                // [1/sqrt(gamma)]
  declareProperty( "DeltaTheta", dt = 0.5 );              //[rad]
  declareProperty( "DeltaPhi", dp = 0.3 );                // [rad]
  declareProperty( "AverageTheta", at = 1.1 );            //[rad]
  declareProperty( "ThetaDetector", td = 1.1 - 0.5 / 2 ); //[rad]
  declareProperty( "AverageMomentum", am = 0.090 );       //[GeV/c]
  declareProperty( "CrossMomentum", cm = 0.090 );         //[GeV/c]
  declareProperty( "MinimumEnergy", em = 0.050 );         //[GeV]
  declareProperty( "ThetaIntermediate", ti = 0.473 );     // [rad]
  declareProperty( "AlphaScale", al = 1. );
  declareProperty( "ThetaMinus", thm = -1. );    //[rad]
  declareProperty( "ThetaPlus", thp = -1. );     //[rad]
  declareProperty( "AbsoluteError", te = 0.05 ); // absolute error in [nb]
  declareProperty( "RelativeError", re = 0.05 ); // realative error
                                                 //  declareProperty("InitialSeed", Seed = 0);
  declareProperty( "BeamSpread", spread = -1 );  // [MeV]
  //  declareProperty("PhaseShift", phase  = -90);// [deg]
  declareProperty( "Mode5pi", m_fmode5pi = 0 ); // 0 - via b1+/- pi-/+, 1 via b10 pi0
}
 
StatusCode Mcgpj::initialize() {
  MsgStream log( msgSvc(), name() );
 
  log << MSG::INFO << "Mcgpj initialize" << endmsg;
 
  static const bool CREATEIFNOTTHERE( true );
  StatusCode        RndmStatus = service( "BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE );
  if ( !RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc )
  {
    log << MSG::ERROR << " Could not initialize Random Number Service" << endmsg;
    return RndmStatus;
  }
  CLHEP::HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine( "Mcgpj" );
  engine->showStatus();
 
  // GeV -> MeV
  cmE *= 1e3;
  de *= 1e3;
  am *= 1e3;
  cm *= 1e3;
  em *= 1e3;
 
  if ( gRandom ) delete gRandom;
  gRandom = new TRandom3();
  gRandom->SetSeed( engine->getSeed() );
  // gRandom->Dump();
 
  if ( proc < 100 )
  {
 
    gConst  = new TConstants();
    gGlobal = new TRadGlobal();
    gCut    = new TKinemCut();
 
    if ( ( proc % 10 ) == 3 )
    {
      td = 0.025;
      dt = 10;
      dp = 10;
      cm = 0;
      am = 0;
    }
 
    const int MaxList = 20;
    bool      InList[MaxList];
    for ( int j = 0; j < MaxList; j++ ) InList[j] = true;
 
    double EBeam = 0.5 * cmE;
 
    if ( ( EBeam < 100 || EBeam > 2500 ) && !IsNoVacPol )
    {
      log << MSG::ERROR << "Invalid value of beam energy:" << EBeam << endmsg;
      return StatusCode::FAILURE;
    }
 
    gGlobal->Set_TotalError( te );
 
    gGlobal->Set_RelativeError( re );
 
    gGlobal->Set_Type( ( proc % 10 ) );
 
    gGlobal->Set_E( EBeam );
 
    gGlobal->Set_ThetaInt( ti );
 
    if ( 0 > td ) gGlobal->Set_ThetaMin( at - dt / 2 );
    else gGlobal->Set_ThetaMin( td );
 
    if ( 0 > de )
    {
      if ( ( proc % 10 ) == 0 ) gGlobal->Set_dE_Abs( 0.015 * EBeam );
      else gGlobal->Set_dE_Abs( 0.0003 * EBeam );
    }
    else gGlobal->Set_dE_Abs( de );
 
    if ( 0 > nt0 )
    {
      if ( proc > 10 ) gGlobal->Set_Theta0_Rel( 0.0 );
      else gGlobal->Set_Theta0_Rel( 1.5 );
    }
    else gGlobal->Set_Theta0_Rel( nt0 );
 
    if ( 0 < thm )
    {
      gCut->ThetaMinM( thm );
      gGlobal->Set_ThetaMin( thm );
    }
 
    if ( 0 < thp ) gCut->ThetaMinP( thp );
 
    gCut->CosPsi( cos( pc ) );
 
    gCut->DeltaTheta( dt );
 
    gCut->AverageTheta( at );
 
    gCut->DeltaPhi( dp );
 
    gCut->PAverage( am );
 
    gCut->PCross( cm );
 
    gCut->EMin( em );
 
    gConst->SetAlphaScale( al );
 
    gConst->Print();
 
    try
    { gGlobal->Init(); } catch ( std::logic_error lge )
    {
      log << MSG::ERROR << lge.what() << endmsg;
      return StatusCode::FAILURE;
    }
 
    gGlobal->SetDatadir( m_datapath );
    gGlobal->SetVpolFname( m_vpolfname );
    //    gGlobal->SetIntFname("");
 
    gGlobal->Print();
 
    gCut->Init();
    gCut->Print();
 
    log << MSG::INFO << "Cross-section statistical precision will be better than "
        << gGlobal->Get_TotalError() << " nb and " << gGlobal->Get_RelativeError() * 100 << "%"
        << endmsg;
 
    if ( !IsHardPhoton )
      log << MSG::INFO << "Hard photon on big angle is not included!" << endmsg;
 
    if ( IsNoVacPol ) log << MSG::INFO << "Vacuum polarization is not included!" << endmsg;
 
    if ( !IsFSR ) log << MSG::INFO << "Final state radiation is not included!" << endmsg;
 
    if ( proc > 10 ) log << MSG::INFO << "Alpha order generation only!" << endmsg;
 
    log << std::flush;
 
    if ( proc == 0 )
    {
      MatrixElements = new TEPCrossPart();
      InList[18]     = false;
    }
    else if ( proc == 1 || proc == 5 ) MatrixElements = new TMuCrossPart();
    else if ( proc == 2 ) MatrixElements = new TPiCrossPart();
    else if ( proc == 3 ) MatrixElements = new TKnCrossPart();
    else if ( proc == 4 ) MatrixElements = new TKcCrossPart();
    else if ( proc == 6 ) MatrixElements = new TGGCrossPart();
    else if ( proc == 10 )
    {
      MatrixElements = new TEPCrossPart();
      for ( int j = 0; j < MaxList; j++ ) InList[j] = false;
      InList[16] = true;
      InList[17] = true;
      InList[18] = true;
    }
    else return StatusCode::FAILURE;
 
    if ( IsNoVacPol ) MatrixElements->SetZeroVP();
 
    if ( !IsFSR ) MatrixElements->SetNoFSR();
 
    gRad = new TRadCor( MatrixElements );
    gRad->SetNEvents( NRad );
    gRad->SetPartList( InList );
    gRad->Init();
 
    MatrixElements->SetHardPhoton( IsHardPhoton );
    gRad->MakeCrossSection();
    if ( ( proc % 10 ) == 2 ) ( (TPiCrossPart*)MatrixElements )->GetFormFactor()->Print();
 
    if ( ( proc % 10 ) == 0 )
    {
      fpid[0] = 11;
      fpid[1] = -11;
      fM      = 0.51099891;
    }
    if ( ( proc % 10 ) == 1 )
    {
      fpid[0] = 13;
      fpid[1] = -13;
      fM      = 105.658367;
    }
    if ( ( proc % 10 ) == 2 )
    {
      fpid[0] = 211;
      fpid[1] = -211;
      fM      = 139.57018;
    }
    if ( ( proc % 10 ) == 3 )
    {
      fpid[0] = 130;
      fpid[1] = 310;
      fM      = 497.614;
    }
    if ( ( proc % 10 ) == 4 )
    {
      fpid[0] = 321;
      fpid[1] = -321;
      fM      = 493.677;
    }
    if ( ( proc % 10 ) == 5 )
    {
      fpid[0] = 15;
      fpid[1] = -15;
      fM      = 1776.84;
    }
    if ( ( proc % 10 ) == 6 )
    {
      fpid[0] = 22;
      fpid[1] = 22;
      fM      = 0;
    }
  }
  else
  {
    double EBeam = 0.5 * cmE;
    if ( 0 > de ) de = 0.005 * EBeam;
 
    if ( 0 > nt0 ) nt0 = 1;
 
    if ( proc == 100 )
    {
      CS = new T3piCrossPart( EBeam, de, nt0 );
      if ( spread > 0 ) CS->SetBeamSpread( spread );
    }
    else if ( proc == 101 )
    {
      CS = new T4piCrossPart( EBeam, de, nt0 );
      if ( spread > 0 ) CS->SetBeamSpread( spread );
    }
    else if ( proc == 102 )
    {
      CS = new T5piCrossPart( EBeam, de, nt0, m_fmode5pi );
      if ( spread > 0 ) CS->SetBeamSpread( spread );
    }
    CS->MakeParts( te );
  }
 
  log << MSG::INFO << "Mcgpj initialize finished" << endmsg;
 
  return StatusCode::SUCCESS;
}
 
StatusCode Mcgpj::execute() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "Mcgpj executing" << endmsg;
  log << MSG::WARNING << "execute start" << endmsg;
 
  // Fill event information
  GenEvent* evt = new GenEvent( 1, 1 );
 
  GenVertex* prod_vtx = new GenVertex();
 
  evt->add_vertex( prod_vtx );
 
  // incoming beam e+
  GenParticle* p =
      new GenParticle( HepLorentzVector( 0, 0, 0.5 * cmE * 1e-3, 0.5 * cmE * 1e-3 ), -11, 3 );
  p->suggest_barcode( 1 );
  prod_vtx->add_particle_in( p );
 
  // incoming beam e-
  p = new GenParticle( HepLorentzVector( 0, 0, 0.5 * cmE * 1e-3, 0.5 * cmE * 1e-3 ), 11, 3 );
  p->suggest_barcode( 2 );
  prod_vtx->add_particle_in( p );
 
  int npart = 2;
  if ( proc < 100 )
  {
    double mom[4 * 6];
    int    np;
    gRad->MakeEvent( mom, np );
 
    // final particle 1
    for ( int i = 0; i < np; i++ )
    {
      double ptot = mom[i * 4 + 3];
      double px   = ptot * mom[i * 4 + 0];
      double py   = ptot * mom[i * 4 + 1];
      double pz   = ptot * mom[i * 4 + 2];
      double mass = 0;
      int    pid  = 22;
      if ( i < 2 )
      {
        pid  = fpid[i];
        mass = fM;
      }
      double etot = sqrt( ptot * ptot + mass * mass * 1e-6 );
      p           = new GenParticle( HepLorentzVector( px, py, pz, etot ), pid, 1 );
      p->suggest_barcode( i + 3 );
      prod_vtx->add_particle_out( p );
      npart++;
    }
  }
  else
  {
    int    ipart = CS->GenUnWeightedEvent();
    size_t nmax  = CS->GetNfinal() + ( ( ipart == 0 ) ? 1 : 2 );
    for ( size_t i = 0; i < nmax; i++ )
    {
      TLorentzVector& q   = *( CS->GetParticles()[i] );
      int             pid = CS->GetPid( i );
      p                   = new GenParticle(
          HepLorentzVector( q.X() * 1e-3, q.Y() * 1e-3, q.Z() * 1e-3, q.T() * 1e-3 ), pid, 1 );
      p->suggest_barcode( i + 3 );
      prod_vtx->add_particle_out( p );
      npart++;
    }
  }
 
  if ( log.level() < MSG::INFO ) { evt->print(); }
 
  // Check if the McCollection already exists
  SmartDataPtr<McGenEventCol> anMcCol( eventSvc(), "/Event/Gen" );
  if ( anMcCol != 0 )
  {
    // Add event to existing collection
    log << MSG::WARNING << "add event" << endmsg;
    MsgStream log( msgSvc(), name() );
    log << MSG::INFO << "Add McGenEvent to existing collection" << endmsg;
    McGenEvent* mcEvent = new McGenEvent( evt );
    anMcCol->push_back( mcEvent );
  }
  else
  {
    // Create Collection and add  to the transient store
    log << MSG::WARNING << "create collection" << endmsg;
    McGenEventCol* mcColl  = new McGenEventCol;
    McGenEvent*    mcEvent = new McGenEvent( evt );
    mcColl->push_back( mcEvent );
    StatusCode sc = eventSvc()->registerObject( "/Event/Gen", mcColl );
    if ( sc != StatusCode::SUCCESS )
    {
      log << MSG::ERROR << "Could not register McGenEvent" << endmsg;
      delete mcColl;
      delete evt;
      delete mcEvent;
      return StatusCode::FAILURE;
    }
    else
    {
      log << MSG::INFO << "McGenEventCol created and " << npart
          << " particles stored in McGenEvent" << endmsg;
    }
  }
 
  log << MSG::WARNING << "execute end" << endmsg;
  return StatusCode::SUCCESS;
}
 
StatusCode Mcgpj::finalize() {
  MsgStream log( msgSvc(), name() );
 
  delete gRad;
  delete MatrixElements;
  delete gRandom;
  delete gConst;
  delete gGlobal;
  delete gCut;
 
  log << MSG::INFO << "Mcgpj finalized" << endmsg;
 
  return StatusCode::SUCCESS;
}