PUSubtractionMethods.cc

Go to the documentation of this file.

// PUSubtractionMethods.cc
// Authors: Alex Barbieri
//          Kalanand Mishra, Fermilab
//          Inga Bucinskaite, UIC
//
// This file should contain the different algorithms used to perform PU, UE subtraction.


//#include "DataFormats/L1TCalorimeter/interface/CaloRegion.h"
#include "L1Trigger/L1TCalorimeter/interface/PUSubtractionMethods.h"

//#include "DataFormats/L1CaloTrigger/interface/L1CaloRegionDetId.h"
#include <vector>

namespace l1t {

  void HICaloRingSubtraction(const std::vector<l1t::CaloRegion> & regions,
                 std::vector<l1t::CaloRegion> *subRegions)
  {
    int puLevelHI[L1CaloRegionDetId::N_ETA];
    double r_puLevelHI[L1CaloRegionDetId::N_ETA];
    int etaCount[L1CaloRegionDetId::N_ETA];
    for(unsigned i = 0; i < L1CaloRegionDetId::N_ETA; ++i)
    {
      puLevelHI[i] = 0;
      r_puLevelHI[i] = 0.0;
      etaCount[i] = 0;
    }

    for(std::vector<CaloRegion>::const_iterator region = regions.begin();
    region != regions.end(); region++){
      r_puLevelHI[region->hwEta()] += region->hwPt();
      etaCount[region->hwEta()]++;
    }

    for(unsigned i = 0; i < L1CaloRegionDetId::N_ETA; ++i)
    {
      puLevelHI[i] = floor(r_puLevelHI[i]/etaCount[i] + 0.5);
    }

    for(std::vector<CaloRegion>::const_iterator region = regions.begin(); region!= regions.end(); region++){
      int subPt = std::max(0, region->hwPt() - puLevelHI[region->hwEta()]);
      int subEta = region->hwEta();
      int subPhi = region->hwPhi();

      ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > ldummy(0,0,0,0);

      CaloRegion newSubRegion(*&ldummy, 0, 0, subPt, subEta, subPhi, region->hwQual(), region->hwEtEm(), region->hwEtHad());
      subRegions->push_back(newSubRegion);
    }
  }

  void simpleHWSubtraction(const std::vector<l1t::CaloRegion> & regions,
               std::vector<l1t::CaloRegion> *subRegions)
  {
    for(std::vector<CaloRegion>::const_iterator region = regions.begin();
    region != regions.end(); region++){
      int subEta = region->hwEta();
      int subPhi = region->hwPhi();
      int subPt = region->hwPt();

      //std::cout << "pre sub: " << subPt;
      if(subPt != (2<<10)-1)
    subPt = subPt - (10+subEta); // arbitrary value chosen in meeting
      if(subPt < 0)
    subPt = 0;
      //std::cout << " post sub: " << subPt << std::endl;
      ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > ldummy(0,0,0,0);

      CaloRegion newSubRegion(*&ldummy, 0, 0, subPt, subEta, subPhi, region->hwQual(), region->hwEtEm(), region->hwEtHad());
      subRegions->push_back(newSubRegion);
    }
  }



 void RegionCorrection(const std::vector<l1t::CaloRegion> & regions,
               std::vector<l1t::CaloRegion> *subRegions,
               std::vector<double> regionPUSParams,
               std::string regionPUSType)
{


  if(regionPUSType == "None") {
    for(std::vector<CaloRegion>::const_iterator notCorrectedRegion = regions.begin();
    notCorrectedRegion != regions.end(); notCorrectedRegion++){
      CaloRegion newSubRegion= *notCorrectedRegion;
      subRegions->push_back(newSubRegion);
    }
  }

  if (regionPUSType == "HICaloRingSub") {
    HICaloRingSubtraction(regions, subRegions);
  }

  if (regionPUSType == "PUM0") {
    int puMult = 0;

    // ------------ This calulates PUM0 ------------------
    for(std::vector<CaloRegion>::const_iterator notCorrectedRegion = regions.begin();
    notCorrectedRegion != regions.end(); notCorrectedRegion++){
      int regionET = notCorrectedRegion->hwPt();
      // cout << "regionET: " << regionET <<endl;
      if (regionET > 0) {puMult++;}
    }
    int pumbin = (int) puMult/22;
    if(pumbin == 18) pumbin = 17; // if puMult = 396 exactly there is an overflow

    for(std::vector<CaloRegion>::const_iterator notCorrectedRegion = regions.begin();
    notCorrectedRegion != regions.end(); notCorrectedRegion++){

      int regionET = notCorrectedRegion->hwPt();
      int regionEta = notCorrectedRegion->hwEta();
      int regionPhi = notCorrectedRegion->hwPhi();

      int puSub = ceil(regionPUSParams[18*regionEta+pumbin]*2);
      // The values in regionSubtraction are MULTIPLIED by
      // RegionLSB=.5 (physicalRegionEt), so to get back unmultiplied
      // regionSubtraction we want to multiply the number by 2
      // (aka divide by LSB).

      int regionEtCorr = std::max(0, regionET - puSub);

      ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > lorentz(0,0,0,0);
      CaloRegion newSubRegion(*&lorentz, 0, 0, regionEtCorr, regionEta, regionPhi, notCorrectedRegion->hwQual(), notCorrectedRegion->hwEtEm(), notCorrectedRegion->hwEtHad());
      subRegions->push_back(newSubRegion);
    }
  }

}

}