EcalHaloAlgo.cc

Go to the documentation of this file.

#include "RecoMET/METAlgorithms/interface/EcalHaloAlgo.h"
#include "DataFormats/Common/interface/ValueMap.h"

/*
  [class]:  EcalHaloAlgo
  [authors]: R. Remington, The University of Florida
  [description]: See EcalHaloAlgo.h
  [date]: October 15, 2009
*/
namespace {
  constexpr float c_cm_per_ns  = 29.9792458;
};
using namespace std;
using namespace reco;
using namespace edm;

bool CompareTime(const EcalRecHit* x, const EcalRecHit* y){ return x->time() < y->time(); }

EcalHaloAlgo::EcalHaloAlgo()
{
  RoundnessCut = 0 ;
  AngleCut = 0 ;
  EBRecHitEnergyThreshold = 0.;
  EERecHitEnergyThreshold = 0.;
  ESRecHitEnergyThreshold = 0.;
  SumEnergyThreshold = 0.;
  NHitsThreshold =0;

  geo = nullptr;
}

EcalHaloData EcalHaloAlgo::Calculate(const CaloGeometry& TheCaloGeometry, edm::Handle<reco::PhotonCollection>& ThePhotons, edm::Handle<reco::SuperClusterCollection>& TheSuperClusters, edm::Handle<EBRecHitCollection>& TheEBRecHits, edm::Handle<EERecHitCollection>& TheEERecHits, edm::Handle<ESRecHitCollection>& TheESRecHits, edm::Handle<HBHERecHitCollection>& TheHBHERecHits, const edm::EventSetup& TheSetup)
{  

  EcalHaloData TheEcalHaloData;

  // Store energy sum of rechits as a function of iPhi (iphi goes from 1 to 72)       
  float SumE[361];           
  // Store number of rechits as a function of iPhi
  int NumHits[361];               
  // Store minimum time of rechit as a function of iPhi
  float MinTimeHits[361];                             
  // Store maximum time of rechit as a function of iPhi 
  float MaxTimeHits[361];                      

  // initialize
  for(int i = 0 ; i < 361 ; i++ )
    {
      SumE[i] = 0.;
      NumHits[i] = 0 ;
      MinTimeHits[i] = 9999.;
      MaxTimeHits[i] = -9999.;
    }

  // Loop over EB RecHits
  for(EBRecHitCollection::const_iterator hit = TheEBRecHits->begin() ; hit != TheEBRecHits->end() ; hit++ )
    {
      // Arbitrary threshold to kill noise (needs to be optimized with data)
      if (hit->energy() < EBRecHitEnergyThreshold ) continue;
      
      // Get Det Id of the rechit
      DetId id = DetId(hit->id()); 

      // Get EB geometry 
      const CaloSubdetectorGeometry* TheSubGeometry = TheCaloGeometry.getSubdetectorGeometry(DetId::Ecal, 1);
      EBDetId EcalID(id.rawId());
      auto cell = (TheSubGeometry) ? (TheSubGeometry->getGeometry(id)) : nullptr;
  
      if(cell)
    {
      // GlobalPoint globalpos = cell->getPosition();
      //      float r = TMath::Sqrt ( globalpos.y()*globalpos.y() + globalpos.x()*globalpos.x());
      int iPhi = EcalID.iphi();

      if( iPhi < 361 )    // just to be safe
        {
          //iPhi = (iPhi-1)/5 +1;  // convert ecal iphi to phiwedge iphi  (e.g. there are 5 crystal per phi wedge, as in calotowers )
          SumE[iPhi] += hit->energy();
          NumHits[iPhi] ++;

          float time = hit->time();
          MinTimeHits[iPhi] = time < MinTimeHits[iPhi] ? time : MinTimeHits[iPhi];
          MaxTimeHits[iPhi] = time > MaxTimeHits[iPhi] ? time : MaxTimeHits[iPhi];
        }
    }
    }
  
  //for( int iPhi = 1 ; iPhi < 73; iPhi++ )
  for( int iPhi = 1 ; iPhi < 361; iPhi++ )
    {
      if( SumE[iPhi] >= SumEnergyThreshold && NumHits[iPhi] > NHitsThreshold )
    {
      // Build PhiWedge and store to EcalHaloData if energy or #hits pass thresholds
      PhiWedge wedge(SumE[iPhi], iPhi, NumHits[iPhi], MinTimeHits[iPhi], MaxTimeHits[iPhi]);
      
      // Loop over rechits again to calculate direction based on timing info
      
      // Loop over EB RecHits
      std::vector<const EcalRecHit*>  Hits;
      for(EBRecHitCollection::const_iterator hit = TheEBRecHits->begin() ; hit != TheEBRecHits->end() ; hit++ )
        {
          if (hit->energy() < EBRecHitEnergyThreshold ) continue;
          
          // Get Det Id of the rechit
          DetId id = DetId(hit->id()); 
          EBDetId EcalID(id.rawId());
          int Hit_iPhi = EcalID.iphi();
          //Hit_iPhi = (Hit_iPhi-1)/5 +1; // convert ecal iphi to phiwedge iphi
          if( Hit_iPhi != iPhi ) continue;
          Hits.push_back( &(*hit) );
          
        }
      std::sort( Hits.begin() , Hits.end(), CompareTime);
      float MinusToPlus = 0.;
      float PlusToMinus = 0.;
      for( unsigned int i = 0 ; i < Hits.size() ; i++ )
        {
          DetId id_i = DetId(Hits[i]->id()); 
          EBDetId EcalID_i(id_i.rawId());
          int ieta_i = EcalID_i.ieta();
          for( unsigned int j = (i+1) ; j < Hits.size() ; j++ )
        {
          DetId id_j = DetId(Hits[j]->id() );
          EBDetId EcalID_j(id_j.rawId());
          int ieta_j = EcalID_j.ieta();
          if( ieta_i > ieta_j ) PlusToMinus += TMath::Abs(ieta_j - ieta_i );
          else MinusToPlus += TMath::Abs(ieta_j -ieta_i) ;
        }
        }
      
      float PlusZOriginConfidence = (PlusToMinus+MinusToPlus) ? PlusToMinus / (PlusToMinus+MinusToPlus) : -1.;
      wedge.SetPlusZOriginConfidence(PlusZOriginConfidence);
      TheEcalHaloData.GetPhiWedges().push_back(wedge);
    }
    }

  std::vector<float> vShowerShapes_Roundness;
  std::vector<float> vShowerShapes_Angle ;
  if(TheSuperClusters.isValid()){
  for(reco::SuperClusterCollection::const_iterator cluster = TheSuperClusters->begin() ; cluster != TheSuperClusters->end() ; cluster++ )
    {
      if( abs(cluster->eta()) <= 1.48  )
    { 
      vector<float> shapes = EcalClusterTools::roundnessBarrelSuperClusters( *cluster, (*TheEBRecHits.product()));
      float roundness = shapes[0];
      float angle = shapes[1];
      
      // Check if supercluster belongs to photon and passes the cuts on roundness and angle, if so store the reference to it
      if( (roundness >=0 && roundness < GetRoundnessCut()) &&  angle >= 0 && angle < GetAngleCut() )
        {
          edm::Ref<SuperClusterCollection> TheClusterRef( TheSuperClusters, cluster - TheSuperClusters->begin());
          bool BelongsToPhoton = false;
          if( ThePhotons.isValid() )
        {
          for(reco::PhotonCollection::const_iterator iPhoton = ThePhotons->begin() ; iPhoton != ThePhotons->end() ; iPhoton++ )
            {
              if(iPhoton->isEB())
            if ( TheClusterRef ==  iPhoton->superCluster() )
              {
                BelongsToPhoton = true;
                break;
              }
            }
        }
          //Only store refs to suspicious EB SuperClusters which belong to Photons
          //Showershape variables are more discriminating for these cases  
          if( BelongsToPhoton )
        { 
          TheEcalHaloData.GetSuperClusters().push_back( TheClusterRef ) ; 
        }
        }
      vShowerShapes_Roundness.push_back(shapes[0]);
      vShowerShapes_Angle.push_back(shapes[1]);
    }
      else
    { 
      vShowerShapes_Roundness.push_back(-1.);
      vShowerShapes_Angle.push_back(-1.);
    }
    }
  
  
  edm::ValueMap<float>::Filler TheRoundnessFiller( TheEcalHaloData.GetShowerShapesRoundness() );
  TheRoundnessFiller.insert( TheSuperClusters, vShowerShapes_Roundness.begin(), vShowerShapes_Roundness.end() );
  TheRoundnessFiller.fill();  

  edm::ValueMap<float>::Filler TheAngleFiller( TheEcalHaloData.GetShowerShapesAngle() );
  TheAngleFiller.insert( TheSuperClusters, vShowerShapes_Angle.begin() , vShowerShapes_Angle.end() );
  TheAngleFiller.fill();
  }

  edm::ESHandle<CaloGeometry> pGeo;
  TheSetup.get<CaloGeometryRecord>().get(pGeo);
  geo = pGeo.product();
  
  //Halo cluster building:
  //Various clusters are built, depending on the subdetector.
  //In barrel, one looks for deposits narrow in phi.
  //In endcaps, one looks for localized deposits (dr condition in EE where r =sqrt(dphi*dphi+deta*deta)
  //H/E condition is also applied in EB.
  //The halo cluster building step targets a large efficiency (ideally >99%) for beam halo deposits.
  //These clusters are used as input for the halo pattern finding methods in EcalHaloAlgo and for the CSC-calo matching methods in GlobalHaloAlgo.     

  //Et threshold hardcoded for now. Might one to get it from config
  std::vector<HaloClusterCandidateECAL> haloclustercands_EB;
  haloclustercands_EB=  GetHaloClusterCandidateEB(TheEBRecHits , TheHBHERecHits, 5);

  std::vector<HaloClusterCandidateECAL> haloclustercands_EE;
  haloclustercands_EE=  GetHaloClusterCandidateEE(TheEERecHits , TheHBHERecHits, 10);

  TheEcalHaloData.setHaloClusterCandidatesEB(haloclustercands_EB);
  TheEcalHaloData.setHaloClusterCandidatesEE(haloclustercands_EE);

  return TheEcalHaloData;
}




