CMSSW_4_4_3_patch1/src/Calibration/HcalAlCaRecoProducers/src/AlCaIsoTracksProducer.cc

Go to the documentation of this file.

#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "Calibration/HcalAlCaRecoProducers/interface/AlCaIsoTracksProducer.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h" 
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
#include "DataFormats/HcalDetId/interface/HcalSubdetector.h" 
#include "DataFormats/GeometrySurface/interface/Cylinder.h"
#include "DataFormats/GeometrySurface/interface/Plane.h" 
#include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
#include "MagneticField/Engine/interface/MagneticField.h" 
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"

#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "FWCore/Utilities/interface/Exception.h"

#include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"

#include "DataFormats/TrackReco/interface/TrackExtra.h"
#include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
#include "RecoTracker/TrackProducer/interface/TrackProducerBase.h"

#include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidate.h"
#include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidateFwd.h"

#include "DataFormats/HLTReco/interface/TriggerEvent.h"

#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"

#include "Math/GenVector/VectorUtil.h"
#include "Math/GenVector/PxPyPzE4D.h"

#include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h"
#include "DataFormats/L1Trigger/interface/L1JetParticle.h"

#include "DataFormats/Provenance/interface/ProductID.h"

#include "Calibration/HcalAlCaRecoProducers/plugins/ConeDefinition.h"

#include <boost/regex.hpp> 

double getDistInCM(double eta1, double phi1, double eta2, double phi2)
{
  double dR, Rec;
  double theta1=2*atan(exp(-eta1));
  double theta2=2*atan(exp(-eta2));
  if (fabs(eta1)<1.479) Rec=129; //radius of ECAL barrel
  else Rec=tan(theta1)*317; //distance from IP to ECAL endcap

  //|vect| times tg of acos(scalar product)
  double angle=acos((sin(theta1)*sin(theta2)*(sin(phi1)*sin(phi2)+cos(phi1)*cos(phi2))+cos(theta1)*cos(theta2)));
  if (angle<acos(-1)/2)
    {
      dR=fabs((Rec/sin(theta1))*tan(angle));
      return dR;
    }
  else return 1000;
}

double getDist(double eta1, double phi1, double eta2, double phi2)
{
  double dphi = fabs(phi1 - phi2); 
  if(dphi>acos(-1)) dphi = 2*acos(-1)-dphi;
  double dr = sqrt(dphi*dphi + std::pow(eta1-eta2,2));
  return dr;
}

bool checkHLTMatch(edm::Event& iEvent, edm::InputTag hltEventTag_, std::vector<std::string> hltFilterTag_, double eta, double phi, double hltMatchingCone_)
{
  bool match =false;
  double minDDD=1000;
  
  edm::Handle<trigger::TriggerEvent> trEv;
  iEvent.getByLabel(hltEventTag_,trEv);
  const trigger::TriggerObjectCollection& TOCol(trEv->getObjects());
  
  trigger::Keys KEYS;
  const trigger::size_type nFilt(trEv->sizeFilters());
  for (trigger::size_type iFilt=0; iFilt!=nFilt; iFilt++) 
    {
      for (unsigned l=0; l<hltFilterTag_.size(); l++)
        {
          if ((trEv->filterTag(iFilt).label()).substr(0,27)==hltFilterTag_[l]) 
            {
              KEYS=trEv->filterKeys(iFilt);
            }
        }
    }
  trigger::size_type nReg=KEYS.size();
  for (trigger::size_type iReg=0; iReg<nReg; iReg++)
    {
      const trigger::TriggerObject& TObj(TOCol[KEYS[iReg]]);
      double dHit=getDist(TObj.eta(),TObj.phi(),eta,phi); 
      if (dHit<minDDD) minDDD=dHit;
    }
  if (minDDD>hltMatchingCone_) match=false;
  else match=true;
     
  return match;
}

std::pair<double,double> getL1triggerDirection(edm::Event& iEvent, edm::InputTag hltEventTag_, std::string l1FilterTag_)
{
  edm::Handle<trigger::TriggerEvent> trEv;
  iEvent.getByLabel(hltEventTag_,trEv);
  const trigger::TriggerObjectCollection& TOCol(trEv->getObjects());
  
  trigger::Keys KEYS;
  const trigger::size_type nFilt(trEv->sizeFilters());
  for (trigger::size_type iFilt=0; iFilt!=nFilt; iFilt++)
    {
      if ((trEv->filterTag(iFilt).label()).substr(0,14)==l1FilterTag_) KEYS=trEv->filterKeys(iFilt); 
    }
  trigger::size_type nReg=KEYS.size();
  double etaTrig=-10000;
  double phiTrig=-10000;
  double ptMax=0;
  for (trigger::size_type iReg=0; iReg<nReg; iReg++)
    {
      const trigger::TriggerObject& TObj(TOCol[KEYS[iReg]]);
      if (TObj.pt()>ptMax)
        {
          etaTrig=TObj.eta();
          phiTrig=TObj.phi();
          ptMax=TObj.pt();
        }
    }
  return std::pair<double,double>(etaTrig,phiTrig);
}


