|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
|
// -*- C++ -*-
// Class: IsoPhotonEBSelector
//
// Original Author: Silvia Taroni
// Created: Wed, 29 Nov 2017 18:23:54 GMT
//
//
#include "FWCore/PluginManager/interface/ModuleDef.h"
// system include files
#include <algorithm>
#include <iostream>
#include <memory>
#include <string>
#include <vector>
// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/global/EDFilter.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "DataFormats/EgammaCandidates/interface/Photon.h"
#include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
//
#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "CommonTools/UtilAlgos/interface/SingleObjectSelector.h"
#include "DataFormats/Common/interface/View.h"
using namespace std;
using namespace reco;
namespace edm {
class EventSetup;
}
class IsoPhotonEBSelector {
public:
IsoPhotonEBSelector(const edm::ParameterSet&, edm::ConsumesCollector& iC);
static void fillPSetDescription(edm::ParameterSetDescription& desc);
bool operator()(const reco::Photon&) const;
void newEvent(const edm::Event&, const edm::EventSetup&);
const float getEffectiveArea(float eta) const;
void printEffectiveAreas() const;
const edm::EDGetTokenT<double> theRhoToken_;
edm::Handle<double> rhoHandle_;
const std::vector<double> absEtaMin_; // low limit of the eta range
const std::vector<double> absEtaMax_; // upper limit of the eta range
const std::vector<double> effectiveAreaValues_; // effective area for this eta range
const edm::ParameterSet phIDWP_;
const vector<double> sigmaIEtaIEtaCut_;
const vector<double> hOverECut_;
const vector<double> relCombIso_;
};
void IsoPhotonEBSelector::printEffectiveAreas() const {
printf(" eta_min eta_max effective area\n");
uint nEtaBins = absEtaMin_.size();
for (uint iEta = 0; iEta < nEtaBins; iEta++) {
printf(" %8.4f %8.4f %8.5f\n", absEtaMin_[iEta], absEtaMax_[iEta], effectiveAreaValues_[iEta]);
}
}
const float IsoPhotonEBSelector::getEffectiveArea(float eta) const {
float effArea = 0;
uint nEtaBins = absEtaMin_.size();
for (uint iEta = 0; iEta < nEtaBins; iEta++) {
if (std::abs(eta) >= absEtaMin_[iEta] && std::abs(eta) < absEtaMax_[iEta]) {
effArea = effectiveAreaValues_[iEta];
break;
}
}
return effArea;
}
void IsoPhotonEBSelector::fillPSetDescription(edm::ParameterSetDescription& desc) {
desc.add<edm::InputTag>("rho", edm::InputTag("fixedGridRhoFastjetCentralCalo"));
desc.add<std::vector<double>>("absEtaMin", {0.0000, 1.0000, 1.4790, 2.0000, 2.2000, 2.3000, 2.4000});
desc.add<std::vector<double>>("absEtaMax", {1.0000, 1.4790, 2.0000, 2.2000, 2.3000, 2.4000, 5.0000});
desc.add<std::vector<double>>("effectiveAreaValues", {0.1703, 0.1715, 0.1213, 0.1230, 0.1635, 0.1937, 0.2393});
// Define the photonIDWP PSet
edm::ParameterSetDescription photonIDWP;
photonIDWP.add<std::vector<double>>("full5x5_sigmaIEtaIEtaCut", {0.011, -1.0});
photonIDWP.add<std::vector<double>>("hOverECut", {0.1, -1.0});
photonIDWP.add<std::vector<double>>("relCombIsolationWithEACut", {0.05, -1.0});
// Add the PSet to the main description
desc.add<edm::ParameterSetDescription>("phID", photonIDWP);
}
IsoPhotonEBSelector::IsoPhotonEBSelector(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
: theRhoToken_(iC.consumes<double>(cfg.getParameter<edm::InputTag>("rho"))),
absEtaMin_(cfg.getParameter<std::vector<double>>("absEtaMin")),
absEtaMax_(cfg.getParameter<std::vector<double>>("absEtaMax")),
effectiveAreaValues_(cfg.getParameter<std::vector<double>>("effectiveAreaValues")),
phIDWP_(cfg.getParameter<edm::ParameterSet>("phID")),
sigmaIEtaIEtaCut_(phIDWP_.getParameter<std::vector<double>>("full5x5_sigmaIEtaIEtaCut")),
hOverECut_(phIDWP_.getParameter<std::vector<double>>("hOverECut")),
relCombIso_(phIDWP_.getParameter<std::vector<double>>("relCombIsolationWithEACut")) {
//printEffectiveAreas();
}
void IsoPhotonEBSelector::newEvent(const edm::Event& ev, const edm::EventSetup&) {
ev.getByToken(theRhoToken_, rhoHandle_);
}
bool IsoPhotonEBSelector::operator()(const reco::Photon& ph) const {
float pt_e = ph.pt();
unsigned int ind = 0;
float abseta = fabs((ph.superCluster().get())->position().eta());
if (ph.isEB()) {
if (abseta > 1.479)
return false; // check if it is really needed
}
if (ph.isEE()) {
ind = 1;
if (abseta <= 1.479 || abseta >= 2.5)
return false; // check if it is really needed
}
if (ph.full5x5_sigmaIetaIeta() > sigmaIEtaIEtaCut_[ind])
return false;
if (ph.hadronicOverEm() > hOverECut_[ind])
return false;
const float eA = getEffectiveArea(abseta);
const float rho = rhoHandle_.isValid() ? (float)(*rhoHandle_.product()) : 0;
if ((ph.getPflowIsolationVariables().chargedHadronIso +
std::max(float(0.0),
ph.getPflowIsolationVariables().neutralHadronIso + ph.getPflowIsolationVariables().photonIso -
eA * rho)) > relCombIso_[ind] * pt_e)
return false;
return true;
}
EVENTSETUP_STD_INIT(IsoPhotonEBSelector);
typedef SingleObjectSelector<edm::View<reco::Photon>, IsoPhotonEBSelector> IsoPhotonEBSelectorAndSkim;
DEFINE_FWK_MODULE(IsoPhotonEBSelectorAndSkim);
|