BabayagaNLO.cxx

Go to the documentation of this file.

// #include "HepMC/HEPEVT_Wrapper.h"
#include "HepMC/IO_HEPEVT.h"
// #include "HepMC/IO_Ascii.h"
#include "HepMC/GenEvent.h"
 
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SmartDataPtr.h"
 
#include "BabayagaNLO.h"
#include "BabayagaNLORandom.h"
#include "BesRndmGenSvc/IBesRndmGenSvc.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "GeneratorObject/McGenEvent.h"
 
#include "boost/filesystem.hpp"
 
#define fs boost::filesystem
using namespace HepMC;
 
#ifdef DECLARE_COMPONENT
DECLARE_COMPONENT( BabayagaNLO )
#endif
 
extern "C" {
 
extern void bossinterface_( int* xpari, double* xpard );
extern void init_babayaga_();
extern void generate_event_( bool& use_unweighted );
extern void print_statistics_();
extern void set_ecms_spread_( double& ecms, double& bspread );
 
//    common/momentainitial/pin1,pin2
extern struct {
  double pin1[4];
  double pin2[4];
} momentainitial_;
 
//    common/event_mom/p1,p2,qph
extern struct {
  double p1[4], p2[4], qph[4][40];
} event_mom_;
 
// common/weights/use_weighted,sdif
extern struct {
  double sdif;
  bool   use_unweighted;
} weights_;
 
// common/babayagainitint/in_conf_spin,nphmax,naccepted,nwhenmax,nover,istopsearch,nneg,ng
extern struct {
  int  in_conf_spin, nphmax;
  long naccepted;
  int  nwhenmax, nover, istopsearch, nneg, ng;
} babayagainitint_;
 
// common/vpcnskf/nskfilepath
extern struct {
  char nskfilepath[400];
} vpcnskf_;
}
 
BabayagaNLO::BabayagaNLO( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ) {
  //    enum int_pars
  //{kChannel,kVerbose,kNsearch,kAlphaRun,kPhMode,kModelVP,kPrecision,kUseWeighted,kDarkMode};
  //    enum dbl_pars
  //{kEcm,kEspread,kCutTh0,kCutTh1,kCutAColl,kCutEmin,kCutInvMass0,kCutInvMass1,kDarkMass,kDarkWidth,kDarkVec,kDarkAx,kMaxSdif,kScaleVPerr};
 
  declareProperty( "Channel", m_ch = 0 ); // 0=ee, 1=mumu, 2 = gg
  declareProperty( "Verbose", m_iverbose = 0 );
  declareProperty( "Nsearch", m_nsearch = 4000000 ); // event used for maximum search
  declareProperty( "RunningAlpha", m_arun = 0 );     // alpha running
  declareProperty( "VPparam", m_iteubn = 0 ); // 0= Jegerlehner, 1 = Teubner, 2 = Novosibirsk
  declareProperty( "PhotonNumber",
                   m_photmode = -1 ); // maximum number of photon (-1 = unlimited (40))
  declareProperty( "Precision", m_precision = 0 ); // 0=best, 1=alpha2, 2=born
 
  declareProperty( "CMSEnergy", m_ecmsinput = 3.686 );
  declareProperty( "BeamEnergySpread", m_beamspread = 0.0013 );
 
  declareProperty( "ThetaMin", m_thmin = 10 );   //| Range of charged particles (ee,mumu)
  declareProperty( "ThetaMax", m_thmax = 170 );  //| or one photon pair (gg)
  declareProperty( "AcollMax", m_zmax = 180.0 ); // Maximum acollinearity
  declareProperty( "Emin",
                   m_emin = 0.001 ); // minimum energy of charged tracks/most energetic photons
  declareProperty( "MinInvMass",
                   m_Minv_min = 0.0 ); //| Range of invariant mass to be generated
  declareProperty( "MaxInvMass",
                   m_Minv_max = -1.0 ); //| charged (ee,mumu) or at least one photon pair (gg)
 
  declareProperty( "MaximumVal", m_sdif_max = 1.e-18 ); // starting value for maximum search
  declareProperty( "ScaleVPerr", m_scale_vperr = 0.0 ); // add hadronic vp error scaled by this
                                                        // value ( vp_hadr = vp_hadr +
                                                        // scale*d_hvp, default: no addition)
 
  declareProperty( "UseWeightedEv", m_weighted = 0 ); // enable generation of weighted events
                                                      // (default: unweighted)
 
  declareProperty( "DarkMode", m_darkmode = 0 );              // use dark photon (A') model
  declareProperty( "DarkMass", m_DarkMass = 0.4 );            // A' mass
  declareProperty( "DarkWidth", m_DarkWidth = -1.0 );         // A' width (default: point-like)
  declareProperty( "DarkVectorCoupling", m_DarkVec = 1.e-3 ); //|    A'
  declareProperty( "DarkAxialCoupling", m_DarkAxial = 0 );    //| couplings
 
  std::string path = "./";
  if ( getenv( "BABAYAGANLOROOT" ) != NULL )
  {
    path = getenv( "BABAYAGANLOROOT" );
    path += "/share/";
  }
  path.append( "vpol_novosibirsk.dat" );
  declareProperty( "NskFilePath", m_nsk_filepath = path ); // path to nsk VP input file
}
 
