diff --git a/Configuration/ProcessModifiers/python/seedingDeepCore_cff.py b/Configuration/ProcessModifiers/python/seedingDeepCore_cff.py new file mode 100644 index 0000000000000..9ec8ae8bcdcd9 --- /dev/null +++ b/Configuration/ProcessModifiers/python/seedingDeepCore_cff.py @@ -0,0 +1,5 @@ +import FWCore.ParameterSet.Config as cms + +# This modifier is for activating DeepCore seeding for the JetCore tracking iteration + +seedingDeepCore = cms.Modifier() diff --git a/Configuration/PyReleaseValidation/python/relval_2017.py b/Configuration/PyReleaseValidation/python/relval_2017.py index d7218ad988efd..c4a5da7586441 100644 --- a/Configuration/PyReleaseValidation/python/relval_2017.py +++ b/Configuration/PyReleaseValidation/python/relval_2017.py @@ -35,6 +35,7 @@ # (Patatrack ECAL-only: TTbar - on CPU, on GPU, both, auto) # (Patatrack HCAL-only: TTbar - on CPU, on GPU, both, auto) # (TTbar 0T, TTbar PU 0T) +# (QCD 1.8TeV DeepCore) # 2023 (TTbar, TTbar PU, TTbar PU premix) # 2024 (TTbar, TTbar PU, TTbar PU premix) numWFIB = [10001.0,10002.0,10003.0,10004.0,10005.0,10006.0,10007.0,10008.0,10009.0,10059.0,10071.0, @@ -59,6 +60,7 @@ 11634.511,11634.512, # 11634.513,11634.514, 11634.521,11634.522, # 11634.523,11634.524 11634.24,11834.24, + 11723.17, 12434.0,12634.0,12634.99, 12834.0,13034.0,13034.99] for numWF in numWFIB: diff --git a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py index a9876705bf631..07fc60222b8e6 100644 --- a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py +++ b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py @@ -311,6 +311,30 @@ def condition_(self, fragment, stepList, key, hasHarvest): '--procModifiers': 'trackingMkFit' } +#DeepCore seeding for JetCore iteration workflow +class UpgradeWorkflow_seedingDeepCore(UpgradeWorkflow): + def setup_(self, step, stepName, stepDict, k, properties): + if 'Reco' in step or 'HARVEST' in step: stepDict[stepName][k] = merge([{'--procModifiers': 'seedingDeepCore'}, stepDict[step][k]]) + def condition(self, fragment, stepList, key, hasHarvest): + result = (fragment=="QCD_Pt_1800_2400_14") and ('2021' in key or '2024' in key) and hasHarvest + if result: + # skip ALCA and Nano + skipList = [s for s in stepList if (("ALCA" in s) or ("Nano" in s))] + for skip in skipList: + stepList.remove(skip) + return result +upgradeWFs['seedingDeepCore'] = UpgradeWorkflow_seedingDeepCore( + steps = [ + 'Reco', + 'HARVEST', + 'RecoGlobal', + 'HARVESTGlobal', + ], + PU = [], + suffix = '_seedingDeepCore', + offset = 0.17, +) + # Vector Hits workflows class UpgradeWorkflow_vectorHits(UpgradeWorkflow): def setup_(self, step, stepName, stepDict, k, properties): diff --git a/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronSeedProducer.cc b/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronSeedProducer.cc index 5b2c463e5cc55..26d8d7d9d245f 100644 --- a/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronSeedProducer.cc +++ b/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronSeedProducer.cc @@ -340,6 +340,9 @@ void LowPtGsfElectronSeedProducer::loop(const edm::Handle >& hand // Create ElectronSeed reco::ElectronSeed seed(*(trackRef->seedRef())); + if (seed.nHits() == 0) { //if DeepCore is used in jetCore iteration the seed are hitless, in case skip + continue; + } seed.setCtfTrack(trackRef); // Create PreIds diff --git a/RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.cc b/RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.cc index febe67147f743..240a4a7a86538 100644 --- a/RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.cc +++ b/RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.cc @@ -345,6 +345,10 @@ void GoodSeedProducer::produce(Event& iEvent, const EventSetup& iSetup) { auto tketa = tkmom.eta(); auto tkpt = std::sqrt(tkmom.perp2()); auto const& Seed = (*trackRef->seedRef()); + if (Seed.nHits() == 0) { //if DeepCore is used in jetCore iteration the seed are hitless, in case skip + continue; + } + if (!disablePreId_) { int ipteta = getBin(Tk[i].eta(), Tk[i].pt()); int ibin = ipteta * 9; diff --git a/RecoParticleFlow/PFTracking/python/trackerDrivenElectronSeeds_cfi.py b/RecoParticleFlow/PFTracking/python/trackerDrivenElectronSeeds_cfi.py index c1657934f9f94..cdb6bdbf4357a 100644 --- a/RecoParticleFlow/PFTracking/python/trackerDrivenElectronSeeds_cfi.py +++ b/RecoParticleFlow/PFTracking/python/trackerDrivenElectronSeeds_cfi.py @@ -68,3 +68,4 @@ from Configuration.ProcessModifiers.egamma_lowPt_exclusive_cff import egamma_lowPt_exclusive egamma_lowPt_exclusive.toModify(trackerDrivenElectronSeeds,MinPt = 1.0) + diff --git a/RecoTracker/IterativeTracking/python/JetCoreRegionalStep_cff.py b/RecoTracker/IterativeTracking/python/JetCoreRegionalStep_cff.py index 7df01c115c676..e9f44e537d619 100644 --- a/RecoTracker/IterativeTracking/python/JetCoreRegionalStep_cff.py +++ b/RecoTracker/IterativeTracking/python/JetCoreRegionalStep_cff.py @@ -95,6 +95,14 @@ minPt = 0.1 ) +from Configuration.ProcessModifiers.seedingDeepCore_cff import seedingDeepCore +seedingDeepCore.toModify(jetCoreRegionalStepTrajectoryFilter, + minimumNumberOfHits = 2, + maxConsecLostHits = 2, + maxLostHitsFraction = 1.1, + minPt = 0.9 + ) + from Configuration.Eras.Modifier_pp_on_XeXe_2017_cff import pp_on_XeXe_2017 from Configuration.ProcessModifiers.pp_on_AA_cff import pp_on_AA (pp_on_XeXe_2017 | pp_on_AA).toModify(jetCoreRegionalStepTrajectoryFilter, minPt=5.0) @@ -118,6 +126,28 @@ estimator = 'jetCoreRegionalStepChi2Est', maxDPhiForLooperReconstruction = cms.double(2.0), maxPtForLooperReconstruction = cms.double(0.7) + ) + +seedingDeepCore.toModify(jetCoreRegionalStepTrajectoryBuilder, + maxPtForLooperReconstruction = 0., + keepOriginalIfRebuildFails = True, + lockHits = False, + requireSeedHitsInRebuild = False, +) + +#customized cleaner for DeepCore +from TrackingTools.TrajectoryCleaning.TrajectoryCleanerBySharedHits_cfi import trajectoryCleanerBySharedHits +jetCoreRegionalStepDeepCoreTrajectoryCleaner = trajectoryCleanerBySharedHits.clone( + ComponentName = 'jetCoreRegionalStepDeepCoreTrajectoryCleaner', + fractionShared = 0.45 +) + +import RecoTracker.TkSeedGenerator.deepCoreSeedGenerator_cfi +import Validation.RecoTrack.JetCoreMCtruthSeedGenerator_cfi +seedingDeepCore.toReplaceWith(jetCoreRegionalStepSeeds, + RecoTracker.TkSeedGenerator.deepCoreSeedGenerator_cfi.deepCoreSeedGenerator.clone(#to run MCtruthSeedGenerator clone here from Validation.RecoTrack + vertices="firstStepPrimaryVertices" + ) ) # MAKING OF TRACK CANDIDATES @@ -131,6 +161,10 @@ #numHitsForSeedCleaner = cms.int32(50), #onlyPixelHitsForSeedCleaner = cms.bool(True), ) +seedingDeepCore.toModify(jetCoreRegionalStepTrackCandidates, + TrajectoryCleaner = 'jetCoreRegionalStepDeepCoreTrajectoryCleaner', + doSeedingRegionRebuilding = True, +) # TRACK FITTING @@ -219,6 +253,7 @@ # jetCoreRegionalStepClassifier1,jetCoreRegionalStepClassifier2, jetCoreRegionalStep) JetCoreRegionalStep = cms.Sequence(JetCoreRegionalStepTask) +seedingDeepCore.toReplaceWith(JetCoreRegionalStep,JetCoreRegionalStep.copyAndExclude([jetCoreRegionalStepHitDoublets])) fastSim.toReplaceWith(JetCoreRegionalStepTask, cms.Task(jetCoreRegionalStepTracks, jetCoreRegionalStep)) diff --git a/RecoTracker/TkSeedGenerator/plugins/BuildFile.xml b/RecoTracker/TkSeedGenerator/plugins/BuildFile.xml index abaa5754e68c7..62511823f5ae7 100644 --- a/RecoTracker/TkSeedGenerator/plugins/BuildFile.xml +++ b/RecoTracker/TkSeedGenerator/plugins/BuildFile.xml @@ -1,3 +1,4 @@ + diff --git a/RecoTracker/TkSeedGenerator/plugins/DeepCoreSeedGenerator.cc b/RecoTracker/TkSeedGenerator/plugins/DeepCoreSeedGenerator.cc new file mode 100644 index 0000000000000..6ffb2349d5c00 --- /dev/null +++ b/RecoTracker/TkSeedGenerator/plugins/DeepCoreSeedGenerator.cc @@ -0,0 +1,577 @@ +// -*- C++ -*- +// +// Package: trackJet/DeepCoreSeedGenerator +// Class: DeepCoreSeedGenerator +// +/**\class DeepCoreSeedGenerator DeepCoreSeedGenerator.cc trackJet/DeepCoreSeedGenerator/plugins/DeepCoreSeedGenerator.cc + + Description: [one line class summary] + + Implementation: + [Notes on implementation] +*/ +// +// Original Author: Valerio Bertacchi +// Created: Mon, 18 Dec 2017 16:35:04 GMT +// +// + +// system include files + +#include + +// user include files + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h" +#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h" +#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h" +#include "Geometry/Records/interface/TrackerTopologyRcd.h" + +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/Common/interface/DetSet.h" +#include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/JetReco/interface/Jet.h" +#include "DataFormats/SiPixelDigi/interface/PixelDigi.h" +#include "DataFormats/GeometryVector/interface/VectorUtil.h" +#include "DataFormats/Math/interface/Point3D.h" +#include "DataFormats/Math/interface/Vector3D.h" +#include "DataFormats/Candidate/interface/Candidate.h" + +#include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h" +#include "RecoLocalTracker/Records/interface/TkPixelCPERecord.h" + +#include "TrackingTools/GeomPropagators/interface/StraightLinePlaneCrossing.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" + +#include "MagneticField/Engine/interface/MagneticField.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" + +#include "FWCore/ServiceRegistry/interface/Service.h" + +#include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h" + +#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h" + +#include "PhysicsTools/TensorFlow/interface/TensorFlow.h" + +// +// class declaration +// +struct DeepCoreCache { + const tensorflow::GraphDef* graph_def; +}; + +class DeepCoreSeedGenerator : public edm::stream::EDProducer> { +public: + explicit DeepCoreSeedGenerator(const edm::ParameterSet&, const DeepCoreCache*); + ~DeepCoreSeedGenerator() override; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + // A pointer to a cluster and a list of tracks on it + + // static methods for handling the global cache + static std::unique_ptr initializeGlobalCache(const edm::ParameterSet&); + static void globalEndJob(DeepCoreCache*); + + double jetPt_; + double jetEta_; + double pitchX_ = 0.01; //100 um (pixel pitch in X) + double pitchY_ = 0.015; //150 um (pixel pitch in Y) + static constexpr int jetDimX = 30; //pixel dimension of NN window on layer2 + static constexpr int jetDimY = 30; //pixel dimension of NN window on layer2 + static constexpr int Nlayer = 4; //Number of layer used in DeepCore + static constexpr int Nover = 3; //Max number of tracks recorded per pixel + static constexpr int Npar = 5; //Number of track parameter + +private: + void beginJob(); + void produce(edm::Event&, const edm::EventSetup&) override; + void endJob(); + + // ----------member data --------------------------- + std::string propagatorName_; + edm::ESHandle magfield_; + edm::ESHandle geometry_; + edm::ESHandle propagator_; + + edm::EDGetTokenT> vertices_; + edm::EDGetTokenT> pixelClusters_; + edm::Handle> inputPixelClusters_; + edm::EDGetTokenT> cores_; + const edm::ESGetToken topoToken_; + + double ptMin_; + double deltaR_; + double chargeFracMin_; + double centralMIPCharge_; + std::string pixelCPE_; + std::string weightfilename_; + std::vector inputTensorName_; + std::vector outputTensorName_; + double probThr_; + tensorflow::Session* session_; + + std::pair> findIntersection(const GlobalVector&, + const reco::Candidate::Point&, + const GeomDet*); + + void fillPixelMatrix(const SiPixelCluster&, + int, + Point3DBase, + const GeomDet*, + tensorflow::NamedTensorList); //if not working,: args=2 auto + + std::pair local2Pixel(double, double, const GeomDet*); + + LocalPoint pixel2Local(int, int, const GeomDet*); + + int pixelFlipper(const GeomDet*); + + const GeomDet* DetectorSelector(int, + const reco::Candidate&, + GlobalVector, + const reco::Vertex&, + const TrackerTopology* const, + const edmNew::DetSetVector&); + + std::vector splittedClusterDirections( + const reco::Candidate&, + const TrackerTopology* const, + const PixelClusterParameterEstimator*, + const reco::Vertex&, + int, + const edmNew::DetSetVector&); //if not working,: args=2 auto + + std::pair SeedEvaluation( + tensorflow::NamedTensorList, std::vector); +}; + +DeepCoreSeedGenerator::DeepCoreSeedGenerator(const edm::ParameterSet& iConfig, const DeepCoreCache* cache) + : vertices_(consumes(iConfig.getParameter("vertices"))), + pixelClusters_( + consumes>(iConfig.getParameter("pixelClusters"))), + cores_(consumes>(iConfig.getParameter("cores"))), + topoToken_(esConsumes()), + ptMin_(iConfig.getParameter("ptMin")), + deltaR_(iConfig.getParameter("deltaR")), + chargeFracMin_(iConfig.getParameter("chargeFractionMin")), + centralMIPCharge_(iConfig.getParameter("centralMIPCharge")), + pixelCPE_(iConfig.getParameter("pixelCPE")), + weightfilename_(iConfig.getParameter("weightFile").fullPath()), + inputTensorName_(iConfig.getParameter>("inputTensorName")), + outputTensorName_(iConfig.getParameter>("outputTensorName")), + probThr_(iConfig.getParameter("probThr")), + session_(tensorflow::createSession(cache->graph_def)) + +{ + produces(); + produces(); +} + +DeepCoreSeedGenerator::~DeepCoreSeedGenerator() {} + +void DeepCoreSeedGenerator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + auto result = std::make_unique(); + auto resultTracks = std::make_unique(); + + const tensorflow::TensorShape input_size_eta({1, 1}); + const tensorflow::TensorShape input_size_pt({1, 1}); + const tensorflow::TensorShape input_size_cluster({1, jetDimX, jetDimY, Nlayer}); + std::vector output_names; + output_names.push_back(outputTensorName_[0]); + output_names.push_back(outputTensorName_[1]); + + using namespace edm; + using namespace reco; + + iSetup.get().get(magfield_); + iSetup.get().get(geometry_); + iSetup.get().get("AnalyticalPropagator", propagator_); + + const auto& inputPixelClusters_ = iEvent.get(pixelClusters_); + const auto& vertices = iEvent.get(vertices_); + const auto& cores = iEvent.get(cores_); + + edm::ESHandle pixelCPEhandle; + const PixelClusterParameterEstimator* pixelCPE; + iSetup.get().get(pixelCPE_, pixelCPEhandle); + pixelCPE = pixelCPEhandle.product(); + + const TrackerTopology* const tTopo = &iSetup.getData(topoToken_); + auto output = std::make_unique>(); + + for (const auto& jet : cores) { // jet loop + + if (jet.pt() > ptMin_) { + std::set ids; + const reco::Vertex& jetVertex = vertices[0]; + + std::vector splitClustDirSet = + splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 1, inputPixelClusters_); + bool l2off = (splitClustDirSet.empty()); + if (splitClustDirSet.empty()) { //if layer 1 is broken find direcitons on layer 2 + splitClustDirSet = splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 2, inputPixelClusters_); + } + splitClustDirSet.emplace_back(GlobalVector(jet.px(), jet.py(), jet.pz())); + for (int cc = 0; cc < (int)splitClustDirSet.size(); cc++) { + tensorflow::NamedTensorList input_tensors; //TensorFlow tensors init + input_tensors.resize(3); + input_tensors[0] = + tensorflow::NamedTensor(inputTensorName_[0], tensorflow::Tensor(tensorflow::DT_FLOAT, input_size_eta)); + input_tensors[1] = + tensorflow::NamedTensor(inputTensorName_[1], tensorflow::Tensor(tensorflow::DT_FLOAT, input_size_pt)); + input_tensors[2] = tensorflow::NamedTensor(inputTensorName_[2], + tensorflow::Tensor(tensorflow::DT_FLOAT, {input_size_cluster})); + + //put all the input tensor to 0 + input_tensors[0].second.matrix()(0, 0) = 0.0; + input_tensors[1].second.matrix()(0, 0) = 0.0; + for (int x = 0; x < jetDimX; x++) { + for (int y = 0; y < jetDimY; y++) { + for (int l = 0; l < 4; l++) { + input_tensors[2].second.tensor()(0, x, y, l) = 0.0; + } + } + } //end of TensorFlow tensors init + + GlobalVector bigClustDir = splitClustDirSet[cc]; + + jetEta_ = jet.eta(); + jetPt_ = jet.pt(); + input_tensors[0].second.matrix()(0, 0) = jet.eta(); + input_tensors[1].second.matrix()(0, 0) = jet.pt(); + + const GeomDet* globDet = DetectorSelector( + 2, jet, bigClustDir, jetVertex, tTopo, inputPixelClusters_); //det. aligned to bigClustDir on L2 + + if (globDet == nullptr) //no intersection between bigClustDir and pixel detector modules found + continue; + + const GeomDet* goodDet1 = DetectorSelector(1, jet, bigClustDir, jetVertex, tTopo, inputPixelClusters_); + const GeomDet* goodDet3 = DetectorSelector(3, jet, bigClustDir, jetVertex, tTopo, inputPixelClusters_); + const GeomDet* goodDet4 = DetectorSelector(4, jet, bigClustDir, jetVertex, tTopo, inputPixelClusters_); + + for (const auto& detset : inputPixelClusters_) { + const GeomDet* det = geometry_->idToDet(detset.id()); + + for (const auto& aCluster : detset) { + det_id_type aClusterID = detset.id(); + if (DetId(aClusterID).subdetId() != 1) + continue; + + int lay = tTopo->layer(det->geographicalId()); + + std::pair> interPair = + findIntersection(bigClustDir, (reco::Candidate::Point)jetVertex.position(), det); + if (interPair.first == false) + continue; + Basic3DVector inter = interPair.second; + auto localInter = det->specificSurface().toLocal((GlobalPoint)inter); + + GlobalPoint pointVertex(jetVertex.position().x(), jetVertex.position().y(), jetVertex.position().z()); + + LocalPoint clustPos_local = + pixelCPE->localParametersV(aCluster, (*geometry_->idToDetUnit(detset.id())))[0].first; + + if (std::abs(clustPos_local.x() - localInter.x()) / pitchX_ <= jetDimX / 2 && + std::abs(clustPos_local.y() - localInter.y()) / pitchY_ <= + jetDimY / 2) { //used the baricenter, better description maybe useful + + if (det == goodDet1 || det == goodDet3 || det == goodDet4 || det == globDet) { + fillPixelMatrix(aCluster, lay, localInter, det, input_tensors); + } + } //cluster in ROI + } //cluster + } //detset + + //here the NN produce the seed from the filled input + std::pair seedParamNN = + DeepCoreSeedGenerator::SeedEvaluation(input_tensors, output_names); + + for (int i = 0; i < jetDimX; i++) { + for (int j = 0; j < jetDimY; j++) { + for (int o = 0; o < Nover; o++) { + if (seedParamNN.second[i][j][o] > (probThr_ - o * 0.1 - (l2off ? 0.35 : 0))) { + std::pair> interPair = + findIntersection(bigClustDir, (reco::Candidate::Point)jetVertex.position(), globDet); + auto localInter = globDet->specificSurface().toLocal((GlobalPoint)interPair.second); + + int flip = pixelFlipper(globDet); // 1=not flip, -1=flip + int nx = i - jetDimX / 2; + int ny = j - jetDimY / 2; + nx = flip * nx; + std::pair pixInter = local2Pixel(localInter.x(), localInter.y(), globDet); + nx = nx + pixInter.first; + ny = ny + pixInter.second; + LocalPoint xyLocal = pixel2Local(nx, ny, globDet); + + double xx = xyLocal.x() + seedParamNN.first[i][j][o][0] * 0.01; //0.01=internal normalization of NN + double yy = xyLocal.y() + seedParamNN.first[i][j][o][1] * 0.01; //0.01=internal normalization of NN + LocalPoint localSeedPoint = LocalPoint(xx, yy, 0); + + double track_eta = + seedParamNN.first[i][j][o][2] * 0.01 + bigClustDir.eta(); //pay attention to this 0.01 + double track_theta = 2 * std::atan(std::exp(-track_eta)); + double track_phi = + seedParamNN.first[i][j][o][3] * 0.01 + bigClustDir.phi(); //pay attention to this 0.01 + + double pt = 1. / seedParamNN.first[i][j][o][4]; + double normdirR = pt / sin(track_theta); + + const GlobalVector globSeedDir( + GlobalVector::Polar(Geom::Theta(track_theta), Geom::Phi(track_phi), normdirR)); + LocalVector localSeedDir = globDet->surface().toLocal(globSeedDir); + int64_t seedid = (int64_t(xx * 200.) << 0) + (int64_t(yy * 200.) << 16) + + (int64_t(track_eta * 400.) << 32) + (int64_t(track_phi * 400.) << 48); + if (ids.count(seedid) != 0) { + continue; + } + //to disable jetcore iteration comment from here to probability block end---> + //seeing iteration skipped, useful to total eff and FakeRate comparison + ids.insert(seedid); + + //Covariance matrix, the hadrcoded variances = NN residuals width + //(see https://twiki.cern.ch/twiki/bin/view/CMSPublic/NNJetCoreAtCtD2019) + float em[15] = {0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0}; // (see LocalTrajectoryError for details), order as follow: + em[0] = 0.15 * 0.15; // q/pt + em[2] = 0.5e-5; // dxdz + em[5] = 0.5e-5; // dydz + em[9] = 2e-5; // x + em[14] = 2e-5; // y + long int detId = globDet->geographicalId(); + LocalTrajectoryParameters localParam(localSeedPoint, localSeedDir, TrackCharge(1)); + result->emplace_back(TrajectorySeed(PTrajectoryStateOnDet(localParam, pt, em, detId, /*surfaceSide*/ 0), + edm::OwnVector(), + PropagationDirection::alongMomentum)); + + GlobalPoint globalSeedPoint = globDet->surface().toGlobal(localSeedPoint); + reco::Track::CovarianceMatrix mm; + resultTracks->emplace_back( + reco::Track(1, + 1, + reco::Track::Point(globalSeedPoint.x(), globalSeedPoint.y(), globalSeedPoint.z()), + reco::Track::Vector(globSeedDir.x(), globSeedDir.y(), globSeedDir.z()), + 1, + mm)); + } //probability + } + } + } + } //bigcluster + } //jet > pt + } //jet + iEvent.put(std::move(result)); + iEvent.put(std::move(resultTracks)); +} + +std::pair> DeepCoreSeedGenerator::findIntersection(const GlobalVector& dir, + const reco::Candidate::Point& vertex, + const GeomDet* det) { + StraightLinePlaneCrossing vertexPlane(Basic3DVector(vertex.x(), vertex.y(), vertex.z()), + Basic3DVector(dir.x(), dir.y(), dir.z())); + + std::pair> pos = vertexPlane.position(det->specificSurface()); + + return pos; +} + +std::pair DeepCoreSeedGenerator::local2Pixel(double locX, double locY, const GeomDet* det) { + LocalPoint locXY(locX, locY); + float pixX = (dynamic_cast(det))->specificTopology().pixel(locXY).first; + float pixY = (dynamic_cast(det))->specificTopology().pixel(locXY).second; + std::pair out(pixX, pixY); + return out; +} + +LocalPoint DeepCoreSeedGenerator::pixel2Local(int pixX, int pixY, const GeomDet* det) { + float locX = (dynamic_cast(det))->specificTopology().localX(pixX); + float locY = (dynamic_cast(det))->specificTopology().localY(pixY); + LocalPoint locXY(locX, locY); + return locXY; +} + +int DeepCoreSeedGenerator::pixelFlipper(const GeomDet* det) { + int out = 1; + LocalVector locZdir(0, 0, 1); + GlobalVector globZdir = det->specificSurface().toGlobal(locZdir); + const GlobalPoint& globDetCenter = det->position(); + float direction = + globZdir.x() * globDetCenter.x() + globZdir.y() * globDetCenter.y() + globZdir.z() * globDetCenter.z(); + if (direction < 0) + out = -1; + return out; +} + +void DeepCoreSeedGenerator::fillPixelMatrix(const SiPixelCluster& cluster, + int layer, + Point3DBase inter, + const GeomDet* det, + tensorflow::NamedTensorList input_tensors) { + int flip = pixelFlipper(det); // 1=not flip, -1=flip + + for (int i = 0; i < cluster.size(); i++) { + SiPixelCluster::Pixel pix = cluster.pixel(i); + std::pair pixInter = local2Pixel(inter.x(), inter.y(), det); + int nx = pix.x - pixInter.first; + int ny = pix.y - pixInter.second; + nx = flip * nx; + + if (abs(nx) < jetDimX / 2 && abs(ny) < jetDimY / 2) { + nx = nx + jetDimX / 2; + ny = ny + jetDimY / 2; + //14000 = normalization of ACD counts used in the NN + input_tensors[2].second.tensor()(0, nx, ny, layer - 1) += (pix.adc) / (14000.f); + } + } +} + +std::pair +DeepCoreSeedGenerator::SeedEvaluation(tensorflow::NamedTensorList input_tensors, + std::vector output_names) { + std::vector outputs; + tensorflow::run(session_, input_tensors, output_names, &outputs); + auto matrix_output_par = outputs.at(0).tensor(); + auto matrix_output_prob = outputs.at(1).tensor(); + + std::pair output_combined; + + for (int x = 0; x < jetDimX; x++) { + for (int y = 0; y < jetDimY; y++) { + for (int trk = 0; trk < Nover; trk++) { + output_combined.second[x][y][trk] = + matrix_output_prob(0, x, y, trk, 0); //outputs.at(1).matrix()(0,x,y,trk); + + for (int p = 0; p < Npar; p++) { + output_combined.first[x][y][trk][p] = + matrix_output_par(0, x, y, trk, p); //outputs.at(0).matrix()(0,x,y,trk,p); + } + } + } + } + return output_combined; +} + +const GeomDet* DeepCoreSeedGenerator::DetectorSelector(int llay, + const reco::Candidate& jet, + GlobalVector jetDir, + const reco::Vertex& jetVertex, + const TrackerTopology* const tTopo, + const edmNew::DetSetVector& clusters) { + double minDist = 0.0; + GeomDet* output = (GeomDet*)nullptr; + for (const auto& detset : clusters) { + auto aClusterID = detset.id(); + if (DetId(aClusterID).subdetId() != 1) + continue; + const GeomDet* det = geometry_->idToDet(aClusterID); + int lay = tTopo->layer(det->geographicalId()); + if (lay != llay) + continue; + std::pair> interPair = + findIntersection(jetDir, (reco::Candidate::Point)jetVertex.position(), det); + if (interPair.first == false) + continue; + Basic3DVector inter = interPair.second; + auto localInter = det->specificSurface().toLocal((GlobalPoint)inter); + if ((minDist == 0.0 || std::abs(localInter.x()) < minDist) && std::abs(localInter.y()) < 3.35) { + minDist = std::abs(localInter.x()); + output = (GeomDet*)det; + } + } //detset + return output; +} + +std::vector DeepCoreSeedGenerator::splittedClusterDirections( + const reco::Candidate& jet, + const TrackerTopology* const tTopo, + const PixelClusterParameterEstimator* pixelCPE, + const reco::Vertex& jetVertex, + int layer, + const edmNew::DetSetVector& clusters) { + std::vector clustDirs; + for (const auto& detset_int : clusters) { + const GeomDet* det_int = geometry_->idToDet(detset_int.id()); + int lay = tTopo->layer(det_int->geographicalId()); + if (lay != layer) + continue; //NB: saved bigClusters on all the layers!! + auto detUnit = geometry_->idToDetUnit(detset_int.id()); + for (const auto& aCluster : detset_int) { + GlobalPoint clustPos = det_int->surface().toGlobal(pixelCPE->localParametersV(aCluster, (*detUnit))[0].first); + GlobalPoint vertexPos(jetVertex.position().x(), jetVertex.position().y(), jetVertex.position().z()); + GlobalVector clusterDir = clustPos - vertexPos; + GlobalVector jetDir(jet.px(), jet.py(), jet.pz()); + if (Geom::deltaR(jetDir, clusterDir) < deltaR_) { + clustDirs.emplace_back(clusterDir); + } + } + } + return clustDirs; +} + +std::unique_ptr DeepCoreSeedGenerator::initializeGlobalCache(const edm::ParameterSet& iConfig) { + // this method is supposed to create, initialize and return a DeepCoreCache instance + std::unique_ptr cache = std::make_unique(); + std::string graphPath = iConfig.getParameter("weightFile").fullPath(); + cache->graph_def = tensorflow::loadGraphDef(graphPath); + return cache; +} + +void DeepCoreSeedGenerator::globalEndJob(DeepCoreCache* cache) { delete cache->graph_def; } + +// ------------ method called once each job just before starting event loop ------------ +void DeepCoreSeedGenerator::beginJob() {} + +// ------------ method called once each job just after ending the event loop ------------ +void DeepCoreSeedGenerator::endJob() {} + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +void DeepCoreSeedGenerator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("vertices", edm::InputTag("offlinePrimaryVertices")); + desc.add("pixelClusters", edm::InputTag("siPixelClustersPreSplitting")); + desc.add("cores", edm::InputTag("jetsForCoreTracking")); + desc.add("ptMin", 300); + desc.add("deltaR", 0.1); + desc.add("chargeFractionMin", 18000.0); + desc.add("centralMIPCharge", 2); + desc.add("pixelCPE", "PixelCPEGeneric"); + desc.add( + "weightFile", + edm::FileInPath("RecoTracker/TkSeedGenerator/data/DeepCore/DeepCoreSeedGenerator_TrainedModel_barrel_2017.pb")); + desc.add>("inputTensorName", {"input_1", "input_2", "input_3"}); + desc.add>("outputTensorName", {"output_node0", "output_node1"}); + desc.add("probThr", 0.85); + descriptions.add("deepCoreSeedGenerator", desc); +} + +DEFINE_FWK_MODULE(DeepCoreSeedGenerator); diff --git a/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionImpl.cc b/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionImpl.cc index 0773dc621fbe0..793a75b82aa91 100644 --- a/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionImpl.cc +++ b/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionImpl.cc @@ -116,7 +116,7 @@ double TrackAssociatorByPositionImpl::quality(const TrajectoryStateOnSurface& tr RecoToSimCollection TrackAssociatorByPositionImpl::associateRecoToSim( const edm::RefToBaseVector& tCH, const edm::RefVector& tPCH) const { - RecoToSimCollection outputCollection; + RecoToSimCollection outputCollection(productGetter_); //for each reco track find a matching tracking particle std::pair minPair; const double dQmin_default = 1542543; @@ -168,7 +168,7 @@ RecoToSimCollection TrackAssociatorByPositionImpl::associateRecoToSim( SimToRecoCollection TrackAssociatorByPositionImpl::associateSimToReco( const edm::RefToBaseVector& tCH, const edm::RefVector& tPCH) const { - SimToRecoCollection outputCollection; + SimToRecoCollection outputCollection(productGetter_); //for each tracking particle, find matching tracks. std::pair minPair; diff --git a/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionImpl.h b/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionImpl.h index 252ed66635e74..d03f57c8e34db 100644 --- a/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionImpl.h +++ b/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionImpl.h @@ -21,6 +21,7 @@ #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" +#include "DataFormats/Common/interface/EDProductGetter.h" #include @@ -33,7 +34,8 @@ class TrackAssociatorByPositionImpl : public reco::TrackToTrackingParticleAssoci typedef std::vector SimHitTPAssociationList; enum class Method { chi2, dist, momdr, posdr }; - TrackAssociatorByPositionImpl(const TrackingGeometry* geo, + TrackAssociatorByPositionImpl(edm::EDProductGetter const& productGetter, + const TrackingGeometry* geo, const Propagator* prop, const SimHitTPAssociationList* assocList, double qMinCut, @@ -42,7 +44,8 @@ class TrackAssociatorByPositionImpl : public reco::TrackToTrackingParticleAssoci Method method, bool minIfNoMatch, bool considerAllSimHits) - : theGeometry(geo), + : productGetter_(&productGetter), + theGeometry(geo), thePropagator(prop), theSimHitsTPAssoc(assocList), theQminCut(qMinCut), @@ -65,6 +68,7 @@ class TrackAssociatorByPositionImpl : public reco::TrackToTrackingParticleAssoci private: double quality(const TrajectoryStateOnSurface&, const TrajectoryStateOnSurface&) const; + edm::EDProductGetter const* productGetter_; const TrackingGeometry* theGeometry; const Propagator* thePropagator; const SimHitTPAssociationList* theSimHitsTPAssoc; diff --git a/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionProducer.cc b/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionProducer.cc index 4b4c116c306af..24748ee050e5a 100644 --- a/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionProducer.cc +++ b/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionProducer.cc @@ -121,7 +121,8 @@ void TrackAssociatorByPositionProducer::produce(edm::StreamID, iSetup.get().get(theG); std::unique_ptr impl{ - new TrackAssociatorByPositionImpl(theG.product(), + new TrackAssociatorByPositionImpl(iEvent.productGetter(), + theG.product(), theP.product(), assocList.product(), theQminCut, diff --git a/TrackingTools/GsfTracking/python/CkfElectronCandidateMaker_cff.py b/TrackingTools/GsfTracking/python/CkfElectronCandidateMaker_cff.py index 1f47cce012bb5..69f0587d042ff 100644 --- a/TrackingTools/GsfTracking/python/CkfElectronCandidateMaker_cff.py +++ b/TrackingTools/GsfTracking/python/CkfElectronCandidateMaker_cff.py @@ -41,7 +41,8 @@ updator = 'KFUpdator' ) - +from Configuration.ProcessModifiers.seedingDeepCore_cff import seedingDeepCore +seedingDeepCore.toModify(TrajectoryBuilderForElectrons, maxPtForLooperReconstruction = cms.double(0.0) ) # CKFTrackCandidateMaker from RecoTracker.CkfPattern.CkfTrackCandidates_cff import * diff --git a/Validation/RecoTrack/interface/MTVHistoProducerAlgoForTracker.h b/Validation/RecoTrack/interface/MTVHistoProducerAlgoForTracker.h index 9fb9096a82f4a..99aace914f292 100644 --- a/Validation/RecoTrack/interface/MTVHistoProducerAlgoForTracker.h +++ b/Validation/RecoTrack/interface/MTVHistoProducerAlgoForTracker.h @@ -74,10 +74,6 @@ struct MTVHistoProducerAlgoForTrackerHistograms { std::vector h_reco_dzpvsigcut, h_assoc_dzpvsigcut, h_assoc2_dzpvsigcut, h_simul_dzpvsigcut, h_simul2_dzpvsigcut, h_pileup_dzpvsigcut; - std::vector h_reco_dzpvcut_pt, h_assoc_dzpvcut_pt, h_assoc2_dzpvcut_pt, h_simul_dzpvcut_pt, - h_simul2_dzpvcut_pt, h_pileup_dzpvcut_pt; - std::vector h_reco_dzpvsigcut_pt, h_assoc_dzpvsigcut_pt, h_assoc2_dzpvsigcut_pt, h_simul_dzpvsigcut_pt, - h_simul2_dzpvsigcut_pt, h_pileup_dzpvsigcut_pt; std::vector h_reco_simpvz, h_assoc_simpvz, h_assoc2_simpvz, h_simul_simpvz, h_looper_simpvz, h_pileup_simpvz; std::vector h_reco_seedingLayerSet, h_assoc2_seedingLayerSet, h_looper_seedingLayerSet, diff --git a/Validation/RecoTrack/plugins/JetCoreMCtruthSeedGenerator.cc b/Validation/RecoTrack/plugins/JetCoreMCtruthSeedGenerator.cc new file mode 100644 index 0000000000000..a0f366c42b493 --- /dev/null +++ b/Validation/RecoTrack/plugins/JetCoreMCtruthSeedGenerator.cc @@ -0,0 +1,528 @@ +// -*- C++ -*- +// +// Package: trackJet/JetCoreMCtruthSeedGenerator +// Class: JetCoreMCtruthSeedGenerator +// +/**\class JetCoreMCtruthSeedGenerator JetCoreMCtruthSeedGenerator.cc trackJet/JetCoreMCtruthSeedGenerator/plugins/JetCoreMCtruthSeedGenerator.cc + Description: [one line class summary] + Implementation: + [Notes on implementation] +*/ +// +// Original Author: Valerio Bertacchi +// Created: Mon, 18 Dec 2017 16:35:04 GMT +// +// + +// system include files + +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/one/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h" +#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h" +#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h" +#include "Geometry/Records/interface/TrackerTopologyRcd.h" + +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/Common/interface/DetSet.h" +#include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/JetReco/interface/Jet.h" +#include "DataFormats/SiPixelDigi/interface/PixelDigi.h" +#include "DataFormats/GeometryVector/interface/VectorUtil.h" +#include "DataFormats/Math/interface/Point3D.h" +#include "DataFormats/Math/interface/Vector3D.h" +#include "DataFormats/Candidate/interface/Candidate.h" +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" + +#include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h" +#include "RecoLocalTracker/Records/interface/TkPixelCPERecord.h" + +#include "TrackingTools/GeomPropagators/interface/StraightLinePlaneCrossing.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" + +#include "MagneticField/Engine/interface/MagneticField.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" + +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" + +#include "SimDataFormats/Track/interface/SimTrack.h" +#include "SimDataFormats/Vertex/interface/SimVertex.h" + +#include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h" + +#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h" +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +// +// class declaration +// + +class JetCoreMCtruthSeedGenerator : public edm::one::EDProducer { +public: + explicit JetCoreMCtruthSeedGenerator(const edm::ParameterSet&); + ~JetCoreMCtruthSeedGenerator() override; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + double jetPt_; + double jetEta_; + double pitchX_ = 0.01; //100 um (pixel pitch in X) + double pitchY_ = 0.015; //150 um (pixel pitch in Y) + static constexpr int jetDimX = 30; //pixel dimension of NN window on layer2 + static constexpr int jetDimY = 30; //pixel dimension of NN window on layer2 + bool inclusiveConeSeed_ = + true; //true= fill tracks in a cone of deltaR_, false=fill tracks which have SimHit on globDet + +private: + void beginJob() override; + void produce(edm::Event&, const edm::EventSetup&) override; + void endJob() override; + + // ----------member data --------------------------- + std::string propagatorName_; + edm::ESHandle magfield_; + edm::ESHandle geometry_; + edm::ESHandle propagator_; + + edm::EDGetTokenT> vertices_; + edm::EDGetTokenT> pixelClusters_; + edm::Handle> inputPixelClusters_; + edm::EDGetTokenT> cores_; + edm::EDGetTokenT> simtracksToken_; + edm::EDGetTokenT> simvertexToken_; + edm::EDGetTokenT> PSimHitToken_; + edm::Handle> simhits_; + + double ptMin_; + double deltaR_; + double chargeFracMin_; + double centralMIPCharge_; + std::string pixelCPE_; + + std::pair> findIntersection(const GlobalVector&, + const reco::Candidate::Point&, + const GeomDet*); + + const GeomDet* DetectorSelector(int, + const reco::Candidate&, + GlobalVector, + const reco::Vertex&, + const TrackerTopology* const, + const edmNew::DetSetVector&); + + std::vector splittedClusterDirections( + const reco::Candidate&, + const TrackerTopology* const, + const PixelClusterParameterEstimator*, + const reco::Vertex&, + int, + const edmNew::DetSetVector&); //if not working,: args=2 auto + + std::vector coreHitsFilling(std::vector, + const GeomDet*, + GlobalVector, + const reco::Vertex&); //if not working,: args=0 auto + + std::pair, std::vector> coreTracksFilling( + std::vector, + const std::vector, + const std::vector); //if not working,: args=1,2 auto + + std::vector> seedParFilling(std::pair, std::vector>, + const GeomDet*, + const reco::Candidate&); + + std::pair, std::vector> coreTracksFillingDeltaR( + const std::vector, + const std::vector, + const GeomDet*, + const reco::Candidate&, + const reco::Vertex&); //if not working,: args=0,1 auto +}; + +JetCoreMCtruthSeedGenerator::JetCoreMCtruthSeedGenerator(const edm::ParameterSet& iConfig) + : + + vertices_(consumes(iConfig.getParameter("vertices"))), + pixelClusters_( + consumes>(iConfig.getParameter("pixelClusters"))), + cores_(consumes>(iConfig.getParameter("cores"))), + simtracksToken_(consumes>(iConfig.getParameter("simTracks"))), + simvertexToken_(consumes>(iConfig.getParameter("simVertex"))), + PSimHitToken_(consumes>(iConfig.getParameter("simHit"))), + ptMin_(iConfig.getParameter("ptMin")), + deltaR_(iConfig.getParameter("deltaR")), + chargeFracMin_(iConfig.getParameter("chargeFractionMin")), + centralMIPCharge_(iConfig.getParameter("centralMIPCharge")), + pixelCPE_(iConfig.getParameter("pixelCPE")) + +{ + produces(); + produces(); +} + +JetCoreMCtruthSeedGenerator::~JetCoreMCtruthSeedGenerator() {} + +void JetCoreMCtruthSeedGenerator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + auto result = std::make_unique(); + auto resultTracks = std::make_unique(); + + using namespace edm; + using namespace reco; + + iSetup.get().get(magfield_); + iSetup.get().get(geometry_); + iSetup.get().get("AnalyticalPropagator", propagator_); + + const auto& inputPixelClusters_ = iEvent.get(pixelClusters_); + const auto& simtracksVector = iEvent.get(simtracksToken_); + const auto& simvertexVector = iEvent.get(simvertexToken_); + const auto& simhits_ = iEvent.get(PSimHitToken_); + const auto& vertices = iEvent.get(vertices_); + const auto& cores = iEvent.get(cores_); + + edm::ESHandle pixelCPEhandle; + const PixelClusterParameterEstimator* pixelCPE; + iSetup.get().get(pixelCPE_, pixelCPEhandle); + pixelCPE = pixelCPEhandle.product(); + + edm::ESHandle tTopoHandle; + iSetup.get().get(tTopoHandle); + const TrackerTopology* const tTopo = tTopoHandle.product(); + + auto output = std::make_unique>(); + + for (const auto& jet : cores) { //jet loop + + if (jet.pt() > ptMin_) { + std::set ids; + const reco::Vertex& jetVertex = vertices[0]; + + std::vector splitClustDirSet = + splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 1, inputPixelClusters_); + if (splitClustDirSet.empty()) { //if layer 1 is broken find direcitons on layer 2 + splitClustDirSet = splittedClusterDirections(jet, tTopo, pixelCPE, jetVertex, 2, inputPixelClusters_); + } + if (inclusiveConeSeed_) + splitClustDirSet.clear(); + splitClustDirSet.emplace_back(GlobalVector(jet.px(), jet.py(), jet.pz())); + + for (int cc = 0; cc < (int)splitClustDirSet.size(); cc++) { + GlobalVector bigClustDir = splitClustDirSet[cc]; + + jetEta_ = jet.eta(); + jetPt_ = jet.pt(); + + const auto& jetVert = jetVertex; //trackInfo filling + + std::vector goodSimHit; + + const GeomDet* globDet = DetectorSelector( + 2, jet, bigClustDir, jetVertex, tTopo, inputPixelClusters_); //select detector mostly hitten by the jet + + if (globDet == nullptr) + continue; + + std::pair, std::vector> goodSimTkVx; + + if (inclusiveConeSeed_) { + goodSimTkVx = JetCoreMCtruthSeedGenerator::coreTracksFillingDeltaR( + simtracksVector, simvertexVector, globDet, jet, jetVert); + } else { + std::vector goodSimHit = + JetCoreMCtruthSeedGenerator::coreHitsFilling(simhits_, globDet, bigClustDir, jetVertex); + goodSimTkVx = JetCoreMCtruthSeedGenerator::coreTracksFilling(goodSimHit, simtracksVector, simvertexVector); + } + edm::LogInfo("PerfectSeeder") << "seed number in deltaR cone =" << goodSimTkVx.first.size(); + + std::vector> seedVector = + JetCoreMCtruthSeedGenerator::seedParFilling(goodSimTkVx, globDet, jet); + edm::LogInfo("PerfectSeeder") << "seedVector.size()=" << seedVector.size(); + + for (uint tk = 0; tk < seedVector.size(); tk++) { + for (int pp = 0; pp < 5; pp++) { + edm::LogInfo("PerfectSeeder") + << "seed " << tk << ", int par " << pp << "=" << seedVector[tk][pp] << std::endl; + } + LocalPoint localSeedPoint = LocalPoint(seedVector[tk][0], seedVector[tk][1], 0); + double track_theta = 2 * std::atan(std::exp(-seedVector[tk][2])); + double track_phi = seedVector[tk][3]; + double pt = 1. / seedVector[tk][4]; + + double normdirR = pt / sin(track_theta); + const GlobalVector globSeedDir( + GlobalVector::Polar(Geom::Theta(track_theta), Geom::Phi(track_phi), normdirR)); + LocalVector localSeedDir = globDet->surface().toLocal(globSeedDir); + + int64_t seedid = (int64_t(localSeedPoint.x() * 200.) << 0) + (int64_t(localSeedPoint.y() * 200.) << 16) + + (int64_t(seedVector[tk][2] * 400.) << 32) + (int64_t(track_phi * 400.) << 48); + if (ids.count(seedid) != 0) { + edm::LogInfo("PerfectSeeder") << "seed not removed with DeepCore cleaner"; + } + ids.insert(seedid); + + //Covariance matrix, currently the hadrcoded variances = NN residuals width (see documentation of DeepCore) + //in general: if are not compared with DeepCore but another algo-->to state-of-the art errors + //The "perfect seeds" has no intrinsic error, but the CTF needs errors to propagate... + float em[15] = { + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; // (see LocalTrajectoryError for details), order as follow: + em[0] = 0.15 * 0.15; // q/pt + em[2] = 0.5e-5; // dxdz + em[5] = 0.5e-5; // dydz + em[9] = 2e-5; // x + em[14] = 2e-5; // y + long int detId = globDet->geographicalId(); + LocalTrajectoryParameters localParam(localSeedPoint, localSeedDir, TrackCharge(1)); + result->emplace_back(TrajectorySeed(PTrajectoryStateOnDet(localParam, pt, em, detId, /*surfaceSide*/ 0), + edm::OwnVector(), + PropagationDirection::alongMomentum)); + + GlobalPoint globalSeedPoint = globDet->surface().toGlobal(localSeedPoint); + reco::Track::CovarianceMatrix mm; + resultTracks->emplace_back( + reco::Track(1, + 1, + reco::Track::Point(globalSeedPoint.x(), globalSeedPoint.y(), globalSeedPoint.z()), + reco::Track::Vector(globSeedDir.x(), globSeedDir.y(), globSeedDir.z()), + 1, + mm)); + edm::LogInfo("PerfectSeeder") << "seed " << tk << ", out, pt=" << pt << ", eta=" << globSeedDir.eta() + << ", phi=" << globSeedDir.phi() << std::endl; + } + + } //bigcluster + } //jet > pt + } //jet + iEvent.put(std::move(result)); + iEvent.put(std::move(resultTracks)); +} + +std::pair> JetCoreMCtruthSeedGenerator::findIntersection( + const GlobalVector& dir, const reco::Candidate::Point& vertex, const GeomDet* det) { + StraightLinePlaneCrossing vertexPlane(Basic3DVector(vertex.x(), vertex.y(), vertex.z()), + Basic3DVector(dir.x(), dir.y(), dir.z())); + + std::pair> pos = vertexPlane.position(det->specificSurface()); + + return pos; +} + +const GeomDet* JetCoreMCtruthSeedGenerator::DetectorSelector(int llay, + const reco::Candidate& jet, + GlobalVector jetDir, + const reco::Vertex& jetVertex, + const TrackerTopology* const tTopo, + const edmNew::DetSetVector& clusters) { + struct trkNumCompare { + bool operator()(std::pair x, std::pair y) const { + return x.first > y.first; + } + }; + std::set, trkNumCompare> track4detSet; + + double minDist = 0.0; + GeomDet* output = (GeomDet*)nullptr; + for (const auto& detset : clusters) { + auto aClusterID = detset.id(); + if (DetId(aClusterID).subdetId() != 1) + continue; + const GeomDet* det = geometry_->idToDet(aClusterID); + int lay = tTopo->layer(det->geographicalId()); + if (lay != llay) + continue; + std::pair> interPair = + findIntersection(jetDir, (reco::Candidate::Point)jetVertex.position(), det); + if (interPair.first == false) + continue; + Basic3DVector inter = interPair.second; + auto localInter = det->specificSurface().toLocal((GlobalPoint)inter); + if ((minDist == 0.0 || std::abs(localInter.x()) < minDist) && std::abs(localInter.y()) < 3.35) { + minDist = std::abs(localInter.x()); + output = (GeomDet*)det; + } + } //detset + return output; +} + +std::vector JetCoreMCtruthSeedGenerator::splittedClusterDirections( + const reco::Candidate& jet, + const TrackerTopology* const tTopo, + const PixelClusterParameterEstimator* pixelCPE, + const reco::Vertex& jetVertex, + int layer, + const edmNew::DetSetVector& clusters) { + std::vector clustDirs; + for (const auto& detset_int : clusters) { + const GeomDet* det_int = geometry_->idToDet(detset_int.id()); + int lay = tTopo->layer(det_int->geographicalId()); + if (lay != layer) + continue; //NB: saved bigclusetr on all the layers!! + auto detUnit = *geometry_->idToDetUnit(detset_int.id()); + for (const auto& aCluster : detset_int) { + GlobalPoint clustPos = det_int->surface().toGlobal(pixelCPE->localParametersV(aCluster, detUnit)[0].first); + GlobalPoint vertexPos(jetVertex.position().x(), jetVertex.position().y(), jetVertex.position().z()); + GlobalVector clusterDir = clustPos - vertexPos; + GlobalVector jetDir(jet.px(), jet.py(), jet.pz()); + if (Geom::deltaR(jetDir, clusterDir) < deltaR_) { + clustDirs.emplace_back(clusterDir); + } + } + } + return clustDirs; +} + +std::vector JetCoreMCtruthSeedGenerator::coreHitsFilling(std::vector simhits, + const GeomDet* globDet, + GlobalVector bigClustDir, + const reco::Vertex& jetVertex) { + std::vector goodSimHit; + for (const auto& sh : simhits) { + const GeomDet* det = geometry_->idToDet(sh.detUnitId()); + if (det != globDet) + continue; + std::pair> interPair = + findIntersection(bigClustDir, (reco::Candidate::Point)jetVertex.position(), det); + if (interPair.first == false) + continue; + Basic3DVector inter = interPair.second; + auto localInter = det->specificSurface().toLocal((GlobalPoint)inter); + + if (std::abs((sh.localPosition()).x() - localInter.x()) / pitchX_ <= jetDimX / 2 && + std::abs((sh.localPosition()).y() - localInter.y()) / pitchY_ <= jetDimY / 2) { + goodSimHit.emplace_back(sh); + } + } + return goodSimHit; +} + +std::pair, std::vector> JetCoreMCtruthSeedGenerator::coreTracksFilling( + std::vector goodSimHit, + const std::vector simtracksVector, + const std::vector simvertexVector) { + std::vector goodSimTrk; + std::vector goodSimVtx; + + for (uint j = 0; j < simtracksVector.size(); j++) { + for (std::vector::const_iterator it = goodSimHit.begin(); it != goodSimHit.end(); ++it) { + SimTrack st = simtracksVector[j]; + if (st.trackId() == (*it).trackId()) { + for (uint v = 0; v < simvertexVector.size(); v++) { + SimVertex sv = simvertexVector[v]; + if ((int)sv.vertexId() == (int)st.vertIndex()) { + goodSimTrk.emplace_back(st); + goodSimVtx.emplace_back(sv); + } + } + } + } + } + std::pair, std::vector> output(goodSimTrk, goodSimVtx); + return output; +} + +std::pair, std::vector> JetCoreMCtruthSeedGenerator::coreTracksFillingDeltaR( + const std::vector simtracksVector, + const std::vector simvertexVector, + const GeomDet* globDet, + const reco::Candidate& jet, + const reco::Vertex& jetVertex) { + std::vector goodSimTrk; + std::vector goodSimVtx; + + GlobalVector jetDir(jet.px(), jet.py(), jet.pz()); + + for (uint j = 0; j < simtracksVector.size(); j++) { + SimTrack st = simtracksVector[j]; + GlobalVector trkDir(st.momentum().Px(), st.momentum().Py(), st.momentum().Pz()); + if (st.charge() == 0) + continue; + if (Geom::deltaR(jetDir, trkDir) < deltaR_) { + for (uint v = 0; v < simvertexVector.size(); v++) { + SimVertex sv = simvertexVector[v]; + if ((int)sv.vertexId() == (int)st.vertIndex()) { + goodSimTrk.emplace_back(st); + goodSimVtx.emplace_back(sv); + } + } + } + } + std::pair, std::vector> output(goodSimTrk, goodSimVtx); + return output; +} + +std::vector> JetCoreMCtruthSeedGenerator::seedParFilling( + std::pair, std::vector> goodSimTkVx, + const GeomDet* globDet, + const reco::Candidate& jet) { + std::vector> output; + std::vector goodSimTrk = goodSimTkVx.first; + std::vector goodSimVtx = goodSimTkVx.second; + + edm::LogInfo("PerfectSeeder") << "goodSimTrk size" << goodSimTrk.size(); + for (uint j = 0; j < goodSimTrk.size(); j++) { + SimTrack st = goodSimTrk[j]; + SimVertex sv = goodSimVtx[j]; + GlobalVector trkMom(st.momentum().x(), st.momentum().y(), st.momentum().z()); + GlobalPoint trkPos(sv.position().x(), sv.position().y(), sv.position().z()); + edm::LogInfo("PerfectSeeder") << "seed " << j << ", very int pt" << st.momentum().Pt() + << ", eta=" << st.momentum().Eta() << ", phi=" << st.momentum().Phi() + << "------ internal point=" << trkMom.x() << "," << trkMom.y() << "," << trkMom.z() + << "," << trkPos.x() << "," << trkPos.y() << "," << trkPos.z() << std::endl; + + std::pair> trkInterPair; + trkInterPair = findIntersection(trkMom, (reco::Candidate::Point)trkPos, globDet); + if (trkInterPair.first == false) { + GlobalVector jetDir(jet.px(), jet.py(), jet.pz()); + continue; + } + Basic3DVector trkInter = trkInterPair.second; + + auto localTrkInter = globDet->specificSurface().toLocal((GlobalPoint)trkInter); //trkInter->trkPos if par at vertex + std::array tkPar{ + {localTrkInter.x(), localTrkInter.y(), st.momentum().Eta(), st.momentum().Phi(), 1 / st.momentum().Pt()}}; + output.emplace_back(tkPar); + } + return output; +} + +// ------------ method called once each job just before starting event loop ------------ +void JetCoreMCtruthSeedGenerator::beginJob() {} + +// ------------ method called once each job just after ending the event loop ------------ +void JetCoreMCtruthSeedGenerator::endJob() {} + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +void JetCoreMCtruthSeedGenerator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("vertices", edm::InputTag("offlinePrimaryVertices")); + desc.add("pixelClusters", edm::InputTag("siPixelClustersPreSplitting")); + desc.add("cores", edm::InputTag("jetsForCoreTracking")); + desc.add("ptMin", 300); + desc.add("deltaR", 0.3); + desc.add("chargeFractionMin", 18000.0); + desc.add("simTracks", edm::InputTag("g4SimHits")); + desc.add("simVertex", edm::InputTag("g4SimHits")); + desc.add("simHit", edm::InputTag("g4SimHits", "TrackerHitsPixelBarrelLowTof")); + desc.add("centralMIPCharge", 2.); + desc.add("pixelCPE", "PixelCPEGeneric"); + descriptions.add("JetCoreMCtruthSeedGenerator", desc); +} + +DEFINE_FWK_MODULE(JetCoreMCtruthSeedGenerator); diff --git a/Validation/RecoTrack/plugins/TrackFromSeedProducer.cc b/Validation/RecoTrack/plugins/TrackFromSeedProducer.cc index 091bca71dd0ec..e5627e2811db7 100644 --- a/Validation/RecoTrack/plugins/TrackFromSeedProducer.cc +++ b/Validation/RecoTrack/plugins/TrackFromSeedProducer.cc @@ -40,6 +40,8 @@ #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" #include "DataFormats/BeamSpot/interface/BeamSpot.h" +#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h" +#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h" // // class declaration @@ -59,6 +61,7 @@ class TrackFromSeedProducer : public edm::global::EDProducer<> { edm::EDGetTokenT > seedsToken; edm::EDGetTokenT beamSpotToken; std::string tTRHBuilderName; + const edm::ESGetToken geoToken_; }; // @@ -72,7 +75,8 @@ class TrackFromSeedProducer : public edm::global::EDProducer<> { // // constructors and destructor // -TrackFromSeedProducer::TrackFromSeedProducer(const edm::ParameterSet& iConfig) { +TrackFromSeedProducer::TrackFromSeedProducer(const edm::ParameterSet& iConfig) + : geoToken_(esConsumes()) { //register your products produces(); produces(); @@ -127,14 +131,21 @@ void TrackFromSeedProducer::produce(edm::StreamID, edm::Event& iEvent, const edm iSetup.get().get(httopo); const TrackerTopology& ttopo = *httopo; + const GlobalTrackingGeometry* const geometry_ = &iSetup.getData(geoToken_); + // create tracks from seeds int nfailed = 0; for (size_t iSeed = 0; iSeed < seeds.size(); ++iSeed) { auto const& seed = seeds[iSeed]; // try to create a track - TransientTrackingRecHit::RecHitPointer lastRecHit = tTRHBuilder->build(&*(seed.recHits().end() - 1)); - TrajectoryStateOnSurface state = - trajectoryStateTransform::transientState(seed.startingState(), lastRecHit->surface(), theMF.product()); + TrajectoryStateOnSurface state; + if (seed.nHits() == 0) { //this is for deepCore seeds only + const Surface* deepCore_sruface = &geometry_->idToDet(seed.startingState().detId())->specificSurface(); + state = trajectoryStateTransform::transientState(seed.startingState(), deepCore_sruface, theMF.product()); + } else { + TransientTrackingRecHit::RecHitPointer lastRecHit = tTRHBuilder->build(&*(seed.recHits().end() - 1)); + state = trajectoryStateTransform::transientState(seed.startingState(), lastRecHit->surface(), theMF.product()); + } TrajectoryStateClosestToBeamLine tsAtClosestApproachSeed = tscblBuilder(*state.freeState(), *beamSpot); //as in TrackProducerAlgorithm if (tsAtClosestApproachSeed.isValid()) { diff --git a/Validation/RecoTrack/python/MTVHistoProducerAlgoForTrackerBlock_cfi.py b/Validation/RecoTrack/python/MTVHistoProducerAlgoForTrackerBlock_cfi.py index 9ca553b9aafd3..7ce1eb5fad5c7 100644 --- a/Validation/RecoTrack/python/MTVHistoProducerAlgoForTrackerBlock_cfi.py +++ b/Validation/RecoTrack/python/MTVHistoProducerAlgoForTrackerBlock_cfi.py @@ -83,7 +83,7 @@ # dR_jet mindrj = cms.double(0.001), maxdrj = cms.double(0.5), - nintdrj = cms.int32(250), + nintdrj = cms.int32(100), # # chi2/ndof minChi2 = cms.double(0), diff --git a/Validation/RecoTrack/python/PostProcessorTracker_cfi.py b/Validation/RecoTrack/python/PostProcessorTracker_cfi.py index cec8d3cd3cff9..36a69eca3cccd 100644 --- a/Validation/RecoTrack/python/PostProcessorTracker_cfi.py +++ b/Validation/RecoTrack/python/PostProcessorTracker_cfi.py @@ -138,21 +138,11 @@ def _addNoFlow(module): "fakerate_vs_dzpvcut 'Fake rate vs. dz(PV)' num_assoc(recoToSim)_dzpvcut num_reco_dzpvcut fake", "pileuprate_dzpvcut 'Pileup rate vs. dz(PV)' num_pileup_dzpvcut num_reco_dzpvcut", - "effic_vs_dzpvcut_pt 'Fraction of true p_{T} carried by recoed TPs from PV vs. dz(PV)' num_assoc(simToReco)_dzpvcut_pt num_simul_dzpvcut_pt", - "effic_vs_dzpvcut2_pt 'Fraction of true p_{T} carried by recoed TPs from PV (tracking eff factorized out) vs. dz(PV)' num_assoc(simToReco)_dzpvcut_pt num_simul2_dzpvcut_pt", - "fakerate_vs_dzpvcut_pt 'Fraction of fake p_{T} carried by tracks from PV vs. dz(PV)' num_assoc(recoToSim)_dzpvcut_pt num_reco_dzpvcut_pt fake", - "pileuprate_dzpvcut_pt 'Fraction of pileup p_{T} carried by tracks from PV vs. dz(PV)' num_pileup_dzpvcut_pt num_reco_dzpvcut_pt", - "effic_vs_dzpvsigcut 'Efficiency vs. dz(PV)/dzError' num_assoc(simToReco)_dzpvsigcut num_simul_dzpvsigcut", "effic_vs_dzpvsigcut2 'Efficiency (tracking eff factorized out) vs. dz(PV)/dzError' num_assoc(simToReco)_dzpvsigcut num_simul2_dzpvsigcut", "fakerate_vs_dzpvsigcut 'Fake rate vs. dz(PV)/dzError' num_assoc(recoToSim)_dzpvsigcut num_reco_dzpvsigcut fake", "pileuprate_dzpvsigcut 'Pileup rate vs. dz(PV)/dzError' num_pileup_dzpvsigcut num_reco_dzpvsigcut", - "effic_vs_dzpvsigcut_pt 'Fraction of true p_{T} carried by recoed TPs from PV vs. dz(PV)/dzError' num_assoc(simToReco)_dzpvsigcut_pt num_simul_dzpvsigcut_pt", - "effic_vs_dzpvsigcut2_pt 'Fraction of true p_{T} carried by recoed TPs from PV (tracking eff factorized out) vs. dz(PV)/dzError' num_assoc(simToReco)_dzpvsigcut_pt num_simul2_dzpvsigcut_pt", - "fakerate_vs_dzpvsigcut_pt 'Fraction of fake p_{T} carried by tracks from PV vs. dz(PV)/dzError' num_assoc(recoToSim)_dzpvsigcut_pt num_reco_dzpvsigcut_pt fake", - "pileuprate_dzpvsigcut_pt 'Fraction of pileup p_{T} carried by tracks from PV vs. dz(PV)/dzError' num_pileup_dzpvsigcut_pt num_reco_dzpvsigcut_pt", - "effic_vs_simpvz 'Efficiency vs. sim PV z' num_assoc(simToReco)_simpvz num_simul_simpvz", "fakerate_vs_simpvz 'Fake rate vs. sim PV z' num_assoc(recoToSim)_simpvz num_reco_simpvz fake", "duplicatesRate_simpvz 'Duplicates Rate vs sim PV z' num_duplicate_simpvz num_reco_simpvz", @@ -210,24 +200,12 @@ def _addNoFlow(module): "num_simul_dzpvcut", "num_simul2_dzpvcut", "num_pileup_dzpvcut", - "num_reco_dzpvcut_pt", - "num_assoc(recoToSim)_dzpvcut_pt", - "num_assoc(simToReco)_dzpvcut_pt", - "num_simul_dzpvcut_pt", - "num_simul2_dzpvcut_pt", - "num_pileup_dzpvcut_pt", "num_reco_dzpvsigcut", "num_assoc(recoToSim)_dzpvsigcut", "num_assoc(simToReco)_dzpvsigcut", "num_simul_dzpvsigcut", "num_simul2_dzpvsigcut", "num_pileup_dzpvsigcut", - "num_reco_dzpvsigcut_pt", - "num_assoc(recoToSim)_dzpvsigcut_pt", - "num_assoc(simToReco)_dzpvsigcut_pt", - "num_simul_dzpvsigcut_pt", - "num_simul2_dzpvsigcut_pt", - "num_pileup_dzpvsigcut_pt", "num_reco_mva1cut descending", "num_reco_mva2cut descending", "num_reco_mva2cut_hp descending", @@ -317,6 +295,18 @@ def _addNoFlow(module): postProcessorTrackSummary ) +from Configuration.ProcessModifiers.seedingDeepCore_cff import seedingDeepCore +postProcessorTrackDeepCore = postProcessorTrack.clone() +postProcessorTrackDeepCore.subDirs.extend(["Tracking/JetCore/*"]) +seedingDeepCore.toReplaceWith(postProcessorTrack,postProcessorTrackDeepCore) +postProcessorTrackSummaryDeepCore = postProcessorTrackSummary.clone() +postProcessorTrackSummaryDeepCore.subDirs.extend(["Tracking/JetCore/*"]) +seedingDeepCore.toReplaceWith(postProcessorTrackSummary,postProcessorTrackSummaryDeepCore) +postProcessorTrack2DDeepCore = postProcessorTrack2D.clone() +postProcessorTrack2DDeepCore.subDirs.extend(["Tracking/JetCore/*"]) +seedingDeepCore.toReplaceWith(postProcessorTrack2D,postProcessorTrack2DDeepCore) + + fastSim.toModify(postProcessorTrack, subDirs = [e for e in _defaultSubdirs if e not in ["Tracking/TrackGsf/*","Tracking/TrackConversion/*"]]) fastSim.toModify(postProcessorTrackSummary, subDirs = [e for e in _defaultSubdirsSummary if e not in ["Tracking/TrackGsf","Tracking/TrackConversion"]]) diff --git a/Validation/RecoTrack/python/TrackValidation_cff.py b/Validation/RecoTrack/python/TrackValidation_cff.py index a90f3eac8a605..fa0c3d1e3c578 100644 --- a/Validation/RecoTrack/python/TrackValidation_cff.py +++ b/Validation/RecoTrack/python/TrackValidation_cff.py @@ -353,14 +353,9 @@ def _getMVASelectors(postfix): ptMin = 0, ) -#ByChi2 association (for jetCore usage, not used by default) -MTVTrackAssociationByChi2 = trackingParticleRecoTrackAsssociation.clone( - associator = cms.InputTag('trackAssociatorByChi2') -) - # Select jets for JetCore tracking -highPtJets = cms.EDFilter("CandPtrSelector", src = cms.InputTag("ak4CaloJets"), cut = cms.string("pt()>1000")) -highPtJetsForTrk = highPtJetsForTrk = highPtJets.clone(src = "ak4CaloJetsForTrk") +highPtJets = cms.EDFilter("CandPtrSelector", src = cms.InputTag("ak4CaloJets"), cut = cms.string("pt()>1000")) +highPtJetsForTrk = highPtJets.clone(src = "ak4CaloJetsForTrk") # Select B-hadron TPs trackingParticlesBHadron = _trackingParticleBHadronRefSelector.clone() @@ -369,7 +364,6 @@ def _getMVASelectors(postfix): trackValidator = Validation.RecoTrack.MultiTrackValidator_cfi.multiTrackValidator.clone( useLogPt = cms.untracked.bool(True), dodEdxPlots = True, - # associators=cms.untracked.VInputTag('MTVTrackAssociationByChi2'), #uncomment for byChi2 assoc. for jetcore studies (1/5) doPVAssociationPlots = True #,minpT = cms.double(-1) #,maxpT = cms.double(3) @@ -540,8 +534,6 @@ def _getMVASelectors(postfix): dirName = "Tracking/TrackBuilding/", doMVAPlots = True, doResolutionPlotsForLabels = ['jetCoreRegionalStepTracks'], - # associators = ["trackAssociatorByChi2"], #uncomment for byChi2 assoc. for jetcore studies (2/5) - # UseAssociators = True, #uncomment for byChi2 assoc. for jetcore studies (3/5) ) trackValidatorBuildingPreSplitting = trackValidatorBuilding.clone( associators = ["quickTrackAssociatorByHitsPreSplitting"], @@ -606,7 +598,25 @@ def _uniqueFirstLayers(layerList): for _eraName, _postfix, _era in _relevantEras: _setForEra(trackValidatorGsfTracks.histoProducerAlgoBlock, _eraName, _era, seedingLayerSets=trackValidator.histoProducerAlgoBlock.seedingLayerSets.value()+locals()["_seedingLayerSetsForElectrons"+_postfix]) - +# For jetCore tracks +trackValidatorJetCore = trackValidator.clone(#equivalent to trackBuilding case + dirName = "Tracking/JetCore/", + useLogPt = cms.untracked.bool(True), + dodEdxPlots = False, + associators= ["trackAssociatorByChi2"],#cms.untracked.VInputTag('MTVTrackAssociationByChi2'), + UseAssociators = True, + doPVAssociationPlots = True, +) +for _eraName, _postfix, _era in _relevantEras: + if 'jetCoreRegionalStep' in _cfg.iterationAlgos(_postfix) : + _setForEra(trackValidatorJetCore, _eraName, _era, + label = ["generalTracks", "jetCoreRegionalStepTracks", + "cutsRecoTracksJetCoreRegionalStepByOriginalAlgo","cutsRecoTracksJetCoreRegionalStepByOriginalAlgoHp", + "cutsRecoTracksJetCoreRegionalStep", "cutsRecoTracksJetCoreRegionalStepHp"], + doResolutionPlotsForLabels =["generalTracks", "jetCoreRegionalStepTracks", + "cutsRecoTracksJetCoreRegionalStepByOriginalAlgo","cutsRecoTracksJetCoreRegionalStepByOriginalAlgoHp", + "cutsRecoTracksJetCoreRegionalStep", "cutsRecoTracksJetCoreRegionalStepHp"], + ) # for B-hadrons trackValidatorBHadron = trackValidator.clone( @@ -658,8 +668,7 @@ def _uniqueFirstLayers(layerList): tracksValidationTruth = cms.Task( tpClusterProducer, tpClusterProducerPreSplitting, - # trackAssociatorByChi2, #uncomment for byChi2 assoc. for jetcore studies (4/5) - # MTVTrackAssociationByChi2, #uncomment for byChi2 assoc. for jetcore studies (5/5) + trackAssociatorByChi2, quickTrackAssociatorByHits, quickTrackAssociatorByHitsPreSplitting, trackingParticleRecoTrackAsssociation, @@ -699,13 +708,18 @@ def _uniqueFirstLayers(layerList): tracksPreValidation ) +from Configuration.ProcessModifiers.seedingDeepCore_cff import seedingDeepCore +seedingDeepCore.toReplaceWith(tracksValidation, cms.Sequence(tracksValidation.copy()+trackValidatorJetCore)) + from Configuration.Eras.Modifier_phase2_tracker_cff import phase2_tracker #tracksValidationPhase2 = cms.Sequence(tracksValidation+trackValidatorTPEtaGreater2p7) # it does not work tracksPreValidationPhase2 = tracksPreValidation.copy() tracksPreValidationPhase2.add(trackingParticlesEtaGreater2p7) phase2_tracker.toReplaceWith(tracksPreValidation, tracksPreValidationPhase2) -tracksValidationPhase2 = tracksValidation.copy() +tracksValidationPhase2 = tracksValidation.copyAndExclude([ + trackValidatorJetCore +]) tracksValidationPhase2+=trackValidatorTPEtaGreater2p7 phase2_tracker.toReplaceWith(tracksValidation, tracksValidationPhase2) @@ -849,6 +863,21 @@ def _uniqueFirstLayers(layerList): doSummaryPlots = False, ) + +trackValidatorJetCoreSeedingTrackingOnly = trackValidatorSeedingTrackingOnly.clone( + dirName = "Tracking/JetCore/TrackSeeding/", + associators = ["trackAssociatorByChi2"], + UseAssociators = True, + doSeedPlots = True, +) + +for _eraName, _postfix, _era in _relevantEras: + if 'jetCoreRegionalStep' in _cfg.iterationAlgos(_postfix) : + _setForEra(trackValidatorJetCoreSeedingTrackingOnly, _eraName, _era, + label = [ "seedTracksjetCoreRegionalStepSeeds",], + doResolutionPlotsForLabels = [ "seedTracksjetCoreRegionalStepSeeds",] + ) + for _eraName, _postfix, _era in _relevantErasAndFastSim: _setForEra(trackValidatorSeedingTrackingOnly, _eraName, _era, label = locals()["_seedSelectors"+_postfix]) for _eraName, _postfix, _era in _relevantEras: @@ -881,13 +910,23 @@ def _uniqueFirstLayers(layerList): trackValidatorsTrackingOnly.replace(trackValidatorConversionStandalone, trackValidatorConversionTrackingOnly) trackValidatorsTrackingOnly.remove(trackValidatorGsfTracksStandalone) trackValidatorsTrackingOnly.replace(trackValidatorBHadronStandalone, trackValidatorBHadronTrackingOnly) + +seedingDeepCore.toReplaceWith(trackValidatorsTrackingOnly, cms.Sequence( + trackValidatorsTrackingOnly.copy()+ + trackValidatorJetCore+ + trackValidatorJetCoreSeedingTrackingOnly + ) + ) +phase2_tracker.toReplaceWith(trackValidatorsTrackingOnly, trackValidatorsTrackingOnly.copyAndExclude([ #must be done for each era which does not have jetcore in the iteration + trackValidatorJetCore, + trackValidatorJetCoreSeedingTrackingOnly +])) fastSim.toReplaceWith(trackValidatorsTrackingOnly, trackValidatorsTrackingOnly.copyAndExclude([ trackValidatorBuildingPreSplitting, trackValidatorSeedingPreSplittingTrackingOnly, trackValidatorConversionTrackingOnly, trackValidatorBHadronTrackingOnly ])) - tracksValidationTrackingOnly = cms.Sequence( trackValidatorsTrackingOnly, tracksPreValidationTrackingOnly, diff --git a/Validation/RecoTrack/python/plotting/trackingPlots.py b/Validation/RecoTrack/python/plotting/trackingPlots.py index 374c83cf50ac8..dea68a2a7de13 100644 --- a/Validation/RecoTrack/python/plotting/trackingPlots.py +++ b/Validation/RecoTrack/python/plotting/trackingPlots.py @@ -297,16 +297,6 @@ def _makeMVAPlots(num, hp=False): xtitle="Efficiency vs. cut on dz(PV)/dzError", **_common), Plot(ROC("effic_vs_fakepileup2_dzpvsigcut", "effic_vs_dzpvsigcut", FakeDuplicate("fakepileup_vs_dzpvsigcut", assoc="num_assoc(recoToSim)_dzpvsigcut", reco="num_reco_dzpvsigcut", dup="num_pileup_dzpvsigcut"), zaxis=True), xtitle="Efficiency", ztitle="Cut on dz(PV)/dzError", **_common2), - ## - Plot(ROC("effic_vs_fakepileup_dzpvcut_pt", "effic_vs_dzpvcut_pt", FakeDuplicate("fakepileup_vs_dzpvcut_pt", assoc="num_assoc(recoToSim)_dzpvcut_pt", reco="num_reco_dzpvcut_pt", dup="num_pileup_dzpvcut_pt")), - xtitle="Efficiency (p_{T} weighted) vs. cut on dz(PV)", **_common), - Plot(ROC("effic_vs_fakepileup2_dzpvcut_pt", "effic_vs_dzpvcut_pt", FakeDuplicate("fakepileup_vs_dzpvcut_pt", assoc="num_assoc(recoToSim)_dzpvcut_pt", reco="num_reco_dzpvcut_pt", dup="num_pileup_dzpvcut_pt"), zaxis=True), - xtitle="Efficiency (p_{T} weighted)", ztitle="Cut on dz(PV)", **_common2), - # - Plot(ROC("effic_vs_fakepileup_dzpvsigcut_pt", "effic_vs_dzpvsigcut_pt", FakeDuplicate("fakepileup_vs_dzpvsigcut_pt", assoc="num_assoc(recoToSim)_dzpvsigcut_pt", reco="num_reco_dzpvsigcut_pt", dup="num_pileup_dzpvsigcut_pt")), - xtitle="Efficiency (p_{T} weighted) vs. cut on dz(PV)/dzError", **_common), - Plot(ROC("effic_vs_fakepileup2_dzpvsigcut_pt", "effic_vs_dzpvsigcut_pt", FakeDuplicate("fakepileup_vs_dzpvsigcut_pt", assoc="num_assoc(recoToSim)_dzpvsigcut_pt", reco="num_reco_dzpvsigcut_pt", dup="num_pileup_dzpvsigcut_pt"), zaxis=True), - xtitle="Efficiency (p_{T} weighted)", ztitle="Cut on dz(PV)/dzError", **_common2), ], onlyForPileup=True, legendDy=_legendDy_4rows ) @@ -323,19 +313,6 @@ def _makeMVAPlots(num, hp=False): ], onlyForPileup=True, legendDy=_legendDy_4rows ) -_pvassociation3 = PlotGroup("pvassociation3", [ - Plot("effic_vs_dzpvcut_pt", xtitle="Cut on dz(PV) (cm)", ytitle="Efficiency (p_{T} weighted)", ymax=_maxEff), - Plot("effic_vs_dzpvcut2_pt", xtitle="Cut on dz(PV) (cm)", ytitle="Efficiency (p_{T} weighted, excl. trk eff)", ymax=_maxEff), - Plot("fakerate_vs_dzpvcut_pt", xtitle="Cut on dz(PV) (cm)", ytitle="Fake rate (p_{T} weighted)", ymax=_maxFake), - Plot("pileuprate_dzpvcut_pt", xtitle="Cut on dz(PV) (cm)", ytitle="Pileup rate (p_{T} weighted)", ymax=_maxFake), - # - Plot("effic_vs_dzpvsigcut_pt", xtitle="Cut on dz(PV)/dzError", ytitle="Efficiency (p_{T} weighted)", ymax=_maxEff), - Plot("effic_vs_dzpvsigcut2_pt", xtitle="Cut on dz(PV)/dzError", ytitle="Efficiency (p_{T} weighted, excl. trk eff)", ymax=_maxEff), - Plot("fakerate_vs_dzpvsigcut_pt", xtitle="Cut on dz(PV)/dzError", ytitle="Fake rate (p_{T} weighted)", ymax=_maxFake), - Plot("pileuprate_dzpvsigcut_pt", xtitle="Cut on dz(PV)/dzError", ytitle="Pileup rate (p_{T} weighted)", ymax=_maxFake), -], onlyForPileup=True, - legendDy=_legendDy_4rows -) # These don't exist in FastSim @@ -1223,7 +1200,6 @@ def _trackingFolders(lastDirName="Track"): _dupandfakeSeedingTable, _pvassociation1, _pvassociation2, - _pvassociation3, _dedx, # _chargemisid, _hitsAndPt, diff --git a/Validation/RecoTrack/src/MTVHistoProducerAlgoForTracker.cc b/Validation/RecoTrack/src/MTVHistoProducerAlgoForTracker.cc index 8664c23ed5abb..e3370837a2a4c 100644 --- a/Validation/RecoTrack/src/MTVHistoProducerAlgoForTracker.cc +++ b/Validation/RecoTrack/src/MTVHistoProducerAlgoForTracker.cc @@ -584,23 +584,6 @@ void MTVHistoProducerAlgoForTracker::bookSimTrackPVAssociationHistos(DQMStore::I histograms.h_simul2_dzpvcut.push_back(ibook.book1D( "num_simul2_dzpvcut", "N of simulated tracks (associated to any track) from sim PV", nintDzpvCum, 0, maxDzpvCum)); - histograms.h_assoc_dzpvcut_pt.push_back(ibook.book1D("num_assoc(simToReco)_dzpvcut_pt", - "#sump_{T} of associated tracks (simToReco) vs dz(PV)", - nintDzpvCum, - 0, - maxDzpvCum)); - histograms.h_simul_dzpvcut_pt.push_back( - ibook.book1D("num_simul_dzpvcut_pt", "#sump_{T} of simulated tracks from sim PV", nintDzpvCum, 0, maxDzpvCum)); - histograms.h_simul2_dzpvcut_pt.push_back( - ibook.book1D("num_simul2_dzpvcut_pt", - "#sump_{T} of simulated tracks (associated to any track) from sim PV", - nintDzpvCum, - 0, - maxDzpvCum)); - histograms.h_assoc_dzpvcut_pt.back()->enableSumw2(); - histograms.h_simul_dzpvcut_pt.back()->enableSumw2(); - histograms.h_simul2_dzpvcut_pt.back()->enableSumw2(); - histograms.h_assoc_dzpvsigcut.push_back(ibook.book1D("num_assoc(simToReco)_dzpvsigcut", "N of associated tracks (simToReco) vs dz(PV)/dzError", nintDzpvsigCum, @@ -614,24 +597,6 @@ void MTVHistoProducerAlgoForTracker::bookSimTrackPVAssociationHistos(DQMStore::I nintDzpvsigCum, 0, maxDzpvsigCum)); - - histograms.h_assoc_dzpvsigcut_pt.push_back( - ibook.book1D("num_assoc(simToReco)_dzpvsigcut_pt", - "#sump_{T} of associated tracks (simToReco) vs dz(PV)/dzError", - nintDzpvsigCum, - 0, - maxDzpvsigCum)); - histograms.h_simul_dzpvsigcut_pt.push_back(ibook.book1D( - "num_simul_dzpvsigcut_pt", "#sump_{T} of simulated tracks from sim PV/dzError", nintDzpvsigCum, 0, maxDzpvsigCum)); - histograms.h_simul2_dzpvsigcut_pt.push_back( - ibook.book1D("num_simul2_dzpvsigcut_pt", - "#sump_{T} of simulated tracks (associated to any track) from sim PV/dzError", - nintDzpvsigCum, - 0, - maxDzpvsigCum)); - histograms.h_assoc_dzpvsigcut_pt.back()->enableSumw2(); - histograms.h_simul_dzpvsigcut_pt.back()->enableSumw2(); - histograms.h_simul2_dzpvsigcut_pt.back()->enableSumw2(); } void MTVHistoProducerAlgoForTracker::bookRecoHistos(DQMStore::IBooker& ibook, @@ -1533,22 +1498,6 @@ void MTVHistoProducerAlgoForTracker::bookRecoPVAssociationHistos(DQMStore::IBook histograms.h_pileup_dzpvcut.push_back(ibook.book1D( "num_pileup_dzpvcut", "N of associated (recoToSim) pileup tracks vs dz(PV)", nintDzpvCum, 0, maxDzpvCum)); - histograms.h_reco_dzpvcut_pt.push_back( - ibook.book1D("num_reco_dzpvcut_pt", "#sump_{T} of reco track vs dz(PV)", nintDzpvCum, 0, maxDzpvCum)); - histograms.h_assoc2_dzpvcut_pt.push_back(ibook.book1D("num_assoc(recoToSim)_dzpvcut_pt", - "#sump_{T} of associated (recoToSim) tracks vs dz(PV)", - nintDzpvCum, - 0, - maxDzpvCum)); - histograms.h_pileup_dzpvcut_pt.push_back(ibook.book1D("num_pileup_dzpvcut_pt", - "#sump_{T} of associated (recoToSim) pileup tracks vs dz(PV)", - nintDzpvCum, - 0, - maxDzpvCum)); - histograms.h_reco_dzpvcut_pt.back()->enableSumw2(); - histograms.h_assoc2_dzpvcut_pt.back()->enableSumw2(); - histograms.h_pileup_dzpvcut_pt.back()->enableSumw2(); - histograms.h_reco_dzpvsigcut.push_back( ibook.book1D("num_reco_dzpvsigcut", "N of reco track vs dz(PV)/dzError", nintDzpvsigCum, 0, maxDzpvsigCum)); histograms.h_assoc2_dzpvsigcut.push_back(ibook.book1D("num_assoc(recoToSim)_dzpvsigcut", @@ -1561,24 +1510,6 @@ void MTVHistoProducerAlgoForTracker::bookRecoPVAssociationHistos(DQMStore::IBook nintDzpvsigCum, 0, maxDzpvsigCum)); - - histograms.h_reco_dzpvsigcut_pt.push_back(ibook.book1D( - "num_reco_dzpvsigcut_pt", "#sump_{T} of reco track vs dz(PV)/dzError", nintDzpvsigCum, 0, maxDzpvsigCum)); - histograms.h_assoc2_dzpvsigcut_pt.push_back( - ibook.book1D("num_assoc(recoToSim)_dzpvsigcut_pt", - "#sump_{T} of associated (recoToSim) tracks vs dz(PV)/dzError", - nintDzpvsigCum, - 0, - maxDzpvsigCum)); - histograms.h_pileup_dzpvsigcut_pt.push_back( - ibook.book1D("num_pileup_dzpvsigcut_pt", - "#sump_{T} of associated (recoToSim) pileup tracks vs dz(PV)/dzError", - nintDzpvsigCum, - 0, - maxDzpvsigCum)); - histograms.h_reco_dzpvsigcut_pt.back()->enableSumw2(); - histograms.h_assoc2_dzpvsigcut_pt.back()->enableSumw2(); - histograms.h_pileup_dzpvsigcut_pt.back()->enableSumw2(); } void MTVHistoProducerAlgoForTracker::bookRecodEdxHistos(DQMStore::IBooker& ibook, Histograms& histograms) { @@ -1941,8 +1872,6 @@ void MTVHistoProducerAlgoForTracker::fill_recoAssociated_simTrack_histos( histograms.h_simul_dzpvcut[count]->Fill(0); histograms.h_simul_dzpvsigcut[count]->Fill(0); - histograms.h_simul_dzpvcut_pt[count]->Fill(0, pt); - histograms.h_simul_dzpvsigcut_pt[count]->Fill(0, pt); if (isMatched) { histograms.h_assocdzpv[count]->Fill(dzPVSim); @@ -1950,14 +1879,10 @@ void MTVHistoProducerAlgoForTracker::fill_recoAssociated_simTrack_histos( histograms.h_simul2_dzpvcut[count]->Fill(0); histograms.h_simul2_dzpvsigcut[count]->Fill(0); - histograms.h_simul2_dzpvcut_pt[count]->Fill(0, pt); - histograms.h_simul2_dzpvsigcut_pt[count]->Fill(0, pt); const double dzpvcut = std::abs(track->dz(*pvPosition)); const double dzpvsigcut = dzpvcut / track->dzError(); histograms.h_assoc_dzpvcut[count]->Fill(dzpvcut); histograms.h_assoc_dzpvsigcut[count]->Fill(dzpvsigcut); - histograms.h_assoc_dzpvcut_pt[count]->Fill(dzpvcut, pt); - histograms.h_assoc_dzpvsigcut_pt[count]->Fill(dzpvsigcut, pt); } } if (simPVPosition) { @@ -2068,8 +1993,6 @@ void MTVHistoProducerAlgoForTracker::fill_generic_recoTrack_histos(const Histogr histograms.h_reco_dzpvcut[count]->Fill(std::abs(dzpv)); histograms.h_reco_dzpvsigcut[count]->Fill(std::abs(dzpvsig)); - histograms.h_reco_dzpvcut_pt[count]->Fill(std::abs(dzpv), pt); - histograms.h_reco_dzpvsigcut_pt[count]->Fill(std::abs(dzpvsig), pt); } if (simPVPosition) { histograms.h_reco_simpvz[count]->Fill(simpvz); @@ -2139,8 +2062,6 @@ void MTVHistoProducerAlgoForTracker::fill_generic_recoTrack_histos(const Histogr histograms.h_assoc2_dzpvcut[count]->Fill(std::abs(dzpv)); histograms.h_assoc2_dzpvsigcut[count]->Fill(std::abs(dzpvsig)); - histograms.h_assoc2_dzpvcut_pt[count]->Fill(std::abs(dzpv), pt); - histograms.h_assoc2_dzpvsigcut_pt[count]->Fill(std::abs(dzpvsig), pt); } if (simPVPosition) { histograms.h_assoc2_simpvz[count]->Fill(simpvz); @@ -2264,8 +2185,6 @@ void MTVHistoProducerAlgoForTracker::fill_generic_recoTrack_histos(const Histogr histograms.h_pileup_dzpvcut[count]->Fill(std::abs(dzpv)); histograms.h_pileup_dzpvsigcut[count]->Fill(std::abs(dzpvsig)); - histograms.h_pileup_dzpvcut_pt[count]->Fill(std::abs(dzpv), pt); - histograms.h_pileup_dzpvsigcut_pt[count]->Fill(std::abs(dzpvsig), pt); } if (simPVPosition) { histograms.h_pileup_simpvz[count]->Fill(simpvz);