Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Simple HGC Clustering and HGC GsfElectrons #4418

Merged
merged 4 commits into from Jul 1, 2014
Merged
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
14 changes: 14 additions & 0 deletions RecoEcal/EgammaClusterAlgos/src/PFECALSuperClusterAlgo.cc
Expand Up @@ -251,6 +251,10 @@ loadAndSortPFClusters(const edm::Event &iEvent) {
_clustersEE.push_back(calib_cluster);
}
break;
case PFLayer::HGC_ECAL:
if( calib_cluster->energy() > threshPFClusterEndcap_ ) {
_clustersEE.push_back(calib_cluster);
}
default:
break;
}
Expand Down Expand Up @@ -309,6 +313,13 @@ buildSuperCluster(CalibClusterPtr& seed,
<< " in the ECAL endcap!" << std::endl;
isEE = true;
break;
case PFLayer::HGC_ECAL:
IsClusteredWithSeed.phiwidthSuperCluster_ = phiwidthSuperClusterEndcap_;
IsClusteredWithSeed.etawidthSuperCluster_ = etawidthSuperClusterEndcap_;
edm::LogInfo("PFClustering") << "Building HGC SC number "
<< superClustersEE_->size() + 1
<< " in the HGC ECAL!" << std::endl;
break;
default:
break;
}
Expand Down Expand Up @@ -491,6 +502,9 @@ buildSuperCluster(CalibClusterPtr& seed,
case PFLayer::ECAL_ENDCAP:
superClustersEE_->push_back(new_sc);
break;
case PFLayer::HGC_ECAL:
superClustersEE_->push_back(new_sc);
break;
default:
break;
}
Expand Down
7 changes: 7 additions & 0 deletions RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc
Expand Up @@ -21,6 +21,7 @@
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "DataFormats/ForwardDetId/interface/HGCEEDetId.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/GeometryVector/interface/GlobalVector.h"

Expand Down Expand Up @@ -1032,10 +1033,12 @@ void GsfElectronAlgo::setCutBasedPreselectionFlag( GsfElectron * ele, const reco
LogTrace("GsfElectronAlgo") << "HoE1 : " << ele->hcalDepth1OverEcal() << ", HoE2 : " << ele->hcalDepth2OverEcal();
double had = ele->hcalOverEcal()*ele->superCluster()->energy() ;
const reco::CaloCluster & seedCluster = *(ele->superCluster()->seed()) ;
int component = seedCluster.hitsAndFractions()[0].first.det();
int detector = seedCluster.hitsAndFractions()[0].first.subdetId() ;
bool HoEveto = false ;
if (detector==EcalBarrel && (had<cfg->maxHBarrel || (had/ele->superCluster()->energy())<cfg->maxHOverEBarrel)) HoEveto=true;
else if (detector==EcalEndcap && (had<cfg->maxHEndcaps || (had/ele->superCluster()->energy())<cfg->maxHOverEEndcaps)) HoEveto=true;
else if (component== DetId::Forward && detector==HGCEE && (had<cfg->maxHEndcaps || (had/ele->superCluster()->energy())<cfg->maxHOverEEndcaps)) HoEveto=true;
if ( !HoEveto ) return ;
LogTrace("GsfElectronAlgo") << "H/E criteria are satisfied";

Expand Down Expand Up @@ -1198,6 +1201,7 @@ void GsfElectronAlgo::createElectron()
//====================================================

reco::GsfElectron::FiducialFlags fiducialFlags ;
int component = seedXtalId.det();
int detector = seedXtalId.subdetId() ;
double feta=std::abs(electronData_->superClusterRef->position().eta()) ;
if (detector==EcalBarrel)
Expand Down Expand Up @@ -1232,6 +1236,9 @@ void GsfElectronAlgo::createElectron()
{
fiducialFlags.isEE = true;
}
else if( component == DetId::Forward && detector==HGCEE ) {
fiducialFlags.isEE = true;
}
else
{ throw cms::Exception("GsfElectronAlgo|UnknownXtalRegion")<<"createElectron(): do not know if it is a barrel or endcap seed cluster !!!!" ; }

Expand Down
Expand Up @@ -37,6 +37,8 @@
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"

#include "DataFormats/ForwardDetId/interface/HGCEEDetId.h"

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
Expand Down Expand Up @@ -224,9 +226,11 @@ void ElectronSeedProducer::filterClusters
had2 = hcalHelper_->hcalESumDepth2(scl);
had = had1+had2 ;
scle = scl.energy() ;
int component = scl.seed()->hitsAndFractions()[0].first.det() ;
int detector = scl.seed()->hitsAndFractions()[0].first.subdetId() ;
if (detector==EcalBarrel && (had<maxHBarrel_ || had/scle<maxHOverEBarrel_)) HoeVeto=true;
else if (detector==EcalEndcap && (had<maxHEndcaps_ || had/scle<maxHOverEEndcaps_)) HoeVeto=true;
else if (component==DetId::Forward && detector==HGCEE && (had<maxHEndcaps_ || had/scle<maxHOverEEndcaps_)) HoeVeto=true;

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess this "if" chain should first select appropriate det() and only then test subdetId()

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi Fedor, I agree that logic would be cleaner, but given that this code operates on super-clusters (where the data volume is significantly reduced) it makes no difference in performance.

if (HoeVeto)
{
sclRefs.push_back(edm::Ref<reco::SuperClusterCollection>(superClusters,i)) ;
Expand Down
Expand Up @@ -32,7 +32,10 @@ class InitialClusteringStepBase {
{"HCAL_BARREL2_RING1",100*(int)PFLayer::HCAL_BARREL2},
{"HCAL_ENDCAP",(int)PFLayer::HCAL_ENDCAP},
{"HF_EM",(int)PFLayer::HF_EM},
{"HF_HAD",(int)PFLayer::HF_HAD} }),
{"HF_HAD",(int)PFLayer::HF_HAD},
{"HGC_ECAL",(int)PFLayer::HGC_ECAL},
{"HGC_HCALF",(int)PFLayer::HGC_HCALF},
{"HGC_HCALB",(int)PFLayer::HGC_HCALB},}),
_algoName(conf.getParameter<std::string>("algoName")) {
const std::vector<edm::ParameterSet>& thresholds =
conf.getParameterSetVector("thresholdsByDetector");
Expand Down
Expand Up @@ -78,6 +78,8 @@ template <typename DET,PFLayer::Layer Layer,ForwardSubdetector subdet>
position.x(), position.y(), position.z(),
0, 0, 0 );

