CMS 3D CMS Logo

AMPTHadronizer.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <cmath>
3 
4 #include "boost/lexical_cast.hpp"
5 
10 
13 
17 
18 #include "HepMC/GenEvent.h"
19 #include "HepMC/HeavyIon.h"
20 #include "HepMC/SimpleVector.h"
21 #include "CLHEP/Random/RandomEngine.h"
22 
23 static const double pi = 3.14159265358979;
24 
25 using namespace edm;
26 using namespace std;
27 using namespace gen;
28 
29 static CLHEP::HepRandomEngine* amptRandomEngine;
30 
31 extern "C"
32 {
33  float gen::ranart_(int *idummy)
34  {
35  if(false) idummy = idummy;
36  float rannum = amptRandomEngine->flat();
37  return rannum;
38  }
39 }
40 
41 extern "C"
42 {
43  float gen::ran1_(int *idummy)
44  {
45  if(false) idummy = idummy;
46  return amptRandomEngine->flat();
47  }
48 }
49 
50 AMPTHadronizer::AMPTHadronizer(const ParameterSet &pset) :
51  BaseHadronizer(pset),
52  evt(nullptr),
53  pset_(pset),
54  bmax_(pset.getParameter<double>("bMax")),
55  bmin_(pset.getParameter<double>("bMin")),
56  efrm_(pset.getParameter<double>("comEnergy")),
57  frame_(pset.getParameter<string>("frame")),
58  proj_(pset.getParameter<string>("proj")),
59  targ_(pset.getParameter<string>("targ")),
60  iap_(pset.getParameter<int>("iap")),
61  izp_(pset.getParameter<int>("izp")),
62  iat_(pset.getParameter<int>("iat")),
63  izt_(pset.getParameter<int>("izt")),
64  amptmode_(pset.getParameter<int>("amptmode")),
65  ntmax_(pset.getParameter<int>("ntmax")),
66  dt_(pset.getParameter<double>("dt")),
67  stringFragA_(pset.getParameter<double>("stringFragA")),
68  stringFragB_(pset.getParameter<double>("stringFragB")),
69  popcornmode_(pset.getParameter<bool>("popcornmode")),
70  popcornpar_(pset.getParameter<double>("popcornpar")),
71  shadowingmode_(pset.getParameter<bool>("shadowingmode")),
72  quenchingmode_(pset.getParameter<bool>("quenchingmode")),
73  quenchingpar_(pset.getParameter<double>("quenchingpar")),
74  pthard_(pset.getParameter<double>("pthard")),
75  mu_(pset.getParameter<double>("mu")),
76  izpc_(pset.getParameter<int>("izpc")),
77  alpha_(pset.getParameter<double>("alpha")),
78  dpcoal_(pset.getParameter<double>("dpcoal")),
79  drcoal_(pset.getParameter<double>("drcoal")),
80  ks0decay_(pset.getParameter<bool>("ks0decay")),
81  phidecay_(pset.getParameter<bool>("phidecay")),
82  deuteronmode_(pset.getParameter<int>("deuteronmode")),
83  deuteronfactor_(pset.getParameter<int>("deuteronfactor")),
84  deuteronxsec_(pset.getParameter<int>("deuteronxsec")),
85  minijetpt_(pset.getParameter<double>("minijetpt")),
86  maxmiss_(pset.getParameter<int>("maxmiss")),
87  doInitialAndFinalRadiation_(pset.getParameter<int>("doInitialAndFinalRadiation")),
88  ktkick_(pset.getParameter<int>("ktkick")),
89  diquarkembedding_(pset.getParameter<int>("diquarkembedding")),
90  diquarkpx_(pset.getParameter<double>("diquarkpx")),
91  diquarkpy_(pset.getParameter<double>("diquarkpy")),
92  diquarkx_(pset.getParameter<double>("diquarkx")),
93  diquarky_(pset.getParameter<double>("diquarky")),
94  phi0_(0.),
95  sinphi0_(0.),
96  cosphi0_(1.),
97  rotate_(pset.getParameter<bool>("rotateEventPlane"))
98 {
99 }
100 
101 //_____________________________________________________________________
103 {
104 }
105 
106 //_____________________________________________________________________
107 void AMPTHadronizer::doSetRandomEngine(CLHEP::HepRandomEngine* v)
108 {
110 }
111 
112 //_____________________________________________________________________
114 {
115  // heavy ion record in the final CMSSW Event
116  HepMC::HeavyIon* hi = new HepMC::HeavyIon(
117  hmain1.jatt, // Ncoll_hard/N of SubEvents
118  hmain1.np, // Npart_proj
119  hmain1.nt, // Npart_targ
120  hmain1.n0+hmain1.n01+hmain1.n10+hmain1.n11, // Ncoll
121  0, // spectator_neutrons
122  0, // spectator_protons
123  hmain1.n01, // N_Nwounded_collisions
124  hmain1.n10, // Nwounded_N_collisions
125  hmain1.n11, // Nwounded_Nwounded_collisions
126  hparnt.hint1[18], // impact_parameter in [fm]
127  phi0_, // event_plane_angle
128  0, // eccentricity
129  hparnt.hint1[11] // sigma_inel_NN
130  );
131  evt->set_heavy_ion(*hi);
132  delete hi;
133 }
134 
135 //___________________________________________________________________
137 {
138  // Build particle object corresponding to index in ampt
139 
140  float px0 = hbt.plast[index][0];
141  float py0 = hbt.plast[index][1];
142  float pz0 = hbt.plast[index][2];
143  float m = hbt.plast[index][3];
144 
145  float px = px0*cosphi0_-py0*sinphi0_;
146  float py = py0*cosphi0_+px0*sinphi0_;
147  float pz = pz0;
148  float e = sqrt(px*px+py*py+pz*pz+m*m);
149  int status = 1;
150 
152  HepMC::FourVector(px,
153  py,
154  pz,
155  e),
156  INVFLV(hbt.lblast[index]),// id
157  status // status
158  );
159 
160  p->suggest_barcode(barcode);
161  return p;
162 }
163 
164 //___________________________________________________________________
165 HepMC::GenVertex* AMPTHadronizer::build_ampt_vertex(int i,int id)
166 {
167  // build verteces for the ampt stored events
168  HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(0,0,0,0),id);
169  return vertex;
170 }
171 //_____________________________________________________________________
173 {
174  // generate single event
175  if(rotate_) rotateEvtPlane();
176 
177  // generate a AMPT event
178  AMPT(frame_.data(), bmin_, bmax_, strlen(frame_.data()));
179 
180  // event information
181  HepMC::GenEvent *evt = new HepMC::GenEvent();
182  get_particles(evt);
183 
184  add_heavy_ion_rec(evt);
185 
186  event().reset(evt);
187 
188 
189  return true;
190 }
191 
192 //_____________________________________________________________________
193 bool AMPTHadronizer::get_particles(HepMC::GenEvent *evt )
194 {
195  HepMC::GenVertex* vertice;
196 
197  vector<HepMC::GenParticle*> particles;
198  vector<int> mother_ids;
199  vector<HepMC::GenVertex*> prods;
200 
201  vertice = new HepMC::GenVertex(HepMC::FourVector(0,0,0,0),0);
202  evt->add_vertex(vertice);
203  if(!evt->signal_process_vertex()) evt->set_signal_process_vertex(vertice);
204 
205  const unsigned int knumpart = hbt.nlast;
206  for (unsigned int ipart = 0; ipart<knumpart; ipart++) {
207  int mid = 0;
208  particles.push_back(build_ampt(ipart,ipart+1));
209  prods.push_back(build_ampt_vertex(ipart,0));
210  mother_ids.push_back(mid);
211  LogDebug("DecayChain")<<"Mother index : "<<mid;
212  }
213 
214  LogDebug("AMPT")<<"Number of particles in vector "<<particles.size();
215 
216  for (unsigned int ipart = 0; ipart<particles.size(); ipart++) {
217  HepMC::GenParticle* part = particles[ipart];
218 
219  int mid = mother_ids[ipart];
220  LogDebug("DecayChain")<<"Particle "<<ipart;
221  LogDebug("DecayChain")<<"Mother's ID "<<mid;
222  LogDebug("DecayChain")<<"Particle's PDG ID "<<part->pdg_id();
223 
224  if(mid <= 0){
225  vertice->add_particle_out(part);
226  continue;
227  }
228 
229  if(mid > 0){
230  HepMC::GenParticle* mother = particles[mid];
231  LogDebug("DecayChain")<<"Mother's PDG ID "<<mother->pdg_id();
232 
233  HepMC::GenVertex* prod_vertex = mother->end_vertex();
234  if(!prod_vertex){
235  prod_vertex = prods[ipart];
236  prod_vertex->add_particle_in(mother);
237  evt->add_vertex(prod_vertex);
238  prods[ipart]=nullptr; // mark to protect deletion
239 
240  }
241  prod_vertex->add_particle_out(part);
242  }
243  }
244 
245  // cleanup vertices not assigned to evt
246  for (unsigned int i = 0; i<prods.size(); i++) {
247  if(prods[i]) delete prods[i];
248  }
249 
250  return true;
251 }
252 
253 //_____________________________________________________________________
254 bool AMPTHadronizer::call_amptset(double efrm, std::string frame, std::string proj, std::string targ, int iap, int izp, int iat, int izt)
255 {
256  // initialize hydjet
257  AMPTSET(efrm,frame.data(),proj.data(),targ.data(),iap,izp,iat,izt,strlen(frame.data()),strlen(proj.data()),strlen(targ.data()));
258  return true;
259 }
260 //______________________________________________________________________
262 {
263  anim.isoft=amptmode_;
264  input2.ntmax=ntmax_;
265  input1.dt=dt_;
266  ludat1.parj[40]=stringFragA_;
267  ludat1.parj[41]=stringFragB_;
268  popcorn.ipop=popcornmode_;
269  ludat1.parj[4]=popcornpar_;
270  hparnt.ihpr2[5]=shadowingmode_;
271  hparnt.ihpr2[3]=quenchingmode_;
272  hparnt.hipr1[13]=quenchingpar_;
273  hparnt.hipr1[7]=pthard_;
274  para2.xmu=mu_;
275  anim.izpc=izpc_;
276  para2.alpha=alpha_;
277  coal.dpcoal=dpcoal_;
278  coal.drcoal=drcoal_;
279  resdcy.iksdcy=ks0decay_;
280  phidcy.iphidcy=phidecay_;
281  para8.idpert=deuteronmode_;
282  para8.npertd=deuteronfactor_;
283  para8.idxsec=deuteronxsec_;
284  phidcy.pttrig=minijetpt_;
285  phidcy.maxmiss=maxmiss_;
287  hparnt.ihpr2[4]=ktkick_;
288  embed.iembed=diquarkembedding_;
289  embed.pxqembd=diquarkpx_;
290  embed.pyqembd=diquarkpy_;
291  embed.xembd=diquarkx_;
292  embed.yembd=diquarky_;
293 
294  return true;
295 }
296 
297 //_____________________________________________________________________
299 
300  // ampt running options
301  ampt_init(pset_);
302 
303  // initialize ampt
304  LogInfo("AMPTinAction") << "##### Calling AMPTSET(" << efrm_ << "," <<frame_<<","<<proj_<<","<<targ_<<","<<iap_<<","<<izp_<<","<<iat_<<","<<izt_<<") ####";
305 
307 
308  return true;
309 }
310 
311 bool AMPTHadronizer::declareStableParticles( const std::vector<int>& pdg )
312 {
313  return true;
314 }
315 
316 //________________________________________________________________
318  int zero = 0;
319  double test = (double)gen::ranart_(&zero);
320  phi0_ = 2.*pi*test - pi;
321  sinphi0_ = sin(phi0_);
322  cosphi0_ = cos(phi0_);
323 }
324 
325 //________________________________________________________________
327 {
328  return false;
329 }
330 
332 {
333  return true;
334 }
335 
337 {
338  return true;
339 }
340 
342  return;
343 }
344 
346  return;
347 }
348 
349 const char* AMPTHadronizer::classname() const
350 {
351  return "gen::AMPTHadronizer";
352 }
353 
#define LogDebug(id)
#define popcorn
Definition: AMPTWrapper.h:185
float ran1_(int *)
#define ludat1
Definition: AMPTWrapper.h:70
bool initializeForInternalPartons()
#define phidcy
Definition: AMPTWrapper.h:167
#define AMPT
Definition: AMPTWrapper.h:19
HepMC::GenEvent * evt
void add_heavy_ion_rec(HepMC::GenEvent *evt)
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
bool generatePartonsAndHadronize()
#define hparnt
Definition: AMPTWrapper.h:60
static const double pi
#define anim
Definition: AMPTWrapper.h:80
#define nullptr
double v[5][pyjets_maxn]
#define input2
Definition: AMPTWrapper.h:149
std::auto_ptr< HepMC::GenEvent > & event()
bool declareStableParticles(const std::vector< int > &)
#define embed
Definition: AMPTWrapper.h:178
const char * classname() const
#define AMPTSET
Definition: AMPTWrapper.h:14
#define resdcy
Definition: AMPTWrapper.h:157
bool call_amptset(double efrm, std::string frame, std::string proj, std::string targ, int iap, int izp, int iat, int izt)
#define hbt
Definition: AMPTWrapper.h:50
double p[5][pyjets_maxn]
T sqrt(T t)
Definition: SSEVec.h:18
~AMPTHadronizer() override
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
#define hmain1
Definition: AMPTWrapper.h:40
#define para2
Definition: AMPTWrapper.h:100
#define input1
Definition: AMPTWrapper.h:129
edm::ParameterSet pset_
static CLHEP::HepRandomEngine * amptRandomEngine
float ranart_(int *)
part
Definition: HCALResponse.h:20
HLT enums.
#define para8
Definition: AMPTWrapper.h:118
void doSetRandomEngine(CLHEP::HepRandomEngine *v) override
HepMC::GenParticle * build_ampt(int index, int barcode)
bool ampt_init(const edm::ParameterSet &pset)
#define INVFLV
Definition: AMPTWrapper.h:24
HepMC::GenVertex * build_ampt_vertex(int i, int id)
#define coal
Definition: AMPTWrapper.h:89
bool get_particles(HepMC::GenEvent *evt)