CMS 3D CMS Logo

CaloTowersMerger.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: CaloTowersMerger
4 // Class: CaloTowersMerger
5 //
13 //
14 // Original Author: Jean-Roch Vlimant,40 3-A28,+41227671209,
15 // Created: Thu Nov 4 16:36:30 CET 2010
16 //
17 //
18 
19 // Anton Anastassov (Northwestern):
20 // Add code for actual tower merging, define two distinct inputs of
21 // "regular" and "extra" towers
22 // This functionality is subject to some restrictions
23 // (see notes below)
24 
31 
32 #include <memory>
33 
34 //
35 // class declaration
36 //
37 
39 public:
40  explicit CaloTowersMerger(const edm::ParameterSet&);
41  ~CaloTowersMerger() override;
42 
43  CaloTower mergedTower(const CaloTower& t1, const CaloTower& t2);
44 
45 private:
46  void produce(edm::Event&, const edm::EventSetup&) override;
47 
48  // ----------member data ---------------------------
49 
53 };
54 
55 //
56 // constants, enums and typedefs
57 //
58 
59 //
60 // static data member definitions
61 //
62 
63 //
64 // constructors and destructor
65 //
67  regularTowerTag = iConfig.getParameter<edm::InputTag>("regularTowerTag");
68  extraTowerTag = iConfig.getParameter<edm::InputTag>("extraTowerTag");
69 
70  // register for data access
71  tok_reg_ = consumes<CaloTowerCollection>(regularTowerTag);
72  tok_ext_ = consumes<CaloTowerCollection>(extraTowerTag);
73 
74  //register your products
75  produces<CaloTowerCollection>();
76 }
77 
79  // do anything here that needs to be done at desctruction time
80  // (e.g. close files, deallocate resources etc.)
81 }
82 
83 //
84 // member functions
85 //
86 
87 // ------------ method called to produce the data ------------
89  edm::Handle<CaloTowerCollection> regTower, extraTower;
90 
91  iEvent.getByToken(tok_reg_, regTower);
92  iEvent.getByToken(tok_ext_, extraTower);
93 
94  if (!regTower.isValid() && !extraTower.isValid()) {
95  edm::LogError("CaloTowersMerger") << "both input tag:" << regularTowerTag << " and " << extraTowerTag
96  << " are invalid. empty merged collection";
97  iEvent.put(std::make_unique<CaloTowerCollection>());
98  return;
99  } else if (!regTower.isValid() || !extraTower.isValid()) {
100  if (!regTower.isValid() && extraTower.isValid())
101  regTower = extraTower;
102  iEvent.put(std::make_unique<CaloTowerCollection>(*regTower));
103  return;
104  } else {
105  //both valid input collections: merging
106  auto output = std::make_unique<CaloTowerCollection>();
107  output->reserve(regTower->size() + extraTower->size());
108 
109  CaloTowerCollection::const_iterator rt_begin = regTower->begin();
110  CaloTowerCollection::const_iterator rt_end = regTower->end();
111  CaloTowerCollection::const_iterator rt_it = rt_begin;
112 
113  //vector of overlapping towers
114  std::vector<CaloTowerCollection::const_iterator> overlappingTowers;
115  overlappingTowers.reserve(extraTower->size());
116 
117  for (; rt_it != rt_end; ++rt_it) {
118  CaloTowerCollection::const_iterator et_it = extraTower->find(rt_it->id());
119  if (et_it != extraTower->end()) {
120  //need to merge the components
121  //FIXME
122 
124  //one needs to merge t1 and t2 into mergedTower
125  //end FIXME
126  CaloTower mt = mergedTower(*rt_it, *et_it);
127 
128  output->push_back(mt);
129  overlappingTowers.push_back(et_it);
130 
131  } else {
132  //just copy the regular tower over
133  output->push_back(*rt_it);
134  }
135  }
136  CaloTowerCollection::const_iterator et_begin = extraTower->begin();
137  CaloTowerCollection::const_iterator et_end = extraTower->end();
138  CaloTowerCollection::const_iterator et_it = et_begin;
139  for (; et_it != et_end; ++et_it) {
140  if (std::find(overlappingTowers.begin(), overlappingTowers.end(), et_it) == overlappingTowers.end())
141  //non overlapping tower
142  //copy the extra tower over
143  output->push_back(*et_it);
144  }
145  iEvent.put(std::move(output));
146  }
147 }
148 
149 // Make a new tower by merging two towers containing exclusive hits
150 // This cannot be done fully consistently and independent of the
151 // creation mechnaism...
152 // This functionlaity it to be used only for testing the effects
153 // of rejected bad hits.
154 
156  double newOuterE = 0;
157 
158  // HO energies are always saved (even if useHO=false)
159  // make sure there is no double counting
160  // crude test if one has HO energy and the other not (possibly due to bad hit cleanup)
161  // However: for |iEta|>15 E_outer has different meening:
162  // it holds either the energy in the outer depths in HCAL, etc
163 
164  if (rt.ietaAbs() < 16 && (fabs(rt.outerEnergy() - et.outerEnergy()) < 0.00001)) {
165  // there is no differnce in the store HO energies
166  newOuterE = rt.outerEnergy();
167  } else {
168  newOuterE = rt.outerEnergy() + et.outerEnergy();
169  }
170 
171  bool rt_hasEcalConstit = false;
172  bool et_hasEcalConstit = false;
173 
174  bool rt_hasHcalConstit = false;
175  bool et_hasHcalConstit = false;
176 
177  // check if there are HCAL/ECAL constituents in the towers
178 
179  std::vector<DetId>::const_iterator rc_begin = rt.constituents().begin();
180  std::vector<DetId>::const_iterator rc_end = rt.constituents().end();
181  std::vector<DetId>::const_iterator rc_it;
182 
183  for (rc_it = rc_begin; rc_it != rc_end; ++rc_it) {
184  if (rc_it->det() == DetId::Hcal)
185  rt_hasHcalConstit = true;
186  break;
187  }
188  for (rc_it = rc_begin; rc_it != rc_end; ++rc_it) {
189  if (rc_it->det() == DetId::Ecal)
190  rt_hasEcalConstit = true;
191  break;
192  }
193 
194  std::vector<DetId>::const_iterator ec_begin = et.constituents().begin();
195  std::vector<DetId>::const_iterator ec_end = et.constituents().end();
196  std::vector<DetId>::const_iterator ec_it;
197 
198  for (ec_it = ec_begin; ec_it != ec_end; ++ec_it) {
199  if (ec_it->det() == DetId::Hcal)
200  et_hasHcalConstit = true;
201  break;
202  }
203  for (ec_it = ec_begin; ec_it != ec_end; ++ec_it) {
204  if (ec_it->det() == DetId::Ecal)
205  et_hasEcalConstit = true;
206  break;
207  }
208 
209  std::vector<DetId> combinedConstituents = rt.constituents();
210  for (ec_it = ec_begin; ec_it != ec_end; ++ec_it) {
211  // if properly resconstructed, the only possible overlap should be for HO hits that
212  // are always listed as constituents if above thereshold (old JetMET request)
213  if (std::find(combinedConstituents.begin(), combinedConstituents.end(), *ec_it) == combinedConstituents.end())
214  combinedConstituents.push_back(*ec_it);
215  }
216 
217  GlobalPoint newEmPosition(0.0, 0.0, 0.0);
218 
219  // The following assumes the current default
220  // momentum reconstruction method (1) and
221  // accepted rechits with some threshod >0
222 
223  if (rt_hasEcalConstit && et_hasEcalConstit) {
224  if (rt.emEnergy() > 0 && et.emEnergy() > 0) {
225  double sumEmE = rt.emEnergy() + et.emEnergy();
226 
227  double x = rt.emEnergy() * rt.emPosition().x() + et.emEnergy() * et.emPosition().x();
228  double y = rt.emEnergy() * rt.emPosition().y() + et.emEnergy() * et.emPosition().y();
229  double z = rt.emEnergy() * rt.emPosition().z() + et.emEnergy() * et.emPosition().z();
230 
231  GlobalPoint weightedEmdPosition(x / sumEmE, y / sumEmE, z / sumEmE);
232  newEmPosition = weightedEmdPosition;
233  }
234 
235  } else if (rt_hasEcalConstit && !et_hasEcalConstit) {
236  newEmPosition = rt.emPosition();
237  } else if (!rt_hasEcalConstit && et_hasEcalConstit) {
238  newEmPosition = et.emPosition();
239  }
240 
241  GlobalPoint newHadPosition(0.0, 0.0, 0.0);
242  // had positions are the same if there is at least one constituent
243  if (rt_hasHcalConstit) {
244  newHadPosition = rt.hadPosition();
245  } else if (et_hasHcalConstit) {
246  newHadPosition = et.hadPosition();
247  }
248 
249  // MAke the new tower and set all values
250 
251  CaloTower mergedTower(rt.id(),
252  rt.emEnergy() + et.emEnergy(),
253  rt.hadEnergy() + et.hadEnergy(),
254  newOuterE,
255  -1,
256  -1,
257  rt.p4() + et.p4(),
258  newEmPosition,
259  newHadPosition);
260 
261  mergedTower.addConstituents(combinedConstituents);
262 
263  (rt.hottestCellE() > et.hottestCellE()) ? mergedTower.setHottestCellE(rt.hottestCellE())
264  : mergedTower.setHottestCellE(et.hottestCellE());
265 
266  unsigned int numBadHcalChan = rt.numBadHcalCells() - et.numProblematicHcalCells() - rt.numRecoveredHcalCells();
267  unsigned int numBadEcalChan = rt.numBadEcalCells() - et.numProblematicEcalCells() - rt.numRecoveredEcalCells();
268 
269  unsigned int numProbHcalChan = rt.numProblematicHcalCells() + et.numProblematicHcalCells();
270  unsigned int numProbEcalChan = rt.numProblematicEcalCells() + et.numProblematicEcalCells();
271 
272  unsigned int numRecHcalChan = rt.numRecoveredHcalCells() + et.numRecoveredHcalCells();
273  unsigned int numRecEcalChan = rt.numRecoveredEcalCells() + et.numRecoveredEcalCells();
274 
276  numBadHcalChan, numBadEcalChan, numRecHcalChan, numRecEcalChan, numProbHcalChan, numProbEcalChan);
277 
278  // use timing from the good tower only (for now, in default reco we use information from good hits only)
279  // time is saved as integer but returned as float in (ns)
280  mergedTower.setEcalTime(int(rt.ecalTime() * 100.0 + 0.5));
281  mergedTower.setHcalTime(int(rt.hcalTime() * 100.0 + 0.5));
282 
283  return mergedTower;
284 }
285 
286 //define this as a plug-in
edm::EDGetTokenT< CaloTowerCollection > tok_reg_
CaloTower mergedTower(const CaloTower &t1, const CaloTower &t2)
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
CaloTowersMerger(const edm::ParameterSet &)
size_type size() const
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
void setCaloTowerStatus(unsigned int numBadHcalChan, unsigned int numBadEcalChan, unsigned int numRecHcalChan, unsigned int numRecEcalChan, unsigned int numProbHcalChan, unsigned int numProbEcalChan)
Definition: CaloTower.cc:221
edm::EDGetTokenT< CaloTowerCollection > tok_ext_
std::vector< CaloTower >::const_iterator const_iterator
edm::InputTag extraTowerTag
Log< level::Error, false > LogError
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
void setHottestCellE(double e)
Definition: CaloTower.h:120
void addConstituents(const std::vector< DetId > &ids)
Definition: CaloTower.cc:199
int iEvent
Definition: GenABIO.cc:224
const_iterator begin() const
const_iterator end() const
void produce(edm::Event &, const edm::EventSetup &) override
edm::InputTag regularTowerTag
~CaloTowersMerger() override
bool isValid() const
Definition: HandleBase.h:70
void setEcalTime(int t)
Definition: CaloTower.h:83
iterator find(key_type k)
void setHcalTime(int t)
Definition: CaloTower.h:84
def move(src, dest)
Definition: eostools.py:511