CMS 3D CMS Logo

Stage2Layer2JetAlgorithmFirmwareImp1.cc
Go to the documentation of this file.
1 //
7 
15 #include "TMath.h"
16 
17 #include <vector>
18 #include <algorithm>
19 #include <cmath>
20 
21 namespace l1t {
22  bool operator > ( const l1t::Jet& a, const l1t::Jet& b ) {
23  return a.hwPt() > b.hwPt();
24  }
25 }
26 
27 // jet mask, needs to be configurable at some point
28 // just a square for now
29 // for 1 do greater than, for 2 do greater than equal to
30 
31 namespace {
32 constexpr int mask_[9][9] = {
33  { 1,2,2,2,2,2,2,2,2 },
34  { 1,1,2,2,2,2,2,2,2 },
35  { 1,1,1,2,2,2,2,2,2 },
36  { 1,1,1,1,2,2,2,2,2 },
37  { 1,1,1,1,0,2,2,2,2 },
38  { 1,1,1,1,1,2,2,2,2 },
39  { 1,1,1,1,1,1,2,2,2 },
40  { 1,1,1,1,1,1,1,2,2 },
41  { 1,1,1,1,1,1,1,1,2 },
42 };
43 
44 }
45 
46 
48  params_(params){}
49 
50 
51 void l1t::Stage2Layer2JetAlgorithmFirmwareImp1::processEvent(const std::vector<l1t::CaloTower> & towers,
52  std::vector<l1t::Jet> & jets,
53  std::vector<l1t::Jet> & alljets) {
54 
55  // find jets
56  create(towers, jets, alljets, params_->jetPUSType());
57 
58  // calibrate all jets
59  calibrate(alljets, 0, true); // pass all jets and the hw threshold above which to calibrate
60 
61  // jets accumulated sort
62  accuSort(jets);
63 
64 }
65 
66 
67 void l1t::Stage2Layer2JetAlgorithmFirmwareImp1::create(const std::vector<l1t::CaloTower> & towers,
68  std::vector<l1t::Jet> & jets,
69  std::vector<l1t::Jet> & alljets,
70  std::string PUSubMethod) {
71 
72 
73  // etaSide=1 is positive eta, etaSide=-1 is negative eta
74  for (int etaSide=1; etaSide>=-1; etaSide-=2) {
75 
76  // the 4 groups of rings
77  std::vector<int> ringGroup1, ringGroup2, ringGroup3, ringGroup4;
78  for (int i=1; i<=CaloTools::mpEta(CaloTools::kHFEnd); i++) {
79  if ( ! ((i-1)%4) ) ringGroup1.push_back( i * etaSide );
80  else if ( ! ((i-2)%4) ) ringGroup2.push_back( i * etaSide );
81  else if ( ! ((i-3)%4) ) ringGroup3.push_back( i * etaSide );
82  else if ( ! ((i-4)%4) ) ringGroup4.push_back( i * etaSide );
83  }
84  std::vector< std::vector<int> > theRings = { ringGroup1, ringGroup2, ringGroup3, ringGroup4 };
85 
86  // the 24 jets in this eta side
87  std::vector<l1t::Jet> jetsHalf;
88 
89  // loop over the 4 groups of rings
90  for ( unsigned ringGroupIt=1; ringGroupIt<=theRings.size(); ringGroupIt++ ) {
91 
92  // the 6 accumulated jets
93  std::vector<l1t::Jet> jetsAccu;
94 
95  // loop over the 10 rings in this group
96  for ( unsigned ringIt=0; ringIt<theRings.at(ringGroupIt-1).size(); ringIt++ ) {
97 
98  int ieta = theRings.at(ringGroupIt-1).at(ringIt);
99 
100  // the jets in this ring
101  std::vector<l1t::Jet> jetsRing;
102 
103  // loop over phi in the ring
104  for ( int iphi=1; iphi<=CaloTools::kHBHENrPhi; ++iphi ) {
105 
106  // no more than 18 jets per ring
107  if (jetsRing.size()==18) break;
108 
109  // seed tower
110  const CaloTower& tow = CaloTools::getTower(towers, CaloTools::caloEta(ieta), iphi);
111 
112  int seedEt = tow.hwPt();
113  int iEt = seedEt;
114  bool satSeed = false;
115  bool vetoCandidate = false;
116 
117  // check it passes the seed threshold
118  if(iEt < floor(params_->jetSeedThreshold()/params_->towerLsbSum())) continue;
119 
120  if (seedEt == CaloTools::kSatHcal ||
121  seedEt == CaloTools::kSatEcal ||
122  seedEt == CaloTools::kSatTower)
123  satSeed=true;
124 
125  // loop over towers in this jet
126  for( int deta = -4; deta < 5; ++deta ) {
127  for( int dphi = -4; dphi < 5; ++dphi ) {
128 
129  int towEt = 0;
130  int ietaTest = ieta+deta;
131  int iphiTest = iphi+dphi;
132 
133  // wrap around phi
134  while ( iphiTest > CaloTools::kHBHENrPhi ) iphiTest -= CaloTools::kHBHENrPhi;
135  while ( iphiTest < 1 ) iphiTest += CaloTools::kHBHENrPhi;
136 
137  // wrap over eta=0
138  if (ieta > 0 && ietaTest <=0) ietaTest -= 1;
139  if (ieta < 0 && ietaTest >=0) ietaTest += 1;
140 
141  // check jet mask and sum tower et
142  const CaloTower& towTest = CaloTools::getTower(towers, CaloTools::caloEta(ietaTest), iphiTest);
143  towEt = towTest.hwPt();
144 
145  if (mask_[8-(dphi+4)][deta+4] == 0) continue;
146  else if (mask_[8-(dphi+4)][deta+4] == 1) vetoCandidate = (seedEt < towEt);
147  else if (mask_[8-(dphi+4)][deta+4] == 2) vetoCandidate = (seedEt <= towEt);
148 
149  if (vetoCandidate) break;
150  else iEt += towEt;
151 
152  }
153  if(vetoCandidate) break;
154  }
155 
156  // add the jet to the list
157  if (!vetoCandidate) {
158 
159  int rawEt = iEt;
160  int puEt(0);
161 
163  int caloEta = CaloTools::caloEta(ieta);
164  l1t::Jet jet( p4, -999, caloEta, iphi, 0);
165 
166  if(!params_->jetBypassPUS()
167  && !satSeed){
168 
169  // comment out for emulator studies
171  PUSubMethod = "PhiRing1";
172  else
173  PUSubMethod = "ChunkyDonut";
174 
175  if (PUSubMethod == "Donut") {
176  puEt = donutPUEstimate(ieta, iphi, 5, towers);
177  iEt -= puEt;
178  }
179 
180  if (PUSubMethod == "ChunkyDonut"){
181  puEt = chunkyDonutPUEstimate(jet, 5, towers);
182  iEt -= puEt;
183  }
184 
185  if(PUSubMethod.find("ChunkySandwich") != std::string::npos){
186  puEt = chunkySandwichPUEstimate(jet, 5, towers, PUSubMethod);
187  iEt -= puEt;
188  }
189 
190 
191  if(PUSubMethod == "PhiRing1"){
192  int size=5;
193  std::map<int,int> SumEtEtaMap=getSumEtEtaMap(towers);
194  // int jetPhi=jet.hwPhi();
195  int jetEta=CaloTools::mpEta(jet.hwEta());
196  int jetEtaLow=jetEta-size+1;
197  int jetEtaHigh=jetEta+size-1;
198  if(jetEta>0 && jetEtaLow<=0) jetEtaLow-=1;
199  if(jetEta<0 && jetEtaHigh>=0) jetEtaHigh+=1;
200  puEt = 0;
201  for(int ieta=jetEtaLow; ieta<=jetEtaHigh; ++ieta){
202  if(ieta==0) continue; // eta start from +- 1
203  auto iter=SumEtEtaMap.find(ieta);
204  puEt+=iter->second;
205  }
206  puEt=(puEt*(size*2-1))/CaloTools::kHBHENrPhi;
207  iEt -= puEt;
208  }
209 
210  if(PUSubMethod == "PhiRing2"){
211  int size=5;
212  std::map<int,int> SumEtEtaMap=getSumEtEtaMap(towers);
213  // int jetPhi=jet.hwPhi();
214  int jetEta=CaloTools::mpEta(jet.hwEta());
215  int jetEtaLow=jetEta-size+1;
216  int jetEtaHigh=jetEta+size-1;
217 
218  if(jetEta>0 && jetEtaLow<=0) jetEtaLow-=1;
219  if(jetEta<0 && jetEtaHigh>=0) jetEtaHigh+=1;
220  puEt = 0;
221 
222  for(int ieta=jetEtaLow; ieta<=jetEtaHigh; ++ieta){
223  if(ieta==0) continue;
224  auto iter=SumEtEtaMap.find(ieta);
225  puEt+=iter->second;
226  }
227  puEt-=iEt;
228  puEt=(puEt*(size*2-1))/(CaloTools::kHBHENrPhi - (size*2-1));
229  iEt -= puEt;
230  }
231 
232 
233 
234  }
235 
236  if (iEt<=0) continue;
237 
238  // if tower Et is saturated, saturate jet Et
239  if (seedEt == CaloTools::kSatHcal || seedEt == CaloTools::kSatEcal || seedEt == CaloTools::kSatTower) iEt = CaloTools::kSatJet;
240 
241  jet.setHwPt(iEt);
242  jet.setRawEt( (short int) rawEt);
243  jet.setSeedEt((short int) seedEt);
244  jet.setTowerIEta((short int) caloEta);
245  jet.setTowerIPhi((short int) iphi);
246  jet.setPUEt((short int) puEt);
247 
248 
249  jetsRing.push_back(jet);
250  alljets.push_back(jet);
251 
252  }
253 
254  }
255 
256  // jet energy corrections
257  calibrate(jetsRing, 0, false); // pass the jet collection and the hw threshold above which to calibrate
258 
259  // sort these jets and keep top 6
260  auto start_ = jetsRing.begin();
261  auto end_ = jetsRing.end();
262  BitonicSort<l1t::Jet>(down, start_, end_);
263  if (jetsRing.size()>6) jetsRing.resize(6);
264 
265  // update jets
266  jets.insert(jets.end(),jetsRing.begin(),jetsRing.end());
267 
268  }
269  }
270  }
271 }
272 
273 
274 //Accumulating sort
276 
278  l1t::Jet tempJet (emptyP4, 0, 0, 0, 0);
279  std::vector< std::vector<l1t::Jet> > jetEtaPos( 41 , std::vector<l1t::Jet>(18, tempJet));
280  std::vector< std::vector<l1t::Jet> > jetEtaNeg( 41 , std::vector<l1t::Jet>(18, tempJet));
281 
282  for (unsigned int iJet = 0; iJet < jets.size(); iJet++)
283  {
284  if (jets.at(iJet).hwEta() > 0) jetEtaPos.at(jets.at(iJet).hwEta()-1).at((72-jets.at(iJet).hwPhi())/4) = jets.at(iJet);
285  else jetEtaNeg.at(-(jets.at(iJet).hwEta()+1)).at((72-jets.at(iJet).hwPhi())/4) = jets.at(iJet);
286  }
287 
288  AccumulatingSort <l1t::Jet> etaPosSorter(7);
289  AccumulatingSort <l1t::Jet> etaNegSorter(7);
290  std::vector<l1t::Jet> accumEtaPos;
291  std::vector<l1t::Jet> accumEtaNeg;
292 
293  for( int ieta = 0 ; ieta < 41 ; ++ieta)
294  {
295  // eta +
296  std::vector<l1t::Jet>::iterator start_, end_;
297  start_ = jetEtaPos.at(ieta).begin();
298  end_ = jetEtaPos.at(ieta).end();
299  BitonicSort<l1t::Jet>(down, start_, end_);
300  etaPosSorter.Merge( jetEtaPos.at(ieta) , accumEtaPos );
301 
302  // eta -
303  start_ = jetEtaNeg.at(ieta).begin();
304  end_ = jetEtaNeg.at(ieta).end();
305  BitonicSort<l1t::Jet>(down, start_, end_);
306  etaNegSorter.Merge( jetEtaNeg.at(ieta) , accumEtaNeg );
307 
308  }
309 
310  //check for 6 & 7th jets with same et and eta. Keep jet with larger phi
311  if(accumEtaPos.at(6).hwPt()==accumEtaPos.at(5).hwPt() && accumEtaPos.at(6).hwEta()==accumEtaPos.at(5).hwEta()
312  && accumEtaPos.at(6).hwPhi() > accumEtaPos.at(5).hwPhi()){
313  accumEtaPos.at(5)=accumEtaPos.at(6);
314  }
315  if(accumEtaNeg.at(6).hwPt()==accumEtaNeg.at(5).hwPt() && accumEtaNeg.at(6).hwEta()==accumEtaNeg.at(5).hwEta()
316  && accumEtaNeg.at(6).hwPhi() > accumEtaNeg.at(5).hwPhi()){
317  accumEtaNeg.at(5)=accumEtaNeg.at(6);
318  }
319 
320  //truncate
321  accumEtaPos.resize(6);
322  accumEtaNeg.resize(6);
323 
324  // put all 12 candidates in the original jet vector, removing zero energy ones
325  jets.clear();
326  for (l1t::Jet accjet : accumEtaPos)
327  {
328  if (accjet.hwPt() > 0) jets.push_back(accjet);
329  }
330  for (l1t::Jet accjet : accumEtaNeg)
331  {
332  if (accjet.hwPt() > 0) jets.push_back(accjet);
333  }
334 
335 
336 }
337 
338 
339 
340 //A function to return the value for donut subtraction around an ieta and iphi position for donut subtraction
341 //Also pass it a vector to store the individual values of the strip for later testing
342 //The size is the number of ieta/iphi units out the ring is (ie for 9x9 jets, we want the 11x11 for PUS therefore we want to go 5 out, so size is 5)
344  int jetPhi,
345  int size,
346  const std::vector<l1t::CaloTower> & towers){
347 
348  //ring is a vector with 4 ring strips, one for each side of the ring
349  std::vector<int> ring(4,0);
350 
351  int iphiUp = jetPhi + size;
352  while ( iphiUp > CaloTools::kHBHENrPhi ) iphiUp -= CaloTools::kHBHENrPhi;
353  int iphiDown = jetPhi - size;
354  while ( iphiDown < 1 ) iphiDown += CaloTools::kHBHENrPhi;
355 
356  int ietaUp = jetEta+size; //(jetEta + size > CaloTools::mpEta(CaloTools::kHFEnd)) ? 999 : jetEta+size;
357  int ietaDown = jetEta-size; //(abs(jetEta - size) > CaloTools::mpEta(CaloTools::kHFEnd)) ? 999 : jetEta-size;
358 
359  for (int ieta = jetEta - size+1; ieta < jetEta + size; ++ieta)
360  {
361 
362  if (abs(ieta) > CaloTools::mpEta(CaloTools::kHFEnd) || abs(ieta) < 1) continue;
363  int towerEta;
364 
365  if (jetEta > 0 && ieta <=0){
366  towerEta = ieta-1;
367  } else if (jetEta < 0 && ieta >=0){
368  towerEta = ieta+1;
369  } else {
370  towerEta=ieta;
371  }
372 
373  const CaloTower& tow = CaloTools::getTower(towers, CaloTools::caloEta(towerEta), iphiUp);
374  int towEt = tow.hwPt();
375  ring[0]+=towEt;
376 
377  const CaloTower& tow2 = CaloTools::getTower(towers, CaloTools::caloEta(towerEta), iphiDown);
378  towEt = tow2.hwPt();
379  ring[1]+=towEt;
380 
381  }
382 
383  for (int iphi = jetPhi - size+1; iphi < jetPhi + size; ++iphi)
384  {
385 
386  int towerPhi = iphi;
387  while ( towerPhi > CaloTools::kHBHENrPhi ) towerPhi -= CaloTools::kHBHENrPhi;
388  while ( towerPhi < 1 ) towerPhi += CaloTools::kHBHENrPhi;
389 
390  const CaloTower& tow = CaloTools::getTower(towers, CaloTools::caloEta(ietaUp), towerPhi);
391  int towEt = tow.hwPt();
392  ring[2]+=towEt;
393 
394  const CaloTower& tow2 = CaloTools::getTower(towers, CaloTools::caloEta(ietaDown), towerPhi);
395  towEt = tow2.hwPt();
396  ring[3]+=towEt;
397  }
398 
399  //for the Donut Subtraction we only use the middle 2 (in energy) ring strips
400  std::sort(ring.begin(), ring.end(), std::greater<int>());
401 
402  return 4*( ring[1]+ring[2] ); // This should really be multiplied by 4.5 not 4.
403 }
404 
405 std::map<int,int> l1t::Stage2Layer2JetAlgorithmFirmwareImp1::getSumEtEtaMap(const std::vector<l1t::CaloTower> & towers){
406 
407 
408  std::map<int,int> SumEtEtaMap;
409  int EtaMin=-41;
410  int EtaMax=41;
411  int phiMin=1;
412  int phiMax=72;
413 
414  for(int iEta=EtaMin ; iEta<=EtaMax; iEta++ )
415  {
416  int SumEt=0;
417  for(int iPhi=phiMin; iPhi<=phiMax; iPhi++){
418  const CaloTower& tow = CaloTools::getTower(towers, CaloTools::caloEta(iEta), iPhi);
419  int towEt=tow.hwPt();
420  SumEt+=towEt;
421 
422  }// end for iPih<=phiMax
423  SumEtEtaMap[iEta]=SumEt;
424 
425  }// end for iEta<=EtaMax
426 
427  return SumEtEtaMap;
428 
429 }
430 
431 
433  const std::vector<l1t::CaloTower> & towers,
434  const std::string chunkyString)
435 {
436  int jetPhi = jet.hwPhi();
437  int jetEta = CaloTools::mpEta(jet.hwEta());
438 
439  // ring is a vector with 4 or 8 eta or strips, one for each side of the ring in ChunkyDonut, else pure phi
440  int ringSize = 4;
441  if(chunkyString == "ChunkySandwich8") ringSize = 8;
442 
443  std::vector<int> ring(ringSize,0);
444 
445  // number of strips in donut - should make this configurable
446  int nStrips = 3;
447 
448  // loop over strips
449  for (int stripIt=0; stripIt<nStrips; stripIt++) {
450 
451  int iphiUp = jetPhi + size + stripIt;
452  int iphiDown = jetPhi - size - stripIt;
453  while ( iphiUp > CaloTools::kHBHENrPhi ) iphiUp -= CaloTools::kHBHENrPhi;
454  while ( iphiDown < 1 ) iphiDown += CaloTools::kHBHENrPhi;
455 
456  // we will always do phiUp and PhiDown
457  // do PhiUp and PhiDown
458  for (int ieta=jetEta-size+1; ieta<jetEta+size; ++ieta) {
459 
460  if (abs(ieta) > CaloTools::mpEta(CaloTools::kHFEnd)) continue;
461 
462  int towEta = ieta;
463  if (jetEta>0 && towEta<=0) towEta-=1;
464  if (jetEta<0 && towEta>=0) towEta+=1;
465 
466  const CaloTower& towPhiUp = CaloTools::getTower(towers, CaloTools::caloEta(towEta), iphiUp);
467  int towEt = towPhiUp.hwPt();
468  ring[0] += towEt;
469 
470  const CaloTower& towPhiDown = CaloTools::getTower(towers, CaloTools::caloEta(towEta), iphiDown);
471  towEt = towPhiDown.hwPt();
472  ring[1] += towEt;
473 
474  }
475 
476  //Only consider eta for chunkydonut
477  if(chunkyString == "ChunkyDonut"){
478  int ietaUp = jetEta + size + stripIt;
479  int ietaDown = jetEta - size - stripIt;
480  if ( jetEta<0 && ietaUp>=0 ) ietaUp += 1;
481  if ( jetEta>0 && ietaDown<=0 ) ietaDown -= 1;
482 
483  // do EtaUp
484  for (int iphi=jetPhi-size+1; iphi<jetPhi+size; ++iphi) {
485 
486  if (abs(ietaUp) <= CaloTools::mpEta(CaloTools::kHFEnd)) {
487  int towPhi = iphi;
488  while ( towPhi > CaloTools::kHBHENrPhi ) towPhi -= CaloTools::kHBHENrPhi;
489  while ( towPhi < 1 ) towPhi += CaloTools::kHBHENrPhi;
490 
491  const CaloTower& towEtaUp = CaloTools::getTower(towers, CaloTools::caloEta(ietaUp), towPhi);
492  int towEt = towEtaUp.hwPt();
493  ring[2] += towEt;
494  }
495  }
496 
497  // do EtaDown
498  for (int iphi=jetPhi-size+1; iphi<jetPhi+size; ++iphi) {
499 
500  if (abs(ietaDown) <= CaloTools::mpEta(CaloTools::kHFEnd)) {
501  int towPhi = iphi;
502  while ( towPhi > CaloTools::kHBHENrPhi ) towPhi -= CaloTools::kHBHENrPhi;
503  while ( towPhi < 1 ) towPhi += CaloTools::kHBHENrPhi;
504 
505  const CaloTower& towEtaDown = CaloTools::getTower(towers, CaloTools::caloEta(ietaDown), towPhi);
506  int towEt = towEtaDown.hwPt();
507  ring[3] += towEt;
508  }
509  }
510  }
511  else if(chunkyString == "ChunkySandwich2"){//simplest of the sandwiches - we will just reuse phi flaps
512  for(int i = 0; i < 2; ++i){ring[i+2] = ring[i];}
513  }
514  else if(chunkyString == "ChunkySandwich2Ave" || chunkyString == "ChunkySandwich"){//simplest of the sandwiches - we will just reuse phi flaps
515  ring[2] = ((ring[0] + ring[1]) >> 1); //bitwise, effective div by 2
516  ring[3] = (ring[0] + ring[1]); // just make sure it is >= max val
517  }
518  else if(chunkyString == "ChunkySandwich4" || chunkyString == "ChunkySandwich8"){//Need to calculate two or siz additional phiflaps. Note that ChunkySandwich defaults to ChunkySandwich4 variant
519 
520  for(int rI = 0; rI < (ringSize-2)/2; ++rI){
521  int iphiUp2 = jetPhi + size + (stripIt + nStrips*(rI+1));
522  int iphiDown2 = jetPhi - size - (stripIt + nStrips*(rI+1));
523  while ( iphiUp2 > CaloTools::kHBHENrPhi ) iphiUp2 -= CaloTools::kHBHENrPhi;
524  while ( iphiDown2 < 1 ) iphiDown2 += CaloTools::kHBHENrPhi;
525 
526  for (int ieta=jetEta-size+1; ieta<jetEta+size; ++ieta) {
527  if (abs(ieta) > CaloTools::mpEta(CaloTools::kHFEnd)) continue;
528 
529  int towEta = ieta;
530  if (jetEta>0 && towEta<=0) towEta-=1;
531  if (jetEta<0 && towEta>=0) towEta+=1;
532 
533  const CaloTower& towPhiUp = CaloTools::getTower(towers, CaloTools::caloEta(towEta), iphiUp2);
534  int towEt = towPhiUp.hwPt();
535  ring[2 + rI*2] += towEt;
536 
537  const CaloTower& towPhiDown = CaloTools::getTower(towers, CaloTools::caloEta(towEta), iphiDown2);
538  towEt = towPhiDown.hwPt();
539  ring[3 + rI*2] += towEt;
540  }
541  }
542  }
543  }
544 
545  //sort our final ring and return;
546  std::sort( ring.begin(), ring.end() );
547 
548  return ring;
549 }
550 
551 
553  const std::vector<l1t::CaloTower> & towers){
554 
555  // ring is a vector with 4 ring strips, one for each side of the ring
556  // PhiUp, PhiDown, EtaUp, EtaDown, sorted w/ call to chunkyring
557  std::vector<int> ring = getChunkyRing(jet, size, towers, "ChunkyDonut");
558  for(unsigned int i=0; i<4; ++i) jet.setPUDonutEt(i, (short int) ring[i]);
559  return ( ring[0] + ring[1] + ring[2] );
560 }
561 
562 //Following example of chunky donut, but considering only phi flaps
563 //Useful for simple handling of zeroed ieta strips and for high PU environments like PbPb
565  const std::vector<l1t::CaloTower> & towers,
566  const std::string chunkySandwichStr){
567  //ring will be handled identically by all variants for now - could consider a different picking w/ ChunkySandwich8, say to exclude a downward fluctuation w/ 1,2,3 from ring, etc.
568  std::vector<int> ring = getChunkyRing(jet, size, towers, chunkySandwichStr);
569  for(unsigned int i=0; i<4; ++i) jet.setPUDonutEt(i, (short int)ring[i]);
570  return (ring[0] + ring[1] + ring[2]);
571 }
572 
573 
574 void l1t::Stage2Layer2JetAlgorithmFirmwareImp1::calibrate(std::vector<l1t::Jet> & jets, int calibThreshold, bool isAllJets) {
575 
576  if( params_->jetCalibrationType() == "function6PtParams22EtaBins" ){ //One eta bin per region
577 
578  //Check the vector of calibration parameters is the correct size
579  //Order the vector in terms of the parameters per eta bin, starting in -ve eta
580  //So first 6 entries are very negative eta, next 6 are the next bin etc.
581 
582  if( params_->jetCalibrationParams().size() != 6*22){
583  edm::LogError("l1t|stage 2") << "Invalid input vector to calo params. Input vector of size: " <<
584  params_->jetCalibrationParams().size() << " Require size: 132 Not calibrating Stage 2 Jets" << std::endl;
585  return;
586  }
587 
588  //Loop over jets and apply corrections
589  for(std::vector<l1t::Jet>::iterator jet = jets.begin(); jet!=jets.end(); jet++){
590 
591  //Check jet is above the calibration threshold, if not do nothing
592  if(jet->hwPt() < calibThreshold) continue;
593  if(jet->hwPt() >= 0xFFFF) continue;
594 
595  int etaBin = CaloTools::regionEta( jet->hwEta() );
596 
597  //Get the parameters from the vector
598  //Each 6 values are the parameters for an eta bin
599  double params[6];
600  for(int i=0; i<6; i++){
601  params[i] = params_->jetCalibrationParams()[etaBin*6 + i];
602  }
603 
604  //Perform the correction based on the calibration function defined
605  //in calibFit
606  //This is derived from the actual physical pt of the jets, not the hwEt
607  //This needs to be addressed in the future
608  double ptPhys = jet->hwPt() * params_->jetLsb();
609  double correction = calibFit(ptPhys, params);
610 
611  jet->setHwPt(correction*jet->hwPt());
612 
613  }
614 
615  }
616  else if( params_->jetCalibrationType() == "function8PtParams22EtaBins" ){
617  // as above but with cap on max correction at low pT
618 
619  if( params_->jetCalibrationParams().size() != 8*22){
620  edm::LogError("l1t|stage 2") << "Invalid input vector to calo params. Input vector of size: " <<
621  params_->jetCalibrationParams().size() << " Require size: 176 Not calibrating Stage 2 Jets" << std::endl;
622  return;
623  }
624 
625  for(std::vector<l1t::Jet>::iterator jet = jets.begin(); jet!=jets.end(); jet++){
626 
627  if(jet->hwPt() < calibThreshold) continue;
628  if(jet->hwPt() >= 0xFFFF) continue;
629 
630  int etaBin = CaloTools::regionEta( jet->hwEta() );
631 
632  double params[8];
633  for(int i=0; i<8; i++){
634  params[i] = params_->jetCalibrationParams()[etaBin*8 + i];
635  }
636 
637  double ptPhys = jet->hwPt() * params_->jetLsb();
638  double correction = params[6];
639  if (ptPhys>params[7]) correction = calibFit(ptPhys, params);
640 
641  jet->setHwPt(correction*jet->hwPt());
642 
643  }
644 
645  }
646  else if( params_->jetCalibrationType() == "functionErf11PtParams16EtaBins" ){
647  // as above but with cap on max correction at low pT
648 
649  if( params_->jetCalibrationParams().size() != 11*16){
650  edm::LogError("l1t|stage 2") << "Invalid input vector to calo params. Input vector of size: " <<
651  params_->jetCalibrationParams().size() << " Require size: 176 Not calibrating Stage 2 Jets" << std::endl;
652  return;
653  }
654 
655  for(std::vector<l1t::Jet>::iterator jet = jets.begin(); jet!=jets.end(); jet++){
656 
657  if(jet->hwPt() < calibThreshold) continue;
658  if(jet->hwPt() >= 0xFFFF) continue;
659 
660  int etaBin = CaloTools::bin16Eta( jet->hwEta() );
661 
662  double params[11];
663  for(int i=0; i<11; i++){
664  params[i] = params_->jetCalibrationParams()[etaBin*11 + i];
665  }
666 
667  double ptPhys = jet->hwPt() * params_->jetLsb();
668  double correction = params[7];
669 
670  if (ptPhys<params[8]) correction = params[7];
671  else if (ptPhys>params[10]) correction = params[9];
672  else correction = calibFitErr(ptPhys, params);
673 
674  jet->setHwPt(correction*jet->hwPt());
675 
676  }
677 
678  }
679  else if( params_->jetCalibrationType() == "LUT" ){
680  // Calibrate using 3 LUTs: pt and eta compression LUTs, and a multiplicand/addend LUT.
681  // The pt and eta are each converted to a compressed scale using individual LUTs
682  // pt : 8 -> 4 bits, eta 6 -> 4 bits
683  // This then forms an address. Using the third LUT, we get a
684  // multiplicand & addend, so we can do y = m*x + c on the original
685  // (i.e. non-compressed) jet pt.
686  // The multiplicand is 10-bit unsigned, addend is 8-bit signed.
687 
688  //Loop over jets and apply corrections
689  for(std::vector<l1t::Jet>::iterator jet = jets.begin(); jet!=jets.end(); jet++){
690 
691  //Check jet is above the calibration threshold, if not do nothing
692  if(jet->hwPt() < calibThreshold) continue;
693 
694  //don't calibrate saturated jets for HT
695  if( isAllJets && (jet->hwPt() == CaloTools::kSatJet) ) continue;
696 
697  // In the firmware, we take bits 1 to 8 of the hwPt.
698  // To avoid getting nonsense by only taking smaller bits,
699  // any values larger than 511 are automatically set to 511.
700  int jetHwPt = jet->hwPt();
701  if (jetHwPt >= 0x200) {
702  jetHwPt = 0x1FF;
703  }
704  unsigned int ptBin = params_->jetCompressPtLUT()->data(jetHwPt>>1);
705  unsigned int etaBin = params_->jetCompressEtaLUT()->data(abs(CaloTools::mpEta(jet->hwEta())));
706  unsigned int compBin = (etaBin<<4) | ptBin;
707 
708  unsigned int addPlusMult = params_->jetCalibrationLUT()->data(compBin);
709  unsigned int multiplier = addPlusMult & 0x3ff;
710  // handles -ve numbers correctly
711  int8_t addend = (addPlusMult>>10);
712  unsigned int jetPtCorr = ((jet->hwPt()*multiplier)>>9) + addend;
713 
714  if(jetPtCorr < 0xFFFF) {
715  jet->setHwPt(jetPtCorr);
716  }
717  else {
718  jet->setHwPt(0xFFFF);
719  }
720  }
721 
722  } else {
723  if(params_->jetCalibrationType() != "None" && params_->jetCalibrationType() != "none")
724  edm::LogError("l1t|stage 2") << "Invalid calibration type in calo params. Not calibrating Stage 2 Jets" << std::endl;
725  return;
726  }
727 
728 
729 
730 }
731 
732 //Function for the calibration, correct as a function of pT in bins of eta
734 
735  double logX = log10(pt);
736 
737  double term1 = par[1] / ( logX * logX + par[2] );
738  double term2 = par[3] * exp( -par[4]*((logX - par[5])*(logX - par[5])) );
739 
740  // Final fitting function, with sanity check
741  double f = par[0] + term1 + term2;
742  if (f < 0)
743  f = 0;
744  if (f != f) // stop NaN
745  f = 1;
746  return f;
747 }
748 
749 //NEW Function for the calibration, correct as a function of pT in bins of eta
751 
752  double f = par[0]+par[1]*TMath::Erf(par[2]*(log10(pt)-par[3])+par[4]*exp(par[5]*(log10(pt)-par[6])*(log10(pt)-par[6])));
753  // sanity check
754  if (f < 0)
755  f = 0;
756  if (f != f) // stop NaN
757  f = 1;
758  return f;
759 }
size
Write out results.
l1t::LUT * jetCalibrationLUT()
static int mpEta(int ieta)
Definition: CaloTools.cc:200
void setRawEt(short int et)
int chunkyDonutPUEstimate(Jet &jet, int pos, const std::vector< l1t::CaloTower > &towers)
std::string jetPUSType() const
static int regionEta(int ieta)
Definition: CaloTools.cc:220
void Merge(const std::vector< T > &aInput, std::vector< T > &aOutput)
std::map< int, int > getSumEtEtaMap(const std::vector< l1t::CaloTower > &towers)
std::string const & jetCalibrationType() const
unsigned jetBypassPUS() const
int hwPhi() const
Definition: L1Candidate.h:50
delete x;
Definition: CaloConfig.h:22
const double EtaMax[kNumberCalorimeter]
static const int kSatHcal
Definition: CaloTools.h:53
void calibrate(std::vector< Jet > &jets, int calibThreshold, bool isAllJets)
void setSeedEt(short int et)
PtEtaPhiMLorentzVectorD PtEtaPhiMLorentzVector
Lorentz vector with cartesian internal representation.
Definition: LorentzVector.h:25
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
Definition: Jet.h:21
std::vector< double > const & jetCalibrationParams() const
unsigned jetPUSUsePhiRing() const
int donutPUEstimate(int jetEta, int jetPhi, int size, const std::vector< l1t::CaloTower > &towers)
double p4[4]
Definition: TauolaWrapper.h:92
double jetLsb() const
vector< PseudoJet > jets
static const int kHFEnd
Definition: CaloTools.h:41
static const int kSatTower
Definition: CaloTools.h:55
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double f[11][100]
const double EtaMin[kNumberCalorimeter]
int hwEta() const
Definition: L1Candidate.h:49
int chunkySandwichPUEstimate(Jet &jet, int pos, const std::vector< l1t::CaloTower > &towers, const std::string chunkySandwichStr)
l1t::LUT * jetCompressPtLUT()
int hwPt() const
Definition: L1Candidate.h:48
static const int kHBHENrPhi
Definition: CaloTools.h:44
void setPUDonutEt(unsigned int i, short int et)
static const int kSatJet
Definition: CaloTools.h:58
double b
Definition: hdecay.h:120
void processEvent(const std::vector< CaloTower > &towers, std::vector< Jet > &jets, std::vector< Jet > &alljets) override
void setTowerIEta(short int ieta)
void create(const std::vector< CaloTower > &towers, std::vector< Jet > &jets, std::vector< Jet > &alljets, std::string PUSubMethod)
static int caloEta(int ietaMP)
Definition: CaloTools.cc:210
static const l1t::CaloTower & getTower(const std::vector< l1t::CaloTower > &towers, int iEta, int iPhi)
Definition: CaloTools.cc:37
double towerLsbSum() const
void setHwPt(int pt)
Definition: L1Candidate.h:41
double a
Definition: hdecay.h:121
void setTowerIPhi(short int iphi)
static int bin16Eta(int ieta)
Definition: CaloTools.cc:246
int data(unsigned int address) const
Definition: LUT.h:46
double jetSeedThreshold() const
void setPUEt(short int et)
bool operator>(const l1t::Jet &a, l1t::Jet &b)
l1t::LUT * jetCompressEtaLUT()
static const int kSatEcal
Definition: CaloTools.h:54
std::vector< int > getChunkyRing(Jet &jet, int pos, const std::vector< l1t::CaloTower > &towers, const std::string chunkyString)
#define constexpr