[HGCal] first version of the HGCalValidator

DataFormats/DetId/src/classes_def.xml

@@ -9,4 +9,8 @@
   <class name="std::pair<DetId,float>"/>               
   <class name="std::map<DetId, std::pair<unsigned int, unsigned int> >"/>
   <class name="std::map<DetId, std::pair<unsigned long, unsigned long> >"/>
+  <class name="std::pair<const DetId,float>"/>
+  <class name="std::map<DetId,float>"/>
+  <class name="edm::Wrapper<std::map<DetId,float> >"/>
+
 </lcgdict>

DataFormats/EgammaReco/src/classes_def.xml

@@ -156,6 +156,5 @@
   <class name="edm::Wrapper<edm::RefVector<edm::AssociationMap<edm::OneToOne<vector<reco::CaloCluster>,vector<reco::ClusterShape>,unsigned int> >,edm::helpers::KeyVal<edm::Ref<vector<reco::CaloCluster>,reco::CaloCluster,edm::refhelper::FindUsingAdvance<vector<reco::CaloCluster>,reco::CaloCluster> >,edm::Ref<vector<reco::ClusterShape>,reco::ClusterShape,edm::refhelper::FindUsingAdvance<vector<reco::ClusterShape>,reco::ClusterShape> > >,edm::AssociationMap<edm::OneToOne<vector<reco::CaloCluster>,vector<reco::ClusterShape>,unsigned int> >::Find> >" />
   <class name="edm::Wrapper<edm::RefVector<edm::AssociationMap<edm::OneToOne<vector<reco::SuperCluster>,vector<reco::HFEMClusterShape>,unsigned int> >,edm::helpers::KeyVal<edm::Ref<vector<reco::SuperCluster>,reco::SuperCluster,edm::refhelper::FindUsingAdvance<vector<reco::SuperCluster>,reco::SuperCluster> >,edm::Ref<vector<reco::HFEMClusterShape>,reco::HFEMClusterShape,edm::refhelper::FindUsingAdvance<vector<reco::HFEMClusterShape>,reco::HFEMClusterShape> > >,edm::AssociationMap<edm::OneToOne<vector<reco::SuperCluster>,vector<reco::HFEMClusterShape>,unsigned int> >::Find> >" />
 
-
 </lcgdict>
 

RecoLocalCalo/HGCalRecAlgos/interface/HGCalImagingAlgo.h

@@ -24,6 +24,8 @@
 
 #include "KDTreeLinkerAlgoT.h"
 
