CMS 3D CMS Logo

EgammaHLTNxNClusterProducer.cc
Go to the documentation of this file.
1 // C/C++ headers
2 //#include <iostream>
3 //#include <vector>
4 
5 // Reconstruction Classes
12 
21 
23 #include "TVector3.h"
24 
25 #include <memory>
26 
28  doBarrel_ (ps.getParameter<bool>("doBarrel")),
29  doEndcaps_ (ps.getParameter<bool>("doEndcaps")),
30  barrelHitProducer_ (consumes<EcalRecHitCollection>(ps.getParameter< edm::InputTag > ("barrelHitProducer"))),
31  endcapHitProducer_ (consumes<EcalRecHitCollection>(ps.getParameter< edm::InputTag > ("endcapHitProducer"))),
32  clusEtaSize_ (ps.getParameter<int> ("clusEtaSize")),
33  clusPhiSize_ (ps.getParameter<int> ("clusPhiSize")),
34  barrelClusterCollection_(ps.getParameter<std::string>("barrelClusterCollection")),
35  endcapClusterCollection_(ps.getParameter<std::string>("endcapClusterCollection")),
36  clusSeedThr_ (ps.getParameter<double> ("clusSeedThr")),
37  clusSeedThrEndCap_ (ps.getParameter<double> ("clusSeedThrEndCap")),
38  useRecoFlag_ (ps.getParameter<bool>("useRecoFlag")),
39  flagLevelRecHitsToUse_ (ps.getParameter<int>("flagLevelRecHitsToUse")),
40  useDBStatus_ (ps.getParameter<bool>("useDBStatus")),
41  statusLevelRecHitsToUse_(ps.getParameter<int>("statusLevelRecHitsToUse")),
42  maxNumberofSeeds_ (ps.getParameter<int> ("maxNumberofSeeds")),
43  maxNumberofClusters_ (ps.getParameter<int> ("maxNumberofClusters")),
44  debug_ (ps.getParameter<int> ("debugLevel")),
45  posCalculator_ (PositionCalc( ps.getParameter<edm::ParameterSet>("posCalcParameters"))) {
46 
47  produces< reco::BasicClusterCollection >(barrelClusterCollection_);
48  produces< reco::BasicClusterCollection >(endcapClusterCollection_);
49 }
50 
51 
53 {}
54 
56 
58  desc.add<bool>(("doBarrel"), true);
59  desc.add<bool>(("doEndcaps"), true);
60  desc.add<edm::InputTag>(("barrelHitProducer"), edm::InputTag("hltEcalRegionalPi0EtaRecHit","EcalRecHitsEB"));
61  desc.add<edm::InputTag>(("endcapHitProducer"), edm::InputTag("hltEcalRegionalPi0EtaRecHit","EcalRecHitsEE"));
62  desc.add<int>(("clusEtaSize"), 3);
63  desc.add<int>(("clusPhiSize"), 3);
64  desc.add<std::string>(("barrelClusterCollection"), "Simple3x3ClustersBarrel");
65  desc.add<std::string>(("endcapClusterCollection"), "Simple3x3ClustersEndcap");
66  desc.add<double>(("clusSeedThr"), 0.5);
67  desc.add<double>(("clusSeedThrEndCap"), 1.0);
68  desc.add<bool>(("useRecoFlag"), false);
69  desc.add<int>(("flagLevelRecHitsToUse"), 1);
70  desc.add<bool>(("useDBStatus"), true);
71  desc.add<int>(("statusLevelRecHitsToUse"), 1);
72 
73  edm::ParameterSetDescription posCalcPSET;
74  posCalcPSET.add<double>("T0_barl", 7.4);
75  posCalcPSET.add<double>("T0_endc", 3.1);
76  posCalcPSET.add<double>("T0_endcPresh", 1.2);
77  posCalcPSET.add<double>("W0", 4.2);
78  posCalcPSET.add<double>("X0", 0.89);
79  posCalcPSET.add<bool>("LogWeighted", true);
80  desc.add<edm::ParameterSetDescription>("posCalcParameters", posCalcPSET);
81 
82  desc.add<int>(("maxNumberofSeeds"), 1000);
83  desc.add<int>(("maxNumberofClusters"), 200);
84  desc.add<int>(("debugLevel"), 0);
85  descriptions.add(("hltEgammaHLTNxNClusterProducer"), desc);
86 }
87 
89 
90  if(doBarrel_) {
91  edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
92  evt.getByToken(barrelHitProducer_, barrelRecHitsHandle);
93  if (!barrelRecHitsHandle.isValid()) {
94  LogDebug("") << "EgammaHLTNxNClusterProducer Error! can't get product eb hit!" << std::endl;
95  }
96 
97  const EcalRecHitCollection *hits_eb = barrelRecHitsHandle.product();
98  if( debug_>=2 )
99  LogDebug("")<<"EgammaHLTNxNClusterProducer nEBrechits: "<< evt.id().run()<<" event "<<evt.id().event() <<" "<< hits_eb->size()<<std::endl;
100 
102 
103  }
104 
105 
106  if(doEndcaps_) {
107  edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
108  evt.getByToken(endcapHitProducer_, endcapRecHitsHandle);
109  if (!endcapRecHitsHandle.isValid()) {
110  LogDebug("") << "EgammaHLTNxNClusterProducer Error! can't get product ee hit!" << std::endl;
111  }
112 
113  const EcalRecHitCollection *hits_ee = endcapRecHitsHandle.product();
114  if( debug_>=2 )
115  LogDebug("")<<"EgammaHLTNxNClusterProducer nEErechits: "<< evt.id().run()<<" event "<<evt.id().event() <<" "<< hits_ee->size()<<std::endl;
117  }
118 }
119 
120 
121 
123 
124  if(useRecoFlag_ ){
125  int flag = rh.recoFlag();
126  if( flagLevelRecHitsToUse_ ==0){
127  if( flag != 0) return false;
128  }
129  else if( flagLevelRecHitsToUse_ ==1){
130  if( flag !=0 && flag != 4 ) return false;
131  }
132  else if( flagLevelRecHitsToUse_ ==2){
133  if( flag !=0 && flag != 4 && flag != 6 && flag != 7) return false;
134  }
135  }
136  if ( useDBStatus_){
137  int status = int(channelStatus[rh.id().rawId()].getStatusCode());
138  if ( status > statusLevelRecHitsToUse_ ) return false;
139  }
140 
141  return true;
142 }
143 
144 
145 
146 
149 {
150 
151 
154  if ( useDBStatus_ ) es.get<EcalChannelStatusRcd>().get(csHandle);
155  const EcalChannelStatus &channelStatus = *csHandle;
156 
157 
158  std::vector<EcalRecHit> seeds;
159 
160  double clusterSeedThreshold ;
161  if (detector == reco::CaloID::DET_ECAL_BARREL){
162  clusterSeedThreshold = clusSeedThr_;
163  }else{
164  clusterSeedThreshold = clusSeedThrEndCap_;
165  }
166 
167 
168  for(EcalRecHitCollection::const_iterator itt = hits->begin(); itt != hits->end(); itt++){
169  double energy = itt->energy();
170  if( ! checkStatusOfEcalRecHit(channelStatus, *itt) ) continue;
171  if (energy > clusterSeedThreshold ) seeds.push_back(*itt);
172 
173  if( int(seeds.size()) > maxNumberofSeeds_){ //too many seeds, like beam splash events
174  seeds.clear();
175  break;
176  }
177  }
178 
179  // get the geometry and topology from the event setup:
180  edm::ESHandle<CaloGeometry> geoHandle;
181  es.get<CaloGeometryRecord>().get(geoHandle);
182 
183  const CaloSubdetectorGeometry *geometry_p;
184  std::unique_ptr<CaloSubdetectorTopology> topology_p;
185  if (detector == reco::CaloID::DET_ECAL_BARREL) {
186  geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
187  topology_p.reset(new EcalBarrelTopology(geoHandle));
188  }else {
189  geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
190  topology_p.reset(new EcalEndcapTopology(geoHandle));
191  }
192 
193  const CaloSubdetectorGeometry *geometryES_p;
194  geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
195 
196 
197 
198  std::vector<reco::BasicCluster> clusters;
199  std::vector<DetId> usedXtals;
200 
201  // sort seed according to Energy
202  sort(seeds.begin(), seeds.end(), ecalRecHitSort());
203 
204 
205 
206  for (std::vector<EcalRecHit>::iterator itseed=seeds.begin(); itseed!=seeds.end(); itseed++) {
207  DetId seed_id = itseed->id();
208  std::vector<DetId>::const_iterator usedIds;
209 
210  std::vector<DetId>::iterator itdet = find(usedXtals.begin(),usedXtals.end(),seed_id);
211  if(itdet != usedXtals.end()) continue;
212 
213  std::vector<DetId> clus_v = topology_p->getWindow(seed_id,clusEtaSize_,clusPhiSize_);
214  std::vector<std::pair<DetId, float> > clus_used;
215 
216  float clus_energy = 0;
217 
218  for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
219  DetId detid = *det;
220 
221  //not yet used
222  std::vector<DetId>::iterator itdet = find(usedXtals.begin(),usedXtals.end(),detid);
223  if(itdet != usedXtals.end()) continue;
224  //inside the collection
226  if( hit == hits->end()) continue;
227 
228  if( ! checkStatusOfEcalRecHit(channelStatus, *hit) ) continue;
229 
230  usedXtals.push_back(detid);
231  clus_used.push_back(std::pair<DetId, float>(detid, 1.) );
232  clus_energy += hit->energy();
233 
234  }
235 
236  if( clus_energy <= 0 ) continue;
237 
238  math::XYZPoint clus_pos = posCalculator_.Calculate_Location(clus_used,hits,geometry_p,geometryES_p);
239 
240  if (debug_>=2 ) LogDebug("")<<"nxn_cluster in run "<< evt.id().run()<<" event "<<evt.id().event()<<" energy: "<<clus_energy <<" eta: "<< clus_pos.Eta()<<" phi: "<< clus_pos.Phi()<<" nRecHits: "<< clus_used.size() <<std::endl;
241 
242  clusters.push_back(reco::BasicCluster(clus_energy, clus_pos, reco::CaloID(detector), clus_used, reco::CaloCluster::island, seed_id));
243  if( int(clusters.size()) > maxNumberofClusters_){
244  clusters.clear();
245  break;
246  }
247 
248  }
249 
250 
251  //Create empty output collections
252  auto clusters_p = std::make_unique<reco::BasicClusterCollection>();
253  clusters_p->assign(clusters.begin(), clusters.end());
254  if (detector == reco::CaloID::DET_ECAL_BARREL){
255  if(debug_>=1) LogDebug("")<<"nxnclusterProducer: "<<clusters_p->size() <<" made in barrel"<<std::endl;
256  evt.put(std::move(clusters_p), barrelClusterCollection_);
257  }
258  else {
259  if(debug_>=1) LogDebug("")<<"nxnclusterProducer: "<<clusters_p->size() <<" made in endcap"<<std::endl;
260  evt.put(std::move(clusters_p), endcapClusterCollection_);
261  }
262 
263 }
264 
265 
#define LogDebug(id)
RunNumber_t run() const
Definition: EventID.h:39
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:44
EventNumber_t event() const
Definition: EventID.h:41
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:136
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:519
std::vector< EcalRecHit >::const_iterator const_iterator
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:20
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
EgammaHLTNxNClusterProducer(const edm::ParameterSet &ps)
const edm::EDGetTokenT< EcalRecHitCollection > endcapHitProducer_
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
ParameterDescriptionBase * add(U const &iLabel, T const &value)
bool isValid() const
Definition: HandleBase.h:74
const_iterator end() const
Definition: DetId.h:18
DetId id() const
get the id
Definition: EcalRecHit.h:77
T const * product() const
Definition: Handle.h:81
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
Flags recoFlag() const
DEPRECATED provided for temporary backward compatibility.
Definition: EcalRecHit.h:206
const T & get() const
Definition: EventSetup.h:58
void add(std::string const &label, ParameterSetDescription const &psetDescription)
void produce(edm::Event &, const edm::EventSetup &) override
edm::EventID id() const
Definition: EventBase.h:60
iterator find(key_type k)
void makeNxNClusters(edm::Event &evt, const edm::EventSetup &es, const EcalRecHitCollection *hits, const reco::CaloID::Detectors detector)
HLT enums.
size_type size() const
math::XYZPoint Calculate_Location(const HitsAndFractions &iDetIds, const edm::SortedCollection< HitType > *iRecHits, const CaloSubdetectorGeometry *iSubGeom, const CaloSubdetectorGeometry *iESGeom=0)
Definition: PositionCalc.h:68
bool checkStatusOfEcalRecHit(const EcalChannelStatus &channelStatus, const EcalRecHit &rh)
const edm::EDGetTokenT< EcalRecHitCollection > barrelHitProducer_
def move(src, dest)
Definition: eostools.py:510
const_iterator begin() const