File indexing completed on 2024-04-06 12:13:43
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021 #include <string>
0022 #include <iostream>
0023
0024
0025 #include "FWCore/Framework/interface/Frameworkfwd.h"
0026 #include "FWCore/Framework/interface/global/EDProducer.h"
0027 #include "FWCore/Framework/interface/EventSetup.h"
0028
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/MakerMacros.h"
0031 #include "FWCore/Utilities/interface/ESGetToken.h"
0032
0033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0034 #include "FWCore/Utilities/interface/InputTag.h"
0035
0036 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0037 #include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
0038 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0039
0040 #include "HepMC/HeavyIon.h"
0041 #include "FWCore/Framework/interface/ESHandle.h"
0042 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0043 using namespace std;
0044
0045
0046
0047
0048
0049 class GenHIEventProducer : public edm::global::EDProducer<> {
0050 public:
0051 explicit GenHIEventProducer(const edm::ParameterSet&);
0052 ~GenHIEventProducer() override = default;
0053
0054 private:
0055 void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0056 const edm::EDGetTokenT<CrossingFrame<edm::HepMCProduct> > hepmcSrc_;
0057 const edm::ESGetToken<ParticleDataTable, edm::DefaultRecord> pdtToken_;
0058 const edm::EDPutTokenT<edm::GenHIEvent> putToken_;
0059
0060 double ptCut_;
0061 const bool doParticleInfo_;
0062 };
0063
0064
0065
0066
0067 GenHIEventProducer::GenHIEventProducer(const edm::ParameterSet& iConfig)
0068 : hepmcSrc_(consumes<CrossingFrame<edm::HepMCProduct> >(iConfig.getParameter<edm::InputTag>("src"))),
0069 pdtToken_(esConsumes<ParticleDataTable, edm::DefaultRecord>()),
0070 putToken_(produces<edm::GenHIEvent>()),
0071 doParticleInfo_(iConfig.getUntrackedParameter<bool>("doParticleInfo", false)) {
0072 if (doParticleInfo_) {
0073 ptCut_ = iConfig.getParameter<double>("ptCut");
0074 }
0075 }
0076
0077
0078
0079
0080
0081
0082 void GenHIEventProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0083 using namespace edm;
0084
0085 const auto& pdt = iSetup.getData(pdtToken_);
0086
0087 double b = -1;
0088 int npart = -1;
0089 int ncoll = 0;
0090 int nhard = 0;
0091 double phi = 0;
0092 double ecc = -1;
0093
0094 int nCharged = 0;
0095 int nChargedMR = 0;
0096 int nChargedPtCut = 0;
0097 int nChargedPtCutMR = 0;
0098
0099 double meanPt = 0;
0100 double meanPtMR = 0;
0101 double EtMR = 0;
0102 double TotEnergy = 0;
0103
0104 const auto& hepmc = iEvent.get(hepmcSrc_);
0105 MixCollection<HepMCProduct> mix(&hepmc);
0106
0107 if (mix.size() < 1) {
0108 throw cms::Exception("MatchVtx") << "Mixing has " << mix.size() << " sub-events, should have been at least 1"
0109 << endl;
0110 }
0111
0112 const HepMCProduct& hievt = mix.getObject(mix.size() - 1);
0113 const HepMC::GenEvent* evt = hievt.GetEvent();
0114 if (doParticleInfo_) {
0115 HepMC::GenEvent::particle_const_iterator begin = evt->particles_begin();
0116 HepMC::GenEvent::particle_const_iterator end = evt->particles_end();
0117 for (HepMC::GenEvent::particle_const_iterator it = begin; it != end; ++it) {
0118 if ((*it)->status() != 1)
0119 continue;
0120 int pdg_id = (*it)->pdg_id();
0121 const ParticleData* part = pdt.particle(pdg_id);
0122 int charge = static_cast<int>(part->charge());
0123
0124 if (charge == 0)
0125 continue;
0126 float pt = (*it)->momentum().perp();
0127 float eta = (*it)->momentum().eta();
0128 float energy = (*it)->momentum().e();
0129
0130 nCharged++;
0131 meanPt += pt;
0132
0133 if (fabs(eta) < 1.0) {
0134 TotEnergy += energy;
0135 }
0136 if (pt > ptCut_) {
0137 nChargedPtCut++;
0138 if (fabs(eta) < 0.5) {
0139 nChargedPtCutMR++;
0140 }
0141 }
0142
0143
0144 if (fabs(eta) > 0.5)
0145 continue;
0146 nChargedMR++;
0147 meanPtMR += pt;
0148 }
0149 }
0150 const HepMC::HeavyIon* hi = evt->heavy_ion();
0151
0152 if (hi) {
0153 ncoll = ncoll + hi->Ncoll();
0154 nhard = nhard + hi->Ncoll_hard();
0155 int np = hi->Npart_proj() + hi->Npart_targ();
0156 if (np >= 0) {
0157 npart = np;
0158 b = hi->impact_parameter();
0159 phi = hi->event_plane_angle();
0160 ecc = hi->eccentricity();
0161 }
0162 }
0163
0164
0165 if (TotEnergy != 0) {
0166 EtMR = TotEnergy / 2;
0167 }
0168
0169 if (nChargedMR != 0) {
0170 meanPtMR /= nChargedMR;
0171 }
0172 if (nCharged != 0) {
0173 meanPt /= nCharged;
0174 }
0175
0176 iEvent.emplace(putToken_,
0177 b,
0178 npart,
0179 ncoll,
0180 nhard,
0181 phi,
0182 ecc,
0183 nCharged,
0184 nChargedMR,
0185 meanPt,
0186 meanPtMR,
0187 EtMR,
0188 nChargedPtCut,
0189 nChargedPtCutMR);
0190 }
0191
0192
0193 DEFINE_FWK_MODULE(GenHIEventProducer);