ReggeGribovPartonMCInterface/src/ReggeGribovPartonMCHadronizer.cc

0001 #include <iostream>
0002 #include <cmath>
0003 #include <string>
0004
0005 #include "GeneratorInterface/ReggeGribovPartonMCInterface/interface/ReggeGribovPartonMCHadronizer.h"
0006 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0007 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0008 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0009
0010 #include "CLHEP/Random/RandomEngine.h"
0011
0012 #include <HepMC/GenCrossSection.h>
0013 #include <HepMC/GenEvent.h>
0014 #include <HepMC/GenVertex.h>
0015 #include <HepMC/GenParticle.h>
0016 #include <HepMC/HeavyIon.h>
0017 #include <HepMC/PdfInfo.h>
0018 #include <HepMC/Units.h>
0019
0020 #include "FWCore/Framework/interface/Event.h"
0021 #include "FWCore/Framework/interface/Run.h"
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030
0031 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0032
0033 using namespace edm;
0034 using namespace std;
0035 using namespace gen;
0036
0037 #include "GeneratorInterface/ReggeGribovPartonMCInterface/interface/IO_EPOS.h"
0038 #include "GeneratorInterface/ReggeGribovPartonMCInterface/interface/EPOS_Wrapper.h"
0039
0040 EPOS::IO_EPOS conv;
0041
0042 static CLHEP::HepRandomEngine* reggeGribovRandomEngine;
0043
0044 extern "C" {
0045 float gen::rangen_() {
0046   float a = float(reggeGribovRandomEngine->flat());
0047   return a;
0048 }
0049
0050 double gen::drangen_(int* idummy) {
0051   double a = reggeGribovRandomEngine->flat();
0052   return a;
0053 }
0054 }
0055
0056 ReggeGribovPartonMCHadronizer::ReggeGribovPartonMCHadronizer(const ParameterSet& pset)
0057     : BaseHadronizer(pset),
0058       pset_(pset),
0059       m_BeamMomentum(pset.getParameter<double>("beammomentum")),
0060       m_TargetMomentum(pset.getParameter<double>("targetmomentum")),
0061       m_BeamID(pset.getParameter<int>("beamid")),
0062       m_TargetID(pset.getParameter<int>("targetid")),
0063       m_HEModel(pset.getParameter<int>("model")),
0064       m_bMin(pset.getParameter<double>("bmin")),
0065       m_bMax(pset.getParameter<double>("bmax")),
0066       m_ParamFileName(pset.getUntrackedParameter<string>("paramFileName")),
0067       m_SkipNuclFrag(pset.getParameter<bool>("skipNuclFrag")),
0068       m_NEvent(0),
0069       m_NParticles(0),
0070       m_ImpactParameter(0.),
0071       m_IsInitialized(false) {
0072   int nevet = 1;      //needed for CS
0073   int noTables = 0;   //don't calculate tables
0074   int LHEoutput = 0;  //no lhe
0075   int dummySeed = 123;
0076   char dummyName[] = "dummy";
0077   crmc_set_f_(nevet,
0078               dummySeed,
0079               m_BeamMomentum,
0080               m_TargetMomentum,
0081               m_BeamID,
0082               m_TargetID,
0083               m_HEModel,
0084               noTables,
0085               LHEoutput,
0086               dummyName,
0087               m_ParamFileName.fullPath().c_str());
0088
0089   //additionally initialise tables
0090   initializeTablePaths();
0091
0092   //change impact paramter
0093   nucl2_.bminim = float(m_bMin);
0094   nucl2_.bmaxim = float(m_bMax);
0095 }
0096
0097 //_____________________________________________________________________
0098 ReggeGribovPartonMCHadronizer::~ReggeGribovPartonMCHadronizer() {
0099   // destructor
0100 }
0101
0102 //_____________________________________________________________________
0103 void ReggeGribovPartonMCHadronizer::doSetRandomEngine(CLHEP::HepRandomEngine* v) { reggeGribovRandomEngine = v; }
0104
0105 //_____________________________________________________________________
0106 bool ReggeGribovPartonMCHadronizer::generatePartonsAndHadronize() {
0107   int iout = 0, ievent = 0;
0108   crmc_f_(iout,
0109           ievent,
0110           m_NParticles,
0111           m_ImpactParameter,
0112           m_PartID[0],
0113           m_PartPx[0],
0114           m_PartPy[0],
0115           m_PartPz[0],
0116           m_PartEnergy[0],
0117           m_PartMass[0],
0118           m_PartStatus[0]);
0119   LogDebug("ReggeGribovPartonMCInterface") << "event generated" << endl;
0120
0121   const bool isHeavyIon = (m_TargetID + m_BeamID > 2);
0122
0123   if (isHeavyIon)
0124     conv.set_trust_beam_particles(false);
0125
0126   conv.set_skip_nuclear_fragments(m_SkipNuclFrag);
0127
0128   HepMC::GenEvent* evt = conv.read_next_event();
0129
0130   evt->set_event_number(m_NEvent++);
0131   int sig_id = -1;
0132   switch (int(c2evt_.typevt))  // if negative typevt mini plasma was created by event (except -4)
0133   {
0134     case 0:
0135       break;  //unknown for qgsjetII
0136     case 1:
0137       sig_id = 101;
0138       break;
0139     case -1:
0140       sig_id = 101;
0141       break;
0142     case 2:
0143       sig_id = 105;
0144       break;
0145     case -2:
0146       sig_id = 105;
0147       break;
0148     case 3:
0149       sig_id = 102;
0150       break;
0151     case -3:
0152       sig_id = 102;
0153       break;
0154     case 4:
0155       sig_id = 103;
0156       break;
0157     case -4:
0158       sig_id = 104;
0159       break;
0160     default:
0161       LogDebug("ReggeGribovPartonMCInterface") << "Signal ID not recognised for setting HEPEVT" << endl;
0162   }
0163   evt->set_signal_process_id(sig_id);  //an integer ID uniquely specifying the signal process (i.e. MSUB in Pythia)
0164
0165 #ifdef HEPMC_HAS_CROSS_SECTION
0166   // set cross section information for this event
0167   HepMC::GenCrossSection theCrossSection;
0168   theCrossSection.set_cross_section(double(hadr5_.sigineaa) * 1e9);  //required in pB
0169   evt->set_cross_section(theCrossSection);
0170 #endif
0171
0172   if (isHeavyIon)  //other than pp
0173   {
0174     HepMC::HeavyIon ion(cevt_.kohevt,                 // Number of hard scatterings
0175                         cevt_.npjevt,                 // Number of projectile participants
0176                         cevt_.ntgevt,                 // Number of target participants
0177                         cevt_.kolevt,                 // Number of NN (nucleon-nucleon) collisions
0178                         cevt_.npnevt + cevt_.ntnevt,  // Number of spectator neutrons
0179                         cevt_.nppevt + cevt_.ntpevt,  // Number of spectator protons
0180                         -1,                           // Number of N-Nwounded collisions
0181                         -1,                           // Number of Nwounded-N collisons
0182                         -1,                           // Number of Nwounded-Nwounded collisions
0183                         cevt_.bimevt,                 // Impact Parameter(fm) of collision
0184                         cevt_.phievt,                 // Azimuthal angle of event plane
0185                         c2evt_.fglevt,                // eccentricity of participating nucleons
0186                         hadr5_.sigine * 1e9);         // nucleon-nucleon inelastic (in pB)
0187     evt->set_heavy_ion(ion);
0188   }
0189
0190   event().reset(evt);
0191   //evt->print();
0192   //EPOS::EPOS_Wrapper::print_hepcom();
0193
0194   return true;
0195 }
0196
0197 //_____________________________________________________________________
0198 bool ReggeGribovPartonMCHadronizer::hadronize() { return false; }
0199
0200 bool ReggeGribovPartonMCHadronizer::decay() { return true; }
0201
0202 bool ReggeGribovPartonMCHadronizer::residualDecay() { return true; }
0203
0204 void ReggeGribovPartonMCHadronizer::finalizeEvent() { return; }
0205
0206 void ReggeGribovPartonMCHadronizer::statistics() { return; }
0207
0208 const char* ReggeGribovPartonMCHadronizer::classname() const { return "gen::ReggeGribovPartonMCHadronizer"; }
0209
0210 bool ReggeGribovPartonMCHadronizer::declareStableParticles(const std::vector<int>&) { return true; }
0211
0212 bool ReggeGribovPartonMCHadronizer::initializeForInternalPartons() {
0213   if (!m_IsInitialized) {
0214     //use set parameters to init models
0215     crmc_init_f_();
0216     m_IsInitialized = true;
0217   }
0218   return true;
0219 }
0220
0221 bool ReggeGribovPartonMCHadronizer::initializeTablePaths() {
0222   //epos
0223   string path_fnii(FileInPath("GeneratorInterface/ReggeGribovPartonMCInterface/data/epos.initl").fullPath());
0224   string path_fnie(FileInPath("GeneratorInterface/ReggeGribovPartonMCInterface/data/epos.iniev").fullPath());
0225   string path_fnrj(FileInPath("GeneratorInterface/ReggeGribovPartonMCInterface/data/epos.inirj").fullPath());
0226   string path_fncs(FileInPath("GeneratorInterface/ReggeGribovPartonMCInterface/data/epos.inics").fullPath());
0227
0228   if (path_fnii.length() >= 500)
0229     LogError("ReggeGribovPartonMCInterface") << "table path too long" << endl;
0230   else
0231     nfname_.nfnii = path_fnii.length();
0232   if (path_fnie.length() >= 500)
0233     LogError("ReggeGribovPartonMCInterface") << "table path too long" << endl;
0234   else
0235     nfname_.nfnie = path_fnie.length();
0236   if (path_fnrj.length() >= 500)
0237     LogError("ReggeGribovPartonMCInterface") << "table path too long" << endl;
0238   else
0239     nfname_.nfnrj = path_fnrj.length();
0240   if (path_fncs.length() >= 500)
0241     LogError("ReggeGribovPartonMCInterface") << "table path too long" << endl;
0242   else
0243     nfname_.nfncs = path_fncs.length();
0244
0245   strcpy(fname_.fnii, path_fnii.c_str());
0246   strcpy(fname_.fnie, path_fnie.c_str());
0247   strcpy(fname_.fnrj, path_fnrj.c_str());
0248   strcpy(fname_.fncs, path_fncs.c_str());
0249
0250   //qgsjet
0251   string path_fndat(FileInPath("GeneratorInterface/ReggeGribovPartonMCInterface/data/qgsjet.dat").fullPath());
0252   string path_fnncs(FileInPath("GeneratorInterface/ReggeGribovPartonMCInterface/data/qgsjet.ncs").fullPath());
0253
0254   if (path_fndat.length() >= 500)
0255     LogError("ReggeGribovPartonMCInterface") << "table path too long" << endl;
0256   else
0257     qgsnfname_.nfndat = path_fndat.length();
0258   if (path_fnncs.length() >= 500)
0259     LogError("ReggeGribovPartonMCInterface") << "table path too long" << endl;
0260   else
0261     qgsnfname_.nfnncs = path_fnncs.length();
0262
0263   strcpy(qgsfname_.fndat, path_fndat.c_str());
0264   strcpy(qgsfname_.fnncs, path_fnncs.c_str());
0265
0266   qgsfname_.ifdat = 1;  //option to overwrite the normal path
0267   qgsfname_.ifncs = 2;
0268
0269   //qgsjetII
0270   string path_fniidat(FileInPath("GeneratorInterface/ReggeGribovPartonMCInterface/data/qgsdat-II-04.lzma").fullPath());
0271   string path_fniincs(FileInPath("GeneratorInterface/ReggeGribovPartonMCInterface/data/sectnu-II-04").fullPath());
0272
0273   if (path_fniidat.length() >= 500)
0274     LogError("ReggeGribovPartonMCInterface") << "table path too long" << endl;
0275   else
0276     qgsiinfname_.nfniidat = path_fniidat.length();
0277   if (path_fniincs.length() >= 500)
0278     LogError("ReggeGribovPartonMCInterface") << "table path too long" << endl;
0279   else
0280     qgsiinfname_.nfniincs = path_fniincs.length();
0281
0282   strcpy(qgsiifname_.fniidat, path_fniidat.c_str());
0283   strcpy(qgsiifname_.fniincs, path_fniincs.c_str());
0284
0285   qgsiifname_.ifiidat = 1;  //option to overwrite the normal path
0286   qgsiifname_.ifiincs = 2;
0287
0288   return true;
0289 }