SQPrimaryParticleGen.h

Go to the documentation of this file.

/*====================================================================
Author: Abinash Pun, Kun Liu
Nov, 2019
Goal: Import the physics generator of E906 experiment (DPSimPrimaryGeneratorAction)
from Kun to E1039 experiment in Fun4All framework
=========================================================================*/
#ifndef __SQPrimaryParticleGen_H__
#define __SQPrimaryParticleGen_H__
 
#ifndef __CINT__
#include <Pythia8/Pythia.h>
#endif
#include <LHAPDF/LHAPDF.h>
#ifndef __CINT__
#include <gsl/gsl_rng.h>
#endif
 
#include <string>
 
#include <TGenPhaseSpace.h>
#include <g4main/PHG4ParticleGeneratorBase.h>
 
class PHCompositeNode;
class PHG4ParticleGeneratorBase;
class PHG4InEvent;
class PHG4Particle;
 
class PHGenIntegral;
class SQEvent;
class SQMCEvent;
class SQDimuon;
class SQDimuonVector;
class SQPrimaryVertexGen;
 
//==========
 
class SQPrimaryParticleGen: public PHG4ParticleGeneratorBase
{
public:
  typedef enum {
    Undef    =  0,
    DrellYan = -1,
    JPsi     = -2,
    PsiPrime = -3,
    Pythia   = -4,
    Custom   = -5
  } ProcessCode_t; 
  
    SQPrimaryParticleGen(const std::string& name = "PrimaryGen");
    virtual ~SQPrimaryParticleGen();
 
    
    int Init(PHCompositeNode* topNode);
    int InitRun(PHCompositeNode* topNode);
    int End(PHCompositeNode* topNode);
    int process_event(PHCompositeNode* topNode);
    //void GeneratePrimaries(G4Event* anEvent);
 
 
    int generateDrellYan(const TVector3& vtx, const double pARatio, double luminosity);
    int generateJPsi(const TVector3& vtx, const double pARatio, double luminosity);
    int generatePsip(const TVector3& vtx, const double pARatio, double luminosity);
    // void generateDarkPhotonFromEta();
    int generatePythia(const TVector3& vtx, const double pARatio);
    //int generateCustomDimuon(PHCompositeNode *topNode, TVector3 vtx, const double pARatio);
  
 
    bool generateDimuon(double mass, double xF);
 
    //swith for the generators; Abi
    //@
    void enablePythia      () { _proc_code = Pythia  ;}
    void enableCustomDimuon() { _proc_code = Custom  ;}
    void enableDrellYanGen () { _proc_code = DrellYan;}
    void enableJPsiGen     () { _proc_code = JPsi    ;}
    void enablePsipGen     () { _proc_code = PsiPrime;}
 
    void set_pdfset(const std::string name) { _pdfset = name; }
 
    void set_pT0DY    (const double val) { _pT0DY     = val; }
    void set_pTpowDY  (const double val) { _pTpowDY   = val; }
    void set_pT0JPsi  (const double val) { _pT0JPsi   = val; }
    void set_pTpowJPsi(const double val) { _pTpowJPsi = val; }
 
    void set_config_file(const char *cfg_file)
    {
     if (cfg_file) _configFile = cfg_file; 
    }
 
    void set_xfRange(const double xmin, const double xmax){
      xfMin = xmin;
      xfMax = xmax;
    }
    void set_massRange(const double mmin, const double mmax){
      massMin = mmin;
      massMax = mmax;
    }
    //@
 
    double CrossSectionDrellYan(const double mass, const double xF, const double pARatio);
    double CrossSectionDrellYan(const double mass, const double xF, const double x1, const double x2, const double pARatio);
    double CrossSectionJPsi(const double xF);
    double CrossSectionPsip(const double xF);
 
 private:
    ProcessCode_t _proc_code;
 
    SQPrimaryVertexGen* _vertexGen;
    PHG4InEvent *ineve;
    SQEvent* _evt; //< An output node
    SQMCEvent* _mcevt; //< An output node
    SQDimuonVector* _vec_dim; //< An output node
    PHGenIntegral *_integral_node; //< An output node
    
    SQDimuon* _dim_gen; //< To hold the kinematics of a dimuon generated
 
    //Pythia generator
    Pythia8::Pythia ppGen;    
    Pythia8::Pythia pnGen;    
    Pythia8::Pythia _Pythia;
    std::string _configFile;
    int read_config(const char *cfg_file = 0);
 
    TGenPhaseSpace phaseGen;
 
    std::string  _pdfset;
    LHAPDF::PDF* _pdf;
 
    // Parameters (being moved from DPGEN)
    double _pT0DY;
    double _pTpowDY;
    double _pT0JPsi;
    double _pTpowJPsi;
  
    //some initializations
    double massMin = 0.22;
    double massMax = 10.;
    double x1Min = 0.;
    double x1Max = 1.;
    double x2Min = 0.;
    double x2Max = 1.;
    double xfMin = -1.;
    double xfMax = 1.;
    double cosThetaMin = -1.;
    double cosThetaMax = 1. ;
    double zOffsetMin = -1.;
    double zOffsetMax = 1.;
 
    void InsertMuonPair(const TVector3& vtx);
    void InsertEventInfo(const double xsec, const double weight, const TVector3& vtx);
 
    double _n_gen_acc_evt; //< N of generator-accepted events
    double _n_proc_evt; //< N of processed events
    double _weight_sum; //< Sum of weights
    double _inte_lumi; //< Integrated luminosity
};
 
//========
 
 
 
 
#endif