EcalHaloAlgo Class Reference

#include <EcalHaloAlgo.h>

Public Member Functions
reco::EcalHaloData	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)
bool	EBClusterShapeandTimeStudy (reco::HaloClusterCandidateECAL hcand, bool ishlt)
	EcalHaloAlgo ()
bool	EEClusterShapeandTimeStudy_ITBH (reco::HaloClusterCandidateECAL hcand, bool ishlt)
bool	EEClusterShapeandTimeStudy_OTBH (reco::HaloClusterCandidateECAL hcand, bool ishlt)
float	GetAngleCut ()
float	GetEBRecHitEnergyThreshold ()
float	GetEERecHitEnergyThreshold ()
float	GetESRecHitEnergyThreshold ()
std::vector < reco::HaloClusterCandidateECAL >	GetHaloClusterCandidateEB (edm::Handle< EcalRecHitCollection > &ecalrechitcoll, edm::Handle< HBHERecHitCollection > &hbherechitcoll, float et_thresh_seedrh)
std::vector < reco::HaloClusterCandidateECAL >	GetHaloClusterCandidateEE (edm::Handle< EcalRecHitCollection > &ecalrechitcoll, edm::Handle< HBHERecHitCollection > &hbherechitcoll, float et_thresh_seedrh)
float	GetPhiWedgeEnergyThreshold ()
int	GetPhiWedgeNHitsThreshold ()
float	GetRoundnessCut ()
void	SetAngleCut (float a=4.)
void	SetPhiWedgeEnergyThreshold (float SumE)
void	SetPhiWedgeNHitsThreshold (int nhits)
void	SetPhiWedgeThresholds (float SumE, int nhits)
void	SetRecHitEnergyThresholds (float EB, float EE, float ES)
void	SetRoundnessCut (float r=100.)
	~EcalHaloAlgo ()

Private Member Functions
math::XYZPoint	getPosition (const DetId &id, reco::Vertex::Point vtx)

Private Attributes
float	AngleCut
float	EBRecHitEnergyThreshold
float	EERecHitEnergyThreshold
float	ESRecHitEnergyThreshold
const CaloGeometry *	geo
int	NHitsThreshold
float	RoundnessCut
float	SumEnergyThreshold

Detailed Description

Definition at line 44 of file EcalHaloAlgo.h.

Constructor & Destructor Documentation

EcalHaloAlgo::EcalHaloAlgo

(

)

Definition at line 19 of file EcalHaloAlgo.cc.

{
  RoundnessCut = 0 ;
  AngleCut = 0 ;
  EBRecHitEnergyThreshold = 0.;
  EERecHitEnergyThreshold = 0.;
  ESRecHitEnergyThreshold = 0.;
  SumEnergyThreshold = 0.;
  NHitsThreshold =0;

  geo = 0;
}

EcalHaloAlgo::~EcalHaloAlgo ( )

inline

Definition at line 50 of file EcalHaloAlgo.h.

{}

Member Function Documentation

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
	)

Definition at line 32 of file EcalHaloAlgo.cc.

References funct::abs(), Abs(), angle(), CompareTime(), DetId::Ecal, edm::EventSetup::get(), reco::EcalHaloData::GetPhiWedges(), reco::EcalHaloData::GetShowerShapesAngle(), reco::EcalHaloData::GetShowerShapesRoundness(), CaloGeometry::getSubdetectorGeometry(), reco::EcalHaloData::GetSuperClusters(), i, hit::id, EBDetId::ieta(), edm::helper::Filler< Map >::insert(), EBDetId::iphi(), edm::HandleBase::isValid(), j, edm::Handle< T >::product(), edm::ESHandle< class >::product(), edm::RefVector< C, T, F >::push_back(), DetId::rawId(), reco::EcalHaloData::setHaloClusterCandidatesEB(), reco::EcalHaloData::setHaloClusterCandidatesEE(), and reco::PhiWedge::SetPlusZOriginConfidence().

Referenced by reco::EcalHaloDataProducer::produce().

{  

  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()); 
      const CaloSubdetectorGeometry* TheSubGeometry = 0;                                                                         
      const CaloCellGeometry* cell = 0 ;                                                                                    

      // Get EB geometry 
      TheSubGeometry = TheCaloGeometry.getSubdetectorGeometry(DetId::Ecal, 1);                                                         
      EBDetId EcalID(id.rawId());
      if( TheSubGeometry )
    cell = TheSubGeometry->getGeometry(id);
  
      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();
  }



  geo = 0;
  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;
}

bool EcalHaloAlgo::EBClusterShapeandTimeStudy	(	reco::HaloClusterCandidateECAL	hcand,
		bool	ishlt
	)

Definition at line 473 of file EcalHaloAlgo.cc.

