Babayaga.cxx

Go to the documentation of this file.

//*****************************************************************************
//
// BABAYAGA.cxx
//
// Algorithm runs BHABHA precess
//
// July 2007 Rong-Gang Ping to migrate BABAYAGA3.5 for BES3
// The accuracy is about 0.1% (see Nucl. Phys. B 758(2006) 227-253)
//*****************************************************************************
 
// Header for this module:-
#include "Babayaga.h"
#include "BabayagaRandom.h"
#include "BesRndmGenSvc/IBesRndmGenSvc.h"
 
// Framework Related Headers:-
#include "HepMC/GenEvent.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
 
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SmartDataPtr.h"
 
#include "CLHEP/Vector/LorentzVector.h"
#include "GeneratorObject/McGenEvent.h"
#include "cfortran/cfortran.h"
 
#include <stdlib.h>
 
// COMMON/MOMSET/PM,PP,PFINEL,PFINPOS,QPH
// dimension PM(0:3),PP(0:3),PFINEL(0:3),PFINPOS(0:3),QPH(40,0:3)  //OUTPUT RESULTS
// PM(0:3),PP(0:3) for beam electron and positron particle
 
typedef struct {
  double q1[4][160001];
  double p1[4][160001];
  double q2[4][160001];
  double p2[4][160001];
  double phot[4][10][160001];
} MOMSET_DEF;
// MOMSET_DEF MOMSET;
#define MOMSET COMMON_BLOCK( MOMSET_DEF, momset )
COMMON_BLOCK_DEF( MOMSET_DEF, MOMSET );
typedef struct {
  double CutNgam, CutEgam;
} RADBB_DEF;
// RADBB_DEF RADBB;
#define RADBB COMMON_BLOCK( RADBB_DEF, radbb )
COMMON_BLOCK_DEF( RADBB_DEF, RADBB );
 
// COMMON/ISRPHOTONS/NCQPH
typedef struct {
  int ncqph[160000];
} ISRPHOTONS_DEF;
#define ISRPHOTONS COMMON_BLOCK( ISRPHOTONS_DEF, isrphotons )
COMMON_BLOCK_DEF( ISRPHOTONS_DEF, ISRPHOTONS );
 
// common/beamEnergy/ebeam
typedef struct {
  double ebeam;
} BEAMENERGY_DEF;
#define BEAMENERGY COMMON_BLOCK( BEAMENERGY_DEF, beamenergy )
COMMON_BLOCK_DEF( BEAMENERGY_DEF, BEAMENERGY );
// BEAMENERGY_DEF,BEAMENERGY
 
// common/expcuts/thmin,thmax,emin,zmax,egmin,thgmin,thgmax
typedef struct {
  double thmin, thmax, emin, zmax, egmin, thgmin, thgmax;
} EXPCUTS_DEF;
#define EXPCUTS COMMON_BLOCK( EXPCUTS_DEF, expcuts )
COMMON_BLOCK_DEF( EXPCUTS_DEF, EXPCUTS );
// EXPCUTS_DEF,EXPCUTS;
 
// common/switchFSR/ON,ILL
typedef struct {
  int on, ill;
} SWITCH_DEF;
#define SWITCH COMMON_BLOCK( SWITCH_DEF, switchfsr )
COMMON_BLOCK_DEF( SWITCH_DEF, SWITCH );
 
// common/switcharun/IARUN
typedef struct {
  int iarun;
} SWITCHARUN_DEF;
#define SWITCHARUN COMMON_BLOCK( SWITCHARUN_DEF, switcharun )
COMMON_BLOCK_DEF( SWITCHARUN_DEF, SWITCHARUN );
 
// COMMON/CHANNEL/ICH
typedef struct {
  int ich;
} CHANNEL_DEF;
#define CHANNEL COMMON_BLOCK( CHANNEL_DEF, channel )
COMMON_BLOCK_DEF( CHANNEL_DEF, CHANNEL );
 