rh.setOriginalRecHit(edm::Ref<HGCRecHitCollection>(recHitHandle,i));

const HGCalGeometry::CornersVec corners( std::move( hgcGeo.getCorners( detid ) ) );
assert( corners.size() == 8 );

Expand Down
Expand Up @@ -9,22 +9,63 @@
algoName = cms.string("PassThruSeedFinder")
)

#topo clusters
#for arbor this is more a pre-clustering step to find little clusters
_arborTopoClusterizer_HGCEE = cms.PSet(
algoName = cms.string("IntraLayerClusteringAlgorithm"),
IntraLayerMaxDistance = cms.double( 19.0 ), # hit separation in mm
ShouldSplitClusterInSingleCaloHitClusters = cms.bool(False), # splitsmall clusters
MaximumSizeForClusterSplitting = cms.uint32( 3 ), #largest of small clusters to split
thresholdsByDetector = cms.VPSet( )
)
_simplePosCalcHGCEE = cms.PSet(
algoName = cms.string("SimplePositionCalc"),
minFractionInCalc = cms.double(0.0)
)

#the real arbor clusterizer
_arborClusterizer_HGCEE = cms.PSet(
algoName = cms.string("SimpleArborClusterizer"),
cellSize = cms.double(2.0),
layerThickness = cms.double(2.0),
thresholdsByDetector = cms.VPSet()
algoName = cms.string("ArborConnectorClusteringAlgorithm"),
# these are taken from the settings for Fine Granularity in ArborPFA
MaximumForwardDistanceForConnection = cms.double(60.0), #in mm
MaximumTransverseDistanceForConnection = cms.double(40.0), #in mm
AllowForwardConnectionForIsolatedObjects = cms.bool(False),
ShouldUseIsolatedObjects = cms.bool(True),
MaximumNumberOfKeptConnectors = cms.uint32(5),
OrderParameterAnglePower = cms.double(1.0),
OrderParameterDistancePower = cms.double(0.5),
minFractionToKeep = cms.double(1e-7)
)

#weights for layers from P.Silva (24 June 2014)
weight_vec = [0.42]
weight_vec.extend([1.00 for x in range(10)])
weight_vec.extend([1.61 for x in range(10)])
weight_vec.extend([2.44 for x in range(10)])

# MIP effective to 1.0/GeV (from fit to data of P. Silva)
#f(x) = a/(1-exp(-bx - c))
# x = cosh(eta)
# a = 168.0
# b = 0.6871
# c = 0.9038

_HGCEE_ElectronEnergy = cms.PSet(
algoName = cms.string("HGCEEElectronEnergyCalibrator"),
weights = cms.vdouble(weight_vec),
effMip_to_InverseGeV_a = cms.double(168.0),
effMip_to_InverseGeV_b = cms.double(0.6871),
effMip_to_InverseGeV_c = cms.double(0.9038),
MipValueInGeV = cms.double(55.1*1e-6)
)

particleFlowClusterHGCEE = cms.EDProducer(
"PFClusterProducer",
recHitsSource = cms.InputTag("particleFlowRecHitHGCEE"),
recHitCleaners = cms.VPSet(),
seedFinder = _noseeds_HGCEE,
initialClusteringStep = _arborClusterizer_HGCEE,
pfClusterBuilder = cms.PSet( ),
positionReCalc = cms.PSet(),
energyCorrector = cms.PSet()
initialClusteringStep = _arborTopoClusterizer_HGCEE,
pfClusterBuilder = cms.PSet( ), #_arborClusterizer_HGCEE,
positionReCalc = cms.PSet( ), #_simplePosCalcHGCEE,
energyCorrector = _HGCEE_ElectronEnergy
)

