CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
HBHETimingShapedFlag.cc
Go to the documentation of this file.
1 #include <iostream>
4 
5 /*
6 v1.0 by Jeff Temple
7 29 August 2009
8 
9 This takes the timing envelope algorithms developed by
10 Phil Dudero and uses them to apply a RecHit Flag.
11 */
12 
14 
16 
17 HBHETimingShapedFlagSetter::HBHETimingShapedFlagSetter(const std::vector<double>& v_userEnvelope) {
18  makeTfilterEnvelope(v_userEnvelope);
19 
20  ignorelowest_ = false;
21  ignorehighest_ = false;
22  win_offset_ = 0.;
23  win_gain_ = 1.;
24 }
25 
26 HBHETimingShapedFlagSetter::HBHETimingShapedFlagSetter(const std::vector<double>& v_userEnvelope,
27  bool ignorelowest,
28  bool ignorehighest,
29  double win_offset,
30  double win_gain) {
31  makeTfilterEnvelope(v_userEnvelope);
32 
33  ignorelowest_ = ignorelowest; // can ignore flagging hits below lowest energy threshold
34  ignorehighest_ = ignorehighest; // can ignore flagging hits above highest energy threshold
35  win_offset_ = win_offset; // timing offset
36  win_gain_ = win_gain; // time gain
37 }
38 
39 void HBHETimingShapedFlagSetter::makeTfilterEnvelope(const std::vector<double>& v_userEnvelope) {
40  // Transform vector of doubles into a map of <energy,lowtime,hitime> triplets
41  // Add extra protection in case vector of doubles is not a multiple of 3
42  if (v_userEnvelope.size() % 3)
43  throw cms::Exception("Invalid tfilterEnvelope definition") << "Must be one energy and two times per point";
44 
45  for (unsigned int i = 0; i < v_userEnvelope.size(); i += 3) {
46  int intGeV = (int)(v_userEnvelope[i] + 0.5);
47  std::pair<double, double> pairOfTimes = std::make_pair(v_userEnvelope[i + 1], v_userEnvelope[i + 2]);
48  if ((pairOfTimes.first > 0) || (pairOfTimes.second < 0))
49  throw cms::Exception("Invalid tfilterEnvelope definition")
50  << "Min and max time values must straddle t=0; use win_offset to shift";
51 
52  tfilterEnvelope_[intGeV] = pairOfTimes;
53  }
54 }
55 
57  std::cout << "Timing Energy/Time parameters are:" << std::endl;
58  TfilterEnvelope_t::const_iterator it;
59  for (it = tfilterEnvelope_.begin(); it != tfilterEnvelope_.end(); ++it)
60  std::cout << "\t" << it->first << " GeV\t" << it->second.first << " ns\t" << it->second.second << " ns"
61  << std::endl;
62 
63  std::cout << "ignorelowest = " << ignorelowest_ << std::endl;
64  std::cout << "ignorehighest = " << ignorehighest_ << std::endl;
65  std::cout << "win_offset = " << win_offset_ << std::endl;
66  std::cout << "win_gain = " << win_gain_ << std::endl;
67 }
68 
70  int status = 0; // 3 bits reserved;status can range from 0-7
71 
72  // tfilterEnvelope stores triples of energy and high/low time limits;
73  // make sure we're checking over an even number of values
74  // energies are guaranteed by std::map to appear in increasing order
75 
76  // need at least two values to make comparison, and must
77  // always have energy, time pair; otherwise, assume "in time" and don't set bits
78  if (tfilterEnvelope_.empty())
79  return 0;
80 
81  double twinmin, twinmax; // min, max 'good' time; values outside this range have flag set
82  double rhtime = hbhe.time();
83  double energy = hbhe.energy();
84 
85  if (energy < (double)tfilterEnvelope_.begin()->first) // less than lowest energy threshold
86  {
87  // Can skip timing test on cells below lowest threshold if so desired
88  if (ignorelowest_)
89  return 0;
90  else {
91  twinmin = tfilterEnvelope_.begin()->second.first;
92  twinmax = tfilterEnvelope_.begin()->second.second;
93  }
94  } else {
95  // Loop over energies in tfilterEnvelope
96  TfilterEnvelope_t::const_iterator it;
97  for (it = tfilterEnvelope_.begin(); it != tfilterEnvelope_.end(); ++it) {
98  // Identify tfilterEnvelope index for this rechit energy
99  if (energy < (double)it->first)
100  break;
101  }
102 
103  if (it == tfilterEnvelope_.end()) {
104  // Skip timing test on cells above highest threshold if so desired
105  if (ignorehighest_)
106  return 0;
107  else {
108  twinmin = tfilterEnvelope_.rbegin()->second.first;
109  twinmax = tfilterEnvelope_.rbegin()->second.second;
110  }
111  } else {
112  // Perform linear interpolation between energy boundaries
113 
114  std::map<int, std::pair<double, double> >::const_iterator prev = it;
115  prev--;
116 
117  // twinmax interpolation
118  double energy1 = prev->first;
119  double lim1 = prev->second.second;
120  double energy2 = it->first;
121  double lim2 = it->second.second;
122 
123  twinmax = lim1 + ((lim2 - lim1) * (energy - energy1) / (energy2 - energy1));
124 
125  // twinmin interpolation
126  lim1 = prev->second.first;
127  lim2 = it->second.first;
128 
129  twinmin = lim1 + ((lim2 - lim1) * (energy - energy1) / (energy2 - energy1));
130  }
131  }
132 
133  // Apply offset, gain
134  twinmin = win_offset_ + twinmin * win_gain_;
135  twinmax = win_offset_ + twinmax * win_gain_;
136 
137  // Set status high if time outside expected range
138  if (rhtime <= twinmin || rhtime >= twinmax)
139  status = 1; // set status to 1
140 
141  return status;
142 }
143 
145  int status = timingStatus(hbhe);
146 
147  // Though we're only using a single bit right now, 3 bits are reserved for these cuts
149 
150  return;
151 }
int timingStatus(const HBHERecHit &hbhe)
TfilterEnvelope_t tfilterEnvelope_
void SetTimingShapedFlags(HBHERecHit &hbhe)
void makeTfilterEnvelope(const std::vector< double > &v_userEnvelope)