HijingInterface/src/HijingHadronizer.cc

0001 #include <iostream>
0002 #include <cmath>
0003
0004 #include "FWCore/Concurrency/interface/SharedResourceNames.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/Run.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Utilities/interface/EDMException.h"
0010
0011 #include "GeneratorInterface/HijingInterface/interface/HijingHadronizer.h"
0012 #include "GeneratorInterface/HijingInterface/interface/HijingPythiaWrapper.h"
0013 #include "GeneratorInterface/HijingInterface/interface/HijingWrapper.h"
0014 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0015 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0016 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0017
0018 #include "HepMC/GenEvent.h"
0019 #include "HepMC/HeavyIon.h"
0020 #include "HepMC/SimpleVector.h"
0021
0022 #include "CLHEP/Random/RandomEngine.h"
0023
0024 static const double pi = 3.14159265358979;
0025
0026 using namespace edm;
0027 using namespace std;
0028 using namespace gen;
0029
0030 static CLHEP::HepRandomEngine* hijRandomEngine;
0031
0032 extern "C" {
0033 float gen::hijran_(int* idummy) { return hijRandomEngine->flat(); }
0034 }
0035
0036 extern "C" {
0037 float ran_(unsigned int* iseed) {
0038   return hijRandomEngine->flat();
0039   //      return ranff_(iseed);
0040   //      return gen::pyr_(0);
0041 }
0042 }
0043
0044 extern "C" {
0045 float rlu_(unsigned int* iseed) {
0046   return hijRandomEngine->flat();
0047   //      return ranff_(iseed);
0048   //      return gen::pyr_(0);
0049 }
0050 }
0051
0052 const std::vector<std::string> HijingHadronizer::theSharedResources = {edm::SharedResourceNames::kPythia6};
0053
0054 HijingHadronizer::HijingHadronizer(const ParameterSet& pset)
0055     : BaseHadronizer(pset),
0056       evt(nullptr),
0057       pset_(pset),
0058       bmax_(pset.getParameter<double>("bMax")),
0059       bmin_(pset.getParameter<double>("bMin")),
0060       efrm_(pset.getParameter<double>("comEnergy")),
0061       frame_(pset.getParameter<string>("frame")),
0062       proj_(pset.getParameter<string>("proj")),
0063       targ_(pset.getParameter<string>("targ")),
0064       iap_(pset.getParameter<int>("iap")),
0065       izp_(pset.getParameter<int>("izp")),
0066       iat_(pset.getParameter<int>("iat")),
0067       izt_(pset.getParameter<int>("izt")),
0068       phi0_(0.),
0069       sinphi0_(0.),
0070       cosphi0_(1.),
0071       rotate_(pset.getParameter<bool>("rotateEventPlane")) {
0072   // Default constructor
0073 }
0074
0075 //_____________________________________________________________________
0076 HijingHadronizer::~HijingHadronizer() {
0077   // destructor
0078 }
0079
0080 //_____________________________________________________________________
0081 void HijingHadronizer::doSetRandomEngine(CLHEP::HepRandomEngine* v) { hijRandomEngine = v; }
0082
0083 //_____________________________________________________________________
0084 void HijingHadronizer::add_heavy_ion_rec(HepMC::GenEvent* evt) {
0085   // heavy ion record in the final CMSSW Event
0086   HepMC::HeavyIon* hi = new HepMC::HeavyIon(himain1.jatt,  // Ncoll_hard/N of SubEvents
0087                                             himain1.np,    // Npart_proj
0088                                             himain1.nt,    // Npart_targ
0089                                             himain1.n0 + himain1.n01 + himain1.n10 + himain1.n11,  // Ncoll
0090                                             0,                                                     // spectator_neutrons
0091                                             0,                                                     // spectator_protons
0092                                             himain1.n01,  // N_Nwounded_collisions
0093                                             himain1.n10,  // Nwounded_N_collisions
0094                                             himain1.n11,  // Nwounded_Nwounded_collisions
0095                                             //gsfs Changed from 19 to 18 (Fortran counts from 1 , not 0)
0096                                             hiparnt.hint1[18],  // impact_parameter in [fm]
0097                                             phi0_,              // event_plane_angle
0098                                             0,                  // eccentricity
0099                                             //gsfs Changed from 12 to 11 (Fortran counts from 1 , not 0)
0100                                             hiparnt.hint1[11]  // sigma_inel_NN
0101   );
0102   evt->set_heavy_ion(*hi);
0103   delete hi;
0104 }
0105
0106 //___________________________________________________________________
0107 HepMC::GenParticle* HijingHadronizer::build_hijing(int index, int barcode) {
0108   // Build particle object corresponding to index in hijing
0109
0110   double x0 = himain2.patt[0][index];
0111   double y0 = himain2.patt[1][index];
0112
0113   double x = x0 * cosphi0_ - y0 * sinphi0_;
0114   double y = y0 * cosphi0_ + x0 * sinphi0_;
0115
0116   // Hijing gives V0's status=4, they need to have status=1 to be decayed in geant
0117   // also change status=11 to status=2
0118   if (himain2.katt[3][index] <= 10 && himain2.katt[3][index] > 0)
0119     himain2.katt[3][index] = 1;
0120   if (himain2.katt[3][index] <= 20 && himain2.katt[3][index] > 10)
0121     himain2.katt[3][index] = 2;
0122
0123   HepMC::GenParticle* p = new HepMC::GenParticle(HepMC::FourVector(x,                        // px
0124                                                                    y,                        // py
0125                                                                    himain2.patt[2][index],   // pz
0126                                                                    himain2.patt[3][index]),  // E
0127                                                  himain2.katt[0][index],                     // id
0128                                                  himain2.katt[3][index]                      // status
0129   );
0130   p->suggest_barcode(barcode);
0131
0132   return p;
0133 }
0134
0135 //___________________________________________________________________
0136 HepMC::GenVertex* HijingHadronizer::build_hijing_vertex(int i, int id) {
0137   // build verteces for the hijing stored events
0138   double x0 = himain2.vatt[0][i];
0139   double y0 = himain2.vatt[1][i];
0140   double x = x0 * cosphi0_ - y0 * sinphi0_;
0141   double y = y0 * cosphi0_ + x0 * sinphi0_;
0142   double z = himain2.vatt[2][i];
0143   double t = himain2.vatt[3][i];
0144
0145   HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(x, y, z, t), id);
0146   return vertex;
0147 }
0148
0149 bool HijingHadronizer::generatePartonsAndHadronize() {
0150   // generate single event
0151   if (rotate_)
0152     rotateEvtPlane();
0153
0154   // generate a HIJING event
0155
0156   float f_bmin = bmin_;
0157   float f_bmax = bmax_;
0158   HIJING(frame_.c_str(), f_bmin, f_bmax, frame_.length());
0159
0160   // event information
0161   HepMC::GenEvent* evt = new HepMC::GenEvent();
0162   get_particles(evt);
0163
0164   //   evt->set_signal_process_id(pypars.msti[0]);      // type of the process
0165   //   evt->set_event_scale(pypars.pari[16]);           // Q^2
0166   add_heavy_ion_rec(evt);
0167
0168   event().reset(evt);
0169
0170   return true;
0171 }
0172
0173 //_____________________________________________________________________
0174 bool HijingHadronizer::get_particles(HepMC::GenEvent* evt) {
0175   HepMC::GenVertex* vertice;
0176
0177   vector<HepMC::GenParticle*> particles;
0178   vector<int> mother_ids;
0179   vector<HepMC::GenVertex*> prods;
0180
0181   vertice = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), 0);
0182   evt->add_vertex(vertice);
0183   if (!evt->signal_process_vertex())
0184     evt->set_signal_process_vertex(vertice);
0185
0186   const unsigned int knumpart = himain1.natt;
0187
0188   for (unsigned int ipart = 0; ipart < knumpart; ipart++) {
0189     int mid = himain2.katt[2][ipart] - 1;  // careful of fortan to c++ array index
0190
0191     particles.push_back(build_hijing(ipart, ipart + 1));
0192     prods.push_back(build_hijing_vertex(ipart, 0));
0193     mother_ids.push_back(mid);
0194     LogDebug("DecayChain") << "Mother index : " << mid;
0195   }
0196
0197   LogDebug("Hijing") << "Number of particles in vector " << particles.size();
0198
0199   for (unsigned int ipart = 0; ipart < particles.size(); ipart++) {
0200     HepMC::GenParticle* part = particles[ipart];
0201
0202     int mid = mother_ids[ipart];
0203     LogDebug("DecayChain") << "Particle " << ipart;
0204     LogDebug("DecayChain") << "Mother's ID " << mid;
0205     LogDebug("DecayChain") << "Particle's PDG ID " << part->pdg_id();
0206
0207     // remove zero pT particles from list, protection for fastJet against pt=0 jets
0208     if (part->status() == 1 &&
0209         sqrt(part->momentum().px() * part->momentum().px() + part->momentum().py() * part->momentum().py()) == 0)
0210       continue;
0211
0212     if (mid <= 0) {
0213       vertice->add_particle_out(part);
0214       continue;
0215     }
0216
0217     if (mid > 0) {
0218       HepMC::GenParticle* mother = particles[mid];
0219       LogDebug("DecayChain") << "Mother's PDG ID " << mother->pdg_id();
0220       HepMC::GenVertex* prod_vertex = mother->end_vertex();
0221       if (!prod_vertex) {
0222         prod_vertex = prods[ipart];
0223         prod_vertex->add_particle_in(mother);
0224
0225         evt->add_vertex(prod_vertex);
0226         prods[ipart] = nullptr;  // mark to protect deletion
0227       }
0228       prod_vertex->add_particle_out(part);
0229     }
0230   }
0231
0232   // cleanup vertices not assigned to evt
0233   for (unsigned int i = 0; i < prods.size(); i++) {
0234     if (prods[i])
0235       delete prods[i];
0236   }
0237
0238   return true;
0239 }
0240
0241 //_____________________________________________________________________
0242 bool HijingHadronizer::call_hijset(
0243     double efrm, std::string frame, std::string proj, std::string targ, int iap, int izp, int iat, int izt) {
0244   float ef = efrm;
0245   // initialize hydjet
0246   frame.resize(4, ' ');
0247   proj.resize(4, ' ');
0248   targ.resize(4, ' ');
0249   HIJSET(ef, frame.data(), proj.data(), targ.data(), iap, izp, iat, izt, frame.length(), proj.length(), targ.length());
0250   return true;
0251 }
0252
0253 //______________________________________________________________
0254 bool HijingHadronizer::initializeForInternalPartons() {
0255   //initialize pythia5
0256
0257   // std::string dumstr = "";
0258   // call_pygive(dumstr);
0259
0260   // initialize hijing
0261   LogInfo("HIJINGinAction") << "##### Calling HIJSET(" << efrm_ << "," << frame_ << "," << proj_ << "," << targ_ << ","
0262                             << iap_ << "," << izp_ << "," << iat_ << "," << izt_ << ") ####";
0263   call_hijset(efrm_, frame_, proj_, targ_, iap_, izp_, iat_, izt_);
0264
0265   return true;
0266 }
0267
0268 bool HijingHadronizer::declareStableParticles(const std::vector<int>& pdg) { return true; }
0269
0270 //________________________________________________________________
0271 void HijingHadronizer::rotateEvtPlane() {
0272   phi0_ = 2. * pi * gen::hijran_(nullptr) - pi;
0273   sinphi0_ = sin(phi0_);
0274   cosphi0_ = cos(phi0_);
0275 }
0276
0277 //________________________________________________________________
0278 bool HijingHadronizer::hadronize() { return false; }
0279
0280 bool HijingHadronizer::decay() { return true; }
0281
0282 bool HijingHadronizer::residualDecay() { return true; }
0283
0284 void HijingHadronizer::finalizeEvent() { return; }
0285
0286 void HijingHadronizer::statistics() { return; }
0287
0288 const char* HijingHadronizer::classname() const { return "gen::HijingHadronizer"; }