// common/resonances/IRES
typedef struct {
  int ires;
} RESONANCES_DEF;
#define RESONANCES COMMON_BLOCK( RESONANCES_DEF, resonances )
COMMON_BLOCK_DEF( RESONANCES_DEF, RESONANCES );
 
// COMMON/DECLAREINT/INT
typedef struct {
  int seed;
} DECLAREINT_DEF;
#define DECLAREINT COMMON_BLOCK( DECLAREINT_DEF, declareint )
COMMON_BLOCK_DEF( DECLAREINT_DEF, DECLAREINT );
 
// COMMON/DECLARESTR/INTUPLE,PHCUT,CUTG
typedef struct {
  char tuple, phcut, cutg;
} DECLARESTR_DEF;
#define DECLARESTR COMMON_BLOCK( DECLARESTR_DEF, declarestr )
COMMON_BLOCK_DEF( DECLARESTR_DEF, DECLARESTR );
 
//--//Unify Babayaga random engine with Bes random servive
extern "C" {
extern float flat_();
}
 
float flat_() {
  float  rr;
  double rd = BabayagaRandom::random();
  //  std::cout<<"BabayagaRandom::random()="<<rd<<endl;
  rr = (float)rd;
  return (float)BabayagaRandom::random();
  //  return rr;
}
 
PROTOCCALLSFSUB1( BABAYAGA, babayaga, INT )
#define BABAYAGA( NEVENTS ) CCALLSFSUB1( BABAYAGA, babayaga, INT, NEVENTS )
 
DECLARE_COMPONENT( Babayaga )
Babayaga::Babayaga( const string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ) {
  //  declareProperty("Seed", m_Int = 62341342); // seed for random number generator
  declareProperty( "Channel", m_Ich = 1 );     // 1: e+e-->e+e-;2:e+e_->mu+mu-;3:e+e-->gamma
                                               // gamma;4:e+e--->pi+pi-
  declareProperty( "Ebeam", m_Ebeam = 1.548 ); // Ecm = 2*Ebeam [GeV]
  declareProperty( "MinThetaAngle", m_Thmin = 70 );       // [degree]
  declareProperty( "MaxThetaAngle", m_Thmax = 178 );      // [degree]
  declareProperty( "MinimumEnergy", m_Emin = 0.4 );       // Minimum Energy(GeV)
  declareProperty( "MaximumAcollinearity", m_Zmax = 10 ); // Maximum acollinearity (degree)
  declareProperty( "RunningAlpha", m_Iarun = 1 );         // running alpha, 0=off, 1=on
  declareProperty( "HadronicResonance", m_Ires = 0 );     // hadronic resoances for ICH=1 or 2
  declareProperty( "FSR_swich", m_on = 0 );         // FSR switch for ICH=2  ( 0=off, 1=on)
  declareProperty( "MinEnerCutG", m_Egmin = 0.01 ); // minimum energy for CUTG=Y (GeV)
  declareProperty( "MinAngCutG", m_Thgmin = 5 );    // minimum angle for cuts =Y (deg.)
  declareProperty( "MaxAngCutG", m_Thgmax = 21 );   // maximum angle for CUTG=Y (deg.)
  declareProperty( "HBOOK", HN = 1 );      // INTUPLE:  if events have to be stored (1/0)
  declareProperty( "PHCUT", m_PHCUT = 0 ); // PHCUT:    to avoid double counting when ICH = 3
                                           // (1/0), for other channels, it set as 0 by default
  declareProperty( "CUTG",
                   m_CUTG = 0 ); // CUTG:    explicit cuts on the generated photons (1/0)
  declareProperty( "CutNgam", m_CutNgam = 0 ); // mininum number of ISR photon cut on Radiative
                                               // Bhabha events
  declareProperty( "CutEgam", m_CutEgam = 0 ); // in GeV, mininum energy of ISR photon cut on
                                               // Radiative Bhabha events
}
 
