#include <CaloExtractorByAssociator.h>

Inheritance diagram for muonisolation::CaloExtractorByAssociator:

Public Member Functions
	CaloExtractorByAssociator ()
	constructors
	CaloExtractorByAssociator (const edm::ParameterSet &par)
virtual reco::IsoDeposit	deposit (const edm::Event &ev, const edm::EventSetup &evSetup, const reco::Track &track) const
	no-op: by design of this extractor the deposits are pulled out all at a time
virtual std::vector < reco::IsoDeposit >	deposits (const edm::Event &ev, const edm::EventSetup &evSetup, const reco::Track &track) const
	return deposits for 3 calorimeter subdetectors (ecal, hcal, ho) -- in this order
virtual void	fillVetos (const edm::Event &ev, const edm::EventSetup &evSetup, const reco::TrackCollection &tracks)
	allows to set extra vetoes (in addition to the muon) -- no-op at this point
virtual	~CaloExtractorByAssociator ()
	destructor
Private Member Functions
double	noiseEcal (const CaloTower &tower) const
	Determine noise for HCAL and ECAL (take some defaults for the time being)
double	noiseHcal (const CaloTower &tower) const
double	noiseHOcal (const CaloTower &tower) const
double	noiseRecHit (const DetId &detId) const
Private Attributes
TrackDetectorAssociator *	theAssociator
TrackAssociatorParameters *	theAssociatorParameters
	associator, its' parameters and the propagator
bool	theCenterConeOnCalIntersection
std::vector< std::string >	theDepositInstanceLabels
	multiple deposits: labels -- expect 3 labels beginning with "e", "h", "ho"
std::string	theDepositLabel
	Label of deposit -- suggest to set to "" (all info is in collection name anyways)
double	theDR_Max
	max cone size in which towers are considered
double	theDR_Veto_E
	cone sizes inside which the Et (towers) are not counted
double	theDR_Veto_H
double	theDR_Veto_HO
double	theNoise_EB
double	theNoise_EE
double	theNoise_HB
double	theNoise_HE
double	theNoise_HO
double	theNoiseTow_EB
double	theNoiseTow_EE
bool	thePrintTimeReport
	flag to turn on/off printing of a time report
std::string	thePropagatorName
	propagator name to feed into the track associator
MuonServiceProxy *	theService
	the event setup proxy, it takes care the services update
double	theThreshold_E
double	theThreshold_H
double	theThreshold_HO
bool	theUseRecHitsFlag
	use towers or rec hits
std::vector< DetId >	theVetoCollection
	Vector of calo Ids to veto -- not used.

Detailed Description

Definition at line 36 of file CaloExtractorByAssociator.h.

Constructor & Destructor Documentation

muonisolation::CaloExtractorByAssociator::CaloExtractorByAssociator ( ) [inline]

constructors

Definition at line 41 of file CaloExtractorByAssociator.h.

{};

CaloExtractorByAssociator::CaloExtractorByAssociator ( const edm::ParameterSet & par )

Definition at line 37 of file CaloExtractorByAssociator.cc.

