CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
List of all members | Public Member Functions | Private Member Functions | Private Attributes | Static Private Attributes
gen::Hydjet2Hadronizer Class Reference

#include <Hydjet2Hadronizer.h>

Inheritance diagram for gen::Hydjet2Hadronizer:
InitialState gen::BaseHadronizer

Public Member Functions

double CharmEnhancementFactor (double, double, double, double)
 
const char * classname () const
 
bool decay ()
 
bool declareSpecialSettings (const std::vector< std::string > &)
 
bool declareStableParticles (const std::vector< int > &)
 
double f (double)
 
double f2 (double, double, double)
 
void finalizeEvent ()
 
bool generatePartonsAndHadronize ()
 
double GetVolEff ()
 
double GetWeakDecayLimit () override
 
bool hadronize ()
 
 Hydjet2Hadronizer (const edm::ParameterSet &)
 
bool IniOfThFreezeoutParameters ()
 
bool initializeForExternalPartons ()
 
bool initializeForInternalPartons ()
 
double MidpointIntegrator2 (double, double, double, double)
 
bool readSettings (int)
 
bool residualDecay ()
 
bool RunDecays () override
 
void SetVolEff (double value)
 
double SimpsonIntegrator (double, double, double, double)
 
double SimpsonIntegrator2 (double, double, double, double)
 
void statistics ()
 
 ~Hydjet2Hadronizer () override
 
- Public Member Functions inherited from InitialState
virtual void Evolve (List_t &secondaries, ParticleAllocator &allocator, double weakDecayLimit)
 
 InitialState ()
 
virtual ~InitialState () noexcept(false)
 
- Public Member Functions inherited from gen::BaseHadronizer
 BaseHadronizer (edm::ParameterSet const &ps)
 
void cleanLHE ()
 
void generateLHE (edm::LuminosityBlock const &lumi, CLHEP::HepRandomEngine *rengine, unsigned int ncpu)
 
edm::EventgetEDMEvent () const
 
std::unique_ptr< HepMC::GenEventgetGenEvent ()
 
std::unique_ptr
< GenEventInfoProduct
getGenEventInfo ()
 
virtual std::unique_ptr
< GenLumiInfoHeader
getGenLumiInfoHeader () const
 
GenRunInfoProductgetGenRunInfo ()
 
std::unique_ptr< lhef::LHEEventgetLHEEvent ()
 
const std::shared_ptr
< lhef::LHERunInfo > & 
getLHERunInfo () const
 
const std::string & gridpackPath () const
 
int randomIndex () const
 
const std::string & randomInitConfigDescription () const
 
void randomizeIndex (edm::LuminosityBlock const &lumi, CLHEP::HepRandomEngine *rengine)
 
void resetEvent (std::unique_ptr< HepMC::GenEvent > event)
 
void resetEventInfo (std::unique_ptr< GenEventInfoProduct > eventInfo)
 
virtual bool select (HepMC::GenEvent *) const
 
void setEDMEvent (edm::Event &event)
 
void setLHEEvent (std::unique_ptr< lhef::LHEEvent > event)
 
void setLHERunInfo (std::unique_ptr< lhef::LHERunInfo > runInfo)
 
void setRandomEngine (CLHEP::HepRandomEngine *v)
 
std::vector< std::string > const & sharedResources () const
 
virtual ~BaseHadronizer () noexcept(false)
 

Private Member Functions

void add_heavy_ion_rec (HepMC::GenEvent *evt)
 
HepMC::GenParticle * build_hyjet2 (int index, int barcode)
 
HepMC::GenVertex * build_hyjet2_vertex (int i, int id)
 
void doSetRandomEngine (CLHEP::HepRandomEngine *v) override
 
std::vector< std::string > const & doSharedResources () const override
 
bool get_particles (HepMC::GenEvent *evt)
 
double nuclear_radius () const
 
void rotateEvtPlane ()
 

Private Attributes

ParticleAllocator allocator
 
float Bgen
 
double cosphi0_
 
float E [500000]
 
bool embedding_
 
int ev
 
HepMC::GenEventevt
 
double fAw
 
double fBfix
 
double fBmax
 
double fBmin
 
int fCharmProd
 
double fCorrC
 
double fCorrS
 
int fDecay
 
double fDelta
 
double fEpsilon
 
int fEtaType
 
int fIanglu
 
int fIenglu
 
int fIfb
 
int fIfDeltaEpsilon
 
int final [500000]
 
int FirstDaughterIndex [500000]
 
int fIshad
 
double fMu_th_pip
 
double fMuB
 
double fMuC
 
double fMuI3
 
double fMuS
 
double fMuTh [1000]
 
double fNccth
 
int fNf
 
int fNhsel
 
double fNocth
 
int fNPartTypes
 
int fPartEnc [1000]
 
double fPartMu [2000]
 
double fPartMult [2000]
 
double fPtmin
 
int fPyhist
 
int fPythDecay
 
double fR
 
edm::Service< TFileServicefs
 
double fSigmaTau
 
double fSqrtS
 
double fT
 
double fT0
 
double fTau
 
double fTau0
 
double fThFO
 
int fTMuType
 
double fUmax
 
double fVolEff
 
double fWeakDecay
 
double fYlmax
 
int Index [500000]
 
int LastDaughterIndex [500000]
 
unsigned int maxEventsToPrint_
 
int MotherIndex [500000]
 
int Mpdg [500000]
 
int Nbcol
 
int NDaughters [500000]
 
int nhard_
 
int Nhyd
 
int Njet
 
int Npart
 
int Npyt
 
int nsoft_
 
int nsub_
 
int Ntot
 
int pdg [500000]
 
double phi0_
 
edm::ParameterSet pset
 
double psiforv3
 
float Px [500000]
 
float Py [500000]
 
Pythia6Servicepythia6Service_
 
unsigned int pythiaPylistVerbosity_
 
int pythiaStatus [500000]
 
float Pz [500000]
 
bool rotate_
 
float Sigin
 
float Sigjet
 
double sinphi0_
 
List_t source
 
edm::InputTag src_
 
int sseed
 
float T [500000]
 
int type [500000]
 
float X [500000]
 
float Y [500000]
 
float Z [500000]
 

Static Private Attributes

static const std::vector
< std::string > 
theSharedResources
 

Additional Inherited Members

- Protected Member Functions inherited from gen::BaseHadronizer
std::unique_ptr
< HepMC::GenEvent > & 
event ()
 
std::unique_ptr
< GenEventInfoProduct > & 
eventInfo ()
 
lhef::LHEEventlheEvent ()
 
lhef::LHERunInfolheRunInfo ()
 
GenRunInfoProductrunInfo ()
 
- Protected Attributes inherited from InitialState
DatabasePDGfDatabase
 
- Protected Attributes inherited from gen::BaseHadronizer
std::string lheFile_
 
int randomIndex_
 

Detailed Description

Definition at line 65 of file Hydjet2Hadronizer.h.

Constructor & Destructor Documentation

Hydjet2Hadronizer::Hydjet2Hadronizer ( const edm::ParameterSet pset)

Definition at line 131 of file Hydjet2Hadronizer.cc.

References embedding_, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), LogDebug, maxEventsToPrint_, pythiaPylistVerbosity_, and src_.

132  : BaseHadronizer(pset),
133  fSqrtS(pset.getParameter<double>("fSqrtS")), // C.m.s. energy per nucleon pair
134  fAw(pset.getParameter<double>("fAw")), // Atomic weigth of nuclei, fAw
135  fIfb(pset.getParameter<int>(
136  "fIfb")), // Flag of type of centrality generation, fBfix (=0 is fixed by fBfix, >0 distributed [fBfmin, fBmax])
137  fBmin(pset.getParameter<double>("fBmin")), // Minimum impact parameter in units of nuclear radius, fBmin
138  fBmax(pset.getParameter<double>("fBmax")), // Maximum impact parameter in units of nuclear radius, fBmax
139  fBfix(pset.getParameter<double>("fBfix")), // Fixed impact parameter in units of nuclear radius, fBfix
140  fT(pset.getParameter<double>("fT")), // Temperature at chemical freeze-out, fT [GeV]
141  fMuB(pset.getParameter<double>("fMuB")), // Chemical baryon potential per unit charge, fMuB [GeV]
142  fMuS(pset.getParameter<double>("fMuS")), // Chemical strangeness potential per unit charge, fMuS [GeV]
143  fMuC(pset.getParameter<double>(
144  "fMuC")), // Chemical charm potential per unit charge, fMuC [GeV] (used if charm production is turned on)
145  fMuI3(pset.getParameter<double>("fMuI3")), // Chemical isospin potential per unit charge, fMuI3 [GeV]
146  fThFO(pset.getParameter<double>("fThFO")), // Temperature at thermal freeze-out, fThFO [GeV]
147  fMu_th_pip(pset.getParameter<double>(
148  "fMu_th_pip")), // Chemical potential of pi+ at thermal freeze-out, fMu_th_pip [GeV]
149  fTau(pset.getParameter<double>(
150  "fTau")), // Proper time proper at thermal freeze-out for central collisions, fTau [fm/c]
151  fSigmaTau(pset.getParameter<double>(
152  "fSigmaTau")), // Duration of emission at thermal freeze-out for central collisions, fSigmaTau [fm/c]
153  fR(pset.getParameter<double>(
154  "fR")), // Maximal transverse radius at thermal freeze-out for central collisions, fR [fm]
155  fYlmax(pset.getParameter<double>("fYlmax")), // Maximal longitudinal flow rapidity at thermal freeze-out, fYlmax
156  fUmax(pset.getParameter<double>(
157  "fUmax")), // Maximal transverse flow rapidity at thermal freeze-out for central collisions, fUmax
158  fDelta(pset.getParameter<double>(
159  "fDelta")), // Momentum azimuthal anizotropy parameter at thermal freeze-out, fDelta
160  fEpsilon(pset.getParameter<double>(
161  "fEpsilon")), // Spatial azimuthal anisotropy parameter at thermal freeze-out, fEpsilon
162  fIfDeltaEpsilon(pset.getParameter<double>(
163  "fIfDeltaEpsilon")), // Flag to specify fDelta and fEpsilon values, fIfDeltaEpsilon (=0 user's ones, >=1 calculated)
164  fDecay(pset.getParameter<int>(
165  "fDecay")), // Flag to switch on/off hadron decays, fDecay (=0 decays off, >=1 decays on)
166  fWeakDecay(pset.getParameter<double>(
167  "fWeakDecay")), // Low decay width threshold fWeakDecay[GeV]: width<fWeakDecay decay off, width>=fDecayWidth decay on; can be used to switch off weak decays
168  fEtaType(pset.getParameter<double>(
169  "fEtaType")), // Flag to choose longitudinal flow rapidity distribution, fEtaType (=0 uniform, >0 Gaussian with the dispersion Ylmax)
170  fTMuType(pset.getParameter<double>(
171  "fTMuType")), // Flag to use calculated T_ch, mu_B and mu_S as a function of fSqrtS, fTMuType (=0 user's ones, >0 calculated)
172  fCorrS(pset.getParameter<double>(
173  "fCorrS")), // Strangeness supression factor gamma_s with fCorrS value (0<fCorrS <=1, if fCorrS <= 0 then it is calculated)
174  fCharmProd(pset.getParameter<int>(
175  "fCharmProd")), // Flag to include thermal charm production, fCharmProd (=0 no charm production, >=1 charm production)
176  fCorrC(pset.getParameter<double>(
177  "fCorrC")), // Charmness enhancement factor gamma_c with fCorrC value (fCorrC >0, if fCorrC<0 then it is calculated)
178  fNhsel(pset.getParameter<int>(
179  "fNhsel")), //Flag to include jet (J)/jet quenching (JQ) and hydro (H) state production, fNhsel (0 H on & J off, 1 H/J on & JQ off, 2 H/J/HQ on, 3 J on & H/JQ off, 4 H off & J/JQ on)
180  fPyhist(pset.getParameter<int>(
181  "fPyhist")), // Flag to suppress the output of particle history from PYTHIA, fPyhist (=1 only final state particles; =0 full particle history from PYTHIA)
182  fIshad(pset.getParameter<int>(
183  "fIshad")), // Flag to switch on/off nuclear shadowing, fIshad (0 shadowing off, 1 shadowing on)
184  fPtmin(pset.getParameter<double>(
185  "fPtmin")), // Minimal pt of parton-parton scattering in PYTHIA event, fPtmin [GeV/c]
186  fT0(pset.getParameter<double>(
187  "fT0")), // Initial QGP temperature for central Pb+Pb collisions in mid-rapidity, fT0 [GeV]
188  fTau0(pset.getParameter<double>("fTau0")), // Proper QGP formation time in fm/c, fTau0 (0.01<fTau0<10)
189  fNf(pset.getParameter<int>("fNf")), // Number of active quark flavours in QGP, fNf (0, 1, 2 or 3)
190  fIenglu(pset.getParameter<int>(
191  "fIenglu")), // Flag to fix type of partonic energy loss, fIenglu (0 radiative and collisional loss, 1 radiative loss only, 2 collisional loss only)
192  fIanglu(pset.getParameter<int>(
193  "fIanglu")), // Flag to fix type of angular distribution of in-medium emitted gluons, fIanglu (0 small-angular, 1 wide-angular, 2 collinear).
194 
195  embedding_(pset.getParameter<bool>("embeddingMode")),
196  rotate_(pset.getParameter<bool>("rotateEventPlane")),
197  evt(nullptr),
198  nsub_(0),
199  nhard_(0),
200  nsoft_(0),
201  phi0_(0.),
202  sinphi0_(0.),
203  cosphi0_(1.),
205 
206 {
207  // constructor
208  // PYLIST Verbosity Level
209  // Valid PYLIST arguments are: 1, 2, 3, 5, 7, 11, 12, 13
210  pythiaPylistVerbosity_ = pset.getUntrackedParameter<int>("pythiaPylistVerbosity", 0);
211  LogDebug("PYLISTverbosity") << "Pythia PYLIST verbosity level = " << pythiaPylistVerbosity_;
212  //Max number of events printed on verbosity level
213  maxEventsToPrint_ = pset.getUntrackedParameter<int>("maxEventsToPrint", 0);
214  LogDebug("Events2Print") << "Number of events to be printed = " << maxEventsToPrint_;
215  if (embedding_)
216  src_ = pset.getParameter<edm::InputTag>("backgroundLabel");
217 }
T getUntrackedParameter(std::string const &, T const &) const
BaseHadronizer(edm::ParameterSet const &ps)
unsigned int pythiaPylistVerbosity_
Pythia6Service * pythia6Service_
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
#define LogDebug(id)
Hydjet2Hadronizer::~Hydjet2Hadronizer ( )
override

