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);