References edm::ParameterSet::getParameter(), MuonServiceProxy_cff::MuonServiceProxy, theAssociator, theAssociatorParameters, theService, and example_cfg::TrackAssociatorParameters.

                                                                            :
  theUseRecHitsFlag(par.getParameter<bool>("UseRecHitsFlag")),
  theDepositLabel(par.getUntrackedParameter<string>("DepositLabel")),
  theDepositInstanceLabels(par.getParameter<std::vector<std::string> >("DepositInstanceLabels")),
  thePropagatorName(par.getParameter<std::string>("PropagatorName")),
  theThreshold_E(par.getParameter<double>("Threshold_E")),
  theThreshold_H(par.getParameter<double>("Threshold_H")),
  theThreshold_HO(par.getParameter<double>("Threshold_HO")),
  theDR_Veto_E(par.getParameter<double>("DR_Veto_E")),
  theDR_Veto_H(par.getParameter<double>("DR_Veto_H")),
  theDR_Veto_HO(par.getParameter<double>("DR_Veto_HO")),
  theCenterConeOnCalIntersection(par.getParameter<bool>("CenterConeOnCalIntersection")),
  theDR_Max(par.getParameter<double>("DR_Max")),
  theNoise_EB(par.getParameter<double>("Noise_EB")),
  theNoise_EE(par.getParameter<double>("Noise_EE")),
  theNoise_HB(par.getParameter<double>("Noise_HB")),
  theNoise_HE(par.getParameter<double>("Noise_HE")),
  theNoise_HO(par.getParameter<double>("Noise_HO")),    
  theNoiseTow_EB(par.getParameter<double>("NoiseTow_EB")),
  theNoiseTow_EE(par.getParameter<double>("NoiseTow_EE")),
  theService(0),
  theAssociator(0),
  thePrintTimeReport(par.getUntrackedParameter<bool>("PrintTimeReport"))
{
  ParameterSet serviceParameters = par.getParameter<ParameterSet>("ServiceParameters");
  theService = new MuonServiceProxy(serviceParameters);

  theAssociatorParameters = new TrackAssociatorParameters(par.getParameter<edm::ParameterSet>("TrackAssociatorParameters"));
  theAssociator = new TrackDetectorAssociator();
}

CaloExtractorByAssociator::~CaloExtractorByAssociator ( ) [virtual]

destructor

Definition at line 68 of file CaloExtractorByAssociator.cc.

References gather_cfg::cout, cond::rpcobimon::current, theAssociator, theAssociatorParameters, thePrintTimeReport, and theService.

                                                     {
  if (thePrintTimeReport) TimingReport::current()->dump(std::cout);
  if (theAssociatorParameters) delete theAssociatorParameters;
  if (theService) delete theService;
  if (theAssociator) delete theAssociator;
}

Member Function Documentation

IsoDeposit CaloExtractorByAssociator::deposit	(	const edm::Event &	ev,
		const edm::EventSetup &	evSetup,
		const reco::Track &	track
	)		const `[virtual]`

no-op: by design of this extractor the deposits are pulled out all at a time

Implements reco::isodeposit::IsoDepositExtractor.

Definition at line 83 of file CaloExtractorByAssociator.cc.

References reco::TrackBase::eta(), and reco::TrackBase::phi().

{
  IsoDeposit::Direction muonDir(muon.eta(), muon.phi());
  IsoDeposit dep(muonDir );

//   LogWarning("CaloExtractorByAssociator")
//     <<"single deposit is not an option here\n"
//     <<"use ::deposits --> extract all and reweight as necessary";

  return dep;

}

std::vector< IsoDeposit > CaloExtractorByAssociator::deposits	(	const edm::Event &	ev,
		const edm::EventSetup &	evSetup,
		const reco::Track &	track
	)		const `[virtual]`

return deposits for 3 calorimeter subdetectors (ecal, hcal, ho) -- in this order

Make separate deposits: for ECAL, HCAL, HO.

check configuration consistency could've been made at construction stage (fix later?)

this should be (eventually) set to the eta-phi of the crossing point of a straight line tangent to a muon at IP and the calorimeter

each deposit type veto is at the point of intersect with that detector

do things based on rec-hits here too much copy-pasting now (refactor later?)

first check if the hit is inside the veto cone by dR-alone

and now pitch those in the crossed list

first check if the hit is inside the veto cone by dR-alone

and now pitch those in the crossed list

first check if the hit is inside the veto cone by dR-alone

and now pitch those in the crossed list

use calo towers

first check if the tower is inside the veto cone by dR-alone

and now pitch those in the crossed list

add the Et of the tower to deposits if it's not a vetoed; put into muonEnergy otherwise

Reimplemented from reco::isodeposit::IsoDepositExtractor.

Definition at line 98 of file CaloExtractorByAssociator.cc.

