HiGenCommon/plugins/GenHIEventProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    GenHIEventProducer
0004 // Class:      GenHIEventProducer
0005 //
0006 /**\class GenHIEventProducer GenHIEventProducer.cc yetkin/GenHIEventProducer/src/GenHIEventProducer.cc
0007
0008 Description: <one line class summary>
0009
0010 Implementation:
0011 <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Yetkin Yilmaz
0015 //         Created:  Thu Aug 13 08:39:51 EDT 2009
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021 #include <string>
0022 #include <iostream>
0023
0024 // user include files
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 // class decleration
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 // constructors and destructor
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 // member functions
0079 //
0080
0081 // ------------ method called to produce the data  ------------
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;    // NchargedPtCut bym
0097   int nChargedPtCutMR = 0;  // NchargedPtCutMR bym
0098
0099   double meanPt = 0;
0100   double meanPtMR = 0;
0101   double EtMR = 0;       // Normalized of total energy bym
0102   double TotEnergy = 0;  // Total energy bym
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();  // energy bym
0129       //float energy = (*it)->momentum().energy(); // energy bym
0130       nCharged++;
0131       meanPt += pt;
0132       // Get the total energy bym
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       // end bym
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   // Get the normalized total energy bym
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 //define this as a plug-in
0193 DEFINE_FWK_MODULE(GenHIEventProducer);