References reco::HaloClusterCandidateECAL::getEtStripIPhiSeedMinus1(), reco::HaloClusterCandidateECAL::getEtStripIPhiSeedPlus1(), reco::HaloClusterCandidateECAL::getNbofCrystalsInEta(), reco::HaloClusterCandidateECAL::getSeedEt(), reco::HaloClusterCandidateECAL::getTimeDiscriminator(), and reco::HaloClusterCandidateECAL::setIsHaloFromPattern().

                                                                                        {
  //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_ITBH	(	reco::HaloClusterCandidateECAL	hcand,
		bool	ishlt
	)

Definition at line 511 of file EcalHaloAlgo.cc.

References reco::HaloClusterCandidateECAL::getClusterSize(), reco::HaloClusterCandidateECAL::getSeedEt(), reco::HaloClusterCandidateECAL::getSeedR(), reco::HaloClusterCandidateECAL::getTimeDiscriminator(), and reco::HaloClusterCandidateECAL::setIsHaloFromPattern().

                                                                                            {
  //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;
}

bool EcalHaloAlgo::EEClusterShapeandTimeStudy_OTBH	(	reco::HaloClusterCandidateECAL	hcand,
		bool	ishlt
	)

Definition at line 496 of file EcalHaloAlgo.cc.

References reco::HaloClusterCandidateECAL::getNbEarlyCrystals(), reco::HaloClusterCandidateECAL::getNbLateCrystals(), reco::HaloClusterCandidateECAL::getSeedEt(), reco::HaloClusterCandidateECAL::getSeedTime(), and reco::HaloClusterCandidateECAL::setIsHaloFromPattern().

                                                                                            {
  //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;
}

float EcalHaloAlgo::GetAngleCut ( )

inline

Definition at line 78 of file EcalHaloAlgo.h.

References AngleCut.

{return AngleCut ;}

float EcalHaloAlgo::GetEBRecHitEnergyThreshold ( )

inline

Definition at line 81 of file EcalHaloAlgo.h.

References EBRecHitEnergyThreshold.

{ return EBRecHitEnergyThreshold;}

float EcalHaloAlgo::GetEERecHitEnergyThreshold ( )

inline

Definition at line 82 of file EcalHaloAlgo.h.

References EERecHitEnergyThreshold.

{ return EERecHitEnergyThreshold;}

float EcalHaloAlgo::GetESRecHitEnergyThreshold ( )

inline

Definition at line 83 of file EcalHaloAlgo.h.

References ESRecHitEnergyThreshold.

{ return ESRecHitEnergyThreshold;}

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

Definition at line 226 of file EcalHaloAlgo.cc.

References funct::abs(), SiPixelRawToDigiRegional_cfi::deltaPhi, CaloRecHit::energy(), EcalRecHit::energy(), eta, HBHERecHit::id(), EcalRecHit::id(), EBDetId::ieta(), phi(), edm::RefVector< C, T, F >::push_back(), reco::HaloClusterCandidateECAL::setBeamHaloRecHitsCandidates(), reco::HaloClusterCandidateECAL::setClusterEt(), reco::HaloClusterCandidateECAL::setEtStripIPhiSeedMinus1(), reco::HaloClusterCandidateECAL::setEtStripIPhiSeedPlus1(), reco::HaloClusterCandidateECAL::setHoverE(), reco::HaloClusterCandidateECAL::setIsHaloFromPattern(), reco::HaloClusterCandidateECAL::setIsHaloFromPattern_HLT(), reco::HaloClusterCandidateECAL::setNbofCrystalsInEta(), reco::HaloClusterCandidateECAL::setSeedEt(), reco::HaloClusterCandidateECAL::setSeedEta(), reco::HaloClusterCandidateECAL::setSeedPhi(), reco::HaloClusterCandidateECAL::setSeedR(), reco::HaloClusterCandidateECAL::setSeedTime(), reco::HaloClusterCandidateECAL::setSeedZ(), reco::HaloClusterCandidateECAL::setTimeDiscriminator(), mathSSE::sqrt(), and EcalRecHit::time().

                                                                                                                                                                                                  {

  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
	)

Definition at line 342 of file EcalHaloAlgo.cc.

References funct::abs(), SiPixelRawToDigiRegional_cfi::deltaPhi, CaloRecHit::energy(), EcalRecHit::energy(), eta, HBHERecHit::id(), EcalRecHit::id(), phi(), funct::pow(), edm::RefVector< C, T, F >::push_back(), reco::HaloClusterCandidateECAL::setBeamHaloRecHitsCandidates(), reco::HaloClusterCandidateECAL::setClusterEt(), reco::HaloClusterCandidateECAL::setClusterSize(), reco::HaloClusterCandidateECAL::setH2overE(), reco::HaloClusterCandidateECAL::setIsHaloFromPattern(), reco::HaloClusterCandidateECAL::setIsHaloFromPattern_HLT(), reco::HaloClusterCandidateECAL::setNbEarlyCrystals(), reco::HaloClusterCandidateECAL::setNbLateCrystals(), reco::HaloClusterCandidateECAL::setSeedEt(), reco::HaloClusterCandidateECAL::setSeedEta(), reco::HaloClusterCandidateECAL::setSeedPhi(), reco::HaloClusterCandidateECAL::setSeedR(), reco::HaloClusterCandidateECAL::setSeedTime(), reco::HaloClusterCandidateECAL::setSeedZ(), reco::HaloClusterCandidateECAL::setTimeDiscriminator(), mathSSE::sqrt(), and EcalRecHit::time().

                                                                                                                                                                                                  {

  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;
} 

