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     PhotonE.clear();
0110     PhotonPt.clear();
0111
0112     for (HepMC::GenVertex::particle_iterator pitr = (*vtx)->particles_begin(HepMC::children);
0113          pitr != (*vtx)->particles_end(HepMC::children);
0114          ++pitr) {
0115       if (abs((*pitr)->pdg_id()) == 15) {
0116         tau_out_vtx_flag = true;
0117       }
0118       if ((*pitr)->pdg_id() == 22) {
0119         NGammaInVtx++;
0120         double ee = (*pitr)->momentum().t();
0121         double px = (*pitr)->momentum().x();
0122         double py = (*pitr)->momentum().y();
0123         double pt = sqrt(px * px + py * py);
0124         PhotonE.push_back(ee);
0125         PhotonPt.push_back(pt);
0126       }
0127     }
0128
0129     if (!tau_in_vtx_flag && !tau_out_vtx_flag)
0130       continue;  // no-tau vtx
0131
0132     if (tau_in_vtx_flag && tau_out_vtx_flag && NGammaInVtx > 0) {
0133       if (fVerbosity > 0)
0134         (*vtx)->print();
0135       fTauPhotonVtxCounter++;
0136       fTauPhotonCounter += NGammaInVtx;
0137       for (int i1 = 0; i1 < NGammaInVtx; i1++) {
0138         fPhotonEnergy->Fill(PhotonE[i1]);
0139         fPhotonPt->Fill(PhotonPt[i1]);
0140         if (PhotonPt[i1] > 0.01) {
0141           fPhotonEnergyGt10MeV->Fill(PhotonE[i1]);
0142           fPhotonPtGt10MeV->Fill(PhotonPt[i1]);
0143         }
0144       }
0145     }
0146
0147     if (!tau_in_vtx_flag && tau_out_vtx_flag) {
0148       if (fVerbosity > 0 && NGammaInVtx > 0)
0149         (*vtx)->print();
0150       if (NGammaInVtx > 0) {
0151         fTauPhotonCounter += NGammaInVtx;
0152         for (int i2 = 0; i2 < NGammaInVtx; i2++) {
0153           fPhotonEnergy->Fill(PhotonE[i2]);
0154           fPhotonPt->Fill(PhotonPt[i2]);
0155           if (PhotonPt[i2] > 0.01) {
0156             fPhotonEnergyGt10MeV->Fill(PhotonE[i2]);
0157             fPhotonPtGt10MeV->Fill(PhotonPt[i2]);
0158           }
0159         }
0160       }
0161     }
0162
0163     if (tau_in_vtx_flag && !tau_out_vtx_flag)  // real tau decay vtx
0164     {
0165       if (fVerbosity > 0 && NGammaInVtx > 0)
0166         (*vtx)->print();
0167       if (NGammaInVtx > 0) {
0168         fTauDecPhotonCounter += NGammaInVtx;
0169         for (int i3 = 0; i3 < NGammaInVtx; i3++) {
0170           fPhotonEnergy->Fill(PhotonE[i3]);
0171           fPhotonPt->Fill(PhotonPt[i3]);
0172           if (PhotonPt[i3] > 0.01) {
0173             fPhotonEnergyGt10MeV->Fill(PhotonE[i3]);
0174             fPhotonPtGt10MeV->Fill(PhotonPt[i3]);
0175           }
0176         }
0177       }
0178     }
0179   }
0180
0181   return;
0182 }
0183
0184 void TauPhotonTester::endJob() {
0185   std::cout << "Tau-Photon Vtx Counter: " << fTauPhotonVtxCounter << std::endl;
0186   std::cout << "Tau-Photon Counter: " << fTauPhotonCounter << std::endl;
0187   std::cout << "TauDec-Photon Counter: " << fTauDecPhotonCounter << std::endl;
0188
0189   return;
0190 }
0191
0192 DEFINE_FWK_MODULE(TauPhotonTester);