CMS 3D CMS Logo

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