CalIsolationAlgo.h

Go to the documentation of this file.

#ifndef IsolationUtils_CalIsolationAlgo_h
#define IsolationUtils_CalIsolationAlgo_h
/* \class CalIsolationAlgo<T1, C2>
 *
 * \author Christian Autermann, U Hamburg
 *
 * template class to calculate calorimeter isolation, the extrapolation
 * to the calorimeter surface is optional.
 *
 */
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/GeometryVector/interface/GlobalVector.h"
#include "DataFormats/Math/interface/deltaR.h"
#include "PhysicsTools/IsolationAlgos/interface/PropagateToCal.h"

template <typename T1, typename C2>
class CalIsolationAlgo {
public:
  typedef double value_type;
  CalIsolationAlgo( ) { }
  CalIsolationAlgo(double dRMin, double dRMax, bool do_propagation,
                   double radius, double minZ, double maxZ, bool theIgnoreMaterial):
     dRMin_( dRMin ), dRMax_( dRMax ), do_propagation_( do_propagation ),
     SrcAtCal(radius, minZ, maxZ, theIgnoreMaterial) { }
  ~CalIsolationAlgo();
  
  void setBfield( const MagneticField * bField ) {
       bField_ = bField;  }
  double operator()(const T1 &, const C2 &) const;

private:
  double dRMin_, dRMax_;
  bool   do_propagation_;
  const MagneticField * bField_;
  PropagateToCal SrcAtCal;
};


template <typename T1, typename C2>
CalIsolationAlgo<T1,C2>::~CalIsolationAlgo() {
}

template <typename T1, typename C2> double CalIsolationAlgo<T1,C2>::
operator()(const T1 & cand, const C2 & elements) const {
  const GlobalPoint Vertex(cand.vx(), cand.vy(), cand.vz());//@@check if this is [cm]!
  //GlobalVector Cand(cand.pt(), cand.eta(), cand.phi()); 
  GlobalVector Cand(cand.px(), cand.py(), cand.pz()); 

  if (do_propagation_ && cand.charge()!=0) 
     SrcAtCal.propagate(Vertex, Cand, cand.charge(), bField_);

  double etSum = 0;
  for( typename C2::const_iterator elem = elements.begin(); 
       elem != elements.end(); ++elem ) {
    double dR = deltaR( elem->eta(), elem->phi(), 
                        (double)Cand.eta(), (double)Cand.phi() );
    if ( dR < dRMax_ && dR > dRMin_ ) {
      etSum += elem->et();
    }
  }
  return etSum;
}

#endif