{
  theService->update(eventSetup);
  theAssociator->setPropagator(&*(theService->propagator(thePropagatorName)));

  if (theDepositInstanceLabels.size() != 3){
    LogError("MuonIsolation")<<"Configuration is inconsistent: Need 3 deposit instance labels";
  }
  if (! theDepositInstanceLabels[0].compare(0,1, std::string("e")) == 0
      || ! theDepositInstanceLabels[1].compare(0,1, std::string("h")) == 0
      || ! theDepositInstanceLabels[2].compare(0,2, std::string("ho")) == 0){
    LogWarning("MuonIsolation")<<"Deposit instance labels do not look like  (e*, h*, ho*):"
                               <<"proceed at your own risk. The extractor interprets lab0=from ecal; lab1=from hcal; lab2=from ho";
  }

  typedef IsoDeposit::Veto Veto;
  IsoDeposit::Direction muonDir(muon.eta(), muon.phi());
  
  IsoDeposit depEcal(muonDir);
  IsoDeposit depHcal(muonDir);
  IsoDeposit depHOcal(muonDir);

  edm::ESHandle<MagneticField> bField;
  eventSetup.get<IdealMagneticFieldRecord>().get(bField);


  reco::TransientTrack tMuon(muon, &*bField);
  FreeTrajectoryState iFTS = tMuon.initialFreeState();
  TrackDetMatchInfo mInfo = theAssociator->associate(event, eventSetup, iFTS, *theAssociatorParameters);

  depEcal.setVeto(Veto(reco::isodeposit::Direction(mInfo.trkGlobPosAtEcal.eta(), mInfo.trkGlobPosAtEcal.phi()),
                       theDR_Veto_E));
  depHcal.setVeto(Veto(reco::isodeposit::Direction(mInfo.trkGlobPosAtHcal.eta(), mInfo.trkGlobPosAtHcal.phi()),
                       theDR_Veto_H));
  depHOcal.setVeto(Veto(reco::isodeposit::Direction(mInfo.trkGlobPosAtHO.eta(), mInfo.trkGlobPosAtHO.phi()),
                        theDR_Veto_HO));

  if (theCenterConeOnCalIntersection){
    reco::isodeposit::Direction dirTmp = depEcal.veto().vetoDir;
    double dRtmp = depEcal.veto().dR;
    depEcal = IsoDeposit(dirTmp); depEcal.setVeto(Veto(dirTmp, dRtmp));

    dirTmp = depHcal.veto().vetoDir;
    dRtmp = depHcal.veto().dR;
    depHcal = IsoDeposit(dirTmp); depHcal.setVeto(Veto(dirTmp, dRtmp));

    dirTmp = depHOcal.veto().vetoDir;
    dRtmp = depHOcal.veto().dR;
    depHOcal = IsoDeposit(dirTmp); depHOcal.setVeto(Veto(dirTmp, dRtmp));
  }

  if (theUseRecHitsFlag){
    edm::ESHandle<CaloGeometry> caloGeom;
    eventSetup.get<CaloGeometryRecord>().get(caloGeom);

    //Ecal
    std::vector<const EcalRecHit*>::const_iterator eHitCI = mInfo.ecalRecHits.begin();
    for (; eHitCI != mInfo.ecalRecHits.end(); ++eHitCI){
      const EcalRecHit* eHitCPtr = *eHitCI;
      GlobalPoint eHitPos = caloGeom->getPosition(eHitCPtr->detid());
      double deltar0 = reco::deltaR(muon, eHitPos);
      double cosTheta = 1./cosh(eHitPos.eta());
      double energy = eHitCPtr->energy();
      double et = energy*cosTheta; 
      if (deltar0 > theDR_Max 
          || ! (et > theThreshold_E && energy > 3*noiseRecHit(eHitCPtr->detid()))) continue;

      bool vetoHit = false;
      double deltar = reco::deltaR(mInfo.trkGlobPosAtEcal, eHitPos);
      if (deltar < theDR_Veto_E ){
        LogDebug("RecoMuon|CaloExtractorByAssociator")
          << " >>> Veto ECAL hit: Calo deltaR= " << deltar;
        LogDebug("RecoMuon|CaloExtractorByAssociator")
          << " >>> Calo eta phi ethcal: " << eHitPos.eta() << " " << eHitPos.phi() << " " << et;
        vetoHit = true;
      }
      if (! vetoHit){
        for (unsigned int iH = 0; iH< mInfo.crossedEcalIds.size() && ! vetoHit; ++iH){
          if (mInfo.crossedEcalIds[iH].rawId() == eHitCPtr->detid().rawId()) vetoHit = true;
        }
      }

      if (vetoHit ){
        depEcal.addCandEnergy(et);
      } else {
        depEcal.addDeposit(reco::isodeposit::Direction(eHitPos.eta(), eHitPos.phi()), et);      
      }
    }

    //Hcal
    std::vector<const HBHERecHit*>::const_iterator hHitCI = mInfo.hcalRecHits.begin();
    for (; hHitCI != mInfo.hcalRecHits.end(); ++hHitCI){
      const HBHERecHit* hHitCPtr = *hHitCI;
      GlobalPoint hHitPos = caloGeom->getPosition(hHitCPtr->detid());
      double deltar0 = reco::deltaR(muon, hHitPos);
      double cosTheta = 1./cosh(hHitPos.eta());
      double energy = hHitCPtr->energy();
      double et = energy*cosTheta;
      if (deltar0 > theDR_Max 
          || ! (et > theThreshold_H && energy > 3*noiseRecHit(hHitCPtr->detid()))) continue;

      bool vetoHit = false;
      double deltar = reco::deltaR(mInfo.trkGlobPosAtHcal, hHitPos);
      if (deltar < theDR_Veto_H ){
        LogDebug("RecoMuon|CaloExtractorByAssociator")
          << " >>> Veto HBHE hit: Calo deltaR= " << deltar;
        LogDebug("RecoMuon|CaloExtractorByAssociator")
          << " >>> Calo eta phi ethcal: " << hHitPos.eta() << " " << hHitPos.phi() << " " << et;
        vetoHit = true;
      }
      if (! vetoHit){
        for (unsigned int iH = 0; iH< mInfo.crossedHcalIds.size() && ! vetoHit; ++iH){
          if (mInfo.crossedHcalIds[iH].rawId() == hHitCPtr->detid().rawId()) vetoHit = true;
        }
      }

      if (vetoHit ){
        depHcal.addCandEnergy(et);
      } else {
        depHcal.addDeposit(reco::isodeposit::Direction(hHitPos.eta(), hHitPos.phi()), et);      
      }
    }

    //HOcal
    std::vector<const HORecHit*>::const_iterator hoHitCI = mInfo.hoRecHits.begin();
    for (; hoHitCI != mInfo.hoRecHits.end(); ++hoHitCI){
      const HORecHit* hoHitCPtr = *hoHitCI;
      GlobalPoint hoHitPos = caloGeom->getPosition(hoHitCPtr->detid());
      double deltar0 = reco::deltaR(muon, hoHitPos);
      double cosTheta = 1./cosh(hoHitPos.eta());
      double energy = hoHitCPtr->energy();
      double et = energy*cosTheta;
      if (deltar0 > theDR_Max 
          || ! (et > theThreshold_HO && energy > 3*noiseRecHit(hoHitCPtr->detid()))) continue;

      bool vetoHit = false;
      double deltar = reco::deltaR(mInfo.trkGlobPosAtHO, hoHitPos);
      if (deltar < theDR_Veto_HO ){
        LogDebug("RecoMuon|CaloExtractorByAssociator")
          << " >>> Veto HO hit: Calo deltaR= " << deltar;
        LogDebug("RecoMuon|CaloExtractorByAssociator")
          << " >>> Calo eta phi ethcal: " << hoHitPos.eta() << " " << hoHitPos.phi() << " " << et;
        vetoHit = true;
      }
      if (! vetoHit){
        for (unsigned int iH = 0; iH< mInfo.crossedHOIds.size() && ! vetoHit; ++iH){
          if (mInfo.crossedHOIds[iH].rawId() == hoHitCPtr->detid().rawId()) vetoHit = true;
        }
      }

      if (vetoHit ){
        depHOcal.addCandEnergy(et);
      } else {
        depHOcal.addDeposit(reco::isodeposit::Direction(hoHitPos.eta(), hoHitPos.phi()), et);           
      }
    }


  } else {
    std::vector<const CaloTower*>::const_iterator calCI = mInfo.towers.begin();
    for (; calCI != mInfo.towers.end(); ++calCI){
      const CaloTower* calCPtr = *calCI;
      double deltar0 = reco::deltaR(muon,*calCPtr);
      if (deltar0>theDR_Max) continue;
    
      //even more copy-pasting .. need to refactor
      double etecal = calCPtr->emEt();
      double eecal = calCPtr->emEnergy();
      bool doEcal = etecal>theThreshold_E && eecal>3*noiseEcal(*calCPtr);
      double ethcal = calCPtr->hadEt();
      double ehcal = calCPtr->hadEnergy();
      bool doHcal = ethcal>theThreshold_H && ehcal>3*noiseHcal(*calCPtr);
      double ethocal = calCPtr->outerEt();
      double ehocal = calCPtr->outerEnergy();
      bool doHOcal = ethocal>theThreshold_HO && ehocal>3*noiseHOcal(*calCPtr);
      if ((!doEcal) && (!doHcal) && (!doHcal)) continue;
    
      bool vetoTowerEcal = false;
      double deltarEcal = reco::deltaR(mInfo.trkGlobPosAtEcal, *calCPtr);
      if (deltarEcal < theDR_Veto_E ){
        LogDebug("RecoMuon|CaloExtractorByAssociator")
          << " >>> Veto ecal tower: Calo deltaR= " << deltarEcal;
        LogDebug("RecoMuon|CaloExtractorByAssociator")
          << " >>> Calo eta phi ethcal: " << calCPtr->eta() << " " << calCPtr->phi() << " " << ethcal;
        vetoTowerEcal = true;
      }
      bool vetoTowerHcal = false;
      double deltarHcal = reco::deltaR(mInfo.trkGlobPosAtHcal, *calCPtr);
      if (deltarHcal < theDR_Veto_H ){
        LogDebug("RecoMuon|CaloExtractorByAssociator")
          << " >>> Veto hcal tower: Calo deltaR= " << deltarHcal;
        LogDebug("RecoMuon|CaloExtractorByAssociator")
          << " >>> Calo eta phi ethcal: " << calCPtr->eta() << " " << calCPtr->phi() << " " << ethcal;
        vetoTowerHcal = true;
      }
      bool vetoTowerHOCal = false;
      double deltarHOcal = reco::deltaR(mInfo.trkGlobPosAtHO, *calCPtr);
      if (deltarHOcal < theDR_Veto_HO ){
        LogDebug("RecoMuon|CaloExtractorByAssociator")
          << " >>> Veto HO tower: Calo deltaR= " << deltarHOcal;
        LogDebug("RecoMuon|CaloExtractorByAssociator")
          << " >>> Calo eta phi ethcal: " << calCPtr->eta() << " " << calCPtr->phi() << " " << ethcal;
        vetoTowerHOCal = true;
      }

      if (! (vetoTowerHOCal && vetoTowerHcal &&  vetoTowerEcal )){
        for (unsigned int iH = 0; iH< mInfo.crossedTowerIds.size(); ++iH){
          if (mInfo.crossedTowerIds[iH].rawId() == calCPtr->id().rawId()){
            vetoTowerEcal = true;
            vetoTowerHcal = true;
            vetoTowerHOCal = true;
            break;
          }
        }
      }

      reco::isodeposit::Direction towerDir(calCPtr->eta(), calCPtr->phi());
      if (doEcal){
        if (vetoTowerEcal) depEcal.addCandEnergy(etecal);
        else depEcal.addDeposit(towerDir, etecal);
      }
      if (doHcal){
        if (vetoTowerHcal) depHcal.addCandEnergy(ethcal);
        else depHcal.addDeposit(towerDir, ethcal);
      }
      if (doHOcal){
        if (vetoTowerHOCal) depHOcal.addCandEnergy(ethocal);
        else depHOcal.addDeposit(towerDir, ethocal);
      }
    }
  }

  std::vector<IsoDeposit> resultDeps;    
  resultDeps.push_back(depEcal);
  resultDeps.push_back(depHcal);
  resultDeps.push_back(depHOcal);

  return resultDeps;

}