Definition at line 220 of file Hydjet2Hadronizer.cc.

References pythia6Service_.

220  {
221  // destructor
222  call_pystat(1);
223  delete pythia6Service_;
224 }
Pythia6Service * pythia6Service_

Member Function Documentation

void Hydjet2Hadronizer::add_heavy_ion_rec ( HepMC::GenEvent evt)
private

Definition at line 1212 of file Hydjet2Hadronizer.cc.

References Bgen, Nbcol, npart, Npart, nsub_, nuclear_radius(), phi0_, and Sigin.

Referenced by generatePartonsAndHadronize().

1212  {
1213  // heavy ion record in the final CMSSW Event
1214  double npart = Npart;
1215  int nproj = static_cast<int>(npart / 2);
1216  int ntarg = static_cast<int>(npart - nproj);
1217 
1218  HepMC::HeavyIon* hi = new HepMC::HeavyIon(nsub_, // Ncoll_hard/N of SubEvents
1219  nproj, // Npart_proj
1220  ntarg, // Npart_targ
1221  Nbcol, // Ncoll
1222  0, // spectator_neutrons
1223  0, // spectator_protons
1224  0, // N_Nwounded_collisions
1225  0, // Nwounded_N_collisions
1226  0, // Nwounded_Nwounded_collisions
1227  Bgen * nuclear_radius(), // impact_parameter in [fm]
1228  phi0_, // event_plane_angle
1229  0, // eccentricity
1230  Sigin // sigma_inel_NN
1231  );
1232 
1233  evt->set_heavy_ion(*hi);
1234  delete hi;
1235 }
double npart
Definition: HydjetWrapper.h:46
double nuclear_radius() const
HepMC::GenParticle * Hydjet2Hadronizer::build_hyjet2 ( int  index,
int  barcode 
)
private

Definition at line 1174 of file Hydjet2Hadronizer.cc.

References cosphi0_, E, GenParticle::GenParticle, gen::p, pdg, Px, Py, Pz, and sinphi0_.

Referenced by get_particles().

1174  {
1175  // Build particle object corresponding to index in hyjets (soft+hard)
1176 
1177  double px0 = Px[index];
1178  double py0 = Py[index];
1179 
1180  double px = px0 * cosphi0_ - py0 * sinphi0_;
1181  double py = py0 * cosphi0_ + px0 * sinphi0_;
1182  // cout<< "status: "<<convertStatus(final[index], type[index])<<endl;
1183  HepMC::GenParticle* p = new HepMC::GenParticle(HepMC::FourVector(px, // px
1184  py, // py
1185  Pz[index], // pz
1186  E[index]), // E
1187  pdg[index], // id
1188  convertStatusForComponents(final[index], type[index]) // status
1189 
1190  );
1191 
1192  p->suggest_barcode(barcode);
1193  return p;
1194 }
double p[5][pyjets_maxn]
HepMC::GenVertex * Hydjet2Hadronizer::build_hyjet2_vertex ( int  i,
int  id 
)
private

Definition at line 1197 of file Hydjet2Hadronizer.cc.

References cosphi0_, mps_fire::i, sinphi0_, submitPVValidationJobs::t, X, Y, and Z.

Referenced by get_particles().

1197  {
1198  // build verteces for the hyjets stored events
1199 
1200  double x0 = X[i];
1201  double y0 = Y[i];
1202  double x = x0 * cosphi0_ - y0 * sinphi0_;
1203  double y = y0 * cosphi0_ + x0 * sinphi0_;
1204  double z = Z[i];
1205  double t = T[i];
1206 
1207  HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(x, y, z, t), id);
1208  return vertex;
1209 }
long double T
double Hydjet2Hadronizer::CharmEnhancementFactor ( double  Ncc,
double  Ndth,
double  NJPsith,
double  Epsilon 
)

Definition at line 1055 of file Hydjet2Hadronizer.cc.

References BesselI0(), BesselI1(), and mps_fire::i.

Referenced by generatePartonsAndHadronize().

1055  {
1056  double gammaC = 100.;
1057  double x1 = gammaC * Ndth;
1058  double var1 = Ncc - 0.5 * gammaC * Ndth * TMath::BesselI1(x1) / TMath::BesselI0(x1) - gammaC * gammaC * NJPsith;
1059  LogInfo("Charam") << "gammaC 20"
1060  << " var " << var1 << endl;
1061  gammaC = 1.;
1062  double x0 = gammaC * Ndth;
1063  double var0 = Ncc - 0.5 * gammaC * Ndth * TMath::BesselI1(x0) / TMath::BesselI0(x0) - gammaC * gammaC * NJPsith;
1064  LogInfo("Charam") << "gammaC 1"
1065  << " var " << var0;
1066 
1067  for (int i = 1; i < 1000; i++) {
1068  if (var1 * var0 < 0) {
1069  gammaC = gammaC + 0.01 * i;
1070  double x = gammaC * Ndth;
1071  var0 = Ncc - 0.5 * gammaC * Ndth * TMath::BesselI1(x) / TMath::BesselI0(x) - gammaC * gammaC * NJPsith;
1072  } else {
1073  LogInfo("Charam") << "gammaC " << gammaC << " var0 " << var0;
1074  return gammaC;
1075  }
1076  }
1077  LogInfo("Charam") << "gammaC not found ? " << gammaC << " var0 " << var0;
1078  return -100;
1079 }
double BesselI0(double x)
Log< level::Info, false > LogInfo
double BesselI1(double x)
const char * Hydjet2Hadronizer::classname ( ) const

Definition at line 968 of file Hydjet2Hadronizer.cc.

968 { return "gen::Hydjet2Hadronizer"; }
bool Hydjet2Hadronizer::decay ( )

Definition at line 964 of file Hydjet2Hadronizer.cc.

964 { return true; }
bool gen::Hydjet2Hadronizer::declareSpecialSettings ( const std::vector< std::string > &  )
inline

Definition at line 71 of file Hydjet2Hadronizer.h.

71 { return true; }
bool Hydjet2Hadronizer::declareStableParticles ( const std::vector< int > &  _pdg)

Definition at line 951 of file Hydjet2Hadronizer.cc.

References gen::call_pygive(), gather_cfg::cout, mps_fire::i, pdg, and gen::pycomp_().

951  {
952  std::vector<int> pdg = _pdg;
953  for (size_t i = 0; i < pdg.size(); i++) {
954  int pyCode = pycomp_(pdg[i]);
955  std::ostringstream pyCard;
956  pyCard << "MDCY(" << pyCode << ",1)=0";
957  std::cout << pyCard.str() << std::endl;
958  call_pygive(pyCard.str());
959  }
960  return true;
961 }
bool call_pygive(const std::string &line)
int pycomp_(int &)
tuple cout
Definition: gather_cfg.py:144
void Hydjet2Hadronizer::doSetRandomEngine ( CLHEP::HepRandomEngine *  v)
overrideprivatevirtual

Reimplemented from gen::BaseHadronizer.

Definition at line 227 of file Hydjet2Hadronizer.cc.

References hjRandomEngine, pythia6Service_, gen::Pythia6Service::setRandomEngine(), and gen::v.

