Hydjet2Hadronizer.h

Go to the documentation of this file.

/*                                                                            
 Based on class InitalStateHydjet:                                          
 Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna
 amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru 
 November. 2, 2005                                

*/

#ifndef Hydjet2Hadronizer_h
#define Hydjet2Hadronizer_h

#include "DatabasePDG.h"
#include "Particle.h"
#include "InitialState.h"

#define PYCOMP pycomp_

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "GeneratorInterface/Core/interface/BaseHadronizer.h"

#include "CLHEP/Random/RandomEngine.h"
#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/RandPoisson.h"
#include "CLHEP/Random/RandGauss.h"

#include <map>
#include <string>
#include <vector>
#include <cmath>

#include "HYJET_COMMONS.h"
extern HYIPARCommon HYIPAR;
extern HYFPARCommon HYFPAR;
extern HYJPARCommon HYJPAR;
extern HYPARTCommon HYPART;
extern SERVICECommon SERVICE;

#define kMax 500000

namespace CLHEP {
  class HepRandomEngine;
  class RandFlat;
  class RandPoisson;
  class RandGauss;
}  // namespace CLHEP

CLHEP::HepRandomEngine* hjRandomEngine;

namespace HepMC {
  class GenEvent;
  class GenParticle;
  class GenVertex;
}  // namespace HepMC

namespace gen {
  class Pythia6Service;
  class Hydjet2Hadronizer : public InitialState, public BaseHadronizer {
  public:
    Hydjet2Hadronizer(const edm::ParameterSet&);
    ~Hydjet2Hadronizer() override;

    bool readSettings(int);
    bool declareSpecialSettings(const std::vector<std::string>&) { return true; }
    bool initializeForInternalPartons();
    bool initializeForExternalPartons();  //0
    bool generatePartonsAndHadronize();
    bool declareStableParticles(const std::vector<int>&);

    bool hadronize();  //0
    bool decay();      //0
    bool residualDecay();
    void finalizeEvent();
    void statistics();
    const char* classname() const;

    //________________________________________

    void SetVolEff(double value) { fVolEff = value; }
    double GetVolEff() { return fVolEff; }
    bool RunDecays() override { return (fDecay > 0 ? kTRUE : kFALSE); }
    double GetWeakDecayLimit() override { return fWeakDecay; }

    bool IniOfThFreezeoutParameters();

    double f(double);
    double f2(double, double, double);

    double SimpsonIntegrator(double, double, double, double);
    double SimpsonIntegrator2(double, double, double, double);
    double MidpointIntegrator2(double, double, double, double);
    double CharmEnhancementFactor(double, double, double, double);

  private:
    void doSetRandomEngine(CLHEP::HepRandomEngine* v) override;
    void rotateEvtPlane();
    bool get_particles(HepMC::GenEvent* evt);
    HepMC::GenParticle* build_hyjet2(int index, int barcode);
    HepMC::GenVertex* build_hyjet2_vertex(int i, int id);
    void add_heavy_ion_rec(HepMC::GenEvent* evt);

    std::vector<std::string> const& doSharedResources() const override { return theSharedResources; }
    static const std::vector<std::string> theSharedResources;

    inline double nuclear_radius() const;

    double fVolEff;  // the effective volume

    // the list of initial state parameters

    //int fNevnt;                      // number of events
    double fSqrtS;  // cms energy per nucleon
    double fAw;     // atomic number of colliding nuclei
    int fIfb;       // flag of type of centrality generation (=0 is fixed by fBfix, not 0
                    // impact parameter is generated in each event between fBfmin
                    // and fBmax according with Glauber model (f-la 30)
    double fBmin;   // minimum impact parameter in units of nuclear radius RA
    double fBmax;   // maximum impact parameter in units of nuclear radius RA
    double fBfix;   // fix impact parameter in units of nuclear radius RA

    double fT;          // chemical freeze-out temperature in GeV
    double fMuB;        // baryon potential
    double fMuS;        // strangeness potential
    double fMuC;        // charm potential
    double fMuI3;       // isospin potential
    double fThFO;       // thermal freeze-out temperature T^th in GeV
    double fMu_th_pip;  // effective chemical potential of positivly charged pions at thermal in GeV

    double fTau;          // proper time value
    double fSigmaTau;     // its standart deviation (emission duration)
    double fR;            // maximal transverse radius
    double fYlmax;        // maximal longitudinal rapidity
    double fUmax;         // maximal transverse velocity multiplaed on \gamma_r
    double fDelta;        // momentum asymmetry parameter
    double fEpsilon;      // coordinate asymmetry parameter
    int fIfDeltaEpsilon;  // flag to specify fDelta and fEpsilon values(=0 user's ones, >=1 parametrized)