AlCaIsoTracksProducer::AlCaIsoTracksProducer(const edm::ParameterSet& iConfig)
{ 
  
  m_inputTrackLabel_ = iConfig.getParameter<edm::InputTag>("InputTracksLabel");

  hoLabel_ = iConfig.getParameter<edm::InputTag>("HOInput");

  ecalLabels_=iConfig.getParameter<std::vector<edm::InputTag> >("ECALInputs");

  hbheLabel_= iConfig.getParameter<edm::InputTag>("HBHEInput");
  
  m_dvCut = iConfig.getParameter<double>("vtxCut");
  m_ddirCut = iConfig.getParameter<double>("RIsolAtHCALSurface");
  useConeCorr_=iConfig.getParameter<bool>("UseLowPtConeCorrection");
  m_pCut = iConfig.getParameter<double>("MinTrackP");
  m_ptCut = iConfig.getParameter<double>("MinTrackPt");
  m_ecalCut = iConfig.getUntrackedParameter<double>("NeutralIsolCut",8.);

  taECALCone_=iConfig.getUntrackedParameter<double>("TrackAssociatorECALCone",0.5);
  taHCALCone_=iConfig.getUntrackedParameter<double>("TrackAssociatorHCALCone",1.0);

  skipNeutrals_=iConfig.getUntrackedParameter<bool>("SkipNeutralIsoCheck",false);

  nHitsMinCore_=iConfig.getParameter<int>("MinNumberOfHitsCoreTrack");
  nHitsMinIso_=iConfig.getParameter<int>("MinNumberOfHitsInConeTracks");

  isolE_ = iConfig.getParameter<double>("MaxNearbyTrackEnergy");
  etaMax_= iConfig.getParameter<double>("MaxTrackEta");
  cluRad_ = iConfig.getParameter<double>("ECALClusterRadius");
  ringOutRad_ = iConfig.getParameter<double>("ECALRingOuterRadius");
  ringInnRad_=iConfig.getParameter<double>("ECALRingInnerRadius");  

  useECALCluMatrix_=iConfig.getParameter<bool>("ClusterECALasMatrix");
  matrixSize_=iConfig.getParameter<int>("ECALMatrixFullSize");

  checkHLTMatch_=iConfig.getParameter<bool>("CheckHLTMatch");
  hltEventTag_=iConfig.getParameter<edm::InputTag>("hltTriggerEventLabel");  
  hltFiltTag_=iConfig.getParameter<std::vector<std::string> >("hltL3FilterLabels");
  hltMatchingCone_=iConfig.getParameter<double>("hltMatchingCone");
  l1FilterTag_=iConfig.getParameter<std::string>("l1FilterLabel");
  l1jetVetoCone_=iConfig.getParameter<double>("l1JetVetoCone");
  
//
// Parameters for track associator   ===========================
//  
  edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
  parameters_.loadParameters( parameters );
  trackAssociator_.useDefaultPropagator();
// ===============================================================

  //create also IsolatedPixelTrackCandidateCollection which contains isolation info and reference to primary track
  produces<reco::IsolatedPixelTrackCandidateCollection>("HcalIsolatedTrackCollection");

  produces<reco::TrackCollection>("IsoTrackTracksCollection");
  produces<reco::TrackExtraCollection>("IsoTrackExtraTracksCollection");

  produces<EcalRecHitCollection>("IsoTrackEcalRecHitCollection");
  produces<EcalRecHitCollection>("IsoTrackPSEcalRecHitCollection");
  
  produces<HBHERecHitCollection>("IsoTrackHBHERecHitCollection");
  produces<HORecHitCollection>("IsoTrackHORecHitCollection");

}


AlCaIsoTracksProducer::~AlCaIsoTracksProducer() { }


void AlCaIsoTracksProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
{
  std::auto_ptr<reco::IsolatedPixelTrackCandidateCollection> outputHcalIsoTrackColl(new reco::IsolatedPixelTrackCandidateCollection);
  std::auto_ptr<reco::TrackCollection> outputTColl(new reco::TrackCollection);
  std::auto_ptr<reco::TrackExtraCollection> outputExTColl(new reco::TrackExtraCollection);
  
  reco::TrackExtraRefProd rTrackExtras = iEvent.getRefBeforePut<reco::TrackExtraCollection>("IsoTrackExtraTracksCollection"); 
  std::auto_ptr<EcalRecHitCollection> outputEColl(new EcalRecHitCollection);
  std::auto_ptr<EcalRecHitCollection> outputESColl(new EcalRecHitCollection);
  std::auto_ptr<HBHERecHitCollection> outputHColl(new HBHERecHitCollection);
  std::auto_ptr<HORecHitCollection> outputHOColl(new HORecHitCollection);
  
  edm::ESHandle<CaloGeometry> pG;
  iSetup.get<CaloGeometryRecord>().get(pG);
  geo = pG.product();
  
  edm::Handle<reco::TrackCollection> trackCollection;
  iEvent.getByLabel(m_inputTrackLabel_,trackCollection);
  
  
  // temporary collection of EB+EE recHits
  std::auto_ptr<EcalRecHitCollection> tmpEcalRecHitCollection(new EcalRecHitCollection);
  
  std::vector<edm::InputTag>::const_iterator i;
  for (i=ecalLabels_.begin(); i!=ecalLabels_.end(); i++) 
    {
      edm::Handle<EcalRecHitCollection> ec;
      iEvent.getByLabel(*i,ec);
      for(EcalRecHitCollection::const_iterator recHit = (*ec).begin(); recHit != (*ec).end(); ++recHit)
        {
          tmpEcalRecHitCollection->push_back(*recHit);
        }
    }      
  
  edm::Handle<HBHERecHitCollection> hbheRHcol;
  iEvent.getByLabel(hbheLabel_, hbheRHcol);

  const reco::TrackCollection tC = *(trackCollection.product());
  
  int itrk=0;
  int nisotr=0;
  edm::Ref<reco::TrackExtraCollection>::key_type  idx = 0;
  
  //   Parameters for TrackDetAssociator ================================
  // For Low momentum tracks need to look for larger cone for search ====
  // ====================================================================
  
  parameters_.useEcal = true;
  parameters_.useHcal = true;
  parameters_.useCalo = false;
  parameters_.useMuon = false;
  parameters_.dREcal = taECALCone_;
  parameters_.dRHcal = taHCALCone_;
  
  
  std::vector<int> usedHitsHC;
  std::vector<int> usedHitsEC;

  // main loop over input tracks
  for (reco::TrackCollection::const_iterator track=tC.begin(); track!=tC.end(); track++) {
    bool noChargedTracks = true;
    int itrk1=0;
    itrk++;
    double px = track->px();
    double py = track->py();
    double pz = track->pz();
    double ptrack = sqrt(px*px+py*py+pz*pz);
    
    if (ptrack < m_pCut || track->pt() < m_ptCut ) continue; 
 
    if (track->hitPattern().numberOfValidHits() < nHitsMinCore_) continue;

    // check that track is not in the region of L1 jet
    double l1jDR=deltaR(track->eta(), track->phi(), getL1triggerDirection(iEvent,hltEventTag_,l1FilterTag_).first,getL1triggerDirection(iEvent,hltEventTag_,l1FilterTag_).second);
    if (l1jDR<l1jetVetoCone_) continue;
            
    TrackDetMatchInfo info = trackAssociator_.associate(iEvent, iSetup, trackAssociator_.getFreeTrajectoryState(iSetup, *track), parameters_);
    
    GlobalPoint gPointHcal(info.trkGlobPosAtHcal.x(),info.trkGlobPosAtHcal.y(),info.trkGlobPosAtHcal.z());

    GlobalVector trackMomAtHcal = info.trkMomAtHcal;
    
    double etaecal=info.trkGlobPosAtEcal.eta();
    double phiecal=info.trkGlobPosAtEcal.phi();

    if (etaecal==0&&phiecal==0) continue;    

    //check matching to HLT object (optional)

    if (checkHLTMatch_&&!checkHLTMatch(iEvent, hltEventTag_, hltFiltTag_, etaecal, phiecal,hltMatchingCone_)) continue;
    
    if (fabs(track->eta())>etaMax_) continue;
    
    // check charged isolation from all other tracks
    double maxPNearby=-10;
    double sumPNearby=0;
    for (reco::TrackCollection::const_iterator track1=tC.begin(); track1!=tC.end(); track1++)
      {
        itrk1++;
        if (track == track1) continue;
        if (track->hitPattern().numberOfValidHits() < nHitsMinIso_) continue;
        double ptrack1 = sqrt(track1->px()*track1->px()+track1->py()*track1->py()+track1->pz()*track1->pz());
        
        TrackDetMatchInfo info1 = trackAssociator_.associate(iEvent, iSetup, trackAssociator_.getFreeTrajectoryState(iSetup, *track1), parameters_);
        
        GlobalPoint gPointHcal1(info1.trkGlobPosAtHcal.x(),info1.trkGlobPosAtHcal.y(),info1.trkGlobPosAtHcal.z());

        double etaecal1=info1.trkGlobPosAtEcal.eta();
        double phiecal1=info1.trkGlobPosAtEcal.phi();
   
        if (etaecal1==0&&phiecal1==0) continue; 

        double hcDist=getDistInPlaneTrackDir(gPointHcal, trackMomAtHcal, gPointHcal1);

//        double hcDist=getDistInCM(etaecal,phiecal,etaecal1,phiecal1);

        // increase required separation for low momentum tracks
        double factor=1.;
        double factor1=1.;

        if (useConeCorr_)
           {
             if(ptrack<10.)factor+=(10.-ptrack)/20.;
             if(ptrack1<10.)factor1+=(10.-ptrack1)/20.;
           }
        
        if( hcDist <  m_ddirCut*factor*factor1 ) 
          {
            //calculate maximum P and sum P near seed track
            if (track1->p()>maxPNearby)
              {
                maxPNearby=track1->p();
              }
            sumPNearby+=track1->p();
            
            //apply loose isolation criteria
            if (track1->p()>isolE_) 
              {
                noChargedTracks = false;
                break;
              }
          }
        
      } //second track loop

    bool noNeutrals = false;
    
    // we have a good charge-isolated track, so check neutral isolation and write it out
        
    if(noChargedTracks) 
      {
        //find ecal cluster energy and write ecal recHits
        double ecClustR=0;
        double ecClustN=0;
        double ecOutRingR=0;

        //get index of ECAL crystal hit by track
        std::vector<const EcalRecHit*> crossedECids=info.crossedEcalRecHits;
        int etaIDcenter=-10000;
        int phiIDcenter=-10000;
        double enMax=0;
        for (unsigned int i=0; i<crossedECids.size(); i++)
          {
            if ((*crossedECids[i]).id().subdetId()==EcalEndcap)
              {
                EEDetId did(crossedECids[i]->id());
                if (crossedECids[i]->energy()>enMax)
                  {
                    enMax=crossedECids[i]->energy();
                    etaIDcenter=did.iy();
                    phiIDcenter=did.ix();
                  }
              }
            if ((*crossedECids[i]).id().subdetId()==EcalBarrel)
              {
                EBDetId did(crossedECids[i]->id());
                if (crossedECids[i]->energy()>enMax)
                  {
                    enMax=crossedECids[i]->energy();
                    etaIDcenter=did.ieta();
                    phiIDcenter=did.iphi();
                  }
              }
          }
        for (std::vector<EcalRecHit>::const_iterator ehit=tmpEcalRecHitCollection->begin(); ehit!=tmpEcalRecHitCollection->end(); ehit++) 
          {
            // R scheme of ECAL CLUSTERIZATION
            GlobalPoint posH = geo->getPosition((*ehit).detid());
            double phihit = posH.phi();
            double etahit = posH.eta();
            
            double dHit=deltaR(etaecal,phiecal,etahit,phihit);
            
            double dHitCM=getDistInPlaneTrackDir(gPointHcal, trackMomAtHcal, posH);
            
            if (dHitCM<cluRad_)
              {
                ecClustR+=ehit->energy();
              }
            
            if (dHitCM>ringInnRad_&&dHitCM<ringOutRad_)
              {
                ecOutRingR+=ehit->energy();
              }

            //NxN scheme & check whether hit was used or not, if not used push into usedHits
            bool hitIsUsed=false;
            int hitHashedIndex=-10000;
            if (ehit->id().subdetId()==EcalBarrel)
              {
                EBDetId did(ehit->id());
                hitHashedIndex=did.hashedIndex();
                if (fabs(did.ieta()-etaIDcenter)<=matrixSize_/2&&fabs(did.iphi()-phiIDcenter)<=matrixSize_/2) ecClustN+=ehit->energy();
              }
            
            if (ehit->id().subdetId()==EcalEndcap)
              {
                EEDetId did(ehit->id());
                hitHashedIndex=did.hashedIndex();
                if (fabs(did.iy()-etaIDcenter)<=matrixSize_/2&&fabs(did.ix()-phiIDcenter)<=matrixSize_/2) ecClustN+=ehit->energy();
              }
            for (uint32_t i=0; i<usedHitsEC.size(); i++)
              {
                if (usedHitsEC[i]==hitHashedIndex) hitIsUsed=true;
              }

            if (hitIsUsed) continue; //skip if used
            usedHitsEC.push_back(hitHashedIndex);
            
            if(dHit<1.)  
              {
                outputEColl->push_back(*ehit);
              }   
          }

        //check neutrals 
        if (ecOutRingR<m_ecalCut) noNeutrals=true;
        else noNeutrals=false;

        if (noNeutrals||skipNeutrals_)
        {

        // Take info on the track extras and keep it in the outercollection

        reco::TrackExtraRef myextra = (*track).extra();
        reco::TrackExtraRef teref= reco::TrackExtraRef ( rTrackExtras, idx ++ );

        outputTColl->push_back(*track);
        reco::Track & mytrack = outputTColl->back();

        mytrack.setExtra( teref );
        outputExTColl->push_back(*myextra);
        //      reco::TrackExtra & tx = outputExTColl->back();

        //Create IsolatedPixelTrackCandidate (will change naming in future release)
        reco::IsolatedPixelTrackCandidate newHITCandidate(reco::Candidate::LorentzVector(track->px(),track->py(),track->pz(),track->p()));
        newHITCandidate.SetSumPtPxl(sumPNearby);
        newHITCandidate.SetMaxPtPxl(maxPNearby);

        //set cluster energy deposition and ring energy deposition and push_back
        if (!useECALCluMatrix_) newHITCandidate.SetEnergyIn(ecClustR);
        else newHITCandidate.SetEnergyIn(ecClustN);
        newHITCandidate.SetEnergyOut(ecOutRingR);
        outputHcalIsoTrackColl->push_back(newHITCandidate);

        //save hcal recHits
        for (std::vector<HBHERecHit>::const_iterator hhit=hbheRHcol->begin(); hhit!=hbheRHcol->end(); hhit++) 
          {
            //check that this hit was not considered before and push it into usedHits
            bool hitIsUsed=false;
            int hitHashedIndex=hhit->id().hashed_index();
            for (uint32_t i=0; i<usedHitsHC.size(); i++)
              {
                if (usedHitsHC[i]==hitHashedIndex) hitIsUsed=true;
              }
            if (hitIsUsed) continue;
            usedHitsHC.push_back(hitHashedIndex);
            
            GlobalPoint posH = geo->getPosition((*hhit).detid());
            double phihit = posH.phi();
            double etahit = posH.eta();
            
            double dHit=deltaR(etaecal,phiecal,etahit,phihit);
            
            if(dHit<1.)  

              {
                outputHColl->push_back(*hhit);
              }
            }
          }
        
        nisotr++;
        
      } //if (noNeutrals....
  
  } // end of main track loop
  
  if(outputTColl->size() > 0)
    {
      //   Take HO collection
      edm::Handle<HORecHitCollection> ho;
      iEvent.getByLabel(hoLabel_,ho);
      
      const HORecHitCollection Hitho = *(ho.product());
      for(HORecHitCollection::const_iterator hoItr=Hitho.begin(); hoItr!=Hitho.end(); hoItr++)
        {
          outputHOColl->push_back(*hoItr);
        }
           
      // Take Preshower collection      
      
      // get the ps geometry
      edm::ESHandle<CaloGeometry> geoHandle;
      iSetup.get<CaloGeometryRecord>().get(geoHandle);
      
      // get the ps topology
      EcalPreshowerTopology psTopology(geoHandle);
      
      // process rechits
      edm::Handle<EcalRecHitCollection> pRecHits;
      iEvent.getByLabel("ecalPreshowerRecHit","EcalRecHitsES",pRecHits);
      const EcalRecHitCollection& psrechits = *(pRecHits.product());

      typedef EcalRecHitCollection::const_iterator IT;
      
      for(IT i=psrechits.begin(); i!=psrechits.end(); i++) 
        {
          outputESColl->push_back( *i );
        }
    }  

  iEvent.put( outputHcalIsoTrackColl, "HcalIsolatedTrackCollection");
  iEvent.put( outputTColl, "IsoTrackTracksCollection");
  iEvent.put( outputExTColl, "IsoTrackExtraTracksCollection");
  iEvent.put( outputEColl, "IsoTrackEcalRecHitCollection");
  iEvent.put( outputESColl, "IsoTrackPSEcalRecHitCollection");
  iEvent.put( outputHColl, "IsoTrackHBHERecHitCollection");
  iEvent.put( outputHOColl, "IsoTrackHORecHitCollection");
}

void AlCaIsoTracksProducer::endJob(void) {}