diff --git a/.gitignore b/.gitignore deleted file mode 100644 index 89543e26bf3d1..0000000000000 --- a/.gitignore +++ /dev/null @@ -1,15 +0,0 @@ -__init__.py -*.pyc -.*.swp -.#* -#*# -*~ -*.pb - -# ignore partially applied patches and their backup -*.rej -*.orig - -# ignore files under the top level $CMSSW_BASE/src directory, but not its subdirectories -/* -!/*/ diff --git a/Configuration/ProcessModifiers/python/kf_cff.py b/Configuration/ProcessModifiers/python/kf_cff.py new file mode 100644 index 0000000000000..0d00a11cc6baa --- /dev/null +++ b/Configuration/ProcessModifiers/python/kf_cff.py @@ -0,0 +1,5 @@ +import FWCore.ParameterSet.Config as cms + +# This modifier is for injecting KF-based iterations in TICL. + +kf = cms.Modifier() diff --git a/DataFormats/HGCTrackingRecHit/BuildFile.xml b/DataFormats/HGCTrackingRecHit/BuildFile.xml new file mode 100644 index 0000000000000..e69de29bb2d1d diff --git a/DataFormats/HGCTrackingRecHit/interface/HGCTrackingRecHit.h b/DataFormats/HGCTrackingRecHit/interface/HGCTrackingRecHit.h new file mode 100644 index 0000000000000..fa2cfbdd4058b --- /dev/null +++ b/DataFormats/HGCTrackingRecHit/interface/HGCTrackingRecHit.h @@ -0,0 +1,27 @@ +#ifndef RecoHGCal_HGCTracking_HGCTrackingRecHit_h +#define RecoHGCal_HGCTracking_HGCTrackingRecHit_h + +/// Wrapper for TrackingRecHits for PatternRecognitionByKF + +#include +#include "DataFormats/TrackingRecHit/interface/RecHit2DLocalPos.h" + +class HGCTrackingRecHit : public RecHit2DLocalPos { + public: + HGCTrackingRecHit() {} + HGCTrackingRecHit(DetId id, const LocalPoint &pos, const LocalError &err, const float energy): + RecHit2DLocalPos(id), + pos_(pos), err_(err), energy_(energy) {} + + virtual HGCTrackingRecHit * clone() const override { return new HGCTrackingRecHit(*this); } + virtual LocalPoint localPosition() const override { return pos_; } + virtual LocalError localPositionError() const override { return err_; } + float energy() const { return energy_; } + + protected: + LocalPoint pos_; + LocalError err_; + float energy_; +}; + +#endif \ No newline at end of file diff --git a/DataFormats/HGCTrackingRecHit/src/classes.h b/DataFormats/HGCTrackingRecHit/src/classes.h new file mode 100644 index 0000000000000..e69de29bb2d1d diff --git a/DataFormats/HGCTrackingRecHit/src/classes_def.xml b/DataFormats/HGCTrackingRecHit/src/classes_def.xml new file mode 100644 index 0000000000000..e69de29bb2d1d diff --git a/DataFormats/HGCalReco/interface/KFHit.h b/DataFormats/HGCalReco/interface/KFHit.h new file mode 100644 index 0000000000000..a86c1e72b212c --- /dev/null +++ b/DataFormats/HGCalReco/interface/KFHit.h @@ -0,0 +1,43 @@ +// Author: Mark Matthewman - mark.matthewman@cern.ch +// Date: 04/2023 + +// Note!!!! +// KFHit is a temporary class provided for easier handling of the PatternRecognitionByKF. +// Once the result format for the Algo has been decided, then the KFHit is obsolete and needs to be deleated. + +#ifndef DataFormats_HGCalReco_KFHit_h +#define DataFormats_HGCalReco_KFHit_h + +#include +#include +#include "DataFormats/Provenance/interface/ProductID.h" +#include "DataFormats/Math/interface/Vector3D.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" +#include + +struct KFHit +{ + explicit KFHit(TrajectoryStateOnSurface tsos, uint32_t id): + center(tsos.globalPosition()), + localError(tsos.localError().matrix()), + cartesianError(tsos.cartesianError().matrix()), + curvilinearError(tsos.curvilinearError().matrix()), + xx(tsos.localError().positionError().xx()), + xy(tsos.localError().positionError().xy()), + yy(tsos.localError().positionError().yy()), + charge(tsos.charge()), + detid(id) + {} + KFHit(){} + + GlobalPoint center; + AlgebraicSymMatrix55 localError; + AlgebraicSymMatrix66 cartesianError; + AlgebraicSymMatrix55 curvilinearError; + float xx; + float xy; + float yy; + int charge; + uint32_t detid; +}; +#endif diff --git a/DataFormats/HGCalReco/src/classes.h b/DataFormats/HGCalReco/src/classes.h index d871bfb485a71..d11063f9be659 100644 --- a/DataFormats/HGCalReco/src/classes.h +++ b/DataFormats/HGCalReco/src/classes.h @@ -1,6 +1,7 @@ #include #include #include "DataFormats/HGCalReco/interface/Trackster.h" +#include "DataFormats/HGCalReco/interface/KFHit.h" #include "DataFormats/HGCalReco/interface/TICLLayerTile.h" #include "DataFormats/HGCalReco/interface/TICLSeedingRegion.h" #include "DataFormats/HGCalReco/interface/TICLCandidate.h" diff --git a/DataFormats/HGCalReco/src/classes_def.xml b/DataFormats/HGCalReco/src/classes_def.xml index 9ac5b29b4f2e2..e20d75553d37e 100644 --- a/DataFormats/HGCalReco/src/classes_def.xml +++ b/DataFormats/HGCalReco/src/classes_def.xml @@ -54,4 +54,10 @@ + + + + + + diff --git a/Geometry/CommonTopologies/interface/HGCDiskGeomDet.h b/Geometry/CommonTopologies/interface/HGCDiskGeomDet.h new file mode 100644 index 0000000000000..2a7930fc9917a --- /dev/null +++ b/Geometry/CommonTopologies/interface/HGCDiskGeomDet.h @@ -0,0 +1,35 @@ +#ifndef RecoHGCal_HGCTracking_HGCDiskGeomDet_h +#define RecoHGCal_HGCTracking_HGCDiskGeomDet_h + +// (mmatthew) Wrapper class for GeomDet corresponding to one layer of HGCal. Used by PatternRecognitionByKF. + +#include "Geometry/CommonDetUnit/interface/GeomDet.h" +#include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h" +#include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h" +#include "FWCore/Utilities/interface/Exception.h" + +class HGCDiskGeomDet : public GeomDet { + public: + HGCDiskGeomDet(int subdet, int zside, int layer, float z, float rmin, float rmax, float radlen, float xi) ; + + int subdet() const { return subdet_; } + int zside() const { return zside_; } + int layer() const { return layer_; } + float rmin() const { return rmin_; } + float rmax() const { return rmax_; } + bool isSilicon() const { + if(subdet_ == 8 || subdet_ == 9){ + return true; + } + else if (subdet_ == 10){ + return false; + } + else throw cms::Exception("LogicError", "Subdetector not defined"); + } + + protected: + const int subdet_, zside_, layer_; + const float rmin_, rmax_; +}; + +#endif \ No newline at end of file diff --git a/Geometry/CommonTopologies/src/HGCGeomDet.cc b/Geometry/CommonTopologies/src/HGCGeomDet.cc new file mode 100644 index 0000000000000..aeb611446fdd3 --- /dev/null +++ b/Geometry/CommonTopologies/src/HGCGeomDet.cc @@ -0,0 +1,13 @@ +#include "Geometry/CommonTopologies/interface/HGCDiskGeomDet.h" +#include "DataFormats/GeometrySurface/interface/MediumProperties.h" +#include "DataFormats/GeometrySurface/interface/BoundDisk.h" + +HGCDiskGeomDet::HGCDiskGeomDet(int subdet, int zside, int layer, float z, float rmin, float rmax, float radlen, float xi) : + GeomDet( Disk::build(Disk::PositionType(0,0,z), Disk::RotationType(), SimpleDiskBounds(rmin, rmax, -20, 20)).get() ), + subdet_(subdet), zside_(zside), layer_(layer), rmin_(rmin), rmax_(rmax) +{ + if(subdet == 8 || subdet == 9) setDetId(DetId::Detector(subdet)); + if (radlen > 0) { + (const_cast(surface())).setMediumProperties(MediumProperties(radlen,xi)); + } +} diff --git a/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py b/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py index 9c5ea39db66ed..22e97a75cceb8 100644 --- a/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py +++ b/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py @@ -18,8 +18,9 @@ 'keep *_ticlTrackstersHFNoseTrk_*_*', 'keep *_ticlTrackstersHFNoseMIP_*_*', 'keep *_ticlTrackstersHFNoseHAD_*_*', - 'keep *_ticlTrackstersHFNoseMerge_*_*',] + - ['keep *_pfTICL_*_*'] + 'keep *_ticlTrackstersHFNoseMerge_*_*'] + + ['keep *_pfTICL_*_*'] + + ['keep *_ticlRecHitTile_*_*'] ) ) TICL_RECO.outputCommands.extend(TICL_AOD.outputCommands) @@ -33,7 +34,10 @@ ) TICL_FEVT.outputCommands.extend(TICL_RECO.outputCommands) +print("RecoHGCalEventContents!!!!!!!!!!!!") + def customiseHGCalOnlyEventContent(process): + print("customiseHGCalOnlyEventContent!!!!!!!!!!!!") def cleanOutputAndSet(outputModule, ticl_outputCommads): outputModule.outputCommands = ['drop *_*_*_*'] outputModule.outputCommands.extend(ticl_outputCommads) diff --git a/RecoHGCal/TICL/BuildFile.xml b/RecoHGCal/TICL/BuildFile.xml index 29a694d95b8eb..f38188dc45aa4 100644 --- a/RecoHGCal/TICL/BuildFile.xml +++ b/RecoHGCal/TICL/BuildFile.xml @@ -1,8 +1,10 @@ + + diff --git a/RecoHGCal/TICL/interface/PatternRecognitionAlgoBase.h b/RecoHGCal/TICL/interface/PatternRecognitionAlgoBase.h index 507777de13d29..82934521a6a49 100644 --- a/RecoHGCal/TICL/interface/PatternRecognitionAlgoBase.h +++ b/RecoHGCal/TICL/interface/PatternRecognitionAlgoBase.h @@ -8,6 +8,7 @@ #include #include "DataFormats/CaloRecHit/interface/CaloCluster.h" #include "DataFormats/HGCalReco/interface/Trackster.h" +#include "DataFormats/HGCalReco/interface/KFHit.h" #include "DataFormats/HGCalReco/interface/TICLLayerTile.h" #include "DataFormats/HGCalReco/interface/TICLSeedingRegion.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" @@ -16,6 +17,7 @@ #include "RecoHGCal/TICL/interface/commons.h" #include "FWCore/Framework/interface/ConsumesCollector.h" #include "PhysicsTools/TensorFlow/interface/TensorFlow.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" namespace edm { class Event; @@ -50,10 +52,17 @@ namespace ticl { const tensorflow::Session* tS) : ev(eV), es(eS), layerClusters(lC), mask(mS), layerClustersTime(lT), tiles(tL), regions(rG), tfSession(tS) {} }; + + // (mmatthew): makeTrajectories used only by PatternRecognitionByKF. + // TODO: combine with makeTracksters as pure virtual function + virtual void makeTrajectories(const Inputs& input, + std::vector& kfhits, + std::vector& prophits, + float& abs_fail){}; virtual void makeTracksters(const Inputs& input, std::vector& result, - std::unordered_map>& seedToTracksterAssociation) = 0; + std::unordered_map>& seedToTracksterAssociation){}; protected: int algo_verbosity_; diff --git a/RecoHGCal/TICL/plugins/BuildFile.xml b/RecoHGCal/TICL/plugins/BuildFile.xml index 8c631381e21b4..d1a12881cbcef 100644 --- a/RecoHGCal/TICL/plugins/BuildFile.xml +++ b/RecoHGCal/TICL/plugins/BuildFile.xml @@ -8,6 +8,7 @@ + @@ -27,9 +28,11 @@ + + diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionPluginFactory.cc b/RecoHGCal/TICL/plugins/PatternRecognitionPluginFactory.cc index ac6296db5bf3b..db69634feac8c 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionPluginFactory.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionPluginFactory.cc @@ -1,6 +1,7 @@ #include "RecoHGCal/TICL/plugins/PatternRecognitionPluginFactory.h" #include "PatternRecognitionbyCA.h" #include "PatternRecognitionbyCLUE3D.h" +#include "PatternRecognitionbyKF.h" #include "PatternRecognitionbyFastJet.h" #include "FWCore/ParameterSet/interface/ValidatedPluginFactoryMacros.h" #include "FWCore/ParameterSet/interface/ValidatedPluginMacros.h" @@ -9,5 +10,6 @@ EDM_REGISTER_VALIDATED_PLUGINFACTORY(PatternRecognitionFactory, "PatternRecognit EDM_REGISTER_VALIDATED_PLUGINFACTORY(PatternRecognitionHFNoseFactory, "PatternRecognitionHFNoseFactory"); DEFINE_EDM_VALIDATED_PLUGIN(PatternRecognitionFactory, ticl::PatternRecognitionbyCA, "CA"); DEFINE_EDM_VALIDATED_PLUGIN(PatternRecognitionFactory, ticl::PatternRecognitionbyCLUE3D, "CLUE3D"); +DEFINE_EDM_VALIDATED_PLUGIN(PatternRecognitionFactory, ticl::PatternRecognitionbyKF, "KF"); DEFINE_EDM_VALIDATED_PLUGIN(PatternRecognitionFactory, ticl::PatternRecognitionbyFastJet, "FastJet"); -DEFINE_EDM_VALIDATED_PLUGIN(PatternRecognitionHFNoseFactory, ticl::PatternRecognitionbyCA, "CA"); +DEFINE_EDM_VALIDATED_PLUGIN(PatternRecognitionHFNoseFactory, ticl::PatternRecognitionbyCA, "CA"); \ No newline at end of file diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionPluginFactory.h b/RecoHGCal/TICL/plugins/PatternRecognitionPluginFactory.h index 2b6308a216085..5220ee2881507 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionPluginFactory.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionPluginFactory.h @@ -14,4 +14,6 @@ typedef edmplugin::PluginFactory PatternRecognitionHFNoseFactory; + + #endif diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc index 8470bbea04675..542409669d320 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc @@ -153,6 +153,7 @@ void PatternRecognitionbyCLUE3D::makeTracksters( std::vector &result, std::unordered_map> &seedToTracksterAssociation) { // Protect from events with no seeding regions + if (input.regions.empty()) return; diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyFastJet.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyFastJet.h index 2524a29921019..e7d8668f5db9c 100644 --- a/RecoHGCal/TICL/plugins/PatternRecognitionbyFastJet.h +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyFastJet.h @@ -50,4 +50,4 @@ namespace ticl { }; } // namespace ticl -#endif +#endif \ No newline at end of file diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyKF.cc b/RecoHGCal/TICL/plugins/PatternRecognitionbyKF.cc new file mode 100644 index 0000000000000..3c5945b792902 --- /dev/null +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyKF.cc @@ -0,0 +1,481 @@ +// Author: Mark Matthewman - mark.matthewman@cern.ch +// Date: 05/2022 +#include +#include +#include +#include + +#include "PatternRecognitionbyKF.h" + +#include "DataFormats/TrajectorySeed/interface/PropagationDirection.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h" +#include "DataFormats/HGCTrackingRecHit/interface/HGCTrackingRecHit.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Utilities/interface/Exception.h" +#include "FWCore/Framework/interface/EventSetup.h" + +#include "Geometry/Records/interface/CaloGeometryRecord.h" +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" +#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h" +#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h" + +#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" + +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h" +#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h" +#include "TrackingTools/KalmanUpdators/interface/KFUpdator.h" +#include "TrackingTools/PatternTools/interface/TrajectoryStateUpdator.h" +#include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h" +#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" + +using namespace ticl; + +template +PatternRecognitionbyKF::PatternRecognitionbyKF(const edm::ParameterSet &conf, edm::ConsumesCollector iC) + : PatternRecognitionAlgoBaseT(conf, iC), + caloGeomToken_(iC.esConsumes()), + radlen_(conf.getParameter>("radlen")), + xi_(conf.getParameter>("xi")), + propName_(conf.getParameter("propagator")), + propNameOppo_(conf.getParameter("propagatorOpposite")), + bfieldtoken_(iC.esConsumes()), + propagatortoken_(iC.esConsumes(edm::ESInputTag("",propName_))), + propagatorOppoToken_(iC.esConsumes(edm::ESInputTag("",propNameOppo_))), + estimatorToken_(iC.esConsumes(edm::ESInputTag("","Chi2"))), + updatorToken_(iC.esConsumes(edm::ESInputTag("","KFUpdator"))), + trackToken_(iC.consumes(conf.getParameter("tracks"))), + hgcalRecHitsEEToken_(iC.consumes(conf.getParameter("HGCEEInput"))), + hgcalRecHitsFHToken_(iC.consumes(conf.getParameter("HGCFHInput"))), + hgcalRecHitsBHToken_(iC.consumes(conf.getParameter("HGCBHInput"))), + geomCacheId_(0), + rescale_(conf.getParameter("rescale")){ +}; + +template +void PatternRecognitionbyKF::calculateLocalError(DetId id, const HGCalGeometry* hgcalgeom_){ + // Calculations based on normalized second moment of area + if(rhtools_.isSilicon(id)){ + double A = hgcalgeom_->getArea(id); + double a = sqrt(2*A/(3*sqrt(3))); // side length hexagon + double var = pow(a,4)*5*sqrt(3)/(16*A); + lerr[id] = LocalError(var, 0, var); + } + else{ + const HGCalDDDConstants* ddd = &(hgcalgeom_->topology().dddConstants()); + // Get outer and inner radius + const GlobalPoint &pos = rhtools_.getPosition(id); + double r = sqrt(pos.x()*pos.x() + pos.y()*pos.y()); + auto radiusLayer = ddd->getRadiusLayer(rhtools_.getLayer(id)); + int idx = static_cast(std::lower_bound(radiusLayer.begin(), radiusLayer.end(),r)-radiusLayer.begin()); + double rmax = radiusLayer[idx]; + double rmin = radiusLayer[idx-1]; + // Get angles + double phi = rhtools_.getPhi(id) + M_PI; // set to radians [0, 2pi] + double dphi = rhtools_.getScintDEtaDPhi(id).second; + double phimin = phi - 0.5*dphi; + double phimax = phi + 0.5*dphi; + + // FIXME!!! Using getArea() function results in lower efficiency due to incorrect position of KF hits. To be understood and getArea() to be implemented + // double A = hgcalgeom_->getArea(id); TOD + double A = (rmax*rmax - rmin*rmin)*M_PI*dphi/(2*M_PI); + // Calculate local error + double ex2 = 1/(8*A) * (pow(rmax,4) - pow(rmin,4)) * (-phimin - sin(phimin)*cos(phimin) + phimax + sin(phimax)*cos(phimax)); + double ex = 1/(3*A) * (pow(rmax,3) - pow(rmin,3)) * (sin(phimax) - sin(phimin)); + double varx = ex2 - ex*ex; + double ey2 = 1/(8*A) * (pow(rmax,4) - pow(rmin,4)) * (-phimin + sin(phimin)*cos(phimin) + phimax - sin(phimax)*cos(phimax)); + double ey = 1/(3*A) * (pow(rmax,3) - pow(rmin,3)) * (cos(phimin) - cos(phimax)); + double vary = ey2 - ey*ey; + double varxy = 1/(16*A)*(pow(rmax,4)-pow(rmin,4))*(cos(2*phimin)-cos(2*phimax)) - ex*ey; + lerr[id] = LocalError(varx, varxy, vary); + } +} + +template +const HGCDiskGeomDet * PatternRecognitionbyKF::switchDisk(const HGCDiskGeomDet * from, + const std::vector &vec, + bool isSilicon) const{ + auto it = std::find(vec.begin(), vec.end(),from); + if (it == vec.end()) throw cms::Exception("LogicError", "nextDisk called with invalid starting disk"); + isSilicon? --it: ++it; + return *(it); +} + +template +const HGCDiskGeomDet * PatternRecognitionbyKF::nextDisk(const HGCDiskGeomDet * from, + PropagationDirection direction, + const std::vector &vec, + bool isSilicon) const{ + auto it = std::find(vec.begin(), vec.end(),from); + if (it == vec.end()) throw cms::Exception("LogicError", "nextDisk called with invalid starting disk"); + int currentLayer = (*it)->layer(); + if (direction == alongMomentum){ + if ((*it == vec.back()) || (*it == vec.rbegin()[1])) return nullptr; // if disk last silicon OR last scintillator disk + while ((*it)->layer()isSilicon() == isSilicon) && ((*(it))->layer()==currentLayer+1)) return *(it); + } + } else{ + if (it == vec.begin()) return nullptr; + while ((*it)->layer()>currentLayer-2){ + if ((it) == vec.begin()) break; + --it; + if(((*(it))->isSilicon() == isSilicon) && ((*(it))->layer()==currentLayer-1)) return *(it); + } + } + return nullptr; // Returns nullptr if nextDisk reached end of detector +} + +template +template +std::vector +PatternRecognitionbyKF::advanceOneLayer(const Start &start, + const HGCDiskGeomDet * disk, + const std::vector &disks, + const TILES &tiles, PropagationDirection direction, + bool &isSilicon, + TempTrajectory traj){ + + std::vector ret; + + const Propagator &prop = (direction == alongMomentum ? *propagator_ : *propagatorOppo_); + TrajectoryStateOnSurface tsos = prop.propagate(start, disk->surface()); + if (!tsos.isValid()) return ret; + float r = sqrt(pow(tsos.globalPosition().x(),2)+pow(tsos.globalPosition().y(),2)); + if (((disk->rmin() > r) && (!isSilicon)) || (((r > disk->rmax()) && (isSilicon)))) { + disk = switchDisk(disk, disks, isSilicon); + isSilicon = !isSilicon; + tsos = prop.propagate(start, disk->surface()); + } + + // Collect hits with estimate + int depth = disk->layer()+1; + auto meas = measurements(tsos, *estimator_, tiles, depth); + std::sort(meas.begin(), meas.end(),TrajMeasLessEstim()); + + for (const TrajectoryMeasurement &tm : meas){ + TrajectoryStateOnSurface updated = updator_->update(tm.forwardPredictedState(),*tm.recHit()); + ret.push_back(traj.foundHits() ? traj: TempTrajectory(traj.direction(),0)); + ret.back().push(TrajectoryMeasurement(tm.forwardPredictedState(), + updated, + tm.recHit(), + tm.estimate()), + tm.estimate()); + } + + if (!meas.empty()) return ret; + auto missing = TrackingRecHit::missing; + ret.push_back(traj.foundHits()? traj : TempTrajectory(traj.direction(),0)); + ret.back().push(TrajectoryMeasurement(tsos, std::make_shared(*disk,missing))); + + return ret; +} + +template +void PatternRecognitionbyKF::makeDisks(int subdet, int disks, const CaloGeometry* geom_) { + + const CaloSubdetectorGeometry *subGeom = geom_->getSubdetectorGeometry(DetId::Detector(subdet), ForwardSubdetector::ForwardEmpty); + auto geomEE = static_cast(subGeom); + const HGCalDDDConstants* ddd = &(geomEE->topology().dddConstants()); + + std::vector rmax(disks, 0), rmin(disks, 9e9); + std::vector zsumPos(disks), zsumNeg(disks); + std::vector countPos(disks), countNeg(disks); + const std::vector & ids = subGeom->getValidDetIds(); + + for (auto & i : ids) { + calculateLocalError(i,geomEE); + + const GlobalPoint & pos = rhtools_.getPosition(i); + int layer = rhtools_.getLayer(i)-1; + float z = pos.z(); + float rho = pos.perp(); + int side = z > 0 ? +1 : -1; + + (side > 0 ? zsumPos : zsumNeg)[layer] += z; + (side > 0 ? countPos : countNeg)[layer]++; + if (rho > rmax[layer]) rmax[layer] = rho; + if (rho < rmin[layer]) rmin[layer] = rho; + } + + int layer = ddd->getLayerOffset(); // FIXME: Can be made less stupid + for (int i = 0; i < disks; ++i) { + if (countPos[i]) { + HGCDiskGeomDet* disk = new HGCDiskGeomDet(subdet, +1, layer, zsumPos[i]/countPos[i], rmin[i], rmax[i], radlen_[layer], xi_[layer]); + addDisk(disk); + } + if (countNeg[i]) { + HGCDiskGeomDet* disk = new HGCDiskGeomDet(subdet, -1, layer, zsumNeg[i]/countNeg[i], rmin[i], rmax[i], radlen_[layer], xi_[layer]); + addDisk(disk); + } + layer++; + } +} + + +template +void PatternRecognitionbyKF::fillHitMap(std::map& hitMap, + const HGCRecHitCollection& rechitsEE, + const HGCRecHitCollection& rechitsFH, + const HGCRecHitCollection& rechitsBH) const { + + hitMap.clear(); + for (const auto& hit : rechitsEE) { + hitMap.emplace(hit.detid(), &hit); + } + + for (const auto& hit : rechitsFH) { + hitMap.emplace(hit.detid(), &hit); + } + + for (const auto& hit : rechitsBH) { + hitMap.emplace(hit.detid(), &hit); + } +} + +template +std::vector PatternRecognitionbyKF::measurements( + const TrajectoryStateOnSurface &tsos, + const MeasurementEstimator &mest, + const TILES &tiles, + int depth){ + + std::vector ret; + + // define search window and get bins + + float eta = tsos.globalPosition().eta(); + float phi = tsos.globalPosition().phi(); + + float etaMin = eta - etaBinSize; + float etaMax = eta + etaBinSize; + float phiMin = phi - phiBinSize; + float phiMax = phi + phiBinSize; + + auto bins = tiles[depth].searchBoxEtaPhi(etaMin, etaMax, phiMin, phiMax); + + // loop over candidates + + for (int ieta = bins[0]; ieta < bins[1]; ieta++) { + auto offset = ieta * nPhiBin; + for (int phi = bins[2]; phi < bins[3]; phi++) { + int iphi = ((phi % nPhiBin + nPhiBin) % nPhiBin); + if (!tiles[depth][offset + iphi].empty()) { + for(auto hit: tiles[depth][offset + iphi]) { + const auto rec = hitMap.find(hit)->second; + float energy = rec->energy(); + + GlobalPoint globalpoint = rhtools_.getPosition(hit); + LocalPoint localpoint = tsos.surface().toLocal(globalpoint); + + auto hitptr = std::make_shared(hit,localpoint,lerr[hit],energy); + auto mest_pair = mest.estimate(tsos, *hitptr); + if(mest_pair.first){ + ret.emplace_back(tsos,hitptr,mest_pair.second); + } + } + } + } + } +return ret; +} + +template +void PatternRecognitionbyKF::init( + const edm::Event& ev, const edm::EventSetup& es){ + + //Get Calo Geometry + if (es.get().cacheIdentifier() != geomCacheId_) { + geomCacheId_ = es.get().cacheIdentifier(); + const CaloGeometry* geom = &es.getData(caloGeomToken_); + rhtools_.setGeometry(*geom); + + makeDisks(8, 26, geom); + makeDisks(9, 21, geom); + makeDisks(10,21, geom); + + auto ptrSort = [](const HGCDiskGeomDet *a, const HGCDiskGeomDet *b) -> bool { return (abs(a->position().z())) < (abs(b->position().z())); }; + std::sort(disksPos_.begin(), disksPos_.end(), ptrSort); + std::sort(disksNeg_.begin(), disksNeg_.end(), ptrSort); + } + + bfield_ = es.getHandle(bfieldtoken_); + propagator_ = es.getHandle(propagatortoken_); // Can I move this up? + propagatorOppo_ = es.getHandle(propagatorOppoToken_); // Can I move this up? + estimator_ = es.getHandle(estimatorToken_); // Can I move this up? + updator_ = es.getHandle(updatorToken_); // Can I move this up? + + edm::Handle ee_hits; + edm::Handle fh_hits; + edm::Handle bh_hits; + + ev.getByToken(hgcalRecHitsEEToken_, ee_hits); + ev.getByToken(hgcalRecHitsFHToken_, fh_hits); + ev.getByToken(hgcalRecHitsBHToken_, bh_hits); + fillHitMap(hitMap, *ee_hits, *fh_hits, *bh_hits); +} + + +template +void PatternRecognitionbyKF::makeTrajectories( + const typename PatternRecognitionAlgoBaseT::Inputs &input, + std::vector& kfhits, + std::vector& prophits, + float& abs_fail) { + + edm::EventSetup const &es = input.es; + edm::Event const &ev = input.ev; + init(ev,es); + const TILES &tiles = input.tiles; + + // Build tsos from TrackCollection + edm::Handle tracks_h; + ev.getByToken(trackToken_,tracks_h); + + const reco::TrackCollection& tkx = *tracks_h; + if (tkx.empty()){ + std::cout << "No track found!!!! Exited PatternRecognitionbyKF" << std::endl; + abs_fail+=1; + return; + } + + FreeTrajectoryState fts = trajectoryStateTransform::outerFreeState(tkx.front(),bfield_.product()); + std::cout << "Rescale Factor: None" << std::endl; + // Propagate through all disks + + // Get first disk + + int zside = fts.momentum().eta() > 0 ? +1 : -1; + PropagationDirection direction = alongMomentum; + std::vector disks = (zside > 0? disksPos_ : disksNeg_); + const HGCDiskGeomDet* disk = (zside > 0 ? disksPos_ : disksNeg_).front(); + bool isSilicon = true; + + // Propagation step + + std::vector traj = advanceOneLayer(fts, disk, disks, tiles, direction, isSilicon, TempTrajectory(direction,0)); + if (traj.empty()){ + std::cout << "No track found!!!! Exited PatternRecognitionbyKF" << std::endl; + abs_fail+=1; + return; + } + auto lm = traj.back().lastMeasurement(); + TrajectoryStateOnSurface tsos_prop = lm.predictedState(); + KFHit *prophit = new KFHit(tsos_prop, lm.recHit()->geographicalId()); + prophits.push_back(*prophit); + + TrajectoryStateOnSurface tsos_kf = lm.updatedState(); + KFHit *kfhit = new KFHit(tsos_kf, lm.recHit()->geographicalId()); + kfhits.push_back(*kfhit); + + std::vector traj_prop; + std::vector traj_kf; + traj_prop.push_back(traj.back()); + traj_kf.push_back(traj.back()); + + // Loop over all disks + + unsigned int depth = 2; + for(disk = nextDisk(disk, direction, disks, isSilicon); disk != nullptr; disk = nextDisk(disk, direction, disks, isSilicon), depth++){ + std::vector newcands_kf; + for(TempTrajectory & cand : traj_kf){ + + TrajectoryStateOnSurface start = cand.lastMeasurement().updatedState(); + std::vector hisTrajs = advanceOneLayer(start, disk, disks, tiles, direction, isSilicon, cand); + for(TempTrajectory & t : hisTrajs){ + auto lm = t.lastMeasurement(); + newcands_kf.push_back(t); + + TrajectoryStateOnSurface tsos_kf = lm.updatedState(); + KFHit *kfhit = new KFHit(tsos_kf, lm.recHit()->geographicalId()); + kfhits.push_back(*kfhit); + break; + } + } + traj_kf.swap(newcands_kf); + } + + // Propagator loop used purely for testing purposes + direction = alongMomentum; + disk = (zside > 0 ? disksPos_ : disksNeg_).front(); + isSilicon=1; + depth = 2; + for(disk = nextDisk(disk, direction, disks, isSilicon); disk != nullptr; disk = nextDisk(disk, direction, disks, isSilicon), depth++){ + std::vector newcands_prop; + for(TempTrajectory & cand : traj_prop){ + + TrajectoryStateOnSurface start = cand.lastMeasurement().predictedState(); + std::vector hisTrajs = advanceOneLayer(start, disk, disks, tiles, direction, isSilicon, cand); + + for(TempTrajectory & t : hisTrajs){ + auto lm = t.lastMeasurement(); + newcands_prop.push_back(t); + TrajectoryStateOnSurface tsos_prop = lm.predictedState(); + KFHit *prophit = new KFHit(tsos_prop, lm.recHit()->geographicalId()); + + std::cout << "-------------- Error -----------" << std::endl; + std::cout << "Local Error," < +void PatternRecognitionbyKF::dumpTiles(const TILES &tiles) const { + std::cout << "Entered dumpTiles" << std::endl; + constexpr int nEtaBin = TILES::constants_type_t::nEtaBins; + constexpr int nPhiBin = TILES::constants_type_t::nPhiBins; + std::cout << nEtaBin << "\t" <(rhtools_.lastLayer(false)); + std::cout << lastLayerPerSide << std::endl; + int maxLayer = 2 * lastLayerPerSide - 1; + int count = 0; + for (int layer = 0; layer <= maxLayer; layer++) { + for (int ieta = 0; ieta < nEtaBin; ieta++) { + auto offset = ieta * nPhiBin; + for (int phi = 0; phi < nPhiBin; phi++) { + int iphi = ((phi % nPhiBin + nPhiBin) % nPhiBin); + if (!tiles[layer][offset + iphi].empty()) { + std::cout << "Layer: " << layer << " ieta: " << ieta << " phi: " << phi + << " " << tiles[layer][offset + iphi].size() << std::endl; + count++; + } + } + } + } + std::cout << "Number of RecHits: " << count << std::endl; +} + +template +void PatternRecognitionbyKF::fillPSetDescription(edm::ParameterSetDescription &iDesc) { + iDesc.add("algo_verbosity", 0); + //iDesc.add("propagator", "RungeKuttaTrackerPropagator"); // RungeKutta Propagator + iDesc.add>("radlen",{}); + iDesc.add>("xi",{}); + iDesc.add("propagator", "PropagatorWithMaterial"); // Analytical Propagator + iDesc.add("propagatorOpposite", "PropagatorWithMaterialOpposite"); + iDesc.add("estimator", "Chi2"); + iDesc.add("tracks", edm::InputTag("generalTracks")); + iDesc.add("HGCEEInput", edm::InputTag("HGCalRecHit", "HGCEERecHits")); + iDesc.add("HGCFHInput", edm::InputTag("HGCalRecHit", "HGCHEFRecHits")); + iDesc.add("HGCBHInput", edm::InputTag("HGCalRecHit", "HGCHEBRecHits")); + iDesc.add("rescale",1); +} + +template class ticl::PatternRecognitionbyKF; +template class ticl::PatternRecognitionbyKF; diff --git a/RecoHGCal/TICL/plugins/PatternRecognitionbyKF.h b/RecoHGCal/TICL/plugins/PatternRecognitionbyKF.h new file mode 100644 index 0000000000000..0ef73f53f4174 --- /dev/null +++ b/RecoHGCal/TICL/plugins/PatternRecognitionbyKF.h @@ -0,0 +1,114 @@ +// Author: Marco Rovere - marco.rovere@cern.ch +// Date: 04/2021 + +#ifndef __RecoHGCal_TICL_PRbyKF_H__ +#define __RecoHGCal_TICL_PRbyKF_H__ +#include // unique_ptr + +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" +#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h" + +#include "Geometry/CommonTopologies/interface/HGCDiskGeomDet.h" +#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h" + +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "MagneticField/Engine/interface/MagneticField.h" + +#include "RecoHGCal/TICL/interface/PatternRecognitionAlgoBase.h" + +#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" + +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" +#include "TrackingTools/PatternTools/interface/TrajectoryStateUpdator.h" +#include "TrackingTools/PatternTools/interface/TempTrajectory.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" +#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h" +#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h" + + +namespace ticl { + template + class PatternRecognitionbyKF final : public PatternRecognitionAlgoBaseT { + public: + PatternRecognitionbyKF(const edm::ParameterSet& conf, edm::ConsumesCollector); + ~PatternRecognitionbyKF() override = default; + + void makeTrajectories(const typename PatternRecognitionAlgoBaseT::Inputs& input, + std::vector& kfhits, + std::vector& prophits, + float& abs_fail) override; + + static void fillPSetDescription(edm::ParameterSetDescription& iDesc); + + private: + // Declarations for Constructor + edm::ESGetToken caloGeomToken_; + const std::vector radlen_; + const std::vector xi_; + const std::string propName_; + const std::string propNameOppo_; + edm::ESGetToken bfieldtoken_; + edm::ESGetToken propagatortoken_; + edm::ESGetToken propagatorOppoToken_; + edm::ESGetToken estimatorToken_; + edm::ESGetToken updatorToken_; + edm::EDGetTokenT trackToken_; + edm::EDGetTokenT hgcalRecHitsEEToken_; + edm::EDGetTokenT hgcalRecHitsFHToken_; + edm::EDGetTokenT hgcalRecHitsBHToken_; + edm::ESHandle bfield_; + edm::ESHandle propagator_; + edm::ESHandle estimator_; + edm::ESHandle updator_; + edm::ESHandle propagatorOppo_; + uint32_t geomCacheId_; + int rescale_; + + // Instance Variables + hgcal::RecHitTools rhtools_; + std::map hitMap; + std::map lerr; + std::vector disksPos_, disksNeg_; + + float etaBinSize = (TILES::constants_type_t::maxEta - TILES::constants_type_t::minEta)/TILES::constants_type_t::nEtaBins; + float phiBinSize = 2*M_PI/TILES::constants_type_t::nPhiBins; + int nPhiBin = TILES::constants_type_t::nPhiBins; + int nEtaBin = TILES::constants_type_t::nEtaBins; + + //Member Functions + void dumpTiles(const TILES&) const; + void makeDisks(int subdet, int disks, const CaloGeometry* geom_); + void addDisk(HGCDiskGeomDet *disk) { + (disk->zside() > 0 ? disksPos_ : disksNeg_).push_back(disk); + } + std::vector measurements(const TrajectoryStateOnSurface &tsos, + const MeasurementEstimator &mest, + const TILES &tiles, + int depth); + template + std::vector advanceOneLayer(const Start &start, + const HGCDiskGeomDet * disk, + const std::vector &disks, + const TILES &tiles, + PropagationDirection direction, + bool &isSilicon, + TempTrajectory traj); + const HGCDiskGeomDet * nextDisk(const HGCDiskGeomDet * from, + PropagationDirection direction, + const std::vector &vec, + bool isSilicon) const; + const HGCDiskGeomDet * switchDisk(const HGCDiskGeomDet * from, + const std::vector &vec, + bool isSilicon) const; + virtual void fillHitMap(std::map& hitMap, + const HGCRecHitCollection& recHitsEE, + const HGCRecHitCollection& recHitsFH, + const HGCRecHitCollection& recHitsBH) const; + void calculateLocalError(DetId id, + const HGCalGeometry* hgcalgeom); + void init(const edm::Event& evt, const edm::EventSetup& es); + }; +} // namespace ticl +#endif diff --git a/RecoHGCal/TICL/plugins/TICLLayerTileProducer.cc b/RecoHGCal/TICL/plugins/TICLLayerTileProducer.cc index e3b1c633eddbf..751514e722a2e 100644 --- a/RecoHGCal/TICL/plugins/TICLLayerTileProducer.cc +++ b/RecoHGCal/TICL/plugins/TICLLayerTileProducer.cc @@ -12,8 +12,17 @@ #include "DataFormats/CaloRecHit/interface/CaloCluster.h" #include "DataFormats/HGCalReco/interface/TICLLayerTile.h" +#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h" #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" +#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h" +#include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h" + +#include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h" +#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h" + + + class TICLLayerTileProducer : public edm::stream::EDProducer<> { public: explicit TICLLayerTileProducer(const edm::ParameterSet &ps); @@ -23,27 +32,114 @@ class TICLLayerTileProducer : public edm::stream::EDProducer<> { static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); private: + virtual void fillHitMap(std::map& hitMap, + const HGCRecHitCollection& recHitsEE, + const HGCRecHitCollection& recHitsFH, + const HGCRecHitCollection& recHitsBH) const; + virtual void fillRecHits(std::vector& recHits, + const HGCRecHitCollection& recHitsEE, + const HGCRecHitCollection& recHitsFH, + const HGCRecHitCollection& recHitsBH) const; + //template + //void fillTiles(T& objects, U& results); edm::EDGetTokenT> clusters_token_; edm::EDGetTokenT> clusters_HFNose_token_; + edm::EDGetTokenT hgcalRecHitsEEToken_; + edm::EDGetTokenT hgcalRecHitsFHToken_; + edm::EDGetTokenT hgcalRecHitsBHToken_; edm::ESGetToken geometry_token_; + //edm::EDGetTokenT> caloParticlesToken_; hgcal::RecHitTools rhtools_; std::string detector_; + bool isLC_; bool doNose_; }; +void TICLLayerTileProducer::fillRecHits(std::vector& recHits, + const HGCRecHitCollection& rechitsEE, + const HGCRecHitCollection& rechitsFH, + const HGCRecHitCollection& rechitsBH) const { + recHits.clear(); + for (const auto& hit : rechitsEE) { + recHits.push_back(&hit); + } + + for (const auto& hit : rechitsFH) { + recHits.push_back(&hit); + } + + for (const auto& hit : rechitsBH) { + recHits.push_back(&hit); + } +} // end of EfficiencyStudies::fillHitMap + +void TICLLayerTileProducer::fillHitMap(std::map& hitMap, + const HGCRecHitCollection& rechitsEE, + const HGCRecHitCollection& rechitsFH, + const HGCRecHitCollection& rechitsBH) const { + hitMap.clear(); + for (const auto& hit : rechitsEE) { + hitMap.emplace(hit.detid(), &hit); + } + + for (const auto& hit : rechitsFH) { + hitMap.emplace(hit.detid(), &hit); + } + + for (const auto& hit : rechitsBH) { + hitMap.emplace(hit.detid(), &hit); + } +} // end of EfficiencyStudies::fillHitMap + +/* +template +void TICLLayerTileProducer::fillTiles(T& objects, U& result){ + int objId = 0; + int layer; + float eta, phi; + for (auto const &obj : objects) { + if (typeid(obj)){ + const auto firstHitDetId = obj.hitsAndFractions()[0].first; + layer = rhtools_.getLayerWithOffset(firstHitDetId) + + rhtools_.lastLayer(doNose_) * ((rhtools_.zside(firstHitDetId) + 1) >> 1) - 1; + eta = obj.eta(); + phi = obj.phi(); + } + else { + layer = rhtools_.getLayerWithOffset(obj->detid()); + eta = rhtools_.getEta(obj->detid()); + phi = rhtools_.getPhi(obj->detid()); + } + + + assert(layer >= 0); + result->fill(layer, eta, phi, objId); + LogDebug("TICLLayerTileProducer") << "Adding layerClusterId: " << objId << " into bin [eta,phi]: [ " + << (*result)[layer].etaBin(eta) << ", " << (*result)[layer].phiBin(phi) + << "] for layer: " << layer << std::endl; + objId++; + } +} +*/ TICLLayerTileProducer::TICLLayerTileProducer(const edm::ParameterSet &ps) - : detector_(ps.getParameter("detector")) { + : detector_(ps.getParameter("detector")), + isLC_(ps.getParameter("isLC")) { clusters_HFNose_token_ = consumes>(ps.getParameter("layer_HFNose_clusters")); clusters_token_ = consumes>(ps.getParameter("layer_clusters")); + hgcalRecHitsEEToken_ = consumes(ps.getParameter("HGCEEInput")); + hgcalRecHitsFHToken_ = consumes(ps.getParameter("HGCFHInput")); + hgcalRecHitsBHToken_ = consumes(ps.getParameter("HGCBHInput")); + //caloParticlesToken_ = consumes>(ps.getParameter("caloParticles")); geometry_token_ = esConsumes(); doNose_ = (detector_ == "HFNose"); if (doNose_) - produces(); + produces("Test"); else - produces(); + produces("Test"); + produces>("Test"); } void TICLLayerTileProducer::beginRun(edm::Run const &, edm::EventSetup const &es) { @@ -52,37 +148,158 @@ void TICLLayerTileProducer::beginRun(edm::Run const &, edm::EventSetup const &es } void TICLLayerTileProducer::produce(edm::Event &evt, const edm::EventSetup &) { + + std::cout << "in TICLLayerTileProducer.cc, start putting RecHits in event" << std::endl; + auto result = std::make_unique(); auto resultHFNose = std::make_unique(); - + auto test = std::make_unique>(); + edm::Handle> cluster_h; + edm::Handle ee_hits; + edm::Handle fh_hits; + edm::Handle bh_hits; + std::map hitMap; + std::vector recHits; + + /* + edm::Handle> CaloParticles; + evt.getByToken(caloParticlesToken_, CaloParticles); + const CaloParticleCollection& cps = *CaloParticles; + + */ + /* if (doNose_) evt.getByToken(clusters_HFNose_token_, cluster_h); else evt.getByToken(clusters_token_, cluster_h); + */ + + if (isLC_){ + doNose_ ? evt.getByToken(clusters_HFNose_token_, cluster_h) : evt.getByToken(clusters_token_, cluster_h); + //const auto &objects = *cluster_h; + //doNose_ ? fillTiles(*cluster_h, result) : fillTiles(*cluster_h, resultHFNose); + + int lcId = 0; + for (auto const &lc : *cluster_h) { + const auto firstHitDetId = lc.hitsAndFractions()[0].first; + int layer = rhtools_.getLayerWithOffset(firstHitDetId) + + rhtools_.lastLayer(doNose_) * ((rhtools_.zside(firstHitDetId) + 1) >> 1) - 1; + assert(layer >= 0); + + doNose_ ? resultHFNose->fill(layer, lc.eta(), lc.phi(), lcId) : result->fill(layer, lc.eta(), lc.phi(), lcId); + + LogDebug("TICLLayerTileProducer") << "Adding layerClusterId: " << lcId << " into bin [eta,phi]: [ " + << (*result)[layer].etaBin(lc.eta()) << ", " << (*result)[layer].phiBin(lc.phi()) + << "] for layer: " << layer << std::endl; + lcId++; + } + } + else { + evt.getByToken(hgcalRecHitsEEToken_, ee_hits); + evt.getByToken(hgcalRecHitsFHToken_, fh_hits); + evt.getByToken(hgcalRecHitsBHToken_, bh_hits); + fillHitMap(hitMap, *ee_hits, *fh_hits, *bh_hits); + fillRecHits(recHits, *ee_hits, *fh_hits, *bh_hits); + //doNose_ ? fillTiles(recHits, result) : fillTiles(recHits, resultHFNose); + + //doNose_ ? evt.getByToken(clusters_HFNose_token_, cluster_h) : evt.getByToken(clusters_token_, cluster_h); + //const auto &objects = *cluster_h; + //doNose_ ? fillTiles(*cluster_h, result) : fillTiles(*cluster_h, resultHFNose); + + int objId = 0; + for (auto const &rhit : recHits) { + int layer = rhtools_.getLayerWithOffset(rhit->detid()); + float eta = rhtools_.getEta(rhit->detid()); + float phi = rhtools_.getPhi(rhit->detid()); - const auto &layerClusters = *cluster_h; - int lcId = 0; - for (auto const &lc : layerClusters) { - const auto firstHitDetId = lc.hitsAndFractions()[0].first; - int layer = rhtools_.getLayerWithOffset(firstHitDetId) + - rhtools_.lastLayer(doNose_) * ((rhtools_.zside(firstHitDetId) + 1) >> 1) - 1; + assert(layer >= 0); + + objId = rhit->detid().rawId(); + + if (doNose_){ + resultHFNose->fill(layer, eta, phi, objId); + } + else { + result->fill(layer, eta, phi, objId); + } + LogDebug("TICLLayerTileProducer") << "Adding layerClusterId: " << objId << " into bin [eta,phi]: [ " + << (*result)[layer].etaBin(eta) << ", " << (*result)[layer].phiBin(phi) + << "] for layer: " << layer << std::endl; + objId++; + } + /* + int objId=0; + for (const auto& it_cp : cps ) { + const CaloParticle& cp = ((it_cp)); + const SimClusterRefVector& simclusters = cp.simClusters(); + for (const auto& it_simc : simclusters) { + const SimCluster& simc = (*(it_simc)); + const auto& sc_haf = simc.hits_and_fractions(); + for (const auto& it_sc_haf : sc_haf) { + DetId detid_ = (it_sc_haf.first); + std::map::const_iterator itcheck = hitMap.find(detid_); + if (itcheck != hitMap.end()) { + int layer = rhtools_.getLayerWithOffset(detid_); + float eta = rhtools_.getEta(detid_); + float phi = rhtools_.getPhi(detid_); + + assert(layer >= 0); + + if (doNose_){ + resultHFNose->fill(layer, eta, phi, objId); + } + else { + result->fill(layer, eta, phi, objId); + } + LogDebug("TICLLayerTileProducer") << "Adding recHitId: " << objId << " into bin [eta,phi]: [ " + << (*result)[layer].etaBin(eta) << ", " << (*result)[layer].phiBin(phi) + << "] for layer: " << layer << std::endl; + objId++; + } + } + } + } + */ + } + + //const auto &object = (isLC_? recHits : *cluster_h) + /* + + int objId = 0; + for (auto const &obj : objects) { + if (isLC_){ + const auto firstHitDetId = lc.hitsAndFractions()[0].first; + int layer = rhtools_.getLayerWithOffset(firstHitDetId) + + rhtools_.lastLayer(doNose_) * ((rhtools_.zside(firstHitDetId) + 1) >> 1) - 1; + } + else { + int layer = obj.layer(); + } assert(layer >= 0); - if (doNose_) - resultHFNose->fill(layer, lc.eta(), lc.phi(), lcId); - else - result->fill(layer, lc.eta(), lc.phi(), lcId); - LogDebug("TICLLayerTileProducer") << "Adding layerClusterId: " << lcId << " into bin [eta,phi]: [ " - << (*result)[layer].etaBin(lc.eta()) << ", " << (*result)[layer].phiBin(lc.phi()) + if (doNose_){ + resultHFNose->fill(layer, obj.eta(), obj.phi(), objId); + } + else { + result->fill(layer, obj.eta(), obj.phi(), objId); + } + LogDebug("TICLLayerTileProducer") << "Adding layerClusterId: " << objId << " into bin [eta,phi]: [ " + << (*result)[layer].etaBin(obj.eta()) << ", " << (*result)[layer].phiBin(obj.phi()) << "] for layer: " << layer << std::endl; - lcId++; + objId++; } - if (doNose_) - evt.put(std::move(resultHFNose)); - else - evt.put(std::move(result)); + */ + + std::cout << "put" << std::endl; + if (doNose_){ + evt.put(std::move(resultHFNose),"Test"); + } + else{ + evt.put(std::move(result),"Test"); + } + evt.put(std::move(test),"Test"); } void TICLLayerTileProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { @@ -90,6 +307,11 @@ void TICLLayerTileProducer::fillDescriptions(edm::ConfigurationDescriptions &des desc.add("detector", "HGCAL"); desc.add("layer_clusters", edm::InputTag("hgcalLayerClusters")); desc.add("layer_HFNose_clusters", edm::InputTag("hgcalLayerClustersHFNose")); + desc.add("isLC", true); + desc.add("HGCEEInput", edm::InputTag("HGCalRecHit", "HGCEERecHits")); + desc.add("HGCFHInput", edm::InputTag("HGCalRecHit", "HGCHEFRecHits")); + desc.add("HGCBHInput", edm::InputTag("HGCalRecHit", "HGCHEBRecHits")); + //desc.add("caloParticles", edm::InputTag("mix", "MergedCaloTruth")); descriptions.add("ticlLayerTileProducer", desc); } diff --git a/RecoHGCal/TICL/plugins/TrackstersProducer.cc b/RecoHGCal/TICL/plugins/TrackstersProducer.cc index 005d19ba1ec2d..682d228ec6efd 100644 --- a/RecoHGCal/TICL/plugins/TrackstersProducer.cc +++ b/RecoHGCal/TICL/plugins/TrackstersProducer.cc @@ -18,19 +18,25 @@ #include "DataFormats/ParticleFlowReco/interface/PFCluster.h" #include "DataFormats/HGCalReco/interface/Trackster.h" +#include "DataFormats/HGCalReco/interface/KFHit.h" #include "DataFormats/HGCalReco/interface/TICLLayerTile.h" #include "DataFormats/HGCalReco/interface/TICLSeedingRegion.h" #include "RecoHGCal/TICL/plugins/PatternRecognitionPluginFactory.h" #include "PatternRecognitionbyCA.h" +#include "PatternRecognitionbyKF.h" #include "PatternRecognitionbyMultiClusters.h" #include "PhysicsTools/TensorFlow/interface/TfGraphRecord.h" #include "PhysicsTools/TensorFlow/interface/TensorFlow.h" #include "PhysicsTools/TensorFlow/interface/TfGraphDefWrapper.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" +#include "TrackingTools/PatternTools/interface/TempTrajectory.h" + using namespace ticl; + class TrackstersProducer : public edm::stream::EDProducer<> { public: explicit TrackstersProducer(const edm::ParameterSet&); @@ -57,6 +63,7 @@ class TrackstersProducer : public edm::stream::EDProducer<> { const edm::EDGetTokenT> original_layerclusters_mask_token_; const edm::EDGetTokenT>> clustersTime_token_; edm::EDGetTokenT layer_clusters_tiles_token_; + edm::EDGetTokenT rechit_tiles_token_; edm::EDGetTokenT layer_clusters_tiles_hfnose_token_; const edm::EDGetTokenT> seeding_regions_token_; const std::string itername_; @@ -104,7 +111,10 @@ TrackstersProducer::TrackstersProducer(const edm::ParameterSet& ps) produces>(); produces>(); // Mask to be applied at the next iteration -} + produces("Abs Fail").setBranchAlias("Abs Fail"); + produces>("KFHits").setBranchAlias("KFHits"); + produces>("PropHits").setBranchAlias("PropHits"); + } void TrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { // hgcalMultiClusters @@ -114,8 +124,8 @@ void TrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descri desc.add("filtered_mask", edm::InputTag("filteredLayerClusters", "iterationLabelGoesHere")); desc.add("original_mask", edm::InputTag("hgcalLayerClusters", "InitialLayerClustersMask")); desc.add("time_layerclusters", edm::InputTag("hgcalLayerClusters", "timeLayerCluster")); - desc.add("layer_clusters_tiles", edm::InputTag("ticlLayerTileProducer")); - desc.add("layer_clusters_hfnose_tiles", edm::InputTag("ticlLayerTileHFNose")); + desc.add("layer_clusters_tiles", edm::InputTag("ticlLayerTileProducer", "Test")); + desc.add("layer_clusters_hfnose_tiles", edm::InputTag("ticlLayerTileHFNose", "Test")); desc.add("seeding_regions", edm::InputTag("ticlSeedingRegionProducer")); desc.add("patternRecognitionBy", "CA"); desc.add("itername", "unknown"); @@ -131,6 +141,11 @@ void TrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descri pluginDescClue3D.addNode(edm::PluginDescription("type", "CLUE3D", true)); desc.add("pluginPatternRecognitionByCLUE3D", pluginDescClue3D); + // KF Plugin + edm::ParameterSetDescription pluginDescKF; + pluginDescKF.addNode(edm::PluginDescription("type", "KF", true)); + desc.add("pluginPatternRecognitionByKF", pluginDescKF); + // FastJet Plugin edm::ParameterSetDescription pluginDescFastJet; pluginDescFastJet.addNode(edm::PluginDescription("type", "FastJet", true)); @@ -142,6 +157,11 @@ void TrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descri void TrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) { auto result = std::make_unique>(); auto output_mask = std::make_unique>(); + auto abs_fail = std::make_unique(); + auto kfhits = std::make_unique>(); + auto prophits = std::make_unique>(); + + std::cout<& original_layerclusters_mask = evt.get(original_layerclusters_mask_token_); const auto& layerClusters = evt.get(clusters_token_); @@ -170,7 +190,6 @@ void TrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) { layer_clusters_hfnose_tiles, seeding_regions, tfSession_); - myAlgoHFNose_->makeTracksters(inputHFNose, *result, seedToTrackstersAssociation); } else { @@ -178,7 +197,14 @@ void TrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) { const typename PatternRecognitionAlgoBaseT::Inputs input( evt, es, layerClusters, inputClusterMask, layerClustersTimes, layer_clusters_tiles, seeding_regions, tfSession_); - myAlgo_->makeTracksters(input, *result, seedToTrackstersAssociation); + + // TODO(mmatthew): Delete if conditions once correct function definition for KF is found + if(itername_ == "KF"){ + myAlgo_->makeTrajectories(input,*kfhits, *prophits, *abs_fail); + } else { + myAlgo_->makeTracksters(input, *result, seedToTrackstersAssociation); + } + } // Now update the global mask and put it into the event output_mask->reserve(original_layerclusters_mask.size()); @@ -198,4 +224,7 @@ void TrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) { evt.put(std::move(result)); evt.put(std::move(output_mask)); -} + evt.put(std::move(abs_fail),"Abs Fail"); + evt.put(std::move(kfhits),"KFHits"); + evt.put(std::move(prophits),"PropHits"); + } diff --git a/RecoHGCal/TICL/python/KFStep_cff.py b/RecoHGCal/TICL/python/KFStep_cff.py new file mode 100644 index 0000000000000..631daeb46f6ad --- /dev/null +++ b/RecoHGCal/TICL/python/KFStep_cff.py @@ -0,0 +1,53 @@ +import FWCore.ParameterSet.Config as cms + +from RecoHGCal.TICL.TICLSeedingRegions_cff import ticlSeedingTrk +from RecoHGCal.TICL.trackstersProducer_cfi import trackstersProducer as _trackstersProducer +from RecoHGCal.TICL.filteredLayerClustersProducer_cfi import filteredLayerClustersProducer as _filteredLayerClustersProducer +from RecoHGCal.TICL.multiClustersFromTrackstersProducer_cfi import multiClustersFromTrackstersProducer as _multiClustersFromTrackstersProducer +from RecoHGCal.TICL.ticlLayerTileProducer_cfi import ticlLayerTileProducer + +# CLUSTER FILTERING/MASKING + +filteredLayerClustersKF = _filteredLayerClustersProducer.clone( + clusterFilter = "ClusterFilterByAlgoAndSize", + min_cluster_size = 2, # inclusive + algo_number = 8, + iteration_label = "KF", +) + +# PATTERN RECOGNITION + +ticlTrackstersKF = _trackstersProducer.clone( + filtered_mask = "filteredLayerClustersKF:KF", + layer_clusters_tiles = "ticlRecHitTile:Test", + seeding_regions = "ticlSeedingTrk", + itername = "KF", + patternRecognitionBy = "KF", + pluginPatternRecognitionByKF = dict ( + radlen = [0.562719015399016,0.9551903779795146,1.043875789035635,0.8920189121952249,0.9712515382761755, + 0.9638214213266212,1.0131600728257035,0.931610084339171,0.9762915187242845,0.9673823352145134, + 0.9846489838011926,0.9687667855140328,0.9718771132058326,0.9676763827273792,0.9755766684867133, + 0.968786167714319,1.041570725190613,0.905510318469831,1.5448704774827073,0.9671013911948058, + 1.5527918490549046,0.9644685854108038,1.5537767590772382,0.9681776397486216,1.5531742202624652, + 0.9692846975312643,2.8275418030617567,3.003105180494033,3.005570502299292,3.0076357163745446, + 3.005138681398646,3.0106898519208087,3.005449098667002,3.007544959496376,3.0033570538999084, + 3.005449880885637,3.0056943107449103,4.158164083557213,4.1459435055013865,4.087756410389329, + 4.14190322796051,4.143581858605728,4.135287193893763,4.093154468149641,4.1386581126514566, + 4.145881343241816,4.145663423464052], + xi = [0.0007556128308446044,0.0007601810242278396,0.0007383347392914731,0.0007954300436349306,0.0007079974303217379, + 0.0007587747875111608,0.0007506686955642467,0.000777442404205796,0.0007223198806925819,0.0007536986711079965, + 0.0007687577424509902,0.0007640633376218576,0.0007307246149741174,0.0007518274635026171,0.0007675709426774015, + 0.00075807841348618,0.000743616079227404,0.0007936746547956953,0.0005974640645738568,0.0007508915190685968, + 0.0006382593210233038,0.0007517838017659952,0.0006397287020925379,0.0007579429706626264,0.0006333046151356415, + 0.0007637633509097611,0.0010337707504477352,0.0010108776039545466,0.0010104712217333819,0.0010106368784497104, + 0.0010106003723634636,0.0010158675945647706,0.001010809007077468,0.001009561281210271,0.0010106988945190448, + 0.0010103715066888692,0.0010106495620045056,0.0010210702627516468,0.0010169851234559488,0.0010010856800569403, + 0.0010155575635282914,0.0010146829805806501,0.0010150368006945165,0.0009992234164373278,0.0010145096108228593, + 0.0010126783833862823,0.0010125427722475087], + rescale = 100 + ) +) + +ticlKFStepTask = cms.Task(ticlSeedingTrk + ,filteredLayerClustersKF + ,ticlTrackstersKF) diff --git a/RecoHGCal/TICL/python/iterativeTICL_cff.py b/RecoHGCal/TICL/python/iterativeTICL_cff.py index 69c9989ca8955..26f704027448d 100644 --- a/RecoHGCal/TICL/python/iterativeTICL_cff.py +++ b/RecoHGCal/TICL/python/iterativeTICL_cff.py @@ -3,6 +3,7 @@ from RecoHGCal.TICL.FastJetStep_cff import * from RecoHGCal.TICL.CLUE3DHighStep_cff import * from RecoHGCal.TICL.CLUE3DLowStep_cff import * +from RecoHGCal.TICL.KFStep_cff import * from RecoHGCal.TICL.MIPStep_cff import * from RecoHGCal.TICL.TrkEMStep_cff import * from RecoHGCal.TICL.TrkStep_cff import * @@ -27,6 +28,13 @@ ticlCLUE3DHighStepTask ) +from Configuration.ProcessModifiers.kf_cff import kf +ticlRecHitTile = ticlLayerTileProducer.clone( + isLC = False +) +ticlRecHitTileTask = cms.Task(ticlRecHitTile) +kf.toModify(ticlIterationsTask, func=lambda x : x.add(ticlRecHitTileTask, ticlKFStepTask)) + from Configuration.ProcessModifiers.clue3D_cff import clue3D clue3D.toModify(ticlIterationsTask, func=lambda x : x.add(ticlCLUE3DHighStepTask,ticlCLUE3DLowStepTask)) diff --git a/RecoHGCal/TICL/test/TICLDebugger.log b/RecoHGCal/TICL/test/TICLDebugger.log new file mode 100644 index 0000000000000..0ea155966063a --- /dev/null +++ b/RecoHGCal/TICL/test/TICLDebugger.log @@ -0,0 +1,18 @@ +01-Jun-2022 11:43:57 CEST Initiating request to open file file:/afs/cern.ch/user/m/mmatthew/CMSSW_12_3_0_pre5/src/hgcal_private/production/Eta_16/singlemuon_flatEGun_hgcalCenter/step3/step3_singlemuon_e100GeV_nopu.root +01-Jun-2022 11:44:21 CEST Successfully opened file file:/afs/cern.ch/user/m/mmatthew/CMSSW_12_3_0_pre5/src/hgcal_private/production/Eta_16/singlemuon_flatEGun_hgcalCenter/step3/step3_singlemuon_e100GeV_nopu.root + +TrksIdx: 0 + bary: (10.6407,-194.901,463.461) baryEta: 1.59954 baryPhi: -1.51625 + raw_energy: 1.13398 raw_em_energy: 0 + raw_pt: 0.430689 raw_em_pt: 0 + seedIdx: 0 + Probs: (h):0.6464 (gam):0.2913 (e):0.0604 (pi0):0.0004 (h0):0.0004 (!):0.0004 (mu):0.0004 (?):0.0003 + + time: -99+/--1 vertices: 11 average usage: 1 + link connections: 0 + Seeding Track: + p: 100.642 pt: 39.0459 charge: -1 eta: 1.60003 outerEta: 1.6 phi: -1.53695 outerPhi: -1.51002 + Best CaloParticles Matches: + 0(0.000432394):13 simCl size:1 energy:99.9979 pt:38.7968 momentum:(1.30924,-38.7747,92.1649) + +01-Jun-2022 11:44:48 CEST Closed file file:/afs/cern.ch/user/m/mmatthew/CMSSW_12_3_0_pre5/src/hgcal_private/production/Eta_16/singlemuon_flatEGun_hgcalCenter/step3/step3_singlemuon_e100GeV_nopu.root diff --git a/RecoHGCal/TICL/test/ticlDebugger_cfg.py b/RecoHGCal/TICL/test/ticlDebugger_cfg.py index 44582be049cd3..707fc561b5594 100644 --- a/RecoHGCal/TICL/test/ticlDebugger_cfg.py +++ b/RecoHGCal/TICL/test/ticlDebugger_cfg.py @@ -13,7 +13,7 @@ process.source = cms.Source("PoolSource", # replace 'myfile.root' with the source file you want to use fileNames = cms.untracked.vstring( - 'file:step3.root' + 'file:/afs/cern.ch/user/m/mmatthew/CMSSW_12_3_0_pre5/src/hgcal_private/production/Eta_16/singlemuon_flatEGun_hgcalCenter/step3/step3_singlemuon_e100GeV_nopu.root' ) )