CMS 3D CMS Logo

HijingHadronizer.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <cmath>
3 
4 #include "boost/lexical_cast.hpp"
5 
12 
19 
20 #include "HepMC/GenEvent.h"
21 #include "HepMC/HeavyIon.h"
22 #include "HepMC/SimpleVector.h"
23 
24 #include "CLHEP/Random/RandomEngine.h"
25 
26 static const double pi = 3.14159265358979;
27 
28 using namespace edm;
29 using namespace std;
30 using namespace gen;
31 
32 static CLHEP::HepRandomEngine* hijRandomEngine;
33 
34 extern "C" {
35 float gen::hijran_(int* idummy) { return hijRandomEngine->flat(); }
36 }
37 
38 extern "C" {
39 float ran_(unsigned int* iseed) {
40  return hijRandomEngine->flat();
41  // return ranff_(iseed);
42  // return gen::pyr_(0);
43 }
44 }
45 
46 extern "C" {
47 float rlu_(unsigned int* iseed) {
48  return hijRandomEngine->flat();
49  // return ranff_(iseed);
50  // return gen::pyr_(0);
51 }
52 }
53 
54 const std::vector<std::string> HijingHadronizer::theSharedResources = {edm::SharedResourceNames::kPythia6};
55 
56 HijingHadronizer::HijingHadronizer(const ParameterSet& pset)
58  evt(nullptr),
59  pset_(pset),
60  bmax_(pset.getParameter<double>("bMax")),
61  bmin_(pset.getParameter<double>("bMin")),
62  efrm_(pset.getParameter<double>("comEnergy")),
63  frame_(pset.getParameter<string>("frame")),
64  proj_(pset.getParameter<string>("proj")),
65  targ_(pset.getParameter<string>("targ")),
66  iap_(pset.getParameter<int>("iap")),
67  izp_(pset.getParameter<int>("izp")),
68  iat_(pset.getParameter<int>("iat")),
69  izt_(pset.getParameter<int>("izt")),
70  phi0_(0.),
71  sinphi0_(0.),
72  cosphi0_(1.),
73  rotate_(pset.getParameter<bool>("rotateEventPlane")) {
74  // Default constructor
75 }
76 
77 //_____________________________________________________________________
79  // destructor
80 }
81 
82 //_____________________________________________________________________
83 void HijingHadronizer::doSetRandomEngine(CLHEP::HepRandomEngine* v) { hijRandomEngine = v; }
84 
85 //_____________________________________________________________________
87  // heavy ion record in the final CMSSW Event
88  HepMC::HeavyIon* hi = new HepMC::HeavyIon(himain1.jatt, // Ncoll_hard/N of SubEvents
89  himain1.np, // Npart_proj
90  himain1.nt, // Npart_targ
91  himain1.n0 + himain1.n01 + himain1.n10 + himain1.n11, // Ncoll
92  0, // spectator_neutrons
93  0, // spectator_protons
94  himain1.n01, // N_Nwounded_collisions
95  himain1.n10, // Nwounded_N_collisions
96  himain1.n11, // Nwounded_Nwounded_collisions
97  //gsfs Changed from 19 to 18 (Fortran counts from 1 , not 0)
98  hiparnt.hint1[18], // impact_parameter in [fm]
99  phi0_, // event_plane_angle
100  0, // eccentricity
101  //gsfs Changed from 12 to 11 (Fortran counts from 1 , not 0)
102  hiparnt.hint1[11] // sigma_inel_NN
103  );
104  evt->set_heavy_ion(*hi);
105  delete hi;
106 }
107 
108 //___________________________________________________________________
110  // Build particle object corresponding to index in hijing
111 
112  double x0 = himain2.patt[0][index];
113  double y0 = himain2.patt[1][index];
114 
115  double x = x0 * cosphi0_ - y0 * sinphi0_;
116  double y = y0 * cosphi0_ + x0 * sinphi0_;
117 
118  // Hijing gives V0's status=4, they need to have status=1 to be decayed in geant
119  // also change status=11 to status=2
120  if (himain2.katt[3][index] <= 10 && himain2.katt[3][index] > 0)
121  himain2.katt[3][index] = 1;
122  if (himain2.katt[3][index] <= 20 && himain2.katt[3][index] > 10)
123  himain2.katt[3][index] = 2;
124 
125  HepMC::GenParticle* p = new HepMC::GenParticle(HepMC::FourVector(x, // px
126  y, // py
127  himain2.patt[2][index], // pz
128  himain2.patt[3][index]), // E
129  himain2.katt[0][index], // id
130  himain2.katt[3][index] // status
131  );
132  p->suggest_barcode(barcode);
133 
134  return p;
135 }
136 
137 //___________________________________________________________________
138 HepMC::GenVertex* HijingHadronizer::build_hijing_vertex(int i, int id) {
139  // build verteces for the hijing stored events
140  double x0 = himain2.vatt[0][i];
141  double y0 = himain2.vatt[1][i];
142  double x = x0 * cosphi0_ - y0 * sinphi0_;
143  double y = y0 * cosphi0_ + x0 * sinphi0_;
144  double z = himain2.vatt[2][i];
145  double t = himain2.vatt[3][i];
146 
147  HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(x, y, z, t), id);
148  return vertex;
149 }
150 
152  // generate single event
153  if (rotate_)
154  rotateEvtPlane();
155 
156  // generate a HIJING event
157 
158  float f_bmin = bmin_;
159  float f_bmax = bmax_;
160  HIJING(frame_.data(), f_bmin, f_bmax, strlen(frame_.data()));
161 
162  // event information
165 
166  // evt->set_signal_process_id(pypars.msti[0]); // type of the process
167  // evt->set_event_scale(pypars.pari[16]); // Q^2
169 
170  event().reset(evt);
171 
172  return true;
173 }
174 
175 //_____________________________________________________________________
177  HepMC::GenVertex* vertice;
178 
179  vector<HepMC::GenParticle*> particles;
180  vector<int> mother_ids;
181  vector<HepMC::GenVertex*> prods;
182 
183  vertice = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), 0);
184  evt->add_vertex(vertice);
185  if (!evt->signal_process_vertex())
186  evt->set_signal_process_vertex(vertice);
187 
188  const unsigned int knumpart = himain1.natt;
189 
190  for (unsigned int ipart = 0; ipart < knumpart; ipart++) {
191  int mid = himain2.katt[2][ipart] - 1; // careful of fortan to c++ array index
192 
193  particles.push_back(build_hijing(ipart, ipart + 1));
194  prods.push_back(build_hijing_vertex(ipart, 0));
195  mother_ids.push_back(mid);
196  LogDebug("DecayChain") << "Mother index : " << mid;
197  }
198 
199  LogDebug("Hijing") << "Number of particles in vector " << particles.size();
200 
201  for (unsigned int ipart = 0; ipart < particles.size(); ipart++) {
203 
204  int mid = mother_ids[ipart];
205  LogDebug("DecayChain") << "Particle " << ipart;
206  LogDebug("DecayChain") << "Mother's ID " << mid;
207  LogDebug("DecayChain") << "Particle's PDG ID " << part->pdg_id();
208 
209  // remove zero pT particles from list, protection for fastJet against pt=0 jets
210  if (part->status() == 1 &&
211  sqrt(part->momentum().px() * part->momentum().px() + part->momentum().py() * part->momentum().py()) == 0)
212  continue;
213 
214  if (mid <= 0) {
215  vertice->add_particle_out(part);
216  continue;
217  }
218 
219  if (mid > 0) {
220  HepMC::GenParticle* mother = particles[mid];
221  LogDebug("DecayChain") << "Mother's PDG ID " << mother->pdg_id();
222  HepMC::GenVertex* prod_vertex = mother->end_vertex();
223  if (!prod_vertex) {
224  prod_vertex = prods[ipart];
225  prod_vertex->add_particle_in(mother);
226 
227  evt->add_vertex(prod_vertex);
228  prods[ipart] = nullptr; // mark to protect deletion
229  }
230  prod_vertex->add_particle_out(part);
231  }
232  }
233 
234  // cleanup vertices not assigned to evt
235  for (unsigned int i = 0; i < prods.size(); i++) {
236  if (prods[i])
237  delete prods[i];
238  }
239 
240  return true;
241 }
242 
243 //_____________________________________________________________________
245  double efrm, std::string frame, std::string proj, std::string targ, int iap, int izp, int iat, int izt) {
246  float ef = efrm;
247  // initialize hydjet
248  HIJSET(ef,
249  frame.data(),
250  proj.data(),
251  targ.data(),
252  iap,
253  izp,
254  iat,
255  izt,
256  strlen(frame.data()),
257  strlen(proj.data()),
258  strlen(targ.data()));
259  return true;
260 }
261 
262 //______________________________________________________________
264  //initialize pythia5
265 
266  // std::string dumstr = "";
267  // call_pygive(dumstr);
268 
269  // initialize hijing
270  LogInfo("HIJINGinAction") << "##### Calling HIJSET(" << efrm_ << "," << frame_ << "," << proj_ << "," << targ_ << ","
271  << iap_ << "," << izp_ << "," << iat_ << "," << izt_ << ") ####";
273 
274  return true;
275 }
276 
277 bool HijingHadronizer::declareStableParticles(const std::vector<int>& pdg) { return true; }
278 
279 //________________________________________________________________
281  phi0_ = 2. * pi * gen::hijran_(nullptr) - pi;
282  sinphi0_ = sin(phi0_);
283  cosphi0_ = cos(phi0_);
284 }
285 
286 //________________________________________________________________
287 bool HijingHadronizer::hadronize() { return false; }
288 
289 bool HijingHadronizer::decay() { return true; }
290 
291 bool HijingHadronizer::residualDecay() { return true; }
292 
294 
296 
297 const char* HijingHadronizer::classname() const { return "gen::HijingHadronizer"; }
gen::HijingHadronizer::call_hijset
bool call_hijset(double efrm, std::string frame, std::string proj, std::string targ, int iap, int izp, int iat, int izt)
Definition: HijingHadronizer.cc:244
gen::HijingHadronizer::rotateEvtPlane
void rotateEvtPlane()
Definition: HijingHadronizer.cc:280
amptDefault_cfi.iap
iap
Definition: amptDefault_cfi.py:15
electrons_cff.bool
bool
Definition: electrons_cff.py:393
mps_fire.i
i
Definition: mps_fire.py:428
gen::HijingHadronizer::sinphi0_
double sinphi0_
Definition: HijingHadronizer.h:78
MessageLogger.h
gen::HijingHadronizer::declareStableParticles
bool declareStableParticles(const std::vector< int > &)
Definition: HijingHadronizer.cc:277
gen::HijingHadronizer::build_hijing
HepMC::GenParticle * build_hijing(int index, int barcode)
Definition: HijingHadronizer.cc:109
gen::HijingHadronizer::~HijingHadronizer
~HijingHadronizer() override
Definition: HijingHadronizer.cc:78
gen::HijingHadronizer::izp_
int izp_
Definition: HijingHadronizer.h:70
edm
HLT enums.
Definition: AlignableModifier.h:19
gen::HijingHadronizer::initializeForInternalPartons
bool initializeForInternalPartons()
Definition: HijingHadronizer.cc:263
gen::HijingHadronizer::residualDecay
bool residualDecay()
Definition: HijingHadronizer.cc:291
gen::HijingHadronizer::evt
HepMC::GenEvent * evt
Definition: HijingHadronizer.h:59
amptDefault_cfi.izt
izt
Definition: amptDefault_cfi.py:18
gen::HijingHadronizer::add_heavy_ion_rec
void add_heavy_ion_rec(HepMC::GenEvent *evt)
Definition: HijingHadronizer.cc:86
pi
static const double pi
Definition: HijingHadronizer.cc:26
gen::HijingHadronizer::izt_
int izt_
Definition: HijingHadronizer.h:72
HijingHadronizer.h
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:125
gen::BaseHadronizer
Definition: BaseHadronizer.h:46
GenRunInfoProduct.h
ecalTrigSettings_cff.particles
particles
Definition: ecalTrigSettings_cff.py:11
gen::HijingHadronizer::decay
bool decay()
Definition: HijingHadronizer.cc:289
HepMC::GenEvent
Definition: hepmc_rootio.cc:9
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
rlu_
float rlu_(unsigned int *iseed)
Definition: HijingHadronizer.cc:47
EDMException.h
part
part
Definition: HCALResponse.h:20
HIJSET
#define HIJSET
Definition: HijingWrapper.h:23
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
gen::p
double p[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:76
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
Run.h
gen::HijingHadronizer::phi0_
double phi0_
Definition: HijingHadronizer.h:77
SharedResourceNames.h
gen
Definition: PythiaDecays.h:13
HijingWrapper.h
gen::HijingHadronizer::targ_
std::string targ_
Definition: HijingHadronizer.h:68
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
gen::HijingHadronizer::frame_
std::string frame_
Definition: HijingHadronizer.h:66
gen::HijingHadronizer::get_particles
bool get_particles(HepMC::GenEvent *evt)
Definition: HijingHadronizer.cc:176
amptDefault_cfi.proj
proj
Definition: amptDefault_cfi.py:13
himain2
#define himain2
Definition: HijingWrapper.h:52
bphysicsOniaDQM_cfi.vertex
vertex
Definition: bphysicsOniaDQM_cfi.py:7
HIJING
#define HIJING
Definition: HijingWrapper.h:28
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:223
edm::ParameterSet
Definition: ParameterSet.h:47
iseed
int iseed
Definition: AMPTWrapper.h:134
GenEventInfoProduct.h
Event.h
gen::v
double v[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:76
createfilelist.int
int
Definition: createfilelist.py:10
gen::hijran_
float hijran_(int *)
Definition: HijingHadronizer.cc:35
gen::HijingHadronizer::build_hijing_vertex
HepMC::GenVertex * build_hijing_vertex(int i, int id)
Definition: HijingHadronizer.cc:138
hijRandomEngine
static CLHEP::HepRandomEngine * hijRandomEngine
Definition: HijingHadronizer.cc:32
amptDefault_cfi.iat
iat
Definition: amptDefault_cfi.py:17
gen::HijingHadronizer::iat_
int iat_
Definition: HijingHadronizer.h:71
hiparnt
#define hiparnt
Definition: HijingWrapper.h:62
gen::HijingHadronizer::finalizeEvent
void finalizeEvent()
Definition: HijingHadronizer.cc:293
gen::HijingHadronizer::iap_
int iap_
Definition: HijingHadronizer.h:69
gen::HijingHadronizer::cosphi0_
double cosphi0_
Definition: HijingHadronizer.h:79
hi
Definition: EPCuts.h:4
GenParticle.GenParticle
GenParticle
Definition: GenParticle.py:18
std
Definition: JetResolutionObject.h:76
gen::BaseHadronizer::event
std::unique_ptr< HepMC::GenEvent > & event()
Definition: BaseHadronizer.h:86
amptDefault_cfi.frame
frame
Definition: amptDefault_cfi.py:12
edm::SharedResourceNames::kPythia6
static const std::string kPythia6
Definition: SharedResourceNames.h:26
gen::HijingHadronizer::efrm_
double efrm_
Definition: HijingHadronizer.h:65
gen::HijingHadronizer::classname
const char * classname() const
Definition: HijingHadronizer.cc:297
gen::HijingHadronizer::bmin_
double bmin_
Definition: HijingHadronizer.h:63
pdg
Definition: pdg_functions.h:28
ran_
float ran_(unsigned int *iseed)
Definition: HijingHadronizer.cc:39
AlignmentPI::index
index
Definition: AlignmentPayloadInspectorHelper.h:46
gen::HijingHadronizer::bmax_
double bmax_
Definition: HijingHadronizer.h:61
amptDefault_cfi.izp
izp
Definition: amptDefault_cfi.py:16
gen::HijingHadronizer::proj_
std::string proj_
Definition: HijingHadronizer.h:67
gen::HijingHadronizer::hadronize
bool hadronize()
Definition: HijingHadronizer.cc:287
ParameterSet.h
HepMCProduct.h
himain1
#define himain1
Definition: HijingWrapper.h:43
amptDefault_cfi.targ
targ
Definition: amptDefault_cfi.py:14
gen::HijingHadronizer::statistics
void statistics()
Definition: HijingHadronizer.cc:295
submitPVValidationJobs.t
string t
Definition: submitPVValidationJobs.py:644
gen::HijingHadronizer::doSetRandomEngine
void doSetRandomEngine(CLHEP::HepRandomEngine *v) override
Definition: HijingHadronizer.cc:83
gen::HijingHadronizer::generatePartonsAndHadronize
bool generatePartonsAndHadronize()
Definition: HijingHadronizer.cc:151
gen::HijingHadronizer::rotate_
bool rotate_
Definition: HijingHadronizer.h:80
HijingPythiaWrapper.h
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27