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 },
52 std::vector<l1t::Jet> &
jets,
53 std::vector<l1t::Jet> & alljets) {
68 std::vector<l1t::Jet> &
jets,
69 std::vector<l1t::Jet> & alljets,
74 for (
int etaSide=1; etaSide>=-1; etaSide-=2) {
77 std::vector<int> ringGroup1, ringGroup2, ringGroup3, ringGroup4;
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 );
84 std::vector< std::vector<int> > theRings = { ringGroup1, ringGroup2, ringGroup3, ringGroup4 };
87 std::vector<l1t::Jet> jetsHalf;
90 for (
unsigned ringGroupIt=1; ringGroupIt<=theRings.size(); ringGroupIt++ ) {
93 std::vector<l1t::Jet> jetsAccu;
96 for (
unsigned ringIt=0; ringIt<theRings.at(ringGroupIt-1).size(); ringIt++ ) {
98 int ieta = theRings.at(ringGroupIt-1).at(ringIt);
101 std::vector<l1t::Jet> jetsRing;
107 if (jetsRing.size()==18)
break;
112 int seedEt = tow.
hwPt();
114 bool satSeed =
false;
115 bool vetoCandidate =
false;
126 for(
int deta = -4; deta < 5; ++deta ) {
127 for(
int dphi = -4; dphi < 5; ++dphi ) {
130 int ietaTest = ieta+deta;
131 int iphiTest = iphi+dphi;
138 if (ieta > 0 && ietaTest <=0) ietaTest -= 1;
139 if (ieta < 0 && ietaTest >=0) ietaTest += 1;
143 towEt = towTest.
hwPt();
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);
149 if (vetoCandidate)
break;
153 if(vetoCandidate)
break;
157 if (!vetoCandidate) {
171 PUSubMethod =
"PhiRing1";
173 PUSubMethod =
"ChunkyDonut";
175 if (PUSubMethod ==
"Donut") {
180 if (PUSubMethod ==
"ChunkyDonut"){
185 if(PUSubMethod.find(
"ChunkySandwich") != std::string::npos){
191 if(PUSubMethod ==
"PhiRing1"){
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;
201 for(
int ieta=jetEtaLow; ieta<=jetEtaHigh; ++ieta){
202 if(ieta==0)
continue;
203 auto iter=SumEtEtaMap.find(ieta);
210 if(PUSubMethod ==
"PhiRing2"){
215 int jetEtaLow=jetEta-size+1;
216 int jetEtaHigh=jetEta+size-1;
218 if(jetEta>0 && jetEtaLow<=0) jetEtaLow-=1;
219 if(jetEta<0 && jetEtaHigh>=0) jetEtaHigh+=1;
222 for(
int ieta=jetEtaLow; ieta<=jetEtaHigh; ++ieta){
223 if(ieta==0)
continue;
224 auto iter=SumEtEtaMap.find(ieta);
236 if (iEt<=0)
continue;
249 jetsRing.push_back(jet);
250 alljets.push_back(jet);
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);
266 jets.insert(jets.end(),jetsRing.begin(),jetsRing.end());
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));
282 for (
unsigned int iJet = 0; iJet < jets.size(); iJet++)
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);
290 std::vector<l1t::Jet> accumEtaPos;
291 std::vector<l1t::Jet> accumEtaNeg;
293 for(
int ieta = 0 ; ieta < 41 ; ++ieta)
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 );
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 );
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);
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);
321 accumEtaPos.resize(6);
322 accumEtaNeg.resize(6);
328 if (accjet.hwPt() > 0) jets.push_back(accjet);
332 if (accjet.hwPt() > 0) jets.push_back(accjet);
346 const std::vector<l1t::CaloTower> & towers){
349 std::vector<int>
ring(4,0);
351 int iphiUp = jetPhi +
size;
353 int iphiDown = jetPhi -
size;
356 int ietaUp = jetEta+
size;
357 int ietaDown = jetEta-
size;
359 for (
int ieta = jetEta - size+1; ieta < jetEta +
size; ++ieta)
365 if (jetEta > 0 && ieta <=0){
367 }
else if (jetEta < 0 && ieta >=0){
374 int towEt = tow.
hwPt();
383 for (
int iphi = jetPhi - size+1; iphi < jetPhi +
size; ++iphi)
391 int towEt = tow.
hwPt();
400 std::sort(ring.begin(), ring.end(), std::greater<int>());
402 return 4*( ring[1]+ring[2] );
408 std::map<int,int> SumEtEtaMap;
414 for(
int iEta=EtaMin ; iEta<=
EtaMax; iEta++ )
417 for(
int iPhi=phiMin; iPhi<=
phiMax; iPhi++){
419 int towEt=tow.
hwPt();
423 SumEtEtaMap[iEta]=SumEt;
433 const std::vector<l1t::CaloTower> & towers,
441 if(chunkyString ==
"ChunkySandwich8") ringSize = 8;
443 std::vector<int>
ring(ringSize,0);
449 for (
int stripIt=0; stripIt<nStrips; stripIt++) {
451 int iphiUp = jetPhi + size + stripIt;
452 int iphiDown = jetPhi - size - stripIt;
458 for (
int ieta=jetEta-size+1; ieta<jetEta+
size; ++ieta) {
463 if (jetEta>0 && towEta<=0) towEta-=1;
464 if (jetEta<0 && towEta>=0) towEta+=1;
467 int towEt = towPhiUp.
hwPt();
471 towEt = towPhiDown.
hwPt();
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;
484 for (
int iphi=jetPhi-size+1; iphi<jetPhi+
size; ++iphi) {
492 int towEt = towEtaUp.
hwPt();
498 for (
int iphi=jetPhi-size+1; iphi<jetPhi+
size; ++iphi) {
506 int towEt = towEtaDown.
hwPt();
511 else if(chunkyString ==
"ChunkySandwich2"){
512 for(
int i = 0;
i < 2; ++
i){ring[
i+2] = ring[
i];}
514 else if(chunkyString ==
"ChunkySandwich2Ave" || chunkyString ==
"ChunkySandwich"){
515 ring[2] = ((ring[0] + ring[1]) >> 1);
516 ring[3] = (ring[0] + ring[1]);
518 else if(chunkyString ==
"ChunkySandwich4" || chunkyString ==
"ChunkySandwich8"){
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));
526 for (
int ieta=jetEta-size+1; ieta<jetEta+
size; ++ieta) {
530 if (jetEta>0 && towEta<=0) towEta-=1;
531 if (jetEta<0 && towEta>=0) towEta+=1;
534 int towEt = towPhiUp.
hwPt();
535 ring[2 + rI*2] += towEt;
538 towEt = towPhiDown.
hwPt();
539 ring[3 + rI*2] += towEt;
553 const std::vector<l1t::CaloTower> & towers){
559 return ( ring[0] + ring[1] + ring[2] );
565 const std::vector<l1t::CaloTower> & towers,
570 return (ring[0] + ring[1] + ring[2]);
583 edm::LogError(
"l1t|stage 2") <<
"Invalid input vector to calo params. Input vector of size: " <<
589 for(std::vector<l1t::Jet>::iterator
jet = jets.begin();
jet!=jets.end();
jet++){
592 if(
jet->hwPt() < calibThreshold)
continue;
593 if(
jet->hwPt() >= 0xFFFF)
continue;
600 for(
int i=0;
i<6;
i++){
609 double correction =
calibFit(ptPhys, params);
611 jet->setHwPt(correction*
jet->hwPt());
620 edm::LogError(
"l1t|stage 2") <<
"Invalid input vector to calo params. Input vector of size: " <<
625 for(std::vector<l1t::Jet>::iterator
jet = jets.begin();
jet!=jets.end();
jet++){
627 if(
jet->hwPt() < calibThreshold)
continue;
628 if(
jet->hwPt() >= 0xFFFF)
continue;
633 for(
int i=0;
i<8;
i++){
638 double correction = params[6];
639 if (ptPhys>params[7]) correction =
calibFit(ptPhys, params);
641 jet->setHwPt(correction*
jet->hwPt());
650 edm::LogError(
"l1t|stage 2") <<
"Invalid input vector to calo params. Input vector of size: " <<
655 for(std::vector<l1t::Jet>::iterator
jet = jets.begin();
jet!=jets.end();
jet++){
657 if(
jet->hwPt() < calibThreshold)
continue;
658 if(
jet->hwPt() >= 0xFFFF)
continue;
663 for(
int i=0;
i<11;
i++){
668 double correction = params[7];
670 if (ptPhys<params[8]) correction = params[7];
671 else if (ptPhys>params[10]) correction = params[9];
674 jet->setHwPt(correction*
jet->hwPt());
689 for(std::vector<l1t::Jet>::iterator
jet = jets.begin();
jet!=jets.end();
jet++){
692 if(
jet->hwPt() < calibThreshold)
continue;
700 int jetHwPt =
jet->hwPt();
701 if (jetHwPt >= 0x200) {
706 unsigned int compBin = (etaBin<<4) | ptBin;
709 unsigned int multiplier = addPlusMult & 0x3ff;
711 int8_t addend = (addPlusMult>>10);
712 unsigned int jetPtCorr = ((
jet->hwPt()*multiplier)>>9) + addend;
714 if(jetPtCorr < 0xFFFF) {
715 jet->setHwPt(jetPtCorr);
718 jet->setHwPt(0xFFFF);
724 edm::LogError(
"l1t|stage 2") <<
"Invalid calibration type in calo params. Not calibrating Stage 2 Jets" << std::endl;
735 double logX = log10(pt);
737 double term1 = par[1] / ( logX * logX + par[2] );
738 double term2 = par[3] *
exp( -par[4]*((logX - par[5])*(logX - par[5])) );
741 double f = par[0] + term1 + term2;
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])));
l1t::LUT * jetCalibrationLUT()
void setRawEt(short int et)
int chunkyDonutPUEstimate(Jet &jet, int pos, const std::vector< l1t::CaloTower > &towers)
std::string jetPUSType() const
void Merge(const std::vector< T > &aInput, std::vector< T > &aOutput)
void accuSort(std::vector< Jet > &jets)
std::map< int, int > getSumEtEtaMap(const std::vector< l1t::CaloTower > &towers)
double calibFit(double, double *)
double calibFitErr(double, double *)
std::string const & jetCalibrationType() const
unsigned jetBypassPUS() const
const double EtaMax[kNumberCalorimeter]
CaloParamsHelper const *const params_
void calibrate(std::vector< Jet > &jets, int calibThreshold, bool isAllJets)
void setSeedEt(short int et)
PtEtaPhiMLorentzVectorD PtEtaPhiMLorentzVector
Lorentz vector with cartesian internal representation.
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
std::vector< double > const & jetCalibrationParams() const
unsigned jetPUSUsePhiRing() const
int donutPUEstimate(int jetEta, int jetPhi, int size, const std::vector< l1t::CaloTower > &towers)
Abs< T >::type abs(const T &t)
const double EtaMin[kNumberCalorimeter]
int chunkySandwichPUEstimate(Jet &jet, int pos, const std::vector< l1t::CaloTower > &towers, const std::string chunkySandwichStr)
l1t::LUT * jetCompressPtLUT()
Stage2Layer2JetAlgorithmFirmwareImp1(CaloParamsHelper const *params)
void setPUDonutEt(unsigned int i, short int et)
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)
double towerLsbSum() const
void setTowerIPhi(short int iphi)
int data(unsigned int address) const
double jetSeedThreshold() const
void setPUEt(short int et)
bool operator>(const l1t::Jet &a, l1t::Jet &b)
l1t::LUT * jetCompressEtaLUT()
std::vector< int > getChunkyRing(Jet &jet, int pos, const std::vector< l1t::CaloTower > &towers, const std::string chunkyString)