diff --git a/DataFormats/Math/interface/invertPosDefMatrix.h b/DataFormats/Math/interface/invertPosDefMatrix.h index 42d53e6d7f745..d6b1ad3cc5af3 100644 --- a/DataFormats/Math/interface/invertPosDefMatrix.h +++ b/DataFormats/Math/interface/invertPosDefMatrix.h @@ -17,6 +17,17 @@ inline bool invertPosDefMatrix(ROOT::Math::SMatrix +inline bool invertPosDefMatrix(ROOT::Math::SMatrix > & m) { + m(0,0) = 1./m(0,0); + return true; +} +template<> +inline bool invertPosDefMatrix(ROOT::Math::SMatrix > & m) { + m(0,0) = 1.f/m(0,0); + return true; +} + template inline bool invertPosDefMatrix(ROOT::Math::SMatrix > const & mIn, ROOT::Math::SMatrix > & mOut) { diff --git a/RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.cc b/RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.cc index 8e9b221cd06eb..5c66a4dda7710 100644 --- a/RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.cc +++ b/RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.cc @@ -200,7 +200,7 @@ GoodSeedProducer::produce(Event& iEvent, const EventSetup& iSetup) auto const & Tj=*(tjCollection.product()); LogDebug("GoodSeedProducer")<<"Number of tracks in collection " - < V0coll; iEvent.getByToken(V0list_[il],V0coll); - LogDebug("PFV0Producer")<size()<<" V0 candidates "; + LogDebug("PFV0Producer")<< "V0list_[" << il <<"] contains "<size()<<" V0 candidates "; for (unsigned int iv=0;ivsize();iv++){ VertexCompositeCandidateRef V0(V0coll, iv); vector Tracks; diff --git a/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc b/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc index 85348ce2b0506..29698615027a9 100644 --- a/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc +++ b/RecoTracker/TrackProducer/src/TrackProducerAlgorithm.cc @@ -30,8 +30,10 @@ #include "TrackingTools/TrackFitters/interface/RecHitSorter.h" #include "DataFormats/TrackReco/interface/TrackBase.h" +#include + // #define VI_DEBUG -#define STAT_TSB +// #define STAT_TSB #ifdef VI_DEBUG #define DPRINT(x) std::cout << x << ": " @@ -257,33 +259,36 @@ TrackProducerAlgorithm::buildTrack (const TrajectoryFitter * the auto theTraj = new Trajectory( std::move(trajTmp) ); theTraj->setSeedRef(seedRef); - - // TrajectoryStateOnSurface innertsos; - // TrajectoryStateOnSurface outertsos; - // if (theTraj->direction() == alongMomentum) { - // innertsos = theTraj->firstMeasurement().updatedState(); - // outertsos = theTraj->lastMeasurement().updatedState(); - // } else { - // innertsos = theTraj->lastMeasurement().updatedState(); - // outertsos = theTraj->firstMeasurement().updatedState(); - // } - // std::cout - // << "Nr. of first / last states = " - // << innertsos.components().size() << " " - // << outertsos.components().size() << std::endl; - // std::vector components = - // innertsos.components(); - // double sinTheta = - // sin(innertsos.globalMomentum().theta()); - // for ( std::vector::const_iterator ic=components.begin(); - // ic!=components.end(); ic++ ) { - // std::cout << " comp " << ic-components.begin() << " " - // << (*ic).weight() << " " - // << (*ic).localParameters().vector()[0]/sinTheta << " " - // << sqrt((*ic).localError().matrix()[0][0])/sinTheta << std::endl; - // } - +#ifdef EDM_ML_DEBUG + TrajectoryStateOnSurface innertsos; + TrajectoryStateOnSurface outertsos; + + if (theTraj->direction() == alongMomentum) { + innertsos = theTraj->firstMeasurement().updatedState(); + outertsos = theTraj->lastMeasurement().updatedState(); + } else { + innertsos = theTraj->lastMeasurement().updatedState(); + outertsos = theTraj->firstMeasurement().updatedState(); + } + std::ostringstream ss; + auto dc = [&](TrajectoryStateOnSurface const & tsos){ + std::vector const & components = tsos.components(); + auto sinTheta = std::sin(tsos.globalMomentum().theta()); + for (auto const & ic : components) ss << ic.weight() << "/"; ss << "\n"; + for (auto const & ic : components) ss << ic.localParameters().vector()[0]/sinTheta << "/"; ss << "\n"; + for (auto const & ic : components) ss << std::sqrt(ic.localError().matrix()(0,0))/sinTheta << "/"; + }; + ss << "\ninner comps\n"; + dc(innertsos); + ss << "\nouter comps\n"; + dc(outertsos); + LogDebug("TrackProducer") + << "Nr. of first / last states = " + << innertsos.components().size() << " " + << outertsos.components().size() << ss.str(); +#endif + ndof = 0; for (auto const & tm : theTraj->measurements()) { auto const & h = tm.recHitR(); diff --git a/RecoTracker/TrackProducer/test/TrackRefit.py b/RecoTracker/TrackProducer/test/TrackRefit.py index 64eab2d5c8f37..009ca4aac4e57 100644 --- a/RecoTracker/TrackProducer/test/TrackRefit.py +++ b/RecoTracker/TrackProducer/test/TrackRefit.py @@ -1,80 +1,41 @@ import FWCore.ParameterSet.Config as cms - -process = cms.Process("Refitting") -### standard MessageLoggerConfiguration -process.load("FWCore.MessageService.MessageLogger_cfi") +process = cms.Process("Refitting") ### Standard Configurations -process.load("Configuration.StandardSequences.Services_cff") -process.load('Configuration/StandardSequences/GeometryIdeal_cff') -process.load('Configuration/StandardSequences/Reconstruction_cff') -process.load('Configuration/StandardSequences/MagneticField_AutoFromDBCurrent_cff') - - -## Fitter-smoother: loosen outlier rejection as for first data-taking with LHC "collisions" -process.KFFittingSmootherWithOutliersRejectionAndRK.BreakTrajWith2ConsecutiveMissing = False -process.KFFittingSmootherWithOutliersRejectionAndRK.EstimateCut = 1000 +#process.load("RecoVertex.BeamSpotProducer.BeamSpot_cfi") +#process.load("Configuration.StandardSequences.MagneticField_cff") +#process.load('Configuration.Geometry.GeometryRecoDB_cff') +#process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") + +process.load('Configuration.StandardSequences.Services_cff') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') +process.load('Configuration.StandardSequences.GeometryRecoDB_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('Configuration.StandardSequences.RawToDigi_cff') +process.load('Configuration.StandardSequences.L1Reco_cff') +process.load('Configuration.StandardSequences.Reconstruction_cff') + +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') +from Configuration.AlCa.GlobalTag import GlobalTag +# choose! +process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:run2_data_GRun', '') +# process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:run2_mc_GRun', '') - -### Conditions -process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") -#process.GlobalTag.globaltag = "IDEAL_V5::All" -process.GlobalTag.globaltag = 'GR09_P_V6::All' - ### Track refitter specific stuff -process.load("RecoTracker.TrackProducer.TrackRefitters_cff") #the correct one -#process.load("RecoTracker.TrackProducer.RefitterWithMaterial_cff") #the one for backward compatibility - - -process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) -) - -process.source = cms.Source("PoolSource", -### tracks from collisions - fileNames = cms.untracked.vstring( -'rfio:/castor/cern.ch/user/c/chiochia/09_beam_commissioning/BSCskim_123151_Express.root') -#'/store/relval/CMSSW_2_1_10/RelValTTbar/GEN-SIM-DIGI-RAW-HLTDEBUG-RECO/IDEAL_V9_v1/0001/1E04FC31-F99A-DD11-94EE-0018F3D096DE.root') - -### tracks from cosmics -# fileNames = cms.untracked.vstring( -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/005F51E5-0373-DD11-B6FA-001731AF6B7D.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/005F51E5-0373-DD11-B6FA-001731AF6B7D.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/006F3A6A-0373-DD11-A8E7-00304876A0FF.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/02CF5B1E-6476-DD11-A034-003048769E65.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/02DF31C3-A775-DD11-91C2-001A92971BB8.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/02F71F56-CE74-DD11-9DD0-001A92810AE4.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0446C89C-E072-DD11-A341-0018F3D0960C.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/04750FC3-3E73-DD11-B054-00304876A147.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/04DFD531-0473-DD11-964E-0018F3D096AE.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/067111FB-3873-DD11-AD86-00304875A9C5.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/067982F4-E175-DD11-99F7-001731AF6AC5.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0680EB9B-4F73-DD11-83F8-0018F3D0962E.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/06BF1AF3-E175-DD11-B467-00304876A147.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0A3843F3-E175-DD11-8419-003048767EE7.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0A5AAABA-3973-DD11-B949-003048767FA1.root', -# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0A911B18-0273-DD11-A5A6-001731A283E1.root') - -### tracks from beam halo muons -) - -process.TRACKS = cms.OutputModule("PoolOutputModule", - outputCommands = cms.untracked.vstring('drop *_*_*_*', - 'keep recoTracks_*_*_*', - 'keep recoTrackExtras_*_*_*', - 'keep TrackingRecHitsOwned_*_*_*'), - - fileName = cms.untracked.string('refitting.root') - ) +process.load("RecoTracker.TrackProducer.TrackRefitter_cfi") +process.TrackRefitter.NavigationSchool = '' +process.TrackRefitter.Fitter = 'FlexibleKFFittingSmoother' -process.options = cms.untracked.PSet( wantSummary = cms.untracked.bool(True) ) +process.source = cms.Source ("PoolSource", + fileNames=cms.untracked.vstring('file:pickevents_1.root', + ), + skipEvents=cms.untracked.uint32(0) + ) +process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(-1)) -process.p1 = cms.Path(process.TrackRefitter - #process.TrackRefitterP5 - #process.TrackRefitterBHM -) -process.outpath = cms.EndPath(process.TRACKS) +process.Path = cms.Path(process.TrackRefitter) - diff --git a/RecoTracker/TrackProducer/test/TrackRefitOld.py b/RecoTracker/TrackProducer/test/TrackRefitOld.py new file mode 100644 index 0000000000000..64eab2d5c8f37 --- /dev/null +++ b/RecoTracker/TrackProducer/test/TrackRefitOld.py @@ -0,0 +1,80 @@ +import FWCore.ParameterSet.Config as cms + +process = cms.Process("Refitting") + +### standard MessageLoggerConfiguration +process.load("FWCore.MessageService.MessageLogger_cfi") + +### Standard Configurations +process.load("Configuration.StandardSequences.Services_cff") +process.load('Configuration/StandardSequences/GeometryIdeal_cff') +process.load('Configuration/StandardSequences/Reconstruction_cff') +process.load('Configuration/StandardSequences/MagneticField_AutoFromDBCurrent_cff') + + +## Fitter-smoother: loosen outlier rejection as for first data-taking with LHC "collisions" +process.KFFittingSmootherWithOutliersRejectionAndRK.BreakTrajWith2ConsecutiveMissing = False +process.KFFittingSmootherWithOutliersRejectionAndRK.EstimateCut = 1000 + + + +### Conditions +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +#process.GlobalTag.globaltag = "IDEAL_V5::All" +process.GlobalTag.globaltag = 'GR09_P_V6::All' + +### Track refitter specific stuff +process.load("RecoTracker.TrackProducer.TrackRefitters_cff") #the correct one +#process.load("RecoTracker.TrackProducer.RefitterWithMaterial_cff") #the one for backward compatibility + + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) +) + +process.source = cms.Source("PoolSource", +### tracks from collisions + fileNames = cms.untracked.vstring( +'rfio:/castor/cern.ch/user/c/chiochia/09_beam_commissioning/BSCskim_123151_Express.root') +#'/store/relval/CMSSW_2_1_10/RelValTTbar/GEN-SIM-DIGI-RAW-HLTDEBUG-RECO/IDEAL_V9_v1/0001/1E04FC31-F99A-DD11-94EE-0018F3D096DE.root') + +### tracks from cosmics +# fileNames = cms.untracked.vstring( +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/005F51E5-0373-DD11-B6FA-001731AF6B7D.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/005F51E5-0373-DD11-B6FA-001731AF6B7D.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/006F3A6A-0373-DD11-A8E7-00304876A0FF.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/02CF5B1E-6476-DD11-A034-003048769E65.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/02DF31C3-A775-DD11-91C2-001A92971BB8.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/02F71F56-CE74-DD11-9DD0-001A92810AE4.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0446C89C-E072-DD11-A341-0018F3D0960C.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/04750FC3-3E73-DD11-B054-00304876A147.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/04DFD531-0473-DD11-964E-0018F3D096AE.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/067111FB-3873-DD11-AD86-00304875A9C5.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/067982F4-E175-DD11-99F7-001731AF6AC5.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0680EB9B-4F73-DD11-83F8-0018F3D0962E.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/06BF1AF3-E175-DD11-B467-00304876A147.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0A3843F3-E175-DD11-8419-003048767EE7.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0A5AAABA-3973-DD11-B949-003048767FA1.root', +# '/store/data/CRUZET4_v1/Cosmics/RECO/CRZT210_V1_SuperPointing_v1/0000/0A911B18-0273-DD11-A5A6-001731A283E1.root') + +### tracks from beam halo muons +) + +process.TRACKS = cms.OutputModule("PoolOutputModule", + outputCommands = cms.untracked.vstring('drop *_*_*_*', + 'keep recoTracks_*_*_*', + 'keep recoTrackExtras_*_*_*', + 'keep TrackingRecHitsOwned_*_*_*'), + + fileName = cms.untracked.string('refitting.root') + ) + +process.options = cms.untracked.PSet( wantSummary = cms.untracked.bool(True) ) + + +process.p1 = cms.Path(process.TrackRefitter + #process.TrackRefitterP5 + #process.TrackRefitterBHM +) +process.outpath = cms.EndPath(process.TRACKS) + + diff --git a/RecoTracker/TrackProducer/test/recoDebug.py b/RecoTracker/TrackProducer/test/recoDebug.py new file mode 100644 index 0000000000000..035f13e282e3e --- /dev/null +++ b/RecoTracker/TrackProducer/test/recoDebug.py @@ -0,0 +1,120 @@ +# Auto generated configuration file +# using: +# Revision: 1.19 +# Source: /local/reps/CMSSW/CMSSW/Configuration/Applications/python/ConfigBuilder.py,v +# with command line options: RECO --conditions 76X_dataRun2_v1 -s RAW2DIGI,L1Reco,RECO,EI,PAT,DQM --runUnscheduled --data --eventcontent RECO,AOD,MINIAOD,DQM --scenario pp --datatier RECO,AOD,MINIAOD,DQMIO --customise Configuration/DataProcessing/RecoTLR.customiseDataRun2Common_25ns --filein /store/data/Run2015C/SingleElectron/RAW/v1/000/254/879/00000/8E51CA98-7349-E511-B9AE-02163E01427B.root -n 100 --no_exec --python_filename=RECO_Run2015C.py +import FWCore.ParameterSet.Config as cms + +process = cms.Process('RECO') + +# import of standard configurations +process.load('Configuration.StandardSequences.Services_cff') +process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') +process.load('Configuration.StandardSequences.GeometryRecoDB_cff') +process.load('Configuration.StandardSequences.MagneticField_AutoFromDBCurrent_cff') +process.load('Configuration.StandardSequences.RawToDigi_Data_cff') +process.load('Configuration.StandardSequences.L1Reco_cff') +process.load('Configuration.StandardSequences.Reconstruction_Data_cff') +process.load('CommonTools.ParticleFlow.EITopPAG_cff') +process.load('DQMOffline.Configuration.DQMOffline_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + + +#Adding SimpleMemoryCheck service: +#process.SimpleMemoryCheck=cms.Service("SimpleMemoryCheck", +# ignoreTotal=cms.untracked.int32(1), +# oncePerEventMode=cms.untracked.bool(True)) + + +process.Timing = cms.Service("Timing" +# ,summaryOnly = cms.untracked.bool(True) +) + +# process.add_(cms.Service("Tracer")) + + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(100) +) + +# Input source +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring('file:jetHT256630_RAW.root'), + secondaryFileNames = cms.untracked.vstring() +) + +process.options = cms.untracked.PSet( +) + +# Production Info +process.configurationMetadata = cms.untracked.PSet( + annotation = cms.untracked.string('RECO nevts:100'), + name = cms.untracked.string('Applications'), + version = cms.untracked.string('$Revision: 1.19 $') +) + +# Output definition + +process.RECOoutput = cms.OutputModule("PoolOutputModule", + dataset = cms.untracked.PSet( + dataTier = cms.untracked.string('RECO'), + filterName = cms.untracked.string('') + ), + eventAutoFlushCompressedSize = cms.untracked.int32(5242880), + fileName = cms.untracked.string('RECO_RAW2DIGI_L1Reco_RECO_EI_PAT_DQM.root'), + outputCommands = process.RECOEventContent.outputCommands, + splitLevel = cms.untracked.int32(0) +) + +# Additional output definition + +# Other statements +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, + 'auto:run2_data_GRun', '') # if mc change in mc.... + +# Path and EndPath definitions +process.raw2digi_step = cms.Path(process.RawToDigi) +process.L1Reco_step = cms.Path(process.L1Reco) +process.reconstruction_step = cms.Path(process.reconstruction) +process.eventinterpretaion_step = cms.Path(process.EIsequence) +process.RECOoutput_step = cms.EndPath(process.RECOoutput) + +# Schedule definition +process.schedule = cms.Schedule(process.raw2digi_step,process.L1Reco_step,process.reconstruction_step,process.eventinterpretaion_step,process.RECOoutput_step) + +# customisation of the process. + +# Automatic addition of the customisation function from Configuration.DataProcessing.RecoTLR +from Configuration.DataProcessing.RecoTLR import customiseDataRun2Common_25ns + +#call to customisation function customiseDataRun2Common_25ns imported from Configuration.DataProcessing.RecoTLR +process = customiseDataRun2Common_25ns(process) + +# End of customisation functions + +process.MessageLogger = cms.Service("MessageLogger", + destinations = cms.untracked.vstring("cout"), #1 + debugModules = cms.untracked.vstring("electronGsfTracks"),#initialStepTrackCandidates), #2 + categories = cms.untracked.vstring( + 'TrackProducer', + 'GsfTrackFitters', + # 'GsfMaterialEffectsUpdator', + 'AnalyticalPropagator','RKPropagatorInS', + # 'CkfPattern' + ), #3 + cout = cms.untracked.PSet(threshold = cms.untracked.string("DEBUG"), #4 + DEBUG = cms.untracked.PSet(limit = cms.untracked.int32(0)), #5 + default = cms.untracked.PSet(limit = cms.untracked.int32(0)), #6 + TrackProducer = cms.untracked.PSet(limit = cms.untracked.int32(-1)), #7 + GsfTrackFitters = cms.untracked.PSet(limit = cms.untracked.int32(-1)), #7 + # GsfMaterialEffectsUpdator = cms.untracked.PSet(limit = cms.untracked.int32(-1)), #7 + AnalyticalPropagator = cms.untracked.PSet(limit = cms.untracked.int32(-1)), #7 + RKPropagatorInS = cms.untracked.PSet(limit = cms.untracked.int32(-1)), #7 + # CkfPattern = cms.untracked.PSet(limit = cms.untracked.int32(-1)) #7 + ) + ) + + diff --git a/TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.h b/TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.h index ec29ec4ee86eb..635d238efb513 100644 --- a/TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.h +++ b/TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.h @@ -65,15 +65,12 @@ class MultiGaussianStateAssembler { private: const MultiState theInitialState; -// bool sortStates; double minFractionalWeight; bool combinationDone; double theValidWeightSum; SingleStateContainer theStates; -// static TimingReport::Item * theTimerAdd; -// static TimingReport::Item * theTimerComb; }; diff --git a/TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.icc b/TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.icc index 5e330a1a8f8aa..3db2235eee0a8 100644 --- a/TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.icc +++ b/TrackingTools/GsfTools/interface/MultiGaussianStateAssembler.icc @@ -1,30 +1,12 @@ #include "TrackingTools/GsfTools/interface/GaussianStateLessWeight.h" #include "FWCore/Utilities/interface/Exception.h" + template MultiGaussianStateAssembler::MultiGaussianStateAssembler (const MultiState & state) : theInitialState(state), minFractionalWeight(1.e-16), combinationDone(false), theValidWeightSum(0.) -{ - // - // parameters (could be configurable) - // -// static SimpleConfigurable sortStatesConf(false,"MultiGaussianStateAssembler:sortStates"); -// sortStates = sortStatesConf.value(); -// static SimpleConfigurable minWeightConf(1.e-16,"MultiGaussianStateAssembler:minFractionalWeight"); -// minFractionalWeight = minWeightConf.value(); -// // -// // Timers -// // -// if ( theTimerAdd==0 ) { -// theTimerAdd = -// &(*TimingReport::current())[string("MultiGaussianStateAssembler::addState")]; -// theTimerAdd->switchCPU(false); -// theTimerComb = -// &(*TimingReport::current())[string("MultiGaussianStateAssembler::combinedState")]; -// theTimerComb->switchCPU(false); -// } -} +{} template void MultiGaussianStateAssembler::addState (const SingleStatePtr& sgs) @@ -44,8 +26,6 @@ void MultiGaussianStateAssembler::addState (const SingleStatePtr& sgs) template void MultiGaussianStateAssembler::addState (const MultiState& mgs) { -// // Timer -// TimeMe t(*theTimerAdd,false); // // refuse to add states after combination has been done // @@ -53,12 +33,6 @@ void MultiGaussianStateAssembler::addState (const MultiState& mgs) { throw cms::Exception("LogicError") << "MultiGaussianStateAssembler: trying to add states after combination"; // - // Verify validity of state to be added - //ThS: Is not possible for gaussian states. -// if ( !tsos.isValid() ) -// throw cms::Exception("LogicError") -// << "MultiGaussianStateAssembler: trying to add invalid state"; - // // Add components (i.e. state to be added can be single or multi state) // SingleStateContainer components(mgs.components()); @@ -77,15 +51,8 @@ void MultiGaussianStateAssembler::addStateVector (const SingleStateContainer& // // sum up weights (all components are supposed to be valid!!!) // - double sum(0.); - for ( typename SingleStateContainer::const_iterator i=states.begin(); - i!=states.end(); i++ ) { - //ThS: Is not possible for gaussian states: -// if ( !(i->isValid()) ) -// throw cms::Exception("LogicError") -// << "MultiGaussianStateAssembler: trying to add invalid state"; - sum += (**i).weight(); - } + double sum=0.; + for (auto const & is : states ) sum += (*is).weight(); theValidWeightSum += sum; // // add to vector of states @@ -97,31 +64,19 @@ void MultiGaussianStateAssembler::addStateVector (const SingleStateContainer& template MultiGaussianState MultiGaussianStateAssembler::combinedState () { -// // Timer -// TimeMe t(*theTimerComb,false); // // Prepare resulting state vector prepareCombinedState(); - // ThS: What to return here? I choose a RCMP with empty state vector... -// if ( !prepareCombinedState() ) return TSOS(); - // - // Return new multi state without reweighting - // -// return theInitialState.createNewState(theStates); return MultiState(theStates); } template MultiGaussianState MultiGaussianStateAssembler::combinedState (const float newWeight) { -// // Timer -// TimeMe t(*theTimerComb,false); // // Prepare resulting state vector // prepareCombinedState(); - // ThS: What to return here? I choose a RCMP with empty state vector... -// if ( !prepareCombinedState() ) return TSOS(); // // return reweighted state // @@ -135,19 +90,12 @@ MultiGaussianStateAssembler::prepareCombinedState () { // remaining part to be done only once // if ( combinationDone ) return true; - else combinationDone = true; + combinationDone = true; // // Remove states with negligible weights // removeSmallWeights(); - if ( theStates.empty() ) return false; -// // -// // Sort output by weights? (currently deactivated - make it again configurable?) -// // -// if ( sortStates ) -// sort(theStates.begin(),theStates.end(),GaussianStateLessWeight()); - - return true; + return !theStates.empty(); } template @@ -162,16 +110,11 @@ MultiGaussianStateAssembler::reweightedCombinedState (const double newWeight) // SingleStateContainer reweightedStates; reweightedStates.reserve(theStates.size()); - for ( typename SingleStateContainer::const_iterator i=theStates.begin(); - i!=theStates.end(); i++ ) { - double oldWeight = (**i).weight(); -// reweightedStates.push_back( (**i).createNewState ((**i).mean(), -// (**i).covariance(), factor*oldWeight) ); - reweightedStates.push_back( SingleState((**i).mean(), - (**i).covariance(), factor*oldWeight) ); + for (auto const & is : theStates) { + auto oldWeight = (*is).weight(); + reweightedStates.emplace_back((*is).mean(), + (*is).covariance(), factor*oldWeight); } - -// return theInitialState.createNewState(reweightedStates); return MultiState(reweightedStates); } @@ -186,22 +129,7 @@ MultiGaussianStateAssembler::removeSmallWeights() theStates.clear(); return; } - // - // Loop until no more states are removed - // - bool redo; - do { - redo = false; - for ( typename SingleStateContainer::iterator i=theStates.begin(); - i!=theStates.end(); i++ ) { - if ( (**i).weight()/theValidWeightSum < minFractionalWeight ) { - theStates.erase(i); - redo = true; - break; - } - } - } while (redo); + theStates.erase(std::remove_if(theStates.begin(),theStates.end(), + [&](typename SingleStateContainer::value_type const & s){ return s->weight() < minFractionalWeight*theValidWeightSum;}),theStates.end() + ); } - -// TimingReport::Item * MultiGaussianStateAssembler::theTimerAdd(0); -// TimingReport::Item * MultiGaussianStateAssembler::theTimerComb(0); diff --git a/TrackingTools/GsfTools/src/BasicMultiTrajectoryState.cc b/TrackingTools/GsfTools/src/BasicMultiTrajectoryState.cc index b049dfe66c13d..a6560787b4c69 100644 --- a/TrackingTools/GsfTools/src/BasicMultiTrajectoryState.cc +++ b/TrackingTools/GsfTools/src/BasicMultiTrajectoryState.cc @@ -46,9 +46,7 @@ void BasicMultiTrajectoryState::rescaleError(double factor) { return; } - for (std::vector::iterator it = theStates.begin(); it != theStates.end(); it++) { - it->rescaleError(factor); - } + for (auto & is : theStates) is.rescaleError(factor); combine(); } @@ -88,8 +86,8 @@ BasicMultiTrajectoryState::combine() { AlgebraicVector5 mean; AlgebraicSymMatrix55 covarPart1, covarPart2, covtmp; for (auto it1 = tsos.begin(); it1 != tsos.end(); it1++) { - double weight = it1->weight(); - AlgebraicVector5 param = it1->localParameters().vector(); + auto weight = it1->weight(); + auto const & param = it1->localParameters().vector(); sumw += weight; mean += weight * param; covarPart1 += weight * it1->localError().matrix(); diff --git a/TrackingTools/GsfTools/src/MultiGaussianStateCombiner1D.cc b/TrackingTools/GsfTools/src/MultiGaussianStateCombiner1D.cc index 5f07b8f76c977..6deaa5c892c2e 100644 --- a/TrackingTools/GsfTools/src/MultiGaussianStateCombiner1D.cc +++ b/TrackingTools/GsfTools/src/MultiGaussianStateCombiner1D.cc @@ -69,9 +69,10 @@ MultiGaussianStateCombiner1D::combine(const VSC& theComponents) const measCovar = 0.; weightSum = 0.; } else { - meanMean /= weightSum; - measCovar1 *= (1./weightSum); - measCovar2 *= (1./weightSum/weightSum); + weightSum = 1./weightSum; + meanMean *= weightSum; + measCovar1 *= weightSum; + measCovar2 *= weightSum*weightSum; measCovar = measCovar1 + measCovar2; // measCovar = measCovar1/weightSum + measCovar2/weightSum/weightSum; } diff --git a/TrackingTools/GsfTools/src/MultiStatePropagation.icc b/TrackingTools/GsfTools/src/MultiStatePropagation.icc index e4e0c81a7add2..eb5eabd3bae51 100644 --- a/TrackingTools/GsfTools/src/MultiStatePropagation.icc +++ b/TrackingTools/GsfTools/src/MultiStatePropagation.icc @@ -11,7 +11,7 @@ MultiStatePropagation::propagateWithPath (const TrajectoryStateOnSurface& tso double sWeight(0.); double sWeightPath(0.); // decomposition of input state - MultiTSOS input(tsos.components()); + MultiTSOS const & input = tsos.components(); // vector of result states MultiTrajectoryStateAssembler result; // @@ -19,19 +19,19 @@ MultiStatePropagation::propagateWithPath (const TrajectoryStateOnSurface& tso // bool firstPropagation(true); SurfaceSideDefinition::SurfaceSide firstSide(SurfaceSideDefinition::atCenterOfSurface); - for ( MultiTSOS::const_iterator iTsos=input.begin(); - iTsos!=input.end(); iTsos++ ) { + for ( auto const & iTsos : input) { // // weight of component // - double weight(iTsos->weight()); + double weight(iTsos.weight()); // // geometrical propagation (assumption: only one output state!) // TsosWP newTsosWP = - thePropagator.propagateWithPath(*iTsos,surface); + thePropagator.propagateWithPath(iTsos,surface); // check validity if ( !(newTsosWP.first).isValid() ) { + LogDebug("GsfTrackFitters") << "adding invalid state"; result.addInvalidState(weight); continue; } @@ -44,7 +44,7 @@ MultiStatePropagation::propagateWithPath (const TrajectoryStateOnSurface& tso } else { if ( weightedTsos.surfaceSide()!=firstSide ) { - edm::LogInfo("MultiStatePropagation") + edm::LogInfo("GsfTrackFitters") << "MultiStatePropagation resulted in states " << "on different sides of the target surface"; return TsosWP(TrajectoryStateOnSurface(),0.); @@ -77,14 +77,9 @@ template TrajectoryStateOnSurface MultiStatePropagation::setWeight (const TrajectoryStateOnSurface state, const double weight) const { - // if ( state.hasError() ) return TrajectoryStateOnSurface(weight, state.localParameters(), state.localError(), state.surface(), &(state.globalParameters().magneticField()), state.surfaceSide()); - // else - // return TrajectoryStateOnSurface(weight, state.localParameters(), - // state.surface()); - // return state; } diff --git a/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc b/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc index 071d6fa049582..4ef4f1d1a5fc8 100644 --- a/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc +++ b/TrackingTools/GsfTools/src/MultiTrajectoryStateAssembler.cc @@ -16,23 +16,10 @@ MultiTrajectoryStateAssembler::MultiTrajectoryStateAssembler () : // sortStates = false; minValidFraction = 0.01; - minFractionalWeight = 1.e-6; - // // - // // Timers - // // - // if ( theTimerAdd==0 ) { - // theTimerAdd = - // &(*TimingReport::current())[string("MultiTrajectoryStateAssembler::addState")]; - // theTimerAdd->switchCPU(false); - // theTimerComb = - // &(*TimingReport::current())[string("MultiTrajectoryStateAssembler::combinedState")]; - // theTimerComb->switchCPU(false); - // } + minFractionalWeight = 1.e-6; // 4; } void MultiTrajectoryStateAssembler::addState (const TrajectoryStateOnSurface tsos) { - // // Timer - // TimeMe t(*theTimerAdd,false); // // refuse to add states after combination has been done // @@ -92,8 +79,6 @@ void MultiTrajectoryStateAssembler::addInvalidState (const double weight) { } TrajectoryStateOnSurface MultiTrajectoryStateAssembler::combinedState () { - // // Timer - // TimeMe t(*theTimerComb,false); // // Prepare resulting state vector // @@ -110,8 +95,6 @@ TrajectoryStateOnSurface MultiTrajectoryStateAssembler::combinedState () { } TrajectoryStateOnSurface MultiTrajectoryStateAssembler::combinedState (const float newWeight) { - // // Timer - // TimeMe t(*theTimerComb,false); // // Prepare resulting state vector // @@ -136,12 +119,12 @@ MultiTrajectoryStateAssembler::prepareCombinedState () { // Check for minimum fraction of valid states // double allWeights(theValidWeightSum+theInvalidWeightSum); - if ( theInvalidWeightSum>0. && (theValidWeightSum/allWeights)0. && theValidWeightSum weight(); - reweightedStates.push_back(TrajectoryStateOnSurface(factor*oldWeight, - i->localParameters(), - i->localError(), - i->surface(), - &(i->globalParameters().magneticField()), - i->surfaceSide() - )); + for (auto const & is : theStates) { + auto oldWeight = is.weight(); + reweightedStates.emplace_back(factor*oldWeight, + is.localParameters(), + is.localError(), + is.surface(), + &(is.globalParameters().magneticField()), + is.surfaceSide() + ); } return TSOS((BasicTrajectoryState *)(new BasicMultiTrajectoryState(reweightedStates))); } @@ -191,70 +173,52 @@ MultiTrajectoryStateAssembler::removeSmallWeights() // // check total weight // - double totalWeight(theInvalidWeightSum+theValidWeightSum); + auto totalWeight(theInvalidWeightSum+theValidWeightSum); if ( totalWeight == 0. ) { theStates.clear(); return; } - // - // Loop until no more states are removed - // - bool redo; - do { - redo = false; - for ( MultiTSOS::iterator i=theStates.begin(); - i!=theStates.end(); i++ ) { - if ( (*i).weight()/totalWeight < minFractionalWeight ) { - theStates.erase(i); - redo = true; - break; - } - } - } while (redo); + theStates.erase(std::remove_if(theStates.begin(),theStates.end(), + [&](MultiTSOS::value_type const & s){ return s.weight() < minFractionalWeight*totalWeight;}), + theStates.end()); } void MultiTrajectoryStateAssembler::removeWrongPz () { - // edm::LogDebug("MultiTrajectoryStateAssembler") - // << "MultiTrajectoryStateAssembler: found at least one state with inconsistent pz\n" - // << " #state / weights before cleaning = " << theStates.size() - // << " / " << theValidWeightSum - // << " / " << theInvalidWeightSum; + LogDebug("GsfTrackFitters") + << "MultiTrajectoryStateAssembler: found at least one state with inconsistent pz\n" + << " #state / weights before cleaning = " << theStates.size() + << " / " << theValidWeightSum + << " / " << theInvalidWeightSum; // // Calculate average pz // double meanPz(0.); - for ( MultiTSOS::const_iterator is=theStates.begin(); - is!=theStates.end(); is++ ) { - meanPz += is->weight()*is->localParameters().pzSign(); - // edm::LogDebug("MultiTrajectoryStateAssembler") - // << " weight / pz / global position = " << is->weight() - // << " " << is->localParameters().pzSign() - // << " " << is->globalPosition(); - } + for (auto const & is : theStates) + meanPz += is.weight()*is.localParameters().pzSign(); meanPz /= theValidWeightSum; // // Now keep only states compatible with the average pz // - // double oldValidWeight(theValidWeightSum); theValidWeightSum = 0.; MultiTSOS oldStates(theStates); theStates.clear(); - for ( MultiTSOS::const_iterator is=oldStates.begin(); - is!=oldStates.end(); is++ ) { - if ( meanPz*is->localParameters().pzSign()>=0. ) { - theValidWeightSum += is->weight(); - theStates.push_back(*is); + for (auto const & is :oldStates) { + if ( meanPz*is.localParameters().pzSign()>=0. ) { + theValidWeightSum += is.weight(); + theStates.push_back(std::move(is)); } else { - theInvalidWeightSum += is->weight(); + theInvalidWeightSum += is.weight(); + LogDebug("GsfTrackFitters") + << "removing weight / pz / global position = " << is.weight() + << " " << is.localParameters().pzSign() + << " " << is.globalPosition(); + } } - // edm::LogDebug("MultiTrajectoryStateAssembler") - // << " #state / weights after cleaning = " << theStates.size() - // << " / " << theValidWeightSum - // << " / " << theInvalidWeightSum; + LogDebug("GsfTrackFitters") + << " #state / weights after cleaning = " << theStates.size() + << " / " << theValidWeightSum + << " / " << theInvalidWeightSum; } - -// TimingReport::Item * MultiTrajectoryStateAssembler::theTimerAdd(0); -// TimingReport::Item * MultiTrajectoryStateAssembler::theTimerComb(0); diff --git a/TrackingTools/GsfTools/src/MultiTrajectoryStateMode.cc b/TrackingTools/GsfTools/src/MultiTrajectoryStateMode.cc index 8b3b2f8bdc055..67183474e49a2 100644 --- a/TrackingTools/GsfTools/src/MultiTrajectoryStateMode.cc +++ b/TrackingTools/GsfTools/src/MultiTrajectoryStateMode.cc @@ -278,28 +278,28 @@ MultiTrajectoryStateMode::momentumFromModePPhiEta (const TrajectoryStateOnSurfac ic!=components.end(); ++ic ) { // parameters GlobalVector mom(ic->globalMomentum()); - double px = mom.x(); - double py = mom.y(); - double pz = mom.z(); - double p = mom.mag(); - double pt2 = mom.perp2(); - double phi = mom.phi(); - double eta = mom.eta(); + auto px = mom.x(); + auto py = mom.y(); + auto pz = mom.z(); + auto op = 1./mom.mag(); + auto opt2 = 1./mom.perp2(); + auto phi = mom.phi(); + auto eta = mom.eta(); // jacobian - jacobian(0,0) = px/p; - jacobian(0,1) = py/p; - jacobian(0,2) = pz/p; - jacobian(1,0) = py/pt2; - jacobian(1,1) = -px/pt2; + jacobian(0,0) = px*op; + jacobian(0,1) = py*op; + jacobian(0,2) = pz*op; + jacobian(1,0) = py*opt2; + jacobian(1,1) = -px*opt2; jacobian(1,2) = 0; - jacobian(2,0) = px*pz/(pt2*p); - jacobian(2,1) = py*pz/(pt2*p); - jacobian(2,2) = -1./p; + jacobian(2,0) = px*pz*opt2*op; + jacobian(2,1) = py*pz*opt2*op; + jacobian(2,2) = -op; // extraction of the momentum part from the 6x6 cartesian error matrix // and conversion to p-phi-eta covCart = ic->cartesianError().matrix().Sub(3,3); covPPhiEta = ROOT::Math::Similarity(jacobian,covCart); - pStates.push_back(SingleGaussianState1D(p,covPPhiEta(0,0),ic->weight())); + pStates.push_back(SingleGaussianState1D(1/op,covPPhiEta(0,0),ic->weight())); phiStates.push_back(SingleGaussianState1D(phi,covPPhiEta(1,1),ic->weight())); etaStates.push_back(SingleGaussianState1D(eta,covPPhiEta(2,2),ic->weight())); } @@ -315,13 +315,13 @@ MultiTrajectoryStateMode::momentumFromModePPhiEta (const TrajectoryStateOnSurfac // // parameters from mode (in case of failure: mean) // - double p = pUtils.modeIsValid() ? pUtils.mode().mean() : pUtils.mean(); - double phi = phiUtils.modeIsValid() ? phiUtils.mode().mean() : phiUtils.mean(); - double eta = etaUtils.modeIsValid() ? etaUtils.mode().mean() : etaUtils.mean(); + auto p = pUtils.modeIsValid() ? pUtils.mode().mean() : pUtils.mean(); + auto phi = phiUtils.modeIsValid() ? phiUtils.mode().mean() : phiUtils.mean(); + auto eta = etaUtils.modeIsValid() ? etaUtils.mode().mean() : etaUtils.mean(); // double theta = 2*atan(exp(-eta)); - double tanth2 = exp(-eta); - double pt = p*2*tanth2/(1+tanth2*tanth2); // p*sin(theta) - double pz = p*(1-tanth2*tanth2)/(1+tanth2*tanth2); // p*cos(theta) + auto tanth2 = std::exp(-eta); + auto pt = p*2*tanth2/(1+tanth2*tanth2); // p*sin(theta) + auto pz = p*(1-tanth2*tanth2)/(1+tanth2*tanth2); // p*cos(theta) // conversion to a cartesian momentum vector momentum = GlobalVector(pt*cos(phi),pt*sin(phi),pz); return true; diff --git a/TrackingTools/GsfTracking/interface/GsfBetheHeitlerUpdator.h b/TrackingTools/GsfTracking/interface/GsfBetheHeitlerUpdator.h index a7a9c604ed6d9..906314076b3c5 100644 --- a/TrackingTools/GsfTracking/interface/GsfBetheHeitlerUpdator.h +++ b/TrackingTools/GsfTracking/interface/GsfBetheHeitlerUpdator.h @@ -6,10 +6,6 @@ #include "DataFormats/Math/interface/approx_exp.h" #include "DataFormats/Math/interface/approx_log.h" - - -#include "TrackingTools/GsfTracking/interface/Triplet.h" - #include #include @@ -69,7 +65,7 @@ class GsfBetheHeitlerUpdator final: public GsfMaterialEffectsUpdator { GsfBetheHeitlerUpdator (const std::string fileName, const int correctionFlag); private: - typedef Triplet GSContainer; + struct GSContainer {float *first, *second, *third;}; /// Computation: generates vectors of weights, means and standard deviations virtual void compute (const TrajectoryStateOnSurface&, const PropagationDirection, Effect[]) const; @@ -83,13 +79,13 @@ class GsfBetheHeitlerUpdator final: public GsfMaterialEffectsUpdator { /// Filling of mixture (in terms of z=E/E0) - void getMixtureParameters (const float, GSContainer[]) const; + void getMixtureParameters (const float, GSContainer &) const; /// Correction for weight of component 1 - void correctWeights (GSContainer[]) const; + void correctWeights (GSContainer&) const; /// Correction for mean of component 1 - float correctedFirstMean (const float, const GSContainer[]) const; + float correctedFirstMean (const float, const GSContainer &) const; /// Correction for variance of component 1 - float correctedFirstVar (const float,const GSContainer[]) const; + float correctedFirstVar (const float,const GSContainer &) const; private: diff --git a/TrackingTools/GsfTracking/interface/Triplet.h b/TrackingTools/GsfTracking/interface/Triplet.h deleted file mode 100644 index ee4d889fe9bed..0000000000000 --- a/TrackingTools/GsfTracking/interface/Triplet.h +++ /dev/null @@ -1,47 +0,0 @@ -#ifndef Triplet_H -#define Triplet_H - -/** triplet is identical to stl::pair exceot that it has three members - * instead of two (first, second, third). - */ - -template -struct Triplet { - typedef T1 first_type; - typedef T2 second_type; - typedef T3 third_type; - - T1 first; - T2 second; - T3 third; - Triplet() : first(T1()), second(T2()), third(T3()) {} - Triplet(const T1& a, const T2& b, const T3& c) : - first(a), second(b), third(c) {} - - template - Triplet(const Triplet& p) : - first(p.first), second(p.second), third(p.third) {} -}; - -template -inline bool operator==(const Triplet& x, - const Triplet& y) { - return x.first == y.first && x.second == y.second && x.third==y.third; -} - -template -inline bool operator<(const Triplet& x, - const Triplet& y) { - bool pair_less = - x.first < y.first || (!(y.first < x.first) && x.second < y.second); - - return pair_less || - (!(y.first < x.first) && !(y.second < x.second) && x.third < y.third); -} - -template -inline Triplet make_Triplet(const T1& x, const T2& y, const T3& z) { - return Triplet(x, y, z); -} - -#endif diff --git a/TrackingTools/GsfTracking/src/DebugHelpers.h b/TrackingTools/GsfTracking/src/DebugHelpers.h new file mode 100644 index 0000000000000..c88aee50f0076 --- /dev/null +++ b/TrackingTools/GsfTracking/src/DebugHelpers.h @@ -0,0 +1,96 @@ +#ifdef EDM_ML_DEBUG + +#include "DataFormats/MuonDetId/interface/CSCDetId.h" +#include "DataFormats/MuonDetId/interface/DTWireId.h" +#include "DataFormats/MuonDetId/interface/RPCDetId.h" +#include "DataFormats/MuonDetId/interface/GEMDetId.h" +#include "DataFormats/MuonDetId/interface/ME0DetId.h" +#include "DataFormats/MuonDetId/interface/MuonSubdetId.h" +#include "DataFormats/SiStripDetId/interface/TIBDetId.h" +#include "DataFormats/SiStripDetId/interface/TOBDetId.h" +#include "DataFormats/SiStripDetId/interface/TIDDetId.h" +#include "DataFormats/SiStripDetId/interface/TECDetId.h" +#include "DataFormats/SiPixelDetId/interface/PXBDetId.h" +#include "DataFormats/SiPixelDetId/interface/PXFDetId.h" + +namespace { + inline + void dump(TrackingRecHit const & hit, int hitcounter) { + if (hit.isValid()) { + LogTrace("GsfTrackFitters")<< " ----------------- HIT #" << hitcounter << " (VALID)-----------------------\n" + << " HIT IS AT R " << hit.globalPosition().perp() << "\n" + << " HIT IS AT Z " << hit.globalPosition().z() << "\n" + << " HIT IS AT Phi " << hit.globalPosition().phi() << "\n" + << " HIT IS AT Loc " << hit.localPosition() << "\n" + << " WITH LocError " << hit.localPositionError() << "\n" + << " HIT IS AT Glo " << hit.globalPosition() << "\n" + << "SURFACE POSITION" << "\n" + << hit.surface()->position()<<"\n" + << "SURFACE ROTATION" << "\n" + << hit.surface()->rotation() + << "dimension " << hit.dimension(); + + DetId hitId = hit.geographicalId(); + + LogDebug("GsfTrackFitters") << " hit det=" << hitId.rawId(); + + if(hitId.det() == DetId::Tracker) { + if (hitId.subdetId() == StripSubdetector::TIB ) + LogDebug("GsfTrackFitters") << " I am TIB " << TIBDetId(hitId).layer(); + else if (hitId.subdetId() == StripSubdetector::TOB ) + LogDebug("GsfTrackFitters") << " I am TOB " << TOBDetId(hitId).layer(); + else if (hitId.subdetId() == StripSubdetector::TEC ) + LogDebug("GsfTrackFitters") << " I am TEC " << TECDetId(hitId).wheel(); + else if (hitId.subdetId() == StripSubdetector::TID ) + LogDebug("GsfTrackFitters") << " I am TID " << TIDDetId(hitId).wheel(); + else if (hitId.subdetId() == (int) PixelSubdetector::PixelBarrel ) + LogDebug("GsfTrackFitters") << " I am PixBar " << PXBDetId(hitId).layer(); + else if (hitId.subdetId() == (int) PixelSubdetector::PixelEndcap ) + LogDebug("GsfTrackFitters") << " I am PixFwd " << PXFDetId(hitId).disk(); + else + LogDebug("GsfTrackFitters") << " UNKNOWN TRACKER HIT TYPE "; + } + else if(hitId.det() == DetId::Muon) { + if(hitId.subdetId() == MuonSubdetId::DT) + LogDebug("GsfTrackFitters") << " I am DT " << DTWireId(hitId); + else if (hitId.subdetId() == MuonSubdetId::CSC ) + LogDebug("GsfTrackFitters") << " I am CSC " << CSCDetId(hitId); + else if (hitId.subdetId() == MuonSubdetId::RPC ) + LogDebug("GsfTrackFitters") << " I am RPC " << RPCDetId(hitId); + else if (hitId.subdetId() == MuonSubdetId::GEM ) + LogDebug("GsfTrackFitters") << " I am GEM " << GEMDetId(hitId); + + else if (hitId.subdetId() == MuonSubdetId::ME0 ) + LogDebug("GsfTrackFitters") << " I am ME0 " << ME0DetId(hitId); + else + LogDebug("GsfTrackFitters") << " UNKNOWN MUON HIT TYPE "; + } + else + LogDebug("GsfTrackFitters") << " UNKNOWN HIT TYPE "; + + } else { + LogDebug("GsfTrackFitters") + << " ----------------- INVALID HIT #" << hitcounter << " -----------------------"; + } + } +} +#include + inline void dump(TrajectoryStateOnSurface const & tsos, const char * header) { + std::ostringstream ss; ss<< " weights "; + for (auto const & c : tsos.components()) ss << c.weight() << '/'; + ss << "\nmomentums "; + for (auto const & c : tsos.components()) ss << c.globalMomentum().mag() << '/'; + ss << "\ndeltap/p "; + for (auto const & c : tsos.components()) ss << std::sqrt(tsos.curvilinearError().matrix()(0,0))/c.globalMomentum().mag() << '/'; + LogTrace("GsfTrackFitters") + << header << "! size " << tsos.components().size() << ss.str() << "\n" + <<" with local position " << tsos.localPosition() << "\n" + << tsos; + } +#else +namespace { + inline void dump(TrackingRecHit const &, int) {} + inline void dump(TrajectoryStateOnSurface const &, const char *){} +} +#endif + diff --git a/TrackingTools/GsfTracking/src/FullConvolutionWithMaterial.cc b/TrackingTools/GsfTracking/src/FullConvolutionWithMaterial.cc index f3a8d50d9e9fc..14bc37ab03914 100644 --- a/TrackingTools/GsfTracking/src/FullConvolutionWithMaterial.cc +++ b/TrackingTools/GsfTracking/src/FullConvolutionWithMaterial.cc @@ -1,5 +1,4 @@ #include "TrackingTools/GsfTracking/interface/FullConvolutionWithMaterial.h" - #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateAssembler.h" TrajectoryStateOnSurface @@ -8,8 +7,7 @@ FullConvolutionWithMaterial::operator() (const TrajectoryStateOnSurface& tsos, // // decomposition of input state // - typedef std::vector MultiTSOS; - MultiTSOS input(tsos.components()); + auto && input = tsos.components(); // // vector of result states // @@ -17,20 +15,15 @@ FullConvolutionWithMaterial::operator() (const TrajectoryStateOnSurface& tsos, // // now add material effects to each component // - for ( MultiTSOS::const_iterator iTsos=input.begin(); - iTsos!=input.end(); iTsos++ ) { - // + for (auto const & iTsos : input) { // add material // TrajectoryStateOnSurface updatedTSoS = - theMEUpdator->updateState(*iTsos,propDir); + theMEUpdator->updateState(iTsos,propDir); if ( updatedTSoS.isValid() ) result.addState(updatedTSoS); else - result.addInvalidState(iTsos->weight()); + result.addInvalidState(iTsos.weight()); } return result.combinedState(); } - -// TimingReport::Item* FullConvolutionWithMaterial::theTimer1(0); -// TimingReport::Item* FullConvolutionWithMaterial::theTimer2(0); diff --git a/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc b/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc index 1758715a5fcd2..50bf93c418f48 100644 --- a/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc +++ b/TrackingTools/GsfTracking/src/GsfBetheHeitlerUpdator.cc @@ -7,8 +7,8 @@ #include #include #include -#include +#include namespace { /// Logistic function (needed for transformation of weight and mean) @@ -20,7 +20,12 @@ namespace { /// Second moment of the Bethe-Heitler distribution (in z=E/E0) inline float BetheHeitlerVariance (const float rl) { - const float l3ol2 = std::log(3.)/std::log(2.); +#if __clang__ + const +#else + constexpr +#endif + float l3ol2 = std::log(3.)/std::log(2.); float mean = BetheHeitlerMean(rl); return unsafe_expf<4>(-rl*l3ol2) - mean*mean; } @@ -37,7 +42,12 @@ namespace { /// Second moment of the Bethe-Heitler distribution (in z=E/E0) inline float BetheHeitlerVariance (const float rl) { - constexpr float l3ol2 = std::log(3.)/std::log(2.); +#if __clang__ + const +#else + constexpr +#endif + float l3ol2 = std::log(3.)/std::log(2.); return std::exp(-rl*l3ol2) - std::exp(-2*rl); } } @@ -122,40 +132,37 @@ GsfBetheHeitlerUpdator::compute (const TrajectoryStateOnSurface& TSoS, if ( rl<0.01f ) rl = 0.01f; if ( rl>0.20f ) rl = 0.20f; - #if __clang__ - std::vector mixtureHolder(theNrComponents); - GSContainer *mixture = mixtureHolder.data(); - #else - GSContainer mixture[theNrComponents]; - #endif + float mixtureData[3][theNrComponents]; + GSContainer mixture{mixtureData[0],mixtureData[1],mixtureData[2]}; + getMixtureParameters(rl,mixture); correctWeights(mixture); if ( theCorrectionFlag>=1 ) - mixture[0].second = correctedFirstMean(rl,mixture); + mixture.second[0] = correctedFirstMean(rl,mixture); if ( theCorrectionFlag>=2 ) - mixture[0].third = correctedFirstVar(rl,mixture); + mixture.third[0] = correctedFirstVar(rl,mixture); for ( int i=0; i(vz); - else vz = vz*vz; - - mixture[i]=Triplet(weight,z,vz); + weight[i] = thePolyWeights[i](rl); + z[i] = thePolyMeans[i](rl); + vz[i] = thePolyVars[i](rl); + } + if ( theTransformationCode ) + for ( int i=0; i(vz[i]);; + } + else // theTransformationCode + for ( int i=0; i && predictedComponents = tsos.components(); + auto && predictedComponents = tsos.components(); if (predictedComponents.empty()) { edm::LogError("GsfMultiStateUpdator") << "Trying to update trajectory state with zero components! " ; return TrajectoryStateOnSurface(); } - std::vector && weights = PosteriorWeightsCalculator(predictedComponents).weights(aRecHit); + auto && weights = PosteriorWeightsCalculator(predictedComponents).weights(aRecHit); if ( weights.empty() ) { edm::LogError("GsfMultiStateUpdator") << " no weights could be retreived. invalid updated state !."; return TrajectoryStateOnSurface(); @@ -28,22 +28,20 @@ TrajectoryStateOnSurface GsfMultiStateUpdator::update(const TrajectoryStateOnSur MultiTrajectoryStateAssembler result; int i = 0; - for (std::vector::const_iterator iter = predictedComponents.begin(); - iter != predictedComponents.end(); iter++) { - TrajectoryStateOnSurface updatedTSOS = KFUpdator().update(*iter, aRecHit); + for (auto const & tsosI : predictedComponents) { + TrajectoryStateOnSurface updatedTSOS = KFUpdator().update(tsosI, aRecHit); if (updatedTSOS.isValid()){ result.addState(TrajectoryStateOnSurface(weights[i], updatedTSOS.localParameters(), updatedTSOS.localError(), updatedTSOS.surface(), &(tsos.globalParameters().magneticField()), - (*iter).surfaceSide() + tsosI.surfaceSide() )); - i++; } else{ - edm::LogError("GsfMultiStateUpdator") << "one of the KF updated state is invalid. skipping."; + edm::LogError("GsfMultiStateUpdator") << "KF updated state " << i << " is invalid. skipping."; } - + ++i; } return result.combinedState(); diff --git a/TrackingTools/GsfTracking/src/GsfPropagatorWithMaterial.cc b/TrackingTools/GsfTracking/src/GsfPropagatorWithMaterial.cc index 9995ad824cd84..0189cf74e246d 100644 --- a/TrackingTools/GsfTracking/src/GsfPropagatorWithMaterial.cc +++ b/TrackingTools/GsfTracking/src/GsfPropagatorWithMaterial.cc @@ -15,7 +15,6 @@ GsfPropagatorWithMaterial::GsfPropagatorWithMaterial (const Propagator& aPropaga theConvolutor(new FullConvolutionWithMaterial(aMEUpdator)), theMaterialLocation(atDestination) { - // if ( propWithPathTimer1==0 ) defineTimer(); } GsfPropagatorWithMaterial::GsfPropagatorWithMaterial (const GsfPropagatorAdapter& aGsfPropagator, @@ -25,34 +24,11 @@ GsfPropagatorWithMaterial::GsfPropagatorWithMaterial (const GsfPropagatorAdapter theConvolutor(aConvolutor.clone()), theMaterialLocation(atDestination) { - // if ( propWithPathTimer1==0 ) defineTimer(); } -// void -// GsfPropagatorWithMaterial::defineTimer() -// { -// if ( propWithPathTimer1==0 ) { -// propWithPathTimer1 = -// &(*TimingReport::current())[string("GsfPropagatorWithMaterial:toPlane")]; -// propWithPathTimer2 = -// &(*TimingReport::current())[string("GsfPropagatorWithMaterial:toCylinder")]; -// static SimpleConfigurable timeConf(false,"GsfPropagatorWithMaterial:activateTiming"); -// if ( timeConf.value() ) { -// propWithPathTimer1->switchCPU(false); -// propWithPathTimer2->switchCPU(false); -// } -// else { -// propWithPathTimer1->switchOn(false); -// propWithPathTimer2->switchOn(false); -// } -// } -// } - std::pair GsfPropagatorWithMaterial::propagateWithPath (const TrajectoryStateOnSurface& tsos, const Plane& plane) const { - // TimeMe t1(*propWithPathTimer1,false); - // // add material before propagation? // TrajectoryStateOnSurface stateAtSource; @@ -75,8 +51,6 @@ GsfPropagatorWithMaterial::propagateWithPath (const TrajectoryStateOnSurface& ts std::pair GsfPropagatorWithMaterial::propagateWithPath (const TrajectoryStateOnSurface& tsos, const Cylinder& cylinder) const { - // TimeMe t2(*propWithPathTimer2,false); - // // add material before propagation? // TrajectoryStateOnSurface stateAtSource; @@ -167,5 +141,3 @@ GsfPropagatorWithMaterial::materialAtSource() const { propagationDirection()==alongMomentum); } -// TimingReport::Item* GsfPropagatorWithMaterial::propWithPathTimer1(0); -// TimingReport::Item* GsfPropagatorWithMaterial::propWithPathTimer2(0); diff --git a/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc b/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc index 7c8662789b69a..7631757f19251 100644 --- a/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc +++ b/TrackingTools/GsfTracking/src/GsfTrajectoryFitter.cc @@ -2,17 +2,15 @@ #include "TrackingTools/TrackFitters/interface/TrajectoryStateWithArbitraryError.h" #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h" -// #include "CommonDet/BasicDet/interface/Det.h" #include "DataFormats/GeometrySurface/interface/BoundPlane.h" #include "TrackingTools/GsfTracking/interface/MultiTrajectoryStateMerger.h" -// #include "Utilities/Notification/interface/Verbose.h" -// #include "Utilities/Notification/interface/TimingReport.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" - #include "DataFormats/TrackerRecHit2D/interface/TkCloner.h" #include "DataFormats/TrackerRecHit2D/interface/BaseTrackerRecHit.h" +#include "DebugHelpers.h" + GsfTrajectoryFitter::GsfTrajectoryFitter(const Propagator& aPropagator, const TrajectoryStateUpdator& aUpdator, @@ -26,8 +24,6 @@ GsfTrajectoryFitter::GsfTrajectoryFitter(const Propagator& aPropagator, theGeometry(detLayerGeometry) { if(!theGeometry) theGeometry = &dummyGeometry; - // static SimpleConfigurable timeConf(false,"GsfTrajectoryFitter:activateTiming"); - // theTiming = timeConf.value(); } GsfTrajectoryFitter::~GsfTrajectoryFitter() { @@ -60,22 +56,13 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, const TrajectoryStateOnSurface& firstPredTsos, fitType) const { - // static TimingReport::Item* propTimer = - // &(*TimingReport::current())[string("GsfTrajectoryFitter:propagation")]; - // propTimer->switchCPU(false); - // if ( !theTiming ) propTimer->switchOn(false); - // static TimingReport::Item* updateTimer = - // &(*TimingReport::current())[string("GsfTrajectoryFitter:update")]; - // updateTimer->switchCPU(false); - // if ( !theTiming ) updateTimer->switchOn(false); - if(hits.empty()) return Trajectory(); Trajectory myTraj(aSeed, propagator()->propagationDirection()); TSOS predTsos(firstPredTsos); if(!predTsos.isValid()) { - edm::LogInfo("GsfTrajectoryFitter") + edm::LogInfo("GsfTrackFitter") << "GsfTrajectoryFitter: predicted tsos of first measurement not valid!"; return Trajectory(); } @@ -87,8 +74,8 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, assert( (!(ihit)->canImproveWithTrack()) | (nullptr!=theHitCloner)); assert( (!(ihit)->canImproveWithTrack()) | (nullptr!=dynamic_cast(ihit.get()))); auto preciseHit = theHitCloner->makeShared(ihit,predTsos); + dump(*preciseHit,1); { - // TimeMe t(*updateTimer,false); currTsos = updator()->update(predTsos, *preciseHit); } if (!predTsos.isValid() || !currTsos.isValid()){ @@ -103,36 +90,22 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, edm::LogError("InvalidState")<<"first invalid hit"; return Trajectory(); } - myTraj.push(TM(predTsos, *hits.begin(),0., theGeometry->idToLayer((*hits.begin())->geographicalId()) )); + myTraj.push(TM(predTsos, hits.front(),0., theGeometry->idToLayer((hits.front())->geographicalId()) )); } + int hitcounter = 1; for(RecHitContainer::const_iterator ihit = hits.begin() + 1; - ihit != hits.end(); ihit++) { + ihit != hits.end(); ihit++) { + ++hitcounter; + // // temporary protection copied from KFTrajectoryFitter. // if ((**ihit).isValid() == false && (**ihit).det() == 0) { - LogDebug("GsfTrajectoryFitter") << " Error: invalid hit with no GeomDet attached .... skipping"; + LogDebug("GsfTrackFitters") << " Error: invalid hit with no GeomDet attached .... skipping"; continue; } - //!!! no invalid hits on cylinders anymore?? - // // - // // check type of surface in case of invalid hit - // // (in this version only propagations to planes are - // // supported for multi trajectory states) - // // - // if ( !(**ihit).isValid() ) { - // const BoundPlane* plane = - // dynamic_cast(&(**ihit).det().surface()); - // // - // // no plane: insert invalid measurement - // // - // if ( plane==0 ) { - // myTraj.push(TM(TrajectoryStateOnSurface(),&(**ihit))); - // continue; - // } - // } { // TimeMe t(*propTimer,false); predTsos = propagator()->propagate(currTsos, @@ -140,13 +113,13 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, } if(!predTsos.isValid()) { if ( myTraj.foundHits()>=3 ) { - edm::LogInfo("GsfTrajectoryFitter") + edm::LogInfo("GsfTrackFitters") << "GsfTrajectoryFitter: predicted tsos not valid! \n" << "Returning trajectory with " << myTraj.foundHits() << " found hits."; return myTraj; } else { - edm::LogInfo("GsfTrajectoryFitter") + edm::LogInfo("GsfTrackFitters") << "GsfTrajectoryFitter: predicted tsos not valid after " << myTraj.foundHits() << " hits, discarding candidate!"; return Trajectory(); @@ -156,28 +129,30 @@ Trajectory GsfTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, if((**ihit).isValid()) { //update - assert( (!(*ihit)->canImproveWithTrack()) | (nullptr!=theHitCloner)); - assert( (!(*ihit)->canImproveWithTrack()) | (nullptr!=dynamic_cast((*ihit).get()))); - auto preciseHit = theHitCloner->makeShared(*ihit,predTsos); - { - // TimeMe t(*updateTimer,false); - currTsos = updator()->update(predTsos, *preciseHit); - } + assert( (!(*ihit)->canImproveWithTrack()) | (nullptr!=theHitCloner)); + assert( (!(*ihit)->canImproveWithTrack()) | (nullptr!=dynamic_cast((*ihit).get()))); + auto preciseHit = theHitCloner->makeShared(*ihit,predTsos); + dump(*preciseHit,hitcounter); + currTsos = updator()->update(predTsos, *preciseHit); if (!predTsos.isValid() || !currTsos.isValid()){ edm::LogError("InvalidState")<<"inside hit"; return Trajectory(); } + auto chi2=estimator()->estimate(predTsos, *preciseHit).second; myTraj.push(TM(predTsos, currTsos, preciseHit, - estimator()->estimate(predTsos, *preciseHit).second, + chi2, theGeometry->idToLayer(preciseHit->geographicalId() ))); + LogDebug("GsfTrackFitters") << "added measurement with chi2 " << chi2; } else { currTsos = predTsos; if (!predTsos.isValid()){ - edm::LogError("InvalidState")<<"inside invalid hit"; - return Trajectory(); + edm::LogError("InvalidState")<<"inside invalid hit"; + return Trajectory(); } myTraj.push(TM(predTsos, *ihit,0., theGeometry->idToLayer( (*ihit)->geographicalId()) )); } + dump(predTsos,"predTsos"); + dump(currTsos,"currTsos"); } return myTraj; } diff --git a/TrackingTools/GsfTracking/src/GsfTrajectorySmoother.cc b/TrackingTools/GsfTracking/src/GsfTrajectorySmoother.cc index 6bfca3a78c00e..07c5226873100 100644 --- a/TrackingTools/GsfTracking/src/GsfTrajectorySmoother.cc +++ b/TrackingTools/GsfTracking/src/GsfTrajectorySmoother.cc @@ -6,6 +6,9 @@ #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "DebugHelpers.h" + + GsfTrajectorySmoother::GsfTrajectorySmoother(const GsfPropagatorWithMaterial& aPropagator, const TrajectoryStateUpdator& aUpdator, const MeasurementEstimator& aEstimator, @@ -37,8 +40,6 @@ GsfTrajectorySmoother::GsfTrajectorySmoother(const GsfPropagatorWithMaterial& aP if(!theGeometry) theGeometry = &dummyGeometry; - // static SimpleConfigurable timeConf(false,"GsfTrajectorySmoother:activateTiming"); - // theTiming = timeConf.value(); } GsfTrajectorySmoother::~GsfTrajectorySmoother() { @@ -54,15 +55,6 @@ GsfTrajectorySmoother::~GsfTrajectorySmoother() { Trajectory GsfTrajectorySmoother::trajectory(const Trajectory& aTraj) const { - // static TimingReport::Item* propTimer = - // &(*TimingReport::current())[string("GsfTrajectorySmoother:propagation")]; - // propTimer->switchCPU(false); - // if ( !theTiming ) propTimer->switchOn(false); - // static TimingReport::Item* updateTimer = - // &(*TimingReport::current())[string("GsfTrajectorySmoother:update")]; - // updateTimer->switchCPU(false); - // if ( !theTiming ) updateTimer->switchOn(false); - if(aTraj.empty()) return Trajectory(); const Propagator* usePropagator = theAlongPropagator; @@ -81,7 +73,7 @@ GsfTrajectorySmoother::trajectory(const Trajectory& aTraj) const { predTsos.rescaleError(theErrorRescaling); if(!predTsos.isValid()) { - edm::LogInfo("GsfTrajectorySmoother") + edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: predicted tsos of last measurement not valid!"; return Trajectory(); } @@ -128,49 +120,30 @@ GsfTrajectorySmoother::trajectory(const Trajectory& aTraj) const { TrajectoryStateCombiner combiner; + int hitcounter = avtm.size()-1; for(std::vector::const_reverse_iterator itm = avtm.rbegin() + 1; itm < avtm.rend() - 1; ++itm) { - { - // TimeMe t(*propTimer,false); - // // - // // check type of surface in case of invalid hit - // // (in this version only propagations to planes are - // // supported for multi trajectory states) - // // - // if ( !(*itm).recHit().isValid() ) { - // const BoundPlane* plane = - // dynamic_cast(&(*itm).recHit().det().surface()); - // // - // // no plane: insert invalid - // if ( plane==0 ) { - // myTraj.push(TM(TrajectoryStateOnSurface(), - // (*itm).recHit()); - // continue; - // } - // } - predTsos = usePropagator->propagate(currTsos, + + predTsos = usePropagator->propagate(currTsos, *(*itm).recHit()->surface()); - } if ( predTsos.isValid() && theConvolutor && theMatBeforeUpdate ) predTsos = (*theConvolutor)(predTsos, usePropagator->propagationDirection()); if(!predTsos.isValid()) { - edm::LogInfo("GsfTrajectorySmoother") << "GsfTrajectorySmoother: predicted tsos not valid!"; + edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: predicted tsos not valid!"; return Trajectory(); } if ( theMerger ) predTsos = theMerger->merge(predTsos); if((*itm).recHit()->isValid()) { //update - { - // TimeMe t(*updateTimer,false); currTsos = updator()->update(predTsos, *(*itm).recHit()); - } + if ( currTsos.isValid() && theConvolutor && !theMatBeforeUpdate ) currTsos = (*theConvolutor)(currTsos, usePropagator->propagationDirection()); if(!currTsos.isValid()) { - edm::LogInfo("GsfTrajectorySmoother") + edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: tsos not valid after update / material effects!"; return Trajectory(); } @@ -180,7 +153,7 @@ GsfTrajectorySmoother::trajectory(const Trajectory& aTraj) const { //3: combine bwd-prediction with fwd-filter TSOS combTsos = combiner(predTsos, (*itm).forwardPredictedState()); if(!combTsos.isValid()) { - LogDebug("GsfTrajectorySmoother") << + LogDebug("GsfTrackFitters") << "KFTrajectorySmoother: combined tsos not valid!\n"<< "pred Tsos pos: "<estimate(combTsos, *(*itm).recHit()).second; myTraj.push(TM((*itm).forwardPredictedState(), predTsos, smooTsos, (*itm).recHit(), - estimator()->estimate(combTsos, *(*itm).recHit()).second, + chi2, theGeometry->idToLayer((*itm).recHit()->geographicalId() ) ), (*itm).estimate()); + LogDebug("GsfTrackFitters") << "added measurement #" << hitcounter-- <<" with chi2 " << chi2; + dump(predTsos,"predTsos"); + dump(smooTsos,"smooTsos"); + } else { currTsos = predTsos; TSOS combTsos = combiner(predTsos, (*itm).forwardPredictedState()); if(!combTsos.isValid()) { - LogDebug("GsfTrajectorySmoother") << + LogDebug("GsfTrackFitters") << "KFTrajectorySmoother: combined tsos not valid!"; return Trajectory(); } @@ -231,36 +208,35 @@ GsfTrajectorySmoother::trajectory(const Trajectory& aTraj) const { (*itm).recHit(), 0., theGeometry->idToLayer((*itm).recHit()->geographicalId()) )); + LogDebug("GsfTrackFitters") << "added invalid measurement #" << hitcounter--; + dump(predTsos,"predTsos"); + dump(combTsos,"smooTsos"); } if ( theMerger ) currTsos = theMerger->merge(currTsos); + + dump(currTsos,"currTsos"); } //last smoothed tm is last filtered - { - // TimeMe t(*propTimer,false); predTsos = usePropagator->propagate(currTsos, *avtm.front().recHit()->surface()); - } if ( predTsos.isValid() && theConvolutor && theMatBeforeUpdate ) predTsos = (*theConvolutor)(predTsos, usePropagator->propagationDirection()); if(!predTsos.isValid()) { - edm::LogInfo("GsfTrajectorySmoother") << "GsfTrajectorySmoother: predicted tsos not valid!"; + edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: predicted tsos not valid!"; return Trajectory(); } if ( theMerger ) predTsos = theMerger->merge(predTsos); if(avtm.front().recHit()->isValid()) { //update - { - // TimeMe t(*updateTimer,false); - currTsos = updator()->update(predTsos, *avtm.front().recHit()); - } + currTsos = updator()->update(predTsos, *avtm.front().recHit()); if ( currTsos.isValid() && theConvolutor && !theMatBeforeUpdate ) currTsos = (*theConvolutor)(currTsos, usePropagator->propagationDirection()); if(!currTsos.isValid()) { - edm::LogInfo("GsfTrajectorySmoother") + edm::LogInfo("GsfTrackFitters") << "GsfTrajectorySmoother: tsos not valid after update / material effects!"; return Trajectory(); } @@ -270,14 +246,18 @@ GsfTrajectorySmoother::trajectory(const Trajectory& aTraj) const { return Trajectory(); } + auto chi2 = estimator()->estimate(predTsos, *avtm.front().recHit()).second; myTraj.push(TM(avtm.front().forwardPredictedState(), predTsos, currTsos, avtm.front().recHit(), - estimator()->estimate(predTsos, *avtm.front().recHit()).second, + chi2, theGeometry->idToLayer(avtm.front().recHit()->geographicalId() )), avtm.front().estimate()); - //estimator()->estimate(predTsos, avtm.front().recHit())); + LogDebug("GsfTrackFitters") << "added measurement #" << hitcounter-- <<" with chi2 " << chi2; + dump(predTsos,"predTsos"); + dump(currTsos,"smooTsos"); + } else { if (!avtm.front().forwardPredictedState().isValid()){ @@ -288,6 +268,7 @@ GsfTrajectorySmoother::trajectory(const Trajectory& aTraj) const { avtm.front().recHit(), 0., theGeometry->idToLayer(avtm.front().recHit()->geographicalId()) )); + LogDebug("GsfTrackFitters") << "added invalid measurement #" << hitcounter--; } return myTraj; diff --git a/TrackingTools/GsfTracking/src/MultiTrajectoryStateMerger.cc b/TrackingTools/GsfTracking/src/MultiTrajectoryStateMerger.cc index 24652c09661ae..e79353a437763 100644 --- a/TrackingTools/GsfTracking/src/MultiTrajectoryStateMerger.cc +++ b/TrackingTools/GsfTracking/src/MultiTrajectoryStateMerger.cc @@ -1,5 +1,4 @@ #include "TrackingTools/GsfTracking/interface/MultiTrajectoryStateMerger.h" - #include "TrackingTools/GsfTracking/interface/TsosGaussianStateConversions.h" TrajectoryStateOnSurface diff --git a/TrackingTools/GsfTracking/src/PosteriorWeightsCalculator.cc b/TrackingTools/GsfTracking/src/PosteriorWeightsCalculator.cc index d57f9d86374f8..59bcbacfce910 100644 --- a/TrackingTools/GsfTracking/src/PosteriorWeightsCalculator.cc +++ b/TrackingTools/GsfTracking/src/PosteriorWeightsCalculator.cc @@ -3,6 +3,8 @@ #include "TrackingTools/PatternTools/interface/MeasurementExtractor.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "DataFormats/TrackingRecHit/interface/KfComponentsHolder.h" +#include "DataFormats/Math/interface/invertPosDefMatrix.h" +#include "DataFormats/Math/interface/ProjectMatrix.h" #include @@ -39,7 +41,6 @@ std::vector PosteriorWeightsCalculator::weights(const TrackingRecHit& re VecD r, rMeas; SMatDD V(SMatrixNoInit{}), R(SMatrixNoInit{}); - AlgebraicVector5 x; ProjectMatrix p; // // calculate chi2 and determinant / component and find @@ -49,7 +50,7 @@ std::vector PosteriorWeightsCalculator::weights(const TrackingRecHit& re for ( unsigned int i=0; i(&r, &V, &p, &rMeas, &R, x, predictedComponents[i].localError().matrix()); recHit.getKfComponents(holder); @@ -64,21 +65,23 @@ std::vector PosteriorWeightsCalculator::weights(const TrackingRecHit& re } detRs.push_back(detR); - int ierr = ! R.Invert(); // if (ierr != 0) throw exception; - if ( ierr!=0 ) { + bool ok = invertPosDefMatrix(R); + if ( !ok ) { edm::LogError("PosteriorWeightsCalculator") - << "PosteriorWeightsCalculator: inversion failed, ierr = " << ierr; + << "PosteriorWeightsCalculator: inversion failed"; return std::vector(); } double chi2 = ROOT::Math::Similarity(r,R); chi2s.push_back(chi2); if ( chi2(); } + // // calculate weights (extracting a common factor // exp(-0.5*chi2Min) to avoid numerical problems @@ -96,16 +99,17 @@ std::vector PosteriorWeightsCalculator::weights(const TrackingRecHit& re // Calculation of (non-normalised) weight. Common factors exp(-chi2Norm/2.) and // 1./sqrt(2*pi*recHit.dimension()) have been omitted // - tempWeight = priorWeight * sqrt(1./detRs[i]) * exp(-0.5 * chi2); + tempWeight = priorWeight * std::sqrt(1./detRs[i]) * std::exp(-0.5 * chi2); + } + else { + LogDebug("GsfTrackFitters") << "PosteriorWeightsCalculator: detR < FLT_MIN !!"; } - // else { - // edm::LogInfo("PosteriorWeightsCalculator") << "PosteriorWeightsCalculator: detR < FLT_MIN !!"; - // } weights.push_back(tempWeight); sumWeights += tempWeight; } + if ( sumWeights(); } @@ -114,10 +118,7 @@ std::vector PosteriorWeightsCalculator::weights(const TrackingRecHit& re edm::LogError("PosteriorWeightsCalculator") << "Problem in vector sizes (2)"; return std::vector(); } - for (std::vector::iterator iter = weights.begin(); - iter != weights.end(); iter++) { - (*iter) /= sumWeights; - } - + sumWeights = 1./sumWeights; + for (auto & w : weights) w *= sumWeights; return weights; } diff --git a/TrackingTools/GsfTracking/src/TsosGaussianStateConversions.cc b/TrackingTools/GsfTracking/src/TsosGaussianStateConversions.cc index d11adafc3185f..f1c5373f8131f 100644 --- a/TrackingTools/GsfTracking/src/TsosGaussianStateConversions.cc +++ b/TrackingTools/GsfTracking/src/TsosGaussianStateConversions.cc @@ -12,16 +12,15 @@ namespace GaussianStateConversions { { if ( !tsos.isValid() ) return MultiGaussianState<5>(); - typedef boost::shared_ptr< SingleGaussianState<5> > SingleStatePtr; - const std::vector& components = tsos.components(); + using SingleStatePtr = boost::shared_ptr>; + auto const & components = tsos.components(); MultiGaussianState<5>::SingleStateContainer singleStates; singleStates.reserve(components.size()); - for ( std::vector::const_iterator ic=components.begin(); - ic!=components.end(); ic ++ ) { - if ( ic->isValid() ) { - SingleStatePtr sgs(new SingleGaussianState<5>(ic->localParameters().vector(), - ic->localError().matrix(), - ic->weight())); + for (auto const & ic : components) { + if ( ic.isValid() ) { + SingleStatePtr sgs(new SingleGaussianState<5>(ic.localParameters().vector(), + ic.localError().matrix(), + ic.weight())); singleStates.push_back(sgs); } } @@ -32,23 +31,21 @@ namespace GaussianStateConversions { const TrajectoryStateOnSurface refTsos) { if ( multiState.components().empty() ) return TrajectoryStateOnSurface(); - const Surface& surface = refTsos.surface(); + const Surface & surface = refTsos.surface(); SurfaceSide side = refTsos.surfaceSide(); const MagneticField* field = refTsos.magneticField(); - TrajectoryStateOnSurface refTsos1 = refTsos.components().front(); - double pzSign = refTsos1.localParameters().pzSign(); + auto const & refTsos1 = refTsos.components().front(); + auto pzSign = refTsos1.localParameters().pzSign(); bool charged = refTsos1.charge()!=0; - const MultiGaussianState<5>::SingleStateContainer& singleStates = - multiState.components(); + auto const & singleStates = multiState.components(); std::vector components; components.reserve(singleStates.size()); - for ( MultiGaussianState<5>::SingleStateContainer::const_iterator ic=singleStates.begin(); - ic!=singleStates.end(); ic++ ) { - components.push_back(TrajectoryStateOnSurface((**ic).weight(), LocalTrajectoryParameters((**ic).mean(), - pzSign,charged), - LocalTrajectoryError((**ic).covariance()), - surface,field,side)); + for ( auto const & ic : singleStates ) { + components.emplace_back((*ic).weight(), + LocalTrajectoryParameters((*ic).mean(), pzSign,charged), + LocalTrajectoryError((*ic).covariance()), + surface,field,side); } return TrajectoryStateOnSurface((BasicTrajectoryState*)new BasicMultiTrajectoryState(components)); } diff --git a/TrackingTools/KalmanUpdators/BuildFile.xml b/TrackingTools/KalmanUpdators/BuildFile.xml index 3e9d8f8035384..4fed031aa2979 100644 --- a/TrackingTools/KalmanUpdators/BuildFile.xml +++ b/TrackingTools/KalmanUpdators/BuildFile.xml @@ -1,6 +1,4 @@ - -