HijingHadronizer.cc

Go to the documentation of this file.

#include <iostream>
#include <cmath>

#include "boost/lexical_cast.hpp"

#include "FWCore/Concurrency/interface/SharedResourceNames.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/Run.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/EDMException.h"

#include "GeneratorInterface/HijingInterface/interface/HijingHadronizer.h"
#include "GeneratorInterface/HijingInterface/interface/HijingPythiaWrapper.h"
#include "GeneratorInterface/HijingInterface/interface/HijingWrapper.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"

#include "HepMC/GenEvent.h"
#include "HepMC/HeavyIon.h"
#include "HepMC/SimpleVector.h"

#include "CLHEP/Random/RandomEngine.h"

static const double pi = 3.14159265358979;

using namespace edm;
using namespace std;
using namespace gen;

static CLHEP::HepRandomEngine* hijRandomEngine;

extern "C" {
float gen::hijran_(int* idummy) { return hijRandomEngine->flat(); }
}

extern "C" {
float ran_(unsigned int* iseed) {
  return hijRandomEngine->flat();
  //      return ranff_(iseed);
  //      return gen::pyr_(0);
}
}

extern "C" {
float rlu_(unsigned int* iseed) {
  return hijRandomEngine->flat();
  //      return ranff_(iseed);
  //      return gen::pyr_(0);
}
}

const std::vector<std::string> HijingHadronizer::theSharedResources = {edm::SharedResourceNames::kPythia6};

HijingHadronizer::HijingHadronizer(const ParameterSet& pset)
    : BaseHadronizer(pset),
      evt(nullptr),
      pset_(pset),
      bmax_(pset.getParameter<double>("bMax")),
      bmin_(pset.getParameter<double>("bMin")),
      efrm_(pset.getParameter<double>("comEnergy")),
      frame_(pset.getParameter<string>("frame")),
      proj_(pset.getParameter<string>("proj")),
      targ_(pset.getParameter<string>("targ")),
      iap_(pset.getParameter<int>("iap")),
      izp_(pset.getParameter<int>("izp")),
      iat_(pset.getParameter<int>("iat")),
      izt_(pset.getParameter<int>("izt")),
      phi0_(0.),
      sinphi0_(0.),
      cosphi0_(1.),
      rotate_(pset.getParameter<bool>("rotateEventPlane")) {
  // Default constructor
}

//_____________________________________________________________________
HijingHadronizer::~HijingHadronizer() {
  // destructor
}

//_____________________________________________________________________
void HijingHadronizer::doSetRandomEngine(CLHEP::HepRandomEngine* v) { hijRandomEngine = v; }

//_____________________________________________________________________
void HijingHadronizer::add_heavy_ion_rec(HepMC::GenEvent* evt) {
  // heavy ion record in the final CMSSW Event
  HepMC::HeavyIon* hi = new HepMC::HeavyIon(himain1.jatt,  // Ncoll_hard/N of SubEvents
                                            himain1.np,    // Npart_proj
                                            himain1.nt,    // Npart_targ
                                            himain1.n0 + himain1.n01 + himain1.n10 + himain1.n11,  // Ncoll
                                            0,                                                     // spectator_neutrons
                                            0,                                                     // spectator_protons
                                            himain1.n01,  // N_Nwounded_collisions
                                            himain1.n10,  // Nwounded_N_collisions
                                            himain1.n11,  // Nwounded_Nwounded_collisions
                                            //gsfs Changed from 19 to 18 (Fortran counts from 1 , not 0)
                                            hiparnt.hint1[18],  // impact_parameter in [fm]
                                            phi0_,              // event_plane_angle
                                            0,                  // eccentricity
                                            //gsfs Changed from 12 to 11 (Fortran counts from 1 , not 0)
                                            hiparnt.hint1[11]  // sigma_inel_NN
  );
  evt->set_heavy_ion(*hi);
  delete hi;
}

//___________________________________________________________________
HepMC::GenParticle* HijingHadronizer::build_hijing(int index, int barcode) {
  // Build particle object corresponding to index in hijing

  double x0 = himain2.patt[0][index];
  double y0 = himain2.patt[1][index];

  double x = x0 * cosphi0_ - y0 * sinphi0_;
  double y = y0 * cosphi0_ + x0 * sinphi0_;

  // Hijing gives V0's status=4, they need to have status=1 to be decayed in geant
  // also change status=11 to status=2
  if (himain2.katt[3][index] <= 10 && himain2.katt[3][index] > 0)
    himain2.katt[3][index] = 1;
  if (himain2.katt[3][index] <= 20 && himain2.katt[3][index] > 10)
    himain2.katt[3][index] = 2;

  HepMC::GenParticle* p = new HepMC::GenParticle(HepMC::FourVector(x,                        // px
                                                                   y,                        // py
                                                                   himain2.patt[2][index],   // pz
                                                                   himain2.patt[3][index]),  // E
                                                 himain2.katt[0][index],                     // id
                                                 himain2.katt[3][index]                      // status
  );
  p->suggest_barcode(barcode);

  return p;
}