Expand Up @@ -10,19 +10,21 @@
)

#topo clusters
_arborClusterizer_HGCHEB = cms.PSet(
algoName = cms.string("SimpleArborClusterizer"),
cellSize = cms.double(10.0),
layerThickness = cms.double(7.0),
thresholdsByDetector = cms.VPSet()
#for arbor this is more a pre-clustering step to find little clusters
_arborTopoClusterizer_HGCHEB = cms.PSet(
algoName = cms.string("IntraLayerClusteringAlgorithm"),
IntraLayerMaxDistance = cms.double( 186.0 ), # hit separation in mm
ShouldSplitClusterInSingleCaloHitClusters = cms.bool(False), # splitsmall clusters
MaximumSizeForClusterSplitting = cms.uint32( 3 ), #largest of small clusters to split
thresholdsByDetector = cms.VPSet( )
)

particleFlowClusterHGCHEB = cms.EDProducer(
"PFClusterProducer",
recHitsSource = cms.InputTag("particleFlowRecHitHGCHEB"),
recHitCleaners = cms.VPSet(),
seedFinder = _noseeds_HGCHEB,
initialClusteringStep = _arborClusterizer_HGCHEB,
initialClusteringStep = _arborTopoClusterizer_HGCHEB,
pfClusterBuilder = cms.PSet( ),
positionReCalc = cms.PSet(),
energyCorrector = cms.PSet()
Expand Down
Expand Up @@ -10,19 +10,21 @@
)

#topo clusters
_arborClusterizer_HGCHEF = cms.PSet(
algoName = cms.string("SimpleArborClusterizer"),
cellSize = cms.double(3.0),
layerThickness = cms.double(2.5),
thresholdsByDetector = cms.VPSet()
#for arbor this is more a pre-clustering step to find little clusters
_arborTopoClusterizer_HGCHEF = cms.PSet(
algoName = cms.string("IntraLayerClusteringAlgorithm"),
IntraLayerMaxDistance = cms.double( 96.0 ), # hit separation in mm
ShouldSplitClusterInSingleCaloHitClusters = cms.bool(False), # splitsmall clusters
MaximumSizeForClusterSplitting = cms.uint32( 3 ), #largest of small clusters to split
thresholdsByDetector = cms.VPSet( )
)

particleFlowClusterHGCHEF = cms.EDProducer(
"PFClusterProducer",
recHitsSource = cms.InputTag("particleFlowRecHitHGCHEF"),
recHitCleaners = cms.VPSet(),
seedFinder = _noseeds_HGCHEF,
initialClusteringStep = _arborClusterizer_HGCHEF,
initialClusteringStep = _arborTopoClusterizer_HGCHEF,
pfClusterBuilder = cms.PSet( ),
positionReCalc = cms.PSet(),
energyCorrector = cms.PSet()
Expand Down
Expand Up @@ -13,10 +13,10 @@
qualityTests = cms.VPSet(
cms.PSet(
name = cms.string("PFRecHitQTestThresholdInMIPs"),
thresholdInMIPs = cms.double(0.50),
thresholdInMIPs = cms.double(0.544),
mipValueInkeV = cms.double(55.1),
recHitEnergyIs_keV = cms.bool(True),
recHitEnergyMultiplier = cms.double(12.0)
recHitEnergyIs_keV = cms.bool(False),
recHitEnergyMultiplier = cms.double(1.0)
),
)
)
Expand Down
Expand Up @@ -13,10 +13,10 @@
qualityTests = cms.VPSet(
cms.PSet(
name = cms.string("PFRecHitQTestThresholdInMIPs"),
thresholdInMIPs = cms.double(0.20),
thresholdInMIPs = cms.double(0.6),
mipValueInkeV = cms.double(1498.4),
recHitEnergyIs_keV = cms.bool(True),
recHitEnergyMultiplier = cms.double(12.0)
recHitEnergyIs_keV = cms.bool(False),
recHitEnergyMultiplier = cms.double(1.0)
),
)
)
Expand Down
Expand Up @@ -13,10 +13,10 @@
qualityTests = cms.VPSet(
cms.PSet(
name = cms.string("PFRecHitQTestThresholdInMIPs"),
thresholdInMIPs = cms.double(0.30),
thresholdInMIPs = cms.double(0.50),
mipValueInkeV = cms.double(85.0),
recHitEnergyIs_keV = cms.bool(True),
recHitEnergyMultiplier = cms.double(12.0)
recHitEnergyIs_keV = cms.bool(False),
recHitEnergyMultiplier = cms.double(1.0)
),
)
)
Expand Down