#include <AlCaIsoTracksProducer.h>

Inheritance diagram for AlCaIsoTracksProducer:

Public Member Functions
	AlCaIsoTracksProducer (const edm::ParameterSet &)
void	endJob (void)
virtual void	produce (edm::Event &, const edm::EventSetup &)
	~AlCaIsoTracksProducer ()
Private Attributes
bool	checkHLTMatch_
double	cluRad_
std::vector< edm::InputTag >	ecalLabels_
double	etaMax_
const CaloGeometry *	geo
edm::InputTag	hbheLabel_
edm::InputTag	hltEventTag_
std::vector< std::string >	hltFiltTag_
double	hltMatchingCone_
edm::InputTag	hoLabel_
double	isolE_
std::string	l1FilterTag_
double	l1jetVetoCone_
double	m_ddirCut
double	m_dvCut
double	m_ecalCut
edm::InputTag	m_inputTrackLabel_
double	m_pCut
double	m_ptCut
int	matrixInnerSize_
int	matrixOuterSize_
int	matrixSize_
int	nHitsMinCore_
int	nHitsMinIso_
TrackAssociatorParameters	parameters_
double	ringInnRad_
double	ringOutRad_
bool	skipNeutrals_
double	taECALCone_
double	taHCALCone_
TrackDetectorAssociator	trackAssociator_
bool	useConeCorr_
bool	useECALCluMatrix_

Detailed Description

Definition at line 51 of file AlCaIsoTracksProducer.h.

Constructor & Destructor Documentation

AlCaIsoTracksProducer::AlCaIsoTracksProducer ( const edm::ParameterSet & iConfig ) [explicit]

Definition at line 138 of file AlCaIsoTracksProducer.cc.

References checkHLTMatch_, cluRad_, ecalLabels_, etaMax_, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), hbheLabel_, hltEventTag_, hltFiltTag_, hltMatchingCone_, hoLabel_, isolE_, l1FilterTag_, l1jetVetoCone_, TrackAssociatorParameters::loadParameters(), m_ddirCut, m_dvCut, m_ecalCut, m_inputTrackLabel_, m_pCut, m_ptCut, matrixSize_, nHitsMinCore_, nHitsMinIso_, Parameters::parameters, parameters_, ringInnRad_, ringOutRad_, skipNeutrals_, AlCaHLTBitMon_QueryRunRegistry::string, taECALCone_, taHCALCone_, trackAssociator_, useConeCorr_, TrackDetectorAssociator::useDefaultPropagator(), and useECALCluMatrix_.

{ 
  
  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 ( )

Definition at line 204 of file AlCaIsoTracksProducer.cc.

{ }

Member Function Documentation

void AlCaIsoTracksProducer::endJob ( void ) [virtual]

Reimplemented from edm::EDProducer.

Definition at line 556 of file AlCaIsoTracksProducer.cc.

{}

void AlCaIsoTracksProducer::produce	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)		`[virtual]`

vector of used hits:

Implements edm::EDProducer.

Definition at line 207 of file AlCaIsoTracksProducer.cc.

References TrackDetectorAssociator::associate(), edm::SortedCollection< T, SORT >::begin(), checkHLTMatch(), checkHLTMatch_, cluRad_, TrackDetMatchInfo::crossedEcalRecHits, deltaR(), TrackAssociatorParameters::dREcal, TrackAssociatorParameters::dRHcal, EcalBarrel, EcalEndcap, ecalLabels_, edm::SortedCollection< T, SORT >::end(), relval_parameters_module::energy, PV3DBase< T, PVType, FrameType >::eta(), etaMax_, geo, edm::EventSetup::get(), edm::Event::getByLabel(), getDistInPlaneTrackDir(), TrackDetectorAssociator::getFreeTrajectoryState(), getL1triggerDirection(), CaloGeometry::getPosition(), edm::Event::getRefBeforePut(), EEDetId::hashedIndex(), EBDetId::hashedIndex(), hbheLabel_, hltEventTag_, hltFiltTag_, hltMatchingCone_, hoLabel_, i, customizeTrackingMonitorSeedNumber::idx, EBDetId::ieta(), info, EBDetId::iphi(), isolE_, EEDetId::ix(), EEDetId::iy(), l1FilterTag_, l1jetVetoCone_, m_ddirCut, m_ecalCut, m_inputTrackLabel_, m_ptCut, matrixSize_, nHitsMinCore_, nHitsMinIso_, parameters_, PV3DBase< T, PVType, FrameType >::phi(), edm::ESHandle< T >::product(), edm::Handle< T >::product(), edm::SortedCollection< T, SORT >::push_back(), edm::Event::put(), ringInnRad_, ringOutRad_, reco::IsolatedPixelTrackCandidate::SetEnergyIn(), reco::IsolatedPixelTrackCandidate::SetEnergyOut(), reco::Track::setExtra(), reco::IsolatedPixelTrackCandidate::SetMaxPtPxl(), reco::IsolatedPixelTrackCandidate::SetSumPtPxl(), skipNeutrals_, mathSSE::sqrt(), taECALCone_, taHCALCone_, trackAssociator_, TrackDetMatchInfo::trkGlobPosAtEcal, TrackDetMatchInfo::trkGlobPosAtHcal, TrackDetMatchInfo::trkMomAtHcal, TrackAssociatorParameters::useCalo, useConeCorr_, TrackAssociatorParameters::useEcal, useECALCluMatrix_, TrackAssociatorParameters::useHcal, and TrackAssociatorParameters::useMuon.

{
  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<HcalDetId> 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;
            for (uint32_t i=0; i<usedHitsHC.size(); i++)
              {
                if (usedHitsHC[i]==hhit->id()) hitIsUsed=true;
              }
            if (hitIsUsed) continue;
            usedHitsHC.push_back(hhit->id());
            
            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");
}