CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
hcalCalibUtils.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <fstream>
3 #include <algorithm>
4 #include <numeric>
5 
6 #include "TString.h"
8 
9 //#include "Calibration/HcalCalibAlgos/plugins/CommonUsefulStuff.h"
11 
12 using namespace std;
13 
14 void sumDepths(vector<TCell> &selectCells) {
15 
16  // Assignes teh sum of the energy in cells with the same iEta, iPhi to the cell with depth=1.
17  // All cells with depth>1 are removed form the container. If
18  // the cell at depth=1 is not present: create it and follow the procedure.
19 
20  if (selectCells.size()==0) return;
21 
22  vector<TCell> selectCellsDepth1;
23  vector<TCell> selectCellsHighDepth;
24 
25  //
26  // NB: Here we add depth 3 for iEta==16 in *HE* to the value in the barrel
27  // this approach is reflected in several of the following loops: make sure
28  // to check when making changes.
29  //
30  // In some documents it is described as having depth 1, the mapping in CMSSW uses depth 3.
31 
32  for (vector<TCell>::iterator i_it = selectCells.begin(); i_it != selectCells.end(); ++i_it) {
33  if (HcalDetId(i_it->id()).depth()==1) {
34  selectCellsDepth1.push_back(*i_it);
35  }
36  else {
37  selectCellsHighDepth.push_back(*i_it);
38  }
39  }
40 
41 
42 
43  // case where depth 1 has zero energy, but higher depths with same (iEta, iPhi) have energy.
44  // For iEta<15 there is one depth -> selectCellsHighDepth is empty and we do not get in the loop.
45  for (vector<TCell>::iterator i_it2 = selectCellsHighDepth.begin(); i_it2 != selectCellsHighDepth.end(); ++i_it2) {
46 
47 
48  // protect against corrupt data
49  if (HcalDetId(i_it2->id()).ietaAbs()<15 && HcalDetId(i_it2->id()).depth()>1) {
50  cout << "ERROR!!! there are no HB cells with depth>1 for iEta<15!\n"
51  << "Check the input data..." << endl;
52  cout << "HCalDetId: " << HcalDetId(i_it2->id()) << endl;
53  return;
54  }
55 
56 
57  bool foundDepthOne = false;
58  for (vector<TCell>::iterator i_it = selectCellsDepth1.begin(); i_it != selectCellsDepth1.end(); ++i_it) {
59  if (HcalDetId(i_it->id()).ieta()==HcalDetId(i_it2->id()).ieta() &&
60  HcalDetId(i_it->id()).iphi()==HcalDetId(i_it2->id()).iphi())
61  foundDepthOne = true;
62  continue;
63  }
64  if (!foundDepthOne) { // create entry for depth 1 with 0 energy
65 
66  UInt_t newId;
67  if (abs(HcalDetId(i_it2->id()).ieta())==16)
68  newId = HcalDetId(HcalBarrel, HcalDetId(i_it2->id()).ieta(), HcalDetId(i_it2->id()).iphi(), 1);
69  else
70  newId = HcalDetId(HcalDetId(i_it2->id()).subdet(), HcalDetId(i_it2->id()).ieta(), HcalDetId(i_it2->id()).iphi(), 1);
71 
72  selectCellsDepth1.push_back(TCell(newId, 0.0));
74  }
75  }
76 
77  for (vector<TCell>::iterator i_it = selectCellsDepth1.begin(); i_it != selectCellsDepth1.end(); ++i_it) {
78  for (vector<TCell>::iterator i_it2 = selectCellsHighDepth.begin(); i_it2 != selectCellsHighDepth.end(); ++i_it2) {
79  if (HcalDetId(i_it->id()).ieta()==HcalDetId(i_it2->id()).ieta() &&
80  HcalDetId(i_it->id()).iphi()==HcalDetId(i_it2->id()).iphi()) {
81  i_it->SetE(i_it->e()+i_it2->e());
82  i_it2->SetE(0.0); // paranoid, aren't we...
83  }
84  }
85  }
86 
87  // replace the original vectors with the new ones
88  selectCells = selectCellsDepth1;
89 
90  return;
91 }
92 
93 
94 void combinePhi(vector<TCell> &selectCells) {
95 
96  // Map: NxN -> N cluster
97  // Comine the targetE of cells with the same iEta
98 
99  if (selectCells.size()==0) return;
100 
101  // new container for the TCells
102  // dummy cell id created with iEta; iPhi=1; depth
103  // if combinePhi() is run after combining depths, depth=1
104  vector<TCell> combinedCells;
105 
106  map<UInt_t, vector<Float_t> > etaSliceE; // keyed by id of cell with iEta and **iPhi=1**
107 
108  // map the cells to the eta ring
109  vector<TCell>::iterator i_it = selectCells.begin();
110  for (; i_it != selectCells.end(); ++i_it) {
111  DetId id = HcalDetId(i_it->id());
112  UInt_t thisKey = HcalDetId(HcalDetId(id).subdet(), HcalDetId(id).ieta(), 1, HcalDetId(id).depth() );
113  etaSliceE[thisKey].push_back(i_it->e());
114  }
115 
116  map<UInt_t, vector<Float_t> >::iterator m_it = etaSliceE.begin();
117  for (; m_it != etaSliceE.end(); ++m_it) {
118  combinedCells.push_back(TCell(m_it->first, accumulate(m_it->second.begin(), m_it->second.end(), 0.0) ) );
119  }
120 
121  // replace the original TCell vector with the new one
122  selectCells = combinedCells;
123 
124 }
125 
126 
127 void combinePhi(vector<TCell> &selectCells, vector<TCell> &combinedCells) {
128 
129  // Map: NxN -> N cluster
130  // Comine the targetE of cells with the same iEta
131 
132  if (selectCells.size()==0) return;
133 
134  map<UInt_t, vector<Float_t> > etaSliceE; // keyed by id of cell with iEta and **iPhi=1**
135 
136  // map the cells to the eta ring
137  vector<TCell>::iterator i_it = selectCells.begin();
138  for (; i_it != selectCells.end(); ++i_it) {
139  DetId id = HcalDetId(i_it->id());
140  UInt_t thisKey = HcalDetId(HcalDetId(id).subdet(), HcalDetId(id).ieta(), 1, HcalDetId(id).depth() );
141  etaSliceE[thisKey].push_back(i_it->e());
142  }
143 
144  map<UInt_t, vector<Float_t> >::iterator m_it = etaSliceE.begin();
145  for (; m_it != etaSliceE.end(); ++m_it) {
146  combinedCells.push_back(TCell(m_it->first, accumulate(m_it->second.begin(), m_it->second.end(), 0.0) ) );
147  }
148 
149 }
150 
151 
152 
153 void getIEtaIPhiForHighestE(vector<TCell>& selectCells, Int_t& iEtaMostE, UInt_t& iPhiMostE) {
154 
155 
156  vector<TCell> summedDepthsCells = selectCells;
157 
158  sumDepths(summedDepthsCells);
159  vector<TCell>::iterator highCell = summedDepthsCells.begin();
160 
161  // sum depths locally to get highest energy tower
162 
163  Float_t highE = -999;
164 
165  for (vector<TCell>::iterator it=summedDepthsCells.begin(); it!=summedDepthsCells.end(); ++it) {
166  if (highE < it->e()) {
167  highCell = it;
168  highE = it->e();
169  }
170  }
171 
172  iEtaMostE = HcalDetId(highCell->id()).ieta();
173  iPhiMostE = HcalDetId(highCell->id()).iphi();
174 
175  return;
176 }
177 
178 
179 //
180 // Remove RecHits outside the 3x3 cluster and replace the vector that will
181 // be used in the minimization. Acts on "event" level.
182 // This can not be done for iEta>20 due to segmentation => in principle the result should be restricted
183 // to iEta<20. Attempted to minimize affect at the boundary without a sharp jump.
184 
185 void filterCells3x3(vector<TCell>& selectCells, Int_t iEtaMaxE, UInt_t iPhiMaxE) {
186 
187  vector<TCell> filteredCells;
188 
189  Int_t dEta, dPhi;
190 
191  for (vector<TCell>::iterator it=selectCells.begin(); it!=selectCells.end(); ++it) {
192 
193  Bool_t passDEta = false;
194  Bool_t passDPhi = false;
195 
196  dEta = HcalDetId(it->id()).ieta() - iEtaMaxE;
197  dPhi = HcalDetId(it->id()).iphi() - iPhiMaxE;
198 
199  if (dPhi > 36) dPhi -= 72;
200  if (dPhi < -36) dPhi += 72;
201 
202  if (abs(dEta)<=1 || (iEtaMaxE * HcalDetId(it->id()).ieta() == -1)) passDEta = true;
203 
204  if (abs(iEtaMaxE)<=20) {
205 
206  if (abs(HcalDetId(it->id()).ieta())<=20) {
207  if (abs(dPhi)<=1) passDPhi = true;
208  }
209  else {
210  // iPhi is labelled by odd numbers
211  if (iPhiMaxE%2==0){
212  if (abs(dPhi)<=1) passDPhi = true;
213  }
214  else {
215  if (dPhi== -2 || dPhi==0) passDPhi = true;
216  }
217  }
218 
219  } // if hottest cell with iEta<=20
220 
221  else {
222  if (abs(HcalDetId(it->id()).ieta())<=20) {
223  if (abs(dPhi)<=1 || dPhi==2) passDPhi = true;
224  }
225  else {
226  if (abs(dPhi)<=2) passDPhi = true;
227  }
228  } // if hottest cell with iEta>20
229 
230  if (passDEta && passDPhi) filteredCells.push_back(*it);
231  }
232 
233  selectCells = filteredCells;
234 
235  return;
236 }
237 
238 //
239 // Remove RecHits outside the 5x5 cluster and replace the vector that will
240 // be used in the minimization. Acts on "event" level
241 // In principle the ntuple should be produced with 5x5 already precelected
242 //
243 // Size for iEta>20 is 3x3, but the segmentation changes by x2 in phi.
244 // There is some bias in the selection of towers near the boundary
245 
246 void filterCells5x5(vector<TCell>& selectCells, Int_t iEtaMaxE, UInt_t iPhiMaxE) {
247 
248  vector<TCell> filteredCells;
249 
250  Int_t dEta, dPhi;
251 
252  for (vector<TCell>::iterator it=selectCells.begin(); it!=selectCells.end(); ++it) {
253 
254  dEta = HcalDetId(it->id()).ieta() - iEtaMaxE;
255  dPhi = HcalDetId(it->id()).iphi() - iPhiMaxE;
256 
257  if (dPhi > 36) dPhi -= 72;
258  if (dPhi < -36) dPhi += 72;
259 
260  bool passDPhi = (abs(dPhi)<3);
261 
262  bool passDEta = (abs(dEta)<3 || (iEtaMaxE * HcalDetId(it->id()).ieta() == -2) );
263  // includes +/- eta boundary
264 
265  if (passDPhi && passDEta) filteredCells.push_back(*it);
266 
267  }
268 
269  selectCells = filteredCells;
270 
271  return;
272 }
273 
274 
275 
276 
277 // this is for the problematic layer near the HB/HE boundary
278 // sum depths 1,2 in towers 15,16
279 
280 void sumSmallDepths(vector<TCell> &selectCells) {
281 
282  if (selectCells.size()==0) return;
283 
284  vector<TCell> newCells; // holds unaffected cells to which the modified ones are added
285  vector<TCell> manipulatedCells; // the ones that are combined
286 
287  for (vector<TCell>::iterator i_it = selectCells.begin(); i_it != selectCells.end(); ++i_it) {
288 
289  if ( (HcalDetId(i_it->id()).ietaAbs()==15 && HcalDetId(i_it->id()).depth()<=2) ||
290  (HcalDetId(i_it->id()).ietaAbs()==16 && HcalDetId(i_it->id()).depth()<=2)) {
291  manipulatedCells.push_back(*i_it);
292  }
293  else {
294  newCells.push_back(*i_it);
295  }
296 
297  }
298 
299  // if the list is empty there is nothing to manipulate
300  // leave the original vector unchanged
301 
302  if (manipulatedCells.size()<1) {
303  newCells.clear();
304  return;
305  }
306 
307 
308  // See what cells are needed to hold the combined information:
309  // Make holders for depth=1 for each (iEta,iPhi)
310  // if a cell with these values is present in "manupulatedCells"
311  vector<UInt_t> dummyIds; // to keep track of kreated cells
312  vector<TCell> createdCells; // cells that need to be added or they exists;
313 
314  for (vector<TCell>::iterator i_it = manipulatedCells.begin(); i_it!=manipulatedCells.end(); ++i_it) {
315  UInt_t dummyId = HcalDetId(HcalDetId(i_it->id()).subdet(), HcalDetId(i_it->id()).ieta(), HcalDetId(i_it->id()).iphi(), 1);
316  if (find(dummyIds.begin(), dummyIds.end(), dummyId)==dummyIds.end()) {
317  dummyIds.push_back(dummyId);
318  createdCells.push_back(TCell(dummyId, 0.0));
319  }
320  }
321 
322  for (vector<TCell>::iterator i_it = createdCells.begin(); i_it!=createdCells.end(); ++i_it) {
323  for (vector<TCell>::iterator i_it2 = manipulatedCells.begin(); i_it2!=manipulatedCells.end(); ++i_it2) {
324  if (HcalDetId(i_it->id()).ieta()==HcalDetId(i_it2->id()).ieta() &&
325  HcalDetId(i_it->id()).iphi()==HcalDetId(i_it2->id()).iphi() &&
326  HcalDetId(i_it2->id()).depth()<=2) {
327  i_it->SetE(i_it->e()+i_it2->e());
328  }
329  }
330  }
331 
332  for (vector<TCell>::iterator i_it = createdCells.begin(); i_it!=createdCells.end(); ++i_it) {
333  newCells.push_back(*i_it);
334  }
335 
336 
337  // replace the original vectors with the new ones
338  selectCells = newCells;
339 
340  return;
341 }
342 
343 
344 void filterCellsInCone(std::vector<TCell>& selectCells, const GlobalPoint hitPositionHcal,
345  Float_t maxConeDist, const CaloGeometry* theCaloGeometry) {
346 
347  vector<TCell> filteredCells;
348 
349  for (vector<TCell>::iterator it=selectCells.begin(); it!=selectCells.end(); ++it) {
350 
351  const GlobalPoint recHitPoint = theCaloGeometry->getPosition(it->id());
352 
353  if (getDistInPlaneSimple(hitPositionHcal, recHitPoint)<= maxConeDist)
354  filteredCells.push_back(*it);
355  }
356 
357  selectCells = filteredCells;
358 
359  return;
360 }
361 
362 
363 // From Jim H. => keep till the code is included centrally
364 /*
365 double getDistInPlaneSimple(const GlobalPoint caloPoint, const GlobalPoint rechitPoint) {
366 
367  // Simplified version of getDistInPlane
368  // Assume track direction is origin -> point of hcal intersection
369 
370  const GlobalVector caloIntersectVector(caloPoint.x(),
371  caloPoint.y(),
372  caloPoint.z());
373 
374  const GlobalVector caloIntersectUnitVector = caloIntersectVector.unit();
375 
376  const GlobalVector rechitVector(rechitPoint.x(),
377  rechitPoint.y(),
378  rechitPoint.z());
379 
380  const GlobalVector rechitUnitVector = rechitVector.unit();
381 
382  double dotprod = caloIntersectUnitVector.dot(rechitUnitVector);
383  double rechitdist = caloIntersectVector.mag()/dotprod;
384 
385 
386  const GlobalVector effectiveRechitVector = rechitdist*rechitUnitVector;
387  const GlobalPoint effectiveRechitPoint(effectiveRechitVector.x(),
388  effectiveRechitVector.y(),
389  effectiveRechitVector.z());
390 
391 
392  GlobalVector distance_vector = effectiveRechitPoint-caloPoint;
393 
394  if (dotprod > 0.)
395  {
396  return distance_vector.mag();
397  }
398  else
399  {
400  return 999999.;
401 
402  }
403 
404 
405 }
406 */
void filterCellsInCone(std::vector< TCell > &selectCells, const GlobalPoint hitPositionHcal, Float_t maxConeDist, const CaloGeometry *theCaloGeometry)
tuple sumSmallDepths
Definition: diJetCalib.py:58
Definition: TCell.h:15
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
tuple maxConeDist
Definition: diJetCalib.py:71
tuple combinePhi
Definition: diJetCalib.py:54
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
const GlobalPoint & getPosition(const DetId &id) const
Get the position of a given detector id.
Definition: CaloGeometry.cc:68
void filterCells3x3(std::vector< TCell > &selectCells, Int_t iEta, UInt_t iPhi)
tuple sumDepths
Definition: diJetCalib.py:51
Definition: DetId.h:18
void filterCells5x5(std::vector< TCell > &selectCells, Int_t iEta, UInt_t iPhi)
double getDistInPlaneSimple(const GlobalPoint caloPoint, const GlobalPoint rechitPoint)
Definition: ConeDefinition.h:9
tuple cout
Definition: gather_cfg.py:121
void getIEtaIPhiForHighestE(std::vector< TCell > &selectCells, Int_t &iEta, UInt_t &iPhi)