File indexing completed on 2023-03-17 11:03:58
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
0098 HepMC::HeavyIon* hi = new HepMC::HeavyIon(hmain1.jatt,
0099 hmain1.np,
0100 hmain1.nt,
0101 hmain1.n0 + hmain1.n01 + hmain1.n10 + hmain1.n11,
0102 0,
0103 0,
0104 hmain1.n01,
0105 hmain1.n10,
0106 hmain1.n11,
0107 hparnt.hint1[18],
0108 phi0_,
0109 0,
0110 hparnt.hint1[11]
0111 );
0112 evt->set_heavy_ion(*hi);
0113 delete hi;
0114 }
0115
0116
0117 HepMC::GenParticle* AMPTHadronizer::build_ampt(int index, int barcode) {
0118
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]),
0133 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
0143 HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), id);
0144 return vertex;
0145 }
0146
0147 bool AMPTHadronizer::generatePartonsAndHadronize() {
0148
0149 if (rotate_)
0150 rotateEvtPlane();
0151
0152
0153 AMPT(frame_.c_str(), bmin_, bmax_, frame_.length());
0154
0155
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;
0213 }
0214 prod_vertex->add_particle_out(part);
0215 }
0216 }
0217
0218
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
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
0277 ampt_init(pset_);
0278
0279
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"; }