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/HeavyIon.h>
00016 #include <HepMC/PdfInfo.h>
00017 #include <HepMC/Units.h>
00018 
00019 #include "FWCore/Framework/interface/Event.h"
00020 #include "FWCore/Framework/interface/Run.h"
00021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00023 #include "FWCore/ServiceRegistry/interface/Service.h"
00024 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
00025 #include "FWCore/Utilities/interface/EDMException.h"
00026 #include "FWCore/Framework/interface/Frameworkfwd.h"
00027 #include "FWCore/Framework/interface/EDProducer.h"
00028 #include "FWCore/Framework/interface/ESHandle.h"
00029 
00030 #include "FWCore/Framework/interface/Event.h"
00031 #include "FWCore/Framework/interface/MakerMacros.h"
00032 
00033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00034 
00035 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00036 using namespace edm;
00037 using namespace std;
00038 using namespace gen;
00039 
00040 extern "C"
00041 {
00042   float gen::rangen_()
00043   {
00044     float a = float(gFlatDistribution_->fire());
00045     return a;
00046   }
00047 
00048   double gen::drangen_(int *idummy)
00049   {
00050     double a = gFlatDistribution_->fire();
00051     return a;
00052   }
00053 }
00054 
00055 ReggeGribovPartonMCHadronizer::ReggeGribovPartonMCHadronizer(const ParameterSet &pset) :
00056   BaseHadronizer(pset),
00057   pset_(pset),
00058   m_BeamMomentum(pset.getParameter<double>("beammomentum")),
00059   m_TargetMomentum(pset.getParameter<double>("targetmomentum")),
00060   m_BeamID(pset.getParameter<int>("beamid")),
00061   m_TargetID(pset.getParameter<int>("targetid")),
00062   m_HEModel(pset.getParameter<int>("model")),
00063   m_bMin(pset.getParameter<double>("bmin")),
00064   m_bMax(pset.getParameter<double>("bmax")),
00065   m_ParamFileName(pset.getUntrackedParameter<string>("paramFileName")),
00066   m_NEvent(0),
00067   m_ImpactParameter(0.)
00068 {
00069   Service<RandomNumberGenerator> rng;
00070   if ( ! rng.isAvailable())
00071     {
00072       throw cms::Exception("Configuration")
00073         << "ReggeGribovPartonMCHadronizer requires the RandomNumberGeneratorService\n"
00074         "which is not present in the configuration file.  You must add the service\n"
00075         "in the configuration file or remove the modules that require it.";
00076     }
00077 
00078   gFlatDistribution_.reset(new CLHEP::RandFlat(rng->getEngine(), 0., 1.));
00079 
00080   int  nevet       = 1; //needed for CS
00081   int  noTables    = 0; //don't calculate tables
00082   int  LHEoutput   = 0; //no lhe
00083   int  dummySeed   = 123;
00084   char dummyName[] = "dummy";
00085   crmc_set_f_(nevet,dummySeed,m_BeamMomentum,m_TargetMomentum,m_BeamID,
00086               m_TargetID,m_HEModel,noTables,LHEoutput,dummyName,
00087               m_ParamFileName.fullPath().c_str());
00088 
00089   //additionally initialise tables
00090   initializeTablePaths();
00091 
00092   //change impact paramter
00093   nucl2_.bminim = float(m_bMin);
00094   nucl2_.bmaxim = float(m_bMax);
00095 
00096   //use set parameters to init models
00097   crmc_init_f_();
00098 }
00099 
00100 
00101 //_____________________________________________________________________
00102 ReggeGribovPartonMCHadronizer::~ReggeGribovPartonMCHadronizer()
00103 {
00104   // destructor
00105   gFlatDistribution_.reset();
00106 }
00107 
00108 //_____________________________________________________________________
00109 bool ReggeGribovPartonMCHadronizer::generatePartonsAndHadronize()
00110 {
00111   int iout=0,ievent=0;
00112   crmc_f_(iout,ievent,m_NParticles,m_ImpactParameter,
00113          m_PartID[0],m_PartPx[0],m_PartPy[0],m_PartPz[0],
00114           m_PartEnergy[0],m_PartMass[0],m_PartStatus[0]);
00115   LogDebug("ReggeGribovPartonMCInterface") << "event generated" << endl;
00116 
00117   HepMC::GenEvent* evt = new HepMC::GenEvent();
00118 
00119   evt->set_event_number(m_NEvent++);
00120   int sig_id = -1;
00121   switch (int(c2evt_.typevt)) // if negative typevt mini plasma was created by event (except -4)
00122     {
00123     case  1: sig_id = 101; break;
00124     case -1: sig_id = 101; break;
00125     case  2: sig_id = 105; break;
00126     case -2: sig_id = 105; break;
00127     case  3: sig_id = 102; break;
00128     case -3: sig_id = 102; break;
00129     case  4: sig_id = 103; break;
00130     case -4: sig_id = 104; break;
00131     default: LogDebug("ReggeGribovPartonMCInterface") << "Signal ID not recognised for setting HEPEVT" << endl;
00132     }
00133   evt->set_signal_process_id(sig_id); //an integer ID uniquely specifying the signal process (i.e. MSUB in Pythia)
00134 
00135   //create event structure;
00136   HepMC::GenVertex* theVertex = new HepMC::GenVertex();
00137   evt->add_vertex(theVertex);
00138 
00139   //number of beam particles
00140   for(int i = 0; i < m_NParticles; i++)
00141     {
00142       //consistency check
00143       const double e2 = m_PartEnergy[i] * m_PartEnergy[i];
00144       const double pc2 = m_PartPy[i]*m_PartPy[i] + m_PartPx[i]*m_PartPx[i] + m_PartPz[i]*m_PartPz[i];
00145       if (e2 + 1e-9 < pc2 )
00146         LogWarning("ReggeGribovPartonMCInterface")
00147           << "momentum off  Id:" << m_PartID[i]
00148           << "(" << i << ") "
00149           << sqrt(fabs(e2 - pc2))
00150           << endl;
00151 
00152       //add particle. do not delete. not stored as a copy
00153       HepMC::GenParticle* p = new HepMC::GenParticle(HepMC::FourVector(m_PartPx[i],
00154                                                                        m_PartPy[i],
00155                                                                        m_PartPz[i],
00156                                                                        m_PartEnergy[i]),
00157                                                      m_PartID[i],
00158                                                      m_PartStatus[i]);
00159       theVertex->add_particle_out(p);
00160     }
00161 
00162   if (m_TargetID + m_BeamID > 2) //other than pp
00163     {
00164       HepMC::HeavyIon ion(-1, //cevt_.koievt, // FIXME // Number of hard scatterings
00165                           cevt_.npjevt,                // Number of projectile participants
00166                           cevt_.ntgevt,                // Number of target participants
00167                           cevt_.kolevt,                // Number of NN (nucleon-nucleon) collisions
00168                           cevt_.npnevt + cevt_.ntnevt, // Number of spectator neutrons
00169                           cevt_.nppevt + cevt_.ntpevt, // Number of spectator protons
00170                           -1, //c2evt_.ng1evt, //FIXME // Number of N-Nwounded collisions
00171                           -1, //c2evt_.ng2evt, //FIXME // Number of Nwounded-N collisons
00172                           -1, //c2evt_.ikoevt, //FIXME // Number of Nwounded-Nwounded collisions
00173                           cevt_.bimevt,                // Impact Parameter(fm) of collision
00174                           cevt_.phievt,                // Azimuthal angle of event plane
00175                           c2evt_.fglevt,               // eccentricity of participating nucleons
00176                           hadr5_.sigineaa);            // nucleon-nucleon inelastic
00177       evt->set_heavy_ion(ion);
00178     }
00179   LogDebug("ReggeGribovPartonMCInterface") << "HepEvt and vertex constructed" << endl;
00180   //evt->print();
00181   event().reset(evt);
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_1/src/GeneratorInterface/ReggeGribovPartonMCInterface/src/ReggeGribovPartonMCHadronizer.cc
