diff --git a/Configuration/PyReleaseValidation/python/relval_standard.py b/Configuration/PyReleaseValidation/python/relval_standard.py index 783b425e0c657..a6b085613b8b6 100644 --- a/Configuration/PyReleaseValidation/python/relval_standard.py +++ b/Configuration/PyReleaseValidation/python/relval_standard.py @@ -573,6 +573,10 @@ ### run3-2023 (2023 HI data RawPrime with re-HLT) workflows[142.0] = ['',['RunHIPhysicsRawPrime2023A','HLTDR3_HI2023ARawprime','RECOHIRUN3_reHLT_2023','HARVESTRUN3_HI2023A']] +### run3-2024 (2024 HI UPC data) +workflows[142.901] = ['',['RunUPC2023','RECODR3_UPC','HARVESTDPROMPTR3']] +workflows[142.902] = ['',['RunUPC2023','RECODR3_HIN','HARVESTDPROMPTR3']] + ### fastsim ### workflows[5.1] = ['TTbarFS', ['TTbarFS','HARVESTFS']] workflows[5.2] = ['SingleMuPt10FS', ['SingleMuPt10FS','HARVESTFS']] diff --git a/Configuration/PyReleaseValidation/python/relval_steps.py b/Configuration/PyReleaseValidation/python/relval_steps.py index 0769ec2bd8adf..2addafed96d66 100644 --- a/Configuration/PyReleaseValidation/python/relval_steps.py +++ b/Configuration/PyReleaseValidation/python/relval_steps.py @@ -2689,7 +2689,9 @@ def lhegensim2018ml(fragment,howMuch): steps['RECODR3_reHLT_2023B']=merge([{'--conditions':'auto:run3_data_prompt_relval', '--hltProcess':'reHLT'},steps['RECODR3']]) steps['RECODR3_2023_HIN']=merge([{'--conditions':'auto:run3_data_prompt', '-s':'RAW2DIGI,L1Reco,RECO,DQM:@commonFakeHLT+@standardDQMFakeHLT', '--repacked':'', '-n':1000},steps['RECODR3_2023']]) -steps['RECODR3_2023_UPC']=merge([{'--era':'Run3_2023_UPC', '--conditions':'132X_dataRun3_Prompt_HI_LowPtPhotonReg_v2'},steps['RECODR3_2023_HIN']]) +steps['RECODR3_2023_UPC']=merge([{'--era':'Run3_2023_UPC'},steps['RECODR3_2023_HIN']]) +steps['RECODR3_HIN']=merge([{'--conditions':'auto:run3_data_prompt', '-s':'RAW2DIGI,L1Reco,RECO,DQM:@commonFakeHLT+@standardDQMFakeHLT', '--repacked':'', '-n':1000},steps['RECODR3']]) +steps['RECODR3_UPC']=merge([{'--era':'Run3_UPC'},steps['RECODR3_HIN']]) steps['RECODR3Splash']=merge([{'-n': 2, '-s': 'RAW2DIGI,L1Reco,RECO,PAT,ALCA:SiStripCalZeroBias+SiStripCalMinBias+TkAlMinBias+EcalESAlign,DQM:@standardDQMFakeHLT+@miniAODDQM' diff --git a/DataFormats/TrackReco/interface/DeDxHit.h b/DataFormats/TrackReco/interface/DeDxHit.h index 29f58717a1b91..8f0d3771ad44b 100644 --- a/DataFormats/TrackReco/interface/DeDxHit.h +++ b/DataFormats/TrackReco/interface/DeDxHit.h @@ -12,8 +12,8 @@ namespace reco { public: DeDxHit() {} - DeDxHit(float ch, float mom, float len, uint32_t rawDetId) - : m_charge(ch), m_momentum(mom), m_pathLength(len), m_rawDetId(rawDetId) {} + DeDxHit(float ch, float mom, float len, uint32_t rawDetId, float err = 0) + : m_charge(ch), m_momentum(mom), m_pathLength(len), m_rawDetId(rawDetId), m_error(err) {} /// Return the angle and thick normalized, calibrated energy release float charge() const { return m_charge; } @@ -27,6 +27,9 @@ namespace reco { /// Return the rawDetId uint32_t rawDetId() const { return m_rawDetId; } + /// Return the error of the energy release + float error() const { return m_error; } + bool operator<(const DeDxHit &other) const { return m_charge < other.m_charge; } private: @@ -36,6 +39,7 @@ namespace reco { float m_momentum; float m_pathLength; uint32_t m_rawDetId; + float m_error; }; typedef std::vector DeDxHitCollection; diff --git a/DataFormats/TrackReco/interface/DeDxHitInfo.h b/DataFormats/TrackReco/interface/DeDxHitInfo.h index 1a7a115c2bed9..43781bf587694 100644 --- a/DataFormats/TrackReco/interface/DeDxHitInfo.h +++ b/DataFormats/TrackReco/interface/DeDxHitInfo.h @@ -16,13 +16,15 @@ namespace reco { class DeDxHitInfoContainer { public: DeDxHitInfoContainer() : charge_(0.0f), pathlength_(0.0f) {} - DeDxHitInfoContainer(const float charge, const float pathlength, const DetId& detId, const LocalPoint& pos) - : charge_(charge), pathlength_(pathlength), detId_(detId), pos_(pos) {} + DeDxHitInfoContainer( + const float charge, const float pathlength, const DetId& detId, const LocalPoint& pos, const uint8_t& type) + : charge_(charge), pathlength_(pathlength), detId_(detId), pos_(pos), type_(type) {} float charge() const { return charge_; } float pathlength() const { return pathlength_; } const DetId& detId() const { return detId_; } const LocalPoint& pos() const { return pos_; } + const uint8_t& type() const { return type_; } private: //! total cluster charge @@ -32,8 +34,10 @@ namespace reco { DetId detId_; //! hit position LocalPoint pos_; + uint8_t type_; }; + static constexpr int Complete = 0, Compatible = 1, Calibration = 2; typedef std::vector DeDxHitInfoContainerCollection; public: @@ -43,6 +47,7 @@ namespace reco { float pathlength(size_t i) const { return infos_[i].pathlength(); } DetId detId(size_t i) const { return infos_[i].detId(); } const LocalPoint pos(size_t i) const { return infos_[i].pos(); } + const uint8_t type(size_t i) const { return infos_[i].type(); } const SiPixelCluster* pixelCluster(size_t i) const { size_t P = 0; bool isPixel = false; @@ -90,16 +95,18 @@ namespace reco { const float pathlength, const DetId& detId, const LocalPoint& pos, + const uint8_t& type, const SiStripCluster& stripCluster) { - infos_.push_back(DeDxHitInfoContainer(charge, pathlength, detId, pos)); + infos_.push_back(DeDxHitInfoContainer(charge, pathlength, detId, pos, type)); stripClusters_.push_back(stripCluster); } void addHit(const float charge, const float pathlength, const DetId& detId, const LocalPoint& pos, + const uint8_t& type, const SiPixelCluster& pixelCluster) { - infos_.push_back(DeDxHitInfoContainer(charge, pathlength, detId, pos)); + infos_.push_back(DeDxHitInfoContainer(charge, pathlength, detId, pos, type)); pixelClusters_.push_back(pixelCluster); } diff --git a/DataFormats/TrackReco/src/classes_def.xml b/DataFormats/TrackReco/src/classes_def.xml index 1c74d269b5b4a..e4fbf85b52503 100644 --- a/DataFormats/TrackReco/src/classes_def.xml +++ b/DataFormats/TrackReco/src/classes_def.xml @@ -398,11 +398,13 @@ - + + + @@ -490,8 +493,9 @@ - + + diff --git a/RecoTracker/Configuration/python/RecoTracker_EventContent_cff.py b/RecoTracker/Configuration/python/RecoTracker_EventContent_cff.py index 5137c1f736c6e..eef499b1016bc 100644 --- a/RecoTracker/Configuration/python/RecoTracker_EventContent_cff.py +++ b/RecoTracker/Configuration/python/RecoTracker_EventContent_cff.py @@ -23,6 +23,12 @@ (pp_on_AA | run3_upc).toModify( RecoTrackerAOD.outputCommands, func=lambda outputCommands: outputCommands.extend(['keep recoTracks_hiConformalPixelTracks_*_*']) ) +run3_upc.toModify( RecoTrackerAOD.outputCommands, + func=lambda outputCommands: outputCommands.extend([ + 'keep *_dedxPixelLikelihood_*_*', + 'keep *_dedxStripLikelihood_*_*', + 'keep *_dedxAllLikelihood_*_*' + ])) #RECO content RecoTrackerRECO = cms.PSet( outputCommands = cms.untracked.vstring( diff --git a/RecoTracker/DeDx/BuildFile.xml b/RecoTracker/DeDx/BuildFile.xml index b3b406d3b7a9a..e7b34eb922974 100644 --- a/RecoTracker/DeDx/BuildFile.xml +++ b/RecoTracker/DeDx/BuildFile.xml @@ -8,6 +8,8 @@ + + diff --git a/RecoTracker/DeDx/interface/LikelihoodFitDeDxEstimator.h b/RecoTracker/DeDx/interface/LikelihoodFitDeDxEstimator.h new file mode 100644 index 0000000000000..67b741dff37e9 --- /dev/null +++ b/RecoTracker/DeDx/interface/LikelihoodFitDeDxEstimator.h @@ -0,0 +1,151 @@ +#ifndef RecoTracker_DeDx_LikelihoodFitDeDxEstimator_h +#define RecoTracker_DeDx_LikelihoodFitDeDxEstimator_h + +#include "RecoTracker/DeDx/interface/BaseDeDxEstimator.h" +#include "DataFormats/TrackReco/interface/DeDxHit.h" + +class LikelihoodFitDeDxEstimator : public BaseDeDxEstimator { +public: + LikelihoodFitDeDxEstimator(const edm::ParameterSet& iConfig){}; + + std::pair dedx(const reco::DeDxHitCollection& Hits) override { + if (Hits.empty()) + return {0., 0.}; + // compute original + std::array value; + const auto& chi2 = estimate(Hits, value); + // try to remove lowest dE/dx measurement + const auto& n = Hits.size(); + if (n >= 3 && (chi2 > 1.3 * n + 4 * std::sqrt(1.3 * n))) { + auto hs = Hits; + hs.erase(std::min_element(hs.begin(), hs.end())); + // if got better, accept + std::array v; + if (estimate(hs, v) < chi2 - 12) + value = v; + } + return {value[0], std::sqrt(value[1])}; + } + +private: + void calculate_wrt_epsilon(const reco::DeDxHit&, const double&, std::array&); + void functionEpsilon(const reco::DeDxHitCollection&, const double&, std::array&); + double minimizeAllSaturated(const reco::DeDxHitCollection&, std::array&); + double newtonMethodEpsilon(const reco::DeDxHitCollection&, std::array&); + double estimate(const reco::DeDxHitCollection&, std::array&); +}; + +/*****************************************************************************/ +void LikelihoodFitDeDxEstimator::calculate_wrt_epsilon(const reco::DeDxHit& h, + const double& epsilon, + std::array& value) { + const auto& ls = h.pathLength(); + const auto& sn = h.error(); // energy sigma + const auto y = h.charge() * ls; // = g * y + const auto sD = 2.E-3 + 0.095 * y; + const auto ss = sD * sD + sn * sn; + const auto s = std::sqrt(ss); + const auto delta = epsilon * ls; + const auto dy = delta - y; + constexpr double nu(0.65); + + // calculate derivatives with respect to delta + std::array val{{0.}}; + if (h.rawDetId() == 0) { // normal + if (dy < -nu * s) { + val[0] = -2. * nu * dy / s - nu * nu; + val[1] = -2. * nu / s; + val[2] = 0.; + } else { + val[0] = dy * dy / ss; + val[1] = 2. * dy / ss; + val[2] = 2. / ss; + } + } else { // saturated + if (dy < s) { + val[0] = -dy / s + 1.; + val[1] = -1. / s; + val[2] = 0.; + } else { + val[0] = 0.; + val[1] = 0.; + val[2] = 0.; + } + } + + // d/d delta -> d/d epsilon + val[1] *= ls; + val[2] *= ls * ls; + + // sum + for (size_t k = 0; k < value.size(); k++) + value[k] += val[k]; +} + +/*****************************************************************************/ +void LikelihoodFitDeDxEstimator::functionEpsilon(const reco::DeDxHitCollection& Hits, + const double& epsilon, + std::array& val) { + val = {{0, 0, 0}}; + for (const auto& h : Hits) + calculate_wrt_epsilon(h, epsilon, val); +} + +/*****************************************************************************/ +double LikelihoodFitDeDxEstimator::minimizeAllSaturated(const reco::DeDxHitCollection& Hits, + std::array& value) { + int nStep(0); + double par(3.0); // input MeV/cm + + std::array val{{0}}; + do { + functionEpsilon(Hits, par, val); + if (val[1] != 0) + par += -val[0] / val[1]; + nStep++; + } while (val[0] > 1e-3 && val[1] != 0 && nStep < 10); + + value[0] = par * 1.1; + value[1] = par * par * 0.01; + + return val[0]; +} + +/*****************************************************************************/ +double LikelihoodFitDeDxEstimator::newtonMethodEpsilon(const reco::DeDxHitCollection& Hits, + std::array& value) { + int nStep(0); + double par(3.0); // input MeV/cm + double dpar(0); + + std::array val{{0}}; + do { + functionEpsilon(Hits, par, val); + if (val[2] != 0.) + dpar = -val[1] / std::abs(val[2]); + else + dpar = 1.; // step up, for epsilon + if (par + dpar > 0) + par += dpar; // ok + else + par /= 2.; // half + nStep++; + } while (std::abs(dpar) > 1e-3 && nStep < 50); + + value[0] = par; + value[1] = 2. / val[2]; + + return val[0]; +} + +/*****************************************************************************/ +double LikelihoodFitDeDxEstimator::estimate(const reco::DeDxHitCollection& Hits, std::array& value) { + // use newton method if at least one hit is not saturated + for (const auto& h : Hits) + if (h.rawDetId() == 0) + return newtonMethodEpsilon(Hits, value); + // else use minimize all saturated + return minimizeAllSaturated(Hits, value); +} + +#endif diff --git a/RecoTracker/DeDx/plugins/DeDxEstimatorProducer.cc b/RecoTracker/DeDx/plugins/DeDxEstimatorProducer.cc index e046246b22e09..7aae98746aeb2 100644 --- a/RecoTracker/DeDx/plugins/DeDxEstimatorProducer.cc +++ b/RecoTracker/DeDx/plugins/DeDxEstimatorProducer.cc @@ -30,6 +30,9 @@ void DeDxEstimatorProducer::fillDescriptions(edm::ConfigurationDescriptions& des edm::ParameterSetDescription desc; desc.add("estimator", "generic"); desc.add("tracks", edm::InputTag("generalTracks")); + desc.add("UseDeDxHits", false); + desc.add("pixelDeDxHits", {}); + desc.add("stripDeDxHits", {}); desc.add("UsePixel", false); desc.add("UseStrip", true); desc.add("MeVperADCPixel", 3.61e-06); @@ -62,6 +65,8 @@ DeDxEstimatorProducer::DeDxEstimatorProducer(const edm::ParameterSet& iConfig) m_estimator = std::make_unique(iConfig); else if (estimatorName == "unbinnedFit") m_estimator = std::make_unique(iConfig); + else if (estimatorName == "likelihoodFit") + m_estimator = std::make_unique(iConfig); else if (estimatorName == "productDiscrim") m_estimator = std::make_unique(iConfig, cCollector); else if (estimatorName == "btagDiscrim") @@ -76,6 +81,10 @@ DeDxEstimatorProducer::DeDxEstimatorProducer(const edm::ParameterSet& iConfig) m_tracksTag = consumes(iConfig.getParameter("tracks")); + useDeDxHits = iConfig.getParameter("UseDeDxHits"); + m_pixelDeDxHitsTag = consumes(iConfig.getParameter("pixelDeDxHits")); + m_stripDeDxHitsTag = consumes(iConfig.getParameter("stripDeDxHits")); + usePixel = iConfig.getParameter("UsePixel"); useStrip = iConfig.getParameter("UseStrip"); meVperADCPixel = iConfig.getParameter("MeVperADCPixel"); @@ -110,27 +119,39 @@ void DeDxEstimatorProducer::produce(edm::Event& iEvent, const edm::EventSetup& i edm::Handle trackCollectionHandle; iEvent.getByToken(m_tracksTag, trackCollectionHandle); + const auto& pixelDeDxHits = iEvent.getHandle(m_pixelDeDxHitsTag); + const auto& stripDeDxHits = iEvent.getHandle(m_stripDeDxHitsTag); std::vector dedxEstimate(trackCollectionHandle->size()); for (unsigned int j = 0; j < trackCollectionHandle->size(); j++) { - const reco::TrackRef track = reco::TrackRef(trackCollectionHandle.product(), j); + const reco::TrackRef track = reco::TrackRef(trackCollectionHandle, j); int NClusterSaturating = 0; DeDxHitCollection dedxHits; - auto const& trajParams = track->extra()->trajParams(); - assert(trajParams.size() == track->recHitsSize()); - auto hb = track->recHitsBegin(); dedxHits.reserve(track->recHitsSize() / 2); - for (unsigned int h = 0; h < track->recHitsSize(); h++) { - auto recHit = *(hb + h); - if (!trackerHitRTTI::isFromDet(*recHit)) - continue; - - auto trackDirection = trajParams[h].direction(); - float cosine = trackDirection.z() / trackDirection.mag(); - processHit(recHit, track->p(), cosine, dedxHits, NClusterSaturating); + if (useDeDxHits) { + if (usePixel) + dedxHits = (*pixelDeDxHits)[track]; + if (useStrip) { + const auto& hits = (*stripDeDxHits)[track]; + dedxHits.insert(dedxHits.end(), hits.begin(), hits.end()); + } + } else { + auto const& trajParams = track->extra()->trajParams(); + assert(trajParams.size() == track->recHitsSize()); + auto hb = track->recHitsBegin(); + dedxHits.reserve(track->recHitsSize() / 2); + for (unsigned int h = 0; h < track->recHitsSize(); h++) { + auto recHit = *(hb + h); + if (!trackerHitRTTI::isFromDet(*recHit)) + continue; + + auto trackDirection = trajParams[h].direction(); + float cosine = trackDirection.z() / trackDirection.mag(); + processHit(recHit, track->p(), cosine, dedxHits, NClusterSaturating); + } } sort(dedxHits.begin(), dedxHits.end(), less()); diff --git a/RecoTracker/DeDx/plugins/DeDxEstimatorProducer.h b/RecoTracker/DeDx/plugins/DeDxEstimatorProducer.h index a77e88a84eaf1..b6029c587105b 100644 --- a/RecoTracker/DeDx/plugins/DeDxEstimatorProducer.h +++ b/RecoTracker/DeDx/plugins/DeDxEstimatorProducer.h @@ -32,10 +32,10 @@ #include "RecoTracker/DeDx/interface/SmirnovDeDxDiscriminator.h" #include "RecoTracker/DeDx/interface/ASmirnovDeDxDiscriminator.h" #include "RecoTracker/DeDx/interface/BTagLikeDeDxDiscriminator.h" +#include "RecoTracker/DeDx/interface/LikelihoodFitDeDxEstimator.h" #include "RecoTracker/DeDx/interface/DeDxTools.h" -#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h" #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h" // @@ -63,7 +63,10 @@ class DeDxEstimatorProducer : public edm::stream::EDProducer<> { std::unique_ptr m_estimator; edm::EDGetTokenT m_tracksTag; + edm::EDGetTokenT m_pixelDeDxHitsTag; + edm::EDGetTokenT m_stripDeDxHitsTag; + bool useDeDxHits; bool usePixel; bool useStrip; float meVperADCPixel; diff --git a/RecoTracker/DeDx/plugins/DeDxHitCalibrator.cc b/RecoTracker/DeDx/plugins/DeDxHitCalibrator.cc new file mode 100644 index 0000000000000..0953b0146e462 --- /dev/null +++ b/RecoTracker/DeDx/plugins/DeDxHitCalibrator.cc @@ -0,0 +1,427 @@ +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/FileInPath.h" + +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/TrackReco/interface/TrackDeDxHits.h" +#include "DataFormats/TrackReco/interface/DeDxHit.h" +#include "DataFormats/TrackReco/interface/DeDxHitInfo.h" +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "DataFormats/SiStripCommon/interface/ConstantsForHardwareSystems.h" +#include "CondFormats/PhysicsToolsObjects/interface/DeDxCalibration.h" +#include "CondFormats/DataRecord/interface/DeDxCalibrationRcd.h" +#include "CondFormats/SiPixelObjects/interface/PixelIndices.h" + +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "Geometry/Records/interface/TrackerTopologyRcd.h" +#include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h" +#include "CalibTracker/SiPixelESProducers/interface/SiPixelGainCalibrationOfflineService.h" + +#include + +class DeDxHitCalibrator : public edm::stream::EDProducer<> { +public: + static constexpr int kIsNormal = 0, kIsBelow = 1, kIsOver = 2; + static constexpr int PXB = 0, PXF = 1, TIB = 2, TID = 3, TOB = 4, TECThin = 5, TECThick = 6; + typedef std::pair ChipId; + + explicit DeDxHitCalibrator(const edm::ParameterSet&); + ~DeDxHitCalibrator() override = default; + static void fillDescriptions(edm::ConfigurationDescriptions&); + +private: + void beginRun(edm::Run const&, const edm::EventSetup&) override; + void produce(edm::Event&, const edm::EventSetup&) override; + + int getDetId(const DetId&, const float&); + float correctEnergy(const float&, const ChipId&); + void processHitInfo(const reco::DeDxHitInfo&, + const float& trackMomentum, + reco::DeDxHitCollection&, + reco::DeDxHitCollection&); + + double getChi2(const std::vector&, const std::vector >&, const double&, const double&); + void getAlphaBeta(const std::vector&, + const std::vector >&, + CLHEP::HepMatrix&, + CLHEP::HepVector&, + const std::vector&, + const double&, + const double&); + std::pair fitStripCluster(const std::vector >&, const double&, const double&); + + const bool applyGain_; + const double MeVPerElectron_; + const int VCaltoElectronGain_, VCaltoElectronGain_L1_, VCaltoElectronOffset_, VCaltoElectronOffset_L1_; + const int pixelSaturationThr_; + const edm::EDGetTokenT tracksToken_; + const edm::EDGetTokenT dedxHitInfoToken_; + const edm::ESGetToken dedxCalibToken_; + const edm::ESGetToken tkGeomToken_; + const edm::ESGetToken tkTopoToken_; + + SiPixelGainCalibrationOfflineService pixelCalib_; + edm::ESHandle dedxCalib_; + edm::ESHandle tkGeom_; + edm::ESHandle tkTopo_; +}; + +DeDxHitCalibrator::DeDxHitCalibrator(const edm::ParameterSet& iConfig) + : applyGain_(iConfig.getParameter("applyGain")), + MeVPerElectron_(iConfig.getParameter("MeVPerElectron")), + VCaltoElectronGain_(iConfig.getParameter("VCaltoElectronGain")), + VCaltoElectronGain_L1_(iConfig.getParameter("VCaltoElectronGain_L1")), + VCaltoElectronOffset_(iConfig.getParameter("VCaltoElectronOffset")), + VCaltoElectronOffset_L1_(iConfig.getParameter("VCaltoElectronOffset_L1")), + pixelSaturationThr_(iConfig.getParameter("pixelSaturationThr")), + tracksToken_(consumes(iConfig.getParameter("trackProducer"))), + dedxHitInfoToken_(consumes(iConfig.getParameter("dedxHitInfo"))), + dedxCalibToken_(esConsumes()), + tkGeomToken_(esConsumes()), + tkTopoToken_(esConsumes()), + pixelCalib_(iConfig, consumesCollector()) { + produces("PixelHits"); + produces("StripHits"); +} + +void DeDxHitCalibrator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("applyGain", true); + desc.add("MeVPerElectron", 3.61e-06); + desc.add("VCaltoElectronGain", 65); + desc.add("VCaltoElectronGain_L1", 65); + desc.add("VCaltoElectronOffset", -414); + desc.add("VCaltoElectronOffset_L1", -414); + desc.add("pixelSaturationThr", 254); + desc.add("trackProducer", edm::InputTag("generalTracks")); + desc.add("dedxHitInfo", edm::InputTag("dedxHitInfo")); + descriptions.add("dedxHitCalibrator", desc); +} + +void DeDxHitCalibrator::beginRun(edm::Run const&, const edm::EventSetup& iSetup) { + dedxCalib_ = iSetup.getHandle(dedxCalibToken_); + tkGeom_ = iSetup.getHandle(tkGeomToken_); + tkTopo_ = iSetup.getHandle(tkTopoToken_); +} + +void DeDxHitCalibrator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + const auto& tracks = iEvent.getHandle(tracksToken_); + const auto& dedxHitInfo = iEvent.get(dedxHitInfoToken_); + pixelCalib_.setESObjects(iSetup); + + // creates the output collection + auto pixelTrackDeDxHitAss = std::make_unique(reco::TrackRefProd(tracks)); + auto stripTrackDeDxHitAss = std::make_unique(reco::TrackRefProd(tracks)); + + for (size_t i = 0; i < tracks->size(); i++) { + const auto& track = reco::TrackRef(tracks, i); + const auto& dedxHits = dedxHitInfo[track]; + if (dedxHits.isNull()) + continue; + reco::DeDxHitCollection pixelHits, stripHits; + processHitInfo(*dedxHits, track->p(), pixelHits, stripHits); + pixelTrackDeDxHitAss->setValue(i, pixelHits); + stripTrackDeDxHitAss->setValue(i, stripHits); + } + + iEvent.put(std::move(pixelTrackDeDxHitAss), "PixelHits"); + iEvent.put(std::move(stripTrackDeDxHitAss), "StripHits"); +} + +/*****************************************************************************/ +void DeDxHitCalibrator::processHitInfo(const reco::DeDxHitInfo& info, + const float& trackMomentum, + reco::DeDxHitCollection& pixelHits, + reco::DeDxHitCollection& stripHits) { + // Energy loss parametrization from https://doi.org/10.1016/j.nima.2012.06.064 + static constexpr double a = 0.07; //energy loss factor in silicon + static constexpr double l0 = 450e-4; //reference path length in um + + for (size_t i = 0; i < info.size(); i++) { + // Require hits to be complete and compatible + const auto& type = info.type(i); + if (!(type & (1 << reco::DeDxHitInfo::Complete)) || !(type & (1 << reco::DeDxHitInfo::Compatible))) + continue; + + // Effective path length + const auto& pl = info.pathlength(i); + const auto pathLength = pl * (1. + a * std::log(pl / l0)); + const DetId& detId = info.detId(i); + + // Strip + if (const auto& stripCluster = info.stripCluster(i)) { + const auto& thickness = + dynamic_cast(tkGeom_->idToDet(detId))->surface().bounds().thickness(); + const auto& det = getDetId(detId, thickness); + const auto& chipId = ChipId(detId, stripCluster->barycenter() / sistrip::STRIPS_PER_APV); + // Fill measured strip deposits + const auto& thr = dedxCalib_->thr()[det]; + std::vector > b; + b.reserve(stripCluster->amplitudes().size() + 2); + b.emplace_back(thr, kIsBelow); + for (const auto& adc : stripCluster->amplitudes()) { + if (adc > 253) + b.emplace_back(254., kIsOver); + else + b.emplace_back(adc + 0.5, kIsNormal); + } + b.emplace_back(thr, kIsBelow); + // Fit + const auto& result = fitStripCluster(b, dedxCalib_->alpha()[det], 1. / dedxCalib_->sigma()[det]); + // ADC -> e- -> MeV + const auto& energy = correctEnergy(result.first * sistrip::MeVperADCStrip, chipId); + const auto& esigma = correctEnergy(result.second * sistrip::MeVperADCStrip, chipId); + // Compute charge + const auto& charge = pathLength != 0. ? energy / pathLength : 0.; + stripHits.emplace_back(charge, trackMomentum, pathLength, 0, esigma); + } + + // Pixel + if (const auto& pixelCluster = info.pixelCluster(i)) { + const auto& thickness = + dynamic_cast(tkGeom_->idToDet(detId))->surface().bounds().thickness(); + const auto& det = getDetId(detId, thickness); + const auto& chipId = + ChipId(detId, (int(pixelCluster->x() / ROCSizeInX) << 3) + int(pixelCluster->y() / ROCSizeInY)); + // Collect adc + double delta(0); + bool isSaturated(false); + for (size_t j = 0; j < pixelCluster->pixelADC().size(); j++) { + const auto& elec = pixelCluster->pixelADC()[j]; + delta += elec * MeVPerElectron_; + if (isSaturated) + continue; + // ( pos , (x , y) ) + const auto& row = pixelCluster->minPixelRow() + pixelCluster->pixelOffset()[2 * j]; + const auto& col = pixelCluster->minPixelCol() + pixelCluster->pixelOffset()[2 * j + 1]; + // Go back to adc + const auto& DBgain = pixelCalib_.getGain(detId, col, row); + const auto& DBpedestal = pixelCalib_.getPedestal(detId, col, row); + if (elec == std::numeric_limits::max()) + isSaturated = true; + else if (DBgain > 0.) { + double vcal; + const auto& theLayer = (detId.subdetId() == 1) ? tkTopo_->pxbLayer(detId) : 0; + if (theLayer == 1) + vcal = (elec - VCaltoElectronOffset_L1_) / VCaltoElectronGain_L1_; + else + vcal = (elec - VCaltoElectronOffset_) / VCaltoElectronGain_; + const auto adc = std::round(DBpedestal + vcal / DBgain); + if (adc > pixelSaturationThr_) + isSaturated = true; + } + } + // Compute energy + const auto& energy = correctEnergy(delta, chipId); + if (det != PXF && energy <= 50e-3) + continue; + // Estimate sigma for cluster + const unsigned char& nChannels = pixelCluster->pixelADC().size(); + const auto esigma = 10e-3 * std::sqrt(nChannels); + // Compute charge + const auto& charge = pathLength != 0. ? energy / pathLength : 0.; + pixelHits.emplace_back(charge, trackMomentum, pathLength, isSaturated, esigma); + } + } +} + +/*****************************************************************************/ +int DeDxHitCalibrator::getDetId(const DetId& id, const float& thickness) { + const auto& subdet = id.subdetId() - 1; + if (subdet == TECThin && thickness > 400e-4) + return TECThick; + return subdet; +} + +/*****************************************************************************/ +float DeDxHitCalibrator::correctEnergy(const float& energy, const ChipId& detId) { + const auto& g = dedxCalib_->gain().find(detId); + if (applyGain_ && g != dedxCalib_->gain().end()) + return energy * g->second; + return energy; +} + +/*****************************************************************************/ +double DeDxHitCalibrator::getChi2(const std::vector& x, + const std::vector >& b, + const double& coupling, + const double& iSigma) { + const auto& npar = b.size(); + std::vector y(npar, 0.); + for (size_t i = 0; i < npar; i++) { + const auto dx = coupling * x[i]; + if (i >= 1) + y[i - 1] += dx; + y[i] += x[i] - 2 * dx; + if (i + 1 < npar) + y[i + 1] += dx; + } + double chi2(0.); + for (size_t i = 0; i < npar; i++) { + auto q = (b[i].first - y[i]) * iSigma; + if (b[i].second == kIsNormal) + chi2 += q * q; + else if (b[i].second == kIsBelow) { + q -= 2; + if (q < 0) + chi2 += 0.5 * q * q; + } else if (b[i].second == kIsOver) { + q += 2; + if (q > 0) + chi2 += 0.5 * q * q; + } + // Penalty for negatives + if (x[i] < 0) { + q = x[i] * iSigma; + chi2 += q * q; + } + } + return chi2; +} + +/*****************************************************************************/ +void DeDxHitCalibrator::getAlphaBeta(const std::vector& x, + const std::vector >& b, + CLHEP::HepMatrix& alpha, + CLHEP::HepVector& beta, + const std::vector& isFix, + const double& coupling, + const double& iSigma) { + const auto& npar = b.size(); + const auto a0 = coupling * iSigma; + const auto a1 = (1 - 2 * coupling) * iSigma; + const auto a00 = a0 * a0; + const auto a11 = a1 * a1; + const auto a01 = a0 * a1; + std::vector y(npar, 0.); + for (size_t i = 0; i < npar; i++) + if (!isFix[i]) { + const auto dx = coupling * x[i]; + if (i >= 1) + y[i - 1] += dx; + y[i] += x[i] - 2 * dx; + if (i + 1 < npar) + y[i + 1] += dx; + } + for (size_t i = 0; i < npar; i++) { + auto q = (y[i] - b[i].first) * iSigma; + int f(0); + if (b[i].second == kIsNormal) + f = 2; + else if (b[i].second == kIsBelow) { + q += 2; + if (q > 0) + f = 1; + } else if (b[i].second == kIsOver) { + q -= 2; + if (q < 0) + f = 1; + } + if (f > 0) { + if (i >= 1) + if (!isFix[i - 1]) { + alpha[i - 1][i - 1] += f * a00; + if (!isFix[i]) { + alpha[i - 1][i] += f * a01; + alpha[i][i - 1] += f * a01; + } + beta[i - 1] += f * q * a0; + } + if (!isFix[i]) { + alpha[i][i] += f * a11; + beta[i] += f * q * a1; + } + if (i + 1 < npar) + if (!isFix[i + 1]) { + alpha[i + 1][i + 1] += f * a00; + if (!isFix[i]) { + alpha[i + 1][i] += f * a01; + alpha[i][i + 1] += f * a01; + } + beta[i + 1] += f * q * a0; + } + } + // Penalty for negatives + if (!isFix[i]) + if (x[i] < 0) { + alpha[i][i] += 2 * iSigma * iSigma; + q = x[i] * iSigma; + beta[i] += 2 * q * iSigma; + } + } + for (size_t i = 0; i < npar; i++) + if (isFix[i]) { + alpha[i][i] = 1.; + beta[i] = 0.; + } +} + +/*****************************************************************************/ +std::pair DeDxHitCalibrator::fitStripCluster(const std::vector >& b, + const double& coupling, + const double& iSigma) { + const auto& npar = b.size(); + std::vector isFix(npar, false); + std::vector x; + x.reserve(b.size()); + for (const auto& ib : b) + x.emplace_back(ib.first); + + bool ok = false; + CLHEP::HepMatrix hessian(npar, npar); + int iter = 0; + do { + double diff(0); + auto old = getChi2(x, b, coupling, iSigma); + do { + CLHEP::HepMatrix alpha(npar, npar, 0.); + CLHEP::HepVector beta(npar, 0.); + getAlphaBeta(x, b, alpha, beta, isFix, coupling, iSigma); + const auto& delta = CLHEP::solve(alpha, -beta); + double lambda(1.); + std::vector xn(npar); + for (size_t i = 0; i < npar; i++) + xn[i] = x[i] + lambda * delta[i]; + const auto& next = getChi2(xn, b, coupling, iSigma); + diff = old - next; + if (diff > 0) + x = xn; + old = next; + hessian = alpha; + iter++; + } while (diff > 1e-6 && iter < 100); + // Check if we have negatives + const auto& mi = std::min_element(x.begin(), x.end()) - x.begin(); + if (x[mi] < 0) { + x[mi] = 0.; + isFix[mi] = true; + } else + ok = true; + } while (!ok && iter < 100); + hessian /= 2.; + int flag; + hessian.invert(flag); + + double var2(0.); + for (size_t i = 0; i < npar; i++) + for (size_t j = 0; j < npar; j++) + if (!isFix[i] && !isFix[j]) + var2 += hessian[i][j]; + var2 = std::abs(var2); + + double sum(0.); + for (size_t i = 0; i < npar; i++) + if (!isFix[i]) + sum += x[i]; + + return {sum, std::sqrt(var2)}; +} + +//define this as a plug-in +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(DeDxHitCalibrator); diff --git a/RecoTracker/DeDx/plugins/DeDxHitInfoProducer.cc b/RecoTracker/DeDx/plugins/DeDxHitInfoProducer.cc index 2dee2b30b1595..593ba2242b1b8 100644 --- a/RecoTracker/DeDx/plugins/DeDxHitInfoProducer.cc +++ b/RecoTracker/DeDx/plugins/DeDxHitInfoProducer.cc @@ -36,6 +36,8 @@ #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" #include "RecoTracker/DeDx/interface/DeDxTools.h" #include "RecoTracker/DeDx/interface/GenericTruncatedAverageDeDxEstimator.h" +#include "RecoTracker/PixelLowPtUtilities/interface/ClusterShape.h" +#include "RecoTracker/PixelLowPtUtilities/interface/ClusterShapeHitFilter.h" #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h" using namespace reco; @@ -52,9 +54,10 @@ class DeDxHitInfoProducer : public edm::stream::EDProducer<> { void produce(edm::Event&, const edm::EventSetup&) override; void makeCalibrationMap(const TrackerGeometry& tkGeom_); + void processRec(reco::DeDxHitInfo&, const SiStripRecHit2D&, const LocalPoint&, const LocalVector&, const float&); void processHit(const TrackingRecHit* recHit, const float trackMomentum, - const float cosine, + const LocalVector& trackDirection, reco::DeDxHitInfo& hitDeDxInfo, const LocalPoint& hitLocalPos); @@ -79,8 +82,12 @@ class DeDxHitInfoProducer : public edm::stream::EDProducer<> { const bool usePixelForPrescales_; const edm::EDGetTokenT tracksToken_; - edm::ESHandle tkGeom_; + const edm::EDGetTokenT pixShapeCacheToken_; const edm::ESGetToken tkGeomToken_; + const edm::ESGetToken clShapeToken_; + edm::ESHandle tkGeom_; + edm::ESHandle clShape_; + edm::Handle pixShapeCache_; std::vector> calibGains_; unsigned int offsetDU_; @@ -113,7 +120,10 @@ DeDxHitInfoProducer::DeDxHitInfoProducer(const edm::ParameterSet& iConfig) lowPtTracksDeDxThreshold_(iConfig.getParameter("lowPtTracksDeDxThreshold")), usePixelForPrescales_(iConfig.getParameter("usePixelForPrescales")), tracksToken_(consumes(iConfig.getParameter("tracks"))), - tkGeomToken_(esConsumes()) { + pixShapeCacheToken_(consumes(iConfig.getParameter("clusterShapeCache"))), + tkGeomToken_(esConsumes()), + clShapeToken_( + esConsumes(edm::ESInputTag("", "ClusterShapeHitFilter"))) { produces(); produces(); produces>("prescale"); @@ -139,6 +149,9 @@ void DeDxHitInfoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe iEvent.getByToken(tracksToken_, trackCollectionHandle); const TrackCollection& trackCollection(*trackCollectionHandle.product()); + clShape_ = iSetup.getHandle(clShapeToken_); + pixShapeCache_ = iEvent.getHandle(pixShapeCacheToken_); + // creates the output collection auto resultdedxHitColl = std::make_unique(); @@ -172,10 +185,7 @@ void DeDxHitInfoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe if (!trackerHitRTTI::isFromDet(*recHit)) continue; - auto trackDirection = trajParams[h].direction(); - float cosine = trackDirection.z() / trackDirection.mag(); - - processHit(recHit, track.p(), cosine, hitDeDxInfo, trajParams[h].position()); + processHit(recHit, track.p(), trajParams[h].direction(), hitDeDxInfo, trajParams[h].position()); } if (!passPt) { @@ -236,9 +246,35 @@ void DeDxHitInfoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe iEvent.put(std::move(dedxPrescale), "prescale"); } +void DeDxHitInfoProducer::processRec(reco::DeDxHitInfo& hitDeDxInfo, + const SiStripRecHit2D& recHit, + const LocalPoint& lpos, + const LocalVector& ldir, + const float& cos) { + uint8_t type(0); + int meas; + float pred; + const auto& usable = clShape_->getSizes(recHit, {}, ldir, meas, pred); + if (usable && meas <= int(std::abs(pred)) + 4) + type |= (1 << reco::DeDxHitInfo::Complete); + if (clShape_->isCompatible(recHit, ldir)) + type |= (1 << reco::DeDxHitInfo::Compatible); + if (usable) + type |= (1 << reco::DeDxHitInfo::Calibration); + + int NSaturating(0); + const auto& detId = recHit.geographicalId(); + const auto* detUnit = recHit.detUnit(); + if (detUnit == nullptr) + detUnit = tkGeom_->idToDet(detId); + const auto pathLen = detUnit->surface().bounds().thickness() / cos; + float chargeAbs = deDxTools::getCharge(&(recHit.stripCluster()), NSaturating, *detUnit, calibGains_, offsetDU_); + hitDeDxInfo.addHit(chargeAbs, pathLen, detId, lpos, type, recHit.stripCluster()); +} + void DeDxHitInfoProducer::processHit(const TrackingRecHit* recHit, const float trackMomentum, - const float cosine, + const LocalVector& trackDirection, reco::DeDxHitInfo& hitDeDxInfo, const LocalPoint& hitLocalPos) { auto const& thit = static_cast(*recHit); @@ -246,6 +282,7 @@ void DeDxHitInfoProducer::processHit(const TrackingRecHit* recHit, return; //make sure cosine is not 0 + float cosine = trackDirection.z() / trackDirection.mag(); float cosineAbs = std::max(0.00000001f, std::abs(cosine)); auto const& clus = thit.firstClusterRef(); @@ -262,15 +299,25 @@ void DeDxHitInfoProducer::processHit(const TrackingRecHit* recHit, if (!usePixel_) return; + uint8_t type(0); + const auto& pixelDet = *dynamic_cast(detUnit); + const auto& pixelRecHit = *dynamic_cast(recHit); + ClusterData data; + ClusterShape().determineShape(pixelDet, clus.pixelCluster(), data); + if (data.isComplete) + type |= (1 << reco::DeDxHitInfo::Complete); + if (clShape_->isCompatible(pixelRecHit, trackDirection, *pixShapeCache_)) + type |= (1 << reco::DeDxHitInfo::Compatible); + if (data.isComplete && data.isStraight && data.hasBigPixelsOnlyInside) + type |= (1 << reco::DeDxHitInfo::Calibration); + float chargeAbs = clus.pixelCluster().charge(); - hitDeDxInfo.addHit(chargeAbs, pathLen, thit.geographicalId(), hitLocalPos, clus.pixelCluster()); + hitDeDxInfo.addHit(chargeAbs, pathLen, thit.geographicalId(), hitLocalPos, type, clus.pixelCluster()); } else if (clus.isStrip() && !thit.isMatched()) { if (!useStrip_) return; - int NSaturating = 0; - float chargeAbs = deDxTools::getCharge(&(clus.stripCluster()), NSaturating, *detUnit, calibGains_, offsetDU_); - hitDeDxInfo.addHit(chargeAbs, pathLen, thit.geographicalId(), hitLocalPos, clus.stripCluster()); + processRec(hitDeDxInfo, {thit.geographicalId(), clus}, hitLocalPos, trackDirection, cosineAbs); } else if (clus.isStrip() && thit.isMatched()) { if (!useStrip_) return; @@ -278,23 +325,8 @@ void DeDxHitInfoProducer::processHit(const TrackingRecHit* recHit, if (!matchedHit) return; - const auto* detUnitM = matchedHit->monoHit().detUnit(); - if (detUnitM == nullptr) - detUnitM = tkGeom_->idToDet(matchedHit->monoHit().geographicalId()); - int NSaturating = 0; - auto pathLenM = detUnitM->surface().bounds().thickness() / cosineAbs; - float chargeAbs = - deDxTools::getCharge(&(matchedHit->monoHit().stripCluster()), NSaturating, *detUnitM, calibGains_, offsetDU_); - hitDeDxInfo.addHit(chargeAbs, pathLenM, thit.geographicalId(), hitLocalPos, matchedHit->monoHit().stripCluster()); - - const auto* detUnitS = matchedHit->stereoHit().detUnit(); - if (detUnitS == nullptr) - detUnitS = tkGeom_->idToDet(matchedHit->stereoHit().geographicalId()); - NSaturating = 0; - auto pathLenS = detUnitS->surface().bounds().thickness() / cosineAbs; - chargeAbs = - deDxTools::getCharge(&(matchedHit->stereoHit().stripCluster()), NSaturating, *detUnitS, calibGains_, offsetDU_); - hitDeDxInfo.addHit(chargeAbs, pathLenS, thit.geographicalId(), hitLocalPos, matchedHit->stereoHit().stripCluster()); + processRec(hitDeDxInfo, matchedHit->monoHit(), hitLocalPos, trackDirection, cosineAbs); + processRec(hitDeDxInfo, matchedHit->stereoHit(), hitLocalPos, trackDirection, cosineAbs); } } diff --git a/RecoTracker/DeDx/python/dedxEstimators_cff.py b/RecoTracker/DeDx/python/dedxEstimators_cff.py index 5b256fe120cea..c29e32539b8f4 100644 --- a/RecoTracker/DeDx/python/dedxEstimators_cff.py +++ b/RecoTracker/DeDx/python/dedxEstimators_cff.py @@ -16,6 +16,7 @@ useCalibration = cms.bool(False), calibrationPath = cms.string("file:Gains.root"), shapeTest = cms.bool(True), + clusterShapeCache = cms.InputTag("siPixelClusterShapeCache"), lowPtTracksPrescalePass = cms.uint32(100), # prescale factor for low pt tracks above the dEdx cut lowPtTracksPrescaleFail = cms.uint32(2000), # prescale factor for low pt tracks below the dEdx cut @@ -87,3 +88,31 @@ lowPtTracksEstimatorParameters = dict(fraction = 0., exponent = -2.0,truncate = False), usePixelForPrescales = False ) + +# dEdx for Run-3 UPC +from Configuration.Eras.Modifier_run3_upc_cff import run3_upc +run3_upc.toModify(dedxHitInfo, minTrackPt = 0) + +from RecoTracker.DeDx.dedxHitCalibrator_cfi import dedxHitCalibrator as _dedxHitCalibrator +from SimGeneral.MixingModule.SiStripSimParameters_cfi import SiStripSimBlock as _SiStripSimBlock +from RecoLocalTracker.SiPixelClusterizer.SiPixelClusterizer_cfi import siPixelClusters as _siPixelClusters +dedxHitCalibrator = _dedxHitCalibrator.clone( + MeVPerElectron = 1000*_SiStripSimBlock.GevPerElectron.value(), + VCaltoElectronGain = _siPixelClusters.VCaltoElectronGain, + VCaltoElectronGain_L1 = _siPixelClusters.VCaltoElectronGain_L1, + VCaltoElectronOffset = _siPixelClusters.VCaltoElectronOffset, + VCaltoElectronOffset_L1 = _siPixelClusters.VCaltoElectronOffset_L1 +) + +dedxAllLikelihood = _mod.DeDxEstimatorProducer.clone( + UseStrip = True, UsePixel = True, + estimator = 'likelihoodFit', + UseDeDxHits = True, + pixelDeDxHits = 'dedxHitCalibrator:PixelHits', + stripDeDxHits = 'dedxHitCalibrator:StripHits' +) +dedxPixelLikelihood = dedxAllLikelihood.clone(UseStrip = False, UsePixel = True) +dedxStripLikelihood = dedxAllLikelihood.clone(UseStrip = True, UsePixel = False) + +from Configuration.Eras.Modifier_run3_egamma_2023_cff import run3_egamma_2023 +(run3_upc & ~run3_egamma_2023).toReplaceWith(doAlldEdXEstimatorsTask, cms.Task(doAlldEdXEstimatorsTask.copy(), dedxHitCalibrator, dedxStripLikelihood, dedxPixelLikelihood, dedxAllLikelihood))