CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HydjetHadronizer.cc
Go to the documentation of this file.
1 /*
2  * $Id: HydjetHadronizer.cc,v 1.14 2013/05/23 14:40:17 gartung Exp $
3  *
4  * Interface to the HYDJET generator, produces HepMC events
5  *
6  * Original Author: Camelia Mironov
7  */
8 
9 #include <iostream>
10 #include <cmath>
11 
12 #include "boost/lexical_cast.hpp"
13 
22 
27 
28 #include "HepMC/PythiaWrapper6_4.h"
29 #include "HepMC/GenEvent.h"
30 #include "HepMC/HeavyIon.h"
31 #include "HepMC/SimpleVector.h"
32 
37 
38 using namespace edm;
39 using namespace std;
40 using namespace gen;
41 
42 namespace {
43  int convertStatus(int st){
44  if(st<= 0) return 0;
45  if(st<=10) return 1;
46  if(st<=20) return 2;
47  if(st<=30) return 3;
48  else return st;
49  }
50 }
51 
52 
53 //_____________________________________________________________________
54 HydjetHadronizer::HydjetHadronizer(const ParameterSet &pset) :
55  BaseHadronizer(pset),
56  evt(0),
57  pset_(pset),
58  abeamtarget_(pset.getParameter<double>("aBeamTarget")),
59  bfixed_(pset.getParameter<double>("bFixed")),
60  bmax_(pset.getParameter<double>("bMax")),
61  bmin_(pset.getParameter<double>("bMin")),
62  cflag_(pset.getParameter<int>("cFlag")),
63  embedding_(pset.getParameter<bool>("embeddingMode")),
64  comenergy(pset.getParameter<double>("comEnergy")),
65  doradiativeenloss_(pset.getParameter<bool>("doRadiativeEnLoss")),
66  docollisionalenloss_(pset.getParameter<bool>("doCollisionalEnLoss")),
67  fracsoftmult_(pset.getParameter<double>("fracSoftMultiplicity")),
68  hadfreeztemp_(pset.getParameter<double>("hadronFreezoutTemperature")),
69  hymode_(pset.getParameter<string>("hydjetMode")),
70  maxEventsToPrint_(pset.getUntrackedParameter<int>("maxEventsToPrint", 1)),
71  maxlongy_(pset.getParameter<double>("maxLongitudinalRapidity")),
72  maxtrany_(pset.getParameter<double>("maxTransverseRapidity")),
73  nsub_(0),
74  nhard_(0),
75  nmultiplicity_(pset.getParameter<int>("nMultiplicity")),
76  nsoft_(0),
77  nquarkflavor_(pset.getParameter<int>("qgpNumQuarkFlavor")),
78  pythiaPylistVerbosity_(pset.getUntrackedParameter<int>("pythiaPylistVerbosity", 0)),
79  qgpt0_(pset.getParameter<double>("qgpInitialTemperature")),
80  qgptau0_(pset.getParameter<double>("qgpProperTimeFormation")),
81  phi0_(0.),
82  sinphi0_(0.),
83  cosphi0_(1.),
84  rotate_(pset.getParameter<bool>("rotateEventPlane")),
85  shadowingswitch_(pset.getParameter<int>("shadowingSwitch")),
86  signn_(pset.getParameter<double>("sigmaInelNN")),
87  pythia6Service_(new Pythia6Service(pset))
88 {
89  // Default constructor
90 
91  // PYLIST Verbosity Level
92  // Valid PYLIST arguments are: 1, 2, 3, 5, 7, 11, 12, 13
93  pythiaPylistVerbosity_ = pset.getUntrackedParameter<int>("pythiaPylistVerbosity",0);
94  LogDebug("PYLISTverbosity") << "Pythia PYLIST verbosity level = " << pythiaPylistVerbosity_;
95 
96  //Max number of events printed on verbosity level
97  maxEventsToPrint_ = pset.getUntrackedParameter<int>("maxEventsToPrint",0);
98  LogDebug("Events2Print") << "Number of events to be printed = " << maxEventsToPrint_;
99 
100  if(embedding_) src_ = pset.getParameter<edm::InputTag>("backgroundLabel");
101 
102 }
103 
104 
105 //_____________________________________________________________________
107 {
108  // destructor
109  call_pystat(1);
110  delete pythia6Service_;
111 }
112 
113 
114 //_____________________________________________________________________
115 void HydjetHadronizer::add_heavy_ion_rec(HepMC::GenEvent *evt)
116 {
117  // heavy ion record in the final CMSSW Event
118  double npart = hyfpar.npart;
119  int nproj = static_cast<int>(npart / 2);
120  int ntarg = static_cast<int>(npart - nproj);
121 
122  HepMC::HeavyIon* hi = new HepMC::HeavyIon(
123  nsub_, // Ncoll_hard/N of SubEvents
124  nproj, // Npart_proj
125  ntarg, // Npart_targ
126  static_cast<int>(hyfpar.nbcol), // Ncoll
127  0, // spectator_neutrons
128  0, // spectator_protons
129  0, // N_Nwounded_collisions
130  0, // Nwounded_N_collisions
131  0, // Nwounded_Nwounded_collisions
132  hyfpar.bgen * nuclear_radius(), // impact_parameter in [fm]
133  phi0_, // event_plane_angle
134  0, // eccentricity
135  hyjpar.sigin // sigma_inel_NN
136  );
137 
138  evt->set_heavy_ion(*hi);
139  delete hi;
140 }
141 
142 //___________________________________________________________________
144 {
145  // Build particle object corresponding to index in hyjets (soft+hard)
146 
147  double x0 = hyjets.phj[0][index];
148  double y0 = hyjets.phj[1][index];
149 
150  double x = x0*cosphi0_-y0*sinphi0_;
151  double y = y0*cosphi0_+x0*sinphi0_;
152 
154  HepMC::FourVector(x, // px
155  y, // py
156  hyjets.phj[2][index], // pz
157  hyjets.phj[3][index]), // E
158  hyjets.khj[1][index], // id
159  convertStatus(hyjets.khj[0][index] // status
160  )
161  );
162 
163  p->suggest_barcode(barcode);
164  return p;
165 }
166 
167 //___________________________________________________________________
168 HepMC::GenVertex* HydjetHadronizer::build_hyjet_vertex(int i,int id)
169 {
170  // build verteces for the hyjets stored events
171 
172  double x0=hyjets.vhj[0][i];
173  double y0=hyjets.vhj[1][i];
174  double x = x0*cosphi0_-y0*sinphi0_;
175  double y = y0*cosphi0_+x0*sinphi0_;
176  double z=hyjets.vhj[2][i];
177  double t=hyjets.vhj[4][i];
178 
179  HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(x,y,z,t),id);
180  return vertex;
181 }
182 
183 //___________________________________________________________________
184 
186 {
188 
189  // generate single event
190  if(embedding_){
191  cflag_ = 0;
192  const edm::Event& e = getEDMEvent();
194  e.getByLabel(src_,input);
195  const HepMC::GenEvent * inev = input->GetEvent();
196  const HepMC::HeavyIon* hi = inev->heavy_ion();
197  if(hi){
198  bfixed_ = hi->impact_parameter();
199  phi0_ = hi->event_plane_angle();
200  sinphi0_ = sin(phi0_);
201  cosphi0_ = cos(phi0_);
202  }else{
203  LogWarning("EventEmbedding")<<"Background event does not have heavy ion record!";
204  }
205  }else if(rotate_) rotateEvtPlane();
206 
207  nsoft_ = 0;
208  nhard_ = 0;
209 
210  edm::LogInfo("HYDJETmode") << "##### HYDJET nhsel = " << hyjpar.nhsel;
211  edm::LogInfo("HYDJETfpart") << "##### HYDJET fpart = " << hyflow.fpart;
212  edm::LogInfo("HYDJETtf") << "##### HYDJET hadron freez-out temp, Tf = " << hyflow.Tf;
213  edm::LogInfo("HYDJETinTemp") << "##### HYDJET: QGP init temperature, T0 ="<<pyqpar.T0u;
214  edm::LogInfo("HYDJETinTau") << "##### HYDJET: QGP formation time,tau0 ="<<pyqpar.tau0u;
215 
216  // generate a HYDJET event
217  int ntry = 0;
218  while(nsoft_ == 0 && nhard_ == 0){
219  if(ntry > 100){
220  edm::LogError("HydjetEmptyEvent") << "##### HYDJET: No Particles generated, Number of tries ="<<ntry;
221 
222  // Throw an exception. Use the EventCorruption exception since it maps onto SkipEvent
223  // which is what we want to do here.
224 
225  std::ostringstream sstr;
226  sstr << "HydjetHadronizerProducer: No particles generated after " << ntry << " tries.\n";
227  edm::Exception except(edm::errors::EventCorruption, sstr.str());
228  throw except;
229  } else {
230  HYEVNT();
231  nsoft_ = hyfpar.nhyd;
232  nsub_ = hyjpar.njet;
233  nhard_ = hyfpar.npyt;
234  ++ntry;
235  }
236  }
237 
238  if(hyjpar.nhsel < 3) nsub_++;
239 
240  // event information
241  HepMC::GenEvent *evt = new HepMC::GenEvent();
242 
243  if(nhard_>0 || nsoft_>0) get_particles(evt);
244 
245  evt->set_signal_process_id(pypars.msti[0]); // type of the process
246  evt->set_event_scale(pypars.pari[16]); // Q^2
247  add_heavy_ion_rec(evt);
248 
249  event().reset(evt);
250  return true;
251 }
252 
253 
254 //_____________________________________________________________________
255 bool HydjetHadronizer::get_particles(HepMC::GenEvent *evt )
256 {
257  // Hard particles. The first nhard_ lines from hyjets array.
258  // Pythia/Pyquen sub-events (sub-collisions) for a given event
259  // Return T/F if success/failure
260  // Create particles from lujet entries, assign them into vertices and
261  // put the vertices in the GenEvent, for each SubEvent
262  // The SubEvent information is kept by storing indeces of main vertices
263  // of subevents as a vector in GenHIEvent.
264 
265  LogDebug("SubEvent")<< "Number of sub events "<<nsub_;
266  LogDebug("Hydjet")<<"Number of hard events "<<hyjpar.njet;
267  LogDebug("Hydjet")<<"Number of hard particles "<<nhard_;
268  LogDebug("Hydjet")<<"Number of soft particles "<<nsoft_;
269 
270  vector<HepMC::GenVertex*> sub_vertices(nsub_);
271 
272  int ihy = 0;
273  for(int isub=0;isub<nsub_;isub++){
274  LogDebug("SubEvent") <<"Sub Event ID : "<<isub;
275 
276  int sub_up = (isub+1)*50000; // Upper limit in mother index, determining the range of Sub-Event
277  vector<HepMC::GenParticle*> particles;
278  vector<int> mother_ids;
279  vector<HepMC::GenVertex*> prods;
280 
281  sub_vertices[isub] = new HepMC::GenVertex(HepMC::FourVector(0,0,0,0),isub);
282  evt->add_vertex(sub_vertices[isub]);
283  if(!evt->signal_process_vertex()) evt->set_signal_process_vertex(sub_vertices[isub]);
284 
285  while(ihy<nhard_+nsoft_ && (hyjets.khj[2][ihy] < sub_up || ihy > nhard_ )){
286  particles.push_back(build_hyjet(ihy,ihy+1));
287  prods.push_back(build_hyjet_vertex(ihy,isub));
288  mother_ids.push_back(hyjets.khj[2][ihy]);
289  LogDebug("DecayChain")<<"Mother index : "<<hyjets.khj[2][ihy];
290 
291  ihy++;
292  }
293 
294  //Produce Vertices and add them to the GenEvent. Remember that GenParticles are adopted by
295  //GenVertex and GenVertex is adopted by GenEvent.
296 
297  LogDebug("Hydjet")<<"Number of particles in vector "<<particles.size();
298 
299  for (unsigned int i = 0; i<particles.size(); i++) {
300  HepMC::GenParticle* part = particles[i];
301 
302  //The Fortran code is modified to preserve mother id info, by seperating the beginning
303  //mother indices of successive subevents by 5000
304  int mid = mother_ids[i]-isub*50000-1;
305  LogDebug("DecayChain")<<"Particle "<<i;
306  LogDebug("DecayChain")<<"Mother's ID "<<mid;
307  LogDebug("DecayChain")<<"Particle's PDG ID "<<part->pdg_id();
308 
309  if(mid <= 0){
310  sub_vertices[isub]->add_particle_out(part);
311  continue;
312  }
313 
314  if(mid > 0){
315  HepMC::GenParticle* mother = particles[mid];
316  LogDebug("DecayChain")<<"Mother's PDG ID "<<mother->pdg_id();
317 
318  HepMC::GenVertex* prod_vertex = mother->end_vertex();
319  if(!prod_vertex){
320  prod_vertex = prods[i];
321  prod_vertex->add_particle_in(mother);
322  evt->add_vertex(prod_vertex);
323  prods[i]=0; // mark to protect deletion
324  }
325  prod_vertex->add_particle_out(part);
326  }
327  }
328  // cleanup vertices not assigned to evt
329  for (unsigned int i = 0; i<prods.size(); i++) {
330  if(prods[i]) delete prods[i];
331  }
332  }
333  return true;
334 }
335 
336 
337 //______________________________________________________________
338 bool HydjetHadronizer::call_hyinit(double energy,double a, int ifb, double bmin,
339  double bmax,double bfix,int nh)
340 {
341  // initialize hydjet
342 
343  pydatr.mrpy[2]=1;
344  HYINIT(energy,a,ifb,bmin,bmax,bfix,nh);
345  return true;
346 }
347 
348 
349 //______________________________________________________________
351 {
352  // set hydjet options
353 
354  // hydjet running mode mode
355  // kHydroOnly --- nhsel=0 jet production off (pure HYDRO event), nhsel=0
356  // kHydroJets --- nhsle=1 jet production on, jet quenching off (HYDRO+njet*PYTHIA events)
357  // kHydroQJet --- nhsel=2 jet production & jet quenching on (HYDRO+njet*PYQUEN events)
358  // kJetsOnly --- nhsel=3 jet production on, jet quenching off, HYDRO off (njet*PYTHIA events)
359  // kQJetsOnly --- nhsel=4 jet production & jet quenching on, HYDRO off (njet*PYQUEN events)
360 
361  if(hymode_ == "kHydroOnly") hyjpar.nhsel=0;
362  else if ( hymode_ == "kHydroJets") hyjpar.nhsel=1;
363  else if ( hymode_ == "kHydroQJets") hyjpar.nhsel=2;
364  else if ( hymode_ == "kJetsOnly") hyjpar.nhsel=3;
365  else if ( hymode_ == "kQJetsOnly") hyjpar.nhsel=4;
366  else hyjpar.nhsel=2;
367 
368  // fraction of soft hydro induced multiplicity
369  hyflow.fpart = fracsoftmult_;
370 
371  // hadron freez-out temperature
372  hyflow.Tf = hadfreeztemp_;
373 
374  // maximum longitudinal collective rapidity
375  hyflow.ylfl = maxlongy_;
376 
377  // maximum transverse collective rapidity
378  hyflow.ytfl = maxtrany_;
379 
380  // shadowing on=1, off=0
381  hyjpar.ishad = shadowingswitch_;
382 
383  // set inelastic nucleon-nucleon cross section
384  hyjpar.sigin = signn_;
385 
386  // number of active quark flavors in qgp
387  pyqpar.nfu = nquarkflavor_;
388 
389  // initial temperature of QGP
390  pyqpar.T0u = qgpt0_;
391 
392  // proper time of QGP formation
393  pyqpar.tau0u = qgptau0_;
394 
395  // type of medium induced partonic energy loss
397  edm::LogInfo("HydjetEnLoss") << "##### Radiative AND Collisional partonic energy loss ON ####";
398  pyqpar.ienglu = 0;
399  } else if ( doradiativeenloss_ ) {
400  edm::LogInfo("HydjetenLoss") << "##### Only RADIATIVE partonic energy loss ON ####";
401  pyqpar.ienglu = 1;
402  } else if ( docollisionalenloss_ ) {
403  edm::LogInfo("HydjetEnLoss") << "##### Only COLLISIONAL partonic energy loss ON ####";
404  pyqpar.ienglu = 2;
405  } else {
406  edm::LogInfo("HydjetEnLoss") << "##### Radiative AND Collisional partonic energy loss ON ####";
407  pyqpar.ienglu = 0;
408  }
409  return true;
410 }
411 
412 //_____________________________________________________________________
413 
415 
418 
419  return true;
420 
421 }
422 
423 //_____________________________________________________________________
424 
426 
428  // pythia6Service_->setGeneralParams();
429 
430  // the input impact parameter (bxx_) is in [fm]; transform in [fm/RA] for hydjet usage
431  const float ra = nuclear_radius();
432  LogInfo("RAScaling")<<"Nuclear radius(RA) = "<<ra;
433  bmin_ /= ra;
434  bmax_ /= ra;
435  bfixed_ /= ra;
436 
437  // hydjet running options
439  // initialize hydjet
440  LogInfo("HYDJETinAction") << "##### Calling HYINIT("<<comenergy<<","<<abeamtarget_<<","
441  <<cflag_<<","<<bmin_<<","<<bmax_<<","<<bfixed_<<","<<nmultiplicity_<<") ####";
443  return true;
444 }
445 
446 bool HydjetHadronizer::declareStableParticles(const std::vector<int>& _pdg )
447 {
448  std::vector<int> pdg = _pdg;
449  for ( size_t i=0; i < pdg.size(); i++ ) {
450  int pyCode = pycomp_( pdg[i] );
451  std::ostringstream pyCard ;
452  pyCard << "MDCY(" << pyCode << ",1)=0";
453  std::cout << pyCard.str() << std::endl;
454  call_pygive( pyCard.str() );
455  }
456  return true;
457 }
458 
459 //________________________________________________________________
461 {
462  const double pi = 3.14159265358979;
463  phi0_ = 2.*pi*gen::pyr_(0) - pi;
464  sinphi0_ = sin(phi0_);
465  cosphi0_ = cos(phi0_);
466 }
467 
468 
469 //________________________________________________________________
471 {
472  return false;
473 }
474 
476 {
477  return true;
478 }
479 
481 {
482  return true;
483 }
484 
486 {
487 }
488 
490 {
491 }
492 
493 const char* HydjetHadronizer::classname() const
494 {
495  return "gen::HydjetHadronizer";
496 }
#define LogDebug(id)
#define HYEVNT
Definition: HydjetWrapper.h:27
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
HepMC::GenEvent * evt
bool docollisionalenloss_
DEFAULT = true.
bool hydjet_init(const edm::ParameterSet &pset)
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
bool call_pygive(const std::string &line)
double npart
Definition: HydjetWrapper.h:45
void add_heavy_ion_rec(HepMC::GenEvent *evt)
#define hyfpar
Definition: HydjetWrapper.h:50
float float float z
std::auto_ptr< HepMC::GenEvent > & event()
#define HYINIT
Definition: HydjetWrapper.h:19
double p[5][pyjets_maxn]
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
const char * classname() const
unsigned int maxEventsToPrint_
def convertStatus
Definition: crabWrap.py:174
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:361
HepMC::GenParticle * build_hyjet(int index, int barcode)
Pythia6Service * pythia6Service_
#define hyjpar
Definition: HydjetWrapper.h:95
int pycomp_(int &)
#define pypars
double fracsoftmult_
DEFAULT = true.
HepMC::GenVertex * build_hyjet_vertex(int i, int id)
part
Definition: HCALResponse.h:20
#define hyjets
Definition: HydjetWrapper.h:80
double nuclear_radius() const
bool call_hyinit(double energy, double a, int ifb, double bmin, double bmax, double bfix, int nh)
unsigned int pythiaPylistVerbosity_
#define pyqpar
#define hyflow
Definition: HydjetWrapper.h:38
double a
Definition: hdecay.h:121
edm::ParameterSet pset_
unsigned int shadowingswitch_
tuple cout
Definition: gather_cfg.py:121
double pi
Definition: DDAxes.h:10
double pyr_(int *idummy)
bool declareStableParticles(const std::vector< int > &)
bool get_particles(HepMC::GenEvent *evt)
edm::Event & getEDMEvent() const