AMPTInterface/src/AMPTHadronizer.cc

0001 #include <iostream>
0002 #include <cmath>
0003
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/Run.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008
0009 #include "GeneratorInterface/AMPTInterface/interface/AMPTHadronizer.h"
0010 #include "GeneratorInterface/AMPTInterface/interface/AMPTWrapper.h"
0011
0012 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0013 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0014 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0015
0016 #include "HepMC/GenEvent.h"
0017 #include "HepMC/HeavyIon.h"
0018 #include "HepMC/SimpleVector.h"
0019 #include "CLHEP/Random/RandomEngine.h"
0020
0021 static const double pi = 3.14159265358979;
0022
0023 using namespace edm;
0024 using namespace std;
0025 using namespace gen;
0026
0027 static CLHEP::HepRandomEngine* amptRandomEngine;
0028
0029 extern "C" {
0030 float gen::ranart_(int*) {
0031   float rannum = amptRandomEngine->flat();
0032   return rannum;
0033 }
0034 }
0035
0036 extern "C" {
0037 float gen::ran1_(int*) { return amptRandomEngine->flat(); }
0038 }
0039
0040 AMPTHadronizer::AMPTHadronizer(const ParameterSet& pset)
0041     : BaseHadronizer(pset),
0042       evt(nullptr),
0043       pset_(pset),
0044       bmax_(pset.getParameter<double>("bMax")),
0045       bmin_(pset.getParameter<double>("bMin")),
0046       efrm_(pset.getParameter<double>("comEnergy")),
0047       frame_(pset.getParameter<string>("frame")),
0048       proj_(pset.getParameter<string>("proj")),
0049       targ_(pset.getParameter<string>("targ")),
0050       iap_(pset.getParameter<int>("iap")),
0051       izp_(pset.getParameter<int>("izp")),
0052       iat_(pset.getParameter<int>("iat")),
0053       izt_(pset.getParameter<int>("izt")),
0054       amptmode_(pset.getParameter<int>("amptmode")),
0055       ntmax_(pset.getParameter<int>("ntmax")),
0056       dt_(pset.getParameter<double>("dt")),
0057       stringFragA_(pset.getParameter<double>("stringFragA")),
0058       stringFragB_(pset.getParameter<double>("stringFragB")),
0059       popcornmode_(pset.getParameter<bool>("popcornmode")),
0060       popcornpar_(pset.getParameter<double>("popcornpar")),
0061       shadowingmode_(pset.getParameter<bool>("shadowingmode")),
0062       quenchingmode_(pset.getParameter<bool>("quenchingmode")),
0063       quenchingpar_(pset.getParameter<double>("quenchingpar")),
0064       pthard_(pset.getParameter<double>("pthard")),
0065       mu_(pset.getParameter<double>("mu")),
0066       izpc_(pset.getParameter<int>("izpc")),
0067       alpha_(pset.getParameter<double>("alpha")),
0068       dpcoal_(pset.getParameter<double>("dpcoal")),
0069       drcoal_(pset.getParameter<double>("drcoal")),
0070       ks0decay_(pset.getParameter<bool>("ks0decay")),
0071       phidecay_(pset.getParameter<bool>("phidecay")),
0072       deuteronmode_(pset.getParameter<int>("deuteronmode")),
0073       deuteronfactor_(pset.getParameter<int>("deuteronfactor")),
0074       deuteronxsec_(pset.getParameter<int>("deuteronxsec")),
0075       minijetpt_(pset.getParameter<double>("minijetpt")),
0076       maxmiss_(pset.getParameter<int>("maxmiss")),
0077       doInitialAndFinalRadiation_(pset.getParameter<int>("doInitialAndFinalRadiation")),
0078       ktkick_(pset.getParameter<int>("ktkick")),
0079       diquarkembedding_(pset.getParameter<int>("diquarkembedding")),
0080       diquarkpx_(pset.getParameter<double>("diquarkpx")),
0081       diquarkpy_(pset.getParameter<double>("diquarkpy")),
0082       diquarkx_(pset.getParameter<double>("diquarkx")),
0083       diquarky_(pset.getParameter<double>("diquarky")),
0084       phi0_(0.),
0085       sinphi0_(0.),
0086       cosphi0_(1.),
0087       rotate_(pset.getParameter<bool>("rotateEventPlane")) {}
0088
0089 //_____________________________________________________________________
0090 AMPTHadronizer::~AMPTHadronizer() {}
0091
0092 //_____________________________________________________________________
0093 void AMPTHadronizer::doSetRandomEngine(CLHEP::HepRandomEngine* v) { amptRandomEngine = v; }
0094
0095 //_____________________________________________________________________
0096 void AMPTHadronizer::add_heavy_ion_rec(HepMC::GenEvent* evt) {
0097   // heavy ion record in the final CMSSW Event
0098   HepMC::HeavyIon* hi = new HepMC::HeavyIon(hmain1.jatt,  // Ncoll_hard/N of SubEvents
0099                                             hmain1.np,    // Npart_proj
0100                                             hmain1.nt,    // Npart_targ
0101                                             hmain1.n0 + hmain1.n01 + hmain1.n10 + hmain1.n11,  // Ncoll
0102                                             0,                                                 // spectator_neutrons
0103                                             0,                                                 // spectator_protons
0104                                             hmain1.n01,                                        // N_Nwounded_collisions
0105                                             hmain1.n10,                                        // Nwounded_N_collisions
0106                                             hmain1.n11,        // Nwounded_Nwounded_collisions
0107                                             hparnt.hint1[18],  // impact_parameter in [fm]
0108                                             phi0_,             // event_plane_angle
0109                                             0,                 // eccentricity
0110                                             hparnt.hint1[11]   // sigma_inel_NN
0111   );
0112   evt->set_heavy_ion(*hi);
0113   delete hi;
0114 }
0115
0116 //___________________________________________________________________
0117 HepMC::GenParticle* AMPTHadronizer::build_ampt(int index, int barcode) {
0118   // Build particle object corresponding to index in ampt
0119
0120   float px0 = hbt.plast[index][0];
0121   float py0 = hbt.plast[index][1];
0122   float pz0 = hbt.plast[index][2];
0123   float m = hbt.plast[index][3];
0124
0125   float px = px0 * cosphi0_ - py0 * sinphi0_;
0126   float py = py0 * cosphi0_ + px0 * sinphi0_;
0127   float pz = pz0;
0128   float e = sqrt(px * px + py * py + pz * pz + m * m);
0129   int status = 1;
0130
0131   HepMC::GenParticle* p = new HepMC::GenParticle(HepMC::FourVector(px, py, pz, e),
0132                                                  INVFLV(hbt.lblast[index]),  // id
0133                                                  status                      // status
0134   );
0135
0136   p->suggest_barcode(barcode);
0137   return p;
0138 }
0139
0140 //___________________________________________________________________
0141 HepMC::GenVertex* AMPTHadronizer::build_ampt_vertex(int i, int id) {
0142   // build verteces for the ampt stored events
0143   HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), id);
0144   return vertex;
0145 }
0146 //_____________________________________________________________________
0147 bool AMPTHadronizer::generatePartonsAndHadronize() {
0148   // generate single event
0149   if (rotate_)
0150     rotateEvtPlane();
0151
0152   // generate a AMPT event
0153   AMPT(frame_.c_str(), bmin_, bmax_, frame_.length());
0154
0155   // event information
0156   HepMC::GenEvent* evt = new HepMC::GenEvent();
0157   get_particles(evt);
0158
0159   add_heavy_ion_rec(evt);
0160
0161   event().reset(evt);
0162
0163   return true;
0164 }
0165
0166 //_____________________________________________________________________
0167 bool AMPTHadronizer::get_particles(HepMC::GenEvent* evt) {
0168   HepMC::GenVertex* vertice;
0169
0170   vector<HepMC::GenParticle*> particles;
0171   vector<int> mother_ids;
0172   vector<HepMC::GenVertex*> prods;
0173
0174   vertice = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), 0);
0175   evt->add_vertex(vertice);
0176   if (!evt->signal_process_vertex())
0177     evt->set_signal_process_vertex(vertice);
0178
0179   const unsigned int knumpart = hbt.nlast;
0180   for (unsigned int ipart = 0; ipart < knumpart; ipart++) {
0181     int mid = 0;
0182     particles.push_back(build_ampt(ipart, ipart + 1));
0183     prods.push_back(build_ampt_vertex(ipart, 0));
0184     mother_ids.push_back(mid);
0185     LogDebug("DecayChain") << "Mother index : " << mid;
0186   }
0187
0188   LogDebug("AMPT") << "Number of particles in vector " << particles.size();
0189
0190   for (unsigned int ipart = 0; ipart < particles.size(); ipart++) {
0191     HepMC::GenParticle* part = particles[ipart];
0192
0193     int mid = mother_ids[ipart];
0194     LogDebug("DecayChain") << "Particle " << ipart;
0195     LogDebug("DecayChain") << "Mother's ID " << mid;
0196     LogDebug("DecayChain") << "Particle's PDG ID " << part->pdg_id();
0197
0198     if (mid <= 0) {
0199       vertice->add_particle_out(part);
0200       continue;
0201     }
0202
0203     if (mid > 0) {
0204       HepMC::GenParticle* mother = particles[mid];
0205       LogDebug("DecayChain") << "Mother's PDG ID " << mother->pdg_id();
0206
0207       HepMC::GenVertex* prod_vertex = mother->end_vertex();
0208       if (!prod_vertex) {
0209         prod_vertex = prods[ipart];
0210         prod_vertex->add_particle_in(mother);
0211         evt->add_vertex(prod_vertex);
0212         prods[ipart] = nullptr;  // mark to protect deletion
0213       }
0214       prod_vertex->add_particle_out(part);
0215     }
0216   }
0217
0218   // cleanup vertices not assigned to evt
0219   for (unsigned int i = 0; i < prods.size(); i++) {
0220     if (prods[i])
0221       delete prods[i];
0222   }
0223
0224   return true;
0225 }
0226
0227 //_____________________________________________________________________
0228 bool AMPTHadronizer::call_amptset(
0229     double efrm, std::string frame, std::string proj, std::string targ, int iap, int izp, int iat, int izt) {
0230   // initialize hydjet
0231   frame.resize(4, ' ');
0232   proj.resize(4, ' ');
0233   targ.resize(4, ' ');
0234   AMPTSET(
0235       efrm, frame.data(), proj.data(), targ.data(), iap, izp, iat, izt, frame.length(), proj.length(), targ.length());
0236   return true;
0237 }
0238 //______________________________________________________________________
0239 bool AMPTHadronizer::ampt_init(const ParameterSet& pset) {
0240   anim.isoft = amptmode_;
0241   input2.ntmax = ntmax_;
0242   input1.dt = dt_;
0243   ludat1.parj[40] = stringFragA_;
0244   ludat1.parj[41] = stringFragB_;
0245   popcorn.ipop = popcornmode_;
0246   ludat1.parj[4] = popcornpar_;
0247   hparnt.ihpr2[5] = shadowingmode_;
0248   hparnt.ihpr2[3] = quenchingmode_;
0249   hparnt.hipr1[13] = quenchingpar_;
0250   hparnt.hipr1[7] = pthard_;
0251   para2.xmu = mu_;
0252   anim.izpc = izpc_;
0253   para2.alpha = alpha_;
0254   coal.dpcoal = dpcoal_;
0255   coal.drcoal = drcoal_;
0256   resdcy.iksdcy = ks0decay_;
0257   phidcy.iphidcy = phidecay_;
0258   para8.idpert = deuteronmode_;
0259   para8.npertd = deuteronfactor_;
0260   para8.idxsec = deuteronxsec_;
0261   phidcy.pttrig = minijetpt_;
0262   phidcy.maxmiss = maxmiss_;
0263   hparnt.ihpr2[1] = doInitialAndFinalRadiation_;
0264   hparnt.ihpr2[4] = ktkick_;
0265   embed.iembed = diquarkembedding_;
0266   embed.pxqembd = diquarkpx_;
0267   embed.pyqembd = diquarkpy_;
0268   embed.xembd = diquarkx_;
0269   embed.yembd = diquarky_;
0270
0271   return true;
0272 }
0273
0274 //_____________________________________________________________________
0275 bool AMPTHadronizer::initializeForInternalPartons() {
0276   // ampt running options
0277   ampt_init(pset_);
0278
0279   // initialize ampt
0280   LogInfo("AMPTinAction") << "##### Calling AMPTSET(" << efrm_ << "," << frame_ << "," << proj_ << "," << targ_ << ","
0281                           << iap_ << "," << izp_ << "," << iat_ << "," << izt_ << ") ####";
0282
0283   call_amptset(efrm_, frame_, proj_, targ_, iap_, izp_, iat_, izt_);
0284
0285   return true;
0286 }
0287
0288 bool AMPTHadronizer::declareStableParticles(const std::vector<int>& pdg) { return true; }
0289
0290 //________________________________________________________________
0291 void AMPTHadronizer::rotateEvtPlane() {
0292   int zero = 0;
0293   double test = (double)gen::ranart_(&zero);
0294   phi0_ = 2. * pi * test - pi;
0295   sinphi0_ = sin(phi0_);
0296   cosphi0_ = cos(phi0_);
0297 }
0298
0299 //________________________________________________________________
0300 bool AMPTHadronizer::hadronize() { return false; }
0301
0302 bool AMPTHadronizer::decay() { return true; }
0303
0304 bool AMPTHadronizer::residualDecay() { return true; }
0305
0306 void AMPTHadronizer::finalizeEvent() { return; }
0307
0308 void AMPTHadronizer::statistics() { return; }
0309
0310 const char* AMPTHadronizer::classname() const { return "gen::AMPTHadronizer"; }