//___________________________________________________________________
HepMC::GenVertex* HijingHadronizer::build_hijing_vertex(int i, int id) {
  // build verteces for the hijing stored events
  double x0 = himain2.vatt[0][i];
  double y0 = himain2.vatt[1][i];
  double x = x0 * cosphi0_ - y0 * sinphi0_;
  double y = y0 * cosphi0_ + x0 * sinphi0_;
  double z = himain2.vatt[2][i];
  double t = himain2.vatt[3][i];

  HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(x, y, z, t), id);
  return vertex;
}

bool HijingHadronizer::generatePartonsAndHadronize() {
  // generate single event
  if (rotate_)
    rotateEvtPlane();

  // generate a HIJING event

  float f_bmin = bmin_;
  float f_bmax = bmax_;
  HIJING(frame_.data(), f_bmin, f_bmax, strlen(frame_.data()));

  // event information
  HepMC::GenEvent* evt = new HepMC::GenEvent();
  get_particles(evt);

  //   evt->set_signal_process_id(pypars.msti[0]);      // type of the process
  //   evt->set_event_scale(pypars.pari[16]);           // Q^2
  add_heavy_ion_rec(evt);

  event().reset(evt);

  return true;
}

//_____________________________________________________________________
bool HijingHadronizer::get_particles(HepMC::GenEvent* evt) {
  HepMC::GenVertex* vertice;

  vector<HepMC::GenParticle*> particles;
  vector<int> mother_ids;
  vector<HepMC::GenVertex*> prods;

  vertice = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), 0);
  evt->add_vertex(vertice);
  if (!evt->signal_process_vertex())
    evt->set_signal_process_vertex(vertice);

  const unsigned int knumpart = himain1.natt;

  for (unsigned int ipart = 0; ipart < knumpart; ipart++) {
    int mid = himain2.katt[2][ipart] - 1;  // careful of fortan to c++ array index

    particles.push_back(build_hijing(ipart, ipart + 1));
    prods.push_back(build_hijing_vertex(ipart, 0));
    mother_ids.push_back(mid);
    LogDebug("DecayChain") << "Mother index : " << mid;
  }

  LogDebug("Hijing") << "Number of particles in vector " << particles.size();

  for (unsigned int ipart = 0; ipart < particles.size(); ipart++) {
    HepMC::GenParticle* part = particles[ipart];

    int mid = mother_ids[ipart];
    LogDebug("DecayChain") << "Particle " << ipart;
    LogDebug("DecayChain") << "Mother's ID " << mid;
    LogDebug("DecayChain") << "Particle's PDG ID " << part->pdg_id();

    // remove zero pT particles from list, protection for fastJet against pt=0 jets
    if (part->status() == 1 &&
        sqrt(part->momentum().px() * part->momentum().px() + part->momentum().py() * part->momentum().py()) == 0)
      continue;

    if (mid <= 0) {
      vertice->add_particle_out(part);
      continue;
    }

    if (mid > 0) {
      HepMC::GenParticle* mother = particles[mid];
      LogDebug("DecayChain") << "Mother's PDG ID " << mother->pdg_id();
      HepMC::GenVertex* prod_vertex = mother->end_vertex();
      if (!prod_vertex) {
        prod_vertex = prods[ipart];
        prod_vertex->add_particle_in(mother);

        evt->add_vertex(prod_vertex);
        prods[ipart] = nullptr;  // mark to protect deletion
      }
      prod_vertex->add_particle_out(part);
    }
  }

  // cleanup vertices not assigned to evt
  for (unsigned int i = 0; i < prods.size(); i++) {
    if (prods[i])
      delete prods[i];
  }

  return true;
}

//_____________________________________________________________________
bool HijingHadronizer::call_hijset(
    double efrm, std::string frame, std::string proj, std::string targ, int iap, int izp, int iat, int izt) {
  float ef = efrm;
  // initialize hydjet
  HIJSET(ef,
         frame.data(),
         proj.data(),
         targ.data(),
         iap,
         izp,
         iat,
         izt,
         strlen(frame.data()),
         strlen(proj.data()),
         strlen(targ.data()));
  return true;
}

//______________________________________________________________
bool HijingHadronizer::initializeForInternalPartons() {
  //initialize pythia5

  // std::string dumstr = "";
  // call_pygive(dumstr);

  // initialize hijing
  LogInfo("HIJINGinAction") << "##### Calling HIJSET(" << efrm_ << "," << frame_ << "," << proj_ << "," << targ_ << ","
                            << iap_ << "," << izp_ << "," << iat_ << "," << izt_ << ") ####";
  call_hijset(efrm_, frame_, proj_, targ_, iap_, izp_, iat_, izt_);

  return true;
}

bool HijingHadronizer::declareStableParticles(const std::vector<int>& pdg) { return true; }

//________________________________________________________________
void HijingHadronizer::rotateEvtPlane() {
  phi0_ = 2. * pi * gen::hijran_(nullptr) - pi;
  sinphi0_ = sin(phi0_);
  cosphi0_ = cos(phi0_);
}

//________________________________________________________________
bool HijingHadronizer::hadronize() { return false; }

bool HijingHadronizer::decay() { return true; }

bool HijingHadronizer::residualDecay() { return true; }

void HijingHadronizer::finalizeEvent() { return; }

void HijingHadronizer::statistics() { return; }

const char* HijingHadronizer::classname() const { return "gen::HijingHadronizer"; }