void CaloExtractorByAssociator::fillVetos	(	const edm::Event &	ev,
		const edm::EventSetup &	evSetup,
		const reco::TrackCollection &	tracks
	)		`[virtual]`

allows to set extra vetoes (in addition to the muon) -- no-op at this point

Implements reco::isodeposit::IsoDepositExtractor.

Definition at line 75 of file CaloExtractorByAssociator.cc.

{
//   LogWarning("CaloExtractorByAssociator")
//     <<"fillVetos does nothing now: IsoDeposit provides enough functionality\n"
//     <<"to remove a deposit at/around given (eta, phi)";

}

double CaloExtractorByAssociator::noiseEcal ( const CaloTower & tower ) const [private]

Determine noise for HCAL and ECAL (take some defaults for the time being)

Definition at line 358 of file CaloExtractorByAssociator.cc.

References eta(), reco::LeafCandidate::eta(), theNoiseTow_EB, and theNoiseTow_EE.

Referenced by deposits().

                                                                        {
      double noise = theNoiseTow_EB;
      double eta = tower.eta();
      if (fabs(eta)>1.479) noise = theNoiseTow_EE;
      return noise;
}

double CaloExtractorByAssociator::noiseHcal ( const CaloTower & tower ) const [private]

Definition at line 365 of file CaloExtractorByAssociator.cc.