227  {
229  hjRandomEngine = v;
230 }
double v[5][pyjets_maxn]
Pythia6Service * pythia6Service_
CLHEP::HepRandomEngine * hjRandomEngine
void setRandomEngine(CLHEP::HepRandomEngine *v)
std::vector<std::string> const& gen::Hydjet2Hadronizer::doSharedResources ( ) const
inlineoverrideprivatevirtual

Reimplemented from gen::BaseHadronizer.

Definition at line 109 of file Hydjet2Hadronizer.h.

References theSharedResources.

109 { return theSharedResources; }
static const std::vector< std::string > theSharedResources
double gen::Hydjet2Hadronizer::f ( double  )
double Hydjet2Hadronizer::f2 ( double  x,
double  y,
double  Delta 
)

Definition at line 974 of file Hydjet2Hadronizer.cc.

References alignCSCRings::ff, fR, fUmax, and LogDebug.

Referenced by MidpointIntegrator2(), and SimpsonIntegrator().

974  {
975  LogDebug("f2") << "in f2: "
976  << "delta" << Delta;
977  double RsB = fR; //test: podstavit' *coefff_RB
978  double rhou = fUmax * y / RsB;
979  double ff =
980  y * TMath::CosH(rhou) * TMath::Sqrt(1 + Delta * TMath::Cos(2 * x) * TMath::TanH(rhou) * TMath::TanH(rhou));
981  //n_mu u^mu f-la 20
982  return ff;
983 }
#define LogDebug(id)
void Hydjet2Hadronizer::finalizeEvent ( )

Definition at line 966 of file Hydjet2Hadronizer.cc.

966 {}
bool Hydjet2Hadronizer::generatePartonsAndHadronize ( )

Definition at line 480 of file Hydjet2Hadronizer.cc.

References funct::A, funct::abs(), Abs(), add_heavy_ion_rec(), ParticleAllocator::AddParticle(), allocator, HYFPARCommon::bgen, Bgen, gen::C, CharmEnhancementFactor(), funct::cos(), cosphi0_, ztail::d, alignCSCRings::e, E, embedding_, HLT_FULL_cff::Epsilon, ev, gen::BaseHadronizer::event(), edm::errors::EventCorruption, InitialState::Evolve(), evt, fBfix, fCharmProd, fCorrC, InitialState::fDatabase, fDelta, fEpsilon, fEtaType, fIfb, fIfDeltaEpsilon, FirstDaughterIndex, fNccth, fNhsel, fNocth, fNPartTypes, fPartEnc, fPartMu, fPartMult, fR, ParticleAllocator::FreeList(), fSigmaTau, fT, fTau, fThFO, fUmax, fVolEff, fYlmax, get_particles(), edm::Event::getByLabel(), ParticlePDG::GetCharmAQNumber(), ParticlePDG::GetCharmQNumber(), gen::BaseHadronizer::getEDMEvent(), ParticlePDG::GetMass(), DatabasePDG::GetPDGParticle(), ParticlePDG::GetSpin(), GetWeakDecayLimit(), h, hjRandomEngine, hyevnt_(), mps_fire::i, if(), Index, Particle::InitIndexing(), input, dqmiolumiharvest::j, kMax, LastDaughterIndex, ResonanceBuilder::mass, MotherIndex, Mpdg, RPCpg::mu, HYFPARCommon::nbcol, Nbcol, NDaughters, nhard_, Nhyd, HYJPARCommon::njet, Njet, HYPARTCommon::njp, HYFPARCommon::npart, Npart, HYFPARCommon::npart0, Npyt, nsoft_, nsub_, Ntot, pdg, phi0_, Pi, HYPARTCommon::ppart, RandArrayFunction::PrepareTable(), psiforv3, Px, Py, pypars, pythia6Service_, pythiaStatus, Pz, alignCSCRings::r, HYIPARCommon::RA, rotate_, rotateEvtPlane(), RunDecays(), SERVICE, SERVICEEV, HYJPARCommon::sigin, Sigin, HYJPARCommon::sigjet, Sigjet, SimpsonIntegrator2(), funct::sin(), sinphi0_, source, mathSSE::sqrt(), src_, FrontierCondition_GT_autoExpress_cfi::t0, metsig::tau, TwoPi, histoStyle::weight, X, Y, and Z.

