diff --git a/DPGAnalysis/SiStripTools/plugins/TrackerDpgAnalysis.cc b/DPGAnalysis/SiStripTools/plugins/TrackerDpgAnalysis.cc index b5033d206885a..1a10c33e8a3b6 100644 --- a/DPGAnalysis/SiStripTools/plugins/TrackerDpgAnalysis.cc +++ b/DPGAnalysis/SiStripTools/plugins/TrackerDpgAnalysis.cc @@ -89,7 +89,6 @@ #include #include #include -#include #include #include #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" @@ -1104,7 +1103,7 @@ std::vector TrackerDpgAnalysis::onTrackAngles(edm::Handle >& collection,const TransientTrackingRecHit* hit , double tla) { if(!hit) return; - const TSiTrackerMultiRecHit* multihit=dynamic_cast(hit); + const SiTrackerMultiRecHit* multihit=dynamic_cast(hit); const SiStripRecHit2D* singlehit=dynamic_cast(hit); const SiStripRecHit1D* hit1d=dynamic_cast(hit); if(hit1d) { //...->33X diff --git a/DataFormats/TrackerRecHit2D/interface/SiTrackerMultiRecHit.h b/DataFormats/TrackerRecHit2D/interface/SiTrackerMultiRecHit.h index 2b0046ebb33f0..5e64fde8303e8 100644 --- a/DataFormats/TrackerRecHit2D/interface/SiTrackerMultiRecHit.h +++ b/DataFormats/TrackerRecHit2D/interface/SiTrackerMultiRecHit.h @@ -16,14 +16,18 @@ class SiTrackerMultiRecHit : public BaseTrackerRecHit typedef BaseTrackerRecHit Base; SiTrackerMultiRecHit(): theHits(), - theWeights(){} + theWeights(), + annealing_(0){} virtual ~SiTrackerMultiRecHit(){} SiTrackerMultiRecHit(const LocalPoint&, const LocalError&, GeomDet const & idet, - const std::vector< std::pair >&); + const std::vector< std::pair >&, double); virtual SiTrackerMultiRecHit* clone() const {return new SiTrackerMultiRecHit(*this);} +#ifdef NO_DICT + virtual RecHitPointer cloneSH() const { return std::make_shared(*this);} +#endif virtual int dimension() const {return 2;} virtual void getKfComponents( KfComponentsHolder & holder ) const { getKfComponents2D(holder); } @@ -32,11 +36,13 @@ class SiTrackerMultiRecHit : public BaseTrackerRecHit // used by trackMerger (to be improved) virtual OmniClusterRef const & firstClusterRef() const { return static_cast(&theHits.front())->firstClusterRef();} - - //vector of component rechits + /// Access to component RecHits (if any) virtual std::vector recHits() const; - +// virtual void recHitsV(std::vector & ) const; + + /// Non-const access to component RecHits (if any) virtual std::vector recHits() ; +// virtual void recHitsV(std::vector & ); //vector of weights std::vector const & weights() const {return theWeights;} @@ -45,14 +51,18 @@ class SiTrackerMultiRecHit : public BaseTrackerRecHit //returns the weight for the i component float weight(unsigned int i) const {return theWeights[i];} float & weight(unsigned int i) {return theWeights[i];} + + //get the annealing + virtual double getAnnealingFactor() const { return annealing_; } bool sharesInput(const TrackingRecHit* other, SharedInputType what) const; + private: edm::OwnVector theHits; std::vector theWeights; - + double annealing_; }; diff --git a/DataFormats/TrackerRecHit2D/src/SiTrackerMultiRecHit.cc b/DataFormats/TrackerRecHit2D/src/SiTrackerMultiRecHit.cc index f1ece481d14fa..6c0d7afdff491 100644 --- a/DataFormats/TrackerRecHit2D/src/SiTrackerMultiRecHit.cc +++ b/DataFormats/TrackerRecHit2D/src/SiTrackerMultiRecHit.cc @@ -5,13 +5,14 @@ using namespace std; using namespace edm; SiTrackerMultiRecHit::SiTrackerMultiRecHit(const LocalPoint& pos, const LocalError& err, GeomDet const & idet, - const std::vector< std::pair >& aHitMap): + const std::vector< std::pair >& aHitMap, double annealing): BaseTrackerRecHit(pos,err, idet,trackerHitRTTI::multi) { for(std::vector >::const_iterator ihit = aHitMap.begin(); ihit != aHitMap.end(); ihit++){ theHits.push_back(ihit->first->clone()); theWeights.push_back(ihit->second); } + annealing_ = annealing; } @@ -68,3 +69,13 @@ vector SiTrackerMultiRecHit::recHits() { // } return theHits.data(); } + +/* +void TrackingRecHit::recHitsV(std::vector & v) const { + v = recHits(); +} + +void TrackingRecHit::recHitsV(std::vector & v) { + v = recHits(); +} +*/ diff --git a/DataFormats/TrackerRecHit2D/src/classes_def.xml b/DataFormats/TrackerRecHit2D/src/classes_def.xml index 382a328eb9d73..e6f8351caf0aa 100755 --- a/DataFormats/TrackerRecHit2D/src/classes_def.xml +++ b/DataFormats/TrackerRecHit2D/src/classes_def.xml @@ -48,7 +48,8 @@ - + + diff --git a/RecoTracker/SiTrackerMRHTools/BuildFile.xml b/RecoTracker/SiTrackerMRHTools/BuildFile.xml new file mode 100644 index 0000000000000..82cec149978e0 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/BuildFile.xml @@ -0,0 +1,21 @@ + + + + + + + + + + + + + + + + + + + + + diff --git a/RecoTracker/SiTrackerMRHTools/interface/GenericProjectedRecHit2D.h b/RecoTracker/SiTrackerMRHTools/interface/GenericProjectedRecHit2D.h new file mode 100644 index 0000000000000..6e25e52e252c8 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/interface/GenericProjectedRecHit2D.h @@ -0,0 +1,85 @@ +#ifndef RECOTRACKER_TRANSIENTRACKINGRECHIT_GenericProjectedRecHit2D_H +#define RECOTRACKER_TRANSIENTRACKINGRECHIT_GenericProjectedRecHit2D_H + +#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h" +#include "TrackingTools/TransientTrackingRecHit/interface/HelpertRecHit2DLocalPos.h" +#include "TrackingTools/KalmanUpdators/interface/TrackingRecHitPropagator.h" + +class GenericProjectedRecHit2D : public TransientTrackingRecHit { +public: + + GenericProjectedRecHit2D( const LocalPoint& pos, const LocalError& err, + const GeomDet* det, const GeomDet* originaldet, + const TransientTrackingRecHit::ConstRecHitPointer originalHit, + const TrackingRecHitPropagator* propagator); + + virtual AlgebraicSymMatrix parametersError() const { + return HelpertRecHit2DLocalPos().parError( localPositionError(), *det()); + } + + //virtual ~GenericProjectedRecHit2D(){delete theOriginalTransientHit;} + + virtual AlgebraicVector parameters() const ; + + virtual LocalPoint localPosition() const {return theLp;} + + virtual LocalError localPositionError() const {return theLe;} + + virtual AlgebraicMatrix projectionMatrix() const {return theProjectionMatrix;} + + virtual DetId geographicalId() const {return det() ? det()->geographicalId() : DetId();} + + virtual int dimension() const {return theDimension;} + + //this hit lays on the original surface, NOT on the projection surface + virtual const TrackingRecHit * hit() const { return theOriginalTransientHit->hit(); } + + virtual TrackingRecHit * cloneHit() const { return theOriginalTransientHit->cloneHit(); } + + virtual bool isValid() const{return true;} + + virtual std::vector recHits() const { + //return theOriginalTransientHit->hit()->recHits(); + return std::vector(); + } + + virtual std::vector recHits() { + //should it do something different? + return std::vector(); + } + + const TrackingRecHitPropagator* propagator() const {return thePropagator;} + + virtual bool canImproveWithTrack() const {return true;} + + const GeomDet* originalDet() const {return theOriginalDet;} + + static RecHitPointer build( const LocalPoint& pos, const LocalError& err, + const GeomDet* det, const GeomDet* originaldet, + const TransientTrackingRecHit::ConstRecHitPointer originalHit, + const TrackingRecHitPropagator* propagator) { + return RecHitPointer( new GenericProjectedRecHit2D( pos, err, det, originaldet, originalHit, propagator)); + } + + RecHitPointer clone( const TrajectoryStateOnSurface& ts, const TransientTrackingRecHitBuilder*) const; + +private: + + const GeomDet* theOriginalDet; + TransientTrackingRecHit::ConstRecHitPointer theOriginalTransientHit; + LocalPoint theLp; + LocalError theLe; + AlgebraicMatrix theProjectionMatrix; + const TrackingRecHitPropagator* thePropagator; + //const TrackingRecHit* theOriginalHit; + int theDimension; + + virtual GenericProjectedRecHit2D* clone() const { + return new GenericProjectedRecHit2D(*this); + } + +}; + + + +#endif diff --git a/RecoTracker/SiTrackerMRHTools/interface/GroupedDAFHitCollector.h b/RecoTracker/SiTrackerMRHTools/interface/GroupedDAFHitCollector.h new file mode 100644 index 0000000000000..4dd68e0aadcea --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/interface/GroupedDAFHitCollector.h @@ -0,0 +1,56 @@ +/** \class GroupedDAFHitCollector + * Returns a collection of SiTrackerMultiRecHits and InvalidRecHits given a Trajectory. + * Builds a SiTrackerMultiRecHit for each detGroup + * (i.e. a group of detectors mutually exclusive for the track's crossing point) + * + * \author tropiano, genta + * \review in May 2014 by brondolin + */ + +#ifndef SiTrackerMRHTools_GroupedDAFHitCollector_h +#define SiTrackerMRHTools_GroupedDAFHitCollector_h + +#include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h" +#include "RecoTracker/SiTrackerMRHTools/interface/MultiRecHitCollector.h" +#include + +class Propagator; +class MeasurementEstimator; +class MeasurementTracker; +class SiTrackerMultiRecHitUpdator; + +class GroupedDAFHitCollector :public MultiRecHitCollector { + +public: + explicit GroupedDAFHitCollector(const MeasurementTracker* measurementTracker, + const SiTrackerMultiRecHitUpdator* updator, + const MeasurementEstimator* est, + const Propagator* propagator, + const Propagator* reversePropagator, bool debug): + MultiRecHitCollector(measurementTracker), theUpdator(updator), + theEstimator(est), thePropagator(propagator), theReversePropagator(reversePropagator), debug_(debug){} + + + virtual ~GroupedDAFHitCollector(){} + + virtual std::vector recHits(const Trajectory&, + const MeasurementTrackerEvent *theMT) const; + + const SiTrackerMultiRecHitUpdator* getUpdator() const {return theUpdator;} + const MeasurementEstimator* getEstimator() const {return theEstimator;} + const Propagator* getPropagator() const {return thePropagator;} + const Propagator* getReversePropagator() const {return theReversePropagator;} + +private: + void buildMultiRecHits(const std::vector& measgroup, + std::vector& result) const; + + const SiTrackerMultiRecHitUpdator* theUpdator; + const MeasurementEstimator* theEstimator; + const Propagator* thePropagator; + const Propagator* theReversePropagator; + const bool debug_; +}; + + +#endif diff --git a/RecoTracker/SiTrackerMRHTools/interface/MeasurementByLayerGrouper.h b/RecoTracker/SiTrackerMRHTools/interface/MeasurementByLayerGrouper.h new file mode 100644 index 0000000000000..50559453ada36 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/interface/MeasurementByLayerGrouper.h @@ -0,0 +1,39 @@ +#ifndef SiTrackerMRHTools_MeasurementByLayerGrouper_H +#define SiTrackerMRHTools_MeasurementByLayerGrouper_H + +class DetLayer; +class TrajectoryMeasurement; +class GeometricSearchTracker; + +#include +#include + +//groups the TrajectoryMeasurements on a layer by layer basis + +class MeasurementByLayerGrouper { + +private: + + typedef TrajectoryMeasurement TM; + const GeometricSearchTracker* theGeomSearch; + + const DetLayer* getDetLayer(const TM& tm) const; + +public: + + explicit MeasurementByLayerGrouper(const GeometricSearchTracker* search = 0):theGeomSearch(search){}; + + std::vector > > operator()(const std::vector&) const; + + +//to be ported later if needed +/* + vector + operator()(const vector > >&) const; + + vector > > > + operator()(const map >&) const; +*/ + +}; +#endif diff --git a/RecoTracker/SiTrackerMRHTools/interface/MultiRecHitCollector.h b/RecoTracker/SiTrackerMRHTools/interface/MultiRecHitCollector.h new file mode 100644 index 0000000000000..3f6072fe939b3 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/interface/MultiRecHitCollector.h @@ -0,0 +1,26 @@ +#ifndef SiTrackerMRHTools_MultiRecHitCollector_h +#define SiTrackerMRHTools_MultiRecHitCollector_h + +#include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h" +#include + +class Trajectory; +class TrajectoryMeasurement; + +class MultiRecHitCollector { + + public: + MultiRecHitCollector(const MeasurementTracker* meas): theMeasurementTracker(meas){} + + virtual std::vector recHits(const Trajectory&, const MeasurementTrackerEvent *theMTE) const = 0; + + const MeasurementTracker* getMeasurementTracker() const {return theMeasurementTracker;} + + + private: + const MeasurementTracker* theMeasurementTracker; + +}; + +#endif + diff --git a/RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdator.h b/RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdator.h new file mode 100644 index 0000000000000..b4bad7fe1c183 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdator.h @@ -0,0 +1,87 @@ +/** \class SiTrackerMultiRecHitUpdator + * Builds a SiTrackerMultiRecHit out of a vector of TrackingRecHit + * or updates an existing SiTrackerMultiRecHit given a tsos. + * + * \author tropiano, genta + * \review in May 2014 by brondolin + */ + +#ifndef SiTrackerMultiRecHitUpdator_h +#define SiTrackerMultiRecHitUpdator_h + +#include "DataFormats/GeometryVector/interface/LocalPoint.h" +#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h" +#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h" +#include "RecoTracker/TransientTrackingRecHit/interface/TkClonerImpl.h" + +#include + +class SiTrackerMultiRecHit; +class TrajectoryStateOnSurface; +class TrackingRecHit; +class TransientTrackingRecHitBuilder; +class LocalError; +class TrackingRecHitPropagator; + + +class SiTrackerMultiRecHitUpdator{ + +public: + + typedef std::pair LocalParameters; + SiTrackerMultiRecHitUpdator(const TransientTrackingRecHitBuilder* builder, + const TrackingRecHitPropagator* hitpropagator, + const float Chi2Cut, + const std::vector& anAnnealingProgram, bool debug); + virtual ~SiTrackerMultiRecHitUpdator(){}; + + //calls the update method in order to build a SiTrackerMultiRecHit + virtual TransientTrackingRecHit::RecHitPointer buildMultiRecHit(const std::vector& rhv, + const TrajectoryStateOnSurface& tsos, + float annealing=1.) const; + + //updates an existing SiTrackerMultiRecHit + //in case a different kind of rechit is passed it returns clone(tsos) + virtual TransientTrackingRecHit::RecHitPointer update( TransientTrackingRecHit::ConstRecHitPointer original, + const TrajectoryStateOnSurface& tsos, + double annealing=1.) const; + + //returns a SiTrackerMultiRecHit out of the transient components + TransientTrackingRecHit::RecHitPointer update( TransientTrackingRecHit::ConstRecHitContainer& tcomponents, + const TrajectoryStateOnSurface& tsos, + double annealing=1.) const; + + //computes weights or the cut-off value (depending on CutWeight variable) + double ComputeWeight(const TrajectoryStateOnSurface& tsos, const TransientTrackingRecHit& aRecHit, + bool CutWeight, double annealing=1.) const; + template double ComputeWeight(const TrajectoryStateOnSurface& tsos, + const TransientTrackingRecHit& aRecHit, + bool CutWeight, double annealing=1.) const; + + //computes parameters for 2 dim hits + std::pair ComputeParameters2dim(const TrajectoryStateOnSurface& tsos, + const TransientTrackingRecHit& aRecHit) const; + template std::pair ComputeParameters2dim (const TrajectoryStateOnSurface& tsos, + const TransientTrackingRecHit& aRecHit ) const; + + const std::vector& annealingProgram() const {return theAnnealingProgram;} + const std::vector& getAnnealingProgram() const {return theAnnealingProgram;} + + +private: + //Obsolete methods + //LocalError calcParametersError(TransientTrackingRecHit::ConstRecHitContainer& map) const; + //LocalPoint calcParameters(TransientTrackingRecHit::ConstRecHitContainer& map, const LocalError& er) const; + + LocalParameters calcParameters(const TrajectoryStateOnSurface& tsos, + std::vector >& aHitMap) const; + + const TransientTrackingRecHitBuilder* theBuilder; + const TrackingRecHitPropagator* theHitPropagator; + double theChi2Cut; + const std::vector theAnnealingProgram; + TkClonerImpl theHitCloner; + bool debug_; + +}; +#endif diff --git a/RecoTracker/SiTrackerMRHTools/interface/SimpleDAFHitCollector.h b/RecoTracker/SiTrackerMRHTools/interface/SimpleDAFHitCollector.h new file mode 100644 index 0000000000000..4caffea12864e --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/interface/SimpleDAFHitCollector.h @@ -0,0 +1,48 @@ +#ifndef SiTrackerMRHTools_SimpleDAFHitCollector_h +#define SiTrackerMRHTools_SimpleDAFHitCollector_h +#include "RecoTracker/SiTrackerMRHTools/interface/MultiRecHitCollector.h" +#include + +class Propagator; +class MeasurementEstimator; +class SiTrackerMultiRecHitUpdator; + +class SimpleDAFHitCollector :public MultiRecHitCollector { + public: + explicit SimpleDAFHitCollector(const MeasurementTracker* measurementTracker, + const SiTrackerMultiRecHitUpdator* updator, + const MeasurementEstimator* est, + const Propagator* propagator, bool debug + ):MultiRecHitCollector(measurementTracker), theUpdator(updator), theEstimator(est), thePropagator(propagator), debug_(debug){} + + + virtual ~SimpleDAFHitCollector(){} + + //given a trajectory it returns a collection + //of SiTrackerMultiRecHits and InvalidTransientRecHits. + //For each measurement in the trajectory, measurements are looked for according to the + //MeasurementDet::fastMeasurements method only in the detector where the original measurement lays. + //If measurements are found a SiTrackerMultiRecHit is built. + //All the components will lay on the same detector + + virtual std::vector recHits(const Trajectory&, const MeasurementTrackerEvent *theMTE) const; + + const SiTrackerMultiRecHitUpdator* getUpdator() const {return theUpdator;} + const MeasurementEstimator* getEstimator() const {return theEstimator;} + const Propagator* getPropagator() const {return thePropagator;} + + private: + //TransientTrackingRecHit::ConstRecHitContainer buildMultiRecHits(const std::vector& measgroup) const; + void buildMultiRecHits(const std::vector& measgroup, std::vector& result) const; + + private: + const SiTrackerMultiRecHitUpdator* theUpdator; + const MeasurementEstimator* theEstimator; + //this actually is not used in the fastMeasurement method + const Propagator* thePropagator; + const bool debug_; + +}; + + +#endif diff --git a/RecoTracker/SiTrackerMRHTools/plugins/BuildFile.xml b/RecoTracker/SiTrackerMRHTools/plugins/BuildFile.xml new file mode 100644 index 0000000000000..b126ca46b6392 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/plugins/BuildFile.xml @@ -0,0 +1,6 @@ + + + + + + diff --git a/RecoTracker/SiTrackerMRHTools/plugins/MultiRecHitCollectorESProducer.cc b/RecoTracker/SiTrackerMRHTools/plugins/MultiRecHitCollectorESProducer.cc new file mode 100644 index 0000000000000..c43dc61c6b9da --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/plugins/MultiRecHitCollectorESProducer.cc @@ -0,0 +1,75 @@ +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + + +#include "RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdator.h" +#include "RecoTracker/SiTrackerMRHTools/plugins/MultiRecHitCollectorESProducer.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h" +#include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h" +#include "RecoTracker/SiTrackerMRHTools/interface/GroupedDAFHitCollector.h" +#include "RecoTracker/SiTrackerMRHTools/interface/SimpleDAFHitCollector.h" + +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/ModuleFactory.h" +#include "FWCore/Framework/interface/ESProducer.h" + +#include +#include + +using namespace edm; + +MultiRecHitCollectorESProducer::MultiRecHitCollectorESProducer(const edm::ParameterSet& iConfig) +{ + std::string myname = iConfig.getParameter("ComponentName"); + setConf(iConfig); + setWhatProduced(this,myname); +} + +MultiRecHitCollectorESProducer::~MultiRecHitCollectorESProducer() {} + +boost::shared_ptr +MultiRecHitCollectorESProducer::produce(const MultiRecHitRecord& iRecord){ + std::string mode = "Grouped"; + if (conf_.getParameter("Mode")=="Simple") mode = "Simple"; + + std::string mrhupdator = conf_.getParameter("MultiRecHitUpdator"); + std::string propagatorAlongName = conf_.getParameter("propagatorAlong"); + std::string estimatorName = conf_.getParameter("estimator"); + std::string measurementTrackerName = conf_.getParameter("MeasurementTrackerName"); + bool debug = conf_.getParameter("Debug"); + + + ESHandle mrhuhandle; + iRecord.get(mrhupdator, mrhuhandle); + ESHandle propagatorhandle; + iRecord.getRecord().getRecord().get(propagatorAlongName, propagatorhandle); + ESHandle estimatorhandle; + iRecord.getRecord().getRecord().get(estimatorName, estimatorhandle); + ESHandle measurementhandle; + iRecord.getRecord().get(measurementTrackerName, measurementhandle); + + if (mode == "Grouped"){ + std::string propagatorOppositeName = conf_.getParameter("propagatorOpposite"); + ESHandle propagatorOppositehandle; + iRecord.getRecord().getRecord().get(propagatorOppositeName, propagatorOppositehandle); + collector_ = boost::shared_ptr(new GroupedDAFHitCollector(measurementhandle.product(), + mrhuhandle.product(), + estimatorhandle.product(), + propagatorhandle.product(), + propagatorOppositehandle.product(), debug)); + } + else { + collector_ = boost::shared_ptr(new SimpleDAFHitCollector(measurementhandle.product(), + mrhuhandle.product(), + estimatorhandle.product(), + propagatorhandle.product(), debug)); + } + + return collector_; + +} + + diff --git a/RecoTracker/SiTrackerMRHTools/plugins/MultiRecHitCollectorESProducer.h b/RecoTracker/SiTrackerMRHTools/plugins/MultiRecHitCollectorESProducer.h new file mode 100644 index 0000000000000..eb847da40d24a --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/plugins/MultiRecHitCollectorESProducer.h @@ -0,0 +1,39 @@ +#ifndef RecoLocalTracker_ESProducers_MultiRecHitCollectorESProducer_h +#define RecoLocalTracker_ESProducers_ESProducers_MultiRecHitCollectorESProducer_h + +#include "FWCore/Framework/interface/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/InputTag.h" + +#include "FWCore/Framework/interface/ESProducer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "RecoTracker/Record/interface/MultiRecHitRecord.h" +#include "RecoTracker/SiTrackerMRHTools/interface/MultiRecHitCollector.h" + +#include + +class MultiRecHitCollectorESProducer: public edm::ESProducer{ + public: + MultiRecHitCollectorESProducer(const edm::ParameterSet& iConfig); + virtual ~MultiRecHitCollectorESProducer(); + boost::shared_ptr produce(const MultiRecHitRecord &); + + // Set parameter set + void setConf(const edm::ParameterSet& conf){ conf_ = conf; } + + private: + boost::shared_ptr collector_; + edm::ParameterSet conf_; + +}; + + +#endif + + + + diff --git a/RecoTracker/SiTrackerMRHTools/plugins/SiTrackerMultiRecHitUpdatorESProducer.cc b/RecoTracker/SiTrackerMRHTools/plugins/SiTrackerMultiRecHitUpdatorESProducer.cc new file mode 100644 index 0000000000000..93346053e1b13 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/plugins/SiTrackerMultiRecHitUpdatorESProducer.cc @@ -0,0 +1,43 @@ +#include "RecoTracker/SiTrackerMRHTools/plugins/SiTrackerMultiRecHitUpdatorESProducer.h" +#include "TrackingTools/Records/interface/TransientRecHitRecord.h" +#include "RecoTracker/Record/interface/CkfComponentsRecord.h" + +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/ModuleFactory.h" +#include "FWCore/Framework/interface/ESProducer.h" + +#include +#include + +using namespace edm; + +SiTrackerMultiRecHitUpdatorESProducer::SiTrackerMultiRecHitUpdatorESProducer(const edm::ParameterSet & p) +{ + std::string myname = p.getParameter("ComponentName"); + pset_ = p; + setWhatProduced(this,myname); +} + +SiTrackerMultiRecHitUpdatorESProducer::~SiTrackerMultiRecHitUpdatorESProducer() {} + +boost::shared_ptr +SiTrackerMultiRecHitUpdatorESProducer::produce(const MultiRecHitRecord & iRecord){ + std::vector annealingProgram = pset_.getParameter >("AnnealingProgram"); + float Chi2Cut=pset_.getParameter("ChiSquareCut"); + + edm::ESHandle hbuilder; + std::string sname = pset_.getParameter("TTRHBuilder"); + iRecord.getRecord().get(sname, hbuilder); + std::string hitpropagator = pset_.getParameter("HitPropagator"); + edm::ESHandle hhitpropagator; + iRecord.getRecord().getRecord().get(hitpropagator, hhitpropagator); + + bool debug = pset_.getParameter("Debug"); + //_updator = boost::shared_ptr(new SiTrackerMultiRecHitUpdator(pDD.product(), pp, sp, mp, annealingProgram)); + _updator = boost::shared_ptr(new SiTrackerMultiRecHitUpdator(hbuilder.product(),hhitpropagator.product(), Chi2Cut, annealingProgram, debug)); + // _updator = boost::shared_ptr(new SiTrackerMultiRecHitUpdator(hhitpropagator.product(),annealingProgram)); + return _updator; +} + + diff --git a/RecoTracker/SiTrackerMRHTools/plugins/SiTrackerMultiRecHitUpdatorESProducer.h b/RecoTracker/SiTrackerMRHTools/plugins/SiTrackerMultiRecHitUpdatorESProducer.h new file mode 100644 index 0000000000000..d534956b4ec37 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/plugins/SiTrackerMultiRecHitUpdatorESProducer.h @@ -0,0 +1,26 @@ +#ifndef RecoLocalTracker_ESProducers_SiTrackerMultiRecHitUpdatorESProducer_h +#define RecoLocalTracker_ESProducers_ESProducers_SiTrackerMultiRecHitUpdatorESProducer_h + +#include "FWCore/Framework/interface/ESProducer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "RecoTracker/Record/interface/MultiRecHitRecord.h" +#include "RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdator.h" + +#include + +class SiTrackerMultiRecHitUpdatorESProducer: public edm::ESProducer{ + public: + SiTrackerMultiRecHitUpdatorESProducer(const edm::ParameterSet & p); + virtual ~SiTrackerMultiRecHitUpdatorESProducer(); + boost::shared_ptr produce(const MultiRecHitRecord &); + private: + boost::shared_ptr _updator; + edm::ParameterSet pset_; +}; + + +#endif + + + + diff --git a/RecoTracker/SiTrackerMRHTools/plugins/module.cc b/RecoTracker/SiTrackerMRHTools/plugins/module.cc new file mode 100644 index 0000000000000..d6aa149330662 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/plugins/module.cc @@ -0,0 +1,13 @@ +#include "FWCore/PluginManager/interface/ModuleDef.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "RecoTracker/SiTrackerMRHTools/plugins/SiTrackerMultiRecHitUpdatorESProducer.h" +#include "RecoTracker/SiTrackerMRHTools/plugins/MultiRecHitCollectorESProducer.h" + +#include "FWCore/Framework/interface/ModuleFactory.h" +#include "FWCore/Framework/interface/ESProducer.h" +#include "FWCore/Utilities/interface/typelookup.h" + + +DEFINE_FWK_EVENTSETUP_MODULE(SiTrackerMultiRecHitUpdatorESProducer); +DEFINE_FWK_EVENTSETUP_MODULE(MultiRecHitCollectorESProducer); diff --git a/RecoTracker/SiTrackerMRHTools/python/GroupedMultiRecHitCollector_cff.py b/RecoTracker/SiTrackerMRHTools/python/GroupedMultiRecHitCollector_cff.py new file mode 100644 index 0000000000000..60ea5e474732b --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/python/GroupedMultiRecHitCollector_cff.py @@ -0,0 +1,16 @@ +import FWCore.ParameterSet.Config as cms + +import copy +from TrackingTools.KalmanUpdators.Chi2MeasurementEstimatorESProducer_cfi import * +# Chi2MeasurementEstimatorESProducer +RelaxedChi2 = copy.deepcopy(Chi2MeasurementEstimator) +RelaxedChi2.ComponentName = 'RelaxedChi2' +RelaxedChi2.MaxChi2 = 100. +#replace RelaxedChi2.nSigma = 3. +# MeasurementTracker +from RecoTracker.MeasurementDet.MeasurementTrackerESProducer_cfi import * +# MultiRecHitUpdator +from RecoTracker.SiTrackerMRHTools.SiTrackerMultiRecHitUpdator_cff import * +#MultiRecHitCollector +from RecoTracker.SiTrackerMRHTools.GroupedMultiRecHitCollector_cfi import * + diff --git a/RecoTracker/SiTrackerMRHTools/python/GroupedMultiRecHitCollector_cfi.py b/RecoTracker/SiTrackerMRHTools/python/GroupedMultiRecHitCollector_cfi.py new file mode 100644 index 0000000000000..d9c5aec39a8ac --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/python/GroupedMultiRecHitCollector_cfi.py @@ -0,0 +1,14 @@ +import FWCore.ParameterSet.Config as cms + +groupedMultiRecHitCollector = cms.ESProducer("MultiRecHitCollectorESProducer", + propagatorAlong = cms.string('RungeKuttaTrackerPropagator'), + MultiRecHitUpdator = cms.string('SiTrackerMultiRecHitUpdator'), + ComponentName = cms.string('groupedMultiRecHitCollector'), + propagatorOpposite = cms.string('OppositeRungeKuttaTrackerPropagator'), + MeasurementTrackerName = cms.string(''), + estimator = cms.string('RelaxedChi2'), + Mode = cms.string('Grouped'), + Debug = cms.bool(False) +) + + diff --git a/RecoTracker/SiTrackerMRHTools/python/SiTrackerMultiRecHitUpdator_cff.py b/RecoTracker/SiTrackerMRHTools/python/SiTrackerMultiRecHitUpdator_cff.py new file mode 100644 index 0000000000000..a11b565952167 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/python/SiTrackerMultiRecHitUpdator_cff.py @@ -0,0 +1,6 @@ +import FWCore.ParameterSet.Config as cms + +from RecoTracker.TransientTrackingRecHit.TransientTrackingRecHitBuilder_cfi import * +from TrackingTools.KalmanUpdators.TrackingRecHitPropagatorESProducer_cff import * +from RecoTracker.SiTrackerMRHTools.SiTrackerMultiRecHitUpdator_cfi import * + diff --git a/RecoTracker/SiTrackerMRHTools/python/SiTrackerMultiRecHitUpdator_cfi.py b/RecoTracker/SiTrackerMRHTools/python/SiTrackerMultiRecHitUpdator_cfi.py new file mode 100644 index 0000000000000..788d7f456120b --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/python/SiTrackerMultiRecHitUpdator_cfi.py @@ -0,0 +1,12 @@ +import FWCore.ParameterSet.Config as cms + +siTrackerMultiRecHitUpdator = cms.ESProducer("SiTrackerMultiRecHitUpdatorESProducer", + ComponentName = cms.string('SiTrackerMultiRecHitUpdator'), + TTRHBuilder = cms.string('WithAngleAndTemplate'), + HitPropagator = cms.string('trackingRecHitPropagator'), + AnnealingProgram = cms.vdouble(80.0, 9.0, 4.0, 1.0, 1.0, 1.0), + ChiSquareCut =cms.double(15.0), + Debug = cms.bool(False) +) + + diff --git a/RecoTracker/SiTrackerMRHTools/python/SimpleMultiRecHitCollector_cff.py b/RecoTracker/SiTrackerMRHTools/python/SimpleMultiRecHitCollector_cff.py new file mode 100644 index 0000000000000..eabd564f13a09 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/python/SimpleMultiRecHitCollector_cff.py @@ -0,0 +1,20 @@ +import FWCore.ParameterSet.Config as cms + +import copy +from TrackingTools.KalmanUpdators.Chi2MeasurementEstimatorESProducer_cfi import * +# Chi2MeasurementEstimatorESProducer +RelaxedChi2Simple = copy.deepcopy(Chi2MeasurementEstimator) +RelaxedChi2Simple.ComponentName = 'RelaxedChi2Simple' +RelaxedChi2Simple.MaxChi2 = 100. +#replace RelaxedChi2.nSigma = 3. +# PropagatorWithMaterialESProducer +from TrackingTools.MaterialEffects.MaterialPropagator_cfi import * +# PropagatorWithMaterialESProducer +from TrackingTools.MaterialEffects.OppositeMaterialPropagator_cfi import * +# MeasurementTracker +from RecoTracker.MeasurementDet.MeasurementTrackerESProducer_cfi import * +# MultiRecHitUpdator +from RecoTracker.SiTrackerMRHTools.SiTrackerMultiRecHitUpdator_cff import * +#MultiRecHitCollector +from RecoTracker.SiTrackerMRHTools.SimpleMultiRecHitCollector_cfi import * + diff --git a/RecoTracker/SiTrackerMRHTools/python/SimpleMultiRecHitCollector_cfi.py b/RecoTracker/SiTrackerMRHTools/python/SimpleMultiRecHitCollector_cfi.py new file mode 100644 index 0000000000000..16e7f5358cef4 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/python/SimpleMultiRecHitCollector_cfi.py @@ -0,0 +1,13 @@ +import FWCore.ParameterSet.Config as cms + +simpleMultiRecHitCollector = cms.ESProducer("MultiRecHitCollectorESProducer", + propagatorAlong = cms.string('RungeKuttaTrackerPropagator'), + MultiRecHitUpdator = cms.string('SiTrackerMultiRecHitUpdator'), + ComponentName = cms.string('simpleMultiRecHitCollector'), + MeasurementTrackerName = cms.string(''), + estimator = cms.string('RelaxedChi2Simple'), + Mode = cms.string('Simple'), + Debug = cms.bool(False) +) + + diff --git a/RecoTracker/SiTrackerMRHTools/src/GenericProjectedRecHit2D.cc b/RecoTracker/SiTrackerMRHTools/src/GenericProjectedRecHit2D.cc new file mode 100644 index 0000000000000..d34bc07949113 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/src/GenericProjectedRecHit2D.cc @@ -0,0 +1,36 @@ +#include "RecoTracker/SiTrackerMRHTools/interface/GenericProjectedRecHit2D.h" +#include "RecoTracker/TransientTrackingRecHit/interface/TSiPixelRecHit.h" +#include "RecoTracker/TransientTrackingRecHit/interface/TSiStripMatchedRecHit.h" +#include "RecoTracker/TransientTrackingRecHit/interface/TSiStripRecHit2DLocalPos.h" +#include "RecoTracker/TransientTrackingRecHit/interface/ProjectedRecHit2D.h" +#include "FWCore/Utilities/interface/Exception.h" + +GenericProjectedRecHit2D::GenericProjectedRecHit2D( const LocalPoint& pos, const LocalError& err, + const GeomDet* det, const GeomDet* originalDet, + const TransientTrackingRecHit::ConstRecHitPointer originalHit, + const TrackingRecHitPropagator* propagator) : + TrackingRecHit( *det ) //, originalHit->weight(), originalHit->getAnnealingFactor()) +{ + theOriginalDet = originalDet; + thePropagator = propagator; + theOriginalTransientHit = originalHit; + theLp = pos; + theLe = err; + theProjectionMatrix = originalHit->projectionMatrix(); + theDimension = originalHit->dimension(); + //theOriginalHit = originalTransientHit.hit()->clone(); +} + +AlgebraicVector GenericProjectedRecHit2D::parameters() const{ + AlgebraicVector result(2); + result[0] = theLp.x(); + result[1] = theLp.y(); + return result; +} + +TransientTrackingRecHit::RecHitPointer +GenericProjectedRecHit2D::clone( const TrajectoryStateOnSurface& ts, const TransientTrackingRecHitBuilder* builder) const +{ + return thePropagator->project(theOriginalTransientHit, *det(), ts, builder); +} + diff --git a/RecoTracker/SiTrackerMRHTools/src/GroupedDAFHitCollector.cc b/RecoTracker/SiTrackerMRHTools/src/GroupedDAFHitCollector.cc new file mode 100644 index 0000000000000..fba80bf4311a9 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/src/GroupedDAFHitCollector.cc @@ -0,0 +1,220 @@ +#include "RecoTracker/SiTrackerMRHTools/interface/GroupedDAFHitCollector.h" +#include "RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdator.h" +#include "RecoTracker/SiTrackerMRHTools/interface/MeasurementByLayerGrouper.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "TrackingTools/DetLayers/interface/MeasurementEstimator.h" +#include "TrackingTools/DetLayers/interface/DetLayer.h" +#include "TrackingTools/PatternTools/interface/Trajectory.h" +#include "TrackingTools/MeasurementDet/interface/TrajectoryMeasurementGroup.h" +#include "TrackingTools/TransientTrackingRecHit/interface/InvalidTransientRecHit.h" +#include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h" +#include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h" + +#include +#include + +using namespace std; + +vector GroupedDAFHitCollector::recHits(const Trajectory& traj, + const MeasurementTrackerEvent *theMTE) const +{ + + LayerMeasurements theLM (theMTE->measurementTracker(), *theMTE); + + //WARNING: At the moment the trajectories has the measurements with reversed sorting after the track smoothing + const vector& meas = traj.measurements(); + const Propagator* forwardPropagator = getPropagator(); + const Propagator* backwardPropagator = getReversePropagator(); + if (traj.direction() == alongMomentum){ + forwardPropagator = getReversePropagator(); + backwardPropagator = getPropagator(); + } + if (meas.empty()) //return TransientTrackingRecHit::ConstRecHitContainer(); + return vector(); + + //groups the TrajectoryMeasurements on a layer by layer + vector > > mol; + mol = MeasurementByLayerGrouper(getMeasurementTracker()->geometricSearchTracker())(meas); + + vector result; + + //add a protection if all the measurement are on the same layer + if(mol.size()<2)return vector(); + + //it assumes that the measurements are sorted in the smoothing direction + //TrajectoryStateOnSurface current = (*(mol.begin()+1)).second.front().updatedState(); + TrajectoryStateOnSurface current = (*(mol.rbegin()+1)).second.back().updatedState(); + //if (current.isValid()) current.rescaleError(10); + + + //protection for layers with invalid meas with no id associated + //to be fixed + //for the moment no hit are lookerd for in these layers + //remind that: + //groupedMeasurements will return at least a measurement with an invalid hit with no detid + LogDebug("MultiRecHitCollector") << "Layer " << mol.back().first << " has " << mol.back().second.size() << " measurements"; + LogTrace("MultiRecHitCollector") << "Original measurements are:"; + for( unsigned int iLay = 0; iLay < mol.size(); iLay++){ + LogTrace("MultiRecHitCollector") << " Layer " << mol.at(iLay).first << " has " << mol.at(iLay).second.size() << " measurements:"; + vector::const_iterator ibeg = (mol.at(iLay)).second.begin(); + vector::const_iterator iend = (mol.at(iLay)).second.end(); + for (vector::const_iterator imeas = ibeg; imeas != iend; ++imeas){ + if (imeas->recHit()->isValid()){ + LogTrace("MultiRecHitCollector") << " Valid Hit with DetId " << imeas->recHit()->geographicalId().rawId() + << " local position " << imeas->recHit()->hit()->localPosition() + << " global position " << imeas->recHit()->hit()->globalPosition() ; + } else { + LogTrace("MultiRecHitCollector") << " Invalid Hit with DetId " << imeas->recHit()->geographicalId().rawId(); + } + } + } + + //ERICA: I have to understand how are set the TM now. REPLACE THIS PART!! + vector groupedMeas; + if (mol.back().first) + groupedMeas = theLM.groupedMeasurements(*(mol.back().first), current, + *backwardPropagator, *(getEstimator())); + + //Since we have passed the backwardPropagator, we have to sort the detGroups in the opposite way + //(according the forward propagator, not the backward one) + vector sortedgroupedMeas; + for (vector::reverse_iterator iter = groupedMeas.rbegin(); + iter != groupedMeas.rend(); iter++){ + + sortedgroupedMeas.push_back(*iter); + + } + + //for the first layer + buildMultiRecHits(sortedgroupedMeas, result); + + + //for other layers + current = mol.back().second.front().updatedState(); + //if (current.isValid()) current.rescaleError(10); + + for( vector > >::reverse_iterator imol = + mol.rbegin() + 1; imol != mol.rend(); imol++) { + + const DetLayer* lay = (*imol).first; + LogDebug("MultiRecHitCollector") << "Layer " << lay << " has " << (*imol).second.size() << " measurements"; + //debug + vector currentLayerMeas; + if (lay) { + currentLayerMeas = theLM.groupedMeasurements(*lay, current, *forwardPropagator, *(getEstimator())); + } + + buildMultiRecHits(currentLayerMeas, result); + current = (*imol).second.front().updatedState(); + //if (current.isValid()) current.rescaleError(10); + } + + LogTrace("MultiRecHitCollector") << " Ending GroupedDAFHitCollector::recHits >> Original Measurement size " << meas.size() + << "\n >> GroupedDAFHitCollector returned " << result.size() << " measurements"; + //results are sorted in the fitting direction + + // adding a protection against too few hits and invalid hits (due to failed propagation on the same surface of the original hits) + if (result.size()>2) + { + int hitcounter=0; + //check if the vector result has more than 3 valid hits + for (vector::const_iterator iimeas = result.begin(); iimeas != result.end(); ++iimeas) + { + if(iimeas->recHit()->isValid()) hitcounter++; + } + + if(hitcounter>2) + { + return result; + } + + else return vector(); + } + + else{return vector();} + +} + +void GroupedDAFHitCollector::buildMultiRecHits(const vector& measgroup, vector& result) const { + + unsigned int initial_size = result.size(); + + //TransientTrackingRecHit::ConstRecHitContainer rhits; + if (measgroup.empty()) { + LogTrace("MultiRecHitCollector") << "No TrajectoryMeasurementGroups found for this layer\n" ; + //should we do something? + //result.push_back(InvalidTransientRecHit::build(0,TrackingRecHit::missing)); + return; + } + + //we build a MultiRecHit out of each group + //groups are sorted along momentum or opposite to momentum, + //measurements in groups are sorted with increating chi2 + LogTrace("MultiRecHitCollector") << "Found " << measgroup.size() << " groups for this layer"; + + //trajectory state to store the last valid TrajectoryState (if present) to be used + //to add an invalid Measurement in case no valid state or no valid hits are found in any group + for ( vector::const_iterator igroup = measgroup.begin(); + igroup != measgroup.end(); igroup++ ){ + + //the TrajectoryState is the first one + TrajectoryStateOnSurface state = igroup->measurements().front().predictedState(); + if (!state.isValid()){ + LogTrace("MultiRecHitCollector") << "Something wrong! no valid TSOS found in current group "; + continue; + } + + LogTrace("MultiRecHitCollector") << "This group has " << igroup->measurements().size() << " measurements"; + LogTrace("MultiRecHitCollector") << "This group has the following " << igroup->detGroup().size() + << " detector ids: " << endl; + for (DetGroup::const_iterator idet = igroup->detGroup().begin(); idet != igroup->detGroup().end(); ++idet){ + LogTrace("MultiRecHitCollector") << idet->det()->geographicalId().rawId(); + } + + vector hits; + for (vector::const_iterator imeas = igroup->measurements().begin(); + imeas != igroup->measurements().end(); imeas++){ + + //collect the non missing hits to build the MultiRecHits + //we use the recHits method; anyway only simple hits, not MultiHits should be present + if (imeas->recHit()->getType() != TrackingRecHit::missing) { + LogTrace("MultiRecHitCollector") << "This hit is valid "; + hits.push_back(imeas->recHit()->hit()); + } + else{ + LogTrace("MultiRecHitCollector") << " This hit is not valid and will not enter in the MRH. " ; + } + } + + if (hits.empty()){ + LogTrace("MultiRecHitCollector") << "No valid hits found in current group "; + continue; + } + + LogTrace("MultiRecHitCollector") << "The best TSOS in this group is " << state << " it lays on surface located at " << state.surface().position(); + + LogTrace("MultiRecHitCollector") << "For the MRH on this group the following hits will be used"; + for (vector::iterator iter = hits.begin(); iter != hits.end(); iter++){ + string validity = "valid"; + if ((*iter)->getType() == TrackingRecHit::missing ) validity = "missing !should not happen!"; + else if ((*iter)->getType() == TrackingRecHit::inactive) validity = "inactive"; + else if ((*iter)->getType() == TrackingRecHit::bad) validity = "bad"; + LogTrace("MultiRecHitCollector") << "DetId " << (*iter)->geographicalId().rawId() + << " validity: " << validity + << " surface position " << getMeasurementTracker()->geomTracker()->idToDet((*iter)->geographicalId())->position() + << " hit local position " << (*iter)->localPosition(); + } + + result.push_back(TrajectoryMeasurement(state,theUpdator->buildMultiRecHit(hits, state))); + } + //can this happen? it means that the measgroup was not empty but no valid measurement was found inside + //in this case we add an invalid measuremnt for this layer + if (result.size() == initial_size){ + LogTrace("MultiRecHitCollector") << "no valid measuremnt or no valid TSOS in none of the groups"; + //measgroup has been already checked for size != 0 + if (measgroup.back().measurements().size() != 0){ + result.push_back(measgroup.back().measurements().back()); + } + } +} + diff --git a/RecoTracker/SiTrackerMRHTools/src/MeasurementByLayerGrouper.cc b/RecoTracker/SiTrackerMRHTools/src/MeasurementByLayerGrouper.cc new file mode 100644 index 0000000000000..7e7d2a3f14884 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/src/MeasurementByLayerGrouper.cc @@ -0,0 +1,72 @@ +#include "RecoTracker/SiTrackerMRHTools/interface/MeasurementByLayerGrouper.h" +#include "TrackingTools/DetLayers/interface/DetLayer.h" +#include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h" +#include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h" +#include "FWCore/Utilities/interface/Exception.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +using namespace std; + +vector > > MeasurementByLayerGrouper::operator()(const vector& vtm) const{ + + if(vtm.empty()) + return vector > >(); + + vector > > result; + result.reserve(vtm.size()); + + vector::const_iterator start = vtm.begin(); + //here we assume that the TM on the same detLayer are consecutive (as it should) + while(start != vtm.end()) { + vector::const_iterator ipart = start; + do {ipart++;} + while(ipart != vtm.end() && + getDetLayer(*start)==getDetLayer(*ipart) && + getDetLayer(*start) != 0 //the returned pointer will be 0 in case + //the measurement contains an invalid hit with no associated detid. + //This kind of hits are at most one per layer. + //this last condition avoids that 2 consecutive measurements of this kind + //are grouped in the same layer. + //it would be useful if invalid hit out of the active area were + //given the detid reserved for the whole layer instead of 0 + ) ; + + vector group(start, ipart); + result.push_back(pair >(getDetLayer(*start), + group)); + start = ipart; + } +#ifdef debug_MeasurementByLayerGrouper_ + for (vector > >::const_iterator iter = result.begin(); iter != result.end(); iter++){ + LogTrace("MeasurementByLayerGrouper|SiTrackerMultiRecHitUpdator") << "DetLayer " << iter->first << " has " << iter->second.size() << " measurements"; + } +#endif + + + return result; +} + +const DetLayer* MeasurementByLayerGrouper::getDetLayer(const TM& tm) const { + // if the DetLayer is set in the TM... + if (tm.layer()) return tm.layer(); + + //if it corresponds to an invalid hit with no geomdet associated + //we can't retrieve the DetLayer + //because unfortunately the detlayer is not set in these cases + //returns 0 for the moment + //to be revisited + + if (tm.recHit()->det()==0){ + LogDebug("MeasurementByLayerGrouper") <<"This hit has no geomdet associated skipping... "; + return 0; + } + + //now the cases in which the detid is set + + if (!theGeomSearch) { + throw cms::Exception("MeasurementByLayerGrouper") << "Impossible to retrieve the det layer because it's not set in the TM and the pointer to the GeometricSearchTracker is 0 "; + return 0; + } + + return theGeomSearch->detLayer(tm.recHit()->det()->geographicalId()); +} diff --git a/RecoTracker/SiTrackerMRHTools/src/SiTrackerMultiRecHitUpdator.cc b/RecoTracker/SiTrackerMRHTools/src/SiTrackerMultiRecHitUpdator.cc new file mode 100644 index 0000000000000..55c41fa378383 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/src/SiTrackerMultiRecHitUpdator.cc @@ -0,0 +1,367 @@ +#include "RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdator.h" +#include "RecoTracker/SiTrackerMRHTools/interface/GenericProjectedRecHit2D.h" +#include "DataFormats/TrackerRecHit2D/interface/SiTrackerMultiRecHit.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" +#include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h" +#include "TrackingTools/TransientTrackingRecHit/interface/TrackingRecHitProjector.h" +#include "DataFormats/TrackingRecHit/interface/KfComponentsHolder.h" +#include "DataFormats/Math/interface/invertPosDefMatrix.h" +#include "DataFormats/Math/interface/ProjectMatrix.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +SiTrackerMultiRecHitUpdator::SiTrackerMultiRecHitUpdator(const TransientTrackingRecHitBuilder* builder, + const TrackingRecHitPropagator* hitpropagator, + const float Chi2Cut, + const std::vector& anAnnealingProgram, + bool debug): + theBuilder(builder), + theHitPropagator(hitpropagator), + theChi2Cut(Chi2Cut), + theAnnealingProgram(anAnnealingProgram), + debug_(debug){ + theHitCloner = static_cast(builder)->cloner(); + } + + +TransientTrackingRecHit::RecHitPointer SiTrackerMultiRecHitUpdator::buildMultiRecHit(const std::vector& rhv, + const TrajectoryStateOnSurface& tsos, + float annealing) const{ + LogTrace("SiTrackerMultiRecHitUpdator") << "Calling SiTrackerMultiRecHitUpdator::buildMultiRecHit with AnnealingFactor: " << annealing; + TransientTrackingRecHit::ConstRecHitContainer tcomponents; + for (std::vector::const_iterator iter = rhv.begin(); iter != rhv.end(); iter++){ + TransientTrackingRecHit::RecHitPointer transient = theBuilder->build(*iter); + if (transient->isValid()) tcomponents.push_back(transient); + } + return update(tcomponents, tsos, annealing); + +} + +TransientTrackingRecHit::RecHitPointer SiTrackerMultiRecHitUpdator::update( TransientTrackingRecHit::ConstRecHitPointer original, + const TrajectoryStateOnSurface& tsos, + double annealing) const{ + LogTrace("SiTrackerMultiRecHitUpdator") << "Calling SiTrackerMultiRecHitUpdator::update with AnnealingFactor: " << annealing; + + if(!tsos.isValid()) { + //return original->clone(); + throw cms::Exception("SiTrackerMultiRecHitUpdator") << "!!! MultiRecHitUpdator::update(..): tsos NOT valid!!! "; + } + + //check if to clone is the right thing + if(original->isValid()){ + if (original->transientHits().empty()){ + return theHitCloner.makeShared(original,tsos); + } + } else { + return theHitCloner.makeShared(original,tsos); + } + + TransientTrackingRecHit::ConstRecHitContainer tcomponents = original->transientHits(); + return update(tcomponents, tsos, annealing); +} + +/*------------------------------------------------------------------------------------------------------------------------*/ +TransientTrackingRecHit::RecHitPointer SiTrackerMultiRecHitUpdator::update( TransientTrackingRecHit::ConstRecHitContainer& tcomponents, + const TrajectoryStateOnSurface& tsos, + double annealing) const{ + + if (tcomponents.empty()){ + LogTrace("SiTrackerMultiRecHitUpdator") << "Empty components vector passed to SiTrackerMultiRecHitUpdator::update, returning an InvalidTransientRecHit "; + return std::make_shared(); + } + + if(!tsos.isValid()) { + LogTrace("SiTrackerMultiRecHitUpdator")<<"SiTrackerMultiRecHitUpdator::update: tsos NOT valid!!!, returning an InvalidTransientRecHit"; + return std::make_shared(); + } + + std::vector updatedcomponents; + const GeomDet* geomdet = 0; + + //running on all over the MRH components + for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator iter = tcomponents.begin(); iter != tcomponents.end(); iter++){ + + //the first rechit must belong to the same surface of TSOS + if (iter == tcomponents.begin()) { + + if (&((*iter)->det()->surface())!=&(tsos.surface())){ + throw cms::Exception("SiTrackerMultiRecHitUpdator") << "the Trajectory state and the first rechit passed to the SiTrackerMultiRecHitUpdator lay on different surfaces!: state lays on surface " << tsos.surface().position() << " hit with detid " << (*iter)->det()->geographicalId().rawId() << " lays on surface " << (*iter)->det()->surface().position(); + } + + geomdet = (*iter)->det(); + + } + + //if the rechit does not belong to the surface of the tsos + //GenericProjectedRecHit2D is used to prepagate + if (&((*iter)->det()->surface())!=&(tsos.surface())){ + + TransientTrackingRecHit::RecHitPointer cloned = theHitPropagator->project(*iter, *geomdet, tsos, theBuilder); + //if it is used a sensor by sensor grouping this should not appear + if (cloned->isValid()) updatedcomponents.push_back(cloned); + + } else { + TransientTrackingRecHit::RecHitPointer cloned = theHitCloner.makeShared(*iter,tsos); + if (cloned->isValid()){ + updatedcomponents.push_back(cloned); + } + } + } + + std::vector > mymap; + std::vector > normmap; + + double a_sum=0, c_sum=0; + + + for(std::vector::iterator ihit = updatedcomponents.begin(); + ihit != updatedcomponents.end(); ihit++) { + + double a_i = ComputeWeight(tsos, *(*ihit), false, annealing); //exp(-0.5*Chi2)/(2.*M_PI*sqrt(det)); + double c_i = ComputeWeight(tsos, *(*ihit), true, annealing); //exp(-0.5*theChi2Cut/annealing)/(2.*M_PI*sqrt(det)); + mymap.push_back(std::pair((*ihit)->hit(), a_i)); + + a_sum += a_i; + c_sum += c_i; + } + double total_sum = a_sum + c_sum; + + unsigned int counter = 0; + for(std::vector::iterator ihit = updatedcomponents.begin(); + ihit != updatedcomponents.end(); ihit++) { + + double p = ((mymap[counter].second)/total_sum > 1.e-6 ? (mymap[counter].second)/total_sum : 1.e-6); + //ORCA: float p = ((mymap[counter].second)/total_sum > 0.01 ? (mymap[counter].second)/total_sum : 1.e-6); + normmap.push_back(std::pair(mymap[counter].first, p)); + + LogTrace("SiTrackerMultiRecHitUpdator")<< " Component hit type " << typeid(*mymap[counter].first).name() + << " position " << mymap[counter].first->localPosition() + << " error " << mymap[counter].first->localPositionError() + << " with weight " << p; + counter++; + } + + SiTrackerMultiRecHitUpdator::LocalParameters param = calcParameters(tsos, normmap); + + SiTrackerMultiRecHit updated(param.first, param.second, *normmap.front().first->det(), normmap, annealing); + LogTrace("SiTrackerMultiRecHitUpdator") << " Updated Hit position " << updated.localPosition() + << " updated error " << updated.localPositionError() << std::endl; + + return std::make_shared(param.first, param.second, *normmap.front().first->det(), normmap, annealing); +} + + +//--------------------------------------------------------------------------------------------------------------- +double SiTrackerMultiRecHitUpdator::ComputeWeight(const TrajectoryStateOnSurface& tsos, + const TransientTrackingRecHit& aRecHit, bool CutWeight, double annealing) const{ + + switch (aRecHit.dimension()) { + case 1: return ComputeWeight<1>(tsos,aRecHit,CutWeight,annealing); + case 2: return ComputeWeight<2>(tsos,aRecHit,CutWeight,annealing); + case 3: return ComputeWeight<3>(tsos,aRecHit,CutWeight,annealing); + case 4: return ComputeWeight<4>(tsos,aRecHit,CutWeight,annealing); + case 5: return ComputeWeight<5>(tsos,aRecHit,CutWeight,annealing); + } + throw cms::Exception("Rec hit of invalid dimension (not 1,2,3,4,5)") << + "The value was " << aRecHit.dimension() << + ", type is " << typeid(aRecHit).name() << "\n"; +} + +//--------------------------------------------------------------------------------------------------------------- +template +double SiTrackerMultiRecHitUpdator::ComputeWeight(const TrajectoryStateOnSurface& tsos, + const TransientTrackingRecHit& aRecHit, bool CutWeight, double annealing) const { + + //define tsos position depending on the dimension of the hit + typename AlgebraicROOTObject::Vector tsospos; + if( N==1 ){ + tsospos[0]=tsos.localPosition().x(); + if(!CutWeight) + LogTrace("SiTrackerMultiRecHitUpdator")<< " TSOS position " << tsos.localPosition(); + } + else if(N==2){ + tsospos[0]=tsos.localPosition().x(); + tsospos[1]=tsos.localPosition().y(); + if(!CutWeight) + LogTrace("SiTrackerMultiRecHitUpdator")<< " TSOS position " << tsos.localPosition(); + } + else{ + if(!CutWeight) + LogTrace("SiTrackerMultiRecHitUpdator")<< " TSOS position not correct. Rec hit of invalid dimension (not 1, 2) but " + << aRecHit.dimension(); + } + + // define variables that will be used to setup the KfComponentsHolder + ProjectMatrix pf; + typename AlgebraicROOTObject::Matrix H; + typename AlgebraicROOTObject::Vector r, rMeas; + typename AlgebraicROOTObject::SymMatrix V, VMeas, W; + AlgebraicVector5 x = tsos.localParameters().vector(); + const AlgebraicSymMatrix55 &C = (tsos.localError().matrix()); + + // setup the holder with the correct dimensions and get the values + KfComponentsHolder holder; + holder.template setup(&r, &V, &H, &pf, &rMeas, &VMeas, x, C); + aRecHit.getKfComponents(holder); + + typename AlgebraicROOTObject::Vector diff; + diff = r - tsospos; + + //assume that TSOS is smoothed one + V *= annealing; + W = V; + //V += me.measuredError(*ihit);// result = b*V + H*C*H.T() + + //ierr will be set to true when inversion is successfull + bool ierr = invertPosDefMatrix(W); + + //Det2 method will preserve the content of the Matrix + //and return true when the calculation is successfull + double det; + bool ierr2 = V.Det2(det); + + if( !ierr || !ierr2) { + LogTrace("SiTrackerMultiRecHitUpdator")<<"SiTrackerMultiRecHitUpdator::ComputeWeight: W not valid!"< >& aHitMap) const{ + + AlgebraicSymMatrix22 W_sum; + AlgebraicVector2 m_sum; + + for(std::vector >::const_iterator ihit = aHitMap.begin(); ihit != aHitMap.end(); ihit++){ + + std::pair PositionAndError22; + PositionAndError22 = ComputeParameters2dim(tsos, *(ihit->first)); + + W_sum += (ihit->second*PositionAndError22.second); + m_sum += (ihit->second*(PositionAndError22.second*PositionAndError22.first)); + + } + + AlgebraicSymMatrix22 V_sum = W_sum; + bool ierr = invertPosDefMatrix(V_sum); + if( !ierr ) { + edm::LogError("SiTrackerMultiRecHitUpdator")<<"SiTrackerMultiRecHitUpdator::calcParameters: V_sum not valid!"< SiTrackerMultiRecHitUpdator::ComputeParameters2dim(const TrajectoryStateOnSurface& tsos, + const TransientTrackingRecHit& aRecHit) const{ + + switch (aRecHit.dimension()) { + case 2: return ComputeParameters2dim<2>(tsos,aRecHit); + //avoiding the not-2D hit due to the matrix final sum + case ( 1 || 3 || 4 || 5 ):{ + AlgebraicVector2 dummyVector2D (0.0,0.0); + AlgebraicSymMatrix22 dummyMatrix2D; + dummyMatrix2D(0,0) = dummyMatrix2D(1,0) = dummyMatrix2D(1,1) = 0.0; + return std::make_pair(dummyVector2D,dummyMatrix2D); + } + } + throw cms::Exception("Rec hit of invalid dimension (not 1,2,3,4,5)") << + "The value was " << aRecHit.dimension() << + ", type is " << typeid(aRecHit).name() << "\n"; +} + +//--------------------------------------------------------------------------------------------------------------- +template +std::pair SiTrackerMultiRecHitUpdator::ComputeParameters2dim(const TrajectoryStateOnSurface& tsos, + const TransientTrackingRecHit& aRecHit) const { + + // define variables that will be used to setup the KfComponentsHolder + ProjectMatrix pf; + typename AlgebraicROOTObject::Matrix H; + typename AlgebraicROOTObject::Vector r, rMeas; + typename AlgebraicROOTObject::SymMatrix V, VMeas, Wtemp; + AlgebraicVector5 x = tsos.localParameters().vector(); + const AlgebraicSymMatrix55 &C = (tsos.localError().matrix()); + + // setup the holder with the correct dimensions and get the values + KfComponentsHolder holder; + holder.template setup(&r, &V, &H, &pf, &rMeas, &VMeas, x, C); + aRecHit.getKfComponents(holder); + + Wtemp = V; + bool ierr = invertPosDefMatrix(Wtemp); + if( !ierr ) { + edm::LogError("SiTrackerMultiRecHitUpdator")<<"SiTrackerMultiRecHitUpdator::ComputeParameters2dim: W not valid!"<((*ihit)->parametersError())); + AlgebraicSymMatrix22 W(V.Inverse(ierr)); + + if(ierr != 0) { + edm::LogError("SiTrackerMultiRecHitUpdator")<<"MultiRecHit::checkParametersError: W not valid!"<weight()*W); + } + AlgebraicSymMatrix22 parametersError = W_sum.Inverse(ierr); + return LocalError(parametersError(0,0), parametersError(0,1), parametersError(1,1)); +} + +LocalPoint SiTrackerMultiRecHitUpdator::calcParameters(TransientTrackingRecHit::ConstRecHitContainer& map, const LocalError& er) const { + AlgebraicVector2 m_sum; + int ierr; + for( TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit = map.begin(); ihit != map.end(); ihit ++) { + AlgebraicVector2 m(asSVector<2>((*ihit)->parameters())); + AlgebraicSymMatrix22 V(asSMatrix<2>((*ihit)->parametersError())); + AlgebraicSymMatrix22 W(V.Inverse(ierr)); + + if(ierr != 0) { + edm::LogError("SiTrackerMultiRecHitUpdator")<<"MultiRecHit::checkParameters: W not valid!"<weight()*(W*m); + else m_sum += ((*ihit)->weight()*(W*m)); + } + AlgebraicSymMatrix22 V_sum; + + V_sum(0,0) = er.xx(); + V_sum(0,1) = er.xy(); + V_sum(1,1) = er.yy(); + //AlgebraicSymMatrix V_sum(parametersError()); + AlgebraicVector2 parameters = V_sum*m_sum; + return LocalPoint(parameters(0), parameters(1)); +} +*/ diff --git a/RecoTracker/SiTrackerMRHTools/src/SimpleDAFHitCollector.cc b/RecoTracker/SiTrackerMRHTools/src/SimpleDAFHitCollector.cc new file mode 100644 index 0000000000000..ef20c35b309e8 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/src/SimpleDAFHitCollector.cc @@ -0,0 +1,182 @@ +#include "RecoTracker/SiTrackerMRHTools/interface/SimpleDAFHitCollector.h" +#include "RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdator.h" +#include "TrackingTools/DetLayers/interface/MeasurementEstimator.h" +#include "TrackingTools/PatternTools/interface/Trajectory.h" +#include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h" +#include "TrackingTools/MeasurementDet/interface/MeasurementDet.h" +#include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h" + +#include +#include + +#define _debug_SimpleDAFHitCollector_ + +using namespace std; + +vector SimpleDAFHitCollector::recHits(const Trajectory& traj, const MeasurementTrackerEvent *theMTE) const{ + + LogTrace("MultiRecHitCollector") << " Calling SimpleDAFHitCollector::recHits" << std::endl; + + //WARNING: At the moment the trajectories has the measurements + //with reversed sorting after the track smoothing + const vector& meas = traj.measurements(); + + if (meas.empty()) return vector(); + + //debug + LogTrace("MultiRecHitCollector") << "Original measurements are:"; + for(vector::const_iterator itrajmeas = meas.begin(); itrajmeas < meas.end(); + itrajmeas++) { + if (itrajmeas->recHit()->isValid()){ + LogTrace("MultiRecHitCollector") << " Valid Hit with DetId " << itrajmeas->recHit()->geographicalId().rawId() + << " local position " << itrajmeas->recHit()->hit()->localPosition() + << " global position " << itrajmeas->recHit()->hit()->globalPosition() ; + } else { + LogTrace("MultiRecHitCollector") << " Invalid Hit with DetId " << itrajmeas->recHit()->geographicalId().rawId(); + } + } + + //groups the TrajectoryMeasurements on a sensor by sensor + //we have to sort the TrajectoryMeasurements in the opposite way in the fitting direction + vector result; + for(vector::const_reverse_iterator itrajmeas = meas.rbegin(); itrajmeas < meas.rend(); + itrajmeas++) { + + DetId id = itrajmeas->recHit()->geographicalId(); + MeasurementDetWithData measDet = theMTE->idToDet(id); + tracking::TempMeasurements tmps; + std::vector hits; + + TrajectoryStateOnSurface current = itrajmeas->updatedState(); + //the error is scaled in order to take more "compatible" hits + if (current.isValid()) current.rescaleError(10); + + TrajectoryStateOnSurface state = itrajmeas->predictedState(); + if (!state.isValid()){ + LogTrace("MultiRecHitCollector") << "Something wrong! no valid TSOS found in current group "; + continue; + } + //collected hits compatible with the itrajmeas + if( measDet.measurements(current, *(getEstimator()), tmps)){ + LogTrace("MultiRecHitCollector") << " Found " << tmps.size() << " compatible measurements"; + for (std::size_t i=0; i!=tmps.size(); ++i){ + DetId idtemps = tmps.hits[i]->geographicalId(); + + if( idtemps == id && tmps.hits[i]->hit()->isValid() ) { + hits.push_back(tmps.hits[i]->hit()); + } + } + + //I will keep the Invalid hit, IF this is not the first one + if (hits.empty()){ + LogTrace("MultiRecHitCollector") << " -> but no valid hits found in current group."; + + if( result.empty() ) continue; + + result.push_back(TrajectoryMeasurement(state, + std::make_shared(measDet.mdet().geomDet(), TrackingRecHit::missing))); + } else { + //measurements in groups are sorted with increating chi2 + //sort( *hits.begin(), *hits.end(), TrajMeasLessEstim()); + + LogTrace("MultiRecHitCollector") << " -> " << hits.size() << " valid hits for this sensor."; + + //building a MultiRecHit out of each sensor group + result.push_back(TrajectoryMeasurement(state,theUpdator->buildMultiRecHit(hits, state))); + } + } else { + LogTrace("MultiRecHitCollector") << " No measurements found in current group."; + + if( result.empty() ) continue; + + result.push_back(TrajectoryMeasurement(state, + std::make_shared(measDet.mdet().geomDet(), TrackingRecHit::missing))); + + } + + } + LogTrace("MultiRecHitCollector") << " Ending SimpleDAFHitCollector::recHits >> " << result.size(); + + //adding a protection against too few hits and invalid hits + //(due to failed propagation on the same surface of the original hits) + if (result.size()>2) + { + int hitcounter=0; + //check if the vector result has more than 3 valid hits + for (vector::const_iterator iimeas = result.begin(); iimeas != result.end(); ++iimeas) { + if(iimeas->recHit()->isValid()) hitcounter++; + } + + if(hitcounter>2) + return result; + else return vector(); + } + + else{return vector();} + +} + +/*vector SimpleDAFHitCollector::recHits(const Trajectory& traj, const MeasurementTrackerEvent *theMTE) const{ + + + //it assumes that the measurements are sorted in the smoothing direction + vector meas = traj.measurements(); + if (meas.empty()) //return TransientTrackingRecHit::ConstRecHitContainer(); + return vector(); + + + //TransientTrackingRecHit::ConstRecHitContainer result; + vector result; + + for(vector::reverse_iterator imeas = meas.rbegin(); imeas != meas.rend(); imeas++) { + //if the rechit is associated to a valid detId + if(imeas->recHit()->geographicalId().rawId()){ + MeasurementDetWithData md = theMTE->idToDet(imeas->recHit()->geographicalId()); + vector currentLayerMeas = md.fastMeasurements(imeas->updatedState(), imeas->updatedState(), *thePropagator, *(getEstimator())); + buildMultiRecHits(currentLayerMeas, result); + } else { + result.push_back(*imeas); + } + } + LogTrace("MultiRecHitCollector") << "Original Measurement size " << meas.size() << " SimpleDAFHitCollector returned " << result.size() << " rechits"; + //results are sorted in the fitting direction + return result; +} + +void SimpleDAFHitCollector::buildMultiRecHits(const vector& vmeas, vector& result) const { + + if (vmeas.empty()) { + LogTrace("MultiRecHitCollector") << "fastMeasurements returned an empty vector, should not happen " ; + //should we do something? + //result.push_back(InvalidTransientRecHit::build(0,TrackingRecHit::missing)); + return; + } + + TrajectoryStateOnSurface state = vmeas.front().predictedState(); + + if (state.isValid()==false){ + LogTrace("MultiRecHitCollector") << "first state is invalid; skipping "; + return; + } + + vector hits; + // TransientTrackingRecHit::ConstRecHitContainer hits; + + for (vector::const_iterator imeas = vmeas.begin(); imeas != vmeas.end(); imeas++){ + if (imeas->recHit()->getType() != TrackingRecHit::missing) { + LogTrace("MultiRecHitCollector") << "This hit is valid "; + hits.push_back(imeas->recHit()->hit()); + } + } + + if (hits.empty()){ + LogTrace("MultiRecHitCollector") << "No valid hits found "; + return; + } + + result.push_back(TrajectoryMeasurement(state,theUpdator->buildMultiRecHit(hits, state))); + // result.push_back(TrajectoryMeasurement(state,theUpdator->update(hits, state))); + +} + +*/ diff --git a/RecoTracker/SiTrackerMRHTools/src/module.cc b/RecoTracker/SiTrackerMRHTools/src/module.cc new file mode 100644 index 0000000000000..96f633a0da9bd --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/src/module.cc @@ -0,0 +1,10 @@ + +#include "RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdator.h" +//#include "RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdatorMTF.h" +#include "RecoTracker/SiTrackerMRHTools/interface/MultiRecHitCollector.h" +//#include "RecoTracker/SiTrackerMRHTools/interface/MultiTrackFilterHitCollector.h" + +#include "FWCore/Utilities/interface/typelookup.h" + +TYPELOOKUP_DATA_REG(SiTrackerMultiRecHitUpdator); +TYPELOOKUP_DATA_REG(MultiRecHitCollector); diff --git a/RecoTracker/SiTrackerMRHTools/test/DAFRecoTrack_cfg.py b/RecoTracker/SiTrackerMRHTools/test/DAFRecoTrack_cfg.py new file mode 100644 index 0000000000000..01dac466a1bc2 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/test/DAFRecoTrack_cfg.py @@ -0,0 +1,69 @@ +import FWCore.ParameterSet.Config as cms + +# inspired from RecoTracker/TrackProducer/test/TrackRefit.py + +process = cms.Process("Refitting") + +### standard MessageLoggerConfiguration +process.load("FWCore.MessageService.MessageLogger_cfi") + +### Standard Configurations +process.load("Configuration.StandardSequences.Services_cff") +process.load('Configuration/StandardSequences/Geometry_cff') +process.load('Configuration/StandardSequences/Reconstruction_cff') +process.load('Configuration/StandardSequences/MagneticField_cff') + +### Conditions +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +#from Configuration.AlCa.GlobalTag import GlobalTag +#process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:startup', '') +process.GlobalTag.globaltag = 'START71_V1::All'#POSTLS171_V1::All' + +### Track Refitter +process.load("RecoTracker.TrackProducer.TrackRefitters_cff") +process.load("RecoTracker.TrackProducer.CTFFinalFitWithMaterialDAF_cff") +process.ctfWithMaterialTracksDAF.TrajectoryInEvent = True +process.ctfWithMaterialTracksDAF.src = 'TrackRefitter' +process.ctfWithMaterialTracksDAF.TrajAnnealingSaving = True +process.MRHFittingSmoother.EstimateCut = -1 +process.MRHFittingSmoother.MinNumberOfHits = 3 + +process.MessageLogger = cms.Service("MessageLogger", + destinations = cms.untracked.vstring("debugTracking"), #1 + debugModules = cms.untracked.vstring("*"), #2 + categories = cms.untracked.vstring("SiTrackerMultiRecHitUpdator"), #3 + debugTracking = 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 + SiTrackerMultiRecHitUpdator = cms.untracked.PSet(limit = cms.untracked.int32(-1)) #7 + ) + ) + +process.source = cms.Source("PoolSource", +# fileNames = cms.untracked.vstring('file:reco_trk_TTbar_13_5evts.root') + fileNames = cms.untracked.vstring('file:reco_trk_SingleMuPt10_UP15_10evts.root') +) +process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(1)) + +process.out = cms.OutputModule("PoolOutputModule", + outputCommands = cms.untracked.vstring('drop *_*_*_*', + 'keep *_siPixelClusters_*_*', + 'keep *_siStripClusters_*_*', + 'keep *_siPixelDigis_*_*', + 'keep *_siStripDigis_*_*', + 'keep *_offlineBeamSpot_*_*', + 'keep *_pixelVertices_*_*', + 'keep *_siStripMatchedRecHits_*_*', + 'keep *_initialStepSeeds_*_*', + 'keep recoTracks_*_*_*', + 'keep recoTrackExtras_*_*_*', + 'keep TrackingRecHitsOwned_*_*_*'), + fileName = cms.untracked.string('refit_DAF_SingleMuPt10_UP15_100evts.root') + ) + +process.p = cms.Path(process.MeasurementTrackerEvent*process.TrackRefitter*process.ctfWithMaterialTracksDAF) +process.o = cms.EndPath(process.out) + +process.schedule = cms.Schedule(process.p,process.o) + + diff --git a/RecoTracker/SiTrackerMRHTools/test/DAFRecoTrack_debug_cfg.py b/RecoTracker/SiTrackerMRHTools/test/DAFRecoTrack_debug_cfg.py new file mode 100644 index 0000000000000..4981331088651 --- /dev/null +++ b/RecoTracker/SiTrackerMRHTools/test/DAFRecoTrack_debug_cfg.py @@ -0,0 +1,72 @@ +import FWCore.ParameterSet.Config as cms + +# inspired from RecoTracker/TrackProducer/test/TrackRefit.py + +process = cms.Process("Refitting") + +### standard MessageLoggerConfiguration +process.load("FWCore.MessageService.MessageLogger_cfi") + +### Standard Configurations +process.load("Configuration.StandardSequences.Services_cff") +process.load('Configuration/StandardSequences/Geometry_cff') +process.load('Configuration/StandardSequences/Reconstruction_cff') +process.load('Configuration/StandardSequences/MagneticField_cff') + +### Conditions +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +#from Configuration.AlCa.GlobalTag import GlobalTag +#process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:startup', '') +process.GlobalTag.globaltag = 'START71_V1::All'#POSTLS171_V1::All' + +### Track Refitter +process.load("RecoTracker.TrackProducer.TrackRefitters_cff") +process.load("RecoTracker.TrackProducer.CTFFinalFitWithMaterialDAF_cff") +process.ctfWithMaterialTracksDAF.TrajectoryInEvent = True +process.ctfWithMaterialTracksDAF.src = 'TrackRefitter' +process.ctfWithMaterialTracksDAF.TrajAnnealingSaving = True +process.MRHFittingSmoother.EstimateCut = -1 +process.MRHFittingSmoother.MinNumberOfHits = 3 + +# debug +process.MessageLogger = cms.Service("MessageLogger", + destinations = cms.untracked.vstring("debugTracking"), #1 + debugModules = cms.untracked.vstring("*"), #2 + #categories = cms.untracked.vstring("MultiRecHitCollector"),#SiTrackerMultiRecHitUpdator"), #3 + categories = cms.untracked.vstring("DAFTrackProducerAlgorithm"),#SiTrackerMultiRecHitUpdator"), #3 + debugTracking = 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 + #MultiRecHitCollector = cms.untracked.PSet(limit = cms.untracked.int32(-1)) #7 + DAFTrackProducerAlgorithm = cms.untracked.PSet(limit = cms.untracked.int32(-1)) #7 + ) + ) + +process.source = cms.Source("PoolSource", +# fileNames = cms.untracked.vstring('file:reco_trk_TTbar_13_5evts.root') + fileNames = cms.untracked.vstring('file:reco_trk_SingleMuPt10_UP15_10evts.root') +) +process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(1)) + +process.out = cms.OutputModule("PoolOutputModule", + outputCommands = cms.untracked.vstring('drop *_*_*_*', + 'keep *_siPixelClusters_*_*', + 'keep *_siStripClusters_*_*', + 'keep *_siPixelDigis_*_*', + 'keep *_siStripDigis_*_*', + 'keep *_offlineBeamSpot_*_*', + 'keep *_pixelVertices_*_*', + 'keep *_siStripMatchedRecHits_*_*', + 'keep *_initialStepSeeds_*_*', + 'keep recoTracks_*_*_*', + 'keep recoTrackExtras_*_*_*', + 'keep TrackingRecHitsOwned_*_*_*'), + fileName = cms.untracked.string('refit_DAF_SingleMuPt10_UP15_100evts.root') + ) + +process.p = cms.Path(process.MeasurementTrackerEvent*process.TrackRefitter*process.ctfWithMaterialTracksDAF) +process.o = cms.EndPath(process.out) + +process.schedule = cms.Schedule(process.p,process.o) + + diff --git a/RecoTracker/TrackProducer/BuildFile.xml b/RecoTracker/TrackProducer/BuildFile.xml index d598e9aa3ba8a..30732fd9a3385 100644 --- a/RecoTracker/TrackProducer/BuildFile.xml +++ b/RecoTracker/TrackProducer/BuildFile.xml @@ -20,6 +20,7 @@ + diff --git a/RecoTracker/TrackProducer/interface/DAFTrackProducerAlgorithm.h b/RecoTracker/TrackProducer/interface/DAFTrackProducerAlgorithm.h new file mode 100644 index 0000000000000..767590f6ad864 --- /dev/null +++ b/RecoTracker/TrackProducer/interface/DAFTrackProducerAlgorithm.h @@ -0,0 +1,97 @@ +/** \class DAFTrackProducerAlgorithm + * All is needed to run the deterministic annealing algorithm. Ported from ORCA. + * + * \author tropiano, genta + * \review in May 2014 by brondolin + */ + +#ifndef DAFTrackProducerAlgorithm_h +#define DAFTrackProducerAlgorithm_h + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h" +#include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h" +#include "DataFormats/TrajectorySeed/interface/PropagationDirection.h" +#include "DataFormats/BeamSpot/interface/BeamSpot.h" + +class MagneticField; +class TrackingGeometry; +class TrajAnnealing; +class TrajectoryFitter; +class Trajectory; +class TrajectoryStateOnSurface; +class TransientTrackingRecHitBuilder; +class MultiRecHitCollector; +class SiTrackerMultiRecHitUpdator; +class MeasurementTrackerEvent; +namespace reco{ + class Track; +} + +class DAFTrackProducerAlgorithm { + + typedef std::pair > AlgoProduct; + typedef std::vector< AlgoProduct > AlgoProductCollection; + typedef std::vector TrajAnnealingCollection; + + public: + + DAFTrackProducerAlgorithm(const edm::ParameterSet& conf); + ~DAFTrackProducerAlgorithm() {} + + /// Run the Final Fit taking TrackCandidates as input + void runWithCandidate(const TrackingGeometry *, + const MagneticField *, + //const TrackCandidateCollection&, + const std::vector&, + const MeasurementTrackerEvent *measTk, + const TrajectoryFitter *, + const TransientTrackingRecHitBuilder*, + const MultiRecHitCollector* measurementTracker, + const SiTrackerMultiRecHitUpdator*, + const reco::BeamSpot&, + AlgoProductCollection &, + TrajAnnealingCollection &, + bool ) const; + + private: + /// Construct Tracks to be put in the event + bool buildTrack(const Trajectory, + AlgoProductCollection& algoResults, + float, + const reco::BeamSpot&) const; + + /// accomplishes the fitting-smoothing step for each annealing value + Trajectory fit(const std::pair& hits, + const TrajectoryFitter * theFitter, + Trajectory vtraj) const; + + //calculates the ndof according to the DAF prescription + float calculateNdof(const Trajectory vtraj) const; + + //creates MultiRecHits out of a KF trajectory + std::pair collectHits( + const Trajectory vtraj, + const MultiRecHitCollector* measurementCollector, + const MeasurementTrackerEvent *measTk ) const; + + //updates the hits with the specified annealing factor + std::pair updateHits( + const Trajectory vtraj, + const SiTrackerMultiRecHitUpdator* updator, + double annealing) const; + + //removes from the trajectory isolated hits with very low weight + void filter(const TrajectoryFitter* fitter, + std::vector& input, + int minhits, std::vector& output, + const TransientTrackingRecHitBuilder* builder) const; + + int checkHits( Trajectory iInitTraj, const Trajectory iFinalTraj) const; + + edm::ParameterSet conf_; + int minHits_; +}; + +#endif diff --git a/RecoTracker/TrackProducer/plugins/BuildFile.xml b/RecoTracker/TrackProducer/plugins/BuildFile.xml index c8b247d2e4b1e..67f546b8021f8 100644 --- a/RecoTracker/TrackProducer/plugins/BuildFile.xml +++ b/RecoTracker/TrackProducer/plugins/BuildFile.xml @@ -1,3 +1,5 @@ + + diff --git a/RecoTracker/TrackProducer/plugins/DAFTrackProducer.cc b/RecoTracker/TrackProducer/plugins/DAFTrackProducer.cc new file mode 100644 index 0000000000000..2752128ecf6ca --- /dev/null +++ b/RecoTracker/TrackProducer/plugins/DAFTrackProducer.cc @@ -0,0 +1,133 @@ +#include + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "RecoTracker/Record/interface/MultiRecHitRecord.h" +#include "RecoTracker/SiTrackerMRHTools/interface/MultiRecHitCollector.h" +#include "RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdator.h" +#include "RecoTracker/TrackProducer/plugins/DAFTrackProducer.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "TrackingTools/PatternTools/interface/Trajectory.h" +#include "TrackingTools/PatternTools/interface/TrajAnnealing.h" +#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h" +#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" + + + +DAFTrackProducer::DAFTrackProducer(const edm::ParameterSet& iConfig): + KfTrackProducerBase(iConfig.getParameter("TrajectoryInEvent"),false), + theAlgo(iConfig) +{ + setConf(iConfig); + setSrc( consumes(iConfig.getParameter( "src" )), + consumes(iConfig.getParameter( "beamSpot" )), + consumes(iConfig.getParameter( "MeasurementTrackerEvent") )); + src_ = consumes(iConfig.getParameter( "src" )); + setAlias( iConfig.getParameter( "@module_label" ) ); + + //register your products + produces().setBranchAlias( alias_ + "Tracks" ); + produces().setBranchAlias( alias_ + "TrackExtras" ); + produces().setBranchAlias( alias_ + "RecHits" ); + produces >(); + produces(); + produces().setBranchAlias( alias_ + "TrajectoryAnnealing" ); + + TrajAnnSaving_ = iConfig.getParameter("TrajAnnealingSaving"); +} + + +void DAFTrackProducer::produce(edm::Event& theEvent, const edm::EventSetup& setup) +{ + edm::LogInfo("DAFTrackProducer") << "Analyzing event number: " << theEvent.id() << "\n"; + + //empty output collections + std::auto_ptr outputRHColl (new TrackingRecHitCollection); + std::auto_ptr outputTColl(new reco::TrackCollection); + std::auto_ptr outputTEColl(new reco::TrackExtraCollection); + std::auto_ptr > outputTrajectoryColl(new std::vector); + std::auto_ptr outputTrajAnnColl(new TrajAnnealingCollection); + + //declare and get stuff to be retrieved from ES + edm::ESHandle theG; + edm::ESHandle theMF; + edm::ESHandle theFitter; + edm::ESHandle thePropagator; + edm::ESHandle theMeasTk; + edm::ESHandle theBuilder; + getFromES(setup,theG,theMF,theFitter,thePropagator,theMeasTk,theBuilder); + + //get additional es_modules needed by the DAF + edm::ESHandle measurementCollectorHandle; + std::string measurementCollectorName = getConf().getParameter("MeasurementCollector"); + setup.get().get(measurementCollectorName, measurementCollectorHandle); + edm::ESHandle updatorHandle; + std::string updatorName = getConf().getParameter("UpdatorName"); + setup.get().get(updatorName, updatorHandle); + + //get MeasurementTrackerEvent + edm::Handle mte; + theEvent.getByToken(mteSrc_, mte); + + + //declare and get TrackColection + AlgoProductCollection algoResults; + reco::BeamSpot bs; + TrajAnnealingCollection trajannResults; + try{ + + edm::Handle > theTrajectoryCollection; + getFromEvt(theEvent,theTrajectoryCollection,bs); + + //obsolete? + //measurementCollectorHandle->updateEvent(theEvent); + + //run the algorithm + LogDebug("DAFTrackProducer") << "run the DAF algorithm" << "\n"; + theAlgo.runWithCandidate(theG.product(), theMF.product(), *theTrajectoryCollection, &*mte, + theFitter.product(), theBuilder.product(), + measurementCollectorHandle.product(), updatorHandle.product(), bs, + algoResults, trajannResults, TrajAnnSaving_); + + } catch (cms::Exception &e){ + edm::LogInfo("DAFTrackProducer") << "cms::Exception caught!!!" << "\n" << e << "\n"; + throw; + } + + //put everything in the event + putInEvt(theEvent, thePropagator.product(),theMeasTk.product(), + outputRHColl, outputTColl, outputTEColl, + outputTrajectoryColl, algoResults, theBuilder.product()); + putInEvtTrajAnn(theEvent, trajannResults, outputTrajAnnColl); + + LogDebug("DAFTrackProducer") << "end the DAF algorithm." << "\n"; +} + +void DAFTrackProducer::getFromEvt(edm::Event& theEvent,edm::Handle& theTrajectoryCollection, reco::BeamSpot& bs) +{ + + //get the TrajectoryCollection from the event + //WARNING: src has always to be redefined in cfg file + theEvent.getByToken(src_,theTrajectoryCollection ); + + //get the BeamSpot + edm::Handle recoBeamSpotHandle; + theEvent.getByToken(bsSrc_,recoBeamSpotHandle); + bs = *recoBeamSpotHandle; + +} + +void DAFTrackProducer::putInEvtTrajAnn(edm::Event& theEvent, TrajAnnealingCollection & trajannResults, + std::auto_ptr& outputTrajAnnColl){ + const int size = trajannResults.size(); + outputTrajAnnColl->reserve(size); + + for(unsigned int i = 0; i < trajannResults.size() ; i++){ + outputTrajAnnColl->push_back(trajannResults[i]); + } + + theEvent.put( outputTrajAnnColl ); +} + diff --git a/RecoTracker/TrackProducer/plugins/DAFTrackProducer.h b/RecoTracker/TrackProducer/plugins/DAFTrackProducer.h new file mode 100644 index 0000000000000..c96068c706265 --- /dev/null +++ b/RecoTracker/TrackProducer/plugins/DAFTrackProducer.h @@ -0,0 +1,37 @@ +/** \class DAFTrackProducer + * EDProducer for DAFTrackProducerAlgorithm. + * + * \author tropiano, genta + * \review in May 2014 by brondolin + */ + +#ifndef DAFTrackProducer_h +#define DAFTrackProducer_h + +#include "RecoTracker/TrackProducer/interface/KfTrackProducerBase.h" +#include "RecoTracker/TrackProducer/interface/DAFTrackProducerAlgorithm.h" +#include "TrackingTools/PatternTools/interface/Trajectory.h" +#include "TrackingTools/PatternTools/interface/TrajAnnealing.h" + +class DAFTrackProducer : public KfTrackProducerBase, public edm::EDProducer { +public: + + typedef std::vector TrajectoryCollection; +// typedef std::vector TrajAnnealingCollection; + explicit DAFTrackProducer(const edm::ParameterSet& iConfig); + + // Implementation of produce method + virtual void produce(edm::Event&, const edm::EventSetup&); + +private: + DAFTrackProducerAlgorithm theAlgo; + void getFromEvt(edm::Event&, edm::Handle&, reco::BeamSpot&); + void putInEvtTrajAnn(edm::Event& theEvent, TrajAnnealingCollection & trajannResults, + std::auto_ptr& selTrajAnn); + + bool TrajAnnSaving_; + edm::EDGetToken src_; + +}; + +#endif diff --git a/RecoTracker/TrackProducer/plugins/SealModules.cc b/RecoTracker/TrackProducer/plugins/SealModules.cc index 6950872b7c83c..9c87e93069668 100644 --- a/RecoTracker/TrackProducer/plugins/SealModules.cc +++ b/RecoTracker/TrackProducer/plugins/SealModules.cc @@ -1,5 +1,6 @@ #include "FWCore/Framework/interface/MakerMacros.h" +#include "RecoTracker/TrackProducer/plugins/DAFTrackProducer.h" #include "RecoTracker/TrackProducer/plugins/TrackProducer.h" #include "RecoTracker/TrackProducer/plugins/TrackRefitter.h" #include "RecoTracker/TrackProducer/plugins/GsfTrackProducer.h" @@ -7,6 +8,7 @@ #include "RecoTracker/TrackProducer/plugins/ExtraFromSeeds.h" // +DEFINE_FWK_MODULE(DAFTrackProducer); DEFINE_FWK_MODULE(TrackProducer); DEFINE_FWK_MODULE(TrackRefitter); DEFINE_FWK_MODULE(GsfTrackProducer); diff --git a/RecoTracker/TrackProducer/plugins/TrackRefitter.cc b/RecoTracker/TrackProducer/plugins/TrackRefitter.cc index 3818b83a0a3b1..d670911bf3279 100644 --- a/RecoTracker/TrackProducer/plugins/TrackRefitter.cc +++ b/RecoTracker/TrackProducer/plugins/TrackRefitter.cc @@ -74,6 +74,9 @@ void TrackRefitter::produce(edm::Event& theEvent, const edm::EventSetup& setup) { edm::Handle theTCollection; getFromEvt(theEvent,theTCollection,bs); + + LogDebug("TrackRefitter") << "TrackRefitter::produce(none):Number of Trajectories:" << (*theTCollection).size(); + if (theTCollection.failedToGet()){ edm::LogError("TrackRefitter")<<"could not get the reco::TrackCollection."; break;} LogDebug("TrackRefitter") << "run the algorithm" << "\n"; diff --git a/RecoTracker/TrackProducer/python/CTFFinalFitWithMaterialDAF_cff.py b/RecoTracker/TrackProducer/python/CTFFinalFitWithMaterialDAF_cff.py new file mode 100644 index 0000000000000..85d27efa99847 --- /dev/null +++ b/RecoTracker/TrackProducer/python/CTFFinalFitWithMaterialDAF_cff.py @@ -0,0 +1,41 @@ +import FWCore.ParameterSet.Config as cms + + +# KFUpdatoerESProducer +from TrackingTools.KalmanUpdators.KFUpdatorESProducer_cfi import * + +# Chi2MeasurementEstimatorESProducer this is used by the fitting-smoother +from TrackingTools.KalmanUpdators.MRHChi2MeasurementEstimatorESProducer_cfi import * + +import copy +from TrackingTools.MaterialEffects.OppositeMaterialPropagator_cfi import * +# PropagatorWithMaterialESProducer +OppositeRungeKuttaTrackerPropagator = copy.deepcopy(OppositeMaterialPropagator) +OppositeRungeKuttaTrackerPropagator.ComponentName = 'OppositeRungeKuttaTrackerPropagator' +OppositeRungeKuttaTrackerPropagator.useRungeKutta = True + +# KFTrajectoryFitterESProducer +from TrackingTools.TrackFitters.TrackFitters_cff import * + +#MultiMeasurementTracker +from RecoTracker.SiTrackerMRHTools.GroupedMultiRecHitCollector_cff import * +from RecoTracker.SiTrackerMRHTools.SimpleMultiRecHitCollector_cff import * + +#multiRecHitUpdator +from RecoTracker.SiTrackerMRHTools.SiTrackerMultiRecHitUpdator_cff import * + + +# stripCPE +from RecoLocalTracker.SiStripRecHitConverter.StripCPEfromTrackAngle_cfi import * +from RecoLocalTracker.SiStripRecHitConverter.SiStripRecHitMatcher_cfi import * + +# pixelCPE +from RecoLocalTracker.SiPixelRecHits.PixelCPEParmError_cfi import * + +#TransientTrackingBuilder +from RecoTracker.TransientTrackingRecHit.TransientTrackingRecHitBuilder_cfi import * + +# TrackProducer +from RecoTracker.TrackProducer.CTFFinalFitWithMaterialDAF_cfi import * + + diff --git a/RecoTracker/TrackProducer/python/CTFFinalFitWithMaterialDAF_cfi.py b/RecoTracker/TrackProducer/python/CTFFinalFitWithMaterialDAF_cfi.py new file mode 100644 index 0000000000000..f8f6046e049ef --- /dev/null +++ b/RecoTracker/TrackProducer/python/CTFFinalFitWithMaterialDAF_cfi.py @@ -0,0 +1,18 @@ +import FWCore.ParameterSet.Config as cms + +ctfWithMaterialTracksDAF = cms.EDProducer("DAFTrackProducer", + src = cms.InputTag("DAFTrackCandidateMaker"), + UpdatorName = cms.string('SiTrackerMultiRecHitUpdator'), + beamSpot = cms.InputTag("offlineBeamSpot"), + Fitter = cms.string('MRHFittingSmoother'), + MeasurementCollector = cms.string('simpleMultiRecHitCollector'), + NavigationSchool = cms.string(''), + MeasurementTrackerEvent = cms.InputTag('MeasurementTrackerEvent'), + TrajectoryInEvent = cms.bool(False), + TTRHBuilder = cms.string('WithAngleAndTemplate'), + Propagator = cms.string('RungeKuttaTrackerPropagator'), + MinHits = cms.int32(3), + TrajAnnealingSaving = cms.bool(False) +) + + diff --git a/RecoTracker/TrackProducer/src/DAFTrackProducerAlgorithm.cc b/RecoTracker/TrackProducer/src/DAFTrackProducerAlgorithm.cc new file mode 100644 index 0000000000000..70ce70c720dc8 --- /dev/null +++ b/RecoTracker/TrackProducer/src/DAFTrackProducerAlgorithm.cc @@ -0,0 +1,404 @@ +#include "DataFormats/TrackCandidate/interface/TrackCandidate.h" +#include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h" +#include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h" +#include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "Geometry/CommonDetUnit/interface/TrackingGeometry.h" +#include "MagneticField/Engine/interface/MagneticField.h" +#include "RecoTracker/TrackProducer/interface/DAFTrackProducerAlgorithm.h" +#include "RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdator.h" +#include "RecoTracker/SiTrackerMRHTools/interface/MultiRecHitCollector.h" +#include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" +#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h" +#include "RecoTracker/TransientTrackingRecHit/interface/TkClonerImpl.h" +#include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h" +#include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h" +#include "TrackingTools/TrackFitters/interface/TrajectoryStateWithArbitraryError.h" +#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h" +#include "DataFormats/TrackerRecHit2D/interface/SiTrackerMultiRecHit.h" +#include "DataFormats/TrackerRecHit2D/interface/OmniClusterRef.h" +#include "DataFormats/TrackerRecHit2D/interface/TkCloner.h" +#include "TrackingTools/PatternTools/interface/TrajAnnealing.h" + + +DAFTrackProducerAlgorithm::DAFTrackProducerAlgorithm(const edm::ParameterSet& conf): + conf_(conf), + minHits_(conf.getParameter("MinHits")){} + + +void DAFTrackProducerAlgorithm::runWithCandidate(const TrackingGeometry * theG, + const MagneticField * theMF, + const std::vector& theTrajectoryCollection, + const MeasurementTrackerEvent *measTk, + const TrajectoryFitter * theFitter, + const TransientTrackingRecHitBuilder* builder, + const MultiRecHitCollector* measurementCollector, + const SiTrackerMultiRecHitUpdator* updator, + const reco::BeamSpot& bs, + AlgoProductCollection& algoResults, + TrajAnnealingCollection& trajann, + bool TrajAnnSaving_) const +{ + LogDebug("DAFTrackProducerAlgorithm") << "Number of Trajectories: " << theTrajectoryCollection.size() << "\n"; + int cont = 0; + + //running on src trajectory collection + for (std::vector::const_iterator ivtraj = theTrajectoryCollection.begin(); + ivtraj != theTrajectoryCollection.end(); ivtraj++) + { + + float ndof = 0; + Trajectory currentTraj; + + //getting trajectory + //no need to have std::vector vtraj ! + if ( (*ivtraj).isValid() ){ + + LogDebug("DAFTrackProducerAlgorithm") << "The trajectory is valid. \n"; + + //getting the MultiRecHit collection and the trajectory with a first fit-smooth round + std::pair hits = + collectHits(*ivtraj, measurementCollector, &*measTk); + + currentTraj = fit(hits, theFitter, *ivtraj); + + //starting the annealing program + for (std::vector::const_iterator ian = updator->getAnnealingProgram().begin(); + ian != updator->getAnnealingProgram().end(); ian++){ + + if (currentTraj.isValid()){ + + LogDebug("DAFTrackProducerAlgorithm") << "Seed direction is " << currentTraj.seed().direction() + << ".Traj direction is " << currentTraj.direction(); + + //updating MultiRecHits and fit-smooth again + std::pair curiterationhits = + updateHits(currentTraj,updator,*ian); + currentTraj = fit(curiterationhits, theFitter, currentTraj); + + //saving trajectory for each annealing cycle ... + if(TrajAnnSaving_){ + TrajAnnealing temp(currentTraj, *ian); + trajann.push_back(temp); + } + + LogDebug("DAFTrackProducerAlgorithm") << "done annealing value " << (*ian) ; + + } + else break; + + + } //end of annealing program + + LogDebug("DAFTrackProducerAlgorithm") << "Ended annealing program with " << (1.*checkHits(*ivtraj, currentTraj))/(1.*(*ivtraj).measurements().size())*100. << " unchanged." << std::endl; + + //computing the ndof keeping into account the weights + ndof = calculateNdof(currentTraj); + + //checking if the trajectory has the minimum number of valid hits ( weight (>1e-6) ) + //in order to remove tracks with too many outliers. + + //std::vector filtered; + //filter(theFitter, vtraj, minHits_, filtered, builder); + + if(currentTraj.foundHits() >= minHits_) { + + bool ok = buildTrack(currentTraj, algoResults, ndof, bs) ; + if(ok) cont++; + + } + else{ + LogDebug("DAFTrackProducerAlgorithm") << "Rejecting trajectory with " + << currentTraj.foundHits()<<" hits"; + } + } + else + LogDebug("DAFTrackProducerAlgorithm") << "Rejecting empty trajectory" << std::endl; + + } //end run on track collection + + LogDebug("DAFTrackProducerAlgorithm") << "Number of Tracks found: " << cont << "\n"; + +} +/*------------------------------------------------------------------------------------------------------*/ +std::pair +DAFTrackProducerAlgorithm::collectHits(const Trajectory vtraj, + const MultiRecHitCollector* measurementCollector, + const MeasurementTrackerEvent *measTk) const +{ + + LogDebug("DAFTrackProducerAlgorithm") << "Calling DAFTrackProducerAlgorithm::collectHits"; + + //getting the traj measurements from the MeasurementCollector + int nHits = 0; + TransientTrackingRecHit::RecHitContainer hits; + std::vector collectedmeas = measurementCollector->recHits(vtraj, measTk); + + //if the MeasurementCollector is empty, make an "empty" pair + //else taking the collected measured hits and building the pair + if( collectedmeas.empty() ) + return std::make_pair(TransientTrackingRecHit::RecHitContainer(), TrajectoryStateOnSurface()); + + for( std::vector::const_iterator iter = collectedmeas.begin(); + iter!=collectedmeas.end(); iter++ ){ + + nHits++; + hits.push_back(iter->recHit()); + + } + + + //TrajectoryStateWithArbitraryError() == Creates a TrajectoryState with the same parameters + // as the input one, but with "infinite" errors, i.e. errors so big that they don't + // bias a fit starting with this state. + //return std::make_pair(hits,TrajectoryStateWithArbitraryError()(collectedmeas.front().predictedState())); + + // I do not have to rescale the error because it is already rescaled in the fit code + TrajectoryStateOnSurface initialStateFromTrack = collectedmeas.front().predictedState(); + + return std::make_pair(hits, initialStateFromTrack); + +} +/*------------------------------------------------------------------------------------------------------*/ +std::pair +DAFTrackProducerAlgorithm::updateHits(const Trajectory vtraj, + const SiTrackerMultiRecHitUpdator* updator, + double annealing) const +{ + TransientTrackingRecHit::RecHitContainer hits; + std::vector vmeas = vtraj.measurements(); + std::vector::reverse_iterator imeas; + + //I run inversely on the trajectory obtained and update the state + for (imeas = vmeas.rbegin(); imeas != vmeas.rend(); imeas++){ + TransientTrackingRecHit::RecHitPointer updated = updator->update(imeas->recHit(), + imeas->updatedState(), annealing); + hits.push_back(updated); + } + + TrajectoryStateOnSurface updatedStateFromTrack = vmeas.back().predictedState(); + //updatedStateFromTrack.rescaleError(10); + + //return std::make_pair(hits,TrajectoryStateWithArbitraryError()(vmeas.back().updatedState())); + return std::make_pair(hits,updatedStateFromTrack); +} +/*------------------------------------------------------------------------------------------------------*/ +Trajectory DAFTrackProducerAlgorithm::fit(const std::pair& hits, + const TrajectoryFitter * theFitter, + Trajectory vtraj) const { + + //creating a new trajectory starting from the direction of the seed of the input one and the hits + Trajectory newVec = theFitter->fitOne(TrajectorySeed(PTrajectoryStateOnDet(), + BasicTrajectorySeed::recHitContainer(), + vtraj.seed().direction()), + hits.first, hits.second); + + if( newVec.isValid() ) return newVec; + else{ + LogDebug("DAFTrackProducerAlgorithm") << "Fit no valid."; + return Trajectory(); + } + +} +/*------------------------------------------------------------------------------------------------------*/ +bool DAFTrackProducerAlgorithm::buildTrack (const Trajectory vtraj, + AlgoProductCollection& algoResults, + float ndof, + const reco::BeamSpot& bs) const +{ + //variable declarations + reco::Track * theTrack; + Trajectory * theTraj; + + LogDebug("DAFTrackProducerAlgorithm") <<" BUILDER " << std::endl;; + TrajectoryStateOnSurface innertsos; + + if ( vtraj.isValid() ){ + + theTraj = new Trajectory( vtraj ); + + if (vtraj.direction() == alongMomentum) { + //if (theTraj->direction() == oppositeToMomentum) { + innertsos = vtraj.firstMeasurement().updatedState(); + } else { + innertsos = vtraj.lastMeasurement().updatedState(); + } + + TSCBLBuilderNoMaterial tscblBuilder; + TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(*(innertsos.freeState()),bs); + + if (tscbl.isValid()==false) return false; + + GlobalPoint v = tscbl.trackStateAtPCA().position(); + math::XYZPoint pos( v.x(), v.y(), v.z() ); + GlobalVector p = tscbl.trackStateAtPCA().momentum(); + math::XYZVector mom( p.x(), p.y(), p.z() ); + + // LogDebug("TrackProducer") <& input, + int minhits, std::vector& output, + const TransientTrackingRecHitBuilder* builder) const +{ + if (input.empty()) return; + + int ngoodhits = 0; + std::vector vtm = input[0].measurements(); + TransientTrackingRecHit::RecHitContainer hits; + + //count the number of non-outlier and non-invalid hits + for (std::vector::reverse_iterator tm=vtm.rbegin(); tm!=vtm.rend();tm++){ + //if the rechit is valid + if (tm->recHit()->isValid()) { + SiTrackerMultiRecHit const & mHit = dynamic_cast(*tm->recHit()); + std::vector components = mHit.recHits(); + int iComp = 0; + bool isGood = false; + for(std::vector::const_iterator iter = components.begin(); iter != components.end(); iter++, iComp++){ + //if there is at least one component with weight higher than 1e-6 then the hit is not an outlier + if (mHit.weight(iComp)>1e-6) {ngoodhits++; iComp++; isGood = true; break;} + } + if (isGood) { + TkClonerImpl hc = static_cast(builder)->cloner(); + auto tempHit = hc.makeShared(tm->recHit(),tm->updatedState()); + hits.push_back(tempHit); + } + else hits.push_back(std::make_shared(*tm->recHit()->det(), TrackingRecHit::missing)); + } else { + TkClonerImpl hc = static_cast(builder)->cloner(); + auto tempHit = hc.makeShared(tm->recHit(),tm->updatedState()); + hits.push_back(tempHit); + } + } + + + LogDebug("DAFTrackProducerAlgorithm") << "Original number of valid hits " << input[0].foundHits() << "; after filtering " << ngoodhits; + //debug + if (ngoodhits>input[0].foundHits()) edm::LogError("DAFTrackProducerAlgorithm") << "Something wrong: the number of good hits from DAFTrackProducerAlgorithm::filter " << ngoodhits << " is higher than the original one " << input[0].foundHits(); + + if (ngoodhits < minhits) return; + + TrajectoryStateOnSurface curstartingTSOS = input.front().lastMeasurement().updatedState(); + LogDebug("DAFTrackProducerAlgorithm") << "starting tsos for final refitting " << curstartingTSOS ; + //curstartingTSOS.rescaleError(100); + + output = fitter->fit(TrajectorySeed(PTrajectoryStateOnDet(), + BasicTrajectorySeed::recHitContainer(), + input.front().seed().direction()), + hits, + TrajectoryStateWithArbitraryError()(curstartingTSOS)); + LogDebug("DAFTrackProducerAlgorithm") << "After filtering " << output.size() << " trajectories"; + +} +/*------------------------------------------------------------------------------------------------------*/ +float DAFTrackProducerAlgorithm::calculateNdof(const Trajectory vtraj) const +{ + + if (!vtraj.isValid()) return 0; + float ndof = 0; + const std::vector& meas = vtraj.measurements(); + for (std::vector::const_iterator iter = meas.begin(); iter != meas.end(); iter++){ + if (iter->recHit()->isValid()){ + SiTrackerMultiRecHit const & mHit = dynamic_cast(*iter->recHit()); + std::vector components = mHit.recHits(); + int iComp = 0; + for(std::vector::const_iterator iter2 = components.begin(); iter2 != components.end(); iter2++, iComp++){ + if ((*iter2)->isValid()) + ndof += ((*iter2)->dimension())*mHit.weight(iComp); + } + + } + } + + return ndof-5; + +} +//------------------------------------------------------------------------------------------------ +int DAFTrackProducerAlgorithm::checkHits( Trajectory iInitTraj, const Trajectory iFinalTraj) const { + + std::vector initmeasurements = iInitTraj.measurements(); + std::vector finalmeasurements = iFinalTraj.measurements(); + std::vector::iterator imeas; + std::vector::iterator jmeas; + int nSame = 0; + int ihit = 0; + + if( initmeasurements.empty() || finalmeasurements.empty() ){ + LogDebug("DAFTrackProducerAlgorithm") << "Initial or Final Trajectory empty."; + return 0; + } + + if( initmeasurements.size() != finalmeasurements.size() ) { + LogDebug("DAFTrackProducerAlgorithm") << "Initial and Final Trajectory have different size."; + return 0; + } + + for(jmeas = initmeasurements.begin(); jmeas != initmeasurements.end(); jmeas++){ + + const TrackingRecHit* initHit = jmeas->recHit()->hit(); + + if(!initHit->isValid() && ihit == 0 ) continue; + + if(initHit->isValid()){ + + TrajectoryMeasurement imeas = finalmeasurements.at(ihit); + const TrackingRecHit* finalHit = imeas.recHit()->hit(); + const TrackingRecHit* MaxWeightHit=0; + float maxweight = 0; + + const SiTrackerMultiRecHit* mrh = dynamic_cast(finalHit); + if (mrh){ + std::vector components = mrh->recHits(); + std::vector::const_iterator icomp; + int hitcounter=0; + + for (icomp = components.begin(); icomp != components.end(); icomp++) { + if((*icomp)->isValid()) { + double weight = mrh->weight(hitcounter); + if(weight > maxweight) { + MaxWeightHit = *icomp; + maxweight = weight; + } + } + hitcounter++; + } + + auto myref1 = reinterpret_cast(initHit)->firstClusterRef(); + auto myref2 = reinterpret_cast(MaxWeightHit)->firstClusterRef(); + + if( myref1 == myref2 ){ + nSame++; + } + } + } else { + + TrajectoryMeasurement imeas = finalmeasurements.at(ihit); + const TrackingRecHit* finalHit = imeas.recHit()->hit(); + if(!finalHit->isValid()){ + nSame++; + } + } + + ihit++; + } + + return nSame; +} diff --git a/RecoTracker/TransientTrackingRecHit/interface/TSiTrackerMultiRecHit.h b/RecoTracker/TransientTrackingRecHit/interface/TSiTrackerMultiRecHit.h deleted file mode 100644 index 1ece4d2d53f1b..0000000000000 --- a/RecoTracker/TransientTrackingRecHit/interface/TSiTrackerMultiRecHit.h +++ /dev/null @@ -1,94 +0,0 @@ -#ifndef TSiTrackerMultiRecHit_h -#define TSiTrackerMultiRecHit_h - -#include "TrackingTools/TransientTrackingRecHit/interface/TValidTrackingRecHit.h" -#include "DataFormats/TrackerRecHit2D/interface/SiTrackerMultiRecHit.h" -#include "TrackingTools/TransientTrackingRecHit/interface/HelpertRecHit2DLocalPos.h" -#include "Geometry/CommonDetUnit/interface/GeomDetUnit.h" - -/* -A TransientTrackingRecHit for the SiTrackerMultiRecHit -*/ - -class TSiTrackerMultiRecHit GCC11_FINAL : public TValidTrackingRecHit { -public: - //virtual ~TSiTrackerMultiRecHit() {delete theHitData;} - virtual ~TSiTrackerMultiRecHit() {} - - virtual AlgebraicVector parameters() const {return theHitData.parameters();} - virtual AlgebraicSymMatrix parametersError() const { - return HelpertRecHit2DLocalPos().parError( theHitData.localPositionError(), *det()); - //return theHitData.parametersError(); - } - - virtual void getKfComponents( KfComponentsHolder & holder ) const { - HelpertRecHit2DLocalPos().getKfComponents(holder, theHitData, *det()); - } - virtual DetId geographicalId() const {return theHitData.geographicalId();} - virtual AlgebraicMatrix projectionMatrix() const {return theHitData.projectionMatrix();} - virtual int dimension() const {return theHitData.dimension();} - - virtual LocalPoint localPosition() const {return theHitData.localPosition();} - virtual LocalError localPositionError() const {return theHitData.localPositionError();} - - virtual const TrackingRecHit * hit() const {return &theHitData;}; - const SiTrackerMultiRecHit* specificHit() const {return &theHitData;} - - virtual bool isValid() const{return theHitData.isValid();} - - virtual std::vector recHits() const { - return theHitData.recHits(); - } - virtual std::vector recHits() { - return theHitData.recHits(); - } - - - /// interface needed to set and read back an annealing value that has been applied to the current hit error matrix when - /// using it as a component for a composite rec hit (useful for the DAF) - void setAnnealingFactor(float annealing) {annealing_ = annealing;} - float getAnnealingFactor() const {return annealing_;} - - - //vector of weights - std::vector const & weights() const {return theHitData.weights();} - std::vector & weights() {return theHitData.weights();} - - //returns the weight for the i component - float weight(unsigned int i) const {return theHitData.weight(i);} - float & weight(unsigned int i) {return theHitData.weight(i);} - - - virtual const GeomDetUnit* detUnit() const; - - virtual bool canImproveWithTrack() const {return true;} - - virtual RecHitPointer clone(const TrajectoryStateOnSurface& ts) const; - - virtual ConstRecHitContainer transientHits() const {return theComponents;}; - - static RecHitPointer build( const GeomDet * geom, const SiTrackerMultiRecHit* rh, - const ConstRecHitContainer& components, float annealing=1.){ - return RecHitPointer(new TSiTrackerMultiRecHit( geom, rh, components, annealing)); - } - - - -private: - SiTrackerMultiRecHit theHitData; - //holds the TransientTrackingRecHit components of the MultiRecHit - ConstRecHitContainer theComponents; - float annealing_; - - TSiTrackerMultiRecHit(const GeomDet * geom, const SiTrackerMultiRecHit* rh, - const ConstRecHitContainer& components, float annealing): - TValidTrackingRecHit(*geom), theHitData(*rh), theComponents(components), annealing_(annealing){} - - virtual TSiTrackerMultiRecHit* clone() const { - return new TSiTrackerMultiRecHit(*this); - } - - -}; - -#endif diff --git a/RecoTracker/TransientTrackingRecHit/src/TSiTrackerMultiRecHit.cc b/RecoTracker/TransientTrackingRecHit/src/TSiTrackerMultiRecHit.cc deleted file mode 100644 index 6a0d0b134162e..0000000000000 --- a/RecoTracker/TransientTrackingRecHit/src/TSiTrackerMultiRecHit.cc +++ /dev/null @@ -1,36 +0,0 @@ -#include "RecoTracker/TransientTrackingRecHit/interface/TSiTrackerMultiRecHit.h" - -/* -TSiTrackerMultiRecHit::TSiTrackerMultiRecHit(const GeomDet * geom, const std::vector& rhs, const SiTrackerMultiRecHitUpdator* updator, const TrajectoryStateOnSurface& tsos): -TValidTrackingRecHit(geom), theUpdator(updator){ - theHitData = theUpdator->buildMultiRecHit(tsos, rhs, theComponents); - setAnnealingFactor(theUpdator->getCurrentAnnealingValue()); -} -*/ - -const GeomDetUnit* TSiTrackerMultiRecHit::detUnit() const{ - return dynamic_cast(det()); -} - -TransientTrackingRecHit::RecHitPointer TSiTrackerMultiRecHit::clone(const TrajectoryStateOnSurface& ts) const{ -/* - std::vector updatedcomponents = theComponents; - SiTrackerMultiRecHit better = theUpdator->update(ts,&theHitData, updatedcomponents); - RecHitPointer result = TSiTrackerMultiRecHit::build( det(), &better, theUpdator, updatedcomponents ); - return result; -*/ - return RecHitPointer(this->clone()); -} - -/* - -std::vector TSiTrackerMultiRecHit::recHits() const { - std::vector components; - std::vector::const_iterator iter; - for (iter = theComponents.begin(); iter != theComponents.end(); iter++){ - components.push_back(iter->get()); - } - return components; -} - -*/ diff --git a/TrackingTools/KalmanUpdators/interface/MRHChi2MeasurementEstimator.h b/TrackingTools/KalmanUpdators/interface/MRHChi2MeasurementEstimator.h index 4e3052fbaf395..7836d8bd00ec5 100644 --- a/TrackingTools/KalmanUpdators/interface/MRHChi2MeasurementEstimator.h +++ b/TrackingTools/KalmanUpdators/interface/MRHChi2MeasurementEstimator.h @@ -17,8 +17,12 @@ class MRHChi2MeasurementEstimator GCC11_FINAL : public Chi2MeasurementEstimatorB explicit MRHChi2MeasurementEstimator(double maxChi2, double nSigma = 3.) : Chi2MeasurementEstimatorBase( maxChi2, nSigma) {} - virtual std::pair estimate(const TrajectoryStateOnSurface&, - const TrackingRecHit&) const; + std::pair estimate(const TrajectoryStateOnSurface& tsos, + const TrackingRecHit& aRecHit) const ; + template + std::pair estimate(const TrajectoryStateOnSurface& tsos, + const TrackingRecHit& aRecHit) const ; + virtual MRHChi2MeasurementEstimator* clone() const { return new MRHChi2MeasurementEstimator(*this); diff --git a/TrackingTools/KalmanUpdators/interface/TrackingRecHitPropagator.h b/TrackingTools/KalmanUpdators/interface/TrackingRecHitPropagator.h index 7e080291469a2..536d9da400f86 100644 --- a/TrackingTools/KalmanUpdators/interface/TrackingRecHitPropagator.h +++ b/TrackingTools/KalmanUpdators/interface/TrackingRecHitPropagator.h @@ -1,12 +1,16 @@ #ifndef TrackingTools_TrackingRecHitPropagator_h #define TrackingTools_TrackingRecHitPropagator_h +#include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h" +#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h" +#include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h" +#include "RecoTracker/TransientTrackingRecHit/interface/TkClonerImpl.h" + #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h" #include "TrackingTools/KalmanUpdators/interface/TrackingRecHitPropagator.h" #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h" #include "MagneticField/Engine/interface/MagneticField.h" #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" -#include "TrackingTools/TransientTrackingRecHit/interface/InvalidTransientRecHit.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" @@ -23,25 +27,29 @@ class TrackingRecHitPropagator { ~TrackingRecHitPropagator() {delete thePropagator;} - /* + template TrackingRecHit::RecHitPointer project(const TrackingRecHit::ConstRecHitPointer hit, const GeomDet& det, - const TrajectoryStateOnSurface ts) const{ + const TrajectoryStateOnSurface ts, + const TransientTrackingRecHitBuilder* builder) const{ + + TkClonerImpl hc = static_cast(builder)->cloner(); + //1) propagate the best possible track parameters to the surface of the hit you want to "move" using a AnalyticalPropagator ; //2) create LocalTrajectoryParameters with the local x,y of the hit and direction + momentum from the propagated track parameters; //3) create a LocalTrajectoryError matrix which is 0 except for the local x,y submatrix, which is filled with the hit errors; //4) create a TSOS from the result of 2) and 3) and propagate it to the reference surface; //5) create a new hit with the local x,y subspace of the result of 4) - if (!ts.isValid()) return InvalidTransientRecHit::build(hit->det()); + if (!ts.isValid()) return std::make_shared(*hit->det(), TrackingRecHit::missing); // LogTrace("SiTrackerMultiRecHitUpdator") << "the tsos is valid"; //check if the ts lays or not on the destination surface and in case propagate it TrajectoryStateOnSurface propagated =ts; if (hit->surface() != &(ts.surface())) propagated = thePropagator->propagate(ts, *(hit->surface())); - if (!propagated.isValid()) return InvalidTransientRecHit::build(hit->det()); + if (!propagated.isValid()) return std::make_shared(*hit->det(), TrackingRecHit::missing); // LogTrace("SiTrackerMultiRecHitUpdator") << "the propagate tsos is valid"; // LogTrace("SiTrackerMultiRecHitUpdator") << "Original: position: "<parameters()<<" error: "<parametersError()<clone(propagated); + TrackingRecHit::RecHitPointer updatedOriginal = hc.makeShared(hit,propagated); // LogTrace("SiTrackerMultiRecHitUpdator") << "New: position: "<parameters()<<" error: "<parametersError()<propagate(hit_state, det.surface()); - if (!projected_hit_state.isValid()) return InvalidTransientRecHit::build(hit->det()); + if (!projected_hit_state.isValid()) return std::make_shared(*hit->det(), TrackingRecHit::missing); LocalPoint p = projected_hit_state.localPosition(); LocalError e = projected_hit_state.localError().positionError(); // LogTrace("SiTrackerMultiRecHitUpdator") << "position: "<det(), updatedOriginal, this); } - */ + private: const AnalyticalPropagator* thePropagator; }; diff --git a/TrackingTools/KalmanUpdators/src/Chi2MeasurementEstimator.cc b/TrackingTools/KalmanUpdators/src/Chi2MeasurementEstimator.cc index c5513130df7b7..b91841c3c8d5a 100644 --- a/TrackingTools/KalmanUpdators/src/Chi2MeasurementEstimator.cc +++ b/TrackingTools/KalmanUpdators/src/Chi2MeasurementEstimator.cc @@ -21,6 +21,7 @@ Chi2MeasurementEstimator::estimate(const TrajectoryStateOnSurface& tsos, template std::pair Chi2MeasurementEstimator::estimate(const TrajectoryStateOnSurface& tsos, const TrackingRecHit& aRecHit) const { + typedef typename AlgebraicROOTObject::Matrix MatD5; typedef typename AlgebraicROOTObject<5,D>::Matrix Mat5D; typedef typename AlgebraicROOTObject::SymMatrix SMatDD; @@ -37,5 +38,6 @@ Chi2MeasurementEstimator::estimate(const TrajectoryStateOnSurface& tsos, R += RMeas; invertPosDefMatrix(R); double est = ROOT::Math::Similarity(r - rMeas, R); + return returnIt(est); } diff --git a/TrackingTools/KalmanUpdators/src/MRHChi2MeasurementEstimator.cc b/TrackingTools/KalmanUpdators/src/MRHChi2MeasurementEstimator.cc index 627d49ccb16d0..96c8b2903a75b 100644 --- a/TrackingTools/KalmanUpdators/src/MRHChi2MeasurementEstimator.cc +++ b/TrackingTools/KalmanUpdators/src/MRHChi2MeasurementEstimator.cc @@ -1,36 +1,74 @@ -#include "RecoTracker/TransientTrackingRecHit/interface/TSiTrackerMultiRecHit.h" +#include "DataFormats/TrackerRecHit2D/interface/SiTrackerMultiRecHit.h" +#include "RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdator.h" #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" #include "TrackingTools/KalmanUpdators/interface/MRHChi2MeasurementEstimator.h" #include "TrackingTools/PatternTools/interface/MeasurementExtractor.h" #include "DataFormats/GeometrySurface/interface/Plane.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "DataFormats/Math/interface/invertPosDefMatrix.h" +#include "DataFormats/Math/interface/ProjectMatrix.h" +#include "DataFormats/TrackingRecHit/interface/KfComponentsHolder.h" -std::pair -MRHChi2MeasurementEstimator::estimate(const TrajectoryStateOnSurface& tsos, - const TrackingRecHit& aRecHit) const { - if (!aRecHit.isValid()) { - throw cms::Exception("MRHChi2MeasurementEstimator") << "Invalid RecHit passed to the MRHChi2MeasurementEstimator "; - } +std::pair MRHChi2MeasurementEstimator::estimate(const TrajectoryStateOnSurface& tsos, + const TrackingRecHit& aRecHit) const { - typedef AlgebraicROOTObject<2>::Vector Vec; - typedef AlgebraicROOTObject<2>::SymMatrix Mat; + switch (aRecHit.dimension()) { + case 2: return estimate<2>(tsos,aRecHit); + //avoid the not-2D hit due to the final sum + case ( 1 || 3 || 4 || 5 ):{ + LogDebug("MRHChi2MeasurementEstimator") << "WARNING:The hit is not 2D: does not count in the MRH Chi2 estimation." ; + double est = 0.0; + return HitReturnType(false, est); + } + } + throw cms::Exception("Rec hit of invalid dimension (not 1,2,3,4,5)") << + "The value was " << aRecHit.dimension() << + ", type is " << typeid(aRecHit).name() << "\n"; +} - //better be a multihit... - TSiTrackerMultiRecHit const & mHit = dynamic_cast(aRecHit); - MeasurementExtractor me(tsos); +//--------------------------------------------------------------------------------------------------------------- +template +std::pair MRHChi2MeasurementEstimator::estimate(const TrajectoryStateOnSurface& tsos, + const TrackingRecHit& aRecHit) const { + + SiTrackerMultiRecHit const & mHit = dynamic_cast(aRecHit); double est=0; + double annealing = mHit.getAnnealingFactor(); - LogDebug("MRHChi2MeasurementEstimator") << "Current annealing factor is " << annealing; - TransientTrackingRecHit::ConstRecHitContainer components = mHit.transientHits(); - LogDebug("MRHChi2MeasurementEstimator") << "this hit has " << components.size() << " components"; - for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator iter = components.begin(); iter != components.end(); iter++){ - Vec r = asSVector<2>((*iter)->parameters()) - me.measuredParameters<2>(**iter); - Mat R = asSMatrix<2>((*iter)->parametersError())*annealing + me.measuredError<2>(**iter); - //int ierr = ! R.Invert(); // if (ierr != 0) throw exception; // - invertPosDefMatrix(R); - LogDebug("MRHChi2MeasurementEstimator") << "Hit with weight " << (*iter)->weight(); - est += ROOT::Math::Similarity(r, R)*((*iter)->weight()); - } + LogDebug("MRHChi2MeasurementEstimator") << "Current annealing factor is " << annealing; + + std::vector components = mHit.recHits(); + LogDebug("MRHChi2MeasurementEstimator") << "this hit has " << components.size() << " components"; + + int iComp = 0; + for(std::vector::const_iterator iter = components.begin(); iter != components.end(); iter++, iComp++){ + + // define variables that will be used to setup the KfComponentsHolder + ProjectMatrix pf; + typename AlgebraicROOTObject::Matrix H; + typename AlgebraicROOTObject::Vector r, rMeas; + typename AlgebraicROOTObject::SymMatrix V, VMeas; + AlgebraicVector5 x = tsos.localParameters().vector(); + const AlgebraicSymMatrix55 &C = (tsos.localError().matrix()); + + // setup the holder with the correct dimensions and get the values + KfComponentsHolder holder; + holder.template setup(&r, &V, &H, &pf, &rMeas, &VMeas, x, C); + (**iter).getKfComponents(holder); + + r -= rMeas; + V = V*annealing + VMeas; + bool ierr = invertPosDefMatrix(V); + if( !ierr ) { + edm::LogError("SiTrackerMultiRecHitUpdator")<<"SiTrackerMultiRecHitUpdator::ComputeParameters2dim: W not valid!"< + diff --git a/TrackingTools/PatternTools/interface/TrajAnnealing.h b/TrackingTools/PatternTools/interface/TrajAnnealing.h new file mode 100644 index 0000000000000..8ae4fb9f4d6f9 --- /dev/null +++ b/TrackingTools/PatternTools/interface/TrajAnnealing.h @@ -0,0 +1,44 @@ +#ifndef DataFormats_TrajAnnealing_h +#define DataFormats_TrajAnnealing_h + +#include + +/** This class allow to save all the traj info + * for each annealing cycle + * suitable for algorithm like DAF + * (not used in default options) + */ + +class TrajAnnealing +{ +public: + TrajAnnealing(): + traj_(), + annealing_(0), + theWeights(){} + virtual ~TrajAnnealing(){} + + TrajAnnealing( const Trajectory&, float ); + + float getAnnealing() const { return annealing_; } + Trajectory getTraj() const { return traj_; } + + //vector of weights + std::vector const & weights() const {return theWeights;} + std::vector & weights() {return theWeights;} + +private: + + Trajectory traj_; + float annealing_; + std::vector theWeights; + TrackingRecHit::RecHitContainer theHits_; + + std::pair > getAnnealingWeight( const TrackingRecHit& aRecHit ) const ; +}; + +// this is our new product, it is simply a +// collection of TrajAnnealing held in an std::vector +typedef std::vector TrajAnnealingCollection; + +#endif diff --git a/TrackingTools/PatternTools/src/TrajAnnealing.cc b/TrackingTools/PatternTools/src/TrajAnnealing.cc new file mode 100644 index 0000000000000..b3624c45cf50e --- /dev/null +++ b/TrackingTools/PatternTools/src/TrajAnnealing.cc @@ -0,0 +1,53 @@ +#include "TrackingTools/PatternTools/interface/Trajectory.h" +#include "TrackingTools/PatternTools/interface/TrajAnnealing.h" +#include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h" +#include "DataFormats/TrackerRecHit2D/interface/SiTrackerMultiRecHit.h" + +TrajAnnealing::TrajAnnealing( const Trajectory& traj, float ann ){ + traj_ = traj; + annealing_ = ann; + + const Trajectory::DataContainer& measurements = traj_.measurements(); + if( measurements.size() > 2 ){ + + Trajectory::DataContainer::const_iterator ibegin,iend; + int increment(0); + if( traj.direction() == alongMomentum ){ + ibegin = measurements.begin(); + iend = measurements.end(); + increment = 1; + } else { + ibegin = measurements.end(); + iend = measurements.begin(); + increment = -1; + } + + for( Trajectory::DataContainer::const_iterator imeas = ibegin; imeas != iend; imeas += increment ){ + theHits_.push_back(imeas->recHit()); + if (imeas->recHit()->isValid()){ + SiTrackerMultiRecHit const & mHit = dynamic_cast(*imeas->recHit()); + std::vector components = mHit.recHits(); + int iComp = 0; + for(std::vector::const_iterator iter = components.begin(); + iter != components.end(); iter++, iComp++){ + theWeights.push_back(mHit.weight(iComp)); + } + } + } + + } + +} + +std::pair > TrajAnnealing::getAnnealingWeight( const TrackingRecHit& aRecHit ) const { + + if (!aRecHit.isValid()) { + std::vector dumpyVec = {0.0}; + return make_pair(0.0,dumpyVec); + } + + SiTrackerMultiRecHit const & mHit = dynamic_cast(aRecHit); + return make_pair(mHit.getAnnealingFactor(), mHit.weights()); + +} + diff --git a/TrackingTools/PatternTools/src/classes.h b/TrackingTools/PatternTools/src/classes.h index 5565695f5ec20..f237681d1607d 100644 --- a/TrackingTools/PatternTools/src/classes.h +++ b/TrackingTools/PatternTools/src/classes.h @@ -9,6 +9,7 @@ #include "DataFormats/TrackCandidate/interface/TrackCandidate.h" #include "DataFormats/TrackReco/interface/TrackFwd.h" #include "TrackingTools/PatternTools/interface/Trajectory.h" +#include "TrackingTools/PatternTools/interface/TrajAnnealing.h" #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h" #include "DataFormats/GeometrySurface/interface/Surface.h" #include "DataFormats/CLHEP/interface/Migration.h" @@ -29,7 +30,10 @@ namespace TrackingTools_PatternTools { struct dictionary { std::vector kk; edm::Wrapper > trajCollWrapper; - + + std::vector ta; + edm::Wrapper > trajAnnCollaction; + TrajTrackAssociationCollection ttam; edm::Wrapper wttam; TrajTrackAssociation vttam; diff --git a/TrackingTools/PatternTools/src/classes_def.xml b/TrackingTools/PatternTools/src/classes_def.xml index 1ff1f97cc8b96..ee1f939238f25 100755 --- a/TrackingTools/PatternTools/src/classes_def.xml +++ b/TrackingTools/PatternTools/src/classes_def.xml @@ -59,4 +59,8 @@ + + + + diff --git a/TrackingTools/TrackFitters/python/MRHFitters_cff.py b/TrackingTools/TrackFitters/python/MRHFitters_cff.py index 27e7b3d600af5..ecba177950584 100644 --- a/TrackingTools/TrackFitters/python/MRHFitters_cff.py +++ b/TrackingTools/TrackFitters/python/MRHFitters_cff.py @@ -1,5 +1,7 @@ import FWCore.ParameterSet.Config as cms +from TrackingTools.KalmanUpdators.MRHChi2MeasurementEstimatorESProducer_cfi import * + import TrackingTools.TrackFitters.KFTrajectoryFitterESProducer_cfi MRHTrajectoryFitter = TrackingTools.TrackFitters.KFTrajectoryFitterESProducer_cfi.KFTrajectoryFitter.clone( ComponentName = 'MRHFitter', diff --git a/TrackingTools/TrackFitters/python/TrackFitters_cff.py b/TrackingTools/TrackFitters/python/TrackFitters_cff.py index 6ee71422391eb..a2a9b37952606 100644 --- a/TrackingTools/TrackFitters/python/TrackFitters_cff.py +++ b/TrackingTools/TrackFitters/python/TrackFitters_cff.py @@ -5,7 +5,7 @@ from TrackingTools.TrackFitters.KFTrajectoryFitterESProducer_cfi import * from TrackingTools.TrackFitters.KFTrajectorySmootherESProducer_cfi import * -# from TrackingTools.TrackFitters.MRHFitters_cff import * +from TrackingTools.TrackFitters.MRHFitters_cff import * from TrackingTools.TrackFitters.RungeKuttaFitters_cff import * from TrackingTools.TrackFitters.RungeKutta1DFitters_cff import * diff --git a/TrackingTools/TrackFitters/src/KFTrajectoryFitter.cc b/TrackingTools/TrackFitters/src/KFTrajectoryFitter.cc index 7e3ed31e1c951..d7116e28b6403 100644 --- a/TrackingTools/TrackFitters/src/KFTrajectoryFitter.cc +++ b/TrackingTools/TrackFitters/src/KFTrajectoryFitter.cc @@ -86,17 +86,15 @@ Trajectory KFTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, // if unlikely(hit.det() == nullptr) continue; if unlikely( (!hit.isValid()) && hit.surface() == nullptr) { - std::cout << "TrackFitters" << " Error: invalid hit with no GeomDet attached .... skipping" << std::endl; LogDebug("TrackFitters")<< " Error: invalid hit with no GeomDet attached .... skipping"; continue; } - if (hit.det() && hit.geographicalId()<1000U) std::cout << "Problem 0 det id for " << typeid(hit).name() << ' ' << hit.det()->geographicalId() << std::endl; - if (hit.isValid() && hit.geographicalId()<1000U) std::cout << "Problem 0 det id for " << typeid(hit).name() << ' ' << hit.det()->geographicalId() << std::endl; + //if (hit.det() && hit.geographicalId()<1000U) LogDebug("TrackFitters")<< "Problem 0 det id for " << typeid(hit).name() << ' ' << hit.det()->geographicalId() ; + //if (hit.isValid() && hit.geographicalId()<1000U) LogDebug("TrackFitters")<< "Problem 0 det id for " << typeid(hit).name() << ' ' << hit.det()->geographicalId(); #ifdef EDM_ML_DEBUG if (hit.isValid()) { - LogTrace("TrackFitters") - << " ----------------- HIT #" << hitcounter << " (VALID)-----------------------\n" + LogTrace("TrackFitters")<< " ----------------- 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" @@ -197,11 +195,17 @@ Trajectory KFTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, || std::abs(currTsos.localParameters().position().x()) > 1000 ) ) || edm::isNotFinite(currTsos.localParameters().qbp()); if unlikely(badState){ - if (!currTsos.isValid()) edm::LogError("FailedUpdate") - <<"updating with the hit failed. Not updating the trajectory with the hit"; - else if (edm::isNotFinite(currTsos.localParameters().qbp())) edm::LogError("TrajectoryNaN")<<"Trajectory has NaN"; - else LogTrace("FailedUpdate")<<"updated state is valid but pretty bad, skipping. currTsos " - <idToLayer((*ihit)->geographicalId()) )); //There is a no-fail policy here. So, it's time to give up //Keep the traj with invalid TSOS so that it's clear what happened @@ -213,11 +217,15 @@ Trajectory KFTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, return Trajectory(); } } else{ - if (preciseHit->det()) myTraj.push(TM(predTsos, currTsos, preciseHit, + if (preciseHit->det()){ + myTraj.push(TM(predTsos, currTsos, preciseHit, estimator()->estimate(predTsos, *preciseHit).second, theGeometry->idToLayer(preciseHit->geographicalId()) )); - else myTraj.push(TM(predTsos, currTsos, preciseHit, + } + else{ + myTraj.push(TM(predTsos, currTsos, preciseHit, estimator()->estimate(predTsos, *preciseHit).second)); + } } } } else { @@ -228,7 +236,6 @@ Trajectory KFTrajectoryFitter::fitOne(const TrajectorySeed& aSeed, if ((*ihit)->det()) myTraj.push(TM(predTsos, *ihit,0,theGeometry->idToLayer((*ihit)->geographicalId()) )); else myTraj.push(TM(predTsos, *ihit,0)); } - LogTrace("TrackFitters") << "predTsos !" << "\n" << predTsos << "\n" diff --git a/TrackingTools/TrackFitters/src/KFTrajectorySmoother.cc b/TrackingTools/TrackFitters/src/KFTrajectorySmoother.cc index 18d1cfae4670d..31466b0802545 100644 --- a/TrackingTools/TrackFitters/src/KFTrajectorySmoother.cc +++ b/TrackingTools/TrackFitters/src/KFTrajectorySmoother.cc @@ -78,12 +78,10 @@ KFTrajectorySmoother::trajectory(const Trajectory& aTraj) const { //check surface just for safety: should never be ==0 because they are skipped in the fitter // if unlikely(hit->det() == nullptr) continue; if unlikely( hit->surface()==nullptr ) { - LogDebug("TrackFitters")<< " Error: invalid hit with no GeomDet attached .... skipping"; + LogDebug("TrackFitters") << " Error: invalid hit with no GeomDet attached .... skipping"; continue; } - if (hit->det() && hit->geographicalId()<1000U) std::cout << "Problem 0 det id for " << typeid(*hit).name() << ' ' << hit->det()->geographicalId() << std::endl; - if (hit->isValid() && hit->geographicalId()<1000U) std::cout << "Problem 0 det id for " << typeid(*hit).name() << ' ' << hit->det()->geographicalId() << std::endl; if (hitcounter != avtm.size())//no propagation needed for first smoothed (==last fitted) hit