RecoHI/plugins/HiBasicGenTest.cc

0001 #include "HepMC/GenEvent.h"
0002 #include "HepMC/GenParticle.h"
0003 #include "Validation/RecoHI/plugins/HiBasicGenTest.h"
0004 #include "DataFormats/Candidate/interface/Candidate.h"
0005 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0006 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0007
0008 #include "HepMC/GenEvent.h"
0009 #include "HepMC/HeavyIon.h"
0010
0011 #include <TString.h>
0012 #include <TMath.h>
0013
0014 using namespace edm;
0015 using namespace HepMC;
0016
0017 HiBasicGenTest::HiBasicGenTest(const edm::ParameterSet& iPSet) {
0018   generatorToken_ = consumes<edm::HepMCProduct>(iPSet.getParameter<edm::InputTag>("generatorLabel"));
0019   pdtToken_ = esConsumes<edm::Transition::BeginRun>();
0020 }
0021
0022 HiBasicGenTest::~HiBasicGenTest() {}
0023
0024 void HiBasicGenTest::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const&) {
0025   ///Setting the DQM top directories
0026   ibooker.setCurrentFolder("Generator/Particles");
0027
0028   ///Booking the ME's
0029   for (int ibin = 0; ibin < 3; ibin++) {
0030     dnchdeta[ibin] = ibooker.book1D(Form("dnchdeta%d", ibin), ";#eta;dN^{ch}/d#eta", 100, -6.0, 6.0);
0031
0032     dnchdpt[ibin] = ibooker.book1D(Form("dnchdpt%d", ibin), ";p_{T};dN^{ch}/dp_{T}", 200, 0.0, 100.0);
0033
0034     b[ibin] = ibooker.book1D(Form("b%d", ibin), ";b[fm];events", 100, 0.0, 20.0);
0035     dnchdphi[ibin] = ibooker.book1D(Form("dnchdphi%d", ibin), ";#phi;dN^{ch}/d#phi", 100, -3.2, 3.2);
0036   }
0037
0038   rp = ibooker.book1D("phi0", ";#phi_{RP};events", 100, -3.2, 3.2);
0039 }
0040
0041 void HiBasicGenTest::dqmBeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0042   pdt = iSetup.getHandle(pdtToken_);
0043 }
0044
0045 void HiBasicGenTest::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0046   Handle<HepMCProduct> mc;
0047   iEvent.getByToken(generatorToken_, mc);
0048   const HepMC::GenEvent* evt = mc->GetEvent();
0049   const HepMC::HeavyIon* hi = evt->heavy_ion();
0050
0051   int cbin = 0;
0052   double phi0 = 0.;
0053
0054   if (hi) {
0055     double ip = hi->impact_parameter();
0056     phi0 = hi->event_plane_angle();
0057
0058     // fill reaction plane distribution
0059     rp->Fill(phi0);
0060
0061     // if the event is in one of the centrality bins of interest fill hists
0062     int cbin = -1;
0063     if (ip < 5.045)
0064       cbin = 0;
0065     else if (ip < 7.145 && ip > 5.045)
0066       cbin = 1;
0067     else if (ip < 15.202 && ip > 14.283)
0068       cbin = 2;
0069     if (cbin < 0)
0070       return;
0071
0072     // fill impact parameter distributions
0073     b[cbin]->Fill(ip);
0074   }
0075
0076   // loop over particles
0077   HepMC::GenEvent::particle_const_iterator begin = evt->particles_begin();
0078   HepMC::GenEvent::particle_const_iterator end = evt->particles_end();
0079   for (HepMC::GenEvent::particle_const_iterator it = begin; it != end; ++it) {
0080     // only fill hists for status=1 particles
0081     if ((*it)->status() != 1)
0082       continue;
0083
0084     // only fill hists for charged particles
0085     int pdg_id = (*it)->pdg_id();
0086     const ParticleData* part = pdt->particle(pdg_id);
0087     int charge = static_cast<int>(part->charge());
0088     if (charge == 0)
0089       continue;
0090
0091     float eta = (*it)->momentum().eta();
0092     float phi = (*it)->momentum().phi();
0093     float pt = (*it)->momentum().perp();
0094
0095     dnchdeta[cbin]->Fill(eta);
0096     dnchdpt[cbin]->Fill(pt);
0097
0098     double pi = TMath::Pi();
0099     double p = phi - phi0;
0100     if (p > pi)
0101       p = p - 2 * pi;
0102     if (p < -1 * pi)
0103       p = p + 2 * pi;
0104     dnchdphi[cbin]->Fill(p);
0105   }
0106
0107   return;
0108 }
0109
0110 #include "FWCore/Framework/interface/MakerMacros.h"
0111 DEFINE_FWK_MODULE(HiBasicGenTest);