CMS 3D CMS Logo

EgammaIsoDetIdCollectionProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: EgammaIsoDetIdCollectionProducer
4 // Class: EgammaIsoDetIdCollectionProducer
5 //
33 
34 template <class T1>
36 public:
37  typedef std::vector<T1> T1Collection;
41  void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
42 
43 private:
44  // ----------member data ---------------------------
51  double energyCut_;
52  double etCut_;
53  double etCandCut_;
54  double outerRadius_;
55  double innerRadius_;
57 
58  std::vector<int> severitiesexclEB_;
59  std::vector<int> severitiesexclEE_;
60  std::vector<int> flagsexclEB_;
61  std::vector<int> flagsexclEE_;
62 };
63 
64 template <class T1>
66  : recHitsToken_{consumes(iConfig.getParameter<edm::InputTag>("recHitsLabel"))},
67  emObjectToken_{consumes(iConfig.getParameter<edm::InputTag>("emObjectLabel"))},
68  caloGeometryToken_(esConsumes()),
69  sevLvToken_(esConsumes()),
70  recHitsLabel_(iConfig.getParameter<edm::InputTag>("recHitsLabel")),
71  emObjectLabel_(iConfig.getParameter<edm::InputTag>("emObjectLabel")),
72  energyCut_(iConfig.getParameter<double>("energyCut")),
73  etCut_(iConfig.getParameter<double>("etCut")),
74  etCandCut_(iConfig.getParameter<double>("etCandCut")),
75  outerRadius_(iConfig.getParameter<double>("outerRadius")),
76  innerRadius_(iConfig.getParameter<double>("innerRadius")),
77  interestingDetIdCollection_(iConfig.getParameter<std::string>("interestingDetIdCollection")) {
78  auto const& flagnamesEB = iConfig.getParameter<std::vector<std::string>>("RecHitFlagToBeExcludedEB");
79  auto const& flagnamesEE = iConfig.getParameter<std::vector<std::string>>("RecHitFlagToBeExcludedEE");
80 
81  flagsexclEB_ = StringToEnumValue<EcalRecHit::Flags>(flagnamesEB);
82  flagsexclEE_ = StringToEnumValue<EcalRecHit::Flags>(flagnamesEE);
83 
84  auto const& severitynamesEB = iConfig.getParameter<std::vector<std::string>>("RecHitSeverityToBeExcludedEB");
85 
86  severitiesexclEB_ = StringToEnumValue<EcalSeverityLevel::SeverityLevel>(severitynamesEB);
87 
88  auto const& severitynamesEE = iConfig.getParameter<std::vector<std::string>>("RecHitSeverityToBeExcludedEE");
89 
90  severitiesexclEE_ = StringToEnumValue<EcalSeverityLevel::SeverityLevel>(severitynamesEE);
91 
92  //register your products
93  produces<DetIdCollection>(interestingDetIdCollection_);
94 }
95 
96 // ------------ method called to produce the data ------------
97 template <class T1>
100  const edm::EventSetup& iSetup) const {
101  using namespace edm;
102  using namespace std;
103 
104  auto const& emObjects = iEvent.get(emObjectToken_);
105  auto const& ecalRecHits = iEvent.get(recHitsToken_);
106 
107  edm::ESHandle<CaloGeometry> pG = iSetup.getHandle(caloGeometryToken_);
108  const CaloGeometry* caloGeom = pG.product();
109 
110  edm::ESHandle<EcalSeverityLevelAlgo> sevlv = iSetup.getHandle(sevLvToken_);
111  const EcalSeverityLevelAlgo* sevLevel = sevlv.product();
112 
113  std::unique_ptr<CaloDualConeSelector<EcalRecHit>> doubleConeSel_ = nullptr;
114  if (recHitsLabel_.instance() == "EcalRecHitsEB") {
115  doubleConeSel_ =
116  std::make_unique<CaloDualConeSelector<EcalRecHit>>(innerRadius_, outerRadius_, &*pG, DetId::Ecal, EcalBarrel);
117  } else if (recHitsLabel_.instance() == "EcalRecHitsEE") {
118  doubleConeSel_ =
119  std::make_unique<CaloDualConeSelector<EcalRecHit>>(innerRadius_, outerRadius_, &*pG, DetId::Ecal, EcalEndcap);
120  }
121 
122  //Create empty output collections
123  auto detIdCollection = std::make_unique<DetIdCollection>();
124 
125  if (doubleConeSel_) { //if cone selector was created
126  for (auto const& emObj : emObjects) { //Loop over candidates
127 
128  if (emObj.et() < etCandCut_)
129  continue; //don't calculate if object hasn't enough energy
130 
131  GlobalPoint pclu(emObj.caloPosition().x(), emObj.caloPosition().y(), emObj.caloPosition().z());
132  doubleConeSel_->selectCallback(pclu, ecalRecHits, [&](const EcalRecHit& recIt) {
133  if (recIt.energy() < energyCut_)
134  return; //dont fill if below E noise value
135 
136  double et =
137  recIt.energy() * caloGeom->getPosition(recIt.detid()).perp() / caloGeom->getPosition(recIt.detid()).mag();
138 
139  if (et < etCut_)
140  return; //dont fill if below ET noise value
141 
142  bool isBarrel = false;
143  if (fabs(caloGeom->getPosition(recIt.detid()).eta() < 1.479))
144  isBarrel = true;
145 
146  int severityFlag = sevLevel->severityLevel(recIt.detid(), ecalRecHits);
147  if (isBarrel) {
148  auto sit = std::find(severitiesexclEB_.begin(), severitiesexclEB_.end(), severityFlag);
149  if (sit != severitiesexclEB_.end())
150  return;
151  } else {
152  auto sit = std::find(severitiesexclEE_.begin(), severitiesexclEE_.end(), severityFlag);
153  if (sit != severitiesexclEE_.end())
154  return;
155  }
156 
157  if (isBarrel) {
158  // new rechit flag checks
159  if (!recIt.checkFlag(EcalRecHit::kGood)) {
160  if (recIt.checkFlags(flagsexclEB_)) {
161  return;
162  }
163  }
164  } else {
165  // new rechit flag checks
166  if (!recIt.checkFlag(EcalRecHit::kGood)) {
167  if (recIt.checkFlags(flagsexclEE_)) {
168  return;
169  }
170  }
171  }
172 
173  if (std::find(detIdCollection->begin(), detIdCollection->end(), recIt.detid()) == detIdCollection->end())
174  detIdCollection->push_back(recIt.detid());
175  }); //end rechits
176 
177  } //end candidates
178 
179  } //end if cone selector was created
180 
181  iEvent.put(std::move(detIdCollection), interestingDetIdCollection_);
182 }
183 
187 
190 
edm::ESGetToken< CaloGeometry, CaloGeometryRecord > caloGeometryToken_
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
edm::EDGetTokenT< EcalRecHitCollection > recHitsToken_
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
bool checkFlags(const std::vector< int > &flagsvec) const
check if one of the flags in a set is true
Definition: EcalRecHit.h:190
void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override
producer
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
EcalSeverityLevel::SeverityLevel severityLevel(const DetId &id) const
Evaluate status from id use channelStatus from DB.
edm::ESGetToken< EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd > sevLvToken_
int iEvent
Definition: GenABIO.cc:224
T const * product() const
Definition: ESHandle.h:86
GlobalPoint getPosition(const DetId &id) const
Get the position of a given detector id.
Definition: CaloGeometry.cc:50
EgammaIsoDetIdCollectionProducer(const edm::ParameterSet &)
ctor
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:151
T perp() const
Magnitude of transverse component.
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
const DetId & detid() const
Definition: EcalRecHit.h:72
HLT enums.
edm::EDGetTokenT< T1Collection > emObjectToken_
float energy() const
Definition: EcalRecHit.h:68
bool checkFlag(int flag) const
check if the flag is true
Definition: EcalRecHit.h:187
def move(src, dest)
Definition: eostools.py:511