GenHIEventProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    GenHIEventProducer
// Class:      GenHIEventProducer
//
//
// Original Author:  Yetkin Yilmaz
//         Created:  Thu Aug 13 08:39:51 EDT 2009
//
//
 
// system include files
#include <memory>
#include <string>
#include <iostream>
 
// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDProducer.h"
#include "FWCore/Framework/interface/EventSetup.h"
 
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
 
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"
 
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
#include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
 
#include "HepMC/HeavyIon.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
using namespace std;
 
//
// class decleration
//
 
class GenHIEventProducer : public edm::EDProducer {
public:
  explicit GenHIEventProducer(const edm::ParameterSet&);
  ~GenHIEventProducer() override;
 
private:
  void produce(edm::Event&, const edm::EventSetup&) override;
  edm::EDGetTokenT<CrossingFrame<edm::HepMCProduct> > hepmcSrc_;
  edm::ESHandle<ParticleDataTable> pdt;
 
  double ptCut_;
  bool doParticleInfo_;
};
 
//
// constants, enums and typedefs
//
 
//
// static data member definitions
//
 
//
// constructors and destructor
//
GenHIEventProducer::GenHIEventProducer(const edm::ParameterSet& iConfig) {
  produces<edm::GenHIEvent>();
  hepmcSrc_ = consumes<CrossingFrame<edm::HepMCProduct> >(iConfig.getParameter<edm::InputTag>("src"));
  doParticleInfo_ = iConfig.getUntrackedParameter<bool>("doParticleInfo", false);
  if (doParticleInfo_) {
    ptCut_ = iConfig.getParameter<double>("ptCut");
  }
}
 
GenHIEventProducer::~GenHIEventProducer() {
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
}
 
//
// member functions
//
 
// ------------ method called to produce the data  ------------
void GenHIEventProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  using namespace edm;
 
  if (!(pdt.isValid()))
    iSetup.getData(pdt);
 
  double b = -1;
  int npart = -1;
  int ncoll = 0;
  int nhard = 0;
  double phi = 0;
  double ecc = -1;
 
  int nCharged = 0;
  int nChargedMR = 0;
  int nChargedPtCut = 0;    // NchargedPtCut bym
  int nChargedPtCutMR = 0;  // NchargedPtCutMR bym
 
  double meanPt = 0;
  double meanPtMR = 0;
  double EtMR = 0;       // Normalized of total energy bym
  double TotEnergy = 0;  // Total energy bym
 
  Handle<CrossingFrame<edm::HepMCProduct> > hepmc;
  iEvent.getByToken(hepmcSrc_, hepmc);
  MixCollection<HepMCProduct> mix(hepmc.product());
 
  if (mix.size() < 1) {
    throw cms::Exception("MatchVtx") << "Mixing has " << mix.size() << " sub-events, should have been at least 1"
                                     << endl;
  }
 
  const HepMCProduct& hievt = mix.getObject(mix.size() - 1);
  const HepMC::GenEvent* evt = hievt.GetEvent();
  if (doParticleInfo_) {
    HepMC::GenEvent::particle_const_iterator begin = evt->particles_begin();
    HepMC::GenEvent::particle_const_iterator end = evt->particles_end();
    for (HepMC::GenEvent::particle_const_iterator it = begin; it != end; ++it) {
      if ((*it)->status() != 1)
        continue;
      int pdg_id = (*it)->pdg_id();
      const ParticleData* part = pdt->particle(pdg_id);
      int charge = static_cast<int>(part->charge());
 
      if (charge == 0)
        continue;
      float pt = (*it)->momentum().perp();
      float eta = (*it)->momentum().eta();
      float energy = (*it)->momentum().e();  // energy bym
      //float energy = (*it)->momentum().energy(); // energy bym
      nCharged++;
      meanPt += pt;
      // Get the total energy bym
      if (fabs(eta) < 1.0) {
        TotEnergy += energy;
      }
      if (pt > ptCut_) {
        nChargedPtCut++;
        if (fabs(eta) < 0.5) {
          nChargedPtCutMR++;
        }
      }
      // end bym
 
      if (fabs(eta) > 0.5)
        continue;
      nChargedMR++;
      meanPtMR += pt;
    }
  }
  const HepMC::HeavyIon* hi = evt->heavy_ion();
 
  if (hi) {
    ncoll = ncoll + hi->Ncoll();
    nhard = nhard + hi->Ncoll_hard();
    int np = hi->Npart_proj() + hi->Npart_targ();
    if (np >= 0) {
      npart = np;
      b = hi->impact_parameter();
      phi = hi->event_plane_angle();
      ecc = hi->eccentricity();
    }
  }
 
  // Get the normalized total energy bym
  if (TotEnergy != 0) {
    EtMR = TotEnergy / 2;
  }
 
  if (nChargedMR != 0) {
    meanPtMR /= nChargedMR;
  }
  if (nCharged != 0) {
    meanPt /= nCharged;
  }
 
  std::unique_ptr<edm::GenHIEvent> pGenHI(new edm::GenHIEvent(
      b, npart, ncoll, nhard, phi, ecc, nCharged, nChargedMR, meanPt, meanPtMR, EtMR, nChargedPtCut, nChargedPtCutMR));
 
  iEvent.put(std::move(pGenHI));
}
 
//define this as a plug-in
DEFINE_FWK_MODULE(GenHIEventProducer);