+//Density collection
+typedef std::map< DetId, float > Density;
 
 template <typename T>
 std::vector<size_t> sorted_indices(const std::vector<T> &v) {
@@ -41,6 +43,19 @@ std::vector<size_t> sorted_indices(const std::vector<T> &v) {
         return idx;
 }
 
+template <typename T>
+size_t max_index(const std::vector<T> &v) {
+
+        // initialize original index locations
+        std::vector<size_t> idx(v.size(),0);
+        std::iota (std::begin(idx), std::end(idx), 0);
+
+        // take the max index based on comparing values in v
+        auto maxidx = std::max_element(idx.begin(), idx.end(), [&v](size_t i1, size_t i2) {return v[i1].data.rho < v[i2].data.rho;});
+
+        return (*maxidx);
+}
+
 class HGCalImagingAlgo
 {
 
@@ -161,6 +176,9 @@ void reset(){
 }
 void computeThreshold();
 
+//getDensity
+ Density getDensity();
+
 /// point in the space
 typedef math::XYZPoint Point;
 
@@ -195,6 +213,9 @@ hgcal::RecHitTools rhtools_;
 // The algo id
 reco::CaloCluster::AlgoId algoId_;
 
+// For keeping the density per hit
+ Density density_;
+
 // various parameters used for calculating the noise levels for a given sensor (and whether to use them)
 bool dependSensor_;
 std::vector<double> dEdXweights_;
@@ -204,7 +225,7 @@ double fcPerEle_;
 std::vector<double> nonAgedNoises_;
 double noiseMip_;
 std::vector<std::vector<double> > thresholds_;
-std::vector<std::vector<double> > v_sigmaNoise_;
+std::vector<std::vector<double> > sigmaNoise_;
 
 // The verbosity level
 VerbosityLevel verbosity_;
@@ -297,6 +318,9 @@ double calculateDistanceToHigher(std::vector<KDNode> &) const;
 int findAndAssignClusters(std::vector<KDNode> &, KDTree &, double, KDTreeBox &, const unsigned int, std::vector<std::vector<KDNode> >&) const;
 math::XYZPoint calculatePosition(std::vector<KDNode> &) const;
 
+//For keeping the density information
+ void setDensity(const std::vector<KDNode> &nd);
+
 // attempt to find subclusters within a given set of hexels
 std::vector<unsigned> findLocalMaximaInCluster(const std::vector<KDNode>&);
 math::XYZPoint calculatePositionWithFraction(const std::vector<KDNode>&, const std::vector<double>&);

RecoLocalCalo/HGCalRecAlgos/src/HGCalImagingAlgo.cc

@@ -37,7 +37,7 @@ void HGCalImagingAlgo::populate(const HGCRecHitCollection &hits) {
       if (thickness_index == -1)
         thickness_index = 3;
       double storedThreshold = thresholds_[layer - 1][thickness_index];
-      sigmaNoise = v_sigmaNoise_[layer - 1][thickness_index];
+      sigmaNoise = sigmaNoise_[layer - 1][thickness_index];
 
       if (hgrh.energy() < storedThreshold)
         continue; // this sets the ZS threshold at ecut times the sigma noise
@@ -98,6 +98,8 @@ void HGCalImagingAlgo::makeClusters() {
       double maxdensity = calculateLocalDensity(
           points_[i], hit_kdtree, actualLayer); // also stores rho (energy
                                                // density) for each point (node)
+      //Now that we have the density per point we can store it
+      setDensity(points_[i]);
       // calculate distance to nearest point with higher density storing
       // distance (delta) and point's index
       calculateDistanceToHigher(points_[i]);
@@ -115,6 +117,9 @@ std::vector<reco::BasicCluster> HGCalImagingAlgo::getClusters(bool doSharing) {
     for (unsigned int i = 0; i < clsOnLayer.size(); ++i) {
       double energy = 0;
       Point position;
+
+      //Will save the maximum density hit of the cluster
+      size_t rsmax = max_index(clsOnLayer[i]);
 
       if (doSharing) {
 
@@ -158,7 +163,8 @@ std::vector<reco::BasicCluster> HGCalImagingAlgo::getClusters(bool doSharing) {
           }
           clusters_v_.emplace_back(energy, position, caloID, thisCluster,
                                   algoId_);
-          thisCluster.clear();
+	  if (!clusters_v_.empty()){ clusters_v_.back().setSeed( clsOnLayer[i][rsmax].data.detid); }
+	  thisCluster.clear();
         }
       } else {
         position = calculatePosition(clsOnLayer[i]); // energy-weighted position
@@ -178,6 +184,7 @@ std::vector<reco::BasicCluster> HGCalImagingAlgo::getClusters(bool doSharing) {
           std::cout << "*****************************" << std::endl;
         }
         clusters_v_.emplace_back(energy, position, caloID, thisCluster, algoId_);
+	if (!clusters_v_.empty()){ clusters_v_.back().setSeed( clsOnLayer[i][rsmax].data.detid); }
         thisCluster.clear();
       }
     }
@@ -652,19 +659,33 @@ void HGCalImagingAlgo::computeThreshold() {
   const unsigned maxNumberOfThickIndices = 3;
   dummy.resize(maxNumberOfThickIndices + 1, 0); // +1 to accomodate for the Scintillators
   thresholds_.resize(maxlayer, dummy);
-  v_sigmaNoise_.resize(maxlayer, dummy);
+  sigmaNoise_.resize(maxlayer, dummy);
 
   for (unsigned ilayer = 1; ilayer <= maxlayer; ++ilayer) {
     for (unsigned ithick = 0; ithick < maxNumberOfThickIndices; ++ithick) {
       float sigmaNoise =
           0.001f * fcPerEle_ * nonAgedNoises_[ithick] * dEdXweights_[ilayer] /
           (fcPerMip_[ithick] * thicknessCorrection_[ithick]);
       thresholds_[ilayer - 1][ithick] = sigmaNoise * ecut_;
-      v_sigmaNoise_[ilayer - 1][ithick] = sigmaNoise;
+      sigmaNoise_[ilayer - 1][ithick] = sigmaNoise;
     }
     float scintillators_sigmaNoise = 0.001f * noiseMip_ * dEdXweights_[ilayer];
     thresholds_[ilayer - 1][maxNumberOfThickIndices] = ecut_ * scintillators_sigmaNoise;
-    v_sigmaNoise_[ilayer -1][maxNumberOfThickIndices] = scintillators_sigmaNoise;
+    sigmaNoise_[ilayer -1][maxNumberOfThickIndices] = scintillators_sigmaNoise;
   }
 
 }
+
+void HGCalImagingAlgo::setDensity(const std::vector<KDNode> &nd){
+
+  // for each node calculate local density rho and store it
+  for (auto &i : nd){
+    density_[ i.data.detid ] =  i.data.rho ;
+  }   // end loop nodes
+}
+
+//Density 
+Density HGCalImagingAlgo::getDensity(){
+  return density_;
+}
+

RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClusterProducer.cc

@@ -108,6 +108,8 @@ HGCalLayerClusterProducer::HGCalLayerClusterProducer(const edm::ParameterSet &ps
 
   produces<std::vector<reco::BasicCluster> >();
   produces<std::vector<reco::BasicCluster> >("sharing");
+  //density
+  produces< Density >();
   //time for layer clusters
   produces<edm::ValueMap<float> > (timeClname);
 
@@ -168,6 +170,7 @@ void HGCalLayerClusterProducer::produce(edm::Event& evt,
 
   std::unique_ptr<std::vector<reco::BasicCluster> > clusters( new std::vector<reco::BasicCluster> ),
     clusters_sharing( new std::vector<reco::BasicCluster> );
+  auto density = std::make_unique<Density>();
 
   algo->reset();
 
@@ -218,6 +221,11 @@ void HGCalLayerClusterProducer::produce(edm::Event& evt,
 
   auto clusterHandle = evt.put(std::move(clusters));
   auto clusterHandleSharing = evt.put(std::move(clusters_sharing),"sharing");
+
+  //Keep the density 
+  *density = algo->getDensity();
+  evt.put(std::move(density));
+
   edm::PtrVector<reco::BasicCluster> clusterPtrs, clusterPtrsSharing;
 
   std::vector<float> times;

Validation/Configuration/python/autoValidation.py

@@ -10,7 +10,7 @@
                    'miniAODValidation' : ['prevalidationMiniAOD','validationMiniAOD','validationHarvestingMiniAOD'],
                    'standardValidation' : ['prevalidation','validation','validationHarvesting'],
                    'standardValidationNoHLT' : ['prevalidationNoHLT','validationNoHLT','validationHarvestingNoHLT'],
-                   'HGCalValidation' : ['', 'globalValidationHGCal', ''],
+                   'HGCalValidation' : ['', 'globalValidationHGCal', 'hgcalValidatorPostProcessor'],
                    'OuterTrackerValidation' : ['', 'globalValidationOuterTracker', 'postValidationOuterTracker'],
                  }
 

Validation/Configuration/python/hgcalSimValid_cff.py

@@ -5,6 +5,10 @@
 from Validation.HGCalValidation.rechitValidation_cff    import *
 from Validation.HGCalValidation.hgcalHitValidation_cfi  import *
 
+from Validation.HGCalValidation.HGCalValidator_cfi import hgcalValidator
+
+hgcalValidatorSequence = cms.Sequence(hgcalValidator)
+
 hgcalValidation = cms.Sequence(hgcalSimHitValidationEE
                                + hgcalSimHitValidationHEF
                                + hgcalSimHitValidationHEB
@@ -14,4 +18,5 @@
                                + hgcalRecHitValidationEE
                                + hgcalRecHitValidationHEF
                                + hgcalRecHitValidationHEB
-                               + hgcalHitValidationSequence)
+                               + hgcalHitValidationSequence
+                               + hgcalValidatorSequence)

Validation/HGCalValidation/BuildFile.xml

@@ -0,0 +1,33 @@
+<use name="CalibFormats/HcalObjects"/>
+<use name="CondFormats/GeometryObjects"/>
+<use name="DataFormats/ForwardDetId"/>
+<use name="DataFormats/HcalDetId"/>
+<use name="DataFormats/HepMCCandidate"/>
+<use name="DataFormats/ParticleFlowCandidate"/>
+<use name="DetectorDescription/Core"/>
+<use name="DQMServices/Core"/>
+<use name="FWCore/Framework"/>
+<use name="FWCore/ParameterSet"/>
+<use name="FWCore/MessageLogger"/>
+<use name="FWCore/PluginManager"/>
+<use name="Geometry/HcalCommonData"/>
+<use name="Geometry/HGCalCommonData"/>
+<use name="Geometry/HGCalGeometry"/>
+<use name="Geometry/Records"/>
+<use name="SimDataFormats/CaloHit"/>
+<use name="SimDataFormats/CaloTest"/>
+<use name="SimDataFormats/CaloAnalysis"/>
+<use name="SimDataFormats/GeneratorProducts"/>
+<use name="SimDataFormats/Track"/>
+<use name="SimDataFormats/ValidationFormats"/>
+<use name="SimG4Core/Notification"/>
+<use name="SimG4Core/Watcher"/>
+<use name="SimG4CMS/Calo"/>
+<use name="RecoLocalCalo/HGCalRecAlgos"/>
+<use name="geant4core"/>
+<use name="clhep"/>
+<use name="hepmc"/>
+
+<export>
+  <lib name="ValidationHGCalValidation"/>
+</export>

Validation/HGCalValidation/data/D28.cumulative.xo

@@ -0,0 +1,52 @@
+1.0 0.434
+2.0 1.402
+3.0 2.273
+4.0 3.241
+5.0 4.112
+6.0 5.08
+7.0 5.951
+8.0 6.919
+9.0 7.79
+10.0 8.758
+11.0 9.629
+12.0 10.597
+13.0 11.468
+14.0 12.436
+15.0 13.307
+16.0 14.275
+17.0 15.146
+18.0 16.114
+19.0 16.985
+20.0 17.953
+21.0 18.824
+22.0 19.792
+23.0 20.663
+24.0 21.631
+25.0 22.502
+26.0 23.47
+27.0 24.341
+28.0 25.309
+29.0 27.793
+30.0 30.324
+31.0 32.856
+32.0 35.387
+33.0 37.919
+34.0 40.45
+35.0 42.982
+36.0 45.513
+37.0 48.045
+38.0 50.576
+39.0 53.108
+40.0 55.639
+41.0 60.072
+42.0 64.505
+43.0 68.969
+44.0 73.44
+45.0 77.912
+46.0 82.383
+47.0 86.854
+48.0 91.325
+49.0 95.797
+50.0 100.268
+51.0 104.739
+52.0 109.211