std::vector<HaloClusterCandidateECAL> EcalHaloAlgo::GetHaloClusterCandidateEB(edm::Handle<EcalRecHitCollection>& ecalrechitcoll, edm::Handle<HBHERecHitCollection>& hbherechitcoll,float et_thresh_seedrh){

  std::vector<HaloClusterCandidateECAL> TheHaloClusterCandsEB;
  reco::Vertex::Point vtx(0,0,0);

  for(size_t ihit = 0; ihit<ecalrechitcoll->size(); ++ ihit){
    HaloClusterCandidateECAL  clustercand;
    
    const EcalRecHit & rechit = (*ecalrechitcoll)[ ihit ];
    math::XYZPoint rhpos = getPosition(rechit.id(),vtx);
    //Et condition

    double rhet = rechit.energy() * sqrt(rhpos.perp2()/rhpos.mag2());
    if(rhet<et_thresh_seedrh) continue;
    double eta = rhpos.eta();
    double phi = rhpos.phi();
    
    bool isiso = true;
    double etcluster(0);
    int nbcrystalsameeta(0);
    double timediscriminator(0);
    double etstrip_iphiseedplus1(0), etstrip_iphiseedminus1(0);

    //Building the cluster 
    edm::RefVector<EcalRecHitCollection> bhrhcandidates;
    for(size_t jhit = 0; jhit<ecalrechitcoll->size(); ++ jhit){
      const EcalRecHit & rechitj = (*ecalrechitcoll)[ jhit ];
      EcalRecHitRef rhRef(ecalrechitcoll,jhit);
      math::XYZPoint rhposj = getPosition(rechitj.id(),vtx);

      double etaj = rhposj.eta();
      double phij = rhposj.phi();

      double deta = eta - etaj;
      double dphi = deltaPhi(phi,phij);
      if(std::abs(deta)>0.2) continue;//This means +/-11 crystals in eta 
      if(std::abs(dphi)>0.08) continue;//This means +/-4 crystals in phi 

      double rhetj = rechitj.energy()* sqrt(rhposj.perp2()/rhposj.mag2());
      //Rechits with et between 1 and 2 GeV are saved in the rh list but not used in the calculation of the halocluster variables
      if(rhetj<1) continue;
      bhrhcandidates.push_back(rhRef);
      if(rhetj<2) continue;

      if(std::abs(dphi)>0.03){isiso=false;break;}//The strip should be isolated
      if(std::abs(dphi)<0.01) nbcrystalsameeta++;
      if(dphi>0.01) etstrip_iphiseedplus1+=rhetj;
      if(dphi<-0.01) etstrip_iphiseedminus1+=rhetj;
      etcluster+=rhetj;
      //Timing discriminator
      //We assign a weight to the rechit defined as: 
      //Log10(Et)*f(T,R,Z) 
      //where f(T,R,Z) is the separation curve between halo-like and IP-like times.
      //The time difference between a deposit from a outgoing IT halo and a deposit coming from a particle emitted at the IP is given by: 
      //dt= ( - sqrt(R^2+z^2) + |z| )/c
      //Here we take R to be 130 cm. 
      //For EB, the function was parametrized as a function of ieta instead of Z.
      double rhtj = rechitj.time();
      EBDetId detj    = rechitj.id();
      int rhietaj= detj.ieta();
      timediscriminator+= std::log10( rhetj )* ( rhtj +0.5*(sqrt(16900+9*rhietaj*rhietaj)-3*std::abs(rhietaj))/c_cm_per_ns );
      
    }
    //Isolation condition
    if(!isiso) continue;
    
    //Calculate H/E
    double hoe(0);
    for(size_t jhit = 0; jhit<hbherechitcoll->size(); ++ jhit){
      const HBHERecHit & rechitj = (*hbherechitcoll)[ jhit ];
      math::XYZPoint rhposj = getPosition(rechitj.id(),vtx);
      double rhetj = rechitj.energy()* sqrt(rhposj.perp2()/rhposj.mag2());
      if(rhetj<2) continue;
      double etaj = rhposj.eta();
      double phij = rhposj.phi();
      double deta = eta - etaj;
      double dphi = deltaPhi(phi,phij);
      if(std::abs(deta)>0.2) continue;
      if(std::abs(dphi)>0.2) continue;
      hoe+=rhetj/etcluster;
    }
    //H/E condition
    if(hoe>0.1) continue; 
        

    clustercand.setClusterEt(etcluster); 
    clustercand.setSeedEt(rhet); 
    clustercand.setSeedEta(eta); 
    clustercand.setSeedPhi(phi); 
    clustercand.setSeedZ(rhpos.Z());
    clustercand.setSeedR(sqrt(rhpos.perp2())); 
    clustercand.setSeedTime(rechit.time()); 
    clustercand.setHoverE(hoe);
    clustercand.setNbofCrystalsInEta(nbcrystalsameeta);
    clustercand.setEtStripIPhiSeedPlus1(etstrip_iphiseedplus1);
    clustercand.setEtStripIPhiSeedMinus1(etstrip_iphiseedminus1);
    clustercand.setTimeDiscriminator(timediscriminator);
    clustercand.setBeamHaloRecHitsCandidates(bhrhcandidates);


    bool isbeamhalofrompattern = EBClusterShapeandTimeStudy(clustercand,false);
    clustercand.setIsHaloFromPattern(isbeamhalofrompattern);

    bool isbeamhalofrompattern_hlt = EBClusterShapeandTimeStudy(clustercand,true);
    clustercand.setIsHaloFromPattern_HLT(isbeamhalofrompattern_hlt);


    TheHaloClusterCandsEB.push_back(clustercand);
  }

  return  TheHaloClusterCandsEB;
} 




