CMS 3D CMS Logo

Stage1Layer2EtSumAlgorithmImpHI.cc
Go to the documentation of this file.
1 
17 #include <cassert>
18 
20 {
21  //now do what ever initialization is needed
22  for(size_t i=0; i<cordicPhiValues.size(); ++i) {
23  cordicPhiValues[i] = static_cast<int>(pow(2.,16)*(((float) i)-36)*M_PI/36);
24  }
25  for(size_t i=0; i<sines.size(); ++i) {
26  sines[i] = static_cast<long>(pow(2,30)*sin(i*20*M_PI/180));
27  cosines[i] = static_cast<long>(pow(2,30)*cos(i*20*M_PI/180));
28  }
29 }
30 
31 
32 void l1t::Stage1Layer2EtSumAlgorithmImpHI::processEvent(const std::vector<l1t::CaloRegion> & regions,
33  const std::vector<l1t::CaloEmCand> & EMCands,
34  const std::vector<l1t::Jet> * jets,
35  std::vector<l1t::EtSum> * etsums) {
36 
37 
38  //Region Correction will return uncorrected subregions if
39  //regionPUSType is set to None in the config
40  //double jetLsb=params_->jetLsb();
41  double jetLsb = 0.5; // HI O2O does not set this, and it will never change.
42 
43  //int etSumEtaMinEt = params_->etSumEtaMin(0);
44  //int etSumEtaMaxEt = params_->etSumEtaMax(0);
45  //double etSumEtThresholdEt = params_->etSumEtThreshold(0);
46  int etSumEtThresholdEt = (int) (params_->etSumEtThreshold(0) / jetLsb);
47 
48  //int etSumEtaMinHt = params_->etSumEtaMin(1);
49  //int etSumEtaMaxHt = params_->etSumEtaMax(1);
50  //double etSumEtThresholdHt = params_->etSumEtThreshold(1);
51  int etSumEtThresholdHt = (int) (params_->etSumEtThreshold(1) / jetLsb);
52 
53  // These values are not changeable online. O2O code for HI does not set a default,
54  // so previous results were garbage.
55  // Boundaries of 4 and 17 correspond to all non-HF regions.
56  int etSumEtaMinEt = 4;
57  int etSumEtaMaxEt = 17;
58  int etSumEtaMinHt = 4;
59  int etSumEtaMaxHt = 17;
60 
61  //RegionCorrection(regions, subRegions, params_);
62 
63  std::vector<SimpleRegion> regionEtVect;
64  std::vector<SimpleRegion> regionHtVect;
65 
66  // check the un-subtracted regions for overflow
67  bool regionOverflowEt(false);
68  bool regionOverflowHt(false);
69  for (auto& region : regions) {
70  if(region.hwEta() >= etSumEtaMinEt && region.hwEta() <= etSumEtaMaxEt)
71  {
72  if(region.hwPt() >= 1023)
73  {
74  regionOverflowEt = true;
75  }
76  }
77  if ( region.hwEta() >= etSumEtaMinHt && region.hwEta() <= etSumEtaMaxHt)
78  {
79  if(region.hwPt() >= 1023)
80  {
81  regionOverflowHt = true;
82  }
83  }
84  }
85 
86  // hwPt() is the sum ET+HT in region, for stage 1 this will be
87  // the region sum input to MET algorithm
88  // In stage 2, we would move to hwEtEm() and hwEtHad() for separate MET/MHT
89  // Thresholds will be hardware values not physical
90  //for (auto& region : *subRegions) {
91  for (auto& region : regions) {
92  if ( region.hwEta() >= etSumEtaMinEt && region.hwEta() <= etSumEtaMaxEt)
93  {
94  if(region.hwPt() >= etSumEtThresholdEt)
95  {
97  r.ieta = region.hwEta();
98  r.iphi = region.hwPhi();
99  r.et = region.hwPt();
100  regionEtVect.push_back(r);
101  }
102  }
103  if ( region.hwEta() >= etSumEtaMinHt && region.hwEta() <= etSumEtaMaxHt)
104  {
105  if(region.hwPt() >= etSumEtThresholdHt)
106  {
107  SimpleRegion r;
108  r.ieta = region.hwEta();
109  r.iphi = region.hwPhi();
110  r.et = region.hwPt();
111  regionHtVect.push_back(r);
112  }
113  }
114  }
115 
116  int sumET, MET, iPhiET;
117  std::tie(sumET, MET, iPhiET) = doSumAndMET(regionEtVect, ETSumType::kEmSum);
118 
119  int sumHT, MHT, iPhiHT;
120  std::tie(sumHT, MHT, iPhiHT) = doSumAndMET(regionHtVect, ETSumType::kHadronicSum);
121 
122  // Set quality (i.e. overflow) bits appropriately
123  int METqual = 0;
124  int MHTqual = 0;
125  int ETTqual = 0;
126  int HTTqual = 0;
127  if(MET >= 0xfff || regionOverflowEt) // MET 12 bits
128  METqual = 1;
129  if(MHT >= 0x7f || regionOverflowHt) // MHT 7 bits
130  MHTqual = 1;
131  if(sumET >= 0xfff || regionOverflowEt)
132  ETTqual = 1;
133  if(sumHT >= 0xfff || regionOverflowHt)
134  HTTqual = 1;
135 
136  MHT &= 127; // limit MHT to 7 bits as the firmware does, but only after checking for overflow.
137  //MHT is replaced with MHT/HT
138  uint16_t MHToHT=MHToverHT(MHT,sumHT);
139  //iPhiHt is replaced by the dPhi between two most energetic jets
140  iPhiHT = DiJetPhi(jets);
141 
142 
143  const ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > etLorentz(0,0,0,0);
144  l1t::EtSum etMiss(*&etLorentz,EtSum::EtSumType::kMissingEt,MET&0xfff,0,iPhiET,METqual);
145  l1t::EtSum htMiss(*&etLorentz,EtSum::EtSumType::kMissingHt,MHToHT&0x7f,0,iPhiHT,MHTqual);
146  l1t::EtSum etTot (*&etLorentz,EtSum::EtSumType::kTotalEt,sumET&0xfff,0,0,ETTqual);
147  l1t::EtSum htTot (*&etLorentz,EtSum::EtSumType::kTotalHt,sumHT&0xfff,0,0,HTTqual);
148 
149  std::vector<l1t::EtSum> preGtEtSums;
150  preGtEtSums.reserve(4);
151 
152  preGtEtSums.push_back(etMiss);
153  preGtEtSums.push_back(htMiss);
154  preGtEtSums.push_back(etTot);
155  preGtEtSums.push_back(htTot);
156 
157  EtSumToGtScales(params_, &preGtEtSums, etsums);
158 
159 }
160 
161 std::tuple<int, int, int>
162 l1t::Stage1Layer2EtSumAlgorithmImpHI::doSumAndMET(const std::vector<SimpleRegion>& regionEt, ETSumType sumType)
163 {
164  std::array<int, 18> sumEtaPos{};
165  std::array<int, 18> sumEtaNeg{};
166  for (const auto& r : regionEt)
167  {
168  if ( r.ieta < 11 )
169  sumEtaNeg[r.iphi] += r.et;
170  else
171  sumEtaPos[r.iphi] += r.et;
172  }
173 
174  std::array<int, 18> sumEta{};
175  int sumEt(0);
176  for(size_t i=0; i<sumEta.size(); ++i)
177  {
178  assert(sumEtaPos[i] >= 0 && sumEtaNeg[i] >= 0);
179  sumEta[i] = sumEtaPos[i] + sumEtaNeg[i];
180  sumEt += sumEta[i];
181  }
182 
183  // 0, 20, 40, 60, 80 degrees
184  std::array<int, 5> sumsForCos{};
185  std::array<int, 5> sumsForSin{};
186  for(size_t iphi=0; iphi<sumEta.size(); ++iphi)
187  {
188  if ( iphi < 5 )
189  {
190  sumsForCos[iphi] += sumEta[iphi];
191  sumsForSin[iphi] += sumEta[iphi];
192  }
193  else if ( iphi < 9 )
194  {
195  sumsForCos[9-iphi] -= sumEta[iphi];
196  sumsForSin[9-iphi] += sumEta[iphi];
197  }
198  else if ( iphi < 14 )
199  {
200  sumsForCos[iphi-9] -= sumEta[iphi];
201  sumsForSin[iphi-9] -= sumEta[iphi];
202  }
203  else
204  {
205  sumsForCos[18-iphi] += sumEta[iphi];
206  sumsForSin[18-iphi] -= sumEta[iphi];
207  }
208  }
209 
210  long sumX(0l);
211  long sumY(0l);
212  for(int i=0; i<5; ++i)
213  {
214  sumX += sumsForCos[i]*cosines[i];
215  sumY += sumsForSin[i]*sines[i];
216  }
217  assert(abs(sumX)<(1l<<48) && abs(sumY)<(1l<<48));
218  int cordicX = sumX>>25;
219  int cordicY = sumY>>25;
220 
221  uint32_t cordicMag(0);
222  int cordicPhase(0);
223  cordic(cordicX, cordicY, cordicPhase, cordicMag);
224 
225  int met(0);
226  int metPhi(0);
227  if ( sumType == ETSumType::kHadronicSum )
228  {
229  met = (cordicMag % (1<<7)) | ((cordicMag >= (1<<7)) ? (1<<7):0);
230  metPhi = cordicToMETPhi(cordicPhase) >> 2;
231  assert(metPhi >=0 && metPhi < 18);
232  }
233  else
234  {
235  met = (cordicMag % (1<<12)) | ((cordicMag >= (1<<12)) ? (1<<12):0);
236  metPhi = cordicToMETPhi(cordicPhase);
237  assert(metPhi >=0 && metPhi < 72);
238  }
239 
240  return std::make_tuple(sumEt, met, metPhi);
241 }
242 
243 // converts phase from 3Q16 to 0-71
244 // Expects abs(phase) <= 205887 (pi*2^16)
245 int
247 {
248  assert(abs(phase)<=205887);
249  for(size_t i=0; i<cordicPhiValues.size()-1; ++i)
250  if ( phase >= cordicPhiValues[i] && phase < cordicPhiValues[i+1] )
251  return i;
252  // if phase == +205887 (+pi), return zero
253  return 0;
254 }
255 
256 int l1t::Stage1Layer2EtSumAlgorithmImpHI::DiJetPhi(const std::vector<l1t::Jet> * jets) const {
257 
258 
259  int dphi = 10; // initialize to negative physical dphi value
260  if (jets->size()<2) return dphi; // size() not really reliable as we pad the size to 8 (4cen+4for) in the sorter
261  if ((*jets).at(0).hwPt() == 0) return dphi;
262  if ((*jets).at(1).hwPt() == 0) return dphi;
263 
264 
265  int iphi1 = (*jets).at(0).hwPhi();
266  int iphi2 = (*jets).at(1).hwPhi();
267 
268  int difference=abs(iphi1-iphi2);
269 
270  if ( difference > 9 ) difference= L1CaloRegionDetId::N_PHI - difference ; // make Physical dphi always positive
271  return difference;
272 }
273 
274 uint16_t l1t::Stage1Layer2EtSumAlgorithmImpHI::MHToverHT(uint16_t num,uint16_t den) const {
275 
276  uint16_t result;
277  uint32_t numerator(num),denominator(den);
278 
279  if(numerator == denominator)
280  result = 0x7f;
281  else
282  {
283  numerator = numerator << 7;
284  result = numerator/denominator;
285  result = result & 0x7f;
286  }
287 
288  return result;
289 }
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
double etSumEtThreshold(unsigned isum) const
void processEvent(const std::vector< l1t::CaloRegion > &regions, const std::vector< l1t::CaloEmCand > &EMCands, const std::vector< l1t::Jet > *jets, std::vector< l1t::EtSum > *sums) override
std::tuple< int, int, int > doSumAndMET(const std::vector< SimpleRegion > &regionEt, ETSumType sumType)
Stage1Layer2EtSumAlgorithmImpHI(CaloParamsHelper const *params)
int DiJetPhi(const std::vector< l1t::Jet > *jets) const
void EtSumToGtScales(CaloParamsHelper const *params, const std::vector< l1t::EtSum > *input, std::vector< l1t::EtSum > *output)
vector< PseudoJet > jets
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
#define M_PI
met
===> hadronic RAZOR
uint16_t MHToverHT(uint16_t, uint16_t) const
static const unsigned N_PHI
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40