CMS 3D CMS Logo

Stage1Layer2TauAlgorithmImpHW.cc
Go to the documentation of this file.
1 
11 
19 
20 
21 using namespace std;
22 using namespace l1t;
23 
24 
25 Stage1Layer2TauAlgorithmImpHW::Stage1Layer2TauAlgorithmImpHW(CaloParamsHelper const* params) : params_(params),
26  isoTauLut{ make_unique<Stage1TauIsolationLUT>(params_)}
27 {
28 
29 }
30 
31 
32 
33 
34 
35 void l1t::Stage1Layer2TauAlgorithmImpHW::processEvent(const std::vector<l1t::CaloEmCand> & EMCands,
36  const std::vector<l1t::CaloRegion> & regions,
37  std::vector<l1t::Tau> * isoTaus,
38  std::vector<l1t::Tau> * taus) {
39 
40  double towerLsb = params_->towerLsbSum();
41 
42  int tauSeedThreshold= floor( params_->tauSeedThreshold()/towerLsb + 0.5); // convert GeV to HW units
43  int tauNeighbourThreshold= floor( params_->tauNeighbourThreshold()/towerLsb + 0.5); // convert GeV to HW units
44  int tauMaxPtTauVeto = floor( params_->tauMaxPtTauVeto()/towerLsb + 0.5);
45  int isoTauEtaMin = params_->isoTauEtaMin();
46  int isoTauEtaMax = params_->isoTauEtaMax();
47 
48  std::vector<l1t::CaloRegion> subRegions;
49 
50 
51  //Region Correction will return uncorrected subregions if
52  //regionPUSType is set to None in the config
53  RegionCorrection(regions, &subRegions, params_);
54 
55 
56 
57  // ----- need to cluster jets in order to compute jet isolation ----
58  std::vector<l1t::Jet> unCorrJets;
59  TwelveByTwelveFinder(0, &subRegions, &unCorrJets);
60 
61  std::vector<l1t::Tau> preGtTaus;
62  std::vector<l1t::Tau> preSortTaus;
63  std::vector<l1t::Tau> sortedTaus;
64  std::vector<l1t::Tau> preGtIsoTaus;
65  std::vector<l1t::Tau> preSortIsoTaus;
66  std::vector<l1t::Tau> sortedIsoTaus;
67 
68  for(CaloRegionBxCollection::const_iterator region = subRegions.begin();
69  region != subRegions.end(); region++) {
70  if(region->hwEta() < 4 || region->hwEta() > 17)
71  continue;
72 
73  int regionEt = region->hwPt();
74  if(regionEt < tauSeedThreshold) continue;
75 
76  int regionEta = region->hwEta();
77  int regionPhi = region->hwPhi();
78 
79  int tauEt = regionEt;
80  int isoFlag = 0; // is 1 if it passes the relative jet iso requirement
81  int quality = 0; //doesn't really mean anything and isn't used
82 
83  int highestNeighborEt=-1;
84  int highestNeighborEta=-1;
85  int highestNeighborPhi=-1;
86  int highestNeighborTauVeto=-1;
87 
88  //Find neighbor with highest Et
89  for(CaloRegionBxCollection::const_iterator neighbor = subRegions.begin();
90  neighbor != subRegions.end(); neighbor++) {
91 
92  int neighborPhi = neighbor->hwPhi();
93  int neighborEta = neighbor->hwEta();
94  if(neighborEta < 4 || neighborEta > 17)
95  continue;
96 
97  int deltaPhi = regionPhi - neighborPhi;
98  if (std::abs(deltaPhi) == L1CaloRegionDetId::N_PHI-1)
99  deltaPhi = -deltaPhi/std::abs(deltaPhi); //18 regions in phi
100 
101  deltaPhi = std::abs(deltaPhi);
102  int deltaEta = std::abs(regionEta - neighborEta);
103 
104  if (deltaPhi + deltaEta > 0 && deltaPhi + deltaEta < 2) { //nondiagonal neighbors
105  if (neighbor->hwPt() > highestNeighborEt) {
106  highestNeighborEt = neighbor->hwPt();
107  highestNeighborEta = neighbor->hwEta();
108  highestNeighborPhi = neighbor->hwPhi();
109  highestNeighborTauVeto = neighbor->hwQual() & 0x1; // tauVeto should be the first bit of quality integer
110  }
111  }
112  }
113 
114  string NESW = findNESW(regionEta, regionPhi, highestNeighborEta, highestNeighborPhi);
115 
116  if((tauEt > highestNeighborEt && (NESW=="isEast" || NESW=="isNorth"))
117  || (tauEt >= highestNeighborEt && (NESW=="isSouth" || NESW=="isWest"))
118  || highestNeighborEt == 0 ) {
119 
120  if (highestNeighborEt >= tauNeighbourThreshold) tauEt += highestNeighborEt;
121 
122  int regionTauVeto = region->hwQual() & 0x1; // tauVeto should be the first bit of quality integer
123 
124  // compute relative jet isolation
125  if (region->hwEta() >= isoTauEtaMin && region->hwEta() <= isoTauEtaMax ){
126  if ((highestNeighborTauVeto == 0 && regionTauVeto == 0) || tauEt > tauMaxPtTauVeto) {
127  int jetEt=AssociatedJetPt(region->hwEta(), region->hwPhi(),&unCorrJets);
128  if (jetEt>0){
129  unsigned int MAX_LUT_ADDRESS = params_->tauIsolationLUT()->maxSize()-1;
130  // unsigned int lutAddress = isoLutIndex(tauEt,jetEt);
131  unsigned lutAddress = isoTauLut->lutAddress(tauEt,jetEt);
132  if (tauEt >0){
133  if (lutAddress > MAX_LUT_ADDRESS) lutAddress = MAX_LUT_ADDRESS;
134  isoFlag = params_->tauIsolationLUT()->data(lutAddress);
135  }
136 
137  }else{ // no associated jet
138  isoFlag=1;
139  }
140  }
141  }
142 
143  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > tauLorentz(0,0,0,0);
144 
145  l1t::Tau theTau(*&tauLorentz, tauEt, region->hwEta(), region->hwPhi(), quality, isoFlag);
146 
147  preGtTaus.push_back(theTau);
148  if(isoFlag)
149  preGtIsoTaus.push_back(theTau);
150  }
151  }
152 
153  // TauToGtPtScales(params_, preGtTaus, preRankTaus);
154  // TauToGtPtScales(params_, preGtIsoTaus, preRankIsoTaus);
155  calibrateAndRankTaus(params_,&preGtTaus,&preSortTaus);
156  calibrateAndRankTaus(params_,&preGtIsoTaus, &preSortIsoTaus);
157 
158  SortTaus(&preSortTaus, &sortedTaus);
159  SortTaus(&preSortIsoTaus, &sortedIsoTaus);
160 
161  TauToGtEtaScales(params_, &sortedTaus, taus);
162  TauToGtEtaScales(params_, &sortedIsoTaus, isoTaus);
163 
164  // set all filler taus to have isolation bit set
165  for(std::vector<l1t::Tau>::iterator iTau = isoTaus->begin(); iTau != isoTaus->end(); ++iTau)
166  iTau->setHwIso(1);
167 }
168 
169 
170 // Compute jet isolation.
172  const std::vector<l1t::Jet> & jets) const {
173 
174  for(JetBxCollection::const_iterator jet = jets.begin();
175  jet != jets.end(); jet++) {
176 
177  if (ieta==jet->hwEta() && iphi==jet->hwPhi()){
178 
179  double isolation = (double) (jet->hwPt() - et);
180  return isolation/et;
181  }
182  }
183 
184  // set output
185  return 999.;
186 }
187 
188 
189 // Find if the neighbor with the highest Et is N, E, S, or W
190 string l1t::Stage1Layer2TauAlgorithmImpHW::findNESW(int ieta, int iphi, int neta, int nphi) const {
191 
192  int deltaPhi = iphi - nphi;
193  if (std::abs(deltaPhi) == L1CaloRegionDetId::N_PHI-1)
194  deltaPhi = -deltaPhi/std::abs(deltaPhi); //18 regions in phi
195 
196  int deltaEta = ieta - neta;
197 
198  if ((std::abs(deltaPhi) + std::abs(deltaEta)) < 2) {
199  if (deltaEta==-1) {
200  return "isEast";
201  }
202  else if (deltaEta==0) {
203  if (deltaPhi==-1) {
204  return "isNorth";
205  }
206  if (deltaPhi==1) {
207  return "isSouth";
208  }
209  }
210  else {
211  return "isWest";
212  }
213  }
214 
215  return "999";
216 
217 }
218 
220  const std::vector<l1t::Jet> * jets) const {
221 
222  int pt = -1;
223 
224 
225  for(JetBxCollection::const_iterator itJet = jets->begin();
226  itJet != jets->end(); ++itJet){
227 
228  int jetEta = itJet->hwEta();
229  int jetPhi = itJet->hwPhi();
230  if ((jetEta == ieta) && (jetPhi == iphi)){
231  pt = itJet->hwPt();
232  break;
233  }
234  }
235 
236  // set output
237  return pt;
238 }
239 
240 unsigned l1t::Stage1Layer2TauAlgorithmImpHW::isoLutIndex(unsigned int tauPt,unsigned int jetPt) const
241 {
242  //const unsigned int nbitsTau=9; // number of bits used for et in LUT file (needed for left shift operation)
243  //const unsigned int nbitsJet=9;
244 
245  const unsigned int nbitsTau=8; // number of bits used for et in LUT file (needed for left shift operation)
246  const unsigned int nbitsJet=8;
247 
248  const unsigned int maxJet = pow(2,nbitsJet)-1;
249  const unsigned int maxTau = pow(2,nbitsTau)-1;
250 
251  if (nbitsTau < 9)
252  {
253  if (nbitsTau == 6)
254  {
255  tauPt=tauPt>>3;
256  }
257  else if (nbitsTau == 7)
258  {
259  tauPt=tauPt>>2;
260  }
261  else if (nbitsTau == 8)
262  {
263  tauPt=tauPt>>1;
264  }
265  }
266 
267  if (nbitsJet < 9)// no need to do shift if nbits>=9
268  {
269  if (nbitsJet == 6)
270  {
271  jetPt=jetPt>>3;
272  }
273  else if (nbitsJet == 7)
274  {
275  jetPt=jetPt>>2;
276  }
277  else if (nbitsJet == 8)
278  {
279  jetPt=jetPt>>1;
280  }
281  }
282 
283  if (jetPt>maxJet) jetPt=maxJet;
284  if (tauPt>maxTau) tauPt=maxTau;
285 
286  unsigned int address= (jetPt << nbitsTau) + tauPt;
287 
288  return address;
289 }
double JetIsolation(int et, int ieta, int iphi, const std::vector< l1t::Jet > &jets) const
unsigned int maxSize() const
Definition: LUT.h:53
double tauMaxPtTauVeto() const
void calibrateAndRankTaus(CaloParamsHelper const *params, const std::vector< l1t::Tau > *input, std::vector< l1t::Tau > *output)
Definition: Tau.h:21
double tauNeighbourThreshold() const
delete x;
Definition: CaloConfig.h:22
void TwelveByTwelveFinder(const int, const std::vector< l1t::CaloRegion > *regions, std::vector< l1t::Jet > *uncalibjets)
void SortTaus(std::vector< l1t::Tau > *input, std::vector< l1t::Tau > *output)
static const double deltaEta
Definition: CaloConstants.h:8
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) --------— ...
string findNESW(int ieta, int iphi, int neta, int nphi) const
unsigned isoLutIndex(unsigned int tauPt, unsigned int jetPt) const
void processEvent(const std::vector< l1t::CaloEmCand > &EMCands, const std::vector< l1t::CaloRegion > &regions, std::vector< l1t::Tau > *isoTaus, std::vector< l1t::Tau > *taus) override
vector< PseudoJet > jets
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int AssociatedJetPt(int ieta, int iphi, const std::vector< l1t::Jet > *jets) const
void TauToGtEtaScales(CaloParamsHelper const *params, const std::vector< l1t::Tau > *input, std::vector< l1t::Tau > *output)
et
define resolution functions of each parameter
double tauSeedThreshold() const
double towerLsbSum() const
int data(unsigned int address) const
Definition: LUT.h:46
std::unique_ptr< Stage1TauIsolationLUT > isoTauLut
static const unsigned N_PHI
l1t::LUT * tauIsolationLUT()
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
std::vector< T >::const_iterator const_iterator
Definition: BXVector.h:20