std::vector<HaloClusterCandidateECAL> EcalHaloAlgo::GetHaloClusterCandidateEE(edm::Handle<EcalRecHitCollection>& ecalrechitcoll, edm::Handle<HBHERecHitCollection>& hbherechitcoll,float et_thresh_seedrh){

  std::vector<HaloClusterCandidateECAL> TheHaloClusterCandsEE;

  reco::Vertex::Point vtx(0,0,0);

  for(size_t ihit = 0; ihit<ecalrechitcoll->size(); ++ ihit){
    HaloClusterCandidateECAL  clustercand;
    
    const EcalRecHit & rechit = (*ecalrechitcoll)[ ihit ];
    math::XYZPoint rhpos = getPosition(rechit.id(),vtx);
    //Et condition
    double rhet = rechit.energy()* sqrt(rhpos.perp2()/rhpos.mag2());
    if(rhet<et_thresh_seedrh) continue;
    double eta = rhpos.eta();
    double phi = rhpos.phi();
    double rhr = sqrt(rhpos.perp2());
    
    bool isiso = true;
    double etcluster(0);
    double timediscriminator(0);
    int clustersize(0);
    int nbcrystalssmallt(0);
    int nbcrystalshight(0);
    //Building the cluster
    edm::RefVector<EcalRecHitCollection> bhrhcandidates;
    for(size_t jhit = 0; jhit<ecalrechitcoll->size(); ++ jhit){
      const EcalRecHit & rechitj = (*ecalrechitcoll)[ jhit ];
      EcalRecHitRef rhRef(ecalrechitcoll,jhit);
      math::XYZPoint rhposj = getPosition(rechitj.id(),vtx);

      //Ask the hits to be in the same endcap
      if(rhposj.z()*rhpos.z()<0)continue;

      double etaj = rhposj.eta();
      double phij = rhposj.phi();
      double dr = sqrt((eta-etaj)*(eta-etaj)+deltaPhi(phi,phij)*deltaPhi(phi,phij));

      //Outer cone
      if(dr>0.3) continue;
      
      double rhetj = rechitj.energy()* sqrt(rhposj.perp2()/rhposj.mag2());
      //Rechits with et between 1 and 2 GeV are saved in the rh list but not used in the calculation of the halocluster variables
      if(rhetj<1) continue; 
      bhrhcandidates.push_back(rhRef);
      if(rhetj<2) continue;
      
      //Isolation between outer and inner cone
      if(dr>0.05){isiso=false;break;}//The deposit should be isolated

      etcluster+=rhetj;
  
      //Timing infos:
      //Here we target both IT and OT beam halo 
      double rhtj=rechitj.time();

      //Discriminating variables for OT beam halo: 
      if(rhtj>1) nbcrystalshight++;
      if(rhtj<0) nbcrystalssmallt++;
      //Timing test (likelihood ratio), only for seeds with large R (100 cm) and for crystals with et>5, 
      //This targets IT beam halo (t around - 1ns) 
      if(rhtj>5){
      double corrt_j = rhtj + sqrt(rhposj.x()*rhposj.x()+rhposj.y()*rhposj.y() + 320.*320.)/c_cm_per_ns - 320./c_cm_per_ns;
      //BH is modeled by a Gaussian peaking at 0.
      //Collisions is modeled by a Gaussian peaking at 0.3
      //The width is similar and taken to be 0.4
      timediscriminator+= 0.5*(pow( (corrt_j-0.3)/0.4,2)-pow( (corrt_j-0.)/0.4,2)); 
      clustersize++;
      }
      
    }
    //Isolation condition
    if(!isiso) continue;
    
    //Calculate H2/E
    //Only second hcal layer is considered as it can happen that a shower initiated in EE reaches HCAL first layer
    double h2oe(0);
    for(size_t jhit = 0; jhit<hbherechitcoll->size(); ++ jhit){      
      const HBHERecHit & rechitj = (*hbherechitcoll)[ jhit ];
      math::XYZPoint rhposj = getPosition(rechitj.id(),vtx);

      //Ask the hits to be in the same endcap
      if(rhposj.z()*rhpos.z()<0)continue;
      //Selects only second HCAL layer
      if(std::abs(rhposj.z())<425) continue;
      
      double rhetj = rechitj.energy()* sqrt(rhposj.perp2()/rhposj.mag2());
      if(rhetj<2) continue;
      
      double phij = rhposj.phi();
      if(std::abs(deltaPhi(phi,phij))>0.4 ) continue;
      
      double rhrj = sqrt(rhposj.perp2());
      if(std::abs(rhr-rhrj)>50) continue;

      h2oe+=rhetj/etcluster;
    }
    //H/E condition
    if(h2oe>0.1) continue; 
        

    clustercand.setClusterEt(etcluster); 
    clustercand.setSeedEt(rhet); 
    clustercand.setSeedEta(eta); 
    clustercand.setSeedPhi(phi); 
    clustercand.setSeedZ(rhpos.Z());
    clustercand.setSeedR(sqrt(rhpos.perp2())); 
    clustercand.setSeedTime(rechit.time()); 
    clustercand.setH2overE(h2oe);
    clustercand.setNbEarlyCrystals(nbcrystalssmallt);
    clustercand.setNbLateCrystals(nbcrystalshight);
    clustercand.setClusterSize(clustersize);
    clustercand.setTimeDiscriminator(timediscriminator);
    clustercand.setBeamHaloRecHitsCandidates(bhrhcandidates);
    
    bool isbeamhalofrompattern = EEClusterShapeandTimeStudy_ITBH(clustercand,false) || EEClusterShapeandTimeStudy_OTBH(clustercand,false); 
    clustercand.setIsHaloFromPattern(isbeamhalofrompattern);

    bool isbeamhalofrompattern_hlt = EEClusterShapeandTimeStudy_ITBH(clustercand,true) || EEClusterShapeandTimeStudy_OTBH(clustercand,true); 
    clustercand.setIsHaloFromPattern_HLT(isbeamhalofrompattern_hlt);

    TheHaloClusterCandsEE.push_back(clustercand);
  }

  return  TheHaloClusterCandsEE;
} 





