L1TCalorimeter/src/PUSubtractionMethods.cc

0001 // PUSubtractionMethods.cc
0002 // Authors: Alex Barbieri
0003 //          Kalanand Mishra, Fermilab
0004 //          Inga Bucinskaite, UIC
0005 //
0006 // This file should contain the different algorithms used to perform PU, UE subtraction.
0007
0008 //#include "DataFormats/L1TCalorimeter/interface/CaloRegion.h"
0009 #include "L1Trigger/L1TCalorimeter/interface/PUSubtractionMethods.h"
0010 #include "TMath.h"
0011
0012 //#include "DataFormats/L1CaloTrigger/interface/L1CaloRegionDetId.h"
0013 #include <vector>
0014
0015 namespace l1t {
0016
0017   /// --------------- For heavy ion -------------------------------------
0018   void HICaloRingSubtraction(const std::vector<l1t::CaloRegion> &regions,
0019                              std::vector<l1t::CaloRegion> *subRegions,
0020                              CaloParamsHelper const *params) {
0021     int puLevelHI[L1CaloRegionDetId::N_ETA];
0022
0023     for (unsigned i = 0; i < L1CaloRegionDetId::N_ETA; ++i) {
0024       puLevelHI[i] = 0;
0025     }
0026
0027     for (std::vector<CaloRegion>::const_iterator region = regions.begin(); region != regions.end(); region++) {
0028       puLevelHI[region->hwEta()] += region->hwPt();
0029     }
0030
0031     for (unsigned i = 0; i < L1CaloRegionDetId::N_ETA; ++i) {
0032       //puLevelHI[i] = floor(((double)puLevelHI[i] / (double)L1CaloRegionDetId::N_PHI)+0.5);
0033       puLevelHI[i] = (puLevelHI[i] + 9) * 455 / (1 << 13);  // approx equals X/18 +0.5
0034     }
0035
0036     for (std::vector<CaloRegion>::const_iterator region = regions.begin(); region != regions.end(); region++) {
0037       int subPt = std::max(0, region->hwPt() - puLevelHI[region->hwEta()]);
0038       int subEta = region->hwEta();
0039       int subPhi = region->hwPhi();
0040
0041       ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > ldummy(0, 0, 0, 0);
0042
0043       CaloRegion newSubRegion(
0044           *&ldummy, 0, 0, subPt, subEta, subPhi, region->hwQual(), region->hwEtEm(), region->hwEtHad());
0045       subRegions->push_back(newSubRegion);
0046     }
0047   }
0048
0049   void simpleHWSubtraction(const std::vector<l1t::CaloRegion> &regions, std::vector<l1t::CaloRegion> *subRegions) {
0050     for (std::vector<CaloRegion>::const_iterator region = regions.begin(); region != regions.end(); region++) {
0051       int subEta = region->hwEta();
0052       int subPhi = region->hwPhi();
0053       int subPt = region->hwPt();
0054
0055       if (subPt != (2 << 10) - 1)
0056         subPt = subPt - (10 + subEta);  // arbitrary value chosen in meeting
0057       if (subPt < 0)
0058         subPt = 0;
0059       ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > ldummy(0, 0, 0, 0);
0060
0061       CaloRegion newSubRegion(
0062           *&ldummy, 0, 0, subPt, subEta, subPhi, region->hwQual(), region->hwEtEm(), region->hwEtHad());
0063       subRegions->push_back(newSubRegion);
0064     }
0065   }
0066
0067   /// --------- New region correction (PUsub, no response correction at the moment) -----------
0068
0069   void RegionCorrection(const std::vector<l1t::CaloRegion> &regions,
0070                         std::vector<l1t::CaloRegion> *subRegions,
0071                         CaloParamsHelper const *params) {
0072     std::string regionPUSType = params->regionPUSType();
0073
0074     if (regionPUSType == "None") {
0075       for (std::vector<CaloRegion>::const_iterator notCorrectedRegion = regions.begin();
0076            notCorrectedRegion != regions.end();
0077            notCorrectedRegion++) {
0078         const CaloRegion &newSubRegion = *notCorrectedRegion;
0079         subRegions->push_back(newSubRegion);
0080       }
0081     }
0082
0083     if (regionPUSType == "HICaloRingSub") {
0084       HICaloRingSubtraction(regions, subRegions, params);
0085     }
0086
0087     if (regionPUSType == "PUM0") {
0088       int puMult = 0;
0089
0090       // ------------ This calulates PUM0 ------------------
0091       for (std::vector<CaloRegion>::const_iterator notCorrectedRegion = regions.begin();
0092            notCorrectedRegion != regions.end();
0093            notCorrectedRegion++) {
0094         int regionET = notCorrectedRegion->hwPt();
0095         if (regionET > 0) {
0096           puMult++;
0097         }
0098       }
0099       int pumbin = (int)puMult / 22;
0100       if (pumbin == 18)
0101         pumbin = 17;  // if puMult = 396 exactly there is an overflow
0102
0103       for (std::vector<CaloRegion>::const_iterator notCorrectedRegion = regions.begin();
0104            notCorrectedRegion != regions.end();
0105            notCorrectedRegion++) {
0106         int regionET = notCorrectedRegion->hwPt();
0107         int regionEta = notCorrectedRegion->hwEta();
0108         int regionPhi = notCorrectedRegion->hwPhi();
0109
0110         //int puSub = ceil(regionPUSParams[18*regionEta+pumbin]*2);
0111         int puSub = params->regionPUSValue(pumbin, regionEta);
0112         // The values in regionSubtraction are MULTIPLIED by
0113         // RegionLSB=.5 (physicalRegionEt), so to get back unmultiplied
0114         // regionSubtraction we want to multiply the number by 2
0115         // (aka divide by LSB).
0116
0117         int regionEtCorr = std::max(0, regionET - puSub);
0118         if (regionET == 1023)
0119           regionEtCorr = 1023;  // do not subtract overflow regions
0120         if ((regionET == 255) && (regionEta < 4 || regionEta > 17))
0121           regionEtCorr = 255;
0122
0123         ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > lorentz(0, 0, 0, 0);
0124         CaloRegion newSubRegion(*&lorentz,
0125                                 0,
0126                                 0,
0127                                 regionEtCorr,
0128                                 regionEta,
0129                                 regionPhi,
0130                                 notCorrectedRegion->hwQual(),
0131                                 notCorrectedRegion->hwEtEm(),
0132                                 notCorrectedRegion->hwEtHad());
0133         subRegions->push_back(newSubRegion);
0134       }
0135     }
0136   }
0137 }  // namespace l1t