File indexing completed on 2021-10-19 10:28:29
0001 #include <iostream>
0002
0003 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0004
0005
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 {}
0026
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
0077
0078 Handle<HepMCProduct> EvtHandle;
0079
0080
0081
0082 e.getByLabel("VtxSmeared", EvtHandle);
0083
0084 const HepMC::GenEvent* Evt = EvtHandle->GetEvent();
0085
0086
0087
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;
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;
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)
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);