480  {
481  Pythia6Service::InstanceWrapper guard(pythia6Service_);
482 
483  // Initialize the static "last index variable"
485 
486  //----- high-pt part------------------------------
487  TLorentzVector partJMom, partJPos, zeroVec;
488 
489  // generate single event
490  if (embedding_) {
491  fIfb = 0;
492  const edm::Event& e = getEDMEvent();
494  e.getByLabel(src_, input);
495  const HepMC::GenEvent* inev = input->GetEvent();
496  const HepMC::HeavyIon* hi = inev->heavy_ion();
497  if (hi) {
498  fBfix = hi->impact_parameter();
499  phi0_ = hi->event_plane_angle();
500  sinphi0_ = sin(phi0_);
501  cosphi0_ = cos(phi0_);
502  } else {
503  LogWarning("EventEmbedding") << "Background event does not have heavy ion record!";
504  }
505  } else if (rotate_)
506  rotateEvtPlane();
507 
508  nsoft_ = 0;
509  nhard_ = 0;
510  /*
511  edm::LogInfo("HYDJET2mode") << "##### HYDJET2 fNhsel = " << fNhsel;
512  edm::LogInfo("HYDJET2fpart") << "##### HYDJET2 fpart = " << hyflow.fpart;??
513  edm::LogInfo("HYDJET2tf") << "##### HYDJET2 hadron freez-out temp, Tf = " << hyflow.Tf;??
514  edm::LogInfo("HYDJET2tf") << "##### HYDJET2 hadron freez-out temp, Tf = " << hyflow.Tf;??
515  edm::LogInfo("HYDJET2inTemp") << "##### HYDJET2: QGP init temperature, fT0 ="<<fT0;
516  edm::LogInfo("HYDJET2inTau") << "##### HYDJET2: QGP formation time, fTau0 ="<<fTau0;
517  */
518  // generate a HYDJET event
519  int ntry = 0;
520  while (nsoft_ == 0 && nhard_ == 0) {
521  if (ntry > 100) {
522  LogError("Hydjet2EmptyEvent") << "##### HYDJET2: No Particles generated, Number of tries =" << ntry;
523  // Throw an exception. Use the EventCorruption exception since it maps onto SkipEvent
524  // which is what we want to do here.
525  std::ostringstream sstr;
526  sstr << "Hydjet2HadronizerProducer: No particles generated after " << ntry << " tries.\n";
527  edm::Exception except(edm::errors::EventCorruption, sstr.str());
528  throw except;
529  } else {
530  //generate non-equilibrated part event
531  hyevnt_();
532 
534  if (fNhsel != 0) {
535  //get number of particles in jets
536  int numbJetPart = HYPART.njp;
537 
538  for (int i = 0; i < numbJetPart; ++i) {
539  int pythiaStatus = int(HYPART.ppart[i][0]); // PYTHIA status code
540  int pdg = int(HYPART.ppart[i][1]); // PYTHIA species code
541  double px = HYPART.ppart[i][2]; // px
542  double py = HYPART.ppart[i][3]; // py
543  double pz = HYPART.ppart[i][4]; // pz
544  double e = HYPART.ppart[i][5]; // E
545  double vx = HYPART.ppart[i][6]; // x
546  double vy = HYPART.ppart[i][7]; // y
547  double vz = HYPART.ppart[i][8]; // z
548  double vt = HYPART.ppart[i][9]; // t
549  // particle line number in pythia are 1 based while we use a 0 based numbering
550  int mother_index = int(HYPART.ppart[i][10]) - 1; //line number of parent particle
551  int daughter_index1 = int(HYPART.ppart[i][11]) - 1; //line number of first daughter
552  int daughter_index2 = int(HYPART.ppart[i][12]) - 1; //line number of last daughter
553 
554  // For status codes 3, 13, 14 the first and last daughter indexes have a different meaning
555  // used for color flow in PYTHIA. So these indexes will be reset to zero.
556  if (TMath::Abs(daughter_index1) > numbJetPart || TMath::Abs(daughter_index2) > numbJetPart ||
557  TMath::Abs(daughter_index1) > TMath::Abs(daughter_index2)) {
558  daughter_index1 = -1;
559  daughter_index2 = -1;
560  }
561 
562  ParticlePDG* partDef = fDatabase->GetPDGParticle(pdg);
563 
564  int type = 1; //from jet
565  if (partDef) {
566  int motherPdg = int(HYPART.ppart[mother_index][1]);
567  if (motherPdg == 0)
568  motherPdg = -1;
569  partJMom.SetXYZT(px, py, pz, e);
570  partJPos.SetXYZT(vx, vy, vz, vt);
571  Particle particle(partDef, partJPos, partJMom, 0, 0, type, motherPdg, zeroVec, zeroVec);
572  int index = particle.SetIndex();
573  if (index != i) {
574  // LogWarning("Hydjet2Hadronizer") << " Allocated Hydjet2 index is not synchronized with the PYTHIA index !" << endl
575  // << " Collision history information is destroyed! It happens when a PYTHIA code is not" << endl
576  // << " implemented in Hydjet2 particle list particles.data! Check it out!";
577  }
578  particle.SetPythiaStatusCode(pythiaStatus);
579  particle.SetMother(mother_index);
580  particle.SetFirstDaughterIndex(daughter_index1);
581  particle.SetLastDaughterIndex(daughter_index2);
582  if (pythiaStatus != 1)
583  particle.SetDecayed();
584  allocator.AddParticle(particle, source);
585  } else {
586  LogWarning("Hydjet2Hadronizer")
587  << " PYTHIA particle of specie " << pdg << " is not in Hydjet2 particle list" << endl
588  << " Please define it in particles.data, otherwise the history information will be de-synchronized and "
589  "lost!";
590  }
591  }
592  } //nhsel !=0 not only hydro!
593 
594  //----------HYDRO part--------------------------------------
595 
596  // get impact parameter
597  double impactParameter = HYFPAR.bgen;
598 
599  // Sergey psiforv3
600  psiforv3 = 0.; //AS-ML Nov2012 epsilon3 //
601  double e3 = (0.2 / 5.5) * TMath::Power(impactParameter, 1. / 3.);
602  psiforv3 = TMath::TwoPi() * (-0.5 + CLHEP::RandFlat::shoot(hjRandomEngine)) / 3.;
603  SERVICEEV.psiv3 = -psiforv3;
604 
605  if (fNhsel < 3) {
606  const double weightMax = 2 * TMath::CosH(fUmax);
607  const int nBins = 100;
608  double probList[nBins];
609  RandArrayFunction arrayFunctDistE(nBins);
610  RandArrayFunction arrayFunctDistR(nBins);
611  TLorentzVector partPos, partMom, n1, p0;
612  TVector3 vec3;
613  const TLorentzVector zeroVec;
614  //set maximal hadron energy
615  const double eMax = 5.;
616  //-------------------------------------
617  // get impact parameter
618 
619  double Delta = fDelta;
620  double Epsilon = fEpsilon;
621 
622  if (fIfDeltaEpsilon > 0) {
623  double Epsilon0 = 0.5 * impactParameter; //e0=b/2Ra
624  double coeff = (HYIPAR.RA / fR) / 12.; //phenomenological coefficient
625  Epsilon = Epsilon0 * coeff;
626  double C = 5.6;
627  double A = C * Epsilon * (1 - Epsilon * Epsilon);
628  if (TMath::Abs(Epsilon) < 0.0001 || TMath::Abs(A) < 0.0001)
629  Delta = 0.0;
630  if (TMath::Abs(Epsilon) > 0.0001 && TMath::Abs(A) > 0.0001)
631  Delta = 0.5 * (TMath::Sqrt(1 + 4 * A * (Epsilon + A)) - 1) / A;
632  }
633 
634  //effective volume for central
635  double dYl = 2 * fYlmax; //uniform distr. [-Ylmax; Ylmax]
636  if (fEtaType > 0)
637  dYl = TMath::Sqrt(2 * TMath::Pi()) * fYlmax; //Gaussian distr.
638 
639  double VolEffcent = 2 * TMath::Pi() * fTau * dYl * (fR * fR) / TMath::Power((fUmax), 2) *
640  ((fUmax)*TMath::SinH((fUmax)) - TMath::CosH((fUmax)) + 1);
641 
642  //effective volume for non-central Simpson2
643  double VolEffnoncent = fTau * dYl * SimpsonIntegrator2(0., 2. * TMath::Pi(), Epsilon, Delta);
644 
645  fVolEff = VolEffcent * HYFPAR.npart / HYFPAR.npart0;
646 
647  double coeff_RB = TMath::Sqrt(VolEffcent * HYFPAR.npart / HYFPAR.npart0 / VolEffnoncent);
648  double coeff_R1 = HYFPAR.npart / HYFPAR.npart0;
649  coeff_R1 = TMath::Power(coeff_R1, 0.333333);
650 
651  double Veff = fVolEff;
652  //------------------------------------
653  //cycle on particles types
654 
655  double Nbcol = HYFPAR.nbcol;
656  double NccNN = SERVICE.charm;
657  double Ncc = Nbcol * NccNN / dYl;
658  double Nocth = fNocth;
659  double NJPsith = fNccth;
660 
661  double gammaC = 1.0;
662  if (fCorrC <= 0) {
663  gammaC = CharmEnhancementFactor(Ncc, Nocth, NJPsith, 0.001);
664  } else {
665  gammaC = fCorrC;
666  }
667 
668  LogInfo("Hydjet2Hadronizer|Param") << " gammaC = " << gammaC;
669 
670  for (int i = 0; i < fNPartTypes; ++i) {
671  double Mparam = fPartMult[2 * i] * Veff;
672  const int encoding = fPartEnc[i];
673 
674  //ml if(abs(encoding)==443)Mparam = Mparam * gammaC * gammaC;
675  //ml if(abs(encoding)==411 || abs(encoding)==421 ||abs(encoding)==413 || abs(encoding)==423
676  //ml || abs(encoding)==4122 || abs(encoding)==431)
677 
678  ParticlePDG* partDef0 = fDatabase->GetPDGParticle(encoding);
679 
680  if (partDef0->GetCharmQNumber() != 0 || partDef0->GetCharmAQNumber() != 0)
681  Mparam = Mparam * gammaC;
682  if (abs(encoding) == 443)
683  Mparam = Mparam * gammaC;
684 
685  LogInfo("Hydjet2Hadronizer|Param") << encoding << " " << Mparam / dYl;
686 
687  int multiplicity = CLHEP::RandPoisson::shoot(hjRandomEngine, Mparam);
688 
689  LogInfo("Hydjet2Hadronizer|Param")
690  << "specie: " << encoding << "; average mult: = " << Mparam << "; multiplicity = " << multiplicity;
691 
692  if (multiplicity > 0) {
693  ParticlePDG* partDef = fDatabase->GetPDGParticle(encoding);
694  if (!partDef) {
695  LogError("Hydjet2Hadronizer") << "No particle with encoding " << encoding;
696  continue;
697  }
698 
699  if (fCharmProd <= 0 && (partDef->GetCharmQNumber() != 0 || partDef->GetCharmAQNumber() != 0)) {
700  LogInfo("Hydjet2Hadronizer|Param") << "statistical charmed particle not allowed ! " << encoding;
701  continue;
702  }
703  if (partDef->GetCharmQNumber() != 0 || partDef->GetCharmAQNumber() != 0)
704  LogInfo("Hydjet2Hadronizer|Param") << " charm pdg generated " << encoding;
705 
706  //compute chemical potential for single f.o. mu==mu_ch
707  //compute chemical potential for thermal f.o.
708  double mu = fPartMu[2 * i];
709 
710  //choose Bose-Einstein or Fermi-Dirac statistics
711  const double d = !(int(2 * partDef->GetSpin()) & 1) ? -1 : 1;
712  const double mass = partDef->GetMass();
713 
714  //prepare histogram to sample hadron energy:
715  double h = (eMax - mass) / nBins;
716  double x = mass + 0.5 * h;
717  int i;
718  for (i = 0; i < nBins; ++i) {
719  if (x >= mu && fThFO > 0)
720  probList[i] = x * TMath::Sqrt(x * x - mass * mass) / (TMath::Exp((x - mu) / (fThFO)) + d);
721  if (x >= mu && fThFO <= 0)
722  probList[i] = x * TMath::Sqrt(x * x - mass * mass) / (TMath::Exp((x - mu) / (fT)) + d);
723  if (x < mu)
724  probList[i] = 0.;
725  x += h;
726  }
727  arrayFunctDistE.PrepareTable(probList);
728 
729  //prepare histogram to sample hadron transverse radius:
730  h = (fR) / nBins;
731  x = 0.5 * h;
732  double param = (fUmax) / (fR);
733  for (i = 0; i < nBins; ++i) {
734  probList[i] = x * TMath::CosH(param * x);
735  x += h;
736  }
737  arrayFunctDistR.PrepareTable(probList);
738 
739  //loop over hadrons, assign hadron coordinates and momenta
740  double weight = 0., yy = 0., px0 = 0., py0 = 0., pz0 = 0.;
741  double e = 0., x0 = 0., y0 = 0., z0 = 0., t0 = 0., etaF = 0.;
742  double r, RB, phiF;
743 
744  RB = fR * coeff_RB * coeff_R1 * TMath::Sqrt((1 + e3) / (1 - e3));
745 
746  for (int j = 0; j < multiplicity; ++j) {
747  do {
748  fEtaType <= 0 ? etaF = fYlmax * (2. * CLHEP::RandFlat::shoot(hjRandomEngine) - 1.)
749  : etaF = (fYlmax) * (CLHEP::RandGauss::shoot(hjRandomEngine));
750  n1.SetXYZT(0., 0., TMath::SinH(etaF), TMath::CosH(etaF));
751 
752  if (TMath::Abs(etaF) > 5.)
753  continue;
754 
755  //old
756  //double RBold = fR * TMath::Sqrt(1-fEpsilon);
757 
758  //RB = fR * coeff_RB * coeff_R1;
759 
760  //double impactParameter =HYFPAR.bgen;
761  //double e0 = 0.5*impactParameter;
762  //double RBold1 = fR * TMath::Sqrt(1-e0);
763 
764  double rho = TMath::Sqrt(CLHEP::RandFlat::shoot(hjRandomEngine));
765  double phi = TMath::TwoPi() * CLHEP::RandFlat::shoot(hjRandomEngine);
766  double Rx = TMath::Sqrt(1 - Epsilon) * RB;
767  double Ry = TMath::Sqrt(1 + Epsilon) * RB;
768 
769  x0 = Rx * rho * TMath::Cos(phi);
770  y0 = Ry * rho * TMath::Sin(phi);
771  r = TMath::Sqrt(x0 * x0 + y0 * y0);
772  phiF = TMath::Abs(TMath::ATan(y0 / x0));
773 
774  if (x0 < 0 && y0 > 0)
775  phiF = TMath::Pi() - phiF;
776  if (x0 < 0 && y0 < 0)
777  phiF = TMath::Pi() + phiF;
778  if (x0 > 0 && y0 < 0)
779  phiF = 2. * TMath::Pi() - phiF;
780 
781  //new Nov2012 AS-ML
782  if (r > RB * (1 + e3 * TMath::Cos(3 * (phiF + psiforv3))) / (1 + e3))
783  continue;
784 
785  //proper time with emission duration
786  double tau =
787  coeff_R1 * fTau + sqrt(2.) * fSigmaTau * coeff_R1 * (CLHEP::RandGauss::shoot(hjRandomEngine));
788  z0 = tau * TMath::SinH(etaF);
789  t0 = tau * TMath::CosH(etaF);
790  double rhou = fUmax * r / RB;
791  double rhou3 = 0.063 * TMath::Sqrt((0.5 * impactParameter) / 0.67);
792  double rhou4 = 0.023 * ((0.5 * impactParameter) / 0.67);
793  double rrcoeff = 1. / TMath::Sqrt(1. + Delta * TMath::Cos(2 * phiF));
794  //AS-ML Nov.2012
795  rhou3 = 0.; //rhou4=0.;
796  rhou = rhou * (1 + rrcoeff * rhou3 * TMath::Cos(3 * (phiF + psiforv3)) +
797  rrcoeff * rhou4 * TMath::Cos(4 * phiF)); //ML new suggestion of AS mar2012
798  double delta1 = 0.;
799  Delta = Delta * (1.0 + delta1 * TMath::Cos(phiF) - delta1 * TMath::Cos(3 * phiF));
800 
801  double uxf = TMath::SinH(rhou) * TMath::Sqrt(1 + Delta) * TMath::Cos(phiF);
802  double uyf = TMath::SinH(rhou) * TMath::Sqrt(1 - Delta) * TMath::Sin(phiF);
803  double utf = TMath::CosH(etaF) * TMath::CosH(rhou) *
804  TMath::Sqrt(1 + Delta * TMath::Cos(2 * phiF) * TMath::TanH(rhou) * TMath::TanH(rhou));
805  double uzf = TMath::SinH(etaF) * TMath::CosH(rhou) *
806  TMath::Sqrt(1 + Delta * TMath::Cos(2 * phiF) * TMath::TanH(rhou) * TMath::TanH(rhou));
807 
808  vec3.SetXYZ(uxf / utf, uyf / utf, uzf / utf);
809  n1.Boost(-vec3);
810 
811  yy = weightMax * CLHEP::RandFlat::shoot(hjRandomEngine);
812 
813  double php0 = TMath::TwoPi() * CLHEP::RandFlat::shoot(hjRandomEngine);
814  double ctp0 = 2. * CLHEP::RandFlat::shoot(hjRandomEngine) - 1.;
815  double stp0 = TMath::Sqrt(1. - ctp0 * ctp0);
816  e = mass + (eMax - mass) * arrayFunctDistE();
817  double pp0 = TMath::Sqrt(e * e - mass * mass);
818  px0 = pp0 * stp0 * TMath::Sin(php0);
819  py0 = pp0 * stp0 * TMath::Cos(php0);
820  pz0 = pp0 * ctp0;
821  p0.SetXYZT(px0, py0, pz0, e);
822 
823  //weight for rdr
824  weight = (n1 * p0) / e; // weight for rdr gammar: weight = (n1 * p0) / n1[3] / e;
825 
826  } while (yy >= weight);
827 
828  if (abs(z0) > 1000 || abs(x0) > 1000)
829  LogInfo("Hydjet2Hadronizer|Param") << " etaF = " << etaF << std::endl;
830 
831  partMom.SetXYZT(px0, py0, pz0, e);
832  partPos.SetXYZT(x0, y0, z0, t0);
833  partMom.Boost(vec3);
834 
835  int type = 0; //hydro
836  Particle particle(partDef, partPos, partMom, 0., 0, type, -1, zeroVec, zeroVec);
837  particle.SetIndex();
838  allocator.AddParticle(particle, source);
839  } //nhsel==4 , no hydro part
840  }
841  }
842  }
843 
845 
846  Npart = (int)HYFPAR.npart;
847  Bgen = HYFPAR.bgen;
848  Njet = (int)HYJPAR.njet;
849  Nbcol = (int)HYFPAR.nbcol;
850 
851  //if(source.empty()) {
852  // LogError("Hydjet2Hadronizer") << "Source is not initialized!! Trying again...";
853  //return ;
854  //}
855  //Run the decays
856  if (RunDecays())
858 
859  LPIT_t it;
860  LPIT_t e;
861 
862  //Fill the decayed arrays
863  Ntot = 0;
864  Nhyd = 0;
865  Npyt = 0;
866  for (it = source.begin(), e = source.end(); it != e; ++it) {
867  TVector3 pos(it->Pos().Vect());
868  TVector3 mom(it->Mom().Vect());
869  float m1 = it->TableMass();
870  pdg[Ntot] = it->Encoding();
871  Mpdg[Ntot] = it->GetLastMotherPdg();
872  Px[Ntot] = mom[0];
873  Py[Ntot] = mom[1];
874  Pz[Ntot] = mom[2];
875  E[Ntot] = TMath::Sqrt(mom.Mag2() + m1 * m1);
876  X[Ntot] = pos[0];
877  Y[Ntot] = pos[1];
878  Z[Ntot] = pos[2];
879  T[Ntot] = it->T();
880  type[Ntot] = it->GetType();
881  pythiaStatus[Ntot] = convertStatus(it->GetPythiaStatusCode());
882  Index[Ntot] = it->GetIndex();
883  MotherIndex[Ntot] = it->GetMother();
884  NDaughters[Ntot] = it->GetNDaughters();
885  FirstDaughterIndex[Ntot] = -1;
886  LastDaughterIndex[Ntot] = -1;
887  //index of first daughter
888  FirstDaughterIndex[Ntot] = it->GetFirstDaughterIndex();
889  //index of last daughter
890  LastDaughterIndex[Ntot] = it->GetLastDaughterIndex();
891  if (type[Ntot] == 1) { // jets
892  if (pythiaStatus[Ntot] == 1 && NDaughters[Ntot] == 0) { // code for final state particle in pythia
893  final[Ntot] = 1;
894  } else {
895  final[Ntot] = pythiaStatus[Ntot];
896  }
897  }
898  if (type[Ntot] == 0) { // hydro
899  if (NDaughters[Ntot] == 0)
900  final[Ntot] = 1;
901  else
902  final[Ntot] = 2;
903  }
904 
905  if (type[Ntot] == 0)
906  Nhyd++;
907  if (type[Ntot] == 1)
908  Npyt++;
909 
910  Ntot++;
911  if (Ntot > kMax)
912  LogError("Hydjet2Hadronizer") << "Ntot is too large" << Ntot;
913  }
914 
915  nsoft_ = Nhyd;
916  nsub_ = Njet;
917  nhard_ = Npyt;
918 
919  //100 trys
920 
921  ++ntry;
922  }
923  }
924 
925  if (ev == 0) {
926  Sigin = HYJPAR.sigin;
927  Sigjet = HYJPAR.sigjet;
928  }
929  ev = 1;
930 
931  if (fNhsel < 3)
932  nsub_++;
933 
934  // event information
936  if (nhard_ > 0 || nsoft_ > 0)
937  get_particles(evt);
938 
939  evt->set_signal_process_id(pypars.msti[0]); // type of the process
940  evt->set_event_scale(pypars.pari[16]); // Q^2
941  add_heavy_ion_rec(evt);
942 
943  event().reset(evt);
944 
946 
947  return kTRUE;
948 }
const double TwoPi
const double Pi
ParticlePDG * GetPDGParticle(int pdg)
Definition: DatabasePDG.cc:218
ParticleAllocator allocator
#define HYJPAR
Definition: HYJET_COMMONS.h:69
void add_heavy_ion_rec(HepMC::GenEvent *evt)
#define HYPART
double SimpsonIntegrator2(double, double, double, double)
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Log< level::Error, false > LogError
#define SERVICEEV
Definition: HYJET_COMMONS.h:52
double GetCharmAQNumber()
Definition: ParticlePDG.h:81
static std::string const input
Definition: EdmProvDump.cc:47
#define HYFPAR
Definition: HYJET_COMMONS.h:99
tuple d
Definition: ztail.py:151
double GetSpin()
Definition: ParticlePDG.h:73
Pythia6Service * pythia6Service_
std::vector< vec2 > vec3
Definition: HCALResponse.h:17
void FreeList(List_t &list)
Definition: Particle.cc:118
bool get_particles(HepMC::GenEvent *evt)
void hyevnt_()
T sqrt(T t)
Definition: SSEVec.h:19
if(conf_.getParameter< bool >("UseStripCablingDB"))
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
T Abs(T a)
Definition: MathUtil.h:49
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double CharmEnhancementFactor(double, double, double, double)
const int mu
Definition: Constants.h:22
std::unique_ptr< HepMC::GenEvent > & event()
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:500
CLHEP::HepRandomEngine * hjRandomEngine
#define HYIPAR
Definition: HYJET_COMMONS.h:26
virtual void Evolve(List_t &secondaries, ParticleAllocator &allocator, double weakDecayLimit)
Definition: InitialState.cc:13
static void InitIndexing()
Definition: Particle.h:140
#define SERVICE
Definition: HYJET_COMMONS.h:38
Log< level::Info, false > LogInfo
#define pypars
bool RunDecays() override
double GetCharmQNumber()
Definition: ParticlePDG.h:80
DatabasePDG * fDatabase
Definition: InitialState.h:16
double GetWeakDecayLimit() override
double GetMass()
Definition: ParticlePDG.h:70
#define kMax
void AddParticle(const Particle &particle, List_t &list)
Definition: Particle.cc:107
Log< level::Warning, false > LogWarning
int weight
Definition: histoStyle.py:51
std::list< Particle >::iterator LPIT_t
Definition: Particle.h:175
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
long double T
edm::Event & getEDMEvent() const
bool Hydjet2Hadronizer::get_particles ( HepMC::GenEvent evt)
private

