Merge pull request #7530 from sethzenz/topic-configure-hgcal-resolution

Calibration, linking, resolution parameterization for PF for HGCal

DataFormats/ParticleFlowReco/interface/PFCluster.h

@@ -12,6 +12,9 @@
 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
 #include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
 
+#include "DataFormats/TrackReco/interface/Track.h"
+#include "DataFormats/TrackReco/interface/TrackFwd.h"
+
 #include <iostream>
 #include <vector>
 #include <algorithm>
@@ -87,6 +90,10 @@ namespace reco {
     void          setEmEnergy(double em)   { emEnergy_ = em; }
     void          setHadEnergy(double had) { hadEnergy_ = had; }
 
+    //track associated to cluster (HGC)
+    void setTrack(const reco::TrackRef& tk) { track_ = tk; }
+    const reco::TrackRef& track() const { return track_; }
+
     /// cluster time
     double        time() const {return time_;}
     double        depth() const {return depth_;}
@@ -193,6 +200,9 @@ namespace reco {
     //Lindsey: add em/had energies
     double emEnergy_,hadEnergy_;
 
+    //Lindsey: add track reference when doing HGC clustering
+    reco::TrackRef track_;
+
     ///Michalis :Add timing information
     double time_;
     double depth_;

DataFormats/ParticleFlowReco/src/classes_def_2.xml

@@ -1,7 +1,8 @@
 <lcgdict>
 <selection>
 
-  <class name="reco::PFCluster" ClassVersion="15">
+  <class name="reco::PFCluster" ClassVersion="16">
+    <version ClassVersion="16" checksum="1577574542"/>
     <version ClassVersion="15" checksum="1680162253"/>
     <version ClassVersion="14" checksum="907905494"/>
    <version ClassVersion="13" checksum="1305770681"/>

RecoEcal/EgammaClusterAlgos/src/HGCALShowerBasedEmIdentification.cc

@@ -53,7 +53,7 @@ HGCALShowerBasedEmIdentification::HGCALShowerBasedEmIdentification (bool withPil
   corrlnalphalnt1_ = -0.0232;
 
   // cut values, to be moved as configurable parameters
-  cutStartPosition_ = 322.5;
+  cutStartPosition_ = 325.5;
   cutSigmaetaeta_ = 0.0055;
   if (withPileup_) cutSigmaetaeta_ = 0.00480;
   cutHoverem_ = 0.003;

RecoParticleFlow/PFClusterProducer/src/HGCClusterizer.cc

@@ -721,6 +721,7 @@ trackAssistedClustering(const edm::Handle<reco::PFRecHitCollection>& hits_handle
     std::unordered_map<unsigned,unsigned> hits_of_cluster_on_track;
     std::unordered_map<unsigned,bool> clusters_in_track;
     reco::PFCluster temp;
+    temp.setTrack(reco::TrackRef(_tracks,i));
     const reco::Track& tk = tracks[i];
     //std::cout << "got track: " << tk.pt() << ' ' << tk.eta() << ' ' << tk.phi() << std::endl;
     const TrajectoryStateOnSurface myTSOS = trajectoryStateTransform::outerStateOnSurface(tk, *(_tkGeom.product()),_bField.product());

RecoParticleFlow/PFProducer/interface/PFAlgo.h

@@ -67,6 +67,9 @@ class PFAlgo {
                      double nSigmaHCAL, 
                      const boost::shared_ptr<PFEnergyCalibration>& calibration,
 		     const boost::shared_ptr<PFEnergyCalibrationHF>& thepfEnergyCalibrationHF);
+
+  void setHGCalParameters(double HGCalResolutionConst, 
+			  double HGCalResolutionStoch);
 
   void setCandConnectorParameters( const edm::ParameterSet& iCfgCandConnector ){
     connector_.setParameters(iCfgCandConnector);
@@ -249,7 +252,8 @@ class PFAlgo {
 
   /// todo: use PFClusterTools for this
   double neutralHadronEnergyResolution( double clusterEnergy,
-					double clusterEta ) const;
+					double clusterEta,
+					reco::PFBlockElement::Type the_cluster_type ) const;
 
 
   double nSigmaHCAL( double clusterEnergy, 
@@ -309,6 +313,10 @@ class PFAlgo {
   int                algo_;
   bool               debug_;
 
+  // HGCal parameters for clustering configuration
+  double HGCalResolutionConst_;
+  double HGCalResolutionStoch_;
+
   /// Variables for PFElectrons
   std::string mvaWeightFileEleID_;
   std::vector<double> setchi2Values_;

RecoParticleFlow/PFProducer/plugins/PFProducer.cc

@@ -222,6 +222,11 @@ PFProducer::PFProducer(const edm::ParameterSet& iConfig) {
 
   int algoType 
     = iConfig.getParameter<unsigned>("algoType");
+
+  double HGCalResolutionConst
+    = iConfig.getParameter<double>("HGCalResolutionConst");
+  double HGCalResolutionStoch
+    = iConfig.getParameter<double>("HGCalResolutionStoch");
 
   switch(algoType) {
   case 0:
@@ -236,6 +241,8 @@ PFProducer::PFProducer(const edm::ParameterSet& iConfig) {
 			  calibration,
 			  thepfEnergyCalibrationHF);
 
+  pfAlgo_->setHGCalParameters(HGCalResolutionConst, HGCalResolutionStoch);
+
   //PFElectrons: call the method setpfeleparameters
   pfAlgo_->setPFEleParameters(mvaEleCut,
 			      path_mvaWeightFileEleID,

RecoParticleFlow/PFProducer/plugins/importers/ECALClusterImporter.cc

@@ -133,10 +133,14 @@ importToBlock( const edm::Event& e,
 
     // ID charged hadrons by those that pass 
     // shower start cut but fail width and length
+    /*
     if( _emPreID->cutStartPosition(*tempref) && 
 	!( _emPreID->cutSigmaetaeta(*tempref) && 
 	   _emPreID->cutLengthCompatibility(*tempref) )  && 
-	the_type == reco::PFBlockElement::HGC_ECAL ) {    
+	the_type == reco::PFBlockElement::HGC_ECAL ) { 
+    */
+    if( !_emPreID->isEm(*tempref) && 
+	the_type == reco::PFBlockElement::HGC_ECAL) {
       the_type = reco::PFBlockElement::HGC_HCALF;
     }
 

RecoParticleFlow/PFProducer/plugins/linkers/TrackAndHGCEELinker.cc

@@ -62,6 +62,11 @@ testLink( const reco::PFBlockElement* elem1,
   const reco::PFTrajectoryPoint& tkAtHGCEE =
     trackref->extrapolatedPoint( HGCEntrance );
   const reco::PFCluster::REPPoint& tkreppos = tkAtHGCEE.positionREP();
+
+  if( tkelem->trackRef() == clusterref->track() ) {
+    dist = 1e-3;
+    return dist;
+  }
 
   // Check if the linking has been done using the KDTree algo
   // Glowinski & Gouzevitch

RecoParticleFlow/PFProducer/plugins/linkers/TrackAndHGCHEBLinker.cc

@@ -63,6 +63,11 @@ testLink( const reco::PFBlockElement* elem1,
     trackref->extrapolatedPoint( HGCEntrance );
    const reco::PFCluster::REPPoint& tkreppos = tkAtHGCHEB.positionREP();
 
+   if( tkelem->trackRef() == clusterref->track() ) {
+    dist = 1e-3;
+    return dist;
+   }
+
   // Check if the linking has been done using the KDTree algo
   // Glowinski & Gouzevitch
   if ( _useKDTree && tkelem->isMultilinksValide() ) { //KDTree Algo    

RecoParticleFlow/PFProducer/plugins/linkers/TrackAndHGCHEFLinker.cc

@@ -63,6 +63,11 @@ testLink( const reco::PFBlockElement* elem1,
     trackref->extrapolatedPoint( HGCEntrance );
   const reco::PFCluster::REPPoint& tkreppos = tkAtHGCHEF.positionREP();
 
+  if( tkelem->trackRef() == clusterref->track() ) {
+    dist = 1e-3;
+    return dist;
+  }
+
   // Check if the linking has been done using the KDTree algo
   // Glowinski & Gouzevitch
   if ( _useKDTree && tkelem->isMultilinksValide() ) { //KDTree Algo    

RecoParticleFlow/PFProducer/python/particleFlow_cfi.py

@@ -181,7 +181,7 @@
     calibHF_a_EMonly  = cms.vdouble(0.96945,0.96701,0.76309,0.82268,0.87583,0.89718,0.98674,1.4681,1.4580,1.4580),
     calibHF_b_HADonly = cms.vdouble(1.27541,0.85361,0.86333,0.89091,0.94348,0.94348,0.94370,1.0034,1.0444,1.0444),
     calibHF_a_EMHAD   = cms.vdouble(1.42215,1.00496,0.68961,0.81656,0.98504,0.98504,1.00802,1.0593,1.4576,1.4576),
-    calibHF_b_EMHAD   = cms.vdouble(1.27541,0.85361,0.86333,0.89091,0.94348,0.94348,0.94370,1.0034,1.0444,1.0444)
+    calibHF_b_EMHAD   = cms.vdouble(1.27541,0.85361,0.86333,0.89091,0.94348,0.94348,0.94370,1.0034,1.0444,1.0444),
 
 #    toRead = cms.untracked.vstring("PFfa_BARREL",
 #                                   "PFfa_ENDCAP",
@@ -194,6 +194,9 @@
 #                                   "PFfbEta_BARREL",
 #                                   "PFfbEta_ENDCAP") # same strings as fType
 
+    # HGCal neutral hadron resolution
+    HGCalResolutionConst = cms.double(0.15),
+    HGCalResolutionStoch = cms.double(0.92)
 )
 
 

RecoParticleFlow/PFProducer/src/PFAlgo.cc

@@ -148,6 +148,13 @@ PFAlgo::setParameters(double nSigmaECAL,
 
 }
 
+void 
+PFAlgo::setHGCalParameters(double HGCalResolutionConst, 
+			   double HGCalResolutionStoch) {
+  HGCalResolutionConst_ = HGCalResolutionConst;
+  HGCalResolutionStoch_ = HGCalResolutionStoch;
+}
+
 
 PFMuonAlgo* PFAlgo::getPFMuonAlgo() {
   return pfmu_;
@@ -982,15 +989,17 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref,
 
 	double deficit = trackMomentum; 
 	double resol = neutralHadronEnergyResolution(trackMomentum,
-						     elements[iTrack].trackRef()->eta());
+						     elements[iTrack].trackRef()->eta(),
+						     PFBlockElement::NONE); // logically we might not have a connected cluster here
 	resol *= trackMomentum;
 
 	if ( !ecalElems.empty() ) { 
 	  unsigned thisEcal = ecalElems.begin()->second;
 	  reco::PFClusterRef clusterRef = elements[thisEcal].clusterRef();
 	  deficit -= clusterRef->energy();
 	  resol = neutralHadronEnergyResolution(trackMomentum,
-						clusterRef->positionREP().Eta());
+						clusterRef->positionREP().Eta(),
+						elements[thisEcal].type());
 	  resol *= trackMomentum;
 	}
 
@@ -1135,7 +1144,8 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref,
 	if ( !thatIsAMuon && trackRef->ptError() > ptError_) { 
 	  double deficit = trackMomentum + trackRef->p() - clusterRef->energy();
 	  double resol = nSigmaTRACK_*neutralHadronEnergyResolution(trackMomentum+trackRef->p(),
-								    clusterRef->positionREP().Eta());
+								    clusterRef->positionREP().Eta(),
+								    the_ecal_type);
 	  resol *= (trackMomentum+trackRef->p());
 	  if ( deficit > nSigmaTRACK_*resol ) { 
  	    rejectFake = true;
@@ -1360,7 +1370,7 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref,
 	*/
 
 	// Add a photon is the energy excess is large enough
-	double resol = neutralHadronEnergyResolution(trackMomentum,pivotalRef->positionREP().Eta());
+	double resol = neutralHadronEnergyResolution(trackMomentum,pivotalRef->positionREP().Eta(),the_ecal_type);
 	resol *= trackMomentum;
 	if ( neutralEnergy > std::max(0.5,nSigmaECAL_*resol) ) {
 	  neutralEnergy /= slopeEcal;
@@ -2101,7 +2111,7 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref,
     hadronDirection = hadronAtECAL.Unit();
 
     // Determine the expected calo resolution from the total charged momentum
-    double Caloresolution = neutralHadronEnergyResolution( totalChargedMomentum, hclusterref->positionREP().Eta());    
+    double Caloresolution = neutralHadronEnergyResolution( totalChargedMomentum, hclusterref->positionREP().Eta(), the_hcal_type);    
     Caloresolution *= totalChargedMomentum;
     // Account for muons
     Caloresolution = std::sqrt(Caloresolution*Caloresolution + muonHCALError + muonECALError);
@@ -2287,7 +2297,7 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref,
 	  //if ( totalChargedMomentum < caloEnergy ) break;	
 	}
 	// New calo resolution.
-	Caloresolution = neutralHadronEnergyResolution( totalChargedMomentum, hclusterref->positionREP().Eta());    
+	Caloresolution = neutralHadronEnergyResolution( totalChargedMomentum, hclusterref->positionREP().Eta(), the_hcal_type);    
 	Caloresolution *= totalChargedMomentum;
 	Caloresolution = std::sqrt(Caloresolution*Caloresolution + muonHCALError + muonECALError);
       }
@@ -2369,7 +2379,7 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref,
     }
 
     // New determination of the calo and track resolution avec track deletion/rescaling.
-    Caloresolution = neutralHadronEnergyResolution( totalChargedMomentum, hclusterref->positionREP().Eta());    
+    Caloresolution = neutralHadronEnergyResolution( totalChargedMomentum, hclusterref->positionREP().Eta(), the_hcal_type);    
     Caloresolution *= totalChargedMomentum;
     Caloresolution = std::sqrt(Caloresolution*Caloresolution + muonHCALError + muonECALError);
 
@@ -2417,7 +2427,7 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref,
     }
 
     // New determination of the calo and track resolution avec track deletion/rescaling.
-    Caloresolution = neutralHadronEnergyResolution( totalChargedMomentum, hclusterref->positionREP().Eta());    
+    Caloresolution = neutralHadronEnergyResolution( totalChargedMomentum, hclusterref->positionREP().Eta(), the_hcal_type);    
     Caloresolution *= totalChargedMomentum;
     Caloresolution = std::sqrt(Caloresolution*Caloresolution + muonHCALError + muonECALError);
 
@@ -3396,7 +3406,16 @@ PFAlgo::reconstructCluster(const reco::PFCluster& cluster,
 //GMA need the followign two for HO also
 
 double 
-PFAlgo::neutralHadronEnergyResolution(double clusterEnergyHCAL, double eta) const {
+PFAlgo::neutralHadronEnergyResolution(double clusterEnergyHCAL, double eta, PFBlockElement::Type the_cluster_type) const {
+
+  // Use configurable HGCal resolution
+  if (the_cluster_type == PFBlockElement::HGC_HCALF ||
+      the_cluster_type == PFBlockElement::HGC_HCALB ||
+      the_cluster_type == PFBlockElement::HGC_ECAL || 
+      (the_cluster_type == PFBlockElement::NONE && fabs(eta) > 1.48) ) { 
+    if (debug_) cout << "\t\t\tUsing configured HGCal resolution const=" << HGCalResolutionConst_ << " stoch=" << HGCalResolutionStoch_ << endl;
+    return sqrt( HGCalResolutionStoch_*HGCalResolutionStoch_/clusterEnergyHCAL + HGCalResolutionConst_*HGCalResolutionConst_ );
+  }
 
   // Add a protection