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/MsgStream.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/DataSvc.h"
#include "GaudiKernel/SmartDataPtr.h"
 
#include "BabayagaNLO/BabayagaNLO.h"
#include "BabayagaNLO/BabayagaNLORandom.h"
#include "GeneratorObject/McGenEvent.h"
#include "BesKernel/IBesRndmGenSvc.h"
#include "CLHEP/Vector/LorentzVector.h"
 
#include <iostream>
#include "boost/filesystem/path.hpp"
 
#define fs boost::filesystem
using std::cout;
using std::endl;
using namespace HepMC;
 
 
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_;
 
}
 
 
DECLARE_COMPONENT(BabayagaNLO)
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 = getenv("BABAYAGANLOROOT");
    path.append("/share/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()" << endreq;
 
    //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" << endreq;
        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 file Novosibirsk VP file"<<endreq;
    //  return StatusCode::ERROR;
    //}
 
    m_nsk_filepath = p_nskfile.native();
    strcpy(vpcnskf_.nskfilepath,m_nsk_filepath.data());
 
    bossinterface_(xpari,xpard);
 
    init_babayaga_();
 
    // if (!m_weighted) {
    //  log <<MSG::INFO<< "Looking for the maximum" <<endreq;
    //  bool gen_unw = 0;
    //  while  (babayagainitint_.istopsearch==0){
    //      generate_event_(gen_unw);
    //  }
    //  log <<MSG::INFO<< "Search for maximum ended"<<endreq;
    // }
 
    log <<MSG::DEBUG<< "Finish BabayagaNLO initialize()" <<endreq;
    return StatusCode::SUCCESS;
}
 
StatusCode BabayagaNLO::execute() {
    MsgStream log(msgSvc(), name());
    log << MSG::DEBUG << "BabayagaNLO in execute()" << endreq;
 
    // 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" << endreq;
    bool unw = !(m_weighted);
 
    log << MSG::DEBUG << "check point 2" << endreq;
    generate_event_(unw);
 
    log << MSG::DEBUG << "check point 3" << endreq;
    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, 1);
    p->suggest_barcode(3);
    prod_vtx->add_particle_in(p);
 
    //incoming particle +
    p = new GenParticle(CLHEP::HepLorentzVector(p2[1],p2[2],-p2[3],p2[0]), -11, 1);
    p->suggest_barcode(4);
    prod_vtx->add_particle_in(p);
 
    log << MSG::DEBUG << "check point 4" << endreq;
    //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" << endreq;
    //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 = 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" << endreq;
    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"<<endreq;
        MsgStream log(messageService(), name());
        log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
        McGenEvent* mcEvent = new McGenEvent(evt);
        anMcCol->push_back(mcEvent);
    } else {
        // Create Collection and add  to the transient store
        log<<MSG::WARNING<<"create collection"<<endreq;
        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" << endreq;
            delete mcColl;
            delete evt;
            delete mcEvent;
            return StatusCode::FAILURE;
        } else {
            log << MSG::INFO << "McGenEventCol created!" << endreq;
        }
    }
 
 
    log<<MSG::DEBUG<< "before execute() return"<<endreq;
    return StatusCode::SUCCESS;
 
}
 
StatusCode BabayagaNLO::finalize() {
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "BabayagaNLO in finalize()" << endreq;
 
    print_statistics_();
    return StatusCode::SUCCESS;
}