PyquenInterface/test/PyquenAnalyzer.cc

0001
0002 #include <iostream>
0003
0004 #include "PyquenAnalyzer.h"
0005
0006 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0007 #include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
0008 #include "DataFormats/Common/interface/Handle.h"
0009
0010 // essentials !!!
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "HepMC/HeavyIon.h"
0014
0015 #include "FWCore/ServiceRegistry/interface/Service.h"
0016 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0017
0018 #include "TFile.h"
0019 #include "TH1.h"
0020
0021 using namespace edm;
0022 using namespace std;
0023
0024 /************************************************************************
0025  *
0026  *  This code is analyzing pyquen events produced with pyquen+analysis.cfg
0027  *
0028  *  Author: cmironov@lanl.gov
0029  *
0030  *************************************************************************/
0031
0032 PyquenAnalyzer::PyquenAnalyzer(const ParameterSet& pset) : phdNdEta(0), phdNdY(0), phdNdPt(0), phdNdPhi(0) {
0033   // constructor
0034   srcT_ = mayConsume<HepMCProduct>(edm::InputTag("generator", "unsmeared"));
0035 }
0036
0037 //_______________________________________________________________________
0038 void PyquenAnalyzer::beginJob() {
0039   //runs at the begining of the job
0040   edm::Service<TFileService> fs;
0041   TH1::SetDefaultSumw2(true);
0042
0043   phdNdEta = fs->make<TH1D>("phdNdEta", ";#eta;", 100, -10., 10.);
0044   phdNdY = fs->make<TH1D>("phdNdY", ";y;", 100, -10., 10.);
0045   phdNdPt = fs->make<TH1D>("phdNdPt", ";p_{T}(GeV/c);", 100, 0., 10.);
0046   phdNdPhi = fs->make<TH1D>("phdNdPhi", ";d#phi(rad);", 100, -3.15, 3.15);
0047
0048   return;
0049 }
0050
0051 //______________________________________________________________________
0052 void PyquenAnalyzer::analyze(const Event& e, const EventSetup&) {
0053   //runs every event
0054
0055   Handle<HepMCProduct> EvtHandle;
0056   e.getByToken(srcT_, EvtHandle);
0057
0058   //  EvtHandle->GetEvent()->print();
0059
0060   double part_eta, part_y, part_pt, part_phi, part_e, part_pz;
0061   const HepMC::GenEvent* myEvt = EvtHandle->GetEvent();
0062   //HepMC::GenVertex* genvtx = nullptr;
0063   //genvtx = myEvt->signal_process_vertex();
0064   //cout<<" HEPMC Vertex X "<<genvtx->position().x()<<endl;
0065
0066   for (HepMC::GenEvent::particle_const_iterator p = myEvt->particles_begin(); p != myEvt->particles_end(); p++) {
0067     if (!(*p)->end_vertex() && abs((*p)->pdg_id()) == 211) {
0068       part_eta = (*p)->momentum().eta();
0069       part_e = (*p)->momentum().e();
0070       part_pt = (*p)->momentum().perp();
0071       part_phi = (*p)->momentum().phi();
0072       part_pz = (*p)->momentum().z();
0073       part_y = 0.5 * log((part_e + part_pz) / (part_e - part_pz));
0074
0075       phdNdEta->Fill(part_eta);
0076       phdNdY->Fill(part_y);
0077       phdNdPt->Fill(part_pt);
0078       phdNdPhi->Fill(part_phi);
0079     }
0080   }
0081   return;
0082 }
0083
0084 //_____________________________________________________________
0085 void PyquenAnalyzer::endJob() {
0086   // executed at the end of the job
0087
0088   phdNdEta->Scale(phdNdEta->GetBinWidth(0));
0089   phdNdY->Scale(phdNdY->GetBinWidth(0));
0090   phdNdPt->Scale(phdNdPt->GetBinWidth(0));
0091   phdNdPhi->Scale(phdNdPhi->GetBinWidth(0));
0092 }
0093
0094 //define as a plug-in
0095 DEFINE_FWK_MODULE(PyquenAnalyzer);