float EcalHaloAlgo::GetPhiWedgeEnergyThreshold ( )

inline

Definition at line 86 of file EcalHaloAlgo.h.

References SumEnergyThreshold.

{ return SumEnergyThreshold;}

int EcalHaloAlgo::GetPhiWedgeNHitsThreshold ( )

inline

Definition at line 87 of file EcalHaloAlgo.h.

References NHitsThreshold.

{ return NHitsThreshold;}

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

private

Definition at line 533 of file EcalHaloAlgo.cc.

References PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

                                                                          {

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

float EcalHaloAlgo::GetRoundnessCut ( )

inline

Definition at line 76 of file EcalHaloAlgo.h.

References RoundnessCut.

{return RoundnessCut ;}

void EcalHaloAlgo::SetAngleCut ( float  a = 4. )

inline

Definition at line 57 of file EcalHaloAlgo.h.

References a, and AngleCut.

Referenced by reco::EcalHaloDataProducer::produce().

{ AngleCut =a;}     

void EcalHaloAlgo::SetPhiWedgeEnergyThreshold ( float  SumE )

inline

Definition at line 62 of file EcalHaloAlgo.h.

References SumEnergyThreshold.

{ SumEnergyThreshold = SumE ;}

void EcalHaloAlgo::SetPhiWedgeNHitsThreshold ( int  nhits )

inline

Definition at line 63 of file EcalHaloAlgo.h.

References nhits, and NHitsThreshold.

{ NHitsThreshold = nhits ; }

void EcalHaloAlgo::SetPhiWedgeThresholds ( float  SumE, int  nhits  )

inline

Definition at line 64 of file EcalHaloAlgo.h.

References nhits, NHitsThreshold, and SumEnergyThreshold.

Referenced by reco::EcalHaloDataProducer::produce().

{ SumEnergyThreshold = SumE ; NHitsThreshold = nhits ;}

void EcalHaloAlgo::SetRecHitEnergyThresholds ( float  EB, float  EE, float  ES  )

inline

Definition at line 60 of file EcalHaloAlgo.h.

References EB, EBRecHitEnergyThreshold, EE, EERecHitEnergyThreshold, and ESRecHitEnergyThreshold.

Referenced by reco::EcalHaloDataProducer::produce().

{ EBRecHitEnergyThreshold = EB; EERecHitEnergyThreshold = EE ; ESRecHitEnergyThreshold = ES;}

void EcalHaloAlgo::SetRoundnessCut ( float  r = 100. )

inline

Definition at line 55 of file EcalHaloAlgo.h.

References alignCSCRings::r, and RoundnessCut.

Referenced by reco::EcalHaloDataProducer::produce().

{ RoundnessCut=r; }                                       

Member Data Documentation

float EcalHaloAlgo::AngleCut

private

Definition at line 95 of file EcalHaloAlgo.h.

Referenced by GetAngleCut(), and SetAngleCut().

float EcalHaloAlgo::EBRecHitEnergyThreshold

private

Definition at line 98 of file EcalHaloAlgo.h.

Referenced by GetEBRecHitEnergyThreshold(), and SetRecHitEnergyThresholds().

float EcalHaloAlgo::EERecHitEnergyThreshold

private

Definition at line 99 of file EcalHaloAlgo.h.

Referenced by GetEERecHitEnergyThreshold(), and SetRecHitEnergyThresholds().

float EcalHaloAlgo::ESRecHitEnergyThreshold

private

Definition at line 100 of file EcalHaloAlgo.h.

Referenced by GetESRecHitEnergyThreshold(), and SetRecHitEnergyThresholds().

const CaloGeometry* EcalHaloAlgo::geo

private

Definition at line 106 of file EcalHaloAlgo.h.

int EcalHaloAlgo::NHitsThreshold

private

Definition at line 104 of file EcalHaloAlgo.h.

Referenced by GetPhiWedgeNHitsThreshold(), SetPhiWedgeNHitsThreshold(), and SetPhiWedgeThresholds().

float EcalHaloAlgo::RoundnessCut

private

Definition at line 93 of file EcalHaloAlgo.h.

Referenced by GetRoundnessCut(), and SetRoundnessCut().

float EcalHaloAlgo::SumEnergyThreshold

private

Definition at line 103 of file EcalHaloAlgo.h.

Referenced by GetPhiWedgeEnergyThreshold(), SetPhiWedgeEnergyThreshold(), and SetPhiWedgeThresholds().