Definition at line 1093 of file Hydjet2Hadronizer.cc.

References build_hyjet2(), build_hyjet2_vertex(), GenParticle::GenParticle, mps_fire::i, LogDebug, MotherIndex, nhard_, Njet, nsoft_, and nsub_.

Referenced by generatePartonsAndHadronize().

1093  {
1094  // Hard particles. The first nhard_ lines from hyjets array.
1095  // Pythia/Pyquen sub-events (sub-collisions) for a given event
1096  // Return T/F if success/failure
1097  // Create particles from lujet entries, assign them into vertices and
1098  // put the vertices in the GenEvent, for each SubEvent
1099  // The SubEvent information is kept by storing indeces of main vertices
1100  // of subevents as a vector in GenHIEvent.
1101  LogDebug("SubEvent") << "Number of sub events " << nsub_;
1102  LogDebug("Hydjet2") << "Number of hard events " << Njet;
1103  LogDebug("Hydjet2") << "Number of hard particles " << nhard_;
1104  LogDebug("Hydjet2") << "Number of soft particles " << nsoft_;
1105 
1106  vector<HepMC::GenVertex*> sub_vertices(nsub_);
1107 
1108  int ihy = 0;
1109  for (int isub = 0; isub < nsub_; isub++) {
1110  LogDebug("SubEvent") << "Sub Event ID : " << isub;
1111 
1112  int sub_up = (isub + 1) * 150000; // Upper limit in mother index, determining the range of Sub-Event
1113  vector<HepMC::GenParticle*> particles;
1114  vector<int> mother_ids;
1115  vector<HepMC::GenVertex*> prods;
1116 
1117  sub_vertices[isub] = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), isub);
1118  evt->add_vertex(sub_vertices[isub]);
1119 
1120  if (!evt->signal_process_vertex())
1121  evt->set_signal_process_vertex(sub_vertices[isub]);
1122 
1123  while (ihy < nhard_ + nsoft_ && (MotherIndex[ihy] < sub_up || ihy > nhard_)) {
1124  particles.push_back(build_hyjet2(ihy, ihy + 1));
1125  prods.push_back(build_hyjet2_vertex(ihy, isub));
1126  mother_ids.push_back(MotherIndex[ihy]);
1127  LogDebug("DecayChain") << "Mother index : " << MotherIndex[ihy];
1128  ihy++;
1129  }
1130  //Produce Vertices and add them to the GenEvent. Remember that GenParticles are adopted by
1131  //GenVertex and GenVertex is adopted by GenEvent.
1132  LogDebug("Hydjet2") << "Number of particles in vector " << particles.size();
1133 
1134  for (unsigned int i = 0; i < particles.size(); i++) {
1135  HepMC::GenParticle* part = particles[i];
1136  //The Fortran code is modified to preserve mother id info, by seperating the beginning
1137  //mother indices of successive subevents by 5000
1138  int mid = mother_ids[i] - isub * 150000 - 1;
1139  LogDebug("DecayChain") << "Particle " << i;
1140  LogDebug("DecayChain") << "Mother's ID " << mid;
1141  LogDebug("DecayChain") << "Particle's PDG ID " << part->pdg_id();
1142 
1143  if (mid <= 0) {
1144  sub_vertices[isub]->add_particle_out(part);
1145  continue;
1146  }
1147 
1148  if (mid > 0) {
1149  HepMC::GenParticle* mother = particles[mid];
1150  LogDebug("DecayChain") << "Mother's PDG ID " << mother->pdg_id();
1151  HepMC::GenVertex* prod_vertex = mother->end_vertex();
1152  if (!prod_vertex) {
1153  prod_vertex = prods[i];
1154  prod_vertex->add_particle_in(mother);
1155  evt->add_vertex(prod_vertex);
1156  prods[i] = nullptr; // mark to protect deletion
1157  }
1158 
1159  prod_vertex->add_particle_out(part);
1160  }
1161  }
1162 
1163  // cleanup vertices not assigned to evt
1164  for (unsigned int i = 0; i < prods.size(); i++) {
1165  if (prods[i])
1166  delete prods[i];
1167  }
1168  }
1169 
1170  return kTRUE;
1171 }
HepMC::GenVertex * build_hyjet2_vertex(int i, int id)
part
Definition: HCALResponse.h:20
HepMC::GenParticle * build_hyjet2(int index, int barcode)
#define LogDebug(id)
double gen::Hydjet2Hadronizer::GetVolEff ( )
inline