References reco::LeafCandidate::eta(), theNoise_HB, and theNoise_HE.

Referenced by deposits().

                                                                        {
  double noise = fabs(tower.eta())> 1.479 ? theNoise_HE : theNoise_HB;      
  return noise;
}

double CaloExtractorByAssociator::noiseHOcal ( const CaloTower & tower ) const [private]

Definition at line 370 of file CaloExtractorByAssociator.cc.

References theNoise_HO.

Referenced by deposits().

                                                                         {
      double noise = theNoise_HO;
      return noise;
}

double CaloExtractorByAssociator::noiseRecHit ( const DetId & detId ) const [private]

Definition at line 376 of file CaloExtractorByAssociator.cc.

References DetId::det(), DetId::Ecal, EcalBarrel, EcalEndcap, DetId::Hcal, HcalBarrel, HcalEndcap, HcalOuter, DetId::subdetId(), theNoise_EB, theNoise_EE, theNoise_HB, theNoise_HE, and theNoise_HO.

Referenced by deposits().

                                                                      {
  double  noise = 100;
  DetId::Detector det = detId.det();
  if (det == DetId::Ecal){
    EcalSubdetector subDet = (EcalSubdetector)(detId.subdetId());
    if (subDet == EcalBarrel){
      noise = theNoise_EB;
    } else if (subDet == EcalEndcap){
      noise = theNoise_EE;
    }
  } else if (det == DetId::Hcal){
    HcalSubdetector subDet = (HcalSubdetector)(detId.subdetId());
    if (subDet == HcalBarrel){
      noise = theNoise_HB;
    } else if (subDet == HcalEndcap){
      noise = theNoise_HE;      
    } else if (subDet == HcalOuter){
      noise = theNoise_HO;
    }
  }
  return noise;
}