StatusCode Babayaga::initialize() {
 
  MsgStream log( msgSvc(), name() );
 
  log << MSG::INFO << "Babayaga 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( "Babayaga" );
  std::cout << "===============================" << engine << endl;
  BabayagaRandom::setRandomEngine( engine );
  //  *****************
 
  if ( HN == 1 ) { DECLARESTR.tuple = 'Y'; }
  else DECLARESTR.tuple = 'N';
  if ( m_PHCUT == 1 ) { DECLARESTR.phcut = 'Y'; }
  else DECLARESTR.cutg = 'N';
  if ( m_CUTG == 1 ) { DECLARESTR.cutg = 'Y'; }
  else DECLARESTR.cutg = 'N';
  CHANNEL.ich      = m_Ich;
  BEAMENERGY.ebeam = m_Ebeam;
  EXPCUTS.thmin    = m_Thmin;
  EXPCUTS.thmax    = m_Thmax;
  EXPCUTS.emin     = m_Emin;
  EXPCUTS.zmax     = m_Zmax;
  SWITCHARUN.iarun = m_Iarun;
  RESONANCES.ires  = m_Ires;
  SWITCH.on        = m_on;
  EXPCUTS.egmin    = m_Egmin;
  EXPCUTS.thgmin   = m_Thgmin;
  EXPCUTS.thgmax   = m_Thgmax;
 
  //  std::cout<<"m_Ires= "<<m_Ires<<endl;
 
  getMaxEvent().ignore();
  std::cout << "m_evtMax = " << m_evtMax << std::endl;
  DECLAREINT.seed = m_Int;
  BABAYAGA( m_evtMax );
 
  return StatusCode::SUCCESS;
}
 