Definition at line 87 of file Hydjet2Hadronizer.h.

References fVolEff.

87 { return fVolEff; }
double gen::Hydjet2Hadronizer::GetWeakDecayLimit ( )
inlineoverridevirtual

Implements InitialState.

Definition at line 89 of file Hydjet2Hadronizer.h.

References fWeakDecay.

Referenced by generatePartonsAndHadronize().

89 { return fWeakDecay; }
bool Hydjet2Hadronizer::hadronize ( )

Definition at line 963 of file Hydjet2Hadronizer.cc.

963 { return false; }
bool gen::Hydjet2Hadronizer::IniOfThFreezeoutParameters ( )
bool gen::Hydjet2Hadronizer::initializeForExternalPartons ( )
bool Hydjet2Hadronizer::initializeForInternalPartons ( )

Definition at line 247 of file Hydjet2Hadronizer.cc.

References funct::abs(), Abs(), HYIPARCommon::AW, HYIPARCommon::bmaxh, HYIPARCommon::bminh, NAStrangePotential::CalculateStrangePotential(), fAw, fBfix, fBmax, fBmin, fCorrS, InitialState::fDatabase, fEtaType, fIanglu, fIenglu, fIfb, fIshad, fMu_th_pip, fMuB, fMuC, fMuI3, fMuS, fNccth, fNf, fNhsel, fNocth, fNPartTypes, fPartEnc, fPartMu, fPartMult, fPtmin, fPyhist, fR, fSqrtS, fT, fT0, fTau, fTau0, fThFO, fTMuType, fUmax, fVolEff, fYlmax, ParticlePDG::GetBaryonNumber(), ParticlePDG::GetCharmAQNumber(), ParticlePDG::GetCharmness(), ParticlePDG::GetCharmQNumber(), ParticlePDG::GetElectricCharge(), DatabasePDG::GetNParticles(), ParticlePDG::GetPDG(), DatabasePDG::GetPDGParticle(), DatabasePDG::GetPDGParticleByIndex(), ParticlePDG::GetStrangeness(), HYPYIN, HYJPARCommon::iPyhist, HYJPARCommon::ishad, RPCpg::mu, myini_(), HYJPARCommon::nhsel, nuclear_radius(), GrandCanonical::ParticleNumberDensity(), Pi, HYJPARCommon::ptmin, PYQPAR, pythia6Service_, NAStrangePotential::SetBaryonPotential(), and NAStrangePotential::SetTemperature().

247  {
248  Pythia6Service::InstanceWrapper guard(pythia6Service_);
249 
250  // the input impact parameter (bxx_) is in [fm]; transform in [fm/RA] for hydjet usage
251  const float ra = nuclear_radius();
252  LogInfo("Hydjet2Hadronizer|RAScaling") << "Nuclear radius(RA) = " << ra;
253  fBmin /= ra;
254  fBmax /= ra;
255  fBfix /= ra;
256 
257  //check and redefine input parameters
258  if (fTMuType > 0 && fSqrtS > 2.24) {
259  if (fSqrtS < 2.24) {
260  LogError("Hydjet2Hadronizer|sqrtS") << "SqrtS<2.24 not allowed with fTMuType>0";
261  return false;
262  }
263 
264  //sqrt(s) = 2.24 ==> T_kin = 0.8 GeV
265  //see J. Cleymans, H. Oeschler, K. Redlich,S. Wheaton, Phys Rev. C73 034905 (2006)
266  fMuB = 1.308 / (1. + fSqrtS * 0.273);
267  fT = 0.166 - 0.139 * fMuB * fMuB - 0.053 * fMuB * fMuB * fMuB * fMuB;
268  fMuI3 = 0.;
269  fMuS = 0.;
270 
271  //create strange potential object and set strangeness density 0
273  psp->SetBaryonPotential(fMuB);
274  psp->SetTemperature(fT);
275 
276  //compute strangeness potential
277  if (fMuB > 0.01)
279  LogInfo("Hydjet2Hadronizer|Strange") << "fMuS = " << fMuS;
280 
281  //if user choose fYlmax larger then allowed by kinematics at the specified beam energy sqrt(s)
282  if (fYlmax > TMath::Log(fSqrtS / 0.94)) {
283  LogError("Hydjet2Hadronizer|Ylmax") << "fYlmax more then TMath::Log(fSqrtS vs 0.94)!!! ";
284  return false;
285  }
286 
287  if (fCorrS <= 0.) {
288  //see F. Becattini, J. Mannien, M. Gazdzicki, Phys Rev. C73 044905 (2006)
289  fCorrS = 1. - 0.386 * TMath::Exp(-1.23 * fT / fMuB);
290  LogInfo("Hydjet2Hadronizer|Strange")
291  << "The phenomenological f-la F. Becattini et al. PRC73 044905 (2006) for CorrS was used." << endl
292  << "Strangeness suppression parameter = " << fCorrS;
293  }
294  LogInfo("Hydjet2Hadronizer|Strange")
295  << "The phenomenological f-la J. Cleymans et al. PRC73 034905 (2006) for Tch mu_B was used." << endl
296  << "The simulation will be done with the calculated parameters:" << endl
297  << "Baryon chemical potential = " << fMuB << " [GeV]" << endl
298  << "Strangeness chemical potential = " << fMuS << " [GeV]" << endl
299  << "Isospin chemical potential = " << fMuI3 << " [GeV]" << endl
300  << "Strangeness suppression parameter = " << fCorrS << endl
301  << "Eta_max = " << fYlmax;
302  }
303 
304  LogInfo("Hydjet2Hadronizer|Param") << "Used eta_max = " << fYlmax << endl
305  << "maximal allowed eta_max TMath::Log(fSqrtS/0.94)= "
306  << TMath::Log(fSqrtS / 0.94);
307 
308  //initialisation of high-pt part
309  HYJPAR.nhsel = fNhsel;
310  HYJPAR.ptmin = fPtmin;
311  HYJPAR.ishad = fIshad;
312  HYJPAR.iPyhist = fPyhist;
313  HYIPAR.bminh = fBmin;
314  HYIPAR.bmaxh = fBmax;
315  HYIPAR.AW = fAw;
316 
317  HYPYIN.ifb = fIfb;
318  HYPYIN.bfix = fBfix;
319  HYPYIN.ene = fSqrtS;
320 
321  PYQPAR.T0 = fT0;
322  PYQPAR.tau0 = fTau0;
323  PYQPAR.nf = fNf;
324  PYQPAR.ienglu = fIenglu;
325  PYQPAR.ianglu = fIanglu;
326  myini_();
327 
328  // calculation of multiplicities of different particle species
329  // according to the grand canonical approach
330  GrandCanonical gc(15, fT, fMuB, fMuS, fMuI3, fMuC);
331  GrandCanonical gc_ch(15, fT, fMuB, fMuS, fMuI3, fMuC);
332  GrandCanonical gc_pi_th(15, fThFO, 0., 0., fMu_th_pip, fMuC);
333  GrandCanonical gc_th_0(15, fThFO, 0., 0., 0., 0.);
334 
335  // std::ofstream outMult("densities.txt");
336  // outMult<<"encoding particle density chemical potential "<<std::endl;
337 
338  double Nocth = 0; //open charm
339  double NJPsith = 0; //JPsi
340 
341  //effective volume for central
342  double dYl = 2 * fYlmax; //uniform distr. [-Ylmax; Ylmax]
343  if (fEtaType > 0)
344  dYl = TMath::Sqrt(2 * TMath::Pi()) * fYlmax; //Gaussian distr.
345  fVolEff = 2 * TMath::Pi() * fTau * dYl * (fR * fR) / TMath::Power((fUmax), 2) *
346  ((fUmax)*TMath::SinH((fUmax)) - TMath::CosH((fUmax)) + 1);
347  LogInfo("Hydjet2Hadronizer|Param") << "central Effective volume = " << fVolEff << " [fm^3]";
348 
349  double particleDensity_pi_ch = 0;
350  double particleDensity_pi_th = 0;
351  // double particleDensity_th_0=0;
352 
353  if (fThFO != fT && fThFO > 0) {
354  GrandCanonical gc_ch(15, fT, fMuB, fMuS, fMuI3, fMuC);
355  GrandCanonical gc_pi_th(15, fThFO, 0., 0., fMu_th_pip, fMuC);
356  GrandCanonical gc_th_0(15, fThFO, 0., 0., 0., 0.);
357  particleDensity_pi_ch = gc_ch.ParticleNumberDensity(fDatabase->GetPDGParticle(211));
358  particleDensity_pi_th = gc_pi_th.ParticleNumberDensity(fDatabase->GetPDGParticle(211));
359  }
360 
361  for (int particleIndex = 0; particleIndex < fDatabase->GetNParticles(); particleIndex++) {
362  ParticlePDG* currParticle = fDatabase->GetPDGParticleByIndex(particleIndex);
363  int encoding = currParticle->GetPDG();
364 
365  //strangeness supression
366  double gammaS = 1;
367  int S = int(currParticle->GetStrangeness());
368  if (encoding == 333)
369  S = 2;
370  if (fCorrS < 1. && S != 0)
371  gammaS = TMath::Power(fCorrS, -TMath::Abs(S));
372 
373  //average densities
374  double particleDensity = gc.ParticleNumberDensity(currParticle) / gammaS;
375 
376  //compute chemical potential for single f.o. mu==mu_ch
377  double mu = fMuB * int(currParticle->GetBaryonNumber()) + fMuS * int(currParticle->GetStrangeness()) +
378  fMuI3 * int(currParticle->GetElectricCharge()) + fMuC * int(currParticle->GetCharmness());
379 
380  //thermal f.o.
381  if (fThFO != fT && fThFO > 0) {
382  double particleDensity_ch = gc_ch.ParticleNumberDensity(currParticle);
383  double particleDensity_th_0 = gc_th_0.ParticleNumberDensity(currParticle);
384  double numb_dens_bolt = particleDensity_pi_th * particleDensity_ch / particleDensity_pi_ch;
385  mu = fThFO * TMath::Log(numb_dens_bolt / particleDensity_th_0);
386  if (abs(encoding) == 211 || encoding == 111)
387  mu = fMu_th_pip;
388  particleDensity = numb_dens_bolt;
389  }
390 
391  // set particle densities to zero for some particle codes
392  // pythia quark codes
393  if (abs(encoding) <= 9) {
394  particleDensity = 0;
395  }
396  // leptons
397  if (abs(encoding) > 10 && abs(encoding) < 19) {
398  particleDensity = 0;
399  }
400  // exchange bosons
401  if (abs(encoding) > 20 && abs(encoding) < 30) {
402  particleDensity = 0;
403  }
404  // pythia special codes (e.g. strings, clusters ...)
405  if (abs(encoding) > 80 && abs(encoding) < 100) {
406  particleDensity = 0;
407  }
408  // pythia di-quark codes
409  // Note: in PYTHIA all diquark codes have the tens digits equal to zero
410  if (abs(encoding) > 1000 && abs(encoding) < 6000) {
411  int tens = ((abs(encoding) - (abs(encoding) % 10)) / 10) % 10;
412  if (tens == 0) { // its a diquark;
413  particleDensity = 0;
414  }
415  }
416  // K0S and K0L
417  if (abs(encoding) == 130 || abs(encoding) == 310) {
418  particleDensity = 0;
419  }
420  // charmed particles
421 
422  if (encoding == 443)
423  NJPsith = particleDensity * fVolEff / dYl;
424 
425  // We generate thermo-statistically only J/psi(443), D_+(411), D_-(-411), D_0(421),
426  //Dbar_0(-421), D1_+(413), D1_-(-413), D1_0(423), D1bar_0(-423)
427  //Dcs(431) Lambdac(4122)
428  if (currParticle->GetCharmQNumber() != 0 || currParticle->GetCharmAQNumber() != 0) {
429  //ml if(abs(encoding)!=443 &&
430  //ml abs(encoding)!=411 && abs(encoding)!=421 &&
431  //ml abs(encoding)!=413 && abs(encoding)!=423 && abs(encoding)!=4122 && abs(encoding)!=431) {
432  //ml particleDensity=0; }
433 
434  if (abs(encoding) == 441 || abs(encoding) == 10441 || abs(encoding) == 10443 || abs(encoding) == 20443 ||
435  abs(encoding) == 445 || abs(encoding) == 4232 || abs(encoding) == 4322 || abs(encoding) == 4132 ||
436  abs(encoding) == 4312 || abs(encoding) == 4324 || abs(encoding) == 4314 || abs(encoding) == 4332 ||
437  abs(encoding) == 4334) {
438  particleDensity = 0;
439  } else {
440  if (abs(encoding) != 443) { //only open charm
441  Nocth = Nocth + particleDensity * fVolEff / dYl;
442  LogInfo("Hydjet2Hadronizer|Charam") << encoding << " Nochth " << Nocth;
443  // particleDensity=particleDensity*fCorrC;
444  // if(abs(encoding)==443)particleDensity=particleDensity*fCorrC;
445  }
446  }
447  }
448 
449  // bottom mesons
450  if ((abs(encoding) > 500 && abs(encoding) < 600) || (abs(encoding) > 10500 && abs(encoding) < 10600) ||
451  (abs(encoding) > 20500 && abs(encoding) < 20600) || (abs(encoding) > 100500 && abs(encoding) < 100600)) {
452  particleDensity = 0;
453  }
454  // bottom baryons
455  if (abs(encoding) > 5000 && abs(encoding) < 6000) {
456  particleDensity = 0;
457  }
459 
460  if (particleDensity > 0.) {
461  fPartEnc[fNPartTypes] = encoding;
462  fPartMult[2 * fNPartTypes] = particleDensity;
463  fPartMu[2 * fNPartTypes] = mu;
464  ++fNPartTypes;
465  if (fNPartTypes > 1000)
466  LogError("Hydjet2Hadronizer") << "fNPartTypes is too large" << fNPartTypes;
467 
468  //outMult<<encoding<<" "<<particleDensity*fVolEff/dYl <<" "<<mu<<std::endl;
469  }
470  }
471 
472  //put open charm number and cc number in Params
473  fNocth = Nocth;
474  fNccth = NJPsith;
475 
476  return kTRUE;
477 }
const double Pi
void SetBaryonPotential(double value)
ParticlePDG * GetPDGParticle(int pdg)
Definition: DatabasePDG.cc:218
#define HYPYIN
Definition: HYJET_COMMONS.h:84
int GetNParticles(bool all=kFALSE)
Definition: DatabasePDG.cc:582
#define HYJPAR
Definition: HYJET_COMMONS.h:69
#define PYQPAR
int GetPDG()
Definition: ParticlePDG.h:69
double CalculateStrangePotential()
Log< level::Error, false > LogError
double GetCharmAQNumber()
Definition: ParticlePDG.h:81
double nuclear_radius() const
Pythia6Service * pythia6Service_
ParticlePDG * GetPDGParticleByIndex(int index)
Definition: DatabasePDG.cc:198
double GetCharmness()
Definition: ParticlePDG.h:88
double GetElectricCharge()
Definition: ParticlePDG.h:89
T Abs(T a)
Definition: MathUtil.h:49
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
const int mu
Definition: Constants.h:22
double GetBaryonNumber()
Definition: ParticlePDG.h:82
double GetStrangeness()
Definition: ParticlePDG.h:87
#define HYIPAR
Definition: HYJET_COMMONS.h:26
Log< level::Info, false > LogInfo
void myini_()
double GetCharmQNumber()
Definition: ParticlePDG.h:80
DatabasePDG * fDatabase
Definition: InitialState.h:16
void SetTemperature(double value)
double Hydjet2Hadronizer::MidpointIntegrator2 ( double  a,
double  b,
double  Delta,
double  Epsilon 
)

