Pythia6Interface/test/BasicGenTester.cc

0001 #include <iostream>
0002
0003 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0004 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0005 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0006
0007 // essentials !!!
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "DataFormats/Common/interface/Handle.h"
0010 #include "FWCore/Framework/interface/ESHandle.h"
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0016 #include "TH1.h"
0017
0018 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0019
0020 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0021 #include "HepPDT/ParticleDataTable.hh"
0022
0023 class BasicGenTester : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0024 public:
0025   //
0026   explicit BasicGenTester(const edm::ParameterSet&);
0027   virtual ~BasicGenTester() {}  // no need to delete ROOT stuff
0028                                 // as it'll be deleted upon closing TFile
0029
0030   virtual void analyze(const edm::Event&, const edm::EventSetup&);
0031   virtual void beginJob();
0032   virtual void beginRun(const edm::Run&, const edm::EventSetup&);
0033   virtual void endRun(const edm::Run&, const edm::EventSetup&);
0034   virtual void endJob();
0035
0036 private:
0037   TH1D* fNChgPartFinalState;
0038   TH1D* fNNeuPartFinalState;
0039   TH1D* fNPartFinalState;
0040   TH1D* fPtChgPartFinalState;
0041   TH1D* fPtNeuPartFinalState;
0042   TH1D* fEtaChgPartFinalState;
0043   TH1D* fEtaNeuPartFinalState;
0044   int fNPart;
0045   double fPtMin;
0046   double fPtMax;
0047   edm::ESHandle<HepPDT::ParticleDataTable> fPDGTable;
0048   edm::ESGetToken<HepPDT::ParticleDataTable, edm::DefaultRecord> fPDGTableToken;
0049 };
0050
0051 using namespace edm;
0052 using namespace std;
0053
0054 BasicGenTester::BasicGenTester(const ParameterSet& pset)
0055     : fNChgPartFinalState(0),
0056       fNNeuPartFinalState(0),
0057       fNPartFinalState(0),
0058       fPtChgPartFinalState(0),
0059       fPtNeuPartFinalState(0) {
0060   fNPart = pset.getUntrackedParameter<int>("NPartForHisto", 500);
0061   fPtMin = pset.getUntrackedParameter<double>("PtMinForHisto", 0.);
0062   fPtMax = pset.getUntrackedParameter<double>("PtMaxForHisto", 25.);
0063   usesResource(TFileService::kSharedResource);
0064   fPDGTableToken = esConsumes<edm::Transition::BeginRun>();
0065   consumes<HepMCProduct>(edm::InputTag("VtxSmeared"));
0066 }
0067
0068 void BasicGenTester::beginJob() {
0069   Service<TFileService> fs;
0070
0071   fNChgPartFinalState =
0072       fs->make<TH1D>("NChgPartFinalState", "Number of final state charged particles", fNPart, 0., (double)fNPart);
0073   fNNeuPartFinalState =
0074       fs->make<TH1D>("NNeuPartFinalState", "Number of final state neutral particles", fNPart, 0., (double)fNPart);
0075   fNPartFinalState =
0076       fs->make<TH1D>("NPartFinalState", "Total number of final state particles", fNPart, 0., (double)fNPart);
0077   fPtChgPartFinalState =
0078       fs->make<TH1D>("PtChgPartFinalState", "Pt of final state charged particles", 500, fPtMin, fPtMax);
0079   fPtNeuPartFinalState =
0080       fs->make<TH1D>("PtNeuPartFinalState", "Pt of final state neutral particles", 500, fPtMin, fPtMax);
0081   fEtaChgPartFinalState =
0082       fs->make<TH1D>("EtaChgPartFinalState", "Eta of final state charged particles", 100, -5.0, 5.0);
0083   fEtaNeuPartFinalState =
0084       fs->make<TH1D>("EtaNeuPartFinalState", "Eta of final state neutral particles", 100, -5.0, 5.0);
0085   return;
0086 }
0087
0088 void BasicGenTester::beginRun(const edm::Run& r, const edm::EventSetup& es) {
0089   fPDGTable = es.getHandle(fPDGTableToken);
0090
0091   return;
0092 }
0093
0094 void BasicGenTester::analyze(const Event& e, const EventSetup&) {
0095   // here's an example of accessing GenEventInfoProduct
0096   /*
0097    Handle< GenEventInfoProduct > GenInfoHandle;
0098    e.getByLabel( "generator", GenInfoHandle );
0099    double qScale = GenInfoHandle->qScale();
0100    double pthat = ( GenInfoHandle->hasBinningValues() ?
0101                   (GenInfoHandle->binningValues())[0] : 0.0);
0102    cout << " qScale = " << qScale << " pthat = " << pthat << endl;
0103    double evt_weight1 = GenInfoHandle->weights()[0]; // this is "stanrd Py6 evt weight;
0104                                                      // corresponds to PYINT1/VINT(97)
0105    double evt_weight2 = GenInfoHandle->weights()[1]; // in case you run in CSA mode or otherwise
0106                                                      // use PYEVWT routine, this will be weight
0107                              // as returned by PYEVWT, i.e. PYINT1/VINT(99)
0108    //std::cout << " evt_weight1 = " << evt_weight1 << std::endl;
0109    //std::cout << " evt_weight2 = " << evt_weight2 << std::endl;
0110    double weight = GenInfoHandle->weight();
0111    //std::cout << " as returned by the weight() method, integrated event weight = " << weight << std::endl;
0112 */
0113
0114   // here's an example of accessing particles in the event record (HepMCProduct)
0115   //
0116   Handle<HepMCProduct> EvtHandle;
0117
0118   // find initial (unsmeared, unfiltered,...) HepMCProduct
0119   //
0120   e.getByLabel("VtxSmeared", EvtHandle);
0121
0122   const HepMC::GenEvent* Evt = EvtHandle->GetEvent();
0123
0124   double NChgPartFS = 0;
0125   double NNeuPartFS = 0;
0126   for (HepMC::GenEvent::particle_const_iterator part = Evt->particles_begin(); part != Evt->particles_end(); ++part) {
0127     /*
0128       int pid = (*part)->pdg_id();
0129       if ( abs(pid) == 15 )
0130       {
0131          std::cout << "found tau " << std::endl;
0132      int stat = (*part)->status();
0133      if ( (*part)->end_vertex() )
0134      {
0135         (*part)->end_vertex()->print();
0136         std::cout << "done with tau end vertex " << std::endl;
0137      }
0138      std::cout << " end looking at tau" << std::endl;
0139       }
0140 */
0141
0142     if ((*part)->status() == 1 && !((*part)->end_vertex())) {
0143       int PartID = (*part)->pdg_id();
0144
0145       const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0146       double charge = PData->charge();
0147       if (fabs(charge) != 0.) {
0148         NChgPartFS++;
0149         fPtChgPartFinalState->Fill(((*part)->momentum()).perp());
0150         fEtaChgPartFinalState->Fill(((*part)->momentum()).eta());
0151       } else {
0152         NNeuPartFS++;
0153         fPtNeuPartFinalState->Fill(((*part)->momentum()).perp());
0154         fEtaNeuPartFinalState->Fill(((*part)->momentum()).perp());
0155       }
0156     }
0157   }
0158
0159   fNChgPartFinalState->Fill(NChgPartFS);
0160   fNNeuPartFinalState->Fill(NNeuPartFS);
0161   fNPartFinalState->Fill((NChgPartFS + NNeuPartFS));
0162
0163   return;
0164 }
0165
0166 void BasicGenTester::endRun(const edm::Run& r, const edm::EventSetup&) { return; }
0167
0168 void BasicGenTester::endJob() { return; }
0169
0170 DEFINE_FWK_MODULE(BasicGenTester);