Reference Manual

Go to the documentation of this file.
00001 #include <iostream>
00002 #include <cmath>
00003 #include <string>
00004 
00005 #include "GeneratorInterface/ReggeGribovPartonMCInterface/interface/ReggeGribovPartonMCHadronizer.h"
00006 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00007 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
00008 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
00009 
00010 #include "CLHEP/Random/RandomEngine.h"
00011 #include "CLHEP/Random/RandFlat.h"
00012 
00013 #include <HepMC/GenCrossSection.h>
00014 #include <HepMC/GenEvent.h>
00015 #include <HepMC/GenVertex.h>
00016 #include <HepMC/GenParticle.h>
00017 #include <HepMC/HeavyIon.h>
00018 #include <HepMC/PdfInfo.h>
00019 #include <HepMC/Units.h>
00020 
00021 #include "FWCore/Framework/interface/Event.h"
00022 #include "FWCore/Framework/interface/Run.h"
00023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00025 #include "FWCore/ServiceRegistry/interface/Service.h"
00026 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
00027 #include "FWCore/Utilities/interface/EDMException.h"
00028 #include "FWCore/Framework/interface/Frameworkfwd.h"
00029 #include "FWCore/Framework/interface/EDProducer.h"
00030 #include "FWCore/Framework/interface/ESHandle.h"
00031 
00032 #include "FWCore/Framework/interface/Event.h"
00033 #include "FWCore/Framework/interface/MakerMacros.h"
00034 
00035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00036 
00037 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00038 using namespace edm;
00039 using namespace std;
00040 using namespace gen;
00041 
00042 #include "GeneratorInterface/ReggeGribovPartonMCInterface/interface/IO_EPOS.h"
00043 #include "GeneratorInterface/ReggeGribovPartonMCInterface/interface/EPOS_Wrapper.h"
00044 
00045 EPOS::IO_EPOS conv;
00046 
00047 extern "C"
00048 {
00049   float gen::rangen_()
00050   {
00051     float a = float(gFlatDistribution_->fire());
00052     return a;
00053   }
00054 
00055   double gen::drangen_(int *idummy)
00056   {
00057     double a = gFlatDistribution_->fire();
00058     return a;
00059   }
00060 }
00061 
00062 ReggeGribovPartonMCHadronizer::ReggeGribovPartonMCHadronizer(const ParameterSet &pset) :
00063   BaseHadronizer(pset),
00064   pset_(pset),
00065   m_BeamMomentum(pset.getParameter<double>("beammomentum")),
00066   m_TargetMomentum(pset.getParameter<double>("targetmomentum")),
00067   m_BeamID(pset.getParameter<int>("beamid")),
00068   m_TargetID(pset.getParameter<int>("targetid")),
00069   m_HEModel(pset.getParameter<int>("model")),
00070   m_bMin(pset.getParameter<double>("bmin")),
00071   m_bMax(pset.getParameter<double>("bmax")),
00072   m_ParamFileName(pset.getUntrackedParameter<string>("paramFileName")),
00073   m_SkipNuclFrag(pset.getParameter<bool>("skipNuclFrag")),
00074   m_NEvent(0),
00075   m_NParticles(0),
00076   m_ImpactParameter(0.)
00077 {
00078   Service<RandomNumberGenerator> rng;
00079   if ( ! rng.isAvailable())
00080     {
00081       throw cms::Exception("Configuration")
00082         << "ReggeGribovPartonMCHadronizer requires the RandomNumberGeneratorService\n"
00083         "which is not present in the configuration file.  You must add the service\n"
00084         "in the configuration file or remove the modules that require it.";
00085     }
00086 
00087   gFlatDistribution_.reset(new CLHEP::RandFlat(rng->getEngine(), 0., 1.));
00088 
00089   int  nevet       = 1; //needed for CS
00090   int  noTables    = 0; //don't calculate tables
00091   int  LHEoutput   = 0; //no lhe
00092   int  dummySeed   = 123;
00093   char dummyName[] = "dummy";
00094   crmc_set_f_(nevet,dummySeed,m_BeamMomentum,m_TargetMomentum,m_BeamID,
00095               m_TargetID,m_HEModel,noTables,LHEoutput,dummyName,
00096               m_ParamFileName.fullPath().c_str());
00097 
00098   //additionally initialise tables
00099   initializeTablePaths();
00100 
00101   //change impact paramter
00102   nucl2_.bminim = float(m_bMin);
00103   nucl2_.bmaxim = float(m_bMax);
00104 
00105   //use set parameters to init models
00106   crmc_init_f_();
00107 }
00108 
00109 
00110 //_____________________________________________________________________
00111 ReggeGribovPartonMCHadronizer::~ReggeGribovPartonMCHadronizer()
00112 {
00113   // destructor
00114   gFlatDistribution_.reset();
00115 }
00116 
00117 //_____________________________________________________________________
00118 bool ReggeGribovPartonMCHadronizer::generatePartonsAndHadronize()
00119 {
00120   int iout=0,ievent=0;
00121   crmc_f_(iout,ievent,m_NParticles,m_ImpactParameter,
00122          m_PartID[0],m_PartPx[0],m_PartPy[0],m_PartPz[0],
00123           m_PartEnergy[0],m_PartMass[0],m_PartStatus[0]);
00124   LogDebug("ReggeGribovPartonMCInterface") << "event generated" << endl;
00125 
00126   const bool isHeavyIon =  (m_TargetID + m_BeamID > 2);
00127 
00128   if(isHeavyIon)
00129     conv.set_trust_beam_particles(false);
00130 
00131   conv.set_skip_nuclear_fragments(m_SkipNuclFrag);
00132 
00133   HepMC::GenEvent* evt = conv.read_next_event();
00134 
00135   evt->set_event_number(m_NEvent++);
00136   int sig_id = -1;
00137   switch (int(c2evt_.typevt)) // if negative typevt mini plasma was created by event (except -4)
00138     {
00139     case  0: break; //unknown for qgsjetII
00140     case  1: sig_id = 101; break;
00141     case -1: sig_id = 101; break;
00142     case  2: sig_id = 105; break;
00143     case -2: sig_id = 105; break;
00144     case  3: sig_id = 102; break;
00145     case -3: sig_id = 102; break;
00146     case  4: sig_id = 103; break;
00147     case -4: sig_id = 104; break;
00148     default: LogDebug("ReggeGribovPartonMCInterface") << "Signal ID not recognised for setting HEPEVT" << endl;
00149     }
00150   evt->set_signal_process_id(sig_id); //an integer ID uniquely specifying the signal process (i.e. MSUB in Pythia)
00151 
00152 #ifdef HEPMC_HAS_CROSS_SECTION
00153   // set cross section information for this event
00154   HepMC::GenCrossSection theCrossSection;
00155   theCrossSection.set_cross_section(double(hadr5_.sigineaa)*1e9); //required in pB
00156   evt->set_cross_section(theCrossSection);
00157 #endif
00158 
00159   if (isHeavyIon) //other than pp
00160     {
00161       HepMC::HeavyIon ion(cevt_.kohevt,                // Number of hard scatterings
00162                           cevt_.npjevt,                // Number of projectile participants
00163                           cevt_.ntgevt,                // Number of target participants
00164                           cevt_.kolevt,                // Number of NN (nucleon-nucleon) collisions
00165                           cevt_.npnevt + cevt_.ntnevt, // Number of spectator neutrons
00166                           cevt_.nppevt + cevt_.ntpevt, // Number of spectator protons
00167                           -1,                          // Number of N-Nwounded collisions
00168                           -1,                          // Number of Nwounded-N collisons
00169                           -1,                          // Number of Nwounded-Nwounded collisions
00170                           cevt_.bimevt,                // Impact Parameter(fm) of collision
00171                           cevt_.phievt,                // Azimuthal angle of event plane
00172                           c2evt_.fglevt,               // eccentricity of participating nucleons
00173                           hadr5_.sigine*1e9);        // nucleon-nucleon inelastic (in pB)
00174       evt->set_heavy_ion(ion);
00175     }
00176 
00177 
00178   event().reset(evt);
00179   //evt->print();
00180   //EPOS::EPOS_Wrapper::print_hepcom();
00181 
00182   return true;
00183 }
00184 
00185 //_____________________________________________________________________
00186 bool ReggeGribovPartonMCHadronizer::hadronize()
00187 {
00188    return false;
00189 }
00190 
00191 bool ReggeGribovPartonMCHadronizer::decay()
00192 {
00193    return true;
00194 }
00195 
00196 bool ReggeGribovPartonMCHadronizer::residualDecay()
00197 {
00198    return true;
00199 }
00200 
00201 void ReggeGribovPartonMCHadronizer::finalizeEvent(){
00202     return;
00203 }
00204 
00205 void ReggeGribovPartonMCHadronizer::statistics(){
00206     return;
00207 }
00208 
00209 const char* ReggeGribovPartonMCHadronizer::classname() const
00210 {
00211    return "gen::ReggeGribovPartonMCHadronizer";
00212 }
00213 
00214 bool ReggeGribovPartonMCHadronizer::declareStableParticles ( const std::vector<int> )
00215 {
00216  return true;
00217  }
00218 
00219 bool ReggeGribovPartonMCHadronizer::initializeForInternalPartons()
00220 {
00221  return true;
00222  }
00223 
00224 bool ReggeGribovPartonMCHadronizer::initializeTablePaths()
00225 {
00226   //epos
00227   string path_fnii(FileInPath("GeneratorInterface/ReggeGribovPartonMCInterface/data/epos.initl").fullPath());
00228   string path_fnie(FileInPath("GeneratorInterface/ReggeGribovPartonMCInterface/data/epos.iniev").fullPath());
00229   string path_fnrj(FileInPath("GeneratorInterface/ReggeGribovPartonMCInterface/data/epos.inirj").fullPath());
00230   string path_fncs(FileInPath("GeneratorInterface/ReggeGribovPartonMCInterface/data/epos.inics").fullPath());
00231 
00232   if (path_fnii.length() >= 500) LogError("ReggeGribovPartonMCInterface") << "table path too long" << endl;
00233   else nfname_.nfnii = path_fnii.length();
00234   if (path_fnie.length() >= 500) LogError("ReggeGribovPartonMCInterface") << "table path too long" << endl;
00235   else nfname_.nfnie = path_fnie.length();
00236   if (path_fnrj.length() >= 500) LogError("ReggeGribovPartonMCInterface") << "table path too long" << endl;
00237   else nfname_.nfnrj = path_fnrj.length();
00238   if (path_fncs.length() >= 500) LogError("ReggeGribovPartonMCInterface") << "table path too long" << endl;
00239   else nfname_.nfncs = path_fncs.length();
00240 
00241   strcpy(fname_.fnii, path_fnii.c_str());
00242   strcpy(fname_.fnie, path_fnie.c_str());
00243   strcpy(fname_.fnrj, path_fnrj.c_str());
00244   strcpy(fname_.fncs, path_fncs.c_str());
00245 
00246   //qgsjet
00247   string path_fndat(FileInPath("GeneratorInterface/ReggeGribovPartonMCInterface/data/qgsjet.dat").fullPath());
00248   string path_fnncs(FileInPath("GeneratorInterface/ReggeGribovPartonMCInterface/data/qgsjet.ncs").fullPath());
00249 
00250   if (path_fndat.length() >= 500) LogError("ReggeGribovPartonMCInterface") << "table path too long" << endl;
00251   else qgsnfname_.nfndat = path_fndat.length();
00252   if (path_fnncs.length() >= 500) LogError("ReggeGribovPartonMCInterface") << "table path too long" << endl;
00253   else qgsnfname_.nfnncs = path_fnncs.length();
00254 
00255   strcpy(qgsfname_.fndat, path_fndat.c_str());
00256   strcpy(qgsfname_.fnncs, path_fnncs.c_str());
00257 
00258   qgsfname_.ifdat=1; //option to overwrite the normal path
00259   qgsfname_.ifncs=2;
00260 
00261   //qgsjetII
00262   string path_fniidat(FileInPath("GeneratorInterface/ReggeGribovPartonMCInterface/data/qgsdat-II-04.lzma").fullPath());
00263   string path_fniincs(FileInPath("GeneratorInterface/ReggeGribovPartonMCInterface/data/sectnu-II-04").fullPath());
00264 
00265   if (path_fniidat.length() >= 500) LogError("ReggeGribovPartonMCInterface") << "table path too long" << endl;
00266   else qgsiinfname_.nfniidat = path_fniidat.length();
00267   if (path_fniincs.length() >= 500) LogError("ReggeGribovPartonMCInterface") << "table path too long" << endl;
00268   else qgsiinfname_.nfniincs = path_fniincs.length();
00269 
00270   strcpy(qgsiifname_.fniidat, path_fniidat.c_str());
00271   strcpy(qgsiifname_.fniincs, path_fniincs.c_str());
00272 
00273   qgsiifname_.ifiidat=1; //option to overwrite the normal path
00274   qgsiifname_.ifiincs=2;
00275 
00276  return true;
00277 }
/data/refman/pasoursint/CMSSW_6_1_2_SLHC4_patch1/src/GeneratorInterface/ReggeGribovPartonMCInterface/src/ReggeGribovPartonMCHadronizer.cc
