/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoMET/METAlgorithms/src/CaloSpecificAlgo.cc

Go to the documentation of this file.

#include "DataFormats/Math/interface/LorentzVector.h"
#include "RecoMET/METAlgorithms/interface/CaloSpecificAlgo.h"
#include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
#include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"

using namespace reco;
using namespace std;

//-------------------------------------------------------------------------
// This algorithm adds calorimeter specific global event information to 
// the MET object which may be useful/needed for MET Data Quality Monitoring
// and MET cleaning.  This list is not exhaustive and additional 
// information will be added in the future. 
//-------------------------------------

reco::CaloMET CaloSpecificAlgo::addInfo(edm::Handle<edm::View<Candidate> > towers, CommonMETData met, bool noHF, double globalThreshold)
{ 
  // Instantiate the container to hold the calorimeter specific information
  SpecificCaloMETData specific;
  // Initialise the container 
  specific.MaxEtInEmTowers = 0.0;         // Maximum energy in EM towers
  specific.MaxEtInHadTowers = 0.0;        // Maximum energy in HCAL towers
  specific.HadEtInHO = 0.0;          // Hadronic energy fraction in HO
  specific.HadEtInHB = 0.0;          // Hadronic energy in HB
  specific.HadEtInHF = 0.0;          // Hadronic energy in HF
  specific.HadEtInHE = 0.0;          // Hadronic energy in HE
  specific.EmEtInEB = 0.0;           // Em energy in EB
  specific.EmEtInEE = 0.0;           // Em energy in EE
  specific.EmEtInHF = 0.0;           // Em energy in HF
  specific.EtFractionHadronic = 0.0; // Hadronic energy fraction
  specific.EtFractionEm = 0.0;       // Em energy fraction
  specific.CaloSETInpHF = 0.0;        // CaloSET in HF+ 
  specific.CaloSETInmHF = 0.0;        // CaloSET in HF- 
  specific.CaloMETInpHF = 0.0;        // CaloMET in HF+ 
  specific.CaloMETInmHF = 0.0;        // CaloMET in HF- 
  specific.CaloMETPhiInpHF = 0.0;     // CaloMET-phi in HF+ 
  specific.CaloMETPhiInmHF = 0.0;     // CaloMET-phi in HF- 
  specific.METSignificance = 0.0;

  double totalEt = 0.0; 
  double totalEm     = 0.0;
  double totalHad    = 0.0;
  double MaxTowerEm  = 0.0;
  double MaxTowerHad = 0.0;
  double sumEtInpHF = 0.0;
  double sumEtInmHF = 0.0;
  double MExInpHF = 0.0;
  double MEyInpHF = 0.0;
  double MExInmHF = 0.0;
  double MEyInmHF = 0.0;

  if( towers->size() == 0 )  // if there are no towers, return specific = 0
    {
      //   LogDebug("CaloMET") << "Number of Candidate CaloTowers is zero : Unable to calculate calo specific info. " ;
      const LorentzVector p4( met.mex, met.mey, 0.0, met.met );
      const Point vtx( 0.0, 0.0, 0.0 );
      CaloMET specificmet( specific, met.sumet, p4, vtx );
      return specificmet;
    }

  edm::View<Candidate>::const_iterator towerCand = towers->begin();
  for( ; towerCand != towers->end(); towerCand++ ) 
    {
      const Candidate* candidate = &(*towerCand);
      if (candidate) {
        const CaloTower* calotower = dynamic_cast<const CaloTower*> (candidate);
        if (calotower)
          {
            double caloTowerEt=calotower->et();
            double caloTowerHadEt=calotower->hadEt();
            double caloTowerEmEt=calotower->emEt();
            if(caloTowerEt < globalThreshold) continue;
            totalEt  += caloTowerEt;
            totalEm  += caloTowerEmEt;
                   
            //totalHad += caloTowerHadEt + calotower->outerEt() ;
                   
            bool hadIsDone = false;
            bool emIsDone = false;
            int cell = calotower->constituentsSize();
            while ( --cell >= 0 && (!hadIsDone || !emIsDone) ) 
              {
                DetId id = calotower->constituent( cell );
                if( !hadIsDone && id.det() == DetId::Hcal ) 
                  {
                    HcalSubdetector subdet = HcalDetId(id).subdet();
                    if( subdet == HcalBarrel || subdet == HcalOuter )
                      {
                        if( caloTowerHadEt > MaxTowerHad ) MaxTowerHad = caloTowerHadEt;
                        specific.HadEtInHB   += caloTowerHadEt;
                        specific.HadEtInHO   += calotower->outerEt();
                      }
                    else if( subdet == HcalEndcap )
                      {
                        if( caloTowerHadEt > MaxTowerHad ) MaxTowerHad = caloTowerHadEt;
                        specific.HadEtInHE   += caloTowerHadEt;
                      }
                    else if( subdet == HcalForward )
                      {
                        if (!noHF)
                          {
                            if( caloTowerHadEt > MaxTowerHad ) MaxTowerHad = caloTowerHadEt;
                            if( caloTowerEmEt  > MaxTowerEm  ) MaxTowerEm  = caloTowerEmEt;
                            //These quantities should be nonzero only if HF is included, i.e., noHF == false
                            specific.HadEtInHF   += caloTowerHadEt;
                            specific.EmEtInHF    += caloTowerEmEt;
                          }
                        else
                          {
                            //These quantities need to be corrected from above if HF is excluded
                            // totalHad             -= caloTowerHadEt;  
                            totalEm              -= caloTowerEmEt;
                            totalEt              -= caloTowerEt;
                          }
                        // These get calculate regardless of NoHF == true or not.
                        // They are needed below for either case. 
                        if (calotower->eta()>=0)
                          {
                            sumEtInpHF  += caloTowerEt;
                            MExInpHF    -= (caloTowerEt * cos(calotower->phi()));
                            MEyInpHF    -= (caloTowerEt * sin(calotower->phi()));
                          }
                        else
                          {
                            sumEtInmHF  += caloTowerEt;
                            MExInmHF    -= (caloTowerEt * cos(calotower->phi()));
                            MEyInmHF    -= (caloTowerEt * sin(calotower->phi()));
                          }
                      }
                    hadIsDone = true;
                  }
                else if( !emIsDone && id.det() == DetId::Ecal )
                  {
                    EcalSubdetector subdet = EcalSubdetector( id.subdetId() );
                    if( caloTowerEmEt  > MaxTowerEm  ) MaxTowerEm  = caloTowerEmEt;
                    if( subdet == EcalBarrel )
                      {
                        specific.EmEtInEB    += caloTowerEmEt; 
                      }
                    else if( subdet == EcalEndcap ) 
                      {
                        specific.EmEtInEE    += caloTowerEmEt;
                      }
                    emIsDone = true;
                  }
              }
          }
      }
    }
  
  //Following Greg L's suggestion to calculate this quantity outside of the loop and to avoid confusion. 
  //This should work regardless of HO's inclusion / exclusion .
  totalHad += (totalEt - totalEm);
  
  if(!noHF)
    { // Form sub-det specific MET-vectors
      LorentzVector METpHF(MExInpHF, MEyInpHF, 0, sqrt(MExInpHF*MExInpHF + MEyInpHF*MEyInpHF));
      LorentzVector METmHF(MExInmHF, MEyInmHF, 0, sqrt(MExInmHF*MExInmHF + MEyInmHF*MEyInmHF));
      specific.CaloMETInpHF = METpHF.pt();
      specific.CaloMETInmHF = METmHF.pt();
      specific.CaloMETPhiInpHF = METpHF.Phi();
      specific.CaloMETPhiInmHF = METmHF.Phi();
      specific.CaloSETInpHF = sumEtInpHF;
      specific.CaloSETInmHF = sumEtInmHF;
    }
  else
    { // remove HF from MET calculation 
      met.mex   -= (MExInmHF + MExInpHF);
      met.mey   -= (MEyInmHF + MEyInpHF);
      met.sumet -= (sumEtInpHF + sumEtInmHF);
      met.met    = sqrt(met.mex*met.mex + met.mey*met.mey);   
    } 

  specific.MaxEtInEmTowers         = MaxTowerEm;  
  specific.MaxEtInHadTowers        = MaxTowerHad;         
  specific.EtFractionHadronic = totalHad / totalEt; 
  specific.EtFractionEm       =  totalEm / totalEt;       

  const LorentzVector p4( met.mex, met.mey, 0.0, met.met );
  const Point vtx( 0.0, 0.0, 0.0 );
  // Create and return an object of type CaloMET, which is a MET object with 
  // the extra calorimeter specfic information added
  CaloMET specificmet( specific, met.sumet, p4, vtx );
  return specificmet;
}
//-------------------------------------------------------------------------