Definition at line 1026 of file Hydjet2Hadronizer.cc.

References a, alignCSCRings::e, HLT_FULL_cff::Epsilon, f2(), fR, h, mps_fire::i, dqmiolumiharvest::j, and submitPVValidationJobs::t.

1026  {
1027  int nsubIntervals = 2000;
1028  int nsubIntervals2 = 1;
1029  double h = (b - a) / nsubIntervals; //0-pi , phi
1030  double h2 = (fR) / nsubIntervals; //0-R maximal RB ?
1031 
1032  double x = a + 0.5 * h;
1033  double y = 0;
1034 
1035  double t = f2(x, y, Delta);
1036 
1037  double e = Epsilon;
1038 
1039  for (int j = 1; j < nsubIntervals; j++) {
1040  x += h; // integr phi
1041 
1042  double RsB = fR; //test: podstavit' *coefff_RB
1043  double RB = RsB * (TMath::Sqrt(1 - e * e) / TMath::Sqrt(1 + e * TMath::Cos(2 * x))); //f-la7 RB
1044 
1045  nsubIntervals2 = int(RB / h2) + 1;
1046  // integr R
1047  y = 0;
1048  for (int i = 1; i < nsubIntervals2; i++)
1049  t += f2(x, (y += h2), Delta);
1050  }
1051  return t * h * h2;
1052 }
double f2(double, double, double)
double b
Definition: hdecay.h:118
double a
Definition: hdecay.h:119
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
double Hydjet2Hadronizer::nuclear_radius ( ) const
inlineprivate

Definition at line 242 of file Hydjet2Hadronizer.h.

References fAw, and funct::pow().

Referenced by add_heavy_ion_rec(), and initializeForInternalPartons().

242  {
243  // Return the nuclear radius derived from the
244  // beam/target atomic mass number.
245  return 1.15 * pow((double)fAw, 1. / 3.);
246  }
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
bool Hydjet2Hadronizer::readSettings ( int  )

Definition at line 233 of file Hydjet2Hadronizer.cc.

References fNPartTypes, hjRandomEngine, pythia6Service_, SERVICE, and gen::Pythia6Service::setGeneralParams().

233  {
234  Pythia6Service::InstanceWrapper guard(pythia6Service_);
236 
237  SERVICE.iseed_fromC = hjRandomEngine->CLHEP::HepRandomEngine::getSeed();
238  LogInfo("Hydjet2Hadronizer|GenSeed") << "Seed for random number generation: "
239  << hjRandomEngine->CLHEP::HepRandomEngine::getSeed();
240 
241  fNPartTypes = 0; //counter of hadron species
242 
243  return kTRUE;
244 }
Pythia6Service * pythia6Service_
CLHEP::HepRandomEngine * hjRandomEngine
#define SERVICE
Definition: HYJET_COMMONS.h:38
Log< level::Info, false > LogInfo
bool Hydjet2Hadronizer::residualDecay ( )

Definition at line 965 of file Hydjet2Hadronizer.cc.

965 { return true; }
void Hydjet2Hadronizer::rotateEvtPlane ( )
private

Definition at line 1085 of file Hydjet2Hadronizer.cc.

References funct::cos(), cosphi0_, phi0_, pi, gen::pyr_(), funct::sin(), and sinphi0_.

Referenced by generatePartonsAndHadronize().

1085  {
1086  const double pi = 3.14159265358979;
1087  phi0_ = 2. * pi * gen::pyr_(nullptr) - pi;
1088  sinphi0_ = sin(phi0_);
1089  cosphi0_ = cos(phi0_);
1090 }
double pyr_(int *idummy)
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
const Double_t pi
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
bool gen::Hydjet2Hadronizer::RunDecays ( )
inlineoverridevirtual

Implements InitialState.

Definition at line 88 of file Hydjet2Hadronizer.h.

References fDecay.

Referenced by generatePartonsAndHadronize().

88 { return (fDecay > 0 ? kTRUE : kFALSE); }
void gen::Hydjet2Hadronizer::SetVolEff ( double  value)
inline

Definition at line 86 of file Hydjet2Hadronizer.h.

References fVolEff, and relativeConstraints::value.

double Hydjet2Hadronizer::SimpsonIntegrator ( double  a,
double  b,
double  phi,
double  Delta 
)

Definition at line 986 of file Hydjet2Hadronizer.cc.

References a, f2(), h, mps_fire::i, LogDebug, alignCSCRings::s, and submitPVValidationJobs::t.

Referenced by SimpsonIntegrator2().

986  {
987  LogDebug("SimpsonIntegrator") << "in SimpsonIntegrator"
988  << "delta - " << Delta;
989  int nsubIntervals = 100;
990  double h = (b - a) / nsubIntervals;
991  double s = f2(phi, a + 0.5 * h, Delta);
992  double t = 0.5 * (f2(phi, a, Delta) + f2(phi, b, Delta));
993  double x = a;
994  double y = a + 0.5 * h;
995  for (int i = 1; i < nsubIntervals; i++) {
996  x += h;
997  y += h;
998  s += f2(phi, y, Delta);
999  t += f2(phi, x, Delta);
1000  }
1001  t += 2.0 * s;
1002  return t * h / 3.0;
1003 }
double f2(double, double, double)
double b
Definition: hdecay.h:118
double a
Definition: hdecay.h:119
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
#define LogDebug(id)
double Hydjet2Hadronizer::SimpsonIntegrator2 ( double  a,
double  b,
double  Epsilon,
double  Delta 
)

Definition at line 1006 of file Hydjet2Hadronizer.cc.

