Pythia6Interface/test/TauPhotonTester.cc

0001 #include <iostream>
0002
0003 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0004
0005 // essentials !!!
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "FWCore/Framework/interface/ESHandle.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0014 #include "TH1.h"
0015
0016 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0017
0018 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0019 #include "HepPDT/ParticleDataTable.hh"
0020
0021 class TauPhotonTester : public edm::one::EDAnalyzer<edm::one::SharedResources, edm::one::WatchRuns> {
0022 public:
0023   //
0024   explicit TauPhotonTester(const edm::ParameterSet&);
0025   ~TauPhotonTester() override {}  // no need to delete ROOT stuff
0026                                   // as it'll be deleted upon closing TFile
0027
0028   void analyze(const edm::Event&, const edm::EventSetup&) override;
0029   void beginJob() override;
0030   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0031   void endRun(const edm::Run&, const edm::EventSetup&) override {}
0032   void endJob() override;
0033
0034 private:
0035   int fTauPhotonCounter;
0036   int fTauDecPhotonCounter;
0037   int fTauPhotonVtxCounter;
0038   TH1D* fPhotonEnergy;
0039   TH1D* fPhotonPt;
0040   TH1D* fPhotonEnergyGt10MeV;
0041   TH1D* fPhotonPtGt10MeV;
0042   edm::ESHandle<HepPDT::ParticleDataTable> fPDGTable;
0043   edm::ESGetToken<HepPDT::ParticleDataTable, edm::DefaultRecord> fPDGTableToken;
0044   int fVerbosity;
0045 };
0046
0047 using namespace edm;
0048 using namespace std;
0049
0050 TauPhotonTester::TauPhotonTester(const ParameterSet& pset)
0051     : fTauPhotonCounter(0), fTauDecPhotonCounter(0), fTauPhotonVtxCounter(0), fVerbosity(0) {
0052   fVerbosity = pset.getUntrackedParameter<int>("Verbosity", 0);
0053   fPDGTableToken = esConsumes<edm::Transition::BeginRun>();
0054   consumes<HepMCProduct>(edm::InputTag("VtxSmeared"));
0055   usesResource(TFileService::kSharedResource);
0056 }
0057
0058 void TauPhotonTester::beginJob() {
0059   Service<TFileService> fs;
0060   fPhotonEnergy = fs->make<TH1D>("PhotonEnergy", "Energy of the Brem Photon (off tau)", 100, 0.0, 2.0);
0061   fPhotonPt = fs->make<TH1D>("PhotonPt", "Pt of the Brem Photon (off tau)", 100, 0.0, 2.0);
0062   fPhotonEnergyGt10MeV =
0063       fs->make<TH1D>("PhotonEnergyGt10MeV", "Energy of the Brem Photon (off tau), Pt>10MeV", 100, 0.01, 2.01);
0064   fPhotonPtGt10MeV = fs->make<TH1D>("PhotonPtGt10MeV", "Pt of the Brem Photon (off tau), Pt>10MeV", 100, 0.01, 2.01);
0065
0066   return;
0067 }
0068
0069 void TauPhotonTester::beginRun(const edm::Run& r, const edm::EventSetup& es) {
0070   fPDGTable = es.getHandle(fPDGTableToken);
0071
0072   return;
0073 }
0074
0075 void TauPhotonTester::analyze(const Event& e, const EventSetup&) {
0076   // here's an example of accessing particles in the event record (HepMCProduct)
0077   //
0078   Handle<HepMCProduct> EvtHandle;
0079
0080   // find initial (unsmeared, unfiltered,...) HepMCProduct
0081   //
0082   e.getByLabel("VtxSmeared", EvtHandle);
0083
0084   const HepMC::GenEvent* Evt = EvtHandle->GetEvent();
0085
0086   // loop over, find a vertex with outgoing tau of no-null vtx
0087   // see if there's also a photon
0088   //
0089   bool tau_in_vtx_flag = false;
0090   bool tau_out_vtx_flag = false;
0091
0092   vector<double> PhotonE;
0093   vector<double> PhotonPt;
0094
0095   for (HepMC::GenEvent::vertex_const_iterator vtx = Evt->vertices_begin(); vtx != Evt->vertices_end(); ++vtx) {
0096     tau_in_vtx_flag = false;
0097     tau_out_vtx_flag = false;
0098
0099     if ((*vtx)->particles_in_size() == 1 && abs((*((*vtx)->particles_in_const_begin()))->pdg_id()) == 15) {
0100       tau_in_vtx_flag = true;
0101     }
0102
0103     if ((*vtx)->particles_out_size() == 1 && abs((*((*vtx)->particles_out_const_begin()))->pdg_id()) == 15 &&
0104         tau_in_vtx_flag) {
0105       continue;  // "intermediate" tau vtx, skip it
0106     }
0107
0108     int NGammaInVtx = 0;
0109     int NTauInVtx = 0;
0110     PhotonE.clear();
0111     PhotonPt.clear();
0112
0113     for (HepMC::GenVertex::particle_iterator pitr = (*vtx)->particles_begin(HepMC::children);
0114          pitr != (*vtx)->particles_end(HepMC::children);
0115          ++pitr) {
0116       if (abs((*pitr)->pdg_id()) == 15) {
0117         NTauInVtx++;
0118         tau_out_vtx_flag = true;
0119       }
0120       if ((*pitr)->pdg_id() == 22) {
0121         NGammaInVtx++;
0122         double ee = (*pitr)->momentum().t();
0123         double px = (*pitr)->momentum().x();
0124         double py = (*pitr)->momentum().y();
0125         double pt = sqrt(px * px + py * py);
0126         PhotonE.push_back(ee);
0127         PhotonPt.push_back(pt);
0128       }
0129     }
0130
0131     if (!tau_in_vtx_flag && !tau_out_vtx_flag)
0132       continue;  // no-tau vtx
0133
0134     if (tau_in_vtx_flag && tau_out_vtx_flag && NGammaInVtx > 0) {
0135       if (fVerbosity > 0)
0136         (*vtx)->print();
0137       fTauPhotonVtxCounter++;
0138       fTauPhotonCounter += NGammaInVtx;
0139       for (int i1 = 0; i1 < NGammaInVtx; i1++) {
0140         fPhotonEnergy->Fill(PhotonE[i1]);
0141         fPhotonPt->Fill(PhotonPt[i1]);
0142         if (PhotonPt[i1] > 0.01) {
0143           fPhotonEnergyGt10MeV->Fill(PhotonE[i1]);
0144           fPhotonPtGt10MeV->Fill(PhotonPt[i1]);
0145         }
0146       }
0147     }
0148
0149     if (!tau_in_vtx_flag && tau_out_vtx_flag) {
0150       if (fVerbosity > 0 && NGammaInVtx > 0)
0151         (*vtx)->print();
0152       if (NGammaInVtx > 0) {
0153         fTauPhotonCounter += NGammaInVtx;
0154         for (int i2 = 0; i2 < NGammaInVtx; i2++) {
0155           fPhotonEnergy->Fill(PhotonE[i2]);
0156           fPhotonPt->Fill(PhotonPt[i2]);
0157           if (PhotonPt[i2] > 0.01) {
0158             fPhotonEnergyGt10MeV->Fill(PhotonE[i2]);
0159             fPhotonPtGt10MeV->Fill(PhotonPt[i2]);
0160           }
0161         }
0162       }
0163     }
0164
0165     if (tau_in_vtx_flag && !tau_out_vtx_flag)  // real tau decay vtx
0166     {
0167       if (fVerbosity > 0 && NGammaInVtx > 0)
0168         (*vtx)->print();
0169       if (NGammaInVtx > 0) {
0170         fTauDecPhotonCounter += NGammaInVtx;
0171         for (int i3 = 0; i3 < NGammaInVtx; i3++) {
0172           fPhotonEnergy->Fill(PhotonE[i3]);
0173           fPhotonPt->Fill(PhotonPt[i3]);
0174           if (PhotonPt[i3] > 0.01) {
0175             fPhotonEnergyGt10MeV->Fill(PhotonE[i3]);
0176             fPhotonPtGt10MeV->Fill(PhotonPt[i3]);
0177           }
0178         }
0179       }
0180     }
0181   }
0182
0183   return;
0184 }
0185
0186 void TauPhotonTester::endJob() {
0187   std::cout << "Tau-Photon Vtx Counter: " << fTauPhotonVtxCounter << std::endl;
0188   std::cout << "Tau-Photon Counter: " << fTauPhotonCounter << std::endl;
0189   std::cout << "TauDec-Photon Counter: " << fTauDecPhotonCounter << std::endl;
0190
0191   return;
0192 }
0193
0194 DEFINE_FWK_MODULE(TauPhotonTester);