Merge pull request #62 from beaudett/topic_fsimevent2

Got rid of hard-coded values for reachedEE

Definitions.md

@@ -4,7 +4,7 @@
 | ------------- | ------------- | ----- | ----- |
 | event/run/lumi | - | - | general event info |
 | vtx | `g4SimHits` | `std::vector<SimVertex>` | primary vertex position |
-| genpart | `g4SimHits` | `std::vector<SimTrack>` | truth level tracks/particles |
+| genpart | `g4SimHits` | `std::vector<SimTrack>` | truth level tracks/particles and information related to their extrapolation towards HGCAL. In particular `reachedEE==2` indicates that the particles reached HGCAL while `reachedEE==1` is for barrel calorimeter and `reachedEE==0` is for the other cases |
 | simcluster | `mix:MergedCaloTruth` | `std::vector<SimCluster>` | Geant particle and its associated hits (DetIds) in the HGCal |
 | pfcluster | `particleFlowClusterHGCal` | `std::vector<reco::PFCluster>` | mapping of the SimCluster DetIds to the reconstructed hits |
 | cluster2d | `hgcalLayerClusters` | `reco::CaloClusterCollection` | reconstructed layer (2D) clusters - those that are associated to a multicluster have `cluster2d_multicluster >= 0`, which is the index of the `multiclus` in the ntuple |

HGCalAnalysis/plugins/HGCalAnalysis.cc

@@ -150,6 +150,10 @@ class HGCalAnalysis : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one:
   typedef ROOT::Math::Transform3DPJ Transform3D;
   typedef ROOT::Math::Transform3DPJ::Point Point;
 
+  // approximative geometrical values
+  static constexpr float hgcalOuterRadius_ = 160.;
+  static constexpr float hgcalInnerRadius_ = 25.;
+
   HGCalAnalysis();
   explicit HGCalAnalysis(const edm::ParameterSet &);
   ~HGCalAnalysis();
@@ -474,6 +478,7 @@ class HGCalAnalysis : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one:
   std::vector<double> dEdXWeights_;
   std::vector<double> invThicknessCorrection_;
 
+
   // and also the magnetic field
   MagneticField const *aField_;
 
@@ -1125,21 +1130,21 @@ void HGCalAnalysis::analyze(const edm::Event &iEvent, const edm::EventSetup &iSe
 
     if (std::abs(myTrack.vertex().position().z()) >= layerPositions_[0]) continue;
 
-    unsigned nlayers = 40;
+    unsigned nlayers =   recHitTools_.lastLayerFH();
     if (myTrack.noEndVertex())  // || myTrack.genpartIndex()>=0)
     {
       HGCal_helpers::coordinates propcoords;
       bool reachesHGCal = toHGCalPropagator.propagate(
           myTrack.momentum(), myTrack.vertex().position(), myTrack.charge(), propcoords);
       vtx = propcoords.toVector();
 
-      if (reachesHGCal && vtx.Rho() < 160 && vtx.Rho() > 25) {
+      if (reachesHGCal && vtx.Rho() < hgcalOuterRadius_ && vtx.Rho() > hgcalInnerRadius_) {
         reachedEE = 2;
         double dpt = 0;
 
         for (int i = 0; i < myTrack.nDaughters(); ++i) dpt += myTrack.daughter(i).momentum().pt();
         if (abs(myTrack.type()) == 11) fbrem = dpt / myTrack.momentum().pt();
-      } else if (reachesHGCal && vtx.Rho() > 160)
+      } else if (reachesHGCal && vtx.Rho() > hgcalOuterRadius_)
         reachedEE = 1;
 
       HGCal_helpers::simpleTrackPropagator indiv_particleProp(aField_);