bool EcalHaloAlgo::EBClusterShapeandTimeStudy( HaloClusterCandidateECAL hcand, bool ishlt){
  //Conditions on the central strip size in eta.
  //For low size, extra conditions on seed et, isolation and cluster timing 
  //The time condition only targets IT beam halo. 
  //EB rechits from OT beam halos are typically too late (around 5 ns or more) and seem therefore already cleaned by the reconstruction.

  if(hcand.getSeedEt()<5)return false;
  if(hcand.getNbofCrystalsInEta()<4) return false;
  if(hcand.getNbofCrystalsInEta()==4&&hcand.getSeedEt()<10) return false;
  if(hcand.getNbofCrystalsInEta()==4 && hcand.getEtStripIPhiSeedPlus1()>0.1 &&hcand.getEtStripIPhiSeedMinus1()>0.1 ) return false;
  if(hcand.getNbofCrystalsInEta()<=5 &&  hcand.getTimeDiscriminator()>=0.)return false; 
  
    //For HLT, only use conditions without timing and tighten seed et condition
  if(ishlt &&hcand.getNbofCrystalsInEta()<=5)return false;
  if(ishlt && hcand.getSeedEt()<10)return false;

  hcand.setIsHaloFromPattern(true);
  
  return true;
}



bool EcalHaloAlgo::EEClusterShapeandTimeStudy_OTBH(HaloClusterCandidateECAL hcand, bool ishlt){
  //Separate conditions targeting IT and OT beam halos
  //For OT beam halos, just require enough crystals with large T
  if(hcand.getSeedEt()<20)return false;
  if(hcand.getSeedTime()<0.5)return false;
  if(hcand.getNbLateCrystals()-hcand.getNbEarlyCrystals() <2)return false;
  
  //The use of time information does not allow this method to work at HLT
  if(ishlt)return false;

  hcand.setIsHaloFromPattern(true);

  return true;
}

