Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:13:48 

0001 /**
0002  *    \brief Interface to the HYDJET++ (Hydjet2) generator (since core v. 2.4.3), produces HepMC events
0003  *    \version 1.3
0004  *    \author Andrey Belyaev
0005  */
0006 
0007 #include <TLorentzVector.h>
0008 #include <TMath.h>
0009 #include <TVector3.h>
0010 
0011 #include "GeneratorInterface/Hydjet2Interface/interface/Hydjet2Hadronizer.h"
0012 #include <cmath>
0013 #include <fstream>
0014 #include <iostream>
0015 
0016 #include "FWCore/Concurrency/interface/SharedResourceNames.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/Run.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Utilities/interface/EDMException.h"
0022 
0023 #include "GeneratorInterface/Pythia6Interface/interface/Pythia6Declarations.h"
0024 #include "GeneratorInterface/Pythia6Interface/interface/Pythia6Service.h"
0025 
0026 #include "HepMC/GenEvent.h"
0027 #include "HepMC/HeavyIon.h"
0028 #include "HepMC/IO_HEPEVT.h"
0029 #include "HepMC/PythiaWrapper6_4.h"
0030 #include "HepMC/SimpleVector.h"
0031 
0032 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0033 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0034 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0035 #include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
0036 
0037 CLHEP::HepRandomEngine *hjRandomEngine;
0038 
0039 using namespace edm;
0040 using namespace std;
0041 using namespace gen;
0042 
0043 int Hydjet2Hadronizer::convertStatusForComponents(int sta, int typ, int pySta) {
0044   int st = -1;
0045   if (typ == 0)  //soft
0046     st = 2 - sta;
0047   else if (typ == 1)
0048     st = convertStatus(pySta);
0049 
0050   if (st == -1)
0051     throw cms::Exception("ConvertStatus") << "Wrong status code!" << endl;
0052 
0053   if (separateHydjetComponents_) {
0054     if (st == 1 && typ == 0)
0055       return 6;
0056     if (st == 1 && typ == 1)
0057       return 7;
0058     if (st == 2 && typ == 0)
0059       return 16;
0060     if (st == 2 && typ == 1)
0061       return 17;
0062   }
0063   return st;
0064 }
0065 
0066 int Hydjet2Hadronizer::convertStatus(int st) {
0067   if (st <= 0)
0068     return 0;
0069   if (st <= 10)
0070     return 1;
0071   if (st <= 20)
0072     return 2;
0073   if (st <= 30)
0074     return 3;
0075   else
0076     return -1;
0077 }
0078 
0079 const std::vector<std::string> Hydjet2Hadronizer::theSharedResources = {edm::SharedResourceNames::kPythia6};
0080 
0081 //____________________________________________________________________________________________
0082 Hydjet2Hadronizer::Hydjet2Hadronizer(const edm::ParameterSet &pset, edm::ConsumesCollector &&iC)
0083     : BaseHadronizer(pset),
0084       rotate_(pset.getParameter<bool>("rotateEventPlane")),
0085       evt(nullptr),
0086       nsub_(0),
0087       nhard_(0),
0088       nsoft_(0),
0089       phi0_(0.),
0090       sinphi0_(0.),
0091       cosphi0_(1.),
0092       fVertex_(nullptr),
0093       pythia6Service_(new Pythia6Service(pset))
0094 
0095 {
0096   fParams.doPrintInfo = false;
0097   fParams.allowEmptyEvent = false;
0098   fParams.fNevnt = 0;                                      //not used in CMSSW
0099   fParams.femb = pset.getParameter<int>("embeddingMode");  //
0100   fParams.fSqrtS = pset.getParameter<double>("fSqrtS");    // C.m.s. energy per nucleon pair
0101   fParams.fAw = pset.getParameter<double>("fAw");          // Atomic weigth of nuclei, fAw
0102   fParams.fIfb = pset.getParameter<int>(
0103       "fIfb");  // Flag of type of centrality generation, fBfix (=0 is fixed by fBfix, >0 distributed [fBfmin, fBmax])
0104   fParams.fBmin = pset.getParameter<double>("fBmin");  // Minimum impact parameter in units of nuclear radius, fBmin
0105   fParams.fBmax = pset.getParameter<double>("fBmax");  // Maximum impact parameter in units of nuclear radius, fBmax
0106   fParams.fBfix = pset.getParameter<double>("fBfix");  // Fixed impact parameter in units of nuclear radius, fBfix
0107   fParams.fT = pset.getParameter<double>("fT");        // Temperature at chemical freeze-out, fT [GeV]
0108   fParams.fMuB = pset.getParameter<double>("fMuB");    // Chemical baryon potential per unit charge, fMuB [GeV]
0109   fParams.fMuS = pset.getParameter<double>("fMuS");    // Chemical strangeness potential per unit charge, fMuS [GeV]
0110   fParams.fMuC = pset.getParameter<double>(
0111       "fMuC");  // Chemical charm potential per unit charge, fMuC [GeV] (used if charm production is turned on)
0112   fParams.fMuI3 = pset.getParameter<double>("fMuI3");  // Chemical isospin potential per unit charge, fMuI3 [GeV]
0113   fParams.fThFO = pset.getParameter<double>("fThFO");  // Temperature at thermal freeze-out, fThFO [GeV]
0114   fParams.fMu_th_pip =
0115       pset.getParameter<double>("fMu_th_pip");  // Chemical potential of pi+ at thermal freeze-out, fMu_th_pip [GeV]
0116   fParams.fTau = pset.getParameter<double>(
0117       "fTau");  // Proper time proper at thermal freeze-out for central collisions, fTau [fm/c]
0118   fParams.fSigmaTau = pset.getParameter<double>(
0119       "fSigmaTau");  // Duration of emission at thermal freeze-out for central collisions, fSigmaTau [fm/c]
0120   fParams.fR = pset.getParameter<double>(
0121       "fR");  // Maximal transverse radius at thermal freeze-out for central collisions, fR [fm]
0122   fParams.fYlmax =
0123       pset.getParameter<double>("fYlmax");  // Maximal longitudinal flow rapidity at thermal freeze-out, fYlmax
0124   fParams.fUmax = pset.getParameter<double>(
0125       "fUmax");  // Maximal transverse flow rapidity at thermal freeze-out for central collisions, fUmax
0126   fParams.frhou2 = pset.getParameter<double>("fRhou2");  //parameter to swich ON/OFF = 0) rhou2
0127   fParams.frhou3 = pset.getParameter<double>("fRhou3");  //parameter to swich ON/OFF(0) rhou3
0128   fParams.frhou4 = pset.getParameter<double>("fRhou4");  //parameter to swich ON/OFF(0) rhou4
0129   fParams.fDelta =
0130       pset.getParameter<double>("fDelta");  // Momentum azimuthal anizotropy parameter at thermal freeze-out, fDelta
0131   fParams.fEpsilon =
0132       pset.getParameter<double>("fEpsilon");  // Spatial azimuthal anisotropy parameter at thermal freeze-out, fEpsilon
0133   fParams.fv2 = pset.getParameter<double>("fKeps2");  //parameter to swich ON/OFF(0) epsilon2 fluctuations
0134   fParams.fv3 = pset.getParameter<double>("fKeps3");  //parameter to swich ON/OFF(0) epsilon3 fluctuations
0135   fParams.fIfDeltaEpsilon = pset.getParameter<double>(
0136       "fIfDeltaEpsilon");  // Flag to specify fDelta and fEpsilon values, fIfDeltaEpsilon (=0 user's ones, >=1 calculated)
0137   fParams.fDecay =
0138       pset.getParameter<int>("fDecay");  // Flag to switch on/off hadron decays, fDecay (=0 decays off, >=1 decays on)
0139   fParams.fWeakDecay = pset.getParameter<double>(
0140       "fWeakDecay");  // Low decay width threshold fWeakDecay[GeV]: width<fWeakDecay decay off, width>=fDecayWidth decay on; can be used to switch off weak decays
0141   fParams.fEtaType = pset.getParameter<double>(
0142       "fEtaType");  // Flag to choose longitudinal flow rapidity distribution, fEtaType (=0 uniform, >0 Gaussian with the dispersion Ylmax)
0143   fParams.fTMuType = pset.getParameter<double>(
0144       "fTMuType");  // Flag to use calculated T_ch, mu_B and mu_S as a function of fSqrtS, fTMuType (=0 user's ones, >0 calculated)
0145   fParams.fCorrS = pset.getParameter<double>(
0146       "fCorrS");  // Strangeness supression factor gamma_s with fCorrS value (0<fCorrS <=1, if fCorrS <= 0 then it is calculated)
0147   fParams.fCharmProd = pset.getParameter<int>(
0148       "fCharmProd");  // Flag to include thermal charm production, fCharmProd (=0 no charm production, >=1 charm production)
0149   fParams.fCorrC = pset.getParameter<double>(
0150       "fCorrC");  // Charmness enhancement factor gamma_c with fCorrC value (fCorrC >0, if fCorrC<0 then it is calculated)
0151   fParams.fNhsel = pset.getParameter<int>(
0152       "fNhsel");  //Flag to include jet (J)/jet quenching (JQ) and hydro (H) state production, fNhsel (0 H on & J off, 1 H/J on & JQ off, 2 H/J/HQ on, 3 J on & H/JQ off, 4 H off & J/JQ on)
0153   fParams.fPyhist = pset.getParameter<int>(
0154       "fPyhist");  // Flag to suppress the output of particle history from PYTHIA, fPyhist (=1 only final state particles; =0 full particle history from PYTHIA)
0155   fParams.fIshad = pset.getParameter<int>(
0156       "fIshad");  // Flag to switch on/off nuclear shadowing, fIshad (0 shadowing off, 1 shadowing on)
0157   fParams.fPtmin =
0158       pset.getParameter<double>("fPtmin");  // Minimal pt of parton-parton scattering in PYTHIA event, fPtmin [GeV/c]
0159   fParams.fT0 = pset.getParameter<double>(
0160       "fT0");  // Initial QGP temperature for central Pb+Pb collisions in mid-rapidity, fT0 [GeV]
0161   fParams.fTau0 = pset.getParameter<double>("fTau0");  // Proper QGP formation time in fm/c, fTau0 (0.01<fTau0<10)
0162   fParams.fNf = pset.getParameter<int>("fNf");         // Number of active quark flavours in QGP, fNf (0, 1, 2 or 3)
0163   fParams.fIenglu = pset.getParameter<int>(
0164       "fIenglu");  // Flag to fix type of partonic energy loss, fIenglu (0 radiative and collisional loss, 1 radiative loss only, 2 collisional loss only)
0165   fParams.fIanglu = pset.getParameter<int>(
0166       "fIanglu");  // Flag to fix type of angular distribution of in-medium emitted gluons, fIanglu (0 small-angular, 1 wide-angular, 2 collinear).
0167 
0168   edm::FileInPath f1("externals/hydjet2/particles.data");
0169   strcpy(fParams.partDat, (f1.fullPath()).c_str());
0170 
0171   edm::FileInPath f2("externals/hydjet2/tabledecay.txt");
0172   strcpy(fParams.tabDecay, (f2.fullPath()).c_str());
0173 
0174   fParams.fPythiaTune = false;
0175 
0176   if (pset.exists("signalVtx"))
0177     signalVtx_ = pset.getUntrackedParameter<std::vector<double>>("signalVtx");
0178 
0179   if (signalVtx_.size() == 4) {
0180     if (!fVertex_)
0181       fVertex_ = new HepMC::FourVector();
0182     LogDebug("EventSignalVertex") << "Setting event signal vertex "
0183                                   << " x = " << signalVtx_.at(0) << " y = " << signalVtx_.at(1)
0184                                   << "  z= " << signalVtx_.at(2) << " t = " << signalVtx_.at(3) << endl;
0185     fVertex_->set(signalVtx_.at(0), signalVtx_.at(1), signalVtx_.at(2), signalVtx_.at(3));
0186   }
0187 
0188   // PYLIST Verbosity Level
0189   // Valid PYLIST arguments are: 1, 2, 3, 5, 7, 11, 12, 13
0190   pythiaPylistVerbosity_ = pset.getUntrackedParameter<int>("pythiaPylistVerbosity", 0);
0191   LogDebug("PYLISTverbosity") << "Pythia PYLIST verbosity level = " << pythiaPylistVerbosity_;
0192   //Max number of events printed on verbosity level
0193   maxEventsToPrint_ = pset.getUntrackedParameter<int>("maxEventsToPrint", 0);
0194   LogDebug("Events2Print") << "Number of events to be printed = " << maxEventsToPrint_;
0195   if (fParams.femb == 1) {
0196     fParams.fIfb = 0;
0197     src_ = iC.consumes<CrossingFrame<edm::HepMCProduct>>(
0198         pset.getUntrackedParameter<edm::InputTag>("backgroundLabel", edm::InputTag("mix", "generatorSmeared")));
0199   }
0200 
0201   separateHydjetComponents_ = pset.getUntrackedParameter<bool>("separateHydjetComponents", false);
0202 }
0203 //__________________________________________________________________________________________
0204 Hydjet2Hadronizer::~Hydjet2Hadronizer() {
0205   call_pystat(1);
0206   delete pythia6Service_;
0207 }
0208 
0209 //_____________________________________________________________________
0210 void Hydjet2Hadronizer::doSetRandomEngine(CLHEP::HepRandomEngine *v) {
0211   pythia6Service_->setRandomEngine(v);
0212   hjRandomEngine = v;
0213 }
0214 
0215 //______________________________________________________________________________________________________
0216 bool Hydjet2Hadronizer::readSettings(int) {
0217   Pythia6Service::InstanceWrapper guard(pythia6Service_);
0218   pythia6Service_->setGeneralParams();
0219 
0220   fParams.fSeed = hjRandomEngine->CLHEP::HepRandomEngine::getSeed();
0221   LogInfo("Hydjet2Hadronizer|GenSeed") << "Seed for random number generation: "
0222                                        << hjRandomEngine->CLHEP::HepRandomEngine::getSeed();
0223 
0224   return kTRUE;
0225 }
0226 
0227 //______________________________________________________________________________________________________
0228 bool Hydjet2Hadronizer::initializeForInternalPartons() {
0229   Pythia6Service::InstanceWrapper guard(pythia6Service_);
0230 
0231   // the input impact parameter (bxx_) is in [fm]; transform in [fm/RA] for hydjet usage
0232   const double ra = nuclear_radius();
0233   LogInfo("Hydjet2Hadronizer|RAScaling") << "Nuclear radius(RA) =  " << ra;
0234   fParams.fBmin /= ra;
0235   fParams.fBmax /= ra;
0236   fParams.fBfix /= ra;
0237 
0238   hj2 = new Hydjet2(fParams);
0239 
0240   return kTRUE;
0241 }
0242 
0243 //__________________________________________________________________________________________
0244 bool Hydjet2Hadronizer::generatePartonsAndHadronize() {
0245   Pythia6Service::InstanceWrapper guard(pythia6Service_);
0246 
0247   // generate single event
0248   if (fParams.femb == 1) {
0249     const edm::Event &e = getEDMEvent();
0250     HepMC::GenVertex *genvtx = nullptr;
0251     const HepMC::GenEvent *inev = nullptr;
0252     Handle<CrossingFrame<HepMCProduct>> cf;
0253     e.getByToken(src_, cf);
0254     MixCollection<HepMCProduct> mix(cf.product());
0255     if (mix.size() < 1) {
0256       throw cms::Exception("MatchVtx") << "Mixing has " << mix.size() << " sub-events, should have been at least 1"
0257                                        << endl;
0258     }
0259     const HepMCProduct &bkg = mix.getObject(0);
0260     if (!(bkg.isVtxGenApplied())) {
0261       throw cms::Exception("MatchVtx") << "Input background does not have smeared vertex!" << endl;
0262     } else {
0263       inev = bkg.GetEvent();
0264     }
0265 
0266     genvtx = inev->signal_process_vertex();
0267 
0268     if (!genvtx)
0269       throw cms::Exception("MatchVtx") << "Input background does not have signal process vertex!" << endl;
0270 
0271     double aX, aY, aZ, aT;
0272 
0273     aX = genvtx->position().x();
0274     aY = genvtx->position().y();
0275     aZ = genvtx->position().z();
0276     aT = genvtx->position().t();
0277 
0278     if (!fVertex_) {
0279       fVertex_ = new HepMC::FourVector();
0280     }
0281     LogInfo("MatchVtx") << " setting vertex "
0282                         << " aX " << aX << " aY " << aY << " aZ " << aZ << " aT " << aT << endl;
0283     fVertex_->set(aX, aY, aZ, aT);
0284 
0285     const HepMC::HeavyIon *hi = inev->heavy_ion();
0286 
0287     if (hi) {
0288       fParams.fBfix = (hi->impact_parameter()) / nuclear_radius();
0289       phi0_ = hi->event_plane_angle();
0290       sinphi0_ = sin(phi0_);
0291       cosphi0_ = cos(phi0_);
0292     } else {
0293       LogWarning("EventEmbedding") << "Background event does not have heavy ion record!";
0294     }
0295 
0296   } else if (rotate_)
0297     rotateEvtPlane();
0298 
0299   nsoft_ = 0;
0300   nhard_ = 0;
0301 
0302   // generate one HYDJET event
0303   int ntry = 0;
0304 
0305   while (nsoft_ == 0 && nhard_ == 0) {
0306     if (ntry > 100) {
0307       LogError("Hydjet2EmptyEvent") << "##### HYDJET2: No Particles generated, Number of tries =" << ntry;
0308       // Throw an exception. Use the EventCorruption exception since it maps onto SkipEvent
0309       // which is what we want to do here.
0310       std::ostringstream sstr;
0311       sstr << "Hydjet2HadronizerProducer: No particles generated after " << ntry << " tries.\n";
0312       edm::Exception except(edm::errors::EventCorruption, sstr.str());
0313       throw except;
0314     } else {
0315       hj2->GenerateEvent(fParams.fBfix);
0316 
0317       if (hj2->IsEmpty()) {
0318         continue;
0319       }
0320 
0321       nsoft_ = hj2->GetNhyd();
0322       nsub_ = hj2->GetNjet();
0323       nhard_ = hj2->GetNpyt();
0324 
0325       //100 trys
0326       ++ntry;
0327     }
0328   }
0329 
0330   if (ev == 0) {
0331     Sigin = hj2->GetSigin();
0332     Sigjet = hj2->GetSigjet();
0333   }
0334   ev = true;
0335 
0336   if (fParams.fNhsel < 3)
0337     nsub_++;
0338 
0339   // event information
0340   std::unique_ptr<HepMC::GenEvent> evt = std::make_unique<HepMC::GenEvent>();
0341   std::unique_ptr<edm::HepMCProduct> HepMCEvt = std::make_unique<edm::HepMCProduct>();
0342 
0343   if (nhard_ > 0 || nsoft_ > 0)
0344     get_particles(evt.get());
0345 
0346   evt->set_signal_process_id(pypars.msti[0]);  // type of the process
0347   evt->set_event_scale(pypars.pari[16]);       // Q^2
0348   add_heavy_ion_rec(evt.get());
0349 
0350   if (fVertex_) {
0351     // generate new vertex & apply the shift
0352     // Copy the HepMC::GenEvent
0353     HepMCEvt = std::make_unique<edm::HepMCProduct>(evt.get());
0354     HepMCEvt->applyVtxGen(fVertex_);
0355     evt = std::make_unique<HepMC::GenEvent>(*HepMCEvt->GetEvent());
0356   }
0357 
0358   HepMC::HEPEVT_Wrapper::check_hepevt_consistency();
0359   LogDebug("HEPEVT_info") << "Ev numb: " << HepMC::HEPEVT_Wrapper::event_number()
0360                           << " Entries number: " << HepMC::HEPEVT_Wrapper::number_entries() << " Max. entries "
0361                           << HepMC::HEPEVT_Wrapper::max_number_entries() << std::endl;
0362 
0363   event() = std::move(evt);
0364   return kTRUE;
0365 }
0366 
0367 //________________________________________________________________
0368 bool Hydjet2Hadronizer::declareStableParticles(const std::vector<int> &_pdg) {
0369   std::vector<int> pdg = _pdg;
0370   for (size_t i = 0; i < pdg.size(); i++) {
0371     int pyCode = pycomp_(pdg[i]);
0372     std::ostringstream pyCard;
0373     pyCard << "MDCY(" << pyCode << ",1)=0";
0374     std::cout << pyCard.str() << std::endl;
0375     call_pygive(pyCard.str());
0376   }
0377   return true;
0378 }
0379 //________________________________________________________________
0380 bool Hydjet2Hadronizer::hadronize() { return false; }
0381 bool Hydjet2Hadronizer::decay() { return true; }
0382 bool Hydjet2Hadronizer::residualDecay() { return true; }
0383 void Hydjet2Hadronizer::finalizeEvent() {}
0384 void Hydjet2Hadronizer::statistics() {}
0385 const char *Hydjet2Hadronizer::classname() const { return "gen::Hydjet2Hadronizer"; }
0386 
0387 //________________________________________________________________
0388 void Hydjet2Hadronizer::rotateEvtPlane() {
0389   const double pi = 3.14159265358979;
0390   phi0_ = 2. * pi * gen::pyr_(nullptr) - pi;
0391   sinphi0_ = sin(phi0_);
0392   cosphi0_ = cos(phi0_);
0393 }
0394 
0395 //_____________________________________________________________________
0396 bool Hydjet2Hadronizer::get_particles(HepMC::GenEvent *evt) {
0397   LogDebug("Hydjet2") << " Number of sub events " << nsub_;
0398   LogDebug("Hydjet2") << " Number of hard events " << hj2->GetNjet();
0399   LogDebug("Hydjet2") << " Number of hard particles " << nhard_;
0400   LogDebug("Hydjet2") << " Number of soft particles " << nsoft_;
0401   LogDebug("Hydjet2") << " nhard_ + nsoft_ = " << nhard_ + nsoft_ << " Ntot = " << hj2->GetNtot() << endl;
0402 
0403   int ihy = 0;
0404   int isub_l = -1;
0405   int stab = 0;
0406 
0407   vector<HepMC::GenParticle *> particle(hj2->GetNtot());
0408   HepMC::GenVertex *sub_vertices = nullptr;
0409 
0410   while (ihy < hj2->GetNtot()) {
0411     if ((hj2->GetiJet().at(ihy)) != isub_l) {
0412       sub_vertices = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), hj2->GetiJet().at(ihy));
0413       evt->add_vertex(sub_vertices);
0414       if (!evt->signal_process_vertex())
0415         evt->set_signal_process_vertex(sub_vertices);
0416       isub_l = hj2->GetiJet().at(ihy);
0417     }
0418 
0419     if ((convertStatusForComponents(
0420             (hj2->GetFinal()).at(ihy), (hj2->GetType()).at(ihy), (hj2->GetPythiaStatus().at(ihy)))) == 1)
0421       stab++;
0422     LogDebug("Hydjet2_array") << ihy << " MULTin ev.:" << hj2->GetNtot() << " SubEv.#" << hj2->GetiJet().at(ihy)
0423                               << " Part #" << ihy + 1 << ", PDG: " << hj2->GetPdg().at(ihy) << " (st. "
0424                               << convertStatus(hj2->GetPythiaStatus().at(ihy))
0425                               << ") mother=" << hj2->GetMotherIndex().at(ihy) + 1 << ", childs ("
0426                               << hj2->GetFirstDaughterIndex().at(ihy) + 1 << "-"
0427                               << hj2->GetLastDaughterIndex().at(ihy) + 1 << "), vtx (" << hj2->GetX().at(ihy) << ","
0428                               << hj2->GetY().at(ihy) << "," << hj2->GetZ().at(ihy) << ") " << std::endl;
0429 
0430     if ((hj2->GetMotherIndex().at(ihy)) <= 0) {
0431       particle.at(ihy) = build_hyjet2(ihy, ihy + 1);
0432       if (!sub_vertices)
0433         LogError("Hydjet2_array") << "##### HYDJET2: Vertex not initialized!";
0434       else
0435         sub_vertices->add_particle_out(particle.at(ihy));
0436       LogDebug("Hydjet2_array") << " ---> " << ihy + 1 << std::endl;
0437     } else {
0438       particle.at(ihy) = build_hyjet2(ihy, ihy + 1);
0439       int mid = hj2->GetMotherIndex().at(ihy);
0440 
0441       while (((mid + 1) < ihy) && (std::abs(hj2->GetPdg().at(mid)) < 100) &&
0442              ((hj2->GetFirstDaughterIndex().at(mid + 1)) <= ihy)) {
0443         mid++;
0444         LogDebug("Hydjet2_array") << "======== MID changed to " << mid
0445                                   << " ======== PDG(mid) = " << hj2->GetPdg().at(mid) << std::endl;
0446       }
0447 
0448       if (std::abs(hj2->GetPdg().at(mid)) < 100) {
0449         mid = hj2->GetMotherIndex().at(ihy);
0450         LogDebug("Hydjet2_array") << "======== MID changed BACK to " << mid
0451                                   << " ======== PDG(mid) = " << hj2->GetPdg().at(mid) << std::endl;
0452       }
0453 
0454       HepMC::GenParticle *mother = particle.at(mid);
0455       HepMC::GenVertex *prod_vertex = mother->end_vertex();
0456 
0457       if (!prod_vertex) {
0458         prod_vertex = build_hyjet2_vertex(ihy, (hj2->GetiJet().at(ihy)));
0459         prod_vertex->add_particle_in(mother);
0460         LogDebug("Hydjet2_array") << " <--- " << mid + 1 << std::endl;
0461         evt->add_vertex(prod_vertex);
0462       }
0463       prod_vertex->add_particle_out(particle.at(ihy));
0464       LogDebug("Hydjet2_array") << " ---" << mid + 1 << "---> " << ihy + 1 << std::endl;
0465     }
0466     ihy++;
0467   }
0468 
0469   LogDebug("Hydjet2_array") << " MULTin ev.:" << hj2->GetNtot() << ", last index: " << ihy - 1
0470                             << ", stable particles: " << stab << std::endl;
0471   return kTRUE;
0472 }
0473 
0474 //___________________________________________________________________
0475 HepMC::GenParticle *Hydjet2Hadronizer::build_hyjet2(int index, int barcode) {
0476   // Build particle object corresponding to index in hyjets (soft+hard)
0477   double px0 = (hj2->GetPx()).at(index);
0478   double py0 = (hj2->GetPy()).at(index);
0479 
0480   double px = px0 * cosphi0_ - py0 * sinphi0_;
0481   double py = py0 * cosphi0_ + px0 * sinphi0_;
0482 
0483   HepMC::GenParticle *p = new HepMC::GenParticle(
0484       HepMC::FourVector(px,                        // px
0485                         py,                        // py
0486                         (hj2->GetPz()).at(index),  // pz
0487                         (hj2->GetE()).at(index)),  // E
0488       (hj2->GetPdg()).at(index),                   // id
0489       convertStatusForComponents(
0490           (hj2->GetFinal()).at(index), (hj2->GetType()).at(index), (hj2->GetPythiaStatus()).at(index))  // status
0491   );
0492 
0493   p->suggest_barcode(barcode);
0494   return p;
0495 }
0496 
0497 //___________________________________________________________________
0498 HepMC::GenVertex *Hydjet2Hadronizer::build_hyjet2_vertex(int i, int id) {
0499   // build verteces for the hyjets stored events
0500   double x0 = (hj2->GetX()).at(i);
0501   double y0 = (hj2->GetY()).at(i);
0502 
0503   // convert to mm (as in PYTHIA6)
0504   const double fm_to_mm = 1e-12;
0505   double x = fm_to_mm * (x0 * cosphi0_ - y0 * sinphi0_);
0506   double y = fm_to_mm * (y0 * cosphi0_ + x0 * sinphi0_);
0507   double z = fm_to_mm * (hj2->GetZ()).at(i);
0508   double t = fm_to_mm * (hj2->GetT()).at(i);
0509 
0510   HepMC::GenVertex *vertex = new HepMC::GenVertex(HepMC::FourVector(x, y, z, t), id);
0511 
0512   return vertex;
0513 }
0514 
0515 //_____________________________________________________________________
0516 void Hydjet2Hadronizer::add_heavy_ion_rec(HepMC::GenEvent *evt) {
0517   // heavy ion record in the final CMSSW Event
0518   int nproj = static_cast<int>((hj2->GetNpart()) / 2);
0519   int ntarg = static_cast<int>((hj2->GetNpart()) - nproj);
0520 
0521   HepMC::HeavyIon *hi = new HepMC::HeavyIon(nsub_,                              // Ncoll_hard/N of SubEvents
0522                                             nproj,                              // Npart_proj
0523                                             ntarg,                              // Npart_targ
0524                                             hj2->GetNbcol(),                    // Ncoll
0525                                             0,                                  // spectator_neutrons
0526                                             0,                                  // spectator_protons
0527                                             0,                                  // N_Nwounded_collisions
0528                                             0,                                  // Nwounded_N_collisions
0529                                             0,                                  // Nwounded_Nwounded_collisions
0530                                             hj2->GetBgen() * nuclear_radius(),  // impact_parameter in [fm]
0531                                             phi0_,                              // event_plane_angle
0532                                             hj2->GetPsiv3(),                    // eccentricity <<<---- psi for v3!!!
0533                                             Sigin                               // sigma_inel_NN
0534   );
0535 
0536   evt->set_heavy_ion(*hi);
0537   delete hi;
0538 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!