References a, alignCSCRings::e, HLT_FULL_cff::Epsilon, fR, h, dqmiolumiharvest::j, alignCSCRings::s, and SimpsonIntegrator().

Referenced by generatePartonsAndHadronize().

1006  {
1007  LogInfo("SimpsonIntegrator2") << "in SimpsonIntegrator2: epsilon - " << Epsilon << " delta - " << Delta;
1008  int nsubIntervals = 10000;
1009  double h = (b - a) / nsubIntervals; //-1-pi, phi
1010  double s = 0;
1011  // double h2 = (fR)/nsubIntervals; //0-R maximal RB ?
1012 
1013  double x = 0; //phi
1014  for (int j = 1; j < nsubIntervals; j++) {
1015  x += h; // phi
1016  double e = Epsilon;
1017  double RsB = fR; //test: podstavit' *coefff_RB
1018  double RB = RsB * (TMath::Sqrt(1 - e * e) / TMath::Sqrt(1 + e * TMath::Cos(2 * x))); //f-la7 RB
1019  double sr = SimpsonIntegrator(0, RB, x, Delta);
1020  s += sr;
1021  }
1022  return s * h;
1023 }
double SimpsonIntegrator(double, double, double, double)
Log< level::Info, false > LogInfo
double b
Definition: hdecay.h:118
double a
Definition: hdecay.h:119
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
void Hydjet2Hadronizer::statistics ( )

Definition at line 967 of file Hydjet2Hadronizer.cc.

967 {}

Member Data Documentation

ParticleAllocator gen::Hydjet2Hadronizer::allocator
private

Definition at line 239 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

float gen::Hydjet2Hadronizer::Bgen
private

Definition at line 219 of file Hydjet2Hadronizer.h.

Referenced by add_heavy_ion_rec(), and generatePartonsAndHadronize().

double gen::Hydjet2Hadronizer::cosphi0_
private
float gen::Hydjet2Hadronizer::E[500000]
private

Definition at line 223 of file Hydjet2Hadronizer.h.

Referenced by build_hyjet2(), and generatePartonsAndHadronize().

bool gen::Hydjet2Hadronizer::embedding_
private

Definition at line 188 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize(), and Hydjet2Hadronizer().

int gen::Hydjet2Hadronizer::ev
private

Definition at line 217 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

HepMC::GenEvent* gen::Hydjet2Hadronizer::evt
private

Definition at line 190 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

double gen::Hydjet2Hadronizer::fAw
private

Definition at line 120 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons(), and nuclear_radius().

double gen::Hydjet2Hadronizer::fBfix
private
double gen::Hydjet2Hadronizer::fBmax
private

Definition at line 125 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

double gen::Hydjet2Hadronizer::fBmin
private

Definition at line 124 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

int gen::Hydjet2Hadronizer::fCharmProd
private

Definition at line 158 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

double gen::Hydjet2Hadronizer::fCorrC
private

Definition at line 159 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

double gen::Hydjet2Hadronizer::fCorrS
private

Definition at line 157 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

int gen::Hydjet2Hadronizer::fDecay
private

Definition at line 145 of file Hydjet2Hadronizer.h.

Referenced by RunDecays().

double gen::Hydjet2Hadronizer::fDelta
private

Definition at line 141 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

double gen::Hydjet2Hadronizer::fEpsilon
private

Definition at line 142 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

int gen::Hydjet2Hadronizer::fEtaType
private
int gen::Hydjet2Hadronizer::fIanglu
private

Definition at line 185 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

int gen::Hydjet2Hadronizer::fIenglu
private

Definition at line 182 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

int gen::Hydjet2Hadronizer::fIfb
private
int gen::Hydjet2Hadronizer::fIfDeltaEpsilon
private

Definition at line 143 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

int gen::Hydjet2Hadronizer::final[500000]
private

Definition at line 237 of file Hydjet2Hadronizer.h.

int gen::Hydjet2Hadronizer::FirstDaughterIndex[500000]
private

Definition at line 235 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

int gen::Hydjet2Hadronizer::fIshad
private

Definition at line 167 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

double gen::Hydjet2Hadronizer::fMu_th_pip
private

Definition at line 134 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

double gen::Hydjet2Hadronizer::fMuB
private

Definition at line 129 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

double gen::Hydjet2Hadronizer::fMuC
private

Definition at line 131 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

double gen::Hydjet2Hadronizer::fMuI3
private

Definition at line 132 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

double gen::Hydjet2Hadronizer::fMuS
private

Definition at line 130 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

double gen::Hydjet2Hadronizer::fMuTh[1000]
private

Definition at line 208 of file Hydjet2Hadronizer.h.

double gen::Hydjet2Hadronizer::fNccth
private
int gen::Hydjet2Hadronizer::fNf
private

Definition at line 181 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

int gen::Hydjet2Hadronizer::fNhsel
private
double gen::Hydjet2Hadronizer::fNocth
private
int gen::Hydjet2Hadronizer::fNPartTypes
private
int gen::Hydjet2Hadronizer::fPartEnc[1000]
private
double gen::Hydjet2Hadronizer::fPartMu[2000]
private
double gen::Hydjet2Hadronizer::fPartMult[2000]
private
double gen::Hydjet2Hadronizer::fPtmin
private

Definition at line 171 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

int gen::Hydjet2Hadronizer::fPyhist
private

Definition at line 166 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

int gen::Hydjet2Hadronizer::fPythDecay
private

Definition at line 147 of file Hydjet2Hadronizer.h.

double gen::Hydjet2Hadronizer::fR
private
edm::Service<TFileService> gen::Hydjet2Hadronizer::fs
private

Definition at line 215 of file Hydjet2Hadronizer.h.

double gen::Hydjet2Hadronizer::fSigmaTau
private

Definition at line 137 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

double gen::Hydjet2Hadronizer::fSqrtS
private

Definition at line 119 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

double gen::Hydjet2Hadronizer::fT
private
double gen::Hydjet2Hadronizer::fT0
private

Definition at line 176 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

double gen::Hydjet2Hadronizer::fTau
private
double gen::Hydjet2Hadronizer::fTau0
private

Definition at line 180 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

double gen::Hydjet2Hadronizer::fThFO
private
int gen::Hydjet2Hadronizer::fTMuType
private

Definition at line 154 of file Hydjet2Hadronizer.h.

Referenced by initializeForInternalPartons().

double gen::Hydjet2Hadronizer::fUmax
private
double gen::Hydjet2Hadronizer::fVolEff
private
double gen::Hydjet2Hadronizer::fWeakDecay
private

Definition at line 146 of file Hydjet2Hadronizer.h.

Referenced by GetWeakDecayLimit().

double gen::Hydjet2Hadronizer::fYlmax
private
int gen::Hydjet2Hadronizer::Index[500000]
private

Definition at line 232 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

int gen::Hydjet2Hadronizer::LastDaughterIndex[500000]
private

Definition at line 236 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

unsigned int gen::Hydjet2Hadronizer::maxEventsToPrint_
private

Definition at line 200 of file Hydjet2Hadronizer.h.

Referenced by Hydjet2Hadronizer().

int gen::Hydjet2Hadronizer::MotherIndex[500000]
private

Definition at line 233 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize(), and get_particles().

int gen::Hydjet2Hadronizer::Mpdg[500000]
private

Definition at line 229 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

int gen::Hydjet2Hadronizer::Nbcol
private

Definition at line 217 of file Hydjet2Hadronizer.h.

Referenced by add_heavy_ion_rec(), and generatePartonsAndHadronize().

int gen::Hydjet2Hadronizer::NDaughters[500000]
private

Definition at line 234 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

int gen::Hydjet2Hadronizer::nhard_
private

Definition at line 192 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize(), and get_particles().

int gen::Hydjet2Hadronizer::Nhyd
private

Definition at line 217 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

int gen::Hydjet2Hadronizer::Njet
private

Definition at line 217 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize(), and get_particles().

int gen::Hydjet2Hadronizer::Npart
private

Definition at line 217 of file Hydjet2Hadronizer.h.

Referenced by add_heavy_ion_rec(), and generatePartonsAndHadronize().

int gen::Hydjet2Hadronizer::Npyt
private

Definition at line 217 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

int gen::Hydjet2Hadronizer::nsoft_
private

Definition at line 193 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize(), and get_particles().

int gen::Hydjet2Hadronizer::nsub_
private
int gen::Hydjet2Hadronizer::Ntot
private

Definition at line 217 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

int gen::Hydjet2Hadronizer::pdg[500000]
private
double gen::Hydjet2Hadronizer::phi0_
private
edm::ParameterSet gen::Hydjet2Hadronizer::pset
private

Definition at line 214 of file Hydjet2Hadronizer.h.

double gen::Hydjet2Hadronizer::psiforv3
private

Definition at line 218 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

float gen::Hydjet2Hadronizer::Px[500000]
private

Definition at line 220 of file Hydjet2Hadronizer.h.

Referenced by build_hyjet2(), and generatePartonsAndHadronize().

float gen::Hydjet2Hadronizer::Py[500000]
private

Definition at line 221 of file Hydjet2Hadronizer.h.

Referenced by build_hyjet2(), and generatePartonsAndHadronize().

Pythia6Service* gen::Hydjet2Hadronizer::pythia6Service_
private
unsigned int gen::Hydjet2Hadronizer::pythiaPylistVerbosity_
private

Definition at line 199 of file Hydjet2Hadronizer.h.

Referenced by Hydjet2Hadronizer().

int gen::Hydjet2Hadronizer::pythiaStatus[500000]
private

Definition at line 231 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

float gen::Hydjet2Hadronizer::Pz[500000]
private

Definition at line 222 of file Hydjet2Hadronizer.h.

Referenced by build_hyjet2(), and generatePartonsAndHadronize().

bool gen::Hydjet2Hadronizer::rotate_
private

Definition at line 189 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

float gen::Hydjet2Hadronizer::Sigin
private

Definition at line 219 of file Hydjet2Hadronizer.h.

Referenced by add_heavy_ion_rec(), and generatePartonsAndHadronize().

float gen::Hydjet2Hadronizer::Sigjet
private

Definition at line 219 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize().

double gen::Hydjet2Hadronizer::sinphi0_
private
List_t gen::Hydjet2Hadronizer::source
private
edm::InputTag gen::Hydjet2Hadronizer::src_
private

Definition at line 202 of file Hydjet2Hadronizer.h.

Referenced by generatePartonsAndHadronize(), and Hydjet2Hadronizer().

int gen::Hydjet2Hadronizer::sseed
private

Definition at line 217 of file Hydjet2Hadronizer.h.

float gen::Hydjet2Hadronizer::T[500000]
private

Definition at line 227 of file Hydjet2Hadronizer.h.

const std::vector< std::string > Hydjet2Hadronizer::theSharedResources
staticprivate
int gen::Hydjet2Hadronizer::type[500000]
private
float gen::Hydjet2Hadronizer::X[500000]
private
float gen::Hydjet2Hadronizer::Y[500000]
private
float gen::Hydjet2Hadronizer::Z[500000]
private