From 24d8629645aa85ad00101e8b254f4bd3d2cf1c90 Mon Sep 17 00:00:00 2001 From: Valerio Bertacchi Date: Thu, 19 Nov 2020 04:27:59 +0100 Subject: [PATCH] DeepCore PR edit: cleaning, dedicated validation on Run3 only aestethic stuff This is the commit message #1.3: seedingDeepCore modifier This is the commit message #1.4: cleaning DeepCore plugin code This is the commit message #1.5: deepCore module by modifier This is the commit message #1.6: some small fix This is the commit message #1.7: buildfile3 fix auto parameters code format fix validation on relvant eras only, some fix fix eras bug additional plots for Run3 only move to existising object in _cff, cleaning of DeepCore cleaning --- .../python/seedingDeepCore_cff.py | 5 + .../python/upgradeWorkflowComponents.py | 19 + .../plugins/LowPtGsfElectronSeedProducer.cc | 4 +- .../PFTracking/plugins/GoodSeedProducer.cc | 6 +- .../python/trackerDrivenElectronSeeds_cfi.py | 1 + RecoTracker/CkfPattern/BuildFile.xml | 1 - .../src/CkfTrackCandidateMakerBase.cc | 10 - .../src/TransientInitialStateEstimator.cc | 4 - .../python/JetCoreRegionalStep_cff.py | 83 +- .../TkSeedGenerator/plugins/BuildFile.xml | 4 +- .../plugins/DeepCoreSeedGenerator.cc | 638 ++++++++++++++ ...eedGenerator.h => DeepCoreSeedGenerator.h} | 150 ++-- .../plugins/JetCoreDirectSeedGenerator.cc | 792 ------------------ .../plugins/JetCoreMCtruthSeedGenerator.cc | 580 +++++++++++++ .../plugins/JetCoreMCtruthSeedGenerator.h | 193 +++++ .../plugins/JetCorePerfectSeedGenerator.cc | 730 ---------------- .../plugins/JetCorePerfectSeedGenerator.h | 191 ----- .../TkSeedGenerator/plugins/SealModules.cc | 8 +- ...or_cfi.py => DeepCoreSeedGenerator_cfi.py} | 4 +- ....py => JetCoreMCtruthSeedGenerator_cfi.py} | 2 +- ...ackingParticleRecoTrackAsssociation_cff.py | 1 + ...ackingParticleRecoTrackAsssociation_cfi.py | 2 +- .../plugins/TrackAssociatorByChi2Impl.h | 6 +- .../plugins/TrackAssociatorByPositionImpl.cc | 6 +- .../plugins/TrackAssociatorByPositionImpl.h | 2 +- .../python/CkfElectronCandidateMaker_cff.py | 3 +- .../RecoTrack/plugins/MultiTrackValidator.cc | 3 +- .../plugins/TrackFromSeedProducer.cc | 23 +- ...MTVHistoProducerAlgoForTrackerBlock_cfi.py | 4 +- .../python/MultiTrackValidator_cfi.py | 2 - .../python/PostProcessorTracker_cfi.py | 14 +- .../RecoTrack/python/TrackValidation_cff.py | 78 +- .../RecoTrack/python/associators_cff.py | 5 - .../src/MTVHistoProducerAlgoForTracker.cc | 41 +- 34 files changed, 1647 insertions(+), 1968 deletions(-) create mode 100644 Configuration/ProcessModifiers/python/seedingDeepCore_cff.py create mode 100644 RecoTracker/TkSeedGenerator/plugins/DeepCoreSeedGenerator.cc rename RecoTracker/TkSeedGenerator/plugins/{JetCoreDirectSeedGenerator.h => DeepCoreSeedGenerator.h} (52%) delete mode 100644 RecoTracker/TkSeedGenerator/plugins/JetCoreDirectSeedGenerator.cc create mode 100644 RecoTracker/TkSeedGenerator/plugins/JetCoreMCtruthSeedGenerator.cc create mode 100644 RecoTracker/TkSeedGenerator/plugins/JetCoreMCtruthSeedGenerator.h delete mode 100644 RecoTracker/TkSeedGenerator/plugins/JetCorePerfectSeedGenerator.cc delete mode 100644 RecoTracker/TkSeedGenerator/plugins/JetCorePerfectSeedGenerator.h rename RecoTracker/TkSeedGenerator/python/{jetCoreDirectSeedGenerator_cfi.py => DeepCoreSeedGenerator_cfi.py} (85%) rename RecoTracker/TkSeedGenerator/python/{jetCorePerfectSeedGenerator_cfi.py => JetCoreMCtruthSeedGenerator_cfi.py} (89%) 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/upgradeWorkflowComponents.py b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py index a9876705bf631..3e2a64cda7b14 100644 --- a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py +++ b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py @@ -311,6 +311,25 @@ def condition_(self, fragment, stepList, key, hasHarvest): '--procModifiers': 'trackingMkFit' } +#DeepCore seeding for JetCore iteration workflow +class UpgradeWorkflow_seedingDeepCore(UpgradeWorkflowTracking): + def setup_(self, step, stepName, stepDict, k, properties): + if 'Reco' in step: stepDict[stepName][k] = merge([self.step3, stepDict[step][k]]) + def condition_(self, fragment, stepList, key, hasHarvest): + return '2021' in key or '2024' in key +upgradeWFs['seedingDeepCore'] = UpgradeWorkflow_seedingDeepCore( + steps = [ + 'Reco', + 'RecoGlobal', + ], + PU = [], + suffix = '_seedingDeepCore', + offset = 0.13, +) +upgradeWFs['seedingDeepCore'].step3 = { + '--procModifiers': 'seedingDeepCore' +} + # 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 a182ae1c52d63..26d8d7d9d245f 100644 --- a/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronSeedProducer.cc +++ b/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronSeedProducer.cc @@ -337,10 +337,12 @@ void LowPtGsfElectronSeedProducer::loop(const edm::Handle >& hand if (!passThrough_ && (trackRef->pt() < minPtThreshold_)) { continue; } - if(trackRef->algo() == 11) continue; //Skip jetcore tracks because the seeds are hitless // 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 b51ac64ca069d..240a4a7a86538 100644 --- a/RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.cc +++ b/RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.cc @@ -336,8 +336,6 @@ void GoodSeedProducer::produce(Event& iEvent, const EventSetup& iSetup) { for (unsigned int i = 0; i < Tk.size(); ++i) { if (useQuality_ && (!(Tk[i].quality(trackQuality_)))) continue; - if(Tk[i].algo() == 11) - continue; //Skip jetcore tracks because the seeds are hitless reco::PreId myPreId; bool GoodPreId = false; @@ -347,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/CkfPattern/BuildFile.xml b/RecoTracker/CkfPattern/BuildFile.xml index cfcd471fef04f..df7161f709460 100644 --- a/RecoTracker/CkfPattern/BuildFile.xml +++ b/RecoTracker/CkfPattern/BuildFile.xml @@ -19,5 +19,4 @@ - diff --git a/RecoTracker/CkfPattern/src/CkfTrackCandidateMakerBase.cc b/RecoTracker/CkfPattern/src/CkfTrackCandidateMakerBase.cc index 2f67a4f0be373..75217eca7add1 100644 --- a/RecoTracker/CkfPattern/src/CkfTrackCandidateMakerBase.cc +++ b/RecoTracker/CkfPattern/src/CkfTrackCandidateMakerBase.cc @@ -11,8 +11,6 @@ #include "DataFormats/TrackReco/interface/SeedStopInfo.h" #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" -// #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h" - #include "TrackingTools/PatternTools/interface/Trajectory.h" #include "TrackingTools/TrajectoryCleaning/interface/TrajectoryCleanerBySharedHits.h" @@ -28,10 +26,6 @@ #include "RecoTracker/CkfPattern/interface/CachingSeedCleanerBySharedInput.h" #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h" -// #include "RecoTracker/TransientTrackingRecHit/interface/TRecHit2DPosConstraint.h" -// #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h" -// #include "DataFormats/GeometrySurface/interface/SimpleCylinderBounds.h" -// #include "DataFormats/GeometrySurface/interface/Cylinder.h" #include "RecoTracker/Record/interface/NavigationSchoolRecord.h" #include "TrackingTools/DetLayers/interface/NavigationSchool.h" @@ -154,10 +148,6 @@ namespace cms { // set the correct navigation // NavigationSetter setter( *theNavigationSchool); - - //geometry, for jetCore iteration - // edm::ESHandle geometry_; - // es.get().get(geometry_); // propagator edm::ESHandle thePropagator; diff --git a/RecoTracker/CkfPattern/src/TransientInitialStateEstimator.cc b/RecoTracker/CkfPattern/src/TransientInitialStateEstimator.cc index b7554df80c084..39492fbf4dc94 100644 --- a/RecoTracker/CkfPattern/src/TransientInitialStateEstimator.cc +++ b/RecoTracker/CkfPattern/src/TransientInitialStateEstimator.cc @@ -75,8 +75,6 @@ std::pair TransientInitialStateEstimat TSOS startingState = measvec[actualLast].updatedState(); startingState.rescaleError(100.); - // std::cout << "DEBUG DEEPCORE: distance first hit" << measvec[0].recHit()->globalPosition().perp() << std::endl; - // avoid cloning... KFUpdator const aKFUpdator; Chi2MeasurementEstimator const aChi2MeasurementEstimator(100., 3); @@ -120,7 +118,5 @@ std::pair TransientInitialStateEstimat << "\n it's field pointer is: " << firstState.magneticField() << "\n the pointer from the state of the back fit was: " << firstMeas.updatedState().magneticField(); - // std::cout << "DEBUG deepCore: firstState------>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>:\n" << "initial parameters:" << ", inv.Pt=" << firstState.freeState()->parameters().signedInverseTransverseMomentum() << ", trans.Curv=" <transverseCurvature()<< ", p=" << firstState.freeState()->momentum().mag() << ", pt=" << firstState.freeState()->momentum().perp() <<", phi=" <momentum().phi() << ", eta="<momentum().eta() << std::endl; - // std::cout << firstState << std::endl; return std::pair(std::move(firstState), firstMeas.recHit()->det()); } diff --git a/RecoTracker/IterativeTracking/python/JetCoreRegionalStep_cff.py b/RecoTracker/IterativeTracking/python/JetCoreRegionalStep_cff.py index dab262674ac62..53be4141080d2 100644 --- a/RecoTracker/IterativeTracking/python/JetCoreRegionalStep_cff.py +++ b/RecoTracker/IterativeTracking/python/JetCoreRegionalStep_cff.py @@ -90,7 +90,13 @@ # QUALITY CUTS DURING TRACK BUILDING import TrackingTools.TrajectoryFiltering.TrajectoryFilter_cff jetCoreRegionalStepTrajectoryFilter = TrackingTools.TrajectoryFiltering.TrajectoryFilter_cff.CkfBaseTrajectoryFilter_block.clone( - minPt = 0.1, + minimumNumberOfHits = 4, + seedPairPenalty = 0, + minPt = 0.1 +) + +from Configuration.ProcessModifiers.seedingDeepCore_cff import seedingDeepCore +seedingDeepCore.toModify(jetCoreRegionalStepTrajectoryFilter, maxCCCLostHits = cms.int32(9999), maxConsecLostHits = cms.int32(2), maxLostHits = cms.int32(999), @@ -99,10 +105,8 @@ minimumNumberOfHits = cms.int32(2), pixelSeedExtension = cms.bool(False), seedExtension = cms.int32(0), - seedPairPenalty = cms.int32(0), strictSeedExtension = cms.bool(False) - -) + ) 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 @@ -122,75 +126,59 @@ jetCoreRegionalStepTrajectoryBuilder = RecoTracker.CkfPattern.GroupedCkfTrajectoryBuilder_cfi.GroupedCkfTrajectoryBuilder.clone( MeasurementTrackerName = '', trajectoryFilter = cms.PSet(refToPSet_ = cms.string('jetCoreRegionalStepTrajectoryFilter')), + #clustersToSkip = cms.InputTag('jetCoreRegionalStepClusters'), maxCand = 50, estimator = 'jetCoreRegionalStepChi2Est', maxDPhiForLooperReconstruction = cms.double(2.0), + maxPtForLooperReconstruction = cms.double(0.7) + ) +seedingDeepCore.toModify(jetCoreRegionalStepTrajectoryBuilder, maxPtForLooperReconstruction = cms.double(0), keepOriginalIfRebuildFails = True, lockHits = False, requireSeedHitsInRebuild = False, + trajectoryFilter = cms.PSet(refToPSet_ = cms.string('jetCoreRegionalStepDeepCoreTrajectoryFilter')) +) - ) - - -#customized cleaner -trajectoryCleanerBySharedHits_JetCore = cms.ESProducer("TrajectoryCleanerESProducer", - ComponentName = cms.string('jetCoreTrajectoryCleanerBySharedHits'), - ComponentType = cms.string('TrajectoryCleanerBySharedHits'), - MissingHitPenalty = cms.double(20.0), - ValidHitBonus = cms.double(5.0), - allowSharedFirstHit = cms.bool(True), +#customized cleaner for DeepCore +from TrackingTools.TrajectoryCleaning.TrajectoryCleanerBySharedHits_cfi import trajectoryCleanerBySharedHits +jetCoreRegionalStepDeepCoreTrajectoryCleaner = trajectoryCleanerBySharedHits.clone( + ComponentName = cms.string('jetCoreRegionalStepDeepCoreTrajectoryCleaner'), fractionShared = cms.double(0.45) ) -CkfBaseTrajectoryFilter_blockLoose = cms.PSet( - ComponentType = cms.string('CkfBaseTrajectoryFilter'), - chargeSignificance = cms.double(-1.0), - constantValueForLostHitsFractionFilter = cms.double(2.0), - extraNumberOfHitsBeforeTheFirstLoop = cms.int32(4), - maxCCCLostHits = cms.int32(9999), +#DeepCore filter +jetCoreRegionalStepDeepCoreTrajectoryFilter = TrackingTools.TrajectoryFiltering.TrajectoryFilter_cff.CkfBaseTrajectoryFilter_block.clone( #blockloose maxConsecLostHits = cms.int32(2), - maxLostHits = cms.int32(999), maxLostHitsFraction = cms.double(1.1), - maxNumberOfHits = cms.int32(100), - minGoodStripCharge = cms.PSet( - refToPSet_ = cms.string('SiStripClusterChargeCutNone') - ), - minHitsMinPt = cms.int32(3), - minNumberOfHitsForLoopers = cms.int32(13), - minNumberOfHitsPerLoop = cms.int32(4), - minPt = cms.double(0.9), - minimumNumberOfHits = cms.int32(2), - nSigmaMinPt = cms.double(5.0), - pixelSeedExtension = cms.bool(False), - seedExtension = cms.int32(0), - seedPairPenalty = cms.int32(0), - strictSeedExtension = cms.bool(False) + minimumNumberOfHits = cms.int32(2) ) -import RecoTracker.TkSeedGenerator.jetCoreDirectSeedGenerator_cfi -import RecoTracker.TkSeedGenerator.jetCorePerfectSeedGenerator_cfi -jetCoreRegionalStepSeeds = RecoTracker.TkSeedGenerator.jetCoreDirectSeedGenerator_cfi.jetCoreDirectSeedGenerator.clone( -# jetCoreSeeds = RecoTracker.TkSeedGenerator.jetCoreDirectSeedGenerator_cfi.jetCoreDirectSeedGenerator.clone( -# jetCoreSeeds = RecoTracker.TkSeedGenerator.jetCorePerfectSeedGenerator_cfi.JetCorePerfectSeedGenerator.clone( -vertices="firstStepPrimaryVertices" +import RecoTracker.TkSeedGenerator.DeepCoreSeedGenerator_cfi +import RecoTracker.TkSeedGenerator.JetCoreMCtruthSeedGenerator_cfi +seedingDeepCore.toReplaceWith(jetCoreRegionalStepSeeds, + RecoTracker.TkSeedGenerator.DeepCoreSeedGenerator_cfi.DeepCoreSeedGenerator.clone( + # RecoTracker.TkSeedGenerator.JetCoreMCtruthSeedGenerator_cfi.JetCoreMCtruthSeedGenerator.clone( #MCtruthSeedGenerator + vertices="firstStepPrimaryVertices" + ) ) # MAKING OF TRACK CANDIDATES import RecoTracker.CkfPattern.CkfTrackCandidates_cfi jetCoreRegionalStepTrackCandidates = RecoTracker.CkfPattern.CkfTrackCandidates_cfi.ckfTrackCandidates.clone( - # src = 'jetCoreSeeds', src = 'jetCoreRegionalStepSeeds', maxSeedsBeforeCleaning = 10000, TrajectoryBuilderPSet = cms.PSet( refToPSet_ = cms.string('jetCoreRegionalStepTrajectoryBuilder')), - TrajectoryCleaner = 'jetCoreTrajectoryCleanerBySharedHits', NavigationSchool = 'SimpleNavigationSchool', - doSeedingRegionRebuilding = True, ### these two parameters are relevant only for the CachingSeedCleanerBySharedInput #numHitsForSeedCleaner = cms.int32(50), #onlyPixelHitsForSeedCleaner = cms.bool(True), ) +seedingDeepCore.toModify(jetCoreRegionalStepTrackCandidates, + TrajectoryCleaner = 'jetCoreRegionalStepDeepCoreTrajectoryCleaner', + doSeedingRegionRebuilding = True, +) # TRACK FITTING @@ -266,22 +254,19 @@ fastSim.toModify(jetCoreRegionalStep,vertices = 'firstStepPrimaryVerticesBeforeMixing') - - # Final sequence -JetCoreRegionalStepTask = cms.Task(jetsForCoreTracking, +JetCoreRegionalStepTask = cms.Task(jetsForCoreTracking, firstStepGoodPrimaryVertices, #jetCoreRegionalStepClusters, jetCoreRegionalStepSeedLayers, jetCoreRegionalStepTrackingRegions, jetCoreRegionalStepHitDoublets, jetCoreRegionalStepSeeds, - # jetCoreSeeds, jetCoreRegionalStepTrackCandidates, jetCoreRegionalStepTracks, # jetCoreRegionalStepClassifier1,jetCoreRegionalStepClassifier2, jetCoreRegionalStep) JetCoreRegionalStep = cms.Sequence(JetCoreRegionalStepTask) -fastSim.toReplaceWith(JetCoreRegionalStepTask, +fastSim.toReplaceWith(JetCoreRegionalStepTask, cms.Task(jetCoreRegionalStepTracks, jetCoreRegionalStep)) diff --git a/RecoTracker/TkSeedGenerator/plugins/BuildFile.xml b/RecoTracker/TkSeedGenerator/plugins/BuildFile.xml index 5b8614fda81d3..a798fd1009b19 100644 --- a/RecoTracker/TkSeedGenerator/plugins/BuildFile.xml +++ b/RecoTracker/TkSeedGenerator/plugins/BuildFile.xml @@ -1,3 +1,5 @@ + + @@ -41,8 +43,6 @@ - - diff --git a/RecoTracker/TkSeedGenerator/plugins/DeepCoreSeedGenerator.cc b/RecoTracker/TkSeedGenerator/plugins/DeepCoreSeedGenerator.cc new file mode 100644 index 0000000000000..0924fecaa892c --- /dev/null +++ b/RecoTracker/TkSeedGenerator/plugins/DeepCoreSeedGenerator.cc @@ -0,0 +1,638 @@ +// -*- 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 "DeepCoreSeedGenerator.h" + +#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/SiPixelDetId/interface/PXBDetId.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 "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.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 +#include +#include "boost/multi_array.hpp" + +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" + +#include "SimDataFormats/Vertex/interface/SimVertex.h" + +#include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h" + +#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h" + +#include "TTree.h" +#include "PhysicsTools/TensorFlow/interface/TensorFlow.h" + +// +// class declaration +// + +// If the analyzer does not use TFileService, please remove +// the template argument to the base class so the class inherits +// from edm::one::EDAnalyzer<> and also remove the line from +// constructor "usesResource("TFileService");" +// This will improve performance in multithreaded jobs. + +DeepCoreSeedGenerator::DeepCoreSeedGenerator(const edm::ParameterSet& iConfig) + : + + vertices_(consumes(iConfig.getParameter("vertices"))), + pixelClusters_( + consumes>(iConfig.getParameter("pixelClusters"))), + cores_(consumes>(iConfig.getParameter("cores"))), + 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")), + nThreads(iConfig.getParameter("nThreads")), + singleThreadPool(iConfig.getParameter("singleThreadPool")), + probThr(iConfig.getParameter("probThr")) + +{ + produces(); + produces(); +} + +DeepCoreSeedGenerator::~DeepCoreSeedGenerator() {} + +void DeepCoreSeedGenerator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + auto result = std::make_unique(); + auto resultTracks = std::make_unique(); + + //-------------------TensorFlow setup - session (1/2)----------------------// + tensorflow::setLogging("3"); + graph_ = tensorflow::loadGraphDef(weightfilename_); + tensorflow::SessionOptions sessionOptions; + tensorflow::setThreading(sessionOptions, nThreads, singleThreadPool); + session_ = tensorflow::createSession(graph_, sessionOptions); + tensorflow::TensorShape input_size_eta({1, 1}); + tensorflow::TensorShape input_size_pt({1, 1}); + tensorflow::TensorShape input_size_cluster({1, jetDimX, jetDimY, Nlayer}); + //-----------------end of TF setup (1/2)----------------------// + + using namespace edm; + using namespace reco; + + iSetup.get().get(magfield_); + iSetup.get().get(geometry_); + iSetup.get().get("AnalyticalPropagator", propagator_); + + iEvent.getByToken(pixelClusters_, inputPixelClusters); + allSiPixelClusters.clear(); + siPixelDetsWithClusters.clear(); + allSiPixelClusters.reserve( + inputPixelClusters->dataSize()); // this is important, otherwise push_back invalidates the iterators + + Handle> vertices; + iEvent.getByToken(vertices_, vertices); + + Handle> cores; + iEvent.getByToken(cores_, cores); + + //--------------------------debuging lines ---------------------// + edm::ESHandle pe; + const PixelClusterParameterEstimator* pp; + iSetup.get().get(pixelCPE_, pe); + pp = pe.product(); + //--------------------------end ---------------------// + + edm::ESHandle tTopoHandle; + iSetup.get().get(tTopoHandle); + const TrackerTopology* const tTopo = tTopoHandle.product(); + + auto output = std::make_unique>(); + + int jet_number = 0; + for (unsigned int ji = 0; ji < cores->size(); ji++) { //loop jet + jet_number++; + + if ((*cores)[ji].pt() > ptMin_) { + std::set ids; + const reco::Candidate& jet = (*cores)[ji]; + const reco::Vertex& jetVertex = (*vertices)[0]; + + std::vector splitClustDirSet = splittedClusterDirections(jet, tTopo, pp, jetVertex, 1); + bool l2off = (splitClustDirSet.empty()); + if (splitClustDirSet.empty()) { //if layer 1 is broken find direcitons on layer 2 + splitClustDirSet = splittedClusterDirections(jet, tTopo, pp, jetVertex, 2); + } + splitClustDirSet.push_back(GlobalVector(jet.px(), jet.py(), jet.pz())); + for (int cc = 0; cc < (int)splitClustDirSet.size(); cc++) { + //-------------------TensorFlow setup - tensor (2/2)----------------------// + tensorflow::NamedTensorList input_tensors; + 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 TF setup (2/2)----------------------// + + GlobalVector bigClustDir = splitClustDirSet.at(cc); + + LocalPoint jetInter(0, 0, 0); + + jet_eta = jet.eta(); + jet_pt = jet.pt(); + input_tensors[0].second.matrix()(0, 0) = jet.eta(); + input_tensors[1].second.matrix()(0, 0) = jet.pt(); + + edmNew::DetSetVector::const_iterator detIt = inputPixelClusters->begin(); + + const GeomDet* globDet = + DetectorSelector(2, jet, bigClustDir, jetVertex, tTopo); //select detector mostly hitten by the jet + + if (globDet == nullptr) + continue; + + const GeomDet* goodDet1 = DetectorSelector(1, jet, bigClustDir, jetVertex, tTopo); + const GeomDet* goodDet3 = DetectorSelector(3, jet, bigClustDir, jetVertex, tTopo); + const GeomDet* goodDet4 = DetectorSelector(4, jet, bigClustDir, jetVertex, tTopo); + + for (; detIt != inputPixelClusters->end(); detIt++) { //loop deset + const edmNew::DetSet& detset = *detIt; + const GeomDet* det = + geometry_->idToDet(detset.id()); //lui sa il layer con cast a PXBDetId (vedi dentro il layer function) + + for (auto cluster = detset.begin(); cluster != detset.end(); cluster++) { //loop cluster + + const SiPixelCluster& aCluster = *cluster; + 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 cPos_local = pp->localParametersV(aCluster, (*geometry_->idToDetUnit(detIt->id())))[0].first; + + if (std::abs(cPos_local.x() - localInter.x()) / pitchX <= jetDimX / 2 && + std::abs(cPos_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); + + 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]>(0.75-o*0.1-(l2off?0.25:0))){//0.99=probThr (doesn't work the variable, SOLVE THIS ISSUE!!) + if (seedParamNN.second[i][j][o] > + (0.85 - o * 0.1 - + (l2off ? 0.35 : 0))) { //0.99=probThr (doesn't work the variable, SOLVE THIS ISSUE!!) + + 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; + double yy = xyLocal.y() + seedParamNN.first[i][j][o][1] * 0.01; + LocalPoint localSeedPoint = LocalPoint(xx, yy, 0); + + // double jet_theta = 2*std::atan(std::exp(-jet_eta)); + double track_eta = + seedParamNN.first[i][j][o][2] * 0.01 + bigClustDir.eta(); //NOT SURE ABOUT THIS 0.01, only to debug + double track_theta = 2 * std::atan(std::exp(-track_eta)); + double track_phi = + seedParamNN.first[i][j][o][3] * 0.01 + bigClustDir.phi(); //NOT SURE ABOUT THIS 0.01, only to debug + + 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; + } + if (true) { //1 TO JET CORE; 0=NO JET CORE (seeding iteration skipped, useful to total eff and FakeRate comparison) + ids.insert(seedid); + + //seed creation + float em[15] = {0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0}; //sigma**2 of the follwing parameters, LocalTrajectoryError for details + 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 + // [2]=1e-5; + // em[5]=1e-5; + // em[9]=2e-5; + // em[14]=2e-5; + long int detId = globDet->geographicalId(); + LocalTrajectoryParameters localParam(localSeedPoint, localSeedDir, TrackCharge(1)); + result->push_back(TrajectorySeed(PTrajectoryStateOnDet(localParam, pt, em, detId, /*surfaceSide*/ 0), + edm::OwnVector(), + PropagationDirection::alongMomentum)); + + GlobalPoint globalSeedPoint = globDet->surface().toGlobal(localSeedPoint); + reco::Track::CovarianceMatrix mm; + resultTracks->push_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)); + } + } + } + } + } + } //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(); + //float direction = globZdir.dot(globDetCenter); + if (direction < 0) + out = -1; + // out=1; + return out; +} + +void DeepCoreSeedGenerator::fillPixelMatrix( + const SiPixelCluster& cluster, + int layer, + Point3DBase inter, + const GeomDet* det, + tensorflow::NamedTensorList input_tensors) { //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; + input_tensors[2].second.tensor()(0, nx, ny, layer - 1) += (pix.adc) / (float)(14000); + } + } +} + +// std::pair DeepCoreSeedGenerator::SeedEvaluation( +std::pair +DeepCoreSeedGenerator::SeedEvaluation(tensorflow::NamedTensorList input_tensors) { + std::vector outputs; + std::vector output_names; + output_names.push_back(outputTensorName_[0]); + output_names.push_back(outputTensorName_[1]); + 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) { + struct trkNumCompare { + bool operator()(std::pair x, std::pair y) const { + return x.first > y.first; + } + }; + + std::set, trkNumCompare> track4detSet; + + LocalPoint jetInter(0, 0, 0); + + edmNew::DetSetVector::const_iterator detIt = inputPixelClusters->begin(); + + double minDist = 0.0; + GeomDet* output = (GeomDet*)nullptr; + + for (; detIt != inputPixelClusters->end(); detIt++) { //loop deset + + const edmNew::DetSet& detset = *detIt; + const GeomDet* det = geometry_->idToDet(detset.id()); + for (auto cluster = detset.begin(); cluster != detset.end(); cluster++) { //loop cluster + auto aClusterID = detset.id(); + if (DetId(aClusterID).subdetId() != 1) + continue; + 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; + } + } //cluster + } //detset + return output; +} +std::vector DeepCoreSeedGenerator::splittedClusterDirections(const reco::Candidate& jet, + const TrackerTopology* const tTopo, + const PixelClusterParameterEstimator* pp, + const reco::Vertex& jetVertex, + int layer) { + std::vector clustDirs; + + edmNew::DetSetVector::const_iterator detIt_int = inputPixelClusters->begin(); + + for (; detIt_int != inputPixelClusters->end(); detIt_int++) { + const edmNew::DetSet& detset_int = *detIt_int; + 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!! + + for (auto cluster = detset_int.begin(); cluster != detset_int.end(); cluster++) { + const SiPixelCluster& aCluster = *cluster; + GlobalPoint cPos = det_int->surface().toGlobal( + pp->localParametersV(aCluster, (*geometry_->idToDetUnit(detIt_int->id())))[0].first); + GlobalPoint ppv(jetVertex.position().x(), jetVertex.position().y(), jetVertex.position().z()); + GlobalVector clusterDir = cPos - ppv; + GlobalVector jetDir(jet.px(), jet.py(), jet.pz()); + // std::cout <<"deltaR" << Geom::deltaR(jetDir, clusterDir)<<", jetDir="<< jetDir << ", clusterDir=" <("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/DeepCoreSeedGenerator_TrainedModel.pb")); + desc.add>("inputTensorName", {"input_1", "input_2", "input_3"}); + desc.add>("outputTensorName", {"output_node0", "output_node1"}); + desc.add("nThreads", 1); + desc.add("singleThreadPool", "no_threads"); + desc.add("probThr", 0.99); + descriptions.addDefault(desc); +} + +//define this as a plug-in +// DEFINE_FWK_MODULE(DeepCoreSeedGenerator); diff --git a/RecoTracker/TkSeedGenerator/plugins/JetCoreDirectSeedGenerator.h b/RecoTracker/TkSeedGenerator/plugins/DeepCoreSeedGenerator.h similarity index 52% rename from RecoTracker/TkSeedGenerator/plugins/JetCoreDirectSeedGenerator.h rename to RecoTracker/TkSeedGenerator/plugins/DeepCoreSeedGenerator.h index df74a044b6a8e..e29dae8ce0b4e 100644 --- a/RecoTracker/TkSeedGenerator/plugins/JetCoreDirectSeedGenerator.h +++ b/RecoTracker/TkSeedGenerator/plugins/DeepCoreSeedGenerator.h @@ -1,18 +1,11 @@ -#ifndef RecoTracker_TkSeedGenerator_JetCoreDirectSeedGenerator_H -#define RecoTracker_TkSeedGenerator_JetCoreDirectSeedGenerator_H - -#define jetDimX 30 -#define jetDimY 30 -#define Nlayer 4 -#define Nover 3 -#define Npar 5 - +#ifndef RecoTracker_TkSeedGenerator_DeepCoreSeedGenerator_H +#define RecoTracker_TkSeedGenerator_DeepCoreSeedGenerator_H #include // user include files #include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/one/EDProducer.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" @@ -40,7 +33,6 @@ #include "DataFormats/Math/interface/Vector3D.h" #include "DataFormats/Candidate/interface/Candidate.h" - #include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h" #include "RecoLocalTracker/Records/interface/TkPixelCPERecord.h" @@ -53,8 +45,6 @@ #include "MagneticField/Engine/interface/MagneticField.h" #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" - - #include #include #include "boost/multi_array.hpp" @@ -62,92 +52,73 @@ #include "FWCore/ServiceRegistry/interface/Service.h" #include "CommonTools/UtilAlgos/interface/TFileService.h" -// #include "SimG4Core/Application/interface/G4SimTrack.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 "TTree.h" #include "PhysicsTools/TensorFlow/interface/TensorFlow.h" +namespace edm { + class Event; + class EventSetup; +} // namespace edm -namespace edm { class Event; class EventSetup; } +class DeepCoreSeedGenerator : public edm::stream::EDProducer<> { +public: + explicit DeepCoreSeedGenerator(const edm::ParameterSet&); + ~DeepCoreSeedGenerator() override; - -class JetCoreDirectSeedGenerator : public edm::one::EDProducer { - public: - explicit JetCoreDirectSeedGenerator(const edm::ParameterSet&); - ~JetCoreDirectSeedGenerator(); - - static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); // A pointer to a cluster and a list of tracks on it - struct TrackAndState - { - TrackAndState(const reco::Track *aTrack, TrajectoryStateOnSurface aState) : - track(aTrack), state(aState) {} - const reco::Track* track; + struct TrackAndState { + TrackAndState(const reco::Track* aTrack, TrajectoryStateOnSurface aState) : track(aTrack), state(aState) {} + const reco::Track* track; TrajectoryStateOnSurface state; }; - - template - struct ClusterWithTracks - { - ClusterWithTracks(const Cluster &c) : cluster(&c) {} + template + struct ClusterWithTracks { + ClusterWithTracks(const Cluster& c) : cluster(&c) {} const Cluster* cluster; std::vector tracks; }; typedef ClusterWithTracks SiPixelClusterWithTracks; - typedef boost::sub_range > SiPixelClustersWithTracks; + typedef boost::sub_range> SiPixelClustersWithTracks; - TFile* JetCoreDirectSeedGenerator_out; - TTree* JetCoreDirectSeedGeneratorTree; - // static const int jetDimX =30; - // static const int jetDimY =30; - // static const int Nlayer =4; - // static const int Nover = 3; - // static const int Npar = 4; - - // double clusterMeas[jetDimX][jetDimY][Nlayer]; double jet_pt; double jet_eta; - double pitchX = 0.01; - double pitchY = 0.015; - bool print = false; - int evt_counter =0; - - - private: - virtual void beginJob() override; - virtual void produce( edm::Event&, const edm::EventSetup&) override; - virtual void endJob() override; - - - // ----------member data --------------------------- + double pitchX = 0.01; //100 um (pixel pitch in X) + double pitchY = 0.015; //150 um (pixel pitch in Y) + static const int jetDimX = 30; //pixel dimension of NN window on layer2 + static const int jetDimY = 30; //pixel dimension of NN window on layer2 + static const int Nlayer = 4; //Number of layer used in DeepCore + static const int Nover = 3; //Max number of tracks recorded per pixel + static const 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 magfield_; edm::ESHandle geometry_; - edm::ESHandle propagator_; + edm::ESHandle propagator_; - edm::EDGetTokenT > vertices_; - edm::EDGetTokenT > pixelClusters_; + edm::EDGetTokenT> vertices_; + edm::EDGetTokenT> pixelClusters_; std::vector allSiPixelClusters; std::map siPixelDetsWithClusters; - edm::Handle< edm::DetSetVector > pixeldigisimlink; - edm::Handle > inputPixelClusters; - edm::EDGetTokenT< edm::DetSetVector > pixeldigisimlinkToken; - edm::EDGetTokenT > cores_; - // edm::EDGetTokenT > simtracksToken; - // edm::EDGetTokenT > simvertexToken; + edm::Handle> pixeldigisimlink; + edm::Handle> inputPixelClusters; + edm::EDGetTokenT> pixeldigisimlinkToken; + edm::EDGetTokenT> cores_; double ptMin_; double deltaR_; @@ -163,30 +134,39 @@ class JetCoreDirectSeedGenerator : public edm::one::EDProducer> 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> findIntersection(const GlobalVector & , const reco::Candidate::Point & ,const GeomDet*); - - void fillPixelMatrix(const SiPixelCluster &, int, auto, const GeomDet*, tensorflow::NamedTensorList); - - std::pair local2Pixel(double, double, const GeomDet*); + std::pair local2Pixel(double, double, const GeomDet*); LocalPoint pixel2Local(int, int, const GeomDet*); int pixelFlipper(const GeomDet*); - const GeomDet* DetectorSelector(int ,const reco::Candidate& jet, GlobalVector, const reco::Vertex& jetVertex, const TrackerTopology* const); - - std::vector splittedClusterDirectionsOld(const reco::Candidate&, const TrackerTopology* const, auto pp, const reco::Vertex& jetVertex ); - std::vector splittedClusterDirections(const reco::Candidate&, const TrackerTopology* const, auto pp, const reco::Vertex& jetVertex, int ); - - std::pair SeedEvaluation(tensorflow::NamedTensorList); - - - + const GeomDet* DetectorSelector( + int, const reco::Candidate&, GlobalVector, const reco::Vertex&, const TrackerTopology* const); + + std::vector splittedClusterDirectionsOld(const reco::Candidate&, + const TrackerTopology* const, + const PixelClusterParameterEstimator*, + const reco::Vertex&); //if not working,: args=2 auto + std::vector splittedClusterDirections(const reco::Candidate&, + const TrackerTopology* const, + const PixelClusterParameterEstimator*, + const reco::Vertex&, + int); //if not working,: args=2 auto + + std::pair SeedEvaluation( + tensorflow::NamedTensorList); }; #endif diff --git a/RecoTracker/TkSeedGenerator/plugins/JetCoreDirectSeedGenerator.cc b/RecoTracker/TkSeedGenerator/plugins/JetCoreDirectSeedGenerator.cc deleted file mode 100644 index 9d0270dee95b7..0000000000000 --- a/RecoTracker/TkSeedGenerator/plugins/JetCoreDirectSeedGenerator.cc +++ /dev/null @@ -1,792 +0,0 @@ -// -*- C++ -*- -// -// Package: trackJet/JetCoreDirectSeedGenerator -// Class: JetCoreDirectSeedGenerator -// -/**\class JetCoreDirectSeedGenerator JetCoreDirectSeedGenerator.cc trackJet/JetCoreDirectSeedGenerator/plugins/JetCoreDirectSeedGenerator.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 - -#define jetDimX 30 -#define jetDimY 30 -#define Nlayer 4 -#define Nover 3 -#define Npar 5 - -#include "JetCoreDirectSeedGenerator.h" - -#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/SiPixelDetId/interface/PXBDetId.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 "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.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 -#include -#include "boost/multi_array.hpp" - -#include "FWCore/ServiceRegistry/interface/Service.h" -#include "CommonTools/UtilAlgos/interface/TFileService.h" - -// #include "SimG4Core/Application/interface/G4SimTrack.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 "TTree.h" -#include "PhysicsTools/TensorFlow/interface/TensorFlow.h" - - -// -// class declaration -// - -// If the analyzer does not use TFileService, please remove -// the template argument to the base class so the class inherits -// from edm::one::EDAnalyzer<> and also remove the line from -// constructor "usesResource("TFileService");" -// This will improve performance in multithreaded jobs. - - - -JetCoreDirectSeedGenerator::JetCoreDirectSeedGenerator(const edm::ParameterSet& iConfig) : - - vertices_(consumes(iConfig.getParameter("vertices"))), - pixelClusters_(consumes >(iConfig.getParameter("pixelClusters"))), - cores_(consumes >(iConfig.getParameter("cores"))), - 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")), - nThreads(iConfig.getParameter("nThreads")), - singleThreadPool(iConfig.getParameter("singleThreadPool")), - probThr(iConfig.getParameter("probThr")) - -{ - produces(); - produces(); - - - - // edm::Service fileService; - // - // JetCoreDirectSeedGeneratorTree= fileService->make("JetCoreDirectSeedGeneratorTree","JetCoreDirectSeedGeneratorTree"); - // JetCoreDirectSeedGeneratorTree->Branch("cluster_measured",clusterMeas,"cluster_measured[30][30][4]/D"); - // JetCoreDirectSeedGeneratorTree->Branch("jet_eta",&jet_eta); - // JetCoreDirectSeedGeneratorTree->Branch("jet_pt",&jet_pt); - - - // for(int i=0; i(); - auto resultTracks = std::make_unique(); - - //-------------------TensorFlow setup - session (1/2)----------------------// - - tensorflow::setLogging("3"); - graph_=tensorflow::loadGraphDef(weightfilename_); - // output_names_=iConfig.getParameter>("outputNames"); - // for(const auto & s : iConfig.getParameter>("outputFormulas")) { output_formulas_.push_back(StringObjectFunction>(s));} - tensorflow::SessionOptions sessionOptions; - tensorflow::setThreading(sessionOptions, nThreads, singleThreadPool); - session_ = tensorflow::createSession(graph_, sessionOptions); - tensorflow::TensorShape input_size_eta({1,1}) ; - tensorflow::TensorShape input_size_pt({1,1}) ; - tensorflow::TensorShape input_size_cluster({1,jetDimX,jetDimY,Nlayer}); - - // std::cout << "input_size_cluster=" << input_size_cluster.num_elements() << "," << "," << input_size_cluster.dims() << "," << input_size_cluster.dim_size(0) << "," << input_size_cluster.dim_size(1) <<"," << input_size_cluster.dim_size(2) <<"," << input_size_cluster.dim_size(3) << std::endl; - - // input_size_cluster.set_dim(0,1); - // input_size_cluster.set_dim(1,jetDimX); - // input_size_cluster.set_dim(2,jetDimY); - // input_size_cluster.set_dim(3,Nlayer); - ; - - // tensorflow::TensorShape input_size_cluster {1,1,1,1} ; - //-----------------end of TF setup (1/2)----------------------// - - evt_counter++; -// std::cout << "event number (iterative)=" << evt_counter<< ", event number (id)="<< iEvent.id().event() << std::endl; - - - using namespace edm; - using namespace reco; - - - iSetup.get().get( magfield_ ); - iSetup.get().get(geometry_); - iSetup.get().get( "AnalyticalPropagator", propagator_ ); - - iEvent.getByToken(pixelClusters_, inputPixelClusters); - allSiPixelClusters.clear(); siPixelDetsWithClusters.clear(); - allSiPixelClusters.reserve(inputPixelClusters->dataSize()); // this is important, otherwise push_back invalidates the iterators - - // edm::Handle > simtracks; - // iEvent.getByToken(simtracksToken, simtracks); - // edm::Handle > simvertex; - // iEvent.getByToken(simvertexToken, simvertex); - - Handle > vertices; - iEvent.getByToken(vertices_, vertices); - - Handle > cores; - iEvent.getByToken(cores_, cores); - - // iEvent.getByToken(pixeldigisimlinkToken, pixeldigisimlink); - - //--------------------------debuging lines ---------------------// - edm::ESHandle pe; - const PixelClusterParameterEstimator* pp; - iSetup.get().get(pixelCPE_, pe); - pp = pe.product(); - //--------------------------end ---------------------// - - edm::ESHandle tTopoHandle; - iSetup.get().get(tTopoHandle); - const TrackerTopology* const tTopo = tTopoHandle.product(); - - auto output = std::make_unique>(); - -print = false; -int jet_number = 0; - for (unsigned int ji = 0; ji < cores->size(); ji++) { //loop jet - jet_number++; - - if ((*cores)[ji].pt() > ptMin_) { -// std::cout << "|____________________NEW JET_______________________________| jet number=" << jet_number << " " << (*cores)[ji].pt() << " " << (*cores)[ji].eta() << " " << (*cores)[ji].phi() << std::endl; - - std::set ids; - const reco::Candidate& jet = (*cores)[ji]; - const reco::Vertex& jetVertex = (*vertices)[0]; - - std::vector splitClustDirSet = splittedClusterDirections(jet, tTopo, pp, jetVertex, 1); - //std::vector splitClustDirSet = splittedClusterDirections(jet, tTopo, pp, jetVertex); - bool l2off=(splitClustDirSet.size()==0); - if(splitClustDirSet.size()==0) {//if layer 1 is broken find direcitons on layer 2 - splitClustDirSet = splittedClusterDirections(jet, tTopo, pp, jetVertex, 2); - // std::cout << "split on lay2, in numero=" << splitClustDirSet.size() << "+jetDir" << std::endl; - } - splitClustDirSet.push_back(GlobalVector(jet.px(),jet.py(),jet.pz())); - // std::cout << "splitted cluster number=" << splitClustDirSet.size() << std::endl;; - for(int cc=0; cc<(int)splitClustDirSet.size(); cc++){ - - //-------------------TensorFlow setup - tensor (2/2)----------------------// - tensorflow::NamedTensorList input_tensors; - 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()(0,x,y,l) = 0.0; - } - } - } - // auto input_matrix_eta = input_tensors[0].second.tensor(); - // auto input_matrix_pt = input_tensors[1].second.tensor(); - // auto input_matrix_cluster = input_tensors[2].second.tensor(); - - // - // std::vector inputs; - // std::vector input_names; - // - // ouput_names.push_back(inputTensorName_[0]); - // ouput_names.push_back(inputTensorName_[1]); - // ouput_names.push_back(inputTensorName_[2]); - - - - - //-----------------end of TF setup (2/2)----------------------// - - GlobalVector bigClustDir = splitClustDirSet.at(cc); - - LocalPoint jetInter(0,0,0); - - jet_eta = jet.eta(); - jet_pt = jet.pt(); - // input_tensors(0).at(0) = jet.eta(); - // input_tensors[1](0) = jet.pt(); - // input_matrix_eta(0,0) = jet.eta(); - // input_matrix_pt(0,0) = jet.pt(); - input_tensors[0].second.matrix()(0,0) = jet.eta(); - input_tensors[1].second.matrix()(0,0) = jet.pt(); - - auto jetVert = jetVertex; //trackInfo filling - - - - edmNew::DetSetVector::const_iterator detIt = inputPixelClusters->begin(); - - const GeomDet* globDet = DetectorSelector(2, jet, bigClustDir, jetVertex, tTopo); //select detector mostly hitten by the jet - - if(globDet == 0) continue; - - const GeomDet* goodDet1 = DetectorSelector(1, jet, bigClustDir, jetVertex, tTopo); - const GeomDet* goodDet3 = DetectorSelector(3, jet, bigClustDir, jetVertex, tTopo); - const GeomDet* goodDet4 = DetectorSelector(4, jet, bigClustDir, jetVertex, tTopo); - - - - for (; detIt != inputPixelClusters->end(); detIt++) { //loop deset - const edmNew::DetSet& detset = *detIt; - const GeomDet* det = geometry_->idToDet(detset.id()); //lui sa il layer con cast a PXBDetId (vedi dentro il layer function) - - for (auto cluster = detset.begin(); cluster != detset.end(); cluster++) { //loop cluster - - const SiPixelCluster& aCluster = *cluster; - 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()); - - - // GlobalPoint cPos = det->surface().toGlobal(pp->localParametersV(aCluster,(*geometry_->idToDetUnit(detIt->id())))[0].first); - LocalPoint cPos_local = pp->localParametersV(aCluster,(*geometry_->idToDetUnit(detIt->id())))[0].first; - - if(std::abs(cPos_local.x()-localInter.x())/pitchX<=jetDimX/2 && std::abs(cPos_local.y()-localInter.y())/pitchY<=jetDimY/2){ // per ora preso baricentro, da migliorare - - if(det==goodDet1 || det==goodDet3 || det==goodDet4 || det==globDet) { - // fillPixelMatrix(aCluster,lay,localInter, det, input_matrix_cluster); - fillPixelMatrix(aCluster,lay,localInter, det, input_tensors); - - } - } //cluster in ROI - } //cluster - } //detset - - // JetCoreDirectSeedGeneratorTree->Fill(); - - //HERE SOMEHOW THE NN PRODUCE THE SEED FROM THE FILLED INPUT -// std::cout << "Filling complete" << std::endl; - std::pair seedParamNN = JetCoreDirectSeedGenerator::SeedEvaluation(input_tensors); - - for(int i=0; i(0.75-o*0.1-(l2off?0.25:0))){//0.99=probThr (doesn't work the variable, SOLVE THIS ISSUE!!) - if(seedParamNN.second[i][j][o]>(0.85-o*0.1-(l2off?0.35:0))){//0.99=probThr (doesn't work the variable, SOLVE THIS ISSUE!!) - - // std::cout << "prob success=" << seedParamNN.second[i][j][o]<< ", for (x,y)=" << i <<"," <geographicalId(); - LocalTrajectoryParameters localParam(localSeedPoint, localSeedDir, TrackCharge(1)); - result->push_back(TrajectorySeed( PTrajectoryStateOnDet (localParam, pt, em, detId, /*surfaceSide*/ 0), edm::OwnVector< TrackingRecHit >() , PropagationDirection::alongMomentum)); - // LocalTrajectoryParameters localParam2(localSeedPoint, localSeedDir2, TrackCharge(1)); - // result->push_back(TrajectorySeed( PTrajectoryStateOnDet (localParam2, pt2, em, detId, /*surfaceSide*/ 0), edm::OwnVector< TrackingRecHit >() , PropagationDirection::alongMomentum)); - - GlobalPoint globalSeedPoint = globDet->surface().toGlobal(localSeedPoint); - reco::Track::CovarianceMatrix mm; - resultTracks->push_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)); - } - } - } - } - } - - - -// std::cout << "FILL!" << std::endl; - - // for(int i=0; i pt - } //jet - std::cout <<"numero di seed=" << result->size() <<", " << resultTracks->size() << std::endl; -iEvent.put(std::move(result)); -iEvent.put(std::move(resultTracks)); -} - - - - - - - - - - - - - -std::pair> JetCoreDirectSeedGenerator::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 JetCoreDirectSeedGenerator::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 JetCoreDirectSeedGenerator::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 JetCoreDirectSeedGenerator::pixelFlipper(const GeomDet* det){ - int out =1; - LocalVector locZdir(0,0,1); - GlobalVector globZdir = det->specificSurface().toGlobal(locZdir); - GlobalPoint globDetCenter = det->position(); - float direction = globZdir.x()*globDetCenter.x()+ globZdir.y()*globDetCenter.y()+ globZdir.z()*globDetCenter.z(); - //float direction = globZdir.dot(globDetCenter); - if(direction<0) out =-1; - // out=1; - return out; -} - - - -void JetCoreDirectSeedGenerator::fillPixelMatrix(const SiPixelCluster & cluster, int layer, auto inter, const GeomDet* det, tensorflow::NamedTensorList input_tensors ){//tensorflow::NamedTensorList input_tensors){ - - int flip = pixelFlipper(det); // 1=not flip, -1=flip - - for(int i=0; i 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)=0 && nx>=0 && ny>=0) { - // clusterMeas[nx][ny][layer-1] += (pix.adc)/(float)(14000);//std::cout << "clusterMeas[nx][ny][layer-1] += (pix.adc)/(float)(14000) =" << (pix.adc)/(float)(14000) << std::endl;; - // } - } - } - -} - -std::pair JetCoreDirectSeedGenerator::SeedEvaluation(tensorflow::NamedTensorList input_tensors){ - - // tensorflow::TensorShape input_size_cluster {1,jetDimX,jetDimY,Nlayer} ; - // tensorflow::TensorShape input_size_pt {1} ; - // tensorflow::TensorShape input_size_eta {1} ; - // tensorflow::NamedTensorList input_tensors; - // 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)); - - // for(int lay=0; lay()(0,nx,ny,layer-1) = clusterMeas[nx][ny][layer-1]; - // } - // } - // - // // for(size_t j =0; j < values_.size();j++) { - // // input_tensors[0].second.matrix()(0,j) = values_[j]; - // // } - - // debug!!! -/* - for(int x=0; x()(0,x,y,l)!=0){ - std::cout << "input, " << "x=" << x << ", y=" << y <<", lay=" << l << ", val =" << input_tensors[2].second.tensor()(0,x,y,l) << std::endl; - } - } - } - } //end of debug -*/ - std::vector outputs; - std::vector output_names; - output_names.push_back(outputTensorName_[0]); - output_names.push_back(outputTensorName_[1]); - tensorflow::run(session_, input_tensors, output_names, &outputs); - auto matrix_output_par = outputs.at(0).tensor(); - auto matrix_output_prob = outputs.at(1).tensor(); - - // double trackPar[jetDimX][jetDimY][Nover][Npar+1]; //NOFLAG - // double trackProb[jetDimX][jetDimY][Nover]; - - std::pair output_combined; - - - for(int x=0; x()(0,x,y,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()(0,x,y,trk,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); -// if(matrix_output_prob(0,x,y,trk,0)>0.9) std::cout << "internal output, prob= "< + +// 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/SiPixelDetId/interface/PXBDetId.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 "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.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 +#include +#include "boost/multi_array.hpp" + +#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 "TTree.h" + +// +// class declaration +// + +// If the analyzer does not use TFileService, please remove +// the template argument to the base class so the class inherits +// from edm::one::EDAnalyzer<> and also remove the line from +// constructor "usesResource("TFileService");" +// This will improve performance in multithreaded jobs. + +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() {} + +#define foreach BOOST_FOREACH + +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_); + + iEvent.getByToken(pixelClusters_, inputPixelClusters); + allSiPixelClusters.clear(); + siPixelDetsWithClusters.clear(); + allSiPixelClusters.reserve( + inputPixelClusters->dataSize()); // this is important, otherwise push_back invalidates the iterators + + edm::Handle> simtracks; + iEvent.getByToken(simtracksToken, simtracks); + edm::Handle> simvertex; + iEvent.getByToken(simvertexToken, simvertex); + + iEvent.getByToken(PSimHitToken, simhits); + + Handle> vertices; + iEvent.getByToken(vertices_, vertices); + + Handle> cores; + iEvent.getByToken(cores_, cores); + + //--------------------------debuging lines ---------------------// + edm::ESHandle pe; + const PixelClusterParameterEstimator* pp; + iSetup.get().get(pixelCPE_, pe); + pp = pe.product(); + //--------------------------end ---------------------// + + edm::ESHandle tTopoHandle; + iSetup.get().get(tTopoHandle); + const TrackerTopology* const tTopo = tTopoHandle.product(); + + auto output = std::make_unique>(); + + print = false; + int jet_number = 0; + int seed_number = 0; + + for (unsigned int ji = 0; ji < cores->size(); ji++) { //loop jet + jet_number++; + + if ((*cores)[ji].pt() > ptMin_) { + std::set ids; + const reco::Candidate& jet = (*cores)[ji]; + const reco::Vertex& jetVertex = (*vertices)[0]; + + std::vector splitClustDirSet = splittedClusterDirections(jet, tTopo, pp, jetVertex, 1); + if (splitClustDirSet.empty()) { //if layer 1 is broken find direcitons on layer 2 + splitClustDirSet = splittedClusterDirections(jet, tTopo, pp, jetVertex, 2); + } + if (inclusiveConeSeed) + splitClustDirSet.clear(); + splitClustDirSet.push_back(GlobalVector(jet.px(), jet.py(), jet.pz())); + + for (int cc = 0; cc < (int)splitClustDirSet.size(); cc++) { + GlobalVector bigClustDir = splitClustDirSet.at(cc); + + const auto& simtracksVector = simtracks.product(); + const auto& simvertexVector = simvertex.product(); + + LocalPoint jetInter(0, 0, 0); + + jet_eta = jet.eta(); + jet_pt = jet.pt(); + + const auto& jetVert = jetVertex; //trackInfo filling + + std::vector goodSimHit; + + const GeomDet* globDet = + DetectorSelector(2, jet, bigClustDir, jetVertex, tTopo); //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); + } + seed_number = goodSimTkVx.first.size(); + std::cout << "seed number in deltaR cone =" << seed_number << std::endl; + + std::vector> seedVector = + JetCoreMCtruthSeedGenerator::seedParFilling(goodSimTkVx, globDet, jet); + std::cout << "seedVector.size()=" << seedVector.size() << std::endl; + + for (uint tk = 0; tk < seedVector.size(); tk++) { + for (int pp = 0; pp < 5; pp++) { + std::cout << "seed " << tk << ", int par " << pp << "=" << seedVector.at(tk).at(pp) << std::endl; + } + LocalPoint localSeedPoint = LocalPoint(seedVector.at(tk).at(0), seedVector.at(tk).at(1), 0); + double track_theta = 2 * std::atan(std::exp(-seedVector.at(tk).at(2))); + double track_phi = seedVector.at(tk).at(3); + double pt = 1. / seedVector.at(tk).at(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.at(tk).at(2) * 400.) << 32) + (int64_t(track_phi * 400.) << 48); + if (ids.count(seedid) != 0) { + std::cout << "seed not removed with DeepCore cleaner" << std::endl; + } + ids.insert(seedid); + + //seed creation + float em[15] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + em[0] = 0.15 * 0.15; + em[2] = 0.5e-5; + em[5] = 0.5e-5; + em[9] = 2e-5; + em[14] = 2e-5; + long int detId = globDet->geographicalId(); + LocalTrajectoryParameters localParam(localSeedPoint, localSeedDir, TrackCharge(1)); + result->push_back(TrajectorySeed(PTrajectoryStateOnDet(localParam, pt, em, detId, /*surfaceSide*/ 0), + edm::OwnVector(), + PropagationDirection::alongMomentum)); + + GlobalPoint globalSeedPoint = globDet->surface().toGlobal(localSeedPoint); + reco::Track::CovarianceMatrix mm; + resultTracks->push_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)); + std::cout << "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; +} + +std::pair JetCoreMCtruthSeedGenerator::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 JetCoreMCtruthSeedGenerator::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 JetCoreMCtruthSeedGenerator::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 JetCoreMCtruthSeedGenerator::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; + + input_tensors[2].second.tensor()(0, nx, ny, layer - 1) += (pix.adc) / (float)(14000); + } + } +} + +const GeomDet* JetCoreMCtruthSeedGenerator::DetectorSelector(int llay, + const reco::Candidate& jet, + GlobalVector jetDir, + const reco::Vertex& jetVertex, + const TrackerTopology* const tTopo) { + struct trkNumCompare { + bool operator()(std::pair x, std::pair y) const { + return x.first > y.first; + } + }; + + std::set, trkNumCompare> track4detSet; + + LocalPoint jetInter(0, 0, 0); + + edmNew::DetSetVector::const_iterator detIt = inputPixelClusters->begin(); + + double minDist = 0.0; + GeomDet* output = (GeomDet*)nullptr; + + for (; detIt != inputPixelClusters->end(); detIt++) { //loop deset + + const edmNew::DetSet& detset = *detIt; + const GeomDet* det = + geometry_->idToDet(detset.id()); //lui sa il layer con cast a PXBDetId (vedi dentro il layer function) + for (auto cluster = detset.begin(); cluster != detset.end(); cluster++) { //loop cluster + auto aClusterID = detset.id(); + if (DetId(aClusterID).subdetId() != 1) + continue; + 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; + } + } //cluster + } //detset + return output; +} + +std::vector JetCoreMCtruthSeedGenerator::splittedClusterDirections( + const reco::Candidate& jet, + const TrackerTopology* const tTopo, + const PixelClusterParameterEstimator* pp, + const reco::Vertex& jetVertex, + int layer) { + std::vector clustDirs; + + edmNew::DetSetVector::const_iterator detIt_int = inputPixelClusters->begin(); + + for (; detIt_int != inputPixelClusters->end(); detIt_int++) { + const edmNew::DetSet& detset_int = *detIt_int; + 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!! + + for (auto cluster = detset_int.begin(); cluster != detset_int.end(); cluster++) { + const SiPixelCluster& aCluster = *cluster; + GlobalPoint cPos = det_int->surface().toGlobal( + pp->localParametersV(aCluster, (*geometry_->idToDetUnit(detIt_int->id())))[0].first); + GlobalPoint ppv(jetVertex.position().x(), jetVertex.position().y(), jetVertex.position().z()); + GlobalVector clusterDir = cPos - ppv; + GlobalVector jetDir(jet.px(), jet.py(), jet.pz()); + // std::cout <<"deltaR" << Geom::deltaR(jetDir, clusterDir)<<", jetDir="<< jetDir << ", clusterDir=" <> trkInterPair; + trkInterPair = findIntersection(trkMom, (reco::Candidate::Point)trkPos, globDet); + if (trkInterPair.first == false) { + GlobalVector jetDir(jet.px(), jet.py(), jet.pz()); + // double deltar = Geom::deltaR(jetDir, trkMom); + continue; + } + Basic3DVector trkInter = trkInterPair.second; + + auto localTrkInter = globDet->specificSurface().toLocal((GlobalPoint)trkInter); + + std::array tkPar{ + {localTrkInter.x(), localTrkInter.y(), st.momentum().Eta(), st.momentum().Phi(), 1 / st.momentum().Pt()}}; + output.push_back(tkPar); + + //vertex approach-------------------------------- + // auto localPos = globDet->specificSurface().toLocal((GlobalPoint)trkPos); + // std::array tkPar {{localPos.x(), localPos.y(), st.momentum().Eta(), st.momentum().Phi(), 1/st.momentum().Pt()}}; + // output.push_back(tkPar); + //end of vertex approach------------------------ + } + 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) { + //The following says we do not know what parameters are allowed so do no validation + // Please change this to state exactly what you do use, even if it is no parameters + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + +//define this as a plug-in +// DEFINE_FWK_MODULE(JetCoreMCtruthSeedGenerator); diff --git a/RecoTracker/TkSeedGenerator/plugins/JetCoreMCtruthSeedGenerator.h b/RecoTracker/TkSeedGenerator/plugins/JetCoreMCtruthSeedGenerator.h new file mode 100644 index 0000000000000..b1abe5135f172 --- /dev/null +++ b/RecoTracker/TkSeedGenerator/plugins/JetCoreMCtruthSeedGenerator.h @@ -0,0 +1,193 @@ +#ifndef RecoTracker_TkSeedGenerator_JetCoreMCtruthSeedGenerator_H +#define RecoTracker_TkSeedGenerator_JetCoreMCtruthSeedGenerator_H + +#define jetDimX 30 //pixel dimension of NN window on layer2 +#define jetDimY 30 //pixel dimension of NN window on layer2 +#define Nlayer 4 //Number of layer used in DeepCore +#define Nover 3 //Max number of tracks recorded per pixel +#define Npar 5 //Number of track parameter + +#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/SiPixelDetId/interface/PXBDetId.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 "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.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 +#include +#include "boost/multi_array.hpp" + +#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 "TTree.h" +#include "PhysicsTools/TensorFlow/interface/TensorFlow.h" + +namespace edm { + class Event; + class EventSetup; +} // namespace edm + +class JetCoreMCtruthSeedGenerator : public edm::one::EDProducer { +public: + explicit JetCoreMCtruthSeedGenerator(const edm::ParameterSet&); + ~JetCoreMCtruthSeedGenerator() override; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + // A pointer to a cluster and a list of tracks on it + struct TrackAndState { + TrackAndState(const reco::Track* aTrack, TrajectoryStateOnSurface aState) : track(aTrack), state(aState) {} + const reco::Track* track; + TrajectoryStateOnSurface state; + }; + + template + struct ClusterWithTracks { + ClusterWithTracks(const Cluster& c) : cluster(&c) {} + const Cluster* cluster; + std::vector tracks; + }; + + typedef ClusterWithTracks SiPixelClusterWithTracks; + + typedef boost::sub_range> SiPixelClustersWithTracks; + + TFile* JetCoreMCtruthSeedGenerator_out; + TTree* JetCoreMCtruthSeedGeneratorTree; + + double jet_pt; + double jet_eta; + double pitchX = 0.01; //100 um (pixel pitch in X) + double pitchY = 0.015; //150 um (pixel pitch in Y) + bool print = false; + 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_; + std::vector allSiPixelClusters; + std::map siPixelDetsWithClusters; + edm::Handle> pixeldigisimlink; + edm::Handle> inputPixelClusters; + edm::EDGetTokenT> pixeldigisimlinkToken; + 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_; + + tensorflow::GraphDef* graph_; + 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); + + std::vector splittedClusterDirections(const reco::Candidate&, + const TrackerTopology* const, + const PixelClusterParameterEstimator*, + const reco::Vertex&, + int); //if not working,: args=2 auto + + std::vector coreHitsFilling(edm::Handle>, + 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 +}; +#endif diff --git a/RecoTracker/TkSeedGenerator/plugins/JetCorePerfectSeedGenerator.cc b/RecoTracker/TkSeedGenerator/plugins/JetCorePerfectSeedGenerator.cc deleted file mode 100644 index 141b00f459e5c..0000000000000 --- a/RecoTracker/TkSeedGenerator/plugins/JetCorePerfectSeedGenerator.cc +++ /dev/null @@ -1,730 +0,0 @@ -// -*- C++ -*- -// -// Package: trackJet/JetCorePerfectSeedGenerator -// Class: JetCorePerfectSeedGenerator -// -/**\class JetCorePerfectSeedGenerator JetCorePerfectSeedGenerator.cc trackJet/JetCorePerfectSeedGenerator/plugins/JetCorePerfectSeedGenerator.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 - -#define jetDimX 30 -#define jetDimY 30 -#define Nlayer 4 -#define Nover 3 -#define Npar 5 - -#include "JetCorePerfectSeedGenerator.h" - -#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/SiPixelDetId/interface/PXBDetId.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 "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.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 -#include -#include "boost/multi_array.hpp" - -#include "FWCore/ServiceRegistry/interface/Service.h" -#include "CommonTools/UtilAlgos/interface/TFileService.h" - -// #include "SimG4Core/Application/interface/G4SimTrack.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 "TTree.h" - - -// -// class declaration -// - -// If the analyzer does not use TFileService, please remove -// the template argument to the base class so the class inherits -// from edm::one::EDAnalyzer<> and also remove the line from -// constructor "usesResource("TFileService");" -// This will improve performance in multithreaded jobs. - - - -JetCorePerfectSeedGenerator::JetCorePerfectSeedGenerator(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(); - - - - // edm::Service fileService; - // JetCorePerfectSeedGeneratorTree= fileService->make("JetCorePerfectSeedGeneratorTree","JetCorePerfectSeedGeneratorTree"); - // JetCorePerfectSeedGeneratorTree->Branch("cluster_measured",clusterMeas,"cluster_measured[30][30][4]/D"); - // JetCorePerfectSeedGeneratorTree->Branch("jet_eta",&jet_eta); - // JetCorePerfectSeedGeneratorTree->Branch("jet_pt",&jet_pt); - - // for(int i=0; i(); - auto resultTracks = std::make_unique(); - - // evt_counter++; - - - using namespace edm; - using namespace reco; - - - iSetup.get().get( magfield_ ); - iSetup.get().get(geometry_); - iSetup.get().get( "AnalyticalPropagator", propagator_ ); - - iEvent.getByToken(pixelClusters_, inputPixelClusters); - allSiPixelClusters.clear(); siPixelDetsWithClusters.clear(); - allSiPixelClusters.reserve(inputPixelClusters->dataSize()); // this is important, otherwise push_back invalidates the iterators - - edm::Handle > simtracks; - iEvent.getByToken(simtracksToken, simtracks); - edm::Handle > simvertex; - iEvent.getByToken(simvertexToken, simvertex); - - iEvent.getByToken(PSimHitToken, simhits); - - Handle > vertices; - iEvent.getByToken(vertices_, vertices); - - Handle > cores; - iEvent.getByToken(cores_, cores); - - // iEvent.getByToken(pixeldigisimlinkToken, pixeldigisimlink); - - //--------------------------debuging lines ---------------------// - edm::ESHandle pe; - const PixelClusterParameterEstimator* pp; - iSetup.get().get(pixelCPE_, pe); - pp = pe.product(); - //--------------------------end ---------------------// - - edm::ESHandle tTopoHandle; - iSetup.get().get(tTopoHandle); - const TrackerTopology* const tTopo = tTopoHandle.product(); - - auto output = std::make_unique>(); - -print = false; -int jet_number = 0; -int seed_number = 0; - - - for (unsigned int ji = 0; ji < cores->size(); ji++) { //loop jet - jet_number++; - - if ((*cores)[ji].pt() > ptMin_) { -// std::cout << "|____________________NEW JET_______________________________| jet number=" << jet_number << " " << (*cores)[ji].pt() << " " << (*cores)[ji].eta() << " " << (*cores)[ji].phi() << std::endl; - - std::set ids; - const reco::Candidate& jet = (*cores)[ji]; - const reco::Vertex& jetVertex = (*vertices)[0]; - - std::vector splitClustDirSet = splittedClusterDirections(jet, tTopo, pp, jetVertex, 1); - //std::vector splitClustDirSet = splittedClusterDirections(jet, tTopo, pp, jetVertex); - // bool l2off=(splitClustDirSet.size()==0); - if(splitClustDirSet.size()==0) {//if layer 1 is broken find direcitons on layer 2 - splitClustDirSet = splittedClusterDirections(jet, tTopo, pp, jetVertex, 2); - // std::cout << "split on lay2, in numero=" << splitClustDirSet.size() << "+jetDir" << std::endl; - } - if(inclusiveConeSeed) splitClustDirSet.clear(); - splitClustDirSet.push_back(GlobalVector(jet.px(),jet.py(),jet.pz())); - // std::cout << "splitted cluster number=" << splitClustDirSet.size() << std::endl;; - - for(int cc=0; cc<(int)splitClustDirSet.size(); cc++){ - - - GlobalVector bigClustDir = splitClustDirSet.at(cc); - - const auto & simtracksVector = simtracks.product(); - const auto & simvertexVector = simvertex.product(); - - LocalPoint jetInter(0,0,0); - - jet_eta = jet.eta(); - jet_pt = jet.pt(); - - auto jetVert = jetVertex; //trackInfo filling - - std::vector goodSimHit; - - // edmNew::DetSetVector::const_iterator detIt = inputPixelClusters->begin(); - - const GeomDet* globDet = DetectorSelector(2, jet, bigClustDir, jetVertex, tTopo); //select detector mostly hitten by the jet - - if(globDet == 0) continue; - - // const GeomDet* goodDet1 = DetectorSelector(1, jet, bigClustDir, jetVertex, tTopo); - // const GeomDet* goodDet3 = DetectorSelector(3, jet, bigClustDir, jetVertex, tTopo); - // const GeomDet* goodDet4 = DetectorSelector(4, jet, bigClustDir, jetVertex, tTopo); - - - - // for (; detIt != inputPixelClusters->end(); detIt++) { //loop deset //COMMENTATO DA QUIIIII - // const edmNew::DetSet& detset = *detIt; - // const GeomDet* det = geometry_->idToDet(detset.id()); //lui sa il layer con cast a PXBDetId (vedi dentro il layer function) - // - // for (auto cluster = detset.begin(); cluster != detset.end(); cluster++) { //loop cluster - // - // const SiPixelCluster& aCluster = *cluster; - // 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()); - // - // - // // GlobalPoint cPos = det->surface().toGlobal(pp->localParametersV(aCluster,(*geometry_->idToDetUnit(detIt->id())))[0].first); - // LocalPoint cPos_local = pp->localParametersV(aCluster,(*geometry_->idToDetUnit(detIt->id())))[0].first; - // - // if(std::abs(cPos_local.x()-localInter.x())/pitchX<=jetDimX/2 && std::abs(cPos_local.y()-localInter.y())/pitchY<=jetDimY/2){ // per ora preso baricentro, da migliorare - // - // if(det==goodDet1 || det==goodDet3 || det==goodDet4 || det==globDet) { - // // fillPixelMatrix(aCluster,lay,localInter, det, input_matrix_cluster); - // fillPixelMatrix(aCluster,lay,localInter, det, input_tensors); - // } - // } //cluster in ROI - // } //cluster - // } //detset - std::pair,std::vector> goodSimTkVx; - - if(inclusiveConeSeed) { - auto jetVert = jetVertex; - goodSimTkVx = JetCorePerfectSeedGenerator::coreTracksFillingDeltaR(simtracksVector, simvertexVector,globDet,jet,jetVert ); - } - else { - std::vector goodSimHit = JetCorePerfectSeedGenerator::coreHitsFilling(simhits,globDet,bigClustDir,jetVertex); - goodSimTkVx = JetCorePerfectSeedGenerator::coreTracksFilling(goodSimHit,simtracksVector, simvertexVector); - } - seed_number = goodSimTkVx.first.size(); - std::cout << "seed number in deltaR cone =" << seed_number << std::endl; - - std::vector> seedVector = JetCorePerfectSeedGenerator::seedParFilling(goodSimTkVx,globDet, jet); - std::cout << "seedVector.size()=" << seedVector.size()<< std::endl; - - - for(uint tk=0; tkgeographicalId(); - LocalTrajectoryParameters localParam(localSeedPoint, localSeedDir, TrackCharge(1)); - result->push_back(TrajectorySeed( PTrajectoryStateOnDet (localParam, pt, em, detId, /*surfaceSide*/ 0), edm::OwnVector< TrackingRecHit >() , PropagationDirection::alongMomentum)); - - GlobalPoint globalSeedPoint = globDet->surface().toGlobal(localSeedPoint); - reco::Track::CovarianceMatrix mm; - resultTracks->push_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)); - std::cout << "seed " << tk<< ", out, pt=" << pt << ", eta="<< globSeedDir.eta() << ", phi=" << globSeedDir.phi() < pt - } //jet -iEvent.put(std::move(result)); -iEvent.put(std::move(resultTracks)); -} - - - - - - - - - - - - - -std::pair> JetCorePerfectSeedGenerator::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 JetCorePerfectSeedGenerator::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 JetCorePerfectSeedGenerator::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 JetCorePerfectSeedGenerator::pixelFlipper(const GeomDet* det){ - int out =1; - LocalVector locZdir(0,0,1); - GlobalVector globZdir = det->specificSurface().toGlobal(locZdir); - GlobalPoint globDetCenter = det->position(); - float direction = globZdir.x()*globDetCenter.x()+ globZdir.y()*globDetCenter.y()+ globZdir.z()*globDetCenter.z(); - //float direction = globZdir.dot(globDetCenter); - if(direction<0) out =-1; - // out=1; - return out; -} - - - -void JetCorePerfectSeedGenerator::fillPixelMatrix(const SiPixelCluster & cluster, int layer, auto inter, const GeomDet* det, tensorflow::NamedTensorList input_tensors ){//tensorflow::NamedTensorList input_tensors){ - - int flip = pixelFlipper(det); // 1=not flip, -1=flip - - for(int i=0; i 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)=0 && nx>=0 && ny>=0) { - // clusterMeas[nx][ny][layer-1] += (pix.adc)/(float)(14000);//std::cout << "clusterMeas[nx][ny][layer-1] += (pix.adc)/(float)(14000) =" << (pix.adc)/(float)(14000) << std::endl;; - // } - } - } - -} - - - - -const GeomDet* JetCorePerfectSeedGenerator::DetectorSelector(int llay, const reco::Candidate& jet, GlobalVector jetDir, const reco::Vertex& jetVertex, const TrackerTopology* const tTopo){ - - struct trkNumCompare { - bool operator()(std::pair x, std::pair y) const - {return x.first > y.first;} - }; - - std::set, trkNumCompare> track4detSet; - - LocalPoint jetInter(0,0,0); - - edmNew::DetSetVector::const_iterator detIt = inputPixelClusters->begin(); - - double minDist = 0.0; - GeomDet* output = (GeomDet*)0; - - for (; detIt != inputPixelClusters->end(); detIt++) { //loop deset - - const edmNew::DetSet& detset = *detIt; - const GeomDet* det = geometry_->idToDet(detset.id()); //lui sa il layer con cast a PXBDetId (vedi dentro il layer function) - for (auto cluster = detset.begin(); cluster != detset.end(); cluster++) { //loop cluster - // const SiPixelCluster& aCluster = *cluster; - auto aClusterID= detset.id(); - if(DetId(aClusterID).subdetId()!=1) continue; - 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()),std::vector> JetCorePerfectSeedGenerator::coreTracksFilling(std::vector goodSimHit, const auto & simtracksVector, const auto & simvertexVector){ - std::vector goodSimTrk; - std::vector goodSimVtx; - - - - for(uint j=0; jsize(); j++){ - for(std::vector::const_iterator it=goodSimHit.begin(); it!=goodSimHit.end(); ++it) { - SimTrack st = simtracksVector->at(j); - if(st.trackId()==(*it).trackId()) { - // goodSimTrk.push_back(st); - - for(uint v =0; vsize(); v++) { - SimVertex sv = simvertexVector->at(v); - if((int)sv.vertexId()==(int)st.vertIndex()){ - // if(st.vertIndex()==-1) goodSimVtx.push_back((SimVertex)jVert); - //else - goodSimTrk.push_back(st); - goodSimVtx.push_back(sv); - } - } - } - } - } - std::pair,std::vector> output(goodSimTrk,goodSimVtx); - return output; -} - -std::pair,std::vector> JetCorePerfectSeedGenerator::coreTracksFillingDeltaR( const auto & simtracksVector, const auto & simvertexVector,const GeomDet* globDet, const reco::Candidate& jet, auto jetVertex){ - std::vector goodSimTrk; - std::vector goodSimVtx; - - GlobalVector jetDir(jet.px(), jet.py(), jet.pz()); - - for(uint j=0; jsize(); j++){ - SimTrack st = simtracksVector->at(j); - GlobalVector trkDir(st.momentum().Px(), st.momentum().Py(), st.momentum().Pz()); - if(Geom::deltaR(jetDir, trkDir) < deltaR_){ - if(st.charge()==0) continue; - // if((int)st.vertIndex()==-1) { - // goodSimTrk.push_back(st); - // // goodSimVtx.push_back((SimVertex)jetVertex); - // std::cout << "problema" << std::endl; - // } - // else { - for(uint v =0; vsize(); v++) { - SimVertex sv = simvertexVector->at(v); - if((int)sv.vertexId()==(int)st.vertIndex()){ - // if(st.vertIndex()==-1) goodSimVtx.push_back((SimVertex)jVert); - //else - // std::cout << "goodsimtrack " << j<< ", filling good st, pt" << st.momentum().Pt() << ", eta="<< st.momentum().Eta() << ", phi=" << st.momentum().Phi() << std::endl; - goodSimTrk.push_back(st); - goodSimVtx.push_back(sv); - } - } - // } - - } - // else std::cout << "BAD sim track " << j<< ", pt" << st.momentum().Pt() << ", eta="<< st.momentum().Eta() << ", phi=" << st.momentum().Phi() << std::endl; - } - std::pair,std::vector> output(goodSimTrk,goodSimVtx); - return output; -} - - -std::vector> JetCorePerfectSeedGenerator::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; - - std::cout << "goodSimTrk.size()" << goodSimTrk.size() << std::endl; - for(uint j=0; j> trkInterPair; - // std::cout << "sv " << (int)sv.vertexId() << "/// st " << (int)st.vertIndex()<< std::endl; - trkInterPair = findIntersection(trkMom,(reco::Candidate::Point)trkPos, globDet); - if(trkInterPair.first==false) { - GlobalVector jetDir(jet.px(), jet.py(), jet.pz()); - double deltar = Geom::deltaR(jetDir, trkMom); - // std::cout << "not intersection, deltaR=" << deltar << std::endl; - - continue; - } - Basic3DVector trkInter = trkInterPair.second; - - auto localTrkInter = globDet->specificSurface().toLocal((GlobalPoint)trkInter); - // std::cout << ", localtrackInter" << localTrkInter.x() << ", " << localTrkInter.y() << std::endl; - - - - // double tkpar[Npar]; - // // for(uint v =0; vsize(); v++) { - // // SimVertex sv = simvertexVector->at(v); - // // if((int)sv.vertexId()==(int)st.vertIndex()){ - // // - // // } - // // } - // - // tkpar[0] = localTrkInter.x(); - // tkpar[1] = localTrkInter.y(); - // tkpar[2] = st.momentum().Eta(); - // tkpar[3] = st.momentum().Phi(); - // tkpar[4] = 1/st.momentum().Pt(); - - // std::cout << "IN, pt=" < -#include -#include "boost/multi_array.hpp" - -#include "FWCore/ServiceRegistry/interface/Service.h" -#include "CommonTools/UtilAlgos/interface/TFileService.h" - -// #include "SimG4Core/Application/interface/G4SimTrack.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 "TTree.h" -#include "PhysicsTools/TensorFlow/interface/TensorFlow.h" - - -namespace edm { class Event; class EventSetup; } - - -class JetCorePerfectSeedGenerator : public edm::one::EDProducer { - public: - explicit JetCorePerfectSeedGenerator(const edm::ParameterSet&); - ~JetCorePerfectSeedGenerator(); - - static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); - // A pointer to a cluster and a list of tracks on it - struct TrackAndState - { - TrackAndState(const reco::Track *aTrack, TrajectoryStateOnSurface aState) : - track(aTrack), state(aState) {} - const reco::Track* track; - TrajectoryStateOnSurface state; - }; - - - template - struct ClusterWithTracks - { - ClusterWithTracks(const Cluster &c) : cluster(&c) {} - const Cluster* cluster; - std::vector tracks; - }; - - typedef ClusterWithTracks SiPixelClusterWithTracks; - - typedef boost::sub_range > SiPixelClustersWithTracks; - - TFile* JetCorePerfectSeedGenerator_out; - TTree* JetCorePerfectSeedGeneratorTree; - // static const int jetDimX =30; - // static const int jetDimY =30; - // static const int Nlayer =4; - // static const int Nover = 3; - // static const int Npar = 4; - - // double clusterMeas[jetDimX][jetDimY][Nlayer]; - double jet_pt; - double jet_eta; - double pitchX = 0.01; - double pitchY = 0.015; - bool print = false; - int evt_counter =0; - bool inclusiveConeSeed = true; //true= fill tracks in a cone of deltaR_, false=fill tracks which produce SimHit on globDet - - - private: - virtual void beginJob() override; - virtual void produce( edm::Event&, const edm::EventSetup&) override; - virtual void endJob() override; - - - // ----------member data --------------------------- - std::string propagatorName_; - edm::ESHandle magfield_; - edm::ESHandle geometry_; - edm::ESHandle propagator_; - - edm::EDGetTokenT > vertices_; - edm::EDGetTokenT > pixelClusters_; - std::vector allSiPixelClusters; - std::map siPixelDetsWithClusters; - edm::Handle< edm::DetSetVector > pixeldigisimlink; - edm::Handle > inputPixelClusters; - edm::EDGetTokenT< edm::DetSetVector > pixeldigisimlinkToken; - 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_; - - tensorflow::GraphDef* graph_; - tensorflow::Session* session_; - - - - std::pair> findIntersection(const GlobalVector & , const reco::Candidate::Point & ,const GeomDet*); - - void fillPixelMatrix(const SiPixelCluster &, int, auto, const GeomDet*, tensorflow::NamedTensorList); - - std::pair local2Pixel(double, double, const GeomDet*); - - LocalPoint pixel2Local(int, int, const GeomDet*); - - int pixelFlipper(const GeomDet*); - - const GeomDet* DetectorSelector(int ,const reco::Candidate& jet, GlobalVector, const reco::Vertex& jetVertex, const TrackerTopology* const); - - std::vector splittedClusterDirections(const reco::Candidate&, const TrackerTopology* const, auto pp, const reco::Vertex& jetVertex, int ); - - std::vector coreHitsFilling(auto,const GeomDet*,GlobalVector,const reco::Vertex&); - std::pair,std::vector> coreTracksFilling(std::vector, const auto &, const auto &); - - std::vector> seedParFilling(std::pair,std::vector>,const GeomDet*, const reco::Candidate&); - - std::pair,std::vector> coreTracksFillingDeltaR( const auto &, const auto &,const GeomDet* , const reco::Candidate&,auto ); - - -}; -#endif diff --git a/RecoTracker/TkSeedGenerator/plugins/SealModules.cc b/RecoTracker/TkSeedGenerator/plugins/SealModules.cc index 2263747ffba03..4a6427d461680 100644 --- a/RecoTracker/TkSeedGenerator/plugins/SealModules.cc +++ b/RecoTracker/TkSeedGenerator/plugins/SealModules.cc @@ -40,8 +40,8 @@ using SeedCreatorFromRegionConsecutiveHitsTripletOnlyEDProducer = SeedCreatorFromRegionHitsEDProducerT; DEFINE_FWK_MODULE(SeedCreatorFromRegionConsecutiveHitsTripletOnlyEDProducer); -#include "RecoTracker/TkSeedGenerator/plugins/JetCoreDirectSeedGenerator.h" -DEFINE_FWK_MODULE(JetCoreDirectSeedGenerator); +#include "RecoTracker/TkSeedGenerator/plugins/DeepCoreSeedGenerator.h" +DEFINE_FWK_MODULE(DeepCoreSeedGenerator); -#include "RecoTracker/TkSeedGenerator/plugins/JetCorePerfectSeedGenerator.h" -DEFINE_FWK_MODULE(JetCorePerfectSeedGenerator); +#include "RecoTracker/TkSeedGenerator/plugins/JetCoreMCtruthSeedGenerator.h" +DEFINE_FWK_MODULE(JetCoreMCtruthSeedGenerator); diff --git a/RecoTracker/TkSeedGenerator/python/jetCoreDirectSeedGenerator_cfi.py b/RecoTracker/TkSeedGenerator/python/DeepCoreSeedGenerator_cfi.py similarity index 85% rename from RecoTracker/TkSeedGenerator/python/jetCoreDirectSeedGenerator_cfi.py rename to RecoTracker/TkSeedGenerator/python/DeepCoreSeedGenerator_cfi.py index 48106708be6cf..8958bb75160cf 100644 --- a/RecoTracker/TkSeedGenerator/python/jetCoreDirectSeedGenerator_cfi.py +++ b/RecoTracker/TkSeedGenerator/python/DeepCoreSeedGenerator_cfi.py @@ -1,6 +1,6 @@ import FWCore.ParameterSet.Config as cms -jetCoreDirectSeedGenerator = cms.EDProducer("JetCoreDirectSeedGenerator", +DeepCoreSeedGenerator = cms.EDProducer("DeepCoreSeedGenerator", vertices= cms.InputTag("offlinePrimaryVertices"), pixelClusters= cms.InputTag("siPixelClustersPreSplitting"), cores= cms.InputTag("jetsForCoreTracking"), @@ -9,7 +9,7 @@ chargeFractionMin= cms.double(18000.0), centralMIPCharge= cms.double(2), pixelCPE= cms.string( "PixelCPEGeneric" ), - weightFile= cms.FileInPath("RecoTracker/TkSeedGenerator/data/JetCoreDirectSeedGenerator_TrainedModel.pb"), + weightFile= cms.FileInPath("RecoTracker/TkSeedGenerator/data/DeepCoreSeedGenerator_TrainedModel.pb"), inputTensorName= cms.vstring(["input_1","input_2","input_3"]), outputTensorName= cms.vstring(["output_node0","output_node1"]), nThreads= cms.uint32(1), diff --git a/RecoTracker/TkSeedGenerator/python/jetCorePerfectSeedGenerator_cfi.py b/RecoTracker/TkSeedGenerator/python/JetCoreMCtruthSeedGenerator_cfi.py similarity index 89% rename from RecoTracker/TkSeedGenerator/python/jetCorePerfectSeedGenerator_cfi.py rename to RecoTracker/TkSeedGenerator/python/JetCoreMCtruthSeedGenerator_cfi.py index 4ade6668c734b..2b517bbc2c297 100644 --- a/RecoTracker/TkSeedGenerator/python/jetCorePerfectSeedGenerator_cfi.py +++ b/RecoTracker/TkSeedGenerator/python/JetCoreMCtruthSeedGenerator_cfi.py @@ -1,6 +1,6 @@ import FWCore.ParameterSet.Config as cms -JetCorePerfectSeedGenerator = cms.EDProducer("JetCorePerfectSeedGenerator", +JetCoreMCtruthSeedGenerator = cms.EDProducer("JetCoreMCtruthSeedGenerator", vertices= cms.InputTag("offlinePrimaryVertices"), pixelClusters= cms.InputTag("siPixelClustersPreSplitting"), cores= cms.InputTag("jetsForCoreTracking"), diff --git a/SimTracker/TrackAssociation/python/trackingParticleRecoTrackAsssociation_cff.py b/SimTracker/TrackAssociation/python/trackingParticleRecoTrackAsssociation_cff.py index 1d51c2ad6a7b8..6e16b746a54b7 100644 --- a/SimTracker/TrackAssociation/python/trackingParticleRecoTrackAsssociation_cff.py +++ b/SimTracker/TrackAssociation/python/trackingParticleRecoTrackAsssociation_cff.py @@ -14,3 +14,4 @@ assoc2GsfTracks.label_tr = 'electronGsfTracks' assocOutInConversionTracks.label_tr = 'ckfOutInTracksFromConversions' assocInOutConversionTracks.label_tr = 'ckfInOutTracksFromConversions' + diff --git a/SimTracker/TrackAssociation/python/trackingParticleRecoTrackAsssociation_cfi.py b/SimTracker/TrackAssociation/python/trackingParticleRecoTrackAsssociation_cfi.py index 47592623ea727..cdc07e07c9fd8 100644 --- a/SimTracker/TrackAssociation/python/trackingParticleRecoTrackAsssociation_cfi.py +++ b/SimTracker/TrackAssociation/python/trackingParticleRecoTrackAsssociation_cfi.py @@ -2,7 +2,6 @@ trackingParticleRecoTrackAsssociation = cms.EDProducer("TrackAssociatorEDProducer", associator = cms.InputTag('quickTrackAssociatorByHits'), - # associator = cms.InputTag('trackAssociatorByChi2'), label_tp = cms.InputTag("mix","MergedTrackTruth"), label_tr = cms.InputTag("generalTracks"), ignoremissingtrackcollection=cms.untracked.bool(False) @@ -10,3 +9,4 @@ from Configuration.ProcessModifiers.premix_stage2_cff import premix_stage2 premix_stage2.toModify(trackingParticleRecoTrackAsssociation, label_tp = "mixData:MergedTrackTruth") + diff --git a/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByChi2Impl.h b/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByChi2Impl.h index 224542d186d78..448c9393bdd47 100644 --- a/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByChi2Impl.h +++ b/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByChi2Impl.h @@ -2,7 +2,7 @@ #define TrackAssociatorByChi2Impl_h /** \class TrackAssociatorByChi2Impl - * Class that performs the association of reco::Tracks and TrackingParticles evaluating the chi2 of reco tracks parameters and sim tracks parameters. The cut can be tuned from the config file: see data/TrackAssociatorByChi2.cfi. Note that the Association Map is filled with -ch2 and not chi2 because it is ordered using std::greater: the track with the lowest association chi2 will be the first in the output map.It is possible to use only diagonal terms (associator by pulls) seeting onlyDiagonal = true in the PSet + * Class that performs the association of reco::Tracks and TrackingParticles evaluating the chi2 of reco tracks parameters and sim tracks parameters. The cut can be tuned from the config file: see data/TrackAssociatorByChi2.cfi. Note that the Association Map is filled with -ch2 and not chi2 because it is ordered using std::greater: the track with the lowest association chi2 will be the first in the output map.It is possible to use only diagonal terms (associator by pulls) seeting onlyDiagonal = true in the PSet * * \author cerati, magni */ @@ -47,10 +47,10 @@ class TrackAssociatorByChi2Impl : public reco::TrackToTrackingParticleAssociator chi2cut(conf.getParameter("chi2cut")), onlyDiagonal(conf.getParameter("onlyDiagonal")), bsSrc(conf.getParameter("beamSpot")) { - theMF=mF; + theMF=mF; if (onlyDiagonal) edm::LogInfo("TrackAssociator") << " ---- Using Off Diagonal Covariance Terms = 0 ---- " << "\n"; - else + else edm::LogInfo("TrackAssociator") << " ---- Using Off Diagonal Covariance Terms != 0 ---- " << "\n"; } */ diff --git a/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionImpl.cc b/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionImpl.cc index 13edfbbe4b6e9..ec74aab8ca3db 100644 --- a/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionImpl.cc +++ b/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionImpl.cc @@ -116,7 +116,8 @@ double TrackAssociatorByPositionImpl::quality(const TrajectoryStateOnSurface& tr RecoToSimCollection TrackAssociatorByPositionImpl::associateRecoToSim( const edm::RefToBaseVector& tCH, const edm::RefVector& tPCH) const { - RecoToSimCollection outputCollection(productGetter_);; + RecoToSimCollection outputCollection(productGetter_); + ; //for each reco track find a matching tracking particle std::pair minPair; const double dQmin_default = 1542543; @@ -168,7 +169,8 @@ RecoToSimCollection TrackAssociatorByPositionImpl::associateRecoToSim( SimToRecoCollection TrackAssociatorByPositionImpl::associateSimToReco( const edm::RefToBaseVector& tCH, const edm::RefVector& tPCH) const { - SimToRecoCollection outputCollection(productGetter_);; + 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 8b3a804b214c5..be57595be36a4 100644 --- a/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionImpl.h +++ b/SimTracker/TrackAssociatorProducers/plugins/TrackAssociatorByPositionImpl.h @@ -70,7 +70,7 @@ class TrackAssociatorByPositionImpl : public reco::TrackToTrackingParticleAssoci private: double quality(const TrajectoryStateOnSurface&, const TrajectoryStateOnSurface&) const; - + edm::EDProductGetter const* productGetter_; const TrackingGeometry* theGeometry; const Propagator* thePropagator; diff --git a/TrackingTools/GsfTracking/python/CkfElectronCandidateMaker_cff.py b/TrackingTools/GsfTracking/python/CkfElectronCandidateMaker_cff.py index 1f47cce012bb5..b3e6be772535b 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) ) # CKFTrackCandidateMaker from RecoTracker.CkfPattern.CkfTrackCandidates_cff import * diff --git a/Validation/RecoTrack/plugins/MultiTrackValidator.cc b/Validation/RecoTrack/plugins/MultiTrackValidator.cc index c2cc78c99bc18..9d1e404908fee 100644 --- a/Validation/RecoTrack/plugins/MultiTrackValidator.cc +++ b/Validation/RecoTrack/plugins/MultiTrackValidator.cc @@ -314,7 +314,6 @@ void MultiTrackValidator::bookHistograms(DQMStore::IBooker& ibook, ibook.setCurrentFolder(dirName); const bool doResolutionPlots = doResolutionPlots_[www]; - // const bool doResolutionPlots = true; if (doSimTrackPlots_) { histoProducerAlgo_->bookSimTrackHistos(ibook, histograms.histoProducerAlgo, doResolutionPlots); @@ -1163,7 +1162,7 @@ void MultiTrackValidator::dqmAnalyze(const edm::Event& event, TrackingParticle::Vector momentumTP = parametersDefinerTP->momentum(event, setup, tpr); TrackingParticle::Point vertexTP = parametersDefinerTP->vertex(event, setup, tpr); int chargeTP = tpr->charge(); - // if (dirName_=="Tracking/jetCoreRegionalStep/") std::cout << "DEBUG direct JETCORE iteration ----------------------------" << std::endl; + histoProducerAlgo_->fill_ResoAndPull_recoTrack_histos( histograms.histoProducerAlgo, w, momentumTP, vertexTP, chargeTP, *track, bs.position()); } diff --git a/Validation/RecoTrack/plugins/TrackFromSeedProducer.cc b/Validation/RecoTrack/plugins/TrackFromSeedProducer.cc index 98dab2edeff39..ed3745a2644a1 100644 --- a/Validation/RecoTrack/plugins/TrackFromSeedProducer.cc +++ b/Validation/RecoTrack/plugins/TrackFromSeedProducer.cc @@ -128,7 +128,7 @@ void TrackFromSeedProducer::produce(edm::StreamID, edm::Event& iEvent, const edm edm::ESHandle httopo; iSetup.get().get(httopo); const TrackerTopology& ttopo = *httopo; - + edm::ESHandle geometry_; iSetup.get().get(geometry_); @@ -137,17 +137,13 @@ void TrackFromSeedProducer::produce(edm::StreamID, edm::Event& iEvent, const edm 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) { //deepCore seeds (jetCoreDirectSeedGenerator) - // std::cout << "DEBUG: 0 hit seed " << std::endl; - const Surface *deepCore_sruface = &geometry_->idToDet(seed.startingState().detId())->specificSurface(); - state = trajectoryStateTransform::transientState( seed.startingState(), deepCore_sruface, theMF.product()); - } - else { + 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()); + state = trajectoryStateTransform::transientState(seed.startingState(), lastRecHit->surface(), theMF.product()); } TrajectoryStateClosestToBeamLine tsAtClosestApproachSeed = tscblBuilder(*state.freeState(), *beamSpot); //as in TrackProducerAlgorithm @@ -162,13 +158,6 @@ void TrackFromSeedProducer::produce(edm::StreamID, edm::Event& iEvent, const edm PerigeeTrajectoryError seedPerigeeErrors = PerigeeConversions::ftsToPerigeeError(tsAtClosestApproachSeed.trackStateAtPCA()); tracks->emplace_back(0., 0., vSeed1, pSeed, state.charge(), seedPerigeeErrors.covarianceMatrix()); - // std::cout << "DEBUG: SEED VALIDATOR PASSED ------------" << std::endl; - // std::cout << "initial parameters:" << ", inv.Pt=" << state.freeState()->parameters().signedInverseTransverseMomentum() << ", trans.Curv=" <transverseCurvature()<< ", p=" << state.freeState()->momentum().mag() << ", pt=" << state.freeState()->momentum().perp() <<", phi=" <momentum().phi() << ", eta="<momentum().eta() << std::endl; - // std::cout << "initial matrix (diag)=" << std::sqrt(state.freeState()->curvilinearError().matrix()(0, 0)) << " , " << std::sqrt(state.freeState()->curvilinearError().matrix()(1, 1)) << " , " << std::sqrt(state.freeState()->curvilinearError().matrix()(2, 2)) << " , " << std::sqrt(state.freeState()->curvilinearError().matrix()(3, 3)) << " , " << std::sqrt(state.freeState()->curvilinearError().matrix()(4, 4)) << std::endl; - // std::cout << "initial matrix (diag)=" << std::sqrt(state.localError().matrix()(0, 0)) << " , " << std::sqrt(state.localError().matrix()(1, 1)) << " , " << std::sqrt(state.localError().matrix()(2, 2)) << " , " << std::sqrt(state.localError().matrix()(3, 3)) << " , " << std::sqrt(state.localError().matrix()(4, 4)) << std::endl; - // std::cout << "PCA parameters:" << ", inv.Pt=" << tsAtClosestApproachSeed.trackStateAtPCA().parameters().signedInverseTransverseMomentum() << ", trans.Curv=" <1 TeV - ) from Configuration.Eras.Modifier_fastSim_cff import fastSim diff --git a/Validation/RecoTrack/python/PostProcessorTracker_cfi.py b/Validation/RecoTrack/python/PostProcessorTracker_cfi.py index 73c809914c4d7..4444651e54a0e 100644 --- a/Validation/RecoTrack/python/PostProcessorTracker_cfi.py +++ b/Validation/RecoTrack/python/PostProcessorTracker_cfi.py @@ -16,7 +16,7 @@ def _addNoFlow(module): if not tmp[ind-1] in _noflowSeen: module.noFlowDists.append(tmp[ind-1]) -_defaultSubdirs = ["Tracking/Track/*", "Tracking/TrackTPPtLess09/*", "Tracking/TrackFromPV/*", "Tracking/TrackFromPVAllTP/*", "Tracking/TrackAllTPEffic/*", "Tracking/TrackBuilding/*","Tracking/TrackConversion/*", "Tracking/TrackGsf/*", "Tracking/jetCoreRegionalStep/*"] +_defaultSubdirs = ["Tracking/Track/*", "Tracking/TrackTPPtLess09/*", "Tracking/TrackFromPV/*", "Tracking/TrackFromPVAllTP/*", "Tracking/TrackAllTPEffic/*", "Tracking/TrackBuilding/*","Tracking/TrackConversion/*", "Tracking/TrackGsf/*"] _defaultSubdirsSummary = [e.replace("/*","") for e in _defaultSubdirs] postProcessorTrack = DQMEDHarvester("DQMGenericClient", @@ -317,6 +317,18 @@ def _addNoFlow(module): postProcessorTrackSummary ) +from Configuration.Eras.Modifier_run3_common_cff import run3_common +postProcessorTrackRun3 = postProcessorTrack.clone() +postProcessorTrackRun3.subDirs.extend(["Tracking/JetCore/*"]) +run3_common.toReplaceWith(postProcessorTrack,postProcessorTrackRun3) +postProcessorTrackSummaryRun3 = postProcessorTrackSummary.clone() +postProcessorTrackSummaryRun3.subDirs.extend(["Tracking/JetCore/*"]) +run3_common.toReplaceWith(postProcessorTrackSummary,postProcessorTrackSummaryRun3) +postProcessorTrack2DRun3 = postProcessorTrack2D.clone() +postProcessorTrack2DRun3.subDirs.extend(["Tracking/JetCore/*"]) +run3_common.toReplaceWith(postProcessorTrack2D,postProcessorTrack2DRun3) + + 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 7d3c576a1fea2..c55d42ac2244e 100644 --- a/Validation/RecoTrack/python/TrackValidation_cff.py +++ b/Validation/RecoTrack/python/TrackValidation_cff.py @@ -2,7 +2,6 @@ import FWCore.ParameterSet.Config as cms from SimTracker.TrackAssociatorProducers.trackAssociatorByChi2_cfi import * -from SimTracker.TrackAssociatorProducers.trackAssociatorByPosition_cfi import * from SimTracker.TrackAssociatorProducers.quickTrackAssociatorByHits_cfi import * from SimTracker.TrackAssociation.trackingParticleRecoTrackAsssociation_cfi import * import Validation.RecoTrack.MultiTrackValidator_cfi @@ -354,32 +353,17 @@ 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")) #perfectSeeding used in Connecting the Dots -highPtJets = cms.EDFilter("CandPtrSelector", src = cms.InputTag("ak4CaloJets"), cut = cms.string("pt()>1000 && eta()<1.4 && eta()>-1.4")) -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() -# MTVtrackAssociatorByChi2 = SimTracker.TrackAssociatorProducers.trackAssociatorByChi2_cfi.trackAssociatorByChi2.clone() #maybe modification needed - -MTVTrackAssociationByChi2 = trackingParticleRecoTrackAsssociation.clone( - associator = cms.InputTag('trackAssociatorByChi2'), - # UseAssociators = cms.bool(True) -) - - ## MTV instances 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) @@ -550,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"], @@ -617,16 +599,24 @@ def _uniqueFirstLayers(layerList): _setForEra(trackValidatorGsfTracks.histoProducerAlgoBlock, _eraName, _era, seedingLayerSets=trackValidator.histoProducerAlgoBlock.seedingLayerSets.value()+locals()["_seedingLayerSetsForElectrons"+_postfix]) # For jetCore tracks -trackValidatorJetCore = trackValidator.clone( - dirName = "Tracking/jetCoreRegionalStep/", +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, - label = ['jetCoreRegionalStepTracks'], - doResolutionPlotsForLabels = ['jetCoreRegionalStepTracks'], ) +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( @@ -678,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, @@ -719,13 +708,18 @@ def _uniqueFirstLayers(layerList): tracksPreValidation ) +from Configuration.Eras.Modifier_run3_common_cff import run3_common +run3_common.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) @@ -859,8 +853,6 @@ def _uniqueFirstLayers(layerList): trackValidatorSeedingTrackingOnly = _trackValidatorSeedingBuilding.clone( dirName = "Tracking/TrackSeeding/", - associators = ["trackAssociatorByChi2"], - UseAssociators = True, label = _seedSelectors, doSeedPlots = True, doResolutionPlotsForLabels = [ "seedTracksjetCoreRegionalStepSeeds",] @@ -871,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: @@ -900,17 +907,26 @@ def _uniqueFirstLayers(layerList): trackValidatorsTrackingOnly += trackValidatorSeedingPreSplittingTrackingOnly trackValidatorsTrackingOnly += trackValidatorBuilding trackValidatorsTrackingOnly += trackValidatorBuildingPreSplitting -trackValidatorsTrackingOnly += trackValidatorJetCore trackValidatorsTrackingOnly.replace(trackValidatorConversionStandalone, trackValidatorConversionTrackingOnly) trackValidatorsTrackingOnly.remove(trackValidatorGsfTracksStandalone) trackValidatorsTrackingOnly.replace(trackValidatorBHadronStandalone, trackValidatorBHadronTrackingOnly) + +run3_common.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/associators_cff.py b/Validation/RecoTrack/python/associators_cff.py index 6bfd14512b112..ab148f9c9e018 100644 --- a/Validation/RecoTrack/python/associators_cff.py +++ b/Validation/RecoTrack/python/associators_cff.py @@ -19,11 +19,6 @@ hltTrackAssociatorByHits.UseSplitting = cms.bool( False ) hltTrackAssociatorByHits.ThreeHitTracksAreSpecial = cms.bool( False ) -#NOT BY HIT NOW!!! uncomment lines 16-20 to do byhits -# hltTrackAssociatorByHits = SimTracker.TrackAssociatorProducers.trackAssociatorByChi2_cfi.trackAssociatorByChi2.clone() - - - hltTrackAssociatorByDeltaR = SimTracker.TrackAssociatorProducers.trackAssociatorByPosition_cfi.trackAssociatorByPosition.clone() hltTrackAssociatorByDeltaR.method = cms.string('momdr') hltTrackAssociatorByDeltaR.QCut = cms.double(0.5) diff --git a/Validation/RecoTrack/src/MTVHistoProducerAlgoForTracker.cc b/Validation/RecoTrack/src/MTVHistoProducerAlgoForTracker.cc index 22a560bd0ba74..8664c23ed5abb 100644 --- a/Validation/RecoTrack/src/MTVHistoProducerAlgoForTracker.cc +++ b/Validation/RecoTrack/src/MTVHistoProducerAlgoForTracker.cc @@ -237,25 +237,25 @@ MTVHistoProducerAlgoForTracker::MTVHistoProducerAlgoForTracker(const edm::Parame nintMVA = pset.getParameter("nintMVA"); //parameters for resolution plots - ptRes_rangeMin = 10*pset.getParameter("ptRes_rangeMin"); - ptRes_rangeMax = 10*pset.getParameter("ptRes_rangeMax"); - ptRes_nbin = 10*pset.getParameter("ptRes_nbin"); + ptRes_rangeMin = pset.getParameter("ptRes_rangeMin"); + ptRes_rangeMax = pset.getParameter("ptRes_rangeMax"); + ptRes_nbin = pset.getParameter("ptRes_nbin"); - phiRes_rangeMin = 10*pset.getParameter("phiRes_rangeMin"); - phiRes_rangeMax = 10*pset.getParameter("phiRes_rangeMax"); - phiRes_nbin = 10*pset.getParameter("phiRes_nbin"); + phiRes_rangeMin = pset.getParameter("phiRes_rangeMin"); + phiRes_rangeMax = pset.getParameter("phiRes_rangeMax"); + phiRes_nbin = pset.getParameter("phiRes_nbin"); - cotThetaRes_rangeMin = 10*pset.getParameter("cotThetaRes_rangeMin"); - cotThetaRes_rangeMax = 10*pset.getParameter("cotThetaRes_rangeMax"); - cotThetaRes_nbin = 10*pset.getParameter("cotThetaRes_nbin"); + cotThetaRes_rangeMin = pset.getParameter("cotThetaRes_rangeMin"); + cotThetaRes_rangeMax = pset.getParameter("cotThetaRes_rangeMax"); + cotThetaRes_nbin = pset.getParameter("cotThetaRes_nbin"); - dxyRes_rangeMin = 10*pset.getParameter("dxyRes_rangeMin"); - dxyRes_rangeMax = 10*pset.getParameter("dxyRes_rangeMax"); - dxyRes_nbin = 10*pset.getParameter("dxyRes_nbin"); + dxyRes_rangeMin = pset.getParameter("dxyRes_rangeMin"); + dxyRes_rangeMax = pset.getParameter("dxyRes_rangeMax"); + dxyRes_nbin = pset.getParameter("dxyRes_nbin"); - dzRes_rangeMin = 10*pset.getParameter("dzRes_rangeMin"); - dzRes_rangeMax = 10*pset.getParameter("dzRes_rangeMax"); - dzRes_nbin = 10*pset.getParameter("dzRes_nbin"); + dzRes_rangeMin = pset.getParameter("dzRes_rangeMin"); + dzRes_rangeMax = pset.getParameter("dzRes_rangeMax"); + dzRes_nbin = pset.getParameter("dzRes_nbin"); maxDzpvCum = pset.getParameter("maxDzpvCumulative"); nintDzpvCum = pset.getParameter("nintDzpvCumulative"); @@ -2435,17 +2435,6 @@ void MTVHistoProducerAlgoForTracker::fill_ResoAndPull_recoTrack_histos(const His double phiPull = phiRes / phiErrorRec; double dxyPull = dxyRes / track.dxyError(); double dzPull = dzRes / track.dzError(); - - // if(track.algo()==0) - // std::cout << "reco track par:" <<"pt="<