StatusCode BabayagaNLO::initialize() {
  MsgStream log( msgSvc(), name() );
  log << MSG::DEBUG << "BabayagaNLO in initialize()" << endmsg;
 
  // set Bes unified random engine
  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( "BabayagaNLO" );
  BabayagaNLORandom::setRandomEngine( engine );
 
  // start initialize
  int    xpari[10];
  double xpard[14];
 
  xpari[kChannel]     = m_ch;
  xpari[kVerbose]     = m_iverbose;
  xpari[kNsearch]     = m_nsearch;
  xpari[kAlphaRun]    = m_arun;
  xpari[kPhMode]      = m_photmode;
  xpari[kModelVP]     = m_iteubn;
  xpari[kPrecision]   = m_precision;
  xpari[kUseWeighted] = m_weighted;
  xpari[kDarkMode]    = m_darkmode;
 
  xpard[kEcm]         = m_ecmsinput;
  xpard[kEspread]     = m_beamspread;
  xpard[kCutTh0]      = m_thmin;
  xpard[kCutTh1]      = m_thmax;
  xpard[kCutAColl]    = m_zmax;
  xpard[kCutEmin]     = m_emin;
  xpard[kCutInvMass0] = m_Minv_min;
  xpard[kCutInvMass1] = m_Minv_max;
  xpard[kDarkMass]    = m_DarkMass;
  xpard[kDarkWidth]   = m_DarkWidth;
  xpard[kDarkVec]     = m_DarkVec;
  xpard[kDarkAx]      = m_DarkAxial;
  xpard[kMaxSdif]     = m_sdif_max;
  xpard[kScaleVPerr]  = m_scale_vperr;
 
  using namespace boost::filesystem;
  path p_nskfile( m_nsk_filepath.data() );
  if ( !exists( p_nskfile ) && m_arun > 0 && m_iteubn == 2 )
  {
    log << MSG::FATAL << "Cannot find Novosibirsk VP file (" << p_nskfile.native()
        << ")!\n Is the environment variable BABAYAGANLOROOT set?" << endmsg;
    return StatusCode::FAILURE;
  }
 
  m_nsk_filepath = p_nskfile.native();
  strcpy( vpcnskf_.nskfilepath, m_nsk_filepath.data() );
 
  bossinterface_( xpari, xpard );
 
  init_babayaga_();
 
  log << MSG::DEBUG << "Finish BabayagaNLO initialize()" << endmsg;
  return StatusCode::SUCCESS;
}
 
StatusCode BabayagaNLO::execute() {
  MsgStream log( msgSvc(), name() );
  log << MSG::DEBUG << "BabayagaNLO in execute()" << endmsg;
 
  // Fill event information
  GenEvent*  evt      = new GenEvent( 1, 1 );
  GenVertex* prod_vtx = new GenVertex();
  evt->add_vertex( prod_vtx );
 
  log << MSG::DEBUG << "check point 1" << endmsg;
  bool unw = !( m_weighted );
 
  log << MSG::DEBUG << "check point 2" << endmsg;
  generate_event_( unw );
 
  log << MSG::DEBUG << "check point 3" << endmsg;
  int finalpidm, finalpidp;
  if ( m_ch == 0 )
  {
    finalpidm = 11;
    finalpidp = -11;
  }
  else if ( m_ch == 1 )
  {
    finalpidm = 13;
    finalpidp = -13;
  }
  else if ( m_ch == 2 )
  {
    finalpidm = 22;
    finalpidp = 22;
  }
  else
  {
    finalpidm = 11;
    finalpidp = -11;
  }
 
  double* p1 = momentainitial_.pin1;
  double* p2 = momentainitial_.pin2;
  double* p3 = event_mom_.p1;
  double* p4 = event_mom_.p2;
 
  // incoming particle -
  GenParticle* p =
      new GenParticle( CLHEP::HepLorentzVector( p1[1], p1[2], -p1[3], p1[0] ), 11, 3 );
  p->suggest_barcode( 1 );
  prod_vtx->add_particle_in( p );
 
  // incoming particle +
  p = new GenParticle( CLHEP::HepLorentzVector( p2[1], p2[2], -p2[3], p2[0] ), -11, 3 );
  p->suggest_barcode( 2 );
  prod_vtx->add_particle_in( p );
 
  log << MSG::DEBUG << "check point 4" << endmsg;
  // outgoing particle -
  p = new GenParticle( CLHEP::HepLorentzVector( p3[1], p3[2], -p3[3], p3[0] ), finalpidm, 1 );
  p->suggest_barcode( 3 );
  prod_vtx->add_particle_out( p );
 
  // outgoing particle +
  p = new GenParticle( CLHEP::HepLorentzVector( p4[1], p4[2], -p4[3], p4[0] ), finalpidp, 1 );
  p->suggest_barcode( 4 );
  prod_vtx->add_particle_out( p );
 
  log << MSG::DEBUG << "check point 5" << endmsg;
  // photons
  int id_cntr = 5;
  for ( int i = 0; i < babayagainitint_.ng; i++ )
  {
    double px  = event_mom_.qph[1][i];
    double py  = event_mom_.qph[2][i];
    double pz  = -1. * event_mom_.qph[3][i];
    double eph = event_mom_.qph[0][i];
    p          = new GenParticle( CLHEP::HepLorentzVector( px, py, pz, eph ), 22, 1 );
    p->suggest_barcode( id_cntr );
    prod_vtx->add_particle_out( p );
    id_cntr++;
  }
 
  log << MSG::DEBUG << "check point 6" << endmsg;
  if ( !unw )
  { // Add dummy particle with diffential cross section as mass
    p = new GenParticle( CLHEP::HepLorentzVector( 0, 0, 0, weights_.sdif ), 99, 1 );
    p->suggest_barcode( id_cntr );
    prod_vtx->add_particle_out( p );
    id_cntr++;
  }
 
  if ( log.level() <= MSG::DEBUG ) { 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!" << endmsg; }
  }
 
  log << MSG::DEBUG << "before execute() return" << endmsg;
  return StatusCode::SUCCESS;
}
 
StatusCode BabayagaNLO::finalize() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "BabayagaNLO in finalize()" << endmsg;
 
  print_statistics_();
  return StatusCode::SUCCESS;
}