static int evtgen = 1;
StatusCode Babayaga::execute() {
  MsgStream log( msgSvc(), name() );
  // log << MSG::INFO << "BABAYAGA executing" << endmsg;
  //  std::cout<<"BABAYAGA begin executing"<<std::endl;
 
  int npart = 0;
  int pid1, pid2, pst1, pst2;
  if ( m_Ich == 1 )
  {
    pid1 = 11;
    pid2 = -11;
    pst1 = 1;
    pst2 = 1;
  }
  if ( m_Ich == 2 )
  {
    pid1 = 13;
    pid2 = -13;
    pst1 = 1;
    pst2 = 1;
  }
  if ( m_Ich == 3 )
  {
    pid1 = 22;
    pid2 = 22;
    pst1 = 1;
    pst2 = 1;
  }
  if ( m_Ich == 4 )
  {
    pid1 = -211;
    pid2 = 211;
    pst1 = 1;
    pst2 = 1;
  }
 
  // Fill event information
  GenEvent* evt = new GenEvent( 1, 1 );
 
  GenVertex* prod_vtx = new GenVertex();
  //  prod_vtx->add_particle_out( p );
  evt->add_vertex( prod_vtx );
 
  // incoming beam e- HepLorentzVector(px,py,pz,energy) required!
  GenParticle* p =
      new GenParticle( CLHEP::HepLorentzVector( MOMSET.p1[1][evtgen], MOMSET.p1[2][evtgen],
                                                MOMSET.p1[3][evtgen], MOMSET.p1[0][evtgen] ),
                       11, 3 );
  p->suggest_barcode( ++npart );
  prod_vtx->add_particle_in( p );
 
  //  std::cout<<"incoming beam e+"<<endl;
  // incoming beam e+, e+ moving along the z-direction
  p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q1[1][evtgen], MOMSET.q1[2][evtgen],
                                                MOMSET.q1[3][evtgen], MOMSET.q1[4][evtgen] ),
                       -11, 3 );
 
  p->suggest_barcode( ++npart );
  prod_vtx->add_particle_in( p );
 
  // scattered lepton +
  p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.p2[1][evtgen], MOMSET.p2[2][evtgen],
                                                MOMSET.p2[3][evtgen], MOMSET.p2[0][evtgen] ),
                       pid1, pst1 );
  p->suggest_barcode( ++npart );
  prod_vtx->add_particle_out( p );
 
  // debuging pingrg
  //  std::cout<<"evtgen= "<<evtgen<<std::endl;
  // std::cout<<"b:e-: "<<MOMSET.p1[0][evtgen]<<"; "<< MOMSET.p1[1][evtgen]<<";
  // "<<MOMSET.p1[2][evtgen]<<"; "<< MOMSET.p1[3][evtgen]<<std::endl; std::cout<<"b:e+:
  // "<<MOMSET.q1[0][evtgen]<<"; "<< MOMSET.q1[1][evtgen]<<"; "<<MOMSET.q1[2][evtgen]<<"; "<<
  // MOMSET.q1[3][evtgen]<<std::endl; std::cout<<"e-: "<<MOMSET.p2[0][evtgen]<<"; "<<
  // MOMSET.p2[1][evtgen]<<"; "<<MOMSET.p2[2][evtgen]<<"; "<< MOMSET.p2[3][evtgen]<<std::endl;
  // std::cout<<"e+: "<<MOMSET.q2[0][evtgen]<<"; "<< MOMSET.q2[1][evtgen]<<";
  // "<<MOMSET.q2[2][evtgen]<<"; "<< MOMSET.q2[3][evtgen]<<std::endl;
 
  // scattered lepton -
  p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q2[1][evtgen], MOMSET.q2[2][evtgen],
                                                MOMSET.q2[3][evtgen], MOMSET.q2[0][evtgen] ),
                       pid2, pst2 );
  p->suggest_barcode( ++npart );
  prod_vtx->add_particle_out( p );
 
  int iphot = 0;
  // std::cout<<"evtgen, ncqph[events#] "<<evtgen<<"; "<<ISRPHOTONS.ncqph[evtgen-1]<<std::endl;
  for ( iphot = 0; iphot < ISRPHOTONS.ncqph[evtgen - 1]; iphot++ )
  {
    // debuging, pingrg
    //    std::cout<<"evtgen, iphot= "<<evtgen<<"; "<<iphot<<std::endl;
    // std::cout<<MOMSET.phot[0][iphot][evtgen]<<", "<<MOMSET.phot[1][iphot][evtgen]<<",
    // "<<MOMSET.phot[2][iphot][evtgen]<<", "<<MOMSET.phot[3][iphot][evtgen]<<std::endl;
 
    p = new GenParticle( CLHEP::HepLorentzVector(
                             MOMSET.phot[1][iphot][evtgen], MOMSET.phot[2][iphot][evtgen],
                             MOMSET.phot[3][iphot][evtgen], MOMSET.phot[0][iphot][evtgen] ),
                         22, 1 );
    p->suggest_barcode( ++npart );
    prod_vtx->add_particle_out( p );
  }
 
  evtgen++;
 
  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
    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
    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;
    }
  }
 
  return StatusCode::SUCCESS;
}
 
StatusCode Babayaga::finalize() {
  MsgStream log( msgSvc(), name() );
  char      delcmd[300];
  strcpy( delcmd, "cat " );
  strcat( delcmd, "CrossSection.txt" );
  system( delcmd );
 
  std::cout << "BABAYAGA finalized" << endl;
  return StatusCode::SUCCESS;
}
 
StatusCode Babayaga::getMaxEvent() {
  IProperty* appPropMgr = 0;
  StatusCode status =
      serviceLocator()->getService( "ApplicationMgr", IProperty::interfaceID(),
                                    reinterpret_cast<IInterface*&>( appPropMgr ) );
  if ( status.isFailure() ) return status;
 
  IntegerProperty evtMax( "EvtMax", 0 );
  status = appPropMgr->getProperty( &evtMax );
  if ( status.isFailure() ) return status;
 
  m_evtMax = evtMax.value();
  return status;
}