bool EcalHaloAlgo::EEClusterShapeandTimeStudy_ITBH(HaloClusterCandidateECAL hcand, bool ishlt){
  //Separate conditions targeting IT and OT beam halos 
  //For IT beam halos, fakes from collisions are higher => require the cluster size to be small. 
  //Only halos with R>100 cm are considered here.
  //For lower values, the time difference with particles from collisions is too small
  //IT outgoing beam halos that interact in EE at low R is probably the most difficult category to deal with: 
  //Their signature is very close to the one of photon from collisions (similar cluster shape and timing)
  if(hcand.getSeedEt()<20)return false;
  if(hcand.getSeedR()<100)return false;
  if(hcand.getTimeDiscriminator()<1) return false; 
  if(hcand.getClusterSize()<2) return false;
  if(hcand.getClusterSize()>4) return false;
  
  //The use of time information does not allow this method to work at HLT
  if(ishlt)return false;

  hcand.setIsHaloFromPattern(true);

  return true;
}


math::XYZPoint EcalHaloAlgo::getPosition(const DetId &id, reco::Vertex::Point vtx){

  const GlobalPoint& pos=geo->getPosition(id);
  math::XYZPoint posV(pos.x() - vtx.x(),pos.y() - vtx.y(),pos.z() - vtx.z());
  return posV;
}