Merge pull request #4418 from lgray/HGC-Clustering

Simple HGC Clustering and HGC GsfElectrons
cms-sw · Jul 1, 2014 · e5af628 · e5af628
2 parents ddbf7cc + 24a8f96
commit e5af628
Show file tree

Hide file tree

Showing 24 changed files with 1,249 additions and 39 deletions.
diff --git a/Fireworks/Calo/plugins/FWCaloClusterProxyBuilder.cc b/Fireworks/Calo/plugins/FWCaloClusterProxyBuilder.cc
@@ -27,7 +27,7 @@ FWCaloClusterProxyBuilder::build( const reco::CaloCluster& iData, unsigned int i
 
    TEveBoxSet* boxset = new TEveBoxSet();
    boxset->Reset(TEveBoxSet::kBT_FreeBox, true, 64);
-   boxset->UseSingleColor();
+   //boxset->UseSingleColor();
    boxset->SetPickable(1);
 
    for( std::vector<std::pair<DetId, float> >::iterator it = clusterDetIds.begin(), itEnd = clusterDetIds.end();
@@ -40,6 +40,7 @@ FWCaloClusterProxyBuilder::build( const reco::CaloCluster& iData, unsigned int i
       std::vector<float> pnts(24);    
       fireworks::energyTower3DCorners(corners, (*it).second, pnts);
       boxset->AddBox( &pnts[0]);
+      boxset->DigitColor( (4*clusterDetIds.size()+3*iIndex)%50 + 50, 50);     
    }
 
    boxset->RefitPlex();

diff --git a/RecoEcal/EgammaClusterAlgos/src/PFECALSuperClusterAlgo.cc b/RecoEcal/EgammaClusterAlgos/src/PFECALSuperClusterAlgo.cc
@@ -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;
     }
@@ -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;
   }
@@ -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;
     }

diff --git a/RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc b/RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc
@@ -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"
 
@@ -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";
 
@@ -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)
@@ -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 !!!!" ; }
 

diff --git a/RecoEgamma/EgammaElectronProducers/plugins/ElectronSeedProducer.cc b/RecoEgamma/EgammaElectronProducers/plugins/ElectronSeedProducer.cc
@@ -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"
@@ -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;
          if (HoeVeto)
           {
            sclRefs.push_back(edm::Ref<reco::SuperClusterCollection>(superClusters,i)) ;

diff --git a/RecoParticleFlow/PFClusterProducer/interface/InitialClusteringStepBase.h b/RecoParticleFlow/PFClusterProducer/interface/InitialClusteringStepBase.h
@@ -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");

diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFHGCalRecHitCreator.h b/RecoParticleFlow/PFClusterProducer/interface/PFHGCalRecHitCreator.h
@@ -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 );
 

diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterHGCEE_cfi.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterHGCEE_cfi.py
@@ -9,22 +9,71 @@
     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
+_manqiArborClusterizer_HGCEE = cms.PSet(
+    algoName = cms.string("SimpleArborClusterizer"), 
+    # use basic pad sizes in HGCEE
+    cellSize = cms.double(10.0),
+    layerThickness = cms.double(16.0),
+    thresholdsByDetector = cms.VPSet( )
+)
+
 _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 = _manqiArborClusterizer_HGCEE,
+    pfClusterBuilder = cms.PSet( ), #_arborClusterizer_HGCEE,
+    positionReCalc = cms.PSet( ), #_simplePosCalcHGCEE,
+    energyCorrector = _HGCEE_ElectronEnergy
 )
 
diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterHGCHEB_cfi.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterHGCHEB_cfi.py
@@ -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()

diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterHGCHEF_cfi.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterHGCHEF_cfi.py
@@ -10,19 +10,30 @@
 )
 
 #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( )
+)
+
+#the real arbor clusterizer
+_manqiArborClusterizer_HGCHEF = cms.PSet(
+    algoName = cms.string("SimpleArborClusterizer"), 
+    # use basic pad sizes in HGCEE
+    cellSize = cms.double(10.0),
+    layerThickness = cms.double(55.0),
+    thresholdsByDetector = cms.VPSet( )
 )
 
 particleFlowClusterHGCHEF = cms.EDProducer(
     "PFClusterProducer",
     recHitsSource = cms.InputTag("particleFlowRecHitHGCHEF"),
     recHitCleaners = cms.VPSet(),
     seedFinder = _noseeds_HGCHEF,
-    initialClusteringStep = _arborClusterizer_HGCHEF,
+    initialClusteringStep = _manqiArborClusterizer_HGCHEF,
     pfClusterBuilder = cms.PSet( ),
     positionReCalc = cms.PSet(),
     energyCorrector = cms.PSet()

diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHGCEE_cfi.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHGCEE_cfi.py
@@ -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)
                   ),
                 )
            )          

diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHGCHEB_cfi.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHGCHEB_cfi.py
@@ -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)
                   ),
                 )
            )          

diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHGCHEF_cfi.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHGCHEF_cfi.py
@@ -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)
                   ),
                 )
            )          

diff --git a/RecoParticleFlow/PFClusterProducer/src/Arbor.cc b/RecoParticleFlow/PFClusterProducer/src/Arbor.cc
@@ -646,7 +646,7 @@ std::vector< std::vector<int> > Arbor( std::vector<TVector3> inputHits, float Ce
 	BranchBuilding();	
 	BushMerging();
 
-	return LengthSortBranchCollection;
+	return Trees;
 
 }
 }