    int fDecay;         // flag to switch on/off hadron decays (=0 decays off, >=1 decays on), (default: 1)
    double fWeakDecay;  // flag to switch on/off weak hadron decays <0: decays off, >0: decays on, (default: 0)
    int fPythDecay;     // Flag to choose how to decay resonances in high-pt part, fPythDecay: 0 by PYTHIA decayer,
                        // 1 by FASTMC decayer(mstj(21)=0)

    int fEtaType;  // flag to choose rapidity distribution, if fEtaType<=0,
                   // then uniform rapidity distribution in [-fYlmax,fYlmax] if fEtaType>0,
                   // then Gaussian with dispertion = fYlmax

    int fTMuType;  // flag to use calculated chemical freeze-out temperature,
                   // baryon potential and strangeness potential as a function of fSqrtS

    double fCorrS;   // flag and value to include strangeness supression factor
    int fCharmProd;  // flag to include statistical charm production
    double fCorrC;   // flag and value to include charmness supression factor

    int fNhsel;  // flag to switch on/off jet and hydro-state production (0: jet
                 // production off and hydro on, 1: jet production on and jet quenching
                 // off and hydro on, 2: jet production on and jet quenching on and
                 // hydro on, 3: jet production on and jet quenching off and hydro
                 // off, 4: jet production on and jet quenching on and hydro off
    int fPyhist;  // Suppress PYTHIA particle history (=1 only final state particles from hard part; =0 include full particle history)
    int fIshad;  // flag to switch on/off impact parameter dependent nuclear
                 // shadowing for gluons and light sea quarks (u,d,s) (0: shadowing off,
                 // 1: shadowing on for fAw=207, 197, 110, 40, default: 1

    double fPtmin;  // minimal transverse momentum transfer p_T of hard
                    // parton-parton scatterings in GeV (the PYTHIA parameter ckin(3)=fPtmin)

    //  PYQUEN energy loss model parameters:

    double fT0;  // initial temperature (in GeV) of QGP for
                 // central Pb+Pb collisions at mid-rapidity (initial temperature for other
                 // centralities and atomic numbers will be calculated automatically) (allowed range is 0.2<fT0<2)

    double fTau0;  // proper QGP formation time in fm/c (0.01<fTau0<10)
    int fNf;       // number of active quark flavours N_f in QGP fNf=0, 1,2 or 3
    int fIenglu;   // flag to fix type of in-medium partonic energy loss
                   // (0: radiative and collisional loss, 1: radiative loss only, 2:
                   // collisional loss only) (default: 0);
    int fIanglu;   // flag to fix type of angular distribution of in-medium emitted
                   // gluons (0: small-angular, 1: wide-angular, 2:collinear) (default: 0).

    bool embedding_;  // Switch for embedding mode
    bool rotate_;     // Switch to rotate event plane
    HepMC::GenEvent* evt;
    int nsub_;     // number of sub-events
    int nhard_;    // multiplicity of PYTHIA(+PYQUEN)-induced particles in event
    int nsoft_;    // multiplicity of HYDRO-induced particles in event
    double phi0_;  // Event plane angle
    double sinphi0_;
    double cosphi0_;
    Pythia6Service* pythia6Service_;

    unsigned int pythiaPylistVerbosity_;  // pythia verbosity; def=1
    unsigned int maxEventsToPrint_;       // Events to print if verbosity

    edm::InputTag src_;

    int fNPartTypes;         //counter of hadron species
    int fPartEnc[1000];      //Hadron encodings. Maximal number of hadron species is 1000!!!
    double fPartMult[2000];  //Multiplicities of hadron species
    double fPartMu[2000];    //Chemical potentials of hadron species
    double fMuTh[1000];      //Chemical potentials at thermal freezeout of hadron species

    //open charm
    double fNocth;
    double fNccth;

    edm::ParameterSet pset;
    edm::Service<TFileService> fs;

    int ev, sseed, Njet, Nbcol, Npart, Ntot, Npyt, Nhyd;
    double psiforv3;
    float Bgen, Sigin, Sigjet;
    float Px[kMax];
    float Py[kMax];
    float Pz[kMax];
    float E[kMax];
    float X[kMax];
    float Y[kMax];
    float Z[kMax];
    float T[kMax];
    int pdg[kMax];
    int Mpdg[kMax];
    int type[kMax];
    int pythiaStatus[kMax];
    int Index[kMax];
    int MotherIndex[kMax];
    int NDaughters[kMax];
    int FirstDaughterIndex[kMax];
    int LastDaughterIndex[kMax];
    int final[kMax];

    ParticleAllocator allocator;
    List_t source;
  };
  double Hydjet2Hadronizer::nuclear_radius() const {
    // Return the nuclear radius derived from the
    // beam/target atomic mass number.
    return 1.15 * pow((double)fAw, 1. / 3.);
  }
}  // namespace gen
#endif