/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoTauTag/HLTProducers/src/L2TauIsolationProducer.cc

Go to the documentation of this file.

#include "RecoTauTag/HLTProducers/interface/L2TauIsolationProducer.h"
#include "RecoTauTag/HLTProducers/interface/L2TauIsolationAlgs.h"
#include "RecoTauTag/HLTProducers/interface/L2TauSimpleClustering.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"

using namespace reco;

L2TauIsolationProducer::L2TauIsolationProducer(const edm::ParameterSet& iConfig):
  l2CaloJets_(iConfig.getParameter<edm::InputTag>("L2TauJetCollection")),
  EBRecHits_(iConfig.getParameter<edm::InputTag>("EBRecHits")),
  EERecHits_(iConfig.getParameter<edm::InputTag>("EERecHits")),
  crystalThreshold_(iConfig.getParameter<double>("crystalThreshold")),
  towerThreshold_(iConfig.getParameter<double>("towerThreshold"))
 {
        
  //ECAL Isolation
  edm::ParameterSet ECALIsolParams = iConfig.getParameter<edm::ParameterSet>("ECALIsolation") ;
    
  ECALIsolation_innerCone_ =  ECALIsolParams.getParameter<double>( "innerCone" );
  ECALIsolation_outerCone_ =  ECALIsolParams.getParameter<double>( "outerCone" );
  ECALIsolation_run_       =  ECALIsolParams.getParameter<bool>( "runAlgorithm" );
  

  //ECAL Clustering
  edm::ParameterSet ECALClusterParams = iConfig.getParameter<edm::ParameterSet>("ECALClustering") ;

  ECALClustering_run_                 =  ECALClusterParams.getParameter<bool>( "runAlgorithm" );
  ECALClustering_clusterRadius_       =  ECALClusterParams.getParameter<double>( "clusterRadius" );
    

  //Tower Isolation
  edm::ParameterSet TowerIsolParams = iConfig.getParameter<edm::ParameterSet>("TowerIsolation") ;
  TowerIsolation_innerCone_         =  TowerIsolParams.getParameter<double>( "innerCone" );
  TowerIsolation_outerCone_         =  TowerIsolParams.getParameter<double>( "outerCone" );
  TowerIsolation_run_               =  TowerIsolParams.getParameter<bool>( "runAlgorithm" );
 

  //Add the products
  produces<L2TauInfoAssociation>();

}


L2TauIsolationProducer::~L2TauIsolationProducer()
{
  //Destruction

}



void
L2TauIsolationProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
{


   edm::Handle<CaloJetCollection> l2CaloJets; //Handle to the input (L2TauCaloJets);
   iEvent.getByLabel(l2CaloJets_ ,l2CaloJets);//get the handle

   //Create the Association
   std::auto_ptr<L2TauInfoAssociation> l2InfoAssoc( new L2TauInfoAssociation);

   //If the JetCrystalsAssociation exists -> RUN The Producer
   if(l2CaloJets->size()>0)
     {
      CaloJetCollection::const_iterator jcStart = l2CaloJets->begin();

       //Loop on Jets
       for(CaloJetCollection::const_iterator jc = jcStart ;jc!=l2CaloJets->end();++jc)
         {

           //Create Algorithm Object
           L2TauIsolationAlgs alg;
           
           //Create Info Object
           L2TauIsolationInfo info; 

           //get Hits
           math::PtEtaPhiELorentzVectorCollection hitsECAL = getECALHits(*jc,iEvent,iSetup);
           math::PtEtaPhiELorentzVectorCollection hitsHCAL = getHCALHits(*jc);


           //Run ECALIsolation 
           if(ECALIsolation_run_)
             {
               info.setEcalIsolEt( alg.isolatedEt(hitsECAL , jc->p4().Vect(), ECALIsolation_innerCone_,ECALIsolation_outerCone_) );
               if(hitsECAL.size()>0)
                 info.setSeedEcalHitEt(hitsECAL[0].pt());
             }

           //Run ECALClustering 
           if(ECALClustering_run_)
             {
               //load simple clustering algorithm
               L2TauSimpleClustering clustering(ECALClustering_clusterRadius_);
               math::PtEtaPhiELorentzVectorCollection clusters = clustering.clusterize(hitsECAL);
               info.setEcalClusterShape(alg.clusterShape(clusters,jc->p4().Vect(),0,0.5) );
               info.setNEcalHits(clusters.size());
             }

           //Run CaloTower Isolation
           if(TowerIsolation_run_)
             {
               info.setHcalIsolEt( alg.isolatedEt(hitsHCAL , jc->p4().Vect(), TowerIsolation_innerCone_,TowerIsolation_outerCone_) );
               if(hitsHCAL.size()>0)
                 info.setSeedHcalHitEt(hitsHCAL[0].pt());
             }

           //Store the info Class
           edm::Ref<CaloJetCollection> jcRef(l2CaloJets,jc-jcStart);
           l2InfoAssoc->insert(jcRef, info);
         }

     
     } //end of if(*jetCrystalsObj)

    iEvent.put(l2InfoAssoc);
}


void 
L2TauIsolationProducer::beginJob()
{
}


void 
L2TauIsolationProducer::endJob() {
}



math::PtEtaPhiELorentzVectorCollection 
L2TauIsolationProducer::getHCALHits(const CaloJet& jet)
{


  std::vector<CaloTowerPtr> towers = jet.getCaloConstituents();

  math::PtEtaPhiELorentzVectorCollection towers2;

  for(size_t i=0;i<towers.size();++i)
    if(towers[i]->energy()>towerThreshold_)
      towers2.push_back(math::PtEtaPhiELorentzVector(towers[i]->et(),towers[i]->eta(),towers[i]->phi(),towers[i]->energy()));

  return towers2;
}





math::PtEtaPhiELorentzVectorCollection 
L2TauIsolationProducer::getECALHits(const CaloJet& jet,const edm::Event& iEvent,const edm::EventSetup& iSetup)
{

  //Init Geometry
  edm::ESHandle<CaloGeometry> geometry;
  iSetup.get<CaloGeometryRecord>().get(geometry);

  //Create ECAL Geometry
  const CaloSubdetectorGeometry* EB = geometry->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
  const CaloSubdetectorGeometry* EE = geometry->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);

  //Handle To the ECAL
  edm::Handle<EBRecHitCollection> EBRecHits;
  edm::Handle<EERecHitCollection> EERecHits;

  //Read From File
 
  iEvent.getByLabel( EBRecHits_, EBRecHits );
  iEvent.getByLabel( EERecHits_, EERecHits );
  
 

  //Create a container for the hits
  math::PtEtaPhiELorentzVectorCollection jetRecHits;


  //Create the Machinery to put the highest et crystal at the begining
  double ref_pt = 0.;//max et
  int ref_id=0; //id of the max element
  int count=0; //counter

  //Code Copied from JetCrystalsAssociator    
      const std::vector<CaloTowerPtr>  myTowers=jet.getCaloConstituents();
        for (unsigned int iTower = 0; iTower < myTowers.size(); iTower++)
        {
          CaloTowerPtr theTower = myTowers[iTower];
          size_t numRecHits = theTower->constituentsSize();
          // access CaloRecHits
          for (size_t j = 0; j < numRecHits; j++) {
            DetId RecHitDetID=theTower->constituent(j);
            DetId::Detector DetNum=RecHitDetID.det();
            if( DetNum == DetId::Ecal ){
              int EcalNum =  RecHitDetID.subdetId();
              if( EcalNum == 1 ){
                EBDetId EcalID = RecHitDetID;
                EBRecHitCollection::const_iterator theRecHit=EBRecHits->find(EcalID);
                if(theRecHit != EBRecHits->end()){
                  DetId id = theRecHit->detid();
                  const CaloCellGeometry* this_cell = EB->getGeometry(id);
                  if (this_cell) {
                    GlobalPoint posi = this_cell->getPosition();
                    double energy = theRecHit->energy();
                    double eta = posi.eta();
                    double phi = posi.phi();
                    double theta = posi.theta();
                    if(theta > M_PI) theta = 2 * M_PI- theta;
                    double et = energy * sin(theta);
                    //Apply Thresholds Here
                    math::PtEtaPhiELorentzVector p(et, eta, phi, energy);
                     if(p.pt()>crystalThreshold_)
                      { 
                        if(p.pt()>ref_pt)
                         {
                           ref_id=count;
                           ref_pt = p.pt();
                         }
                        jetRecHits.push_back(p);
                        count++;
                      }
                    
                  }
                }
              } else if ( EcalNum == 2 ) {
                EEDetId EcalID = RecHitDetID;
                EERecHitCollection::const_iterator theRecHit=EERecHits->find(EcalID);    
                if(theRecHit != EBRecHits->end()){
                  DetId id = theRecHit->detid();
                  const CaloCellGeometry* this_cell = EE->getGeometry(id);
                  if (this_cell) {
                    GlobalPoint posi = this_cell->getPosition();
                    double energy = theRecHit->energy();
                    double eta = posi.eta();
                    double phi = posi.phi();
                    double theta = posi.theta();
                    if (theta > M_PI) theta = 2 * M_PI - theta;
                    double et = energy * sin(theta);
                    // std::cout <<"Et "<<et<<std::endl;
                    math::PtEtaPhiELorentzVector p(et, eta, phi, energy);
                     if(p.pt()>crystalThreshold_)
                      { 
                        if(p.pt()>ref_pt)
                         {
                           ref_id=count;
                           ref_pt = p.pt();
                         }
                        jetRecHits.push_back(p);
                        count++;
                      }
                  }
                }
              }
            }
          }
        }

        //bring it to the front
        if(jetRecHits.size()>0)
          std::swap(jetRecHits[ref_id],jetRecHits[0]);

        return jetRecHits;
}