diff --git a/CMakeLists.txt b/CMakeLists.txt index 0295f2d..83a9e49 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,8 +11,8 @@ PROJECT( MarlinTrkProcessors ) # project version SET( ${PROJECT_NAME}_VERSION_MAJOR 2 ) -SET( ${PROJECT_NAME}_VERSION_MINOR 15 ) -SET( ${PROJECT_NAME}_VERSION_PATCH 2 ) +SET( ${PROJECT_NAME}_VERSION_MINOR 17 ) +SET( ${PROJECT_NAME}_VERSION_PATCH 0 ) cmake_policy(SET CMP0008 NEW) diff --git a/source/Refitting/include/RefitFinal.h b/source/Refitting/include/RefitFinal.h index f8ebc9a..06ec526 100644 --- a/source/Refitting/include/RefitFinal.h +++ b/source/Refitting/include/RefitFinal.h @@ -15,11 +15,10 @@ class IMarlinTrack; } // namespace MarlinTrk class RefitFinal : public marlin::Processor { - /// Nhits cut struct - struct NHitsCut{ - std::vector detIDs {} ; // list of detectors IDs from which hits will be summed - int nHits_min {}; // minimum number of hits + struct NHitsCut { + std::vector detIDs{}; // list of detectors IDs from which hits will be summed + int nHits_min{}; // minimum number of hits }; public: @@ -80,8 +79,8 @@ class RefitFinal : public marlin::Processor { bool _MSOn = true; bool _ElossOn = true; - StringVec _NHitsCutsStr {}; - std::vector _NHitsCuts {}; + StringVec _NHitsCutsStr{}; + std::vector _NHitsCuts{}; double _ReducedChi2Cut = -1.0; bool _SmoothOn = false; diff --git a/source/Refitting/src/RefitFinal.cc b/source/Refitting/src/RefitFinal.cc index 5db3223..34ccc30 100644 --- a/source/Refitting/src/RefitFinal.cc +++ b/source/Refitting/src/RefitFinal.cc @@ -232,10 +232,10 @@ void RefitFinal::processEvent(LCEvent* evt) { auto lcioTrkPtr = lcio_trk.release(); // if required apply the ReducedChi2 cut - if (_ReducedChi2Cut > 0. && lcioTrkPtr->getChi2()/lcioTrkPtr->getNdf() > _ReducedChi2Cut) { + if (_ReducedChi2Cut > 0. && lcioTrkPtr->getChi2() / lcioTrkPtr->getNdf() > _ReducedChi2Cut) { streamlog_out(DEBUG5) << "Skip track " << lcioTrkPtr->id() << ": " - << "Chi2/ndof " << lcioTrkPtr->getChi2()/lcioTrkPtr->getNdf() << std::endl; - counter++; + << "Chi2/ndof " << lcioTrkPtr->getChi2() / lcioTrkPtr->getNdf() << std::endl; + counter++; continue; } @@ -251,10 +251,9 @@ void RefitFinal::processEvent(LCEvent* evt) { } } // for loop to the tracks - streamlog_out(MESSAGE4) << " Final number of tracks " << trackVec->getNumberOfElements() + streamlog_out(MESSAGE4) << " Final number of tracks " << trackVec->getNumberOfElements() << " Skipped tracks: " << counter << std::endl; - evt->addCollection(trackVec, _output_track_col_name); if (input_rel_col) { evt->addCollection(trackRelationCollection, _output_track_rel_name); @@ -299,4 +298,3 @@ int RefitFinal::FitInit2(Track* track, MarlinTrk::IMarlinTrack* marlinTrk) { return MarlinTrk::IMarlinTrack::success; } - diff --git a/source/Utils/include/FilterClusters.h b/source/Utils/include/FilterClusters.h index 7dd215c..06958f9 100644 --- a/source/Utils/include/FilterClusters.h +++ b/source/Utils/include/FilterClusters.h @@ -1,56 +1,67 @@ #pragma once +#include #include -//#include +// #include -namespace TrackPerf -{ -} +namespace TrackPerf {} -class FilterClusters : public marlin::Processor -{ +class FilterClusters : public marlin::Processor { public: - virtual Processor* newProcessor() { return new FilterClusters ; } + virtual Processor* newProcessor() { return new FilterClusters; } - FilterClusters(const FilterClusters &) = delete ; - FilterClusters& operator =(const FilterClusters &) = delete ; - FilterClusters() ; + FilterClusters(const FilterClusters&) = delete; + FilterClusters& operator=(const FilterClusters&) = delete; + FilterClusters(); - /** Called at the begin of the job before anything is read. - * Use to initialize the processor, e.g. book histograms. - */ - virtual void init() ; + /** Called at the begin of the job before anything is read. + * Use to initialize the processor, e.g. book histograms. + */ + virtual void init(); - /** Called for every run. - */ - virtual void processRunHeader( LCRunHeader* run ) ; + /** Called for every run. + */ + virtual void processRunHeader(LCRunHeader* run); - /** Called for every event - the working horse. - */ - virtual void processEvent(LCEvent* evt) ; + /** Called for every event - the working horse. + */ + virtual void processEvent(LCEvent* evt); - - /** Called after data processing for clean up. - */ - virtual void end() ; + /** Called after data processing for clean up. + */ + virtual void end(); private: - //! Input track collection - std::string _InTrackerHitCollection {}; - std::string _InRelationCollection {}; - - //! Output track collection - std::string _OutTrackerHitCollection {}; - std::string _OutRelationCollection {}; - - //! Ranges for theta - std::vector _ThetaRanges {}; - - //! Cut-offs for cluster size in various theta ranges - std::vector _ClusterSize {}; - - //! Layers to be filtered - std::vector _Layers {}; - -}; + //! Input track collection + std::string _InTrackerHitCollection{}; + std::string _InRelationCollection{}; + std::string _InSimTrackerHitCollection{}; + + //! Output track collection + std::string _OutTrackerHitCollection{}; + std::string _OutRelationCollection{}; + std::string _OutSimTrackerHitCollection{}; + + //! Ranges for theta + std::vector _ThetaRanges{}; + + //! Cut-offs for cluster size in various theta ranges + std::vector _ClusterSize{}; + + //! Layers to be filtered + std::vector _Layers{}; + + //! Number of bins in theta + std::string _ThetaBins{}; + + bool m_fillHistos{}; + + // --- Diagnostic histograms: + TH1F* m_clusterTheta_beforeCut = nullptr; + TH1F* m_clusterTheta_afterCut = nullptr; + TH1F* m_clusterLayer_beforeCut = nullptr; + TH1F* m_clusterLayer_afterCut = nullptr; + TH1F* m_clusterSize_beforeCut = nullptr; + TH1F* m_clusterSize_afterCut = nullptr; +}; \ No newline at end of file diff --git a/source/Utils/include/FilterConeHits.h b/source/Utils/include/FilterConeHits.h index 8b6b36d..c1b31a6 100644 --- a/source/Utils/include/FilterConeHits.h +++ b/source/Utils/include/FilterConeHits.h @@ -1,16 +1,15 @@ #ifndef FilterConeHits_h #define FilterConeHits_h 1 -#include "marlin/Processor.h" #include "lcio.h" +#include "marlin/Processor.h" #include #include #include -using namespace lcio ; -using namespace marlin ; - +using namespace lcio; +using namespace marlin; /** Utility processor that selects and saves the tracker hits that are included in * a DeltaR cone around the MC particle direction along with the corresponding @@ -32,62 +31,53 @@ using namespace marlin ; */ class FilterConeHits : public Processor { - - public: - - virtual Processor* newProcessor() { return new FilterConeHits ; } - - - FilterConeHits() ; +public: + virtual Processor* newProcessor() { return new FilterConeHits; } + + FilterConeHits(); FilterConeHits(const FilterConeHits&) = delete; FilterConeHits& operator=(const FilterConeHits&) = delete; - - virtual void init() ; - - virtual void processRunHeader( LCRunHeader* run ) ; - - virtual void processEvent( LCEvent * evt ) ; - - virtual void check( LCEvent * evt ) ; - - virtual void end() ; - - - protected: - // --- Input/output collection names: - std::string m_inputMCParticlesCollName{} ; - std::vector m_inputTrackerHitsCollNames{} ; - std::vector m_inputTrackerSimHitsCollNames{} ; - std::vector m_inputTrackerHitRelNames{} ; - std::vector m_outputTrackerHitsCollNames{} ; - std::vector m_outputTrackerSimHitsCollNames{} ; - std::vector m_outputTrackerHitRelNames{} ; + virtual void init(); + + virtual void processRunHeader(LCRunHeader* run); + + virtual void processEvent(LCEvent* evt); + + virtual void check(LCEvent* evt); + virtual void end(); + +protected: + // --- Input/output collection names: + std::string m_inputMCParticlesCollName{}; + std::vector m_inputTrackerHitsCollNames{}; + std::vector m_inputTrackerSimHitsCollNames{}; + std::vector m_inputTrackerHitRelNames{}; + std::vector m_outputTrackerHitsCollNames{}; + std::vector m_outputTrackerSimHitsCollNames{}; + std::vector m_outputTrackerHitRelNames{}; // --- Processor parameters: - bool m_fillHistos{} ; - double m_deltaRCut{} ; + bool m_fillHistos{}; + double m_deltaRCut{}; + double m_dist3DCut{}; // --- Diagnostic histograms: - TH1F* m_distXY = nullptr ; - TH1F* m_distZ = nullptr ; - TH1F* m_dist3D = nullptr ; - TH1F* m_angle = nullptr ; - TH1F* m_pathLength = nullptr ; - TH1F* m_time = nullptr ; + TH1F* m_distXY = nullptr; + TH1F* m_distZ = nullptr; + TH1F* m_dist3D = nullptr; + TH1F* m_angle = nullptr; + TH1F* m_pathLength = nullptr; + TH1F* m_time = nullptr; // --- Magneti field value - double m_magneticField{} ; + double m_magneticField{}; const double trackerOuterRadius = 1500.; // mm - - // --- Run and event counters: - int _nRun{} ; - int _nEvt{} ; -} ; + // --- Run and event counters: + int _nRun{}; + int _nEvt{}; +}; #endif - - - diff --git a/source/Utils/include/FilterDoubleLayerHits.h b/source/Utils/include/FilterDoubleLayerHits.h index e14a69c..d396e5e 100644 --- a/source/Utils/include/FilterDoubleLayerHits.h +++ b/source/Utils/include/FilterDoubleLayerHits.h @@ -1,8 +1,8 @@ #ifndef FilterDoubleLayerHits_h #define FilterDoubleLayerHits_h 1 -#include "marlin/Processor.h" #include "marlin/EventModifier.h" +#include "marlin/Processor.h" #include "DDRec/SurfaceManager.h" #include @@ -12,10 +12,8 @@ #include - -using namespace lcio ; -using namespace marlin ; - +using namespace lcio; +using namespace marlin; /** Utility processor that removes tracker hits in double layers if they don't have * a corresponding close-by hit in the other sublayer. @@ -31,12 +29,10 @@ using namespace marlin ; */ class FilterDoubleLayerHits : public Processor { - protected: - static const size_t NHITS_MAX = 10000000; - struct SensorPosition{ + struct SensorPosition { unsigned int layer; unsigned int side; unsigned int ladder; @@ -48,89 +44,79 @@ class FilterDoubleLayerHits : public Processor { }; /// Double layer cut struct - struct DoubleLayerCut{ - unsigned int layer0 ; - unsigned int layer1 ; - double dPhi_max ; - double dTheta_max ; + struct DoubleLayerCut { + unsigned int layer0; + unsigned int layer1; + double dPhi_max; + double dTheta_max; }; +public: + virtual Processor* newProcessor() { return new FilterDoubleLayerHits; } - public: - - virtual Processor* newProcessor() { return new FilterDoubleLayerHits ; } - - - FilterDoubleLayerHits() ; - FilterDoubleLayerHits(const FilterDoubleLayerHits&) = delete ; - FilterDoubleLayerHits& operator=(const FilterDoubleLayerHits&) = delete ; + FilterDoubleLayerHits(); + FilterDoubleLayerHits(const FilterDoubleLayerHits&) = delete; + FilterDoubleLayerHits& operator=(const FilterDoubleLayerHits&) = delete; - virtual const std::string & name() const { return Processor::name() ; } + virtual const std::string& name() const { return Processor::name(); } /** Called at the begin of the job before anything is read. * Use to initialize the processor, e.g. book histograms. */ - virtual void init() ; + virtual void init(); /** Called for every run. */ - virtual void processRunHeader( LCRunHeader* run ) ; + virtual void processRunHeader(LCRunHeader* run); /** Called for every event - the working horse. */ - virtual void processEvent( LCEvent * evt ) ; - - - virtual void check( LCEvent * evt ) ; + virtual void processEvent(LCEvent* evt); + virtual void check(LCEvent* evt); /** Called after data processing for clean up. */ - virtual void end() ; - + virtual void end(); - protected: - - dd4hep::rec::Vector2D globalToLocal(long int cellID, const dd4hep::rec::Vector3D& posGlobal, dd4hep::rec::ISurface** surf) ; +protected: + dd4hep::rec::Vector2D globalToLocal(long int cellID, const dd4hep::rec::Vector3D& posGlobal, + dd4hep::rec::ISurface** surf); ////Input collection name. - std::string _inColName {}; + std::string _inColName{}; ////Output collection name. - std::string _outColName {}; + std::string _outColName{}; ////Maximum time difference between hits in a doublet - double _dtMax {}; + double _dtMax{}; ////Double layer cuts configuration - StringVec _dlCutConfigs {}; + StringVec _dlCutConfigs{}; ////Whether to fill diagnostic histograms - bool _fillHistos {false}; + bool _fillHistos{false}; ////Subdetector name (needed to get the sensor surface manager) - std::string _subDetName {}; + std::string _subDetName{}; - std::vector _dlCuts {}; + std::vector _dlCuts{}; ////Surface map for getting local hit positions at sensor surface - const dd4hep::rec::SurfaceMap* _map {nullptr}; - + const dd4hep::rec::SurfaceMap* _map{nullptr}; ////Array of flags for hits to be accepted - bool _hitAccepted[NHITS_MAX] ; + bool _hitAccepted[NHITS_MAX]; ////Map of vectors of hits grouped by position in the detector - std::map > _hitsGrouped{}; + std::map> _hitsGrouped{}; ////Monitoring histograms - std::map _histos {}; + std::map _histos{}; - int _nRun {}; - int _nEvt {}; -} ; + int _nRun{}; + int _nEvt{}; +}; #endif - - - diff --git a/source/Utils/include/FilterTimeHits.h b/source/Utils/include/FilterTimeHits.h index 3f25c77..32cf74b 100644 --- a/source/Utils/include/FilterTimeHits.h +++ b/source/Utils/include/FilterTimeHits.h @@ -1,8 +1,8 @@ #ifndef FilterTimeHits_h #define FilterTimeHits_h 1 -#include "marlin/Processor.h" #include "lcio.h" +#include "marlin/Processor.h" #include #include @@ -35,47 +35,45 @@ using namespace marlin; * @version $Id: FilterTimeHits.h,v 0.1 2021-06-08 08:58:00 fmeloni Exp $ */ -class FilterTimeHits : public Processor -{ - +class FilterTimeHits : public Processor { public: - virtual Processor *newProcessor() { return new FilterTimeHits; } + virtual Processor* newProcessor() { return new FilterTimeHits; } - FilterTimeHits(); - FilterTimeHits(const FilterTimeHits&) = delete; - FilterTimeHits& operator=(const FilterTimeHits&) = delete; + FilterTimeHits(); + FilterTimeHits(const FilterTimeHits&) = delete; + FilterTimeHits& operator=(const FilterTimeHits&) = delete; - virtual void init(); + virtual void init(); - virtual void processRunHeader(LCRunHeader *run); + virtual void processRunHeader(LCRunHeader* run); - virtual void processEvent(LCEvent *evt); + virtual void processEvent(LCEvent* evt); - virtual void check(LCEvent *evt); + virtual void check(LCEvent* evt); - virtual void end(); + virtual void end(); protected: - // --- Input/output collection names: - std::vector m_inputTrackerHitsCollNames{}; - std::vector m_inputTrackerSimHitsCollNames{}; - std::vector m_inputTrackerHitRelNames{}; - std::vector m_outputTrackerHitsCollNames{}; - std::vector m_outputTrackerSimHitsCollNames{}; - std::vector m_outputTrackerHitRelNames{}; - - // --- Processor parameters: - bool m_fillHistos{}; - double m_beta{}; - double m_time_min{}; - double m_time_max{}; - - // --- Diagnostic histograms: - TH1F *m_corrected_time = nullptr; - - // --- Run and event counters: - int _nRun{}; - int _nEvt{}; + // --- Input/output collection names: + std::vector m_inputTrackerHitsCollNames{}; + std::vector m_inputTrackerSimHitsCollNames{}; + std::vector m_inputTrackerHitRelNames{}; + std::vector m_outputTrackerHitsCollNames{}; + std::vector m_outputTrackerSimHitsCollNames{}; + std::vector m_outputTrackerHitRelNames{}; + + // --- Processor parameters: + bool m_fillHistos{}; + double m_beta{}; + double m_time_min{}; + double m_time_max{}; + + // --- Diagnostic histograms: + TH1F* m_corrected_time = nullptr; + + // --- Run and event counters: + int _nRun{}; + int _nEvt{}; }; #endif diff --git a/source/Utils/include/FilterTracks.h b/source/Utils/include/FilterTracks.h index 92e2950..e1183b5 100644 --- a/source/Utils/include/FilterTracks.h +++ b/source/Utils/include/FilterTracks.h @@ -2,9 +2,7 @@ #include -namespace TrackPerf -{ -} +namespace TrackPerf {} /** Utility processor that removes tracks based on different track properties. * @@ -24,64 +22,61 @@ namespace TrackPerf * @date 7 December 2022 * @version $Id: $ */ -class FilterTracks : public marlin::Processor -{ +class FilterTracks : public marlin::Processor { public: - virtual Processor* newProcessor() { return new FilterTracks ; } + virtual Processor* newProcessor() { return new FilterTracks; } - FilterTracks(const FilterTracks &) = delete ; - FilterTracks& operator =(const FilterTracks &) = delete ; - FilterTracks() ; + FilterTracks(const FilterTracks&) = delete; + FilterTracks& operator=(const FilterTracks&) = delete; + FilterTracks(); - /** Called at the begin of the job before anything is read. - * Use to initialize the processor, e.g. book histograms. - */ - virtual void init() ; + /** Called at the begin of the job before anything is read. + * Use to initialize the processor, e.g. book histograms. + */ + virtual void init(); - /** Called for every run. - */ - virtual void processRunHeader( LCRunHeader* run ) ; + /** Called for every run. + */ + virtual void processRunHeader(LCRunHeader* run); - virtual void buildBfield() ; + virtual void buildBfield(); - /** Called for every event - the working horse. - */ - virtual void processEvent(LCEvent* evt) ; + /** Called for every event - the working horse. + */ + virtual void processEvent(LCEvent* evt); - - /** Called after data processing for clean up. - */ - virtual void end() ; + /** Called after data processing for clean up. + */ + virtual void end(); private: - //! Input track collection - std::string _InputTrackCollection {}; + //! Input track collection + std::string _InputTrackCollection{}; - //! Output track collection - std::string _OutputTrackCollection {}; + //! Output track collection + std::string _OutputTrackCollection{}; - bool _BarrelOnly = false; - bool _HasCaloState = false; + bool _BarrelOnly = false; + bool _HasCaloState = false; - //! Cut off for total number of hits - int _NHitsTotal = 7; - //! Cut off for number of hits in vertex detector (barrel and endcap combined) - int _NHitsVertex = 3; - //! Cut off for number of hits in inner tracker (barrel and endcap combined) - int _NHitsInner = 2; - //! Cut off for number of hits in outer tracker (barrel and endcap combined) - int _NHitsOuter = 1; + //! Cut off for total number of hits + int _NHitsTotal = 7; + //! Cut off for number of hits in vertex detector (barrel and endcap combined) + int _NHitsVertex = 3; + //! Cut off for number of hits in inner tracker (barrel and endcap combined) + int _NHitsInner = 2; + //! Cut off for number of hits in outer tracker (barrel and endcap combined) + int _NHitsOuter = 1; - //! Cut off for maximum number of holes on track - int _MaxHoles = 0; + //! Cut off for maximum number of holes on track + int _MaxHoles = 0; - //! Cut off for momentum (GeV) - float _MinPt = 1.0; //units GeV + //! Cut off for momentum (GeV) + float _MinPt = 1.0; // units GeV - //! Cut off for spatial and temporal chi squared values - float _Chi2Spatial = 0; + //! Cut off for spatial and temporal chi squared values + float _Chi2Spatial = 0; - //! Default magnetic field value (Tesla) - float _Bz = 3.57; + //! Default magnetic field value (Tesla) + float _Bz = 3.57; }; - diff --git a/source/Utils/include/SplitCollectionByLayer.h b/source/Utils/include/SplitCollectionByLayer.h index e2d83e2..7c4710f 100644 --- a/source/Utils/include/SplitCollectionByLayer.h +++ b/source/Utils/include/SplitCollectionByLayer.h @@ -1,16 +1,14 @@ #ifndef SplitCollectionByLayer_h #define SplitCollectionByLayer_h 1 -#include "marlin/Processor.h" #include "marlin/EventModifier.h" +#include "marlin/Processor.h" #include "lcio.h" #include - -using namespace lcio ; -using namespace marlin ; - +using namespace lcio; +using namespace marlin; /** Utility processor that allows to split a collection of Hits into * several collections based on the layer information in the cellID word. @@ -25,21 +23,19 @@ using namespace marlin ; */ class SplitCollectionByLayer : public Processor { - protected: - - ///helper struct - struct OutColInfo{ + /// helper struct + struct OutColInfo { OutColInfo(const OutColInfo&) = default; OutColInfo& operator=(const OutColInfo&) = default; - std::vector layers{} ; + std::vector layers{}; LCCollection* collection{nullptr}; OutColInfo() : layers(std::vector(0)), collection(nullptr) {} }; /// Enum used for hit types - enum HitType{ + enum HitType { SimTrackerHitType = 1, TrackerHitType, TrackerHitPlaneType, @@ -48,58 +44,50 @@ class SplitCollectionByLayer : public Processor { UnkownType }; - public: - - virtual Processor* newProcessor() { return new SplitCollectionByLayer ; } +public: + virtual Processor* newProcessor() { return new SplitCollectionByLayer; } + SplitCollectionByLayer(); + SplitCollectionByLayer(const SplitCollectionByLayer&) = delete; + SplitCollectionByLayer& operator=(const SplitCollectionByLayer&) = delete; - SplitCollectionByLayer() ; - SplitCollectionByLayer(const SplitCollectionByLayer&) = delete ; - SplitCollectionByLayer& operator=(const SplitCollectionByLayer&) = delete ; - - virtual const std::string & name() const { return Processor::name() ; } + virtual const std::string& name() const { return Processor::name(); } /** Called at the begin of the job before anything is read. * Use to initialize the processor, e.g. book histograms. */ - virtual void init() ; + virtual void init(); /** Called for every run. */ - virtual void processRunHeader( LCRunHeader* run ) ; + virtual void processRunHeader(LCRunHeader* run); /** Called for every event - the working horse. */ - virtual void processEvent( LCEvent * evt ) ; - - - virtual void check( LCEvent * evt ) ; + virtual void processEvent(LCEvent* evt); + virtual void check(LCEvent* evt); /** Called after data processing for clean up. */ - virtual void end() ; - - - protected: + virtual void end(); +protected: ////Input collection name. - std::string _colName {}; + std::string _colName{}; /// Output collections and layers: - StringVec _outColAndLayers {}; + StringVec _outColAndLayers{}; /// Whether to add empty collections to the event - bool _addEmptyCollections{false} ; + bool _addEmptyCollections{false}; - std::map _outCols{} ; + std::map _outCols{}; - HitType _type {}; + HitType _type{}; - int _nRun {}; - int _nEvt {}; -} ; + int _nRun{}; + int _nEvt{}; +}; #endif - - diff --git a/source/Utils/include/SplitCollectionByPolarAngle.h b/source/Utils/include/SplitCollectionByPolarAngle.h index 741e5c3..feea833 100644 --- a/source/Utils/include/SplitCollectionByPolarAngle.h +++ b/source/Utils/include/SplitCollectionByPolarAngle.h @@ -1,16 +1,15 @@ #ifndef SplitCollectionByPolarAngle_h #define SplitCollectionByPolarAngle_h 1 -#include "marlin/Processor.h" #include "lcio.h" +#include "marlin/Processor.h" #include #include #include -using namespace lcio ; -using namespace marlin ; - +using namespace lcio; +using namespace marlin; /** Utility processor that selects and saves the tracker hits with a polar angle * within a given range (defined by the lower and upper limit), along with @@ -32,51 +31,43 @@ using namespace marlin ; */ class SplitCollectionByPolarAngle : public Processor { - public: + virtual Processor* newProcessor() { return new SplitCollectionByPolarAngle; } - virtual Processor* newProcessor() { return new SplitCollectionByPolarAngle ; } - - - SplitCollectionByPolarAngle() ; - SplitCollectionByPolarAngle(const SplitCollectionByPolarAngle&) = delete ; - SplitCollectionByPolarAngle& operator=(const SplitCollectionByPolarAngle&) = delete ; + SplitCollectionByPolarAngle(); + SplitCollectionByPolarAngle(const SplitCollectionByPolarAngle&) = delete; + SplitCollectionByPolarAngle& operator=(const SplitCollectionByPolarAngle&) = delete; - virtual void init() ; + virtual void init(); - virtual void processRunHeader( LCRunHeader* run ) ; + virtual void processRunHeader(LCRunHeader* run); - virtual void processEvent( LCEvent * evt ) ; + virtual void processEvent(LCEvent* evt); - virtual void check( LCEvent * evt ) ; - - virtual void end() ; + virtual void check(LCEvent* evt); + virtual void end(); protected: - // --- Input/output collection names: - std::vector m_inputTrackerHitsCollNames{} ; - std::vector m_inputTrackerSimHitsCollNames{} ; - std::vector m_inputTrackerHitRelNames{} ; - std::vector m_outputTrackerHitsCollNames{} ; - std::vector m_outputTrackerSimHitsCollNames{} ; - std::vector m_outputTrackerHitRelNames{} ; - + std::vector m_inputTrackerHitsCollNames{}; + std::vector m_inputTrackerSimHitsCollNames{}; + std::vector m_inputTrackerHitRelNames{}; + std::vector m_outputTrackerHitsCollNames{}; + std::vector m_outputTrackerSimHitsCollNames{}; + std::vector m_outputTrackerHitRelNames{}; // --- Processor parameters: - bool m_fillHistos{} ; - double m_theta_min{} ; - double m_theta_max{} ; + bool m_fillHistos{}; + double m_theta_min{}; + double m_theta_max{}; // --- Diagnostic histograms: - TH1F* m_theta = nullptr ; - + TH1F* m_theta = nullptr; // --- Run and event counters: - int _nRun{} ; - int _nEvt{} ; - -} ; + int _nRun{}; + int _nEvt{}; +}; #endif diff --git a/source/Utils/src/FilterClusters.cc b/source/Utils/src/FilterClusters.cc index 926de76..314a1dc 100644 --- a/source/Utils/src/FilterClusters.cc +++ b/source/Utils/src/FilterClusters.cc @@ -4,13 +4,18 @@ #include +#include #include #include -#include #include -#include -#include +#include #include +#include +#include +#include +#include + +#include #include @@ -18,166 +23,254 @@ #include "marlin/VerbosityLevels.h" +FilterClusters aFilterClusters; -FilterClusters aFilterClusters ; - -FilterClusters::FilterClusters() - : Processor("FilterClusters") -{ +FilterClusters::FilterClusters() : Processor("FilterClusters") { // modify processor description - _description = "FilterClusters processor filters a collection of tracker hits based on cluster size and outputs a filtered collection"; + _description = "FilterClusters processor filters a collection of tracker hits based on cluster size and outputs a " + "filtered collection"; // register steering parameters: name, description, class-variable, default value - registerProcessorParameter("ThetaRanges", - "Divide theta into bins for different cluster size cuts", - _ThetaRanges, - {} - ); - registerProcessorParameter("ClusterSize", - "Maximum cluster size for each theta range", - _ClusterSize, - {} - ); - registerProcessorParameter("Layers", - "Layers to be filtered", - _Layers, - {} - ); - registerInputCollection( LCIO::TRACKERHIT, - "InTrackerHitCollection" , - "Name of the input tracker hit collection", - _InTrackerHitCollection, - _InTrackerHitCollection - ); - - registerInputCollection( LCIO::LCRELATION, - "InRelationCollection" , - "Name of the input relation collection", - _InRelationCollection, - _InRelationCollection - ); - - registerOutputCollection( LCIO::TRACKERHIT, - "OutTrackerHitCollection" , - "Name of output tracker hit collection", - _OutTrackerHitCollection, - std::string("FilteredVBTrackerHits") - ); - - registerOutputCollection( LCIO::LCRELATION, - "OutRelationCollection" , - "Name of output relation collection", - _OutRelationCollection, - std::string("FilteredVBTrackerHitsRelations") - ); + registerProcessorParameter("ThetaRanges", "Divide theta into bins for different cluster size cuts", _ThetaRanges, {}); + registerProcessorParameter("ThetaBins", "Number of bins in theta", _ThetaBins, {}); + registerProcessorParameter("ClusterSize", "Maximum cluster size for each theta range", _ClusterSize, {}); + registerProcessorParameter("Layers", "Layers to be filtered", _Layers, {}); + registerInputCollection(LCIO::SIMTRACKERHIT, "InSimTrackerHitCollection", + "Name of the input sim tracker hit collection", _InSimTrackerHitCollection, + _InSimTrackerHitCollection); + registerInputCollection(LCIO::TRACKERHITPLANE, "InTrackerHitCollection", "Name of the input tracker hit collection", + _InTrackerHitCollection, _InTrackerHitCollection); + + registerInputCollection(LCIO::LCRELATION, "InRelationCollection", "Name of the input relation collection", + _InRelationCollection, _InRelationCollection); + + registerOutputCollection(LCIO::SIMTRACKERHIT, "OutSimTrackerHitCollection", + "Name of output sim tracker hit collection", _OutSimTrackerHitCollection, + _OutSimTrackerHitCollection); + + registerOutputCollection(LCIO::TRACKERHITPLANE, "OutTrackerHitCollection", "Name of output tracker hit collection", + _OutTrackerHitCollection, _OutTrackerHitCollection); + + registerOutputCollection(LCIO::LCRELATION, "OutRelationCollection", "Name of output relation collection", + _OutRelationCollection, _OutRelationCollection); + + registerProcessorParameter("FillHistograms", "Flag to fill the diagnostic histograms", m_fillHistos, false); } -void FilterClusters::init() -{ +void FilterClusters::init() { // Print the initial parameters - printParameters() ; -} + printParameters(); -void FilterClusters::processRunHeader( LCRunHeader* /*run*/) -{ } + marlin::AIDAProcessor::histogramFactory(this); -void FilterClusters::processEvent( LCEvent * evt ) -{ - // Make the output track collection - LCCollectionVec *OutTrackerHitCollection = new LCCollectionVec(LCIO::TRACKERHIT); - OutTrackerHitCollection->setSubset(true); - LCCollectionVec *OutRelationCollection = new LCCollectionVec(LCIO::LCRELATION); - OutRelationCollection->setSubset(true); + m_clusterTheta_beforeCut = new TH1F("m_ClusterTheta_before", "Cluster Theta [radian]", 32, 0., 3.2); + m_clusterTheta_afterCut = new TH1F("m_ClusterTheta_after", "Cluster Theta [radian]", 32, 0., 3.2); + m_clusterLayer_beforeCut = new TH1F("m_ClusterLayer_before", "Cluster Layer Index", 10, 0, 10); + m_clusterLayer_afterCut = new TH1F("m_ClusterLayer_after", "Cluster Layer Index", 10, 0, 10); + m_clusterSize_beforeCut = new TH1F("m_ClusterSize_before", "Cluster hit multiplicity", 20, 0, 20); + m_clusterSize_afterCut = new TH1F("m_ClusterSize_after", "Cluster hit multiplicity", 20, 0, 20); +} +void FilterClusters::processRunHeader(LCRunHeader* /*run*/) {} + +void FilterClusters::processEvent(LCEvent* evt) { // Get input collection - LCCollection* InTrackerHitCollection = evt->getCollection(_InTrackerHitCollection); - LCCollection* InRelationCollection = evt->getCollection(_InRelationCollection); + LCCollection* InTrackerHitCollection = evt->getCollection(_InTrackerHitCollection); + LCCollection* InRelationCollection = evt->getCollection(_InRelationCollection); + LCCollection* InSimTrackerHitCollection = evt->getCollection(_InSimTrackerHitCollection); + + streamlog_out(DEBUG0) << "Starting processing event." << std::endl; - if( InTrackerHitCollection->getTypeName() != lcio::LCIO::TRACKERHITPLANE ) - { throw EVENT::Exception( "Invalid collection type: " + InTrackerHitCollection->getTypeName() ) ; } + if (InTrackerHitCollection->getTypeName() != lcio::LCIO::TRACKERHITPLANE) { + throw EVENT::Exception("Invalid collection type: " + InTrackerHitCollection->getTypeName()); streamlog_out(DEBUG0) << "Wrong collection type for TrackerHitCollection. \n"; - if( InRelationCollection->getTypeName() != lcio::LCIO::LCRELATION ) - { throw EVENT::Exception( "Invalid collection type: " + InRelationCollection->getTypeName() ) ; } + } + if (InRelationCollection->getTypeName() != lcio::LCIO::LCRELATION) { + throw EVENT::Exception("Invalid collection type: " + InRelationCollection->getTypeName()); streamlog_out(DEBUG0) << "Wrong collection type for InRelationCollection. \n"; + } + if (InSimTrackerHitCollection->getTypeName() != lcio::LCIO::SIMTRACKERHIT) { + throw EVENT::Exception("Invalid collection type: " + InSimTrackerHitCollection->getTypeName()); + streamlog_out(DEBUG0) << "Wrong collection type for SimTrackerHitCollection. \n"; + } + + streamlog_out(DEBUG1) << "Number of Elements in Tracker Hits Collection: " + << InTrackerHitCollection->getNumberOfElements() << std::endl; + + // Make the output collections: reco hits, sim hits, reco-sim relationship + std::string encoderString = InTrackerHitCollection->getParameters().getStringVal("CellIDEncoding"); + LCCollectionVec* OutTrackerHitCollection = new LCCollectionVec(InTrackerHitCollection->getTypeName()); + OutTrackerHitCollection->setSubset(true); + OutTrackerHitCollection->parameters().setValue("CellIDEncoding", encoderString); + + // sim hit output collections + LCCollectionVec* OutSimTrackerHitCollection = new LCCollectionVec(InSimTrackerHitCollection->getTypeName()); + OutSimTrackerHitCollection->parameters().setValue("CellIDEncoding", encoderString); + OutSimTrackerHitCollection->setSubset(true); + // reco-sim relation output collections + UTIL::LCRelationNavigator thitNav = UTIL::LCRelationNavigator(LCIO::TRACKERHITPLANE, LCIO::SIMTRACKERHIT); - streamlog_out(DEBUG0) << "Number of Elements in VB Tracker Hits Collection: " << InTrackerHitCollection->getNumberOfElements() <getNumberOfElements(); ++i) //loop through all hits - { - streamlog_out(DEBUG0) << "Loop over hits opened. \n"; - EVENT::TrackerHit *trkhit=static_cast(InTrackerHitCollection->getElementAt(i)); //define trkhit var, pointer to i'th element of tracker hits - EVENT::LCRelation *rel=static_cast(InRelationCollection->getElementAt(i)); - - //Calculating theta - float x = trkhit->getPosition()[0]; - float y = trkhit->getPosition()[1]; - float z = trkhit->getPosition()[2]; - float r = sqrt(pow(x,2)+pow(y,2)); - float incidentTheta = std::atan(r/z); - - if(incidentTheta<0) - incidentTheta += M_PI; - - //Calculating cluster size - const lcio::LCObjectVec &rawHits = trkhit->getRawHits(); - float max = -1000000; - float min = 1000000; - for (size_t j=0; j( rawHits[j] ); - const double *localPos = hitConstituent->getPosition(); - x = localPos[0]; - y = localPos[1]; - if (y < min){ - min = y; - } - else if (y > max){ - max = y; - } + for (int i = 0; i < InTrackerHitCollection->getNumberOfElements(); ++i) // loop through all hits + { + streamlog_out(DEBUG2) << "Loop over hits:" << i << std::endl; + TrackerHitPlane* trkhit = static_cast( + InTrackerHitCollection->getElementAt(i)); // define trkhit var, pointer to i'th element of tracker hits + + if (!trkhit) { + streamlog_out(WARNING) << "Cannot retrieve valid point to cluster. Skipping it" << std::endl; + } + + // Calculating theta + float x = trkhit->getPosition()[0]; + float y = trkhit->getPosition()[1]; + float z = trkhit->getPosition()[2]; + float r = sqrt(pow(x, 2) + pow(y, 2)); + float incidentTheta = std::atan(r / z); + + if (incidentTheta < 0) + incidentTheta += M_PI; + + // Calculating cluster size + const lcio::LCObjectVec& rawHits = trkhit->getRawHits(); + float ymax = -1000000; + float xmax = -1000000; + float ymin = 1000000; + float xmin = 1000000; + streamlog_out(DEBUG1) << "Looping over hits constituents." << std::endl; + for (size_t j = 0; j < rawHits.size(); ++j) { + lcio::SimTrackerHit* hitConstituent = dynamic_cast(rawHits[j]); + const double* localPos = hitConstituent->getPosition(); + float x_local = localPos[0]; + float y_local = localPos[1]; + + if (y_local < ymin) { + ymin = y_local; } - float cluster_size = (max - min)+1; - - //Get hit subdetector/layer - std::string _encoderString = lcio::LCTrackerCellID::encoding_string(); - UTIL::CellIDDecoder decoder(_encoderString); - uint32_t layerID = decoder(trkhit)["layer"]; - bool filter_layer = false; - for (size_t j=0; j<_Layers.size(); ++j){ - if (layerID == std::stof(_Layers[j])) { - filter_layer = true; - } + if (y_local > ymax) { + ymax = y_local; } - streamlog_out(DEBUG0) << "Filter layer: " << filter_layer << std::endl; - - for (size_t j=0; j<_ThetaRanges.size()-1; ++j) { - streamlog_out( DEBUG0 ) << "theta: " << incidentTheta << std::endl; - float min_theta = std::stof(_ThetaRanges[j]); - float max_theta = std::stof(_ThetaRanges[j+1]); - streamlog_out( DEBUG0 ) << "theta range: " << min_theta << ", " << max_theta << std::endl; - - if(incidentTheta > min and incidentTheta <= max and not filter_layer){ - streamlog_out( DEBUG0 ) << "theta in range" << std::endl; - streamlog_out( DEBUG0 ) << "cluster size cut off: " << _ClusterSize[j] << std::endl; - streamlog_out( DEBUG0 ) << "cluster size: " << cluster_size << std::endl; - if(cluster_size < std::stof(_ClusterSize[j])) { - streamlog_out( DEBUG0 ) << "cluster added" << std::endl; - OutTrackerHitCollection->addElement(trkhit); - OutRelationCollection->addElement(rel); } - else { - streamlog_out( DEBUG0 ) << "cluster rejected" << std::endl; + + if (x_local < xmin) { + xmin = x_local; + } + if (x_local > xmax) { + xmax = x_local; + } + } + // float cluster_size_y = (ymax - ymin) + 1; + // float cluster_size_x = (xmax - xmin) + 1; + float cluster_size = rawHits.size(); + + streamlog_out(DEBUG2) << "Cluster size:" << cluster_size << std::endl; + + // Get hit subdetector/layer + std::string _encoderString = lcio::LCTrackerCellID::encoding_string(); + UTIL::CellIDDecoder decoder(_encoderString); + uint32_t systemID = decoder(trkhit)["system"]; + uint32_t layerID = decoder(trkhit)["layer"]; + bool filter_layer = false; + + streamlog_out(DEBUG3) << "Decoded system/layer:" << systemID << "/" << layerID << std::endl; + + size_t rows = _Layers.size(), cols = std::stoi(_ThetaBins); + if ((rows * (cols + 1) != _ThetaRanges.size()) || (rows * cols != _ClusterSize.size())) { + std::cout + << "Either theta cuts or cluster cuts not provided for each layer. Please change the config, exiting now..." + << std::endl; + return; + } + + std::vector> _thetaBins_byLayer; + std::vector> _clusterSizeCuts_byLayer; + + for (size_t k = 0; k < rows; ++k) { + std::vector row; + for (size_t j = 0; j <= cols; ++j) { + row.push_back(std::stof(_ThetaRanges[j + k * (cols + 1)])); + } + _thetaBins_byLayer.push_back(row); + } + + for (size_t k = 0; k < rows; ++k) { + std::vector row; + for (size_t j = 0; j < cols; ++j) { + row.push_back(std::stof(_ClusterSize[j + k * cols])); + } + _clusterSizeCuts_byLayer.push_back(row); + } + + int layerInd = -1; + for (size_t j = 0; j < _Layers.size(); ++j) { + if (layerID == std::stof(_Layers[j])) { + filter_layer = true; + layerInd = j; + break; + } + } + streamlog_out(DEBUG1) << "Filter layer: " << filter_layer << std::endl; + + if (m_fillHistos) { + m_clusterTheta_beforeCut->Fill(incidentTheta); + m_clusterLayer_beforeCut->Fill(layerID); + m_clusterSize_beforeCut->Fill(cluster_size); + } + + bool store_hit = true; + if (filter_layer) { + store_hit = false; + for (size_t j = 0; j < _thetaBins_byLayer[layerInd].size() - 1; ++j) { + streamlog_out(DEBUG0) << "theta: " << incidentTheta << std::endl; + float min_theta = _thetaBins_byLayer[layerInd][j]; + float max_theta = _thetaBins_byLayer[layerInd][j + 1]; + streamlog_out(DEBUG0) << "theta range: " << min_theta << ", " << max_theta << std::endl; + + if (incidentTheta >= min_theta && incidentTheta <= max_theta && filter_layer) { + streamlog_out(DEBUG0) << "theta in range" << std::endl; + streamlog_out(DEBUG0) << "cluster size cut off: " << _clusterSizeCuts_byLayer[layerInd][j] << std::endl; + streamlog_out(DEBUG0) << "cluster size: " << cluster_size << std::endl; + if (cluster_size < _clusterSizeCuts_byLayer[layerInd][j]) { + store_hit = true; + streamlog_out(DEBUG0) << "Adding reco/sim clusters and relation to output collections" << std::endl; + break; } } - else{ - streamlog_out( DEBUG0 ) << "theta out of range or layer filtered" << std::endl; - } } } + if (store_hit) { + if (m_fillHistos) { + m_clusterTheta_afterCut->Fill(incidentTheta); + m_clusterLayer_afterCut->Fill(layerID); + m_clusterSize_afterCut->Fill(cluster_size); + } + + EVENT::LCRelation* rel = static_cast(InRelationCollection->getElementAt(i)); + SimTrackerHit* simhit = dynamic_cast(rel->getTo()); + + OutSimTrackerHitCollection->addElement(simhit); + OutTrackerHitCollection->addElement(trkhit); + thitNav.addRelation(trkhit, simhit); + } else { + streamlog_out(DEBUG0) << "cluster rejected" << std::endl; + } + } + + LCCollection* OutRelationCollection = thitNav.createLCCollection(); + // Save output track collection - evt->addCollection(OutTrackerHitCollection, _OutTrackerHitCollection); - evt->addCollection(OutRelationCollection, _OutRelationCollection); + evt->addCollection(OutTrackerHitCollection, _OutTrackerHitCollection); + evt->addCollection(OutRelationCollection, _OutRelationCollection); + evt->addCollection(OutSimTrackerHitCollection, _OutSimTrackerHitCollection); + + streamlog_out(MESSAGE) << " " << OutTrackerHitCollection->size() + << " reco clusters added to the collection: " << _OutTrackerHitCollection << ", " + << OutSimTrackerHitCollection->size() + << " sim hits added to the collection: " << _OutSimTrackerHitCollection << " and " + << OutRelationCollection->getNumberOfElements() + << " relations added to the collection: " << _OutRelationCollection << std::endl; } -void FilterClusters::end() -{ } +void FilterClusters::end() {} \ No newline at end of file diff --git a/source/Utils/src/FilterConeHits.cc b/source/Utils/src/FilterConeHits.cc index 2265784..6c79c55 100644 --- a/source/Utils/src/FilterConeHits.cc +++ b/source/Utils/src/FilterConeHits.cc @@ -1,405 +1,297 @@ #include "FilterConeHits.h" -#include #include +#include #include #include #include -#include -#include -#include #include +#include +#include +#include #include #include -#include "HelixClass_double.h" - -using namespace lcio ; -using namespace marlin ; +#include +#include "HelixClass_double.h" -FilterConeHits aFilterConeHits ; +using namespace lcio; +using namespace marlin; +FilterConeHits aFilterConeHits; FilterConeHits::FilterConeHits() : Processor("FilterConeHits") { - // --- Processor description: _description = "FilterConeHits selects tracker hits in a cone opened around a MC particle direction"; - // --- Processor parameters: - - registerProcessorParameter("MCParticleCollection", - "Name of the MCParticle collection", - m_inputMCParticlesCollName, - std::string("MCParticle") ); - - registerProcessorParameter("TrackerHitInputCollections", - "Name of the tracker hit input collections", - m_inputTrackerHitsCollNames, - {} ); - - registerProcessorParameter("TrackerSimHitInputCollections", - "Name of the tracker simhit input collections", - m_inputTrackerSimHitsCollNames, - {} ); - - registerProcessorParameter("TrackerHitInputRelations", - "Name of the tracker hit relation collections", - m_inputTrackerHitRelNames, - {} ); - - registerProcessorParameter("TrackerHitOutputCollections", - "Name of the tracker hit output collections", - m_outputTrackerHitsCollNames, - {} ); - - registerProcessorParameter("TrackerSimHitOutputCollections", - "Name of the tracker simhit output collections", - m_outputTrackerSimHitsCollNames, - {} ); - - registerProcessorParameter("TrackerHitOutputRelations", - "Name of the tracker hit relation collections", - m_outputTrackerHitRelNames, - {} ); - - registerProcessorParameter( "DeltaRCut" , - "Maximum angular distance between the hits and the particle direction" , - m_deltaRCut, - double(1.) ); - - registerProcessorParameter( "FillHistograms", - "Flag to fill the diagnostic histograms", - m_fillHistos, - false ); - - -} + registerProcessorParameter("MCParticleCollection", "Name of the MCParticle collection", m_inputMCParticlesCollName, + std::string("MCParticle")); + registerProcessorParameter("TrackerHitInputCollections", "Name of the tracker hit input collections", + m_inputTrackerHitsCollNames, {}); -void FilterConeHits::init() { + registerProcessorParameter("TrackerSimHitInputCollections", "Name of the tracker simhit input collections", + m_inputTrackerSimHitsCollNames, {}); - streamlog_out(DEBUG) << " init called " << std::endl ; - - // --- Print the processor parameters: + registerProcessorParameter("TrackerHitInputRelations", "Name of the tracker hit relation collections", + m_inputTrackerHitRelNames, {}); - printParameters() ; + registerProcessorParameter("TrackerHitOutputCollections", "Name of the tracker hit output collections", + m_outputTrackerHitsCollNames, {}); + registerProcessorParameter("TrackerSimHitOutputCollections", "Name of the tracker simhit output collections", + m_outputTrackerSimHitsCollNames, {}); - // --- Get the value of the magnetic field + registerProcessorParameter("TrackerHitOutputRelations", "Name of the tracker hit relation collections", + m_outputTrackerHitRelNames, {}); - m_magneticField = MarlinUtil::getBzAtOrigin(); + registerProcessorParameter("DeltaRCut", "Maximum angular distance between the hits and the particle direction", + m_deltaRCut, double(-1.)); - - // --- Initialize the run and event counters: + registerProcessorParameter("Dist3DCut", "Maximum distance between the hits and the extrapolated helix", m_dist3DCut, + double(-1.)); - _nRun = 0 ; - _nEvt = 0 ; + registerProcessorParameter("FillHistograms", "Flag to fill the diagnostic histograms", m_fillHistos, false); +} +void FilterConeHits::init() { + streamlog_out(DEBUG) << " init called " << std::endl; - // --- Initialize the AIDAProcessor and book the diagnostic histograms: + // --- Print the processor parameters: - AIDAProcessor::histogramFactory(this); + printParameters(); - m_distXY = new TH1F("m_distXY", "hit-to-helix XY distance;d_{XY} [mm]", 1000, 0., 1000.); - m_distZ = new TH1F("m_distZ", "hit-to-helix Z distance;d_{Z} [mm]", 1000, 0., 1000.); - m_dist3D = new TH1F("m_dist3D", "hit-to-helix 3D distance;d_{3D} [mm]", 1000, 0., 1000.); - m_angle = new TH1F("m_angle", "angle between hit and particle;angle [rad]", 1000, 0., 1.); - m_time = new TH1F("m_time","time at the point of closest approach;T [mm/GeV]", 1000, 0., 2000.); - m_pathLength = new TH1F("m_pathLength","pathlength at the point of closest approach;L [mm]", 1000, 0., 12000.); - -} + // --- Get the value of the magnetic field + m_magneticField = MarlinUtil::getBzAtOrigin(); -void FilterConeHits::processRunHeader( LCRunHeader* ) { + // --- Initialize the run and event counters: - _nRun++ ; + _nRun = 0; + _nEvt = 0; -} + // --- Initialize the AIDAProcessor and book the diagnostic histograms: + AIDAProcessor::histogramFactory(this); + m_distXY = new TH1F("m_distXY", "hit-to-helix XY distance;d_{XY} [mm]", 1000, 0., 1000.); + m_distZ = new TH1F("m_distZ", "hit-to-helix Z distance;d_{Z} [mm]", 1000, 0., 1000.); + m_dist3D = new TH1F("m_dist3D", "hit-to-helix 3D distance;d_{3D} [mm]", 1000, 0., 1000.); + m_angle = new TH1F("m_angle", "angle between hit and particle;angle [rad]", 1000, 0., 1.); + m_time = new TH1F("m_time", "time at the point of closest approach;T [mm/GeV]", 1000, 0., 2000.); + m_pathLength = new TH1F("m_pathLength", "pathlength at the point of closest approach;L [mm]", 1000, 0., 12000.); +} -void FilterConeHits::processEvent( LCEvent * evt ) { +void FilterConeHits::processRunHeader(LCRunHeader*) { _nRun++; } +void FilterConeHits::processEvent(LCEvent* evt) { // --- Check whether the number of input and output collections match - if ( m_inputTrackerHitsCollNames.size() != m_inputTrackerSimHitsCollNames.size() || - m_inputTrackerHitsCollNames.size() != m_inputTrackerHitRelNames.size() ){ - + if (m_inputTrackerHitsCollNames.size() != m_inputTrackerSimHitsCollNames.size() || + m_inputTrackerHitsCollNames.size() != m_inputTrackerHitRelNames.size()) { std::stringstream err_msg; - err_msg << "Mismatch between the reco and sim hits input collections" - << std::endl ; - - throw EVENT::Exception( err_msg.str() ) ; + err_msg << "Mismatch between the reco and sim hits input collections" << std::endl; + throw EVENT::Exception(err_msg.str()); } - if ( m_outputTrackerHitsCollNames.size() != m_outputTrackerSimHitsCollNames.size() || - m_outputTrackerHitsCollNames.size() != m_outputTrackerHitRelNames.size() ){ - + if (m_outputTrackerHitsCollNames.size() != m_outputTrackerSimHitsCollNames.size() || + m_outputTrackerHitsCollNames.size() != m_outputTrackerHitRelNames.size()) { std::stringstream err_msg; - err_msg << "Mismatch between the reco and sim hits output collections" - << std::endl ; - - throw EVENT::Exception( err_msg.str() ) ; + err_msg << "Mismatch between the reco and sim hits output collections" << std::endl; + throw EVENT::Exception(err_msg.str()); } - // --- Get the MC particles collection: LCCollection* m_inputMCParticles = NULL; try { - m_inputMCParticles = evt->getCollection( m_inputMCParticlesCollName ); - } - catch( lcio::DataNotAvailableException& e ) { + m_inputMCParticles = evt->getCollection(m_inputMCParticlesCollName); + } catch (lcio::DataNotAvailableException& e) { streamlog_out(WARNING) << m_inputMCParticlesCollName << " collection not available" << std::endl; return; } - // --- Get the input hit collections and create the corresponding output collections: const unsigned int nTrackerHitCol = m_inputTrackerHitsCollNames.size(); std::vector inputHitColls(nTrackerHitCol); std::vector inputSimHitColls(nTrackerHitCol); - std::vector inputHitRels(nTrackerHitCol); std::vector outputTrackerHitColls(nTrackerHitCol); std::vector outputTrackerSimHitColls(nTrackerHitCol); - std::vector outputTrackerHitRels(nTrackerHitCol); - - for (unsigned int icol=0; icol outputTrackerHitRels(nTrackerHitCol); + for (unsigned int icol = 0; icol < nTrackerHitCol; ++icol) { // get the reco hits try { inputHitColls[icol] = evt->getCollection(m_inputTrackerHitsCollNames[icol]); - } - catch( lcio::DataNotAvailableException& e ) { - streamlog_out(WARNING) << m_inputTrackerHitsCollNames[icol] - << " collection not available" << std::endl; + } catch (lcio::DataNotAvailableException& e) { + streamlog_out(WARNING) << m_inputTrackerHitsCollNames[icol] << " collection not available" << std::endl; continue; } // get the sim hits try { inputSimHitColls[icol] = evt->getCollection(m_inputTrackerSimHitsCollNames[icol]); - } - catch( lcio::DataNotAvailableException& e ) { - streamlog_out(WARNING) << m_inputTrackerSimHitsCollNames[icol] - << " collection not available" << std::endl; - continue; - } - - // get the reco-sim relations - try { - inputHitRels[icol] = evt->getCollection(m_inputTrackerHitRelNames[icol]); - } - catch( lcio::DataNotAvailableException& e ) { - streamlog_out(WARNING) << m_inputTrackerHitRelNames[icol] - << " collection not available" << std::endl; + } catch (lcio::DataNotAvailableException& e) { + streamlog_out(WARNING) << m_inputTrackerSimHitsCollNames[icol] << " collection not available" << std::endl; continue; } // reco hit output collections - std::string encoderString = inputHitColls[icol]->getParameters().getStringVal( "CellIDEncoding" ); - outputTrackerHitColls[icol] = new LCCollectionVec( inputHitColls[icol]->getTypeName() ); - outputTrackerHitColls[icol]->parameters().setValue( "CellIDEncoding", encoderString ); + std::string encoderString = inputHitColls[icol]->getParameters().getStringVal("CellIDEncoding"); + outputTrackerHitColls[icol] = new LCCollectionVec(inputHitColls[icol]->getTypeName()); + outputTrackerHitColls[icol]->parameters().setValue("CellIDEncoding", encoderString); LCFlagImpl lcFlag(inputHitColls[icol]->getFlag()); outputTrackerHitColls[icol]->setFlag(lcFlag.getFlag()); - + outputTrackerHitColls[icol]->setSubset(true); + // sim hit output collections - outputTrackerSimHitColls[icol] = new LCCollectionVec( inputSimHitColls[icol]->getTypeName() ); - outputTrackerSimHitColls[icol]->parameters().setValue( "CellIDEncoding", encoderString ); + outputTrackerSimHitColls[icol] = new LCCollectionVec(inputSimHitColls[icol]->getTypeName()); + outputTrackerSimHitColls[icol]->parameters().setValue("CellIDEncoding", encoderString); LCFlagImpl lcFlag_sim(inputSimHitColls[icol]->getFlag()); outputTrackerSimHitColls[icol]->setFlag(lcFlag_sim.getFlag()); + outputTrackerSimHitColls[icol]->setSubset(true); // reco-sim relation output collections - outputTrackerHitRels[icol] = new LCCollectionVec( inputHitRels[icol]->getTypeName() ); - LCFlagImpl lcFlag_rel(inputHitRels[icol]->getFlag()) ; - outputTrackerHitRels[icol]->setFlag( lcFlag_rel.getFlag() ) ; - + outputTrackerHitRels[icol] = new UTIL::LCRelationNavigator(LCIO::TRACKERHITPLANE, LCIO::SIMTRACKERHIT); } + // Load relations + std::vector> hit2simhits; + for (const std::string& name : m_inputTrackerHitRelNames) { + // Get the collection of tracker hit relations + LCCollection* trackerHitRelationCollection = evt->getCollection(name); + if (trackerHitRelationCollection == nullptr) + continue; + std::shared_ptr hit2simhit = + std::make_shared(trackerHitRelationCollection); + hit2simhits.push_back(hit2simhit); + } // --- Loop over the MC particles: - - std::vector > hits_to_save(nTrackerHitCol); - - for (int ipart=0; ipartgetNumberOfElements(); ++ipart){ - - MCParticle* part = dynamic_cast( m_inputMCParticles->getElementAt(ipart) ); + for (int ipart = 0; ipart < m_inputMCParticles->getNumberOfElements(); ++ipart) { + MCParticle* part = dynamic_cast(m_inputMCParticles->getElementAt(ipart)); // --- Keep only the generator-level particles: - if ( part->getGeneratorStatus() != 1 ) continue; + if (part->getGeneratorStatus() != 1) + continue; - double part_p = sqrt( part->getMomentum()[0]*part->getMomentum()[0] + - part->getMomentum()[1]*part->getMomentum()[1] + - part->getMomentum()[2]*part->getMomentum()[2] ); + double part_p = + sqrt(part->getMomentum()[0] * part->getMomentum()[0] + part->getMomentum()[1] * part->getMomentum()[1] + + part->getMomentum()[2] * part->getMomentum()[2]); HelixClass_double helix; - helix.Initialize_VP( (double*) part->getVertex(), (double*) part->getMomentum(), - (double) part->getCharge(), m_magneticField ); - + helix.Initialize_VP((double*)part->getVertex(), (double*)part->getMomentum(), (double)part->getCharge(), + m_magneticField); // --- Get the intersection point with the barrel outer cylinder. // N.B.: If the particle spirals and doesn't reach the tracker outer cylinder, // getPointOnCircle returns -1e20. double intersectionPoint[3] = {0., 0., 0.}; - double intersectionTime = helix.getPointOnCircle(trackerOuterRadius,(double*) part->getVertex(), intersectionPoint); - + double intersectionTime = helix.getPointOnCircle(trackerOuterRadius, (double*)part->getVertex(), intersectionPoint); // --- Loop over the tracker hits and select hits inside a cone around the particle trajectory: - for (unsigned int icol=0; icolgetNumberOfElements(); ++ihit){ - - TrackerHitPlane* hit = dynamic_cast(hit_col->getElementAt(ihit)); + for (unsigned int icol = 0; icol < inputHitColls.size(); ++icol) { + LCCollection* hit_col = inputHitColls[icol]; + if (!hit_col) + continue; + + for (int ihit = 0; ihit < hit_col->getNumberOfElements(); ++ihit) { + TrackerHitPlane* hit = dynamic_cast(hit_col->getElementAt(ihit)); + + // --- Skip hits that are in the opposite hemisphere w.r.t. the MC particle direction + if ((hit->getPosition()[0] * part->getMomentum()[0] + hit->getPosition()[1] * part->getMomentum()[1] + + hit->getPosition()[2] * part->getMomentum()[2]) < 0.) + continue; + + // --- Get the distance between the hit and the particle trajectory: + double hit_distance[3] = {0., 0., 0.}; + double timeAtPCA = helix.getDistanceToPoint((double*)hit->getPosition(), hit_distance); + + // --- Exclude the opposite side of the helix w.r.t. the production vertex : + if (timeAtPCA < 0.) + continue; + + // --- This is to avoid that the helix of central trajectories reenter the tracker: + if (timeAtPCA > intersectionTime && intersectionTime != -1.e20) + continue; + + double pathLength = part_p * timeAtPCA; + double hit_angle = atan2(hit_distance[2], pathLength); + + if (m_fillHistos) { + m_distXY->Fill(hit_distance[0]); + m_distZ->Fill(hit_distance[1]); + m_dist3D->Fill(hit_distance[2]); + m_angle->Fill(hit_angle); + m_time->Fill(timeAtPCA); + m_pathLength->Fill(pathLength); + } + + bool save = false; + + if (m_deltaRCut > 0.) { + if (hit_angle < m_deltaRCut) + save = true; + } + + if (m_dist3DCut > 0.) { + if (hit_distance[2] < m_dist3DCut) + save = true; + } + + if (save) { + // Find the sim hit + SimTrackerHit* simhit = nullptr; + const LCObjectVec& simHitVector = hit2simhits[icol]->getRelatedToObjects(hit); + if (!simHitVector.empty()) { // Found the sim hit + simhit = dynamic_cast(simHitVector.at(0)); + + outputTrackerSimHitColls[icol]->addElement(simhit); + outputTrackerHitColls[icol]->addElement(hit); + outputTrackerHitRels[icol]->addRelation(hit, simhit); + } + } - // --- Skip hits that are in the opposite hemisphere w.r.t. the MC particle direction - if ( ( hit->getPosition()[0]*part->getMomentum()[0] + - hit->getPosition()[1]*part->getMomentum()[1] + - hit->getPosition()[2]*part->getMomentum()[2] ) < 0. ) continue; - - - // --- Get the distance between the hit and the particle trajectory: - double hit_distance[3] = {0.,0.,0.}; - double timeAtPCA = helix.getDistanceToPoint((double*) hit->getPosition(), hit_distance); - - // --- Exclude the opposite side of the helix w.r.t. the production vertex : - if ( timeAtPCA < 0. ) continue; - - // --- This is to avoid that the helix of central trajectories reenter the tracker: - if ( timeAtPCA > intersectionTime && intersectionTime != -1.e20 ) continue; - - double pathLength = part_p*timeAtPCA; - double hit_angle = atan2(hit_distance[2], pathLength); - - - if ( m_fillHistos ){ - - m_distXY->Fill(hit_distance[0]); - m_distZ->Fill(hit_distance[1]); - m_dist3D->Fill(hit_distance[2]); - m_angle->Fill(hit_angle); - m_time->Fill(timeAtPCA); - m_pathLength->Fill(pathLength); - - } - - if ( hit_angle < m_deltaRCut ) - hits_to_save[icol].insert(ihit); - } // ihit loop - + } // icol loop } // ipart loop + for (unsigned int icol = 0; icol < inputHitColls.size(); ++icol) { + // Save output track collection + evt->addCollection(outputTrackerHitColls[icol], m_outputTrackerHitsCollNames[icol]); + evt->addCollection(outputTrackerSimHitColls[icol], m_outputTrackerSimHitsCollNames[icol]); + evt->addCollection(outputTrackerHitRels[icol]->createLCCollection(), m_outputTrackerHitRelNames[icol]); - // --- Add the filtered hits to the output collections: - - for (unsigned int icol=0; icol(inputHitColls[icol]->getElementAt(ihit)); - TrackerHitPlaneImpl* hit_new = new TrackerHitPlaneImpl(); - - hit_new->setCellID0(hit->getCellID0()); - hit_new->setCellID1(hit->getCellID1()); - hit_new->setType(hit->getType()); - hit_new->setPosition(hit->getPosition()); - hit_new->setU(hit->getU()); - hit_new->setV(hit->getV()); - hit_new->setdU(hit->getdU()); - hit_new->setdV(hit->getdV()); - hit_new->setEDep(hit->getEDep()); - hit_new->setEDepError(hit->getEDepError()); - hit_new->setTime(hit->getTime()); - hit_new->setQuality(hit->getQuality()); - - outputTrackerHitColls[icol]->addElement( hit_new ); - - - LCRelation* rel = dynamic_cast(inputHitRels[icol]->getElementAt(ihit)); - - - SimTrackerHit* simhit = dynamic_cast(rel->getTo()); - SimTrackerHitImpl* simhit_new = new SimTrackerHitImpl(); - - simhit_new->setCellID0(simhit->getCellID0()); - simhit_new->setCellID1(simhit->getCellID1()); - simhit_new->setPosition(simhit->getPosition()); - simhit_new->setEDep(simhit->getEDep()); - simhit_new->setTime(simhit->getTime()); - simhit_new->setMCParticle(simhit->getMCParticle()); - simhit_new->setMomentum(simhit->getMomentum()); - simhit_new->setPathLength(simhit->getPathLength()); - simhit_new->setQuality(simhit->getQuality()); - simhit_new->setOverlay(simhit->isOverlay()); - simhit_new->setProducedBySecondary(simhit->isProducedBySecondary()); - - outputTrackerSimHitColls[icol]->addElement( simhit_new ); - - - LCRelationImpl* rel_new = new LCRelationImpl(); - - rel_new->setFrom(hit_new); - rel_new->setTo(simhit_new); - rel_new->setWeight(rel->getWeight()); - - outputTrackerHitRels[icol]->addElement( rel_new ); - - } // ihit loop - - streamlog_out( MESSAGE ) << " " << hits_to_save[icol].size() << " hits added to the collections: " - << m_outputTrackerHitsCollNames[icol] << ", " - << m_outputTrackerSimHitsCollNames[icol] << ", " - << m_outputTrackerHitRelNames[icol] << std::endl; - - evt->addCollection( outputTrackerHitColls[icol], m_outputTrackerHitsCollNames[icol] ) ; - evt->addCollection( outputTrackerSimHitColls[icol], m_outputTrackerSimHitsCollNames[icol] ) ; - evt->addCollection( outputTrackerHitRels[icol], m_outputTrackerHitRelNames[icol] ) ; - - streamlog_out( DEBUG5 ) << " output collection " << m_outputTrackerHitsCollNames[icol] << " of type " - << outputTrackerHitColls[icol]->getTypeName() << " added to the event \n" - << " output collection " << m_outputTrackerSimHitsCollNames[icol] << " of type " - << outputTrackerSimHitColls[icol]->getTypeName() << " added to the event \n" - << " output collection " << m_outputTrackerHitRelNames[icol] << " of type " - << outputTrackerHitRels[icol]->getTypeName() << " added to the event " - << std::endl ; + streamlog_out(DEBUG5) << " output collection " << m_outputTrackerHitsCollNames[icol] << " of type " + << outputTrackerHitColls[icol]->getTypeName() << " added to the event \n" + << " output collection " << m_outputTrackerSimHitsCollNames[icol] << " of type " + << outputTrackerSimHitColls[icol]->getTypeName() << " added to the event " << std::endl; } // icol loop - streamlog_out(DEBUG) << " processing event: " << evt->getEventNumber() - << " in run: " << evt->getRunNumber() << std::endl ; - - _nEvt ++ ; - -} - - + streamlog_out(DEBUG) << " processing event: " << evt->getEventNumber() << " in run: " << evt->getRunNumber() + << std::endl; -void FilterConeHits::check( LCEvent * ) { + _nEvt++; } +void FilterConeHits::check(LCEvent*) {} - -void FilterConeHits::end(){ - - std::cout << "FilterConeHits::end() " << name() - << " processed " << _nEvt << " events in " << _nRun << " runs " - << std::endl ; - +void FilterConeHits::end() { + std::cout << "FilterConeHits::end() " << name() << " processed " << _nEvt << " events in " << _nRun << " runs " + << std::endl; } diff --git a/source/Utils/src/FilterDoubleLayerHits.cc b/source/Utils/src/FilterDoubleLayerHits.cc index ba430d4..04f82b7 100644 --- a/source/Utils/src/FilterDoubleLayerHits.cc +++ b/source/Utils/src/FilterDoubleLayerHits.cc @@ -1,10 +1,10 @@ #include "FilterDoubleLayerHits.h" -#include #include +#include #include -#include #include +#include #include "marlin/AIDAProcessor.h" @@ -15,182 +15,146 @@ #include #include +using namespace lcio; +using namespace marlin; -using namespace lcio ; -using namespace marlin ; +FilterDoubleLayerHits aFilterDoubleLayerHits; +FilterDoubleLayerHits::FilterDoubleLayerHits() : Processor("FilterDoubleLayerHits"), _dtMax(-1.0), _nRun(0), _nEvt(0) { + // modify processor description + _description = "remove hits in double layers if they don't have a corresponding close-by hit in the other sublayer"; -FilterDoubleLayerHits aFilterDoubleLayerHits ; + // register steering parameters: name, description, class-variable, default value + registerProcessorParameter("InputCollection", "Name of the input collection with hits", _inColName, + std::string("VXDTrackerHitPlanes")); -FilterDoubleLayerHits::FilterDoubleLayerHits() : Processor("FilterDoubleLayerHits") , - _dtMax(-1.0), _nRun(0), _nEvt(0) { + registerProcessorParameter("OutputCollection", "Name of the output collection with filtered hits", _outColName, + std::string("VXDTrackerHitPlanes_DLFiltered")); - // modify processor description - _description = "remove hits in double layers if they don't have a corresponding close-by hit in the other sublayer" ; + registerProcessorParameter("SubDetectorName", "Name of dub detector", _subDetName, std::string("Vertex")); + registerProcessorParameter("FillHistograms", "Whether to fill diagnostic histograms", _fillHistos, false); - // register steering parameters: name, description, class-variable, default value - registerProcessorParameter( "InputCollection" , - "Name of the input collection with hits" , - _inColName , - std::string("VXDTrackerHitPlanes") - ); - - registerProcessorParameter( "OutputCollection" , - "Name of the output collection with filtered hits" , - _outColName , - std::string("VXDTrackerHitPlanes_DLFiltered") - ); - - registerProcessorParameter( "SubDetectorName" , - "Name of dub detector" , - _subDetName , - std::string("Vertex") ); - - registerProcessorParameter( "FillHistograms" , - "Whether to fill diagnostic histograms" , - _fillHistos , - false ); - - registerProcessorParameter( "DeltaTimeMax" , - "Maximum time difference between hits in a doublet [ns]" , - _dtMax , - -1.0 ); - - StringVec dlCutConfigsEx ; - dlCutConfigsEx.push_back("0") ; - dlCutConfigsEx.push_back("1") ; - dlCutConfigsEx.push_back("0.5") ; - dlCutConfigsEx.push_back("0.05") ; - - registerProcessorParameter("DoubleLayerCuts" , - "Layer IDs and angular cuts [mrad] to be applied: " , - _dlCutConfigs , - dlCutConfigsEx - ); + registerProcessorParameter("DeltaTimeMax", "Maximum time difference between hits in a doublet [ns]", _dtMax, -1.0); + StringVec dlCutConfigsEx; + dlCutConfigsEx.push_back("0"); + dlCutConfigsEx.push_back("1"); + dlCutConfigsEx.push_back("0.5"); + dlCutConfigsEx.push_back("0.05"); + registerProcessorParameter("DoubleLayerCuts", + "Layer IDs and angular cuts [mrad] to be applied: ", + _dlCutConfigs, dlCutConfigsEx); } - -dd4hep::rec::Vector2D FilterDoubleLayerHits::globalToLocal(long int cellID, const dd4hep::rec::Vector3D& posGlobal, dd4hep::rec::ISurface** surfptr=nullptr) { +dd4hep::rec::Vector2D FilterDoubleLayerHits::globalToLocal(long int cellID, const dd4hep::rec::Vector3D& posGlobal, + dd4hep::rec::ISurface** surfptr = nullptr) { dd4hep::rec::ISurface* surf; // Using directly the provided surface object if available if (surfptr && *surfptr) { surf = *surfptr; } else { // Finding the surface corresponding to the cellID - dd4hep::rec::SurfaceMap::const_iterator surfIt = _map->find( cellID ); - if( surfIt == _map->end() ){ + dd4hep::rec::SurfaceMap::const_iterator surfIt = _map->find(cellID); + if (surfIt == _map->end()) { std::stringstream err; err << " FilterDoubleLayerHits::processEvent(): no surface found for cellID: " << cellID; - throw Exception ( err.str() ); + throw Exception(err.str()); } surf = surfIt->second; // Saving the surface object outside the function to be reused for the same cellID - if (surfptr) *surfptr = surf; + if (surfptr) + *surfptr = surf; } // Converting global position to local in [cm] - dd4hep::rec::Vector2D posLocal = surf->globalToLocal( dd4hep::mm * posGlobal ); + dd4hep::rec::Vector2D posLocal = surf->globalToLocal(dd4hep::mm * posGlobal); - return dd4hep::rec::Vector2D( posLocal.u() / dd4hep::mm, posLocal.v() / dd4hep::mm ); + return dd4hep::rec::Vector2D(posLocal.u() / dd4hep::mm, posLocal.v() / dd4hep::mm); } - void FilterDoubleLayerHits::init() { - - streamlog_out(DEBUG) << " init called " << std::endl ; - printParameters() ; + streamlog_out(DEBUG) << " init called " << std::endl; + printParameters(); // Extracting double-layer cut configurations - _dlCuts.resize( _dlCutConfigs.size() / 4 ) ; - - unsigned i=0,index=0 ; - while( i < _dlCutConfigs.size() ){ - _dlCuts[index].layer0 = std::atoi( _dlCutConfigs[ i++ ].c_str() ) ; - _dlCuts[index].layer1 = std::atoi( _dlCutConfigs[ i++ ].c_str() ) ; - _dlCuts[index].dPhi_max = std::atof( _dlCutConfigs[ i++ ].c_str() ) / 1e3 ; // converting mrad -> rad - _dlCuts[index].dTheta_max = std::atof( _dlCutConfigs[ i++ ].c_str() ) / 1e3 ; // converting mrad -> rad - ++index ; + _dlCuts.resize(_dlCutConfigs.size() / 4); + + unsigned i = 0, index = 0; + while (i < _dlCutConfigs.size()) { + _dlCuts[index].layer0 = std::atoi(_dlCutConfigs[i++].c_str()); + _dlCuts[index].layer1 = std::atoi(_dlCutConfigs[i++].c_str()); + _dlCuts[index].dPhi_max = std::atof(_dlCutConfigs[i++].c_str()) / 1e3; // converting mrad -> rad + _dlCuts[index].dTheta_max = std::atof(_dlCutConfigs[i++].c_str()) / 1e3; // converting mrad -> rad + ++index; } // Get the surface map from the SurfaceManager dd4hep::Detector& theDetector = dd4hep::Detector::getInstance(); dd4hep::rec::SurfaceManager& surfMan = *theDetector.extension(); - dd4hep::DetElement det = theDetector.detector( _subDetName ); + dd4hep::DetElement det = theDetector.detector(_subDetName); - _map = surfMan.map( det.name() ); + _map = surfMan.map(det.name()); - if( !_map ) { + if (!_map) { std::stringstream err; err << " Could not find surface map for detector: " << _subDetName << " in SurfaceManager"; - throw Exception( err.str() ); + throw Exception(err.str()); } // Booking diagnostic histograms for each configured cut AIDAProcessor::histogramFactory(this); char hname[100]; - for(size_t iCut=0, nCuts=_dlCuts.size() ; iCut layers{cut.layer0, cut.layer1}; for (auto layer : layers) { sprintf(hname, "h2_posUV_rejected_layer_%d", layer); - _histos[ std::string(hname) ] = new TH2F( hname , ";U [mm]; V [mm]", 500, -100, 100, 1000, -200, 200 ); + _histos[std::string(hname)] = new TH2F(hname, ";U [mm]; V [mm]", 500, -100, 100, 1000, -200, 200); } } // Printing the configured cut - streamlog_out( DEBUG5 ) << iCut << ". layers: " << cut.layer0 << " >> " << cut.layer1 << "; dPhi: " - << cut.dPhi_max << " rad; dTheta: " << cut.dTheta_max << " rad" << std::endl; + streamlog_out(DEBUG5) << iCut << ". layers: " << cut.layer0 << " >> " << cut.layer1 << "; dPhi: " << cut.dPhi_max + << " rad; dTheta: " << cut.dTheta_max << " rad" << std::endl; } - _nRun = 0 ; - _nEvt = 0 ; - -} - - -void FilterDoubleLayerHits::processRunHeader( LCRunHeader* ) { - - _nRun++ ; + _nRun = 0; + _nEvt = 0; } +void FilterDoubleLayerHits::processRunHeader(LCRunHeader*) { _nRun++; } - - -void FilterDoubleLayerHits::processEvent( LCEvent * evt ) { - - +void FilterDoubleLayerHits::processEvent(LCEvent* evt) { LCCollection* col(0); try { - col = evt->getCollection( _inColName ); - } - catch( lcio::DataNotAvailableException& e) { - streamlog_out( WARNING ) << " input collection not in event : " << _inColName << " - nothing to do !!! " << std::endl; + col = evt->getCollection(_inColName); + } catch (lcio::DataNotAvailableException& e) { + streamlog_out(WARNING) << " input collection not in event : " << _inColName << " - nothing to do !!! " + << std::endl; return; } - - std::string encoderString = col->getParameters().getStringVal( "CellIDEncoding" ); - UTIL::CellIDDecoder decoder( encoderString ) ; + std::string encoderString = col->getParameters().getStringVal("CellIDEncoding"); + UTIL::CellIDDecoder decoder(encoderString); //---- create the output collection - LCCollectionVec* colOut = new LCCollectionVec( col->getTypeName() ); - colOut->setSubset( true ) ; - colOut->parameters().setValue( "CellIDEncoding", encoderString ); + LCCollectionVec* colOut = new LCCollectionVec(col->getTypeName()); + colOut->setSubset(true); + colOut->parameters().setValue("CellIDEncoding", encoderString); // Set acceptance flags for all hits to FALSE const size_t nHit = col->getNumberOfElements(); @@ -198,12 +162,11 @@ void FilterDoubleLayerHits::processEvent( LCEvent * evt ) { // Splitting hits by sensor ids for faster association _hitsGrouped.clear(); - for (size_t iHit = 0; iHit < nHit ; iHit++) { - - TrackerHitPlane* h = (TrackerHitPlane*)col->getElementAt( iHit ); + for (size_t iHit = 0; iHit < nHit; iHit++) { + TrackerHitPlane* h = (TrackerHitPlane*)col->getElementAt(iHit); - unsigned int layerID = decoder(h)["layer"]; - unsigned int sideID = decoder(h)["side"]; + unsigned int layerID = decoder(h)["layer"]; + unsigned int sideID = decoder(h)["side"]; unsigned int ladderID = decoder(h)["module"]; unsigned int moduleID = decoder(h)["sensor"]; @@ -216,26 +179,28 @@ void FilterDoubleLayerHits::processEvent( LCEvent * evt ) { } //---- loop over hits - for (size_t iHit = 0; iHit < nHit ; iHit++) { - + for (size_t iHit = 0; iHit < nHit; iHit++) { // Skipping hits that are already accepted - if (_hitAccepted[iHit]) continue; + if (_hitAccepted[iHit]) + continue; - TrackerHitPlane* h = (TrackerHitPlane*)col->getElementAt( iHit ); + TrackerHitPlane* h = (TrackerHitPlane*)col->getElementAt(iHit); - unsigned int layerID = decoder(h)["layer"]; - unsigned int sideID = decoder(h)["side"]; + unsigned int layerID = decoder(h)["layer"]; + unsigned int sideID = decoder(h)["side"]; unsigned int ladderID = decoder(h)["module"]; unsigned int moduleID = decoder(h)["sensor"]; - streamlog_out( DEBUG5 ) << " Checking 1st hit " << iHit << " / " << nHit << " at layer: " << layerID << " ladder: " << ladderID << " module: " << moduleID << std::endl ; + streamlog_out(DEBUG5) << " Checking 1st hit " << iHit << " / " << nHit << " at layer: " << layerID + << " ladder: " << ladderID << " module: " << moduleID << std::endl; const SensorPosition sensPos = {layerID, sideID, ladderID, moduleID}; // Checking if the hit is at the inner double layer to be filtered const DoubleLayerCut* dlCut(0); - for (int iCut=0, nCuts=_dlCuts.size(); iCutlayer1) continue; + if (layerID == dlCut->layer1) + continue; // Getting local and global hit positions - dd4hep::rec::Vector3D posGlobal( h->getPosition()[0], h->getPosition()[1], h->getPosition()[2] ); - dd4hep::rec::Vector2D posLocal = globalToLocal( decoder( h ).getValue(), posGlobal ); + dd4hep::rec::Vector3D posGlobal(h->getPosition()[0], h->getPosition()[1], h->getPosition()[2]); + dd4hep::rec::Vector2D posLocal = globalToLocal(decoder(h).getValue(), posGlobal); // Setting the values for closest hits double dR_min(999.0); @@ -264,28 +230,32 @@ void FilterDoubleLayerHits::processEvent( LCEvent * evt ) { size_t nCompatibleHits(0); SensorPosition sensPos2 = sensPos; sensPos2.layer = dlCut->layer1; - dd4hep::rec::ISurface* surf=nullptr; + dd4hep::rec::ISurface* surf = nullptr; // Checking if there are any hits in the corresponding sensor at the other sublayer - if (_hitsGrouped.find(sensPos2) == _hitsGrouped.end()) continue; + if (_hitsGrouped.find(sensPos2) == _hitsGrouped.end()) + continue; for (size_t iHit2 : _hitsGrouped.at(sensPos2)) { - TrackerHitPlane* h2 = (TrackerHitPlane*)col->getElementAt( iHit2 ); + TrackerHitPlane* h2 = (TrackerHitPlane*)col->getElementAt(iHit2); unsigned int layerID2 = decoder(h2)["layer"]; // Checking whether hit is in the time acceptance window double dt = h2->getTime() - h->getTime(); - if (_dtMax >= 0.0 && std::fabs(dt) > _dtMax) continue; + if (_dtMax >= 0.0 && std::fabs(dt) > _dtMax) + continue; // Getting the local and global hit positions - dd4hep::rec::Vector3D posGlobal2( h2->getPosition()[0], h2->getPosition()[1], h2->getPosition()[2] ); - dd4hep::rec::Vector2D posLocal2 = globalToLocal( decoder( h2 ).getValue(), posGlobal2, &surf ); + dd4hep::rec::Vector3D posGlobal2(h2->getPosition()[0], h2->getPosition()[1], h2->getPosition()[2]); + dd4hep::rec::Vector2D posLocal2 = globalToLocal(decoder(h2).getValue(), posGlobal2, &surf); // Checking whether hit is close enough to the 1st one double dU = posLocal2.u() - posLocal.u(); double dTheta = posGlobal2.theta() - posGlobal.theta(); double dPhi = std::fabs(posGlobal2.phi() - posGlobal.phi()); - if (dPhi > dd4hep::pi) dPhi = dd4hep::twopi - dPhi; - double dR = sqrt(dPhi*dPhi + dTheta*dTheta); - streamlog_out( DEBUG0 ) << " Checking 2nd hit at layer: " << layerID2 << "; dPhi: " << dPhi << "; dTheta: " << dTheta << std::endl; + if (dPhi > dd4hep::pi) + dPhi = dd4hep::twopi - dPhi; + double dR = sqrt(dPhi * dPhi + dTheta * dTheta); + streamlog_out(DEBUG0) << " Checking 2nd hit at layer: " << layerID2 << "; dPhi: " << dPhi + << "; dTheta: " << dTheta << std::endl; // Updating the minimal values if (dR < dR_min) { @@ -297,32 +267,35 @@ void FilterDoubleLayerHits::processEvent( LCEvent * evt ) { } // Skipping if the hit is outside the cut window - if (std::fabs(dPhi) > dlCut->dPhi_max) continue; - if (std::fabs(dTheta) > dlCut->dTheta_max) continue; + if (std::fabs(dPhi) > dlCut->dPhi_max) + continue; + if (std::fabs(dTheta) > dlCut->dTheta_max) + continue; nCompatibleHits++; _hitAccepted[iHit2] = true; - streamlog_out( DEBUG5 ) << " Accepted 2nd hit at layer: " << layerID2 << "; dPhi: " << dPhi << "; dTheta: " << dTheta << std::endl; + streamlog_out(DEBUG5) << " Accepted 2nd hit at layer: " << layerID2 << "; dPhi: " << dPhi + << "; dTheta: " << dTheta << std::endl; } // Filling diagnostic histograms if (_fillHistos && dR_min < 998) { - char hname[100]; - sprintf(hname, "h_dU_layers_%d_%d", dlCut->layer0, dlCut->layer1); - _histos[hname]->Fill(dU_closest); - sprintf(hname, "h2_dU_dPhi_layers_%d_%d", dlCut->layer0, dlCut->layer1); - _histos[hname]->Fill(dU_closest, dPhi_closest*1e3); - sprintf(hname, "h_dTheta_layers_%d_%d", dlCut->layer0, dlCut->layer1); - _histos[hname]->Fill(dTheta_closest*1e3); - sprintf(hname, "h_dPhi_layers_%d_%d", dlCut->layer0, dlCut->layer1); - _histos[hname]->Fill(dPhi_closest*1e3); - sprintf(hname, "h_dt_layers_%d_%d", dlCut->layer0, dlCut->layer1); - _histos[hname]->Fill(dt_closest); - } + char hname[100]; + sprintf(hname, "h_dU_layers_%d_%d", dlCut->layer0, dlCut->layer1); + _histos[hname]->Fill(dU_closest); + sprintf(hname, "h2_dU_dPhi_layers_%d_%d", dlCut->layer0, dlCut->layer1); + _histos[hname]->Fill(dU_closest, dPhi_closest * 1e3); + sprintf(hname, "h_dTheta_layers_%d_%d", dlCut->layer0, dlCut->layer1); + _histos[hname]->Fill(dTheta_closest * 1e3); + sprintf(hname, "h_dPhi_layers_%d_%d", dlCut->layer0, dlCut->layer1); + _histos[hname]->Fill(dPhi_closest * 1e3); + sprintf(hname, "h_dt_layers_%d_%d", dlCut->layer0, dlCut->layer1); + _histos[hname]->Fill(dt_closest); + } // Accepting the first hit if it has at least one compatible pair if (nCompatibleHits > 0) { _hitAccepted[iHit] = true; - streamlog_out( DEBUG5 ) << " Accepted 1st hit at layer: " << layerID << std::endl; + streamlog_out(DEBUG5) << " Accepted 1st hit at layer: " << layerID << std::endl; } } @@ -332,12 +305,12 @@ void FilterDoubleLayerHits::processEvent( LCEvent * evt ) { if (!_hitAccepted[iHit]) { // Filling the positions of rejected hits if (_fillHistos) { - TrackerHitPlane* h = (TrackerHitPlane*)col->getElementAt( iHit ); + TrackerHitPlane* h = (TrackerHitPlane*)col->getElementAt(iHit); unsigned int layerID = decoder(h)["layer"]; // Getting local hit position - dd4hep::rec::Vector3D posGlobal( h->getPosition()[0], h->getPosition()[1], h->getPosition()[2] ); - dd4hep::rec::Vector2D posLocal = globalToLocal( decoder( h ).getValue(), posGlobal ); + dd4hep::rec::Vector3D posGlobal(h->getPosition()[0], h->getPosition()[1], h->getPosition()[2]); + dd4hep::rec::Vector2D posLocal = globalToLocal(decoder(h).getValue(), posGlobal); char hname[100]; sprintf(hname, "h2_posUV_rejected_layer_%d", layerID); @@ -345,31 +318,23 @@ void FilterDoubleLayerHits::processEvent( LCEvent * evt ) { } continue; } - colOut->addElement( col->getElementAt( iHit ) ); + colOut->addElement(col->getElementAt(iHit)); nHitsAccepted++; } - streamlog_out( MESSAGE ) << " " << nHitsAccepted << " hits added to collection: " << _outColName << std::endl; - + streamlog_out(MESSAGE) << " " << nHitsAccepted << " hits added to collection: " << _outColName << std::endl; - evt->addCollection( colOut , _outColName ) ; - streamlog_out( DEBUG5 ) << " output collection " << _outColName << " of type " << col->getTypeName() << " added to the event " << std::endl ; + evt->addCollection(colOut, _outColName); + streamlog_out(DEBUG5) << " output collection " << _outColName << " of type " << col->getTypeName() + << " added to the event " << std::endl; - _nEvt ++ ; + _nEvt++; } - - -void FilterDoubleLayerHits::check( LCEvent* ) { +void FilterDoubleLayerHits::check(LCEvent*) { // nothing to check here - could be used to fill checkplots in reconstruction processor } - -void FilterDoubleLayerHits::end(){ - - streamlog_out( MESSAGE ) << "FilterDoubleLayerHits::end() " << name() - << " processed " << _nEvt << " events in " << _nRun << " runs " - << std::endl ; - +void FilterDoubleLayerHits::end() { + streamlog_out(MESSAGE) << "FilterDoubleLayerHits::end() " << name() << " processed " << _nEvt << " events in " + << _nRun << " runs " << std::endl; } - - diff --git a/source/Utils/src/FilterTimeHits.cc b/source/Utils/src/FilterTimeHits.cc index 8e34ca4..1ba2fd1 100644 --- a/source/Utils/src/FilterTimeHits.cc +++ b/source/Utils/src/FilterTimeHits.cc @@ -1,16 +1,18 @@ #include "FilterTimeHits.h" -#include #include +#include #include #include -#include -#include -#include #include +#include +#include +#include + +#include -#include "DD4hep/Detector.h" #include "DD4hep/DD4hepUnits.h" +#include "DD4hep/Detector.h" #include @@ -19,321 +21,207 @@ using namespace marlin; FilterTimeHits aFilterTimeHits; -FilterTimeHits::FilterTimeHits() : Processor("FilterTimeHits") -{ - // --- Processor description: - - _description = "FilterTimeHits selects tracker hits based on their time corrected for the time of flight"; - - // --- Processor parameters: - - registerProcessorParameter("TrackerHitInputCollections", - "Name of the tracker hit input collections", - m_inputTrackerHitsCollNames, - {}); - - registerProcessorParameter("TrackerSimHitInputCollections", - "Name of the tracker simhit input collections", - m_inputTrackerSimHitsCollNames, - {}); - - registerProcessorParameter("TrackerHitInputRelations", - "Name of the tracker hit relation collections", - m_inputTrackerHitRelNames, - {}); - - registerProcessorParameter("TrackerHitOutputCollections", - "Name of the tracker hit output collections", - m_outputTrackerHitsCollNames, - {}); - - registerProcessorParameter("TrackerSimHitOutputCollections", - "Name of the tracker simhit output collections", - m_outputTrackerSimHitsCollNames, - {}); - - registerProcessorParameter("TrackerHitOutputRelations", - "Name of the tracker hit relation collections", - m_outputTrackerHitRelNames, - {}); - - registerProcessorParameter("TargetBeta", - "Target beta=v/c for hit time of flight correction", - m_beta, - double(1.0)); - - registerProcessorParameter("TimeLowerLimit", - "Lower limit on the corrected hit time in ns", - m_time_min, - double(-90.0)); - - registerProcessorParameter("TimeUpperLimit", - "Upper limit on the corrected hit time in ns", - m_time_max, - double(90.0)); - - registerProcessorParameter("FillHistograms", - "Flag to fill the diagnostic histograms", - m_fillHistos, - false); -} - -void FilterTimeHits::init() -{ - - streamlog_out(DEBUG) << " init called " << std::endl; - - // --- Print the processor parameters: - - printParameters(); - - // --- Initialize the run and event counters: - - _nRun = 0; - _nEvt = 0; - - // --- Initialize the AIDAProcessor and book the diagnostic histograms: - - AIDAProcessor::histogramFactory(this); - - m_corrected_time = new TH1F("m_corrected_time", "Corrected time of the hit [ps]", 1000, -250., 250.); -} - -void FilterTimeHits::processRunHeader(LCRunHeader *) -{ - _nRun++; -} - -void FilterTimeHits::processEvent(LCEvent *evt) -{ - - streamlog_out(DEBUG) << " processing event: " << evt->getEventNumber() - << " in run: " << evt->getRunNumber() << std::endl; - - // --- Check whether the number of input and output collections match - - if (m_inputTrackerHitsCollNames.size() != m_inputTrackerSimHitsCollNames.size() || - m_inputTrackerHitsCollNames.size() != m_inputTrackerHitRelNames.size()) - { - - std::stringstream err_msg; - err_msg << "Mismatch between the reco and sim hits input collections" - << std::endl; +FilterTimeHits::FilterTimeHits() : Processor("FilterTimeHits") { + // --- Processor description: - throw EVENT::Exception(err_msg.str()); - } + _description = "FilterTimeHits selects tracker hits based on their time corrected for the time of flight"; - if (m_outputTrackerHitsCollNames.size() != m_outputTrackerSimHitsCollNames.size() || - m_outputTrackerHitsCollNames.size() != m_outputTrackerHitRelNames.size()) - { + // --- Processor parameters: - std::stringstream err_msg; - err_msg << "Mismatch between the reco and sim hits output collections" - << std::endl; + registerProcessorParameter("TrackerHitInputCollections", "Name of the tracker hit input collections", + m_inputTrackerHitsCollNames, {}); - throw EVENT::Exception(err_msg.str()); - } + registerProcessorParameter("TrackerSimHitInputCollections", "Name of the tracker simhit input collections", + m_inputTrackerSimHitsCollNames, {}); - streamlog_out(DEBUG) << " passed container size checks" << std::endl; + registerProcessorParameter("TrackerHitInputRelations", "Name of the tracker hit relation collections", + m_inputTrackerHitRelNames, {}); - // --- Get the input hit collections and create the corresponding output collections: + registerProcessorParameter("TrackerHitOutputCollections", "Name of the tracker hit output collections", + m_outputTrackerHitsCollNames, {}); - const unsigned int nTrackerHitCol = m_inputTrackerHitsCollNames.size(); - std::vector inputHitColls(nTrackerHitCol); - std::vector inputSimHitColls(nTrackerHitCol); - std::vector inputHitRels(nTrackerHitCol); + registerProcessorParameter("TrackerSimHitOutputCollections", "Name of the tracker simhit output collections", + m_outputTrackerSimHitsCollNames, {}); - std::vector outputTrackerHitColls(nTrackerHitCol); - std::vector outputTrackerSimHitColls(nTrackerHitCol); - std::vector outputTrackerHitRels(nTrackerHitCol); + registerProcessorParameter("TrackerHitOutputRelations", "Name of the tracker hit relation collections", + m_outputTrackerHitRelNames, {}); - for (unsigned int icol = 0; icol < nTrackerHitCol; ++icol) - { + registerProcessorParameter("TargetBeta", "Target beta=v/c for hit time of flight correction", m_beta, double(1.0)); - // get the reco hits - try - { - inputHitColls[icol] = evt->getCollection(m_inputTrackerHitsCollNames[icol]); - } - catch (lcio::DataNotAvailableException &e) - { - streamlog_out(WARNING) << m_inputTrackerHitsCollNames[icol] - << " collection not available" << std::endl; - continue; - } + registerProcessorParameter("TimeLowerLimit", "Lower limit on the corrected hit time in ns", m_time_min, + double(-90.0)); - // get the sim hits - try - { - inputSimHitColls[icol] = evt->getCollection(m_inputTrackerSimHitsCollNames[icol]); - } - catch (lcio::DataNotAvailableException &e) - { - streamlog_out(WARNING) << m_inputTrackerSimHitsCollNames[icol] - << " collection not available" << std::endl; - continue; - } + registerProcessorParameter("TimeUpperLimit", "Upper limit on the corrected hit time in ns", m_time_max, double(90.0)); - // get the reco-sim relations - try - { - inputHitRels[icol] = evt->getCollection(m_inputTrackerHitRelNames[icol]); - } - catch (lcio::DataNotAvailableException &e) - { - streamlog_out(WARNING) << m_inputTrackerHitRelNames[icol] - << " collection not available" << std::endl; - continue; - } - - // reco hit output collections - std::string encoderString = inputHitColls[icol]->getParameters().getStringVal("CellIDEncoding"); - outputTrackerHitColls[icol] = new LCCollectionVec(inputHitColls[icol]->getTypeName()); - outputTrackerHitColls[icol]->parameters().setValue("CellIDEncoding", encoderString); - LCFlagImpl lcFlag(inputHitColls[icol]->getFlag()); - outputTrackerHitColls[icol]->setFlag(lcFlag.getFlag()); - - // sim hit output collections - outputTrackerSimHitColls[icol] = new LCCollectionVec(inputSimHitColls[icol]->getTypeName()); - outputTrackerSimHitColls[icol]->parameters().setValue("CellIDEncoding", encoderString); - LCFlagImpl lcFlag_sim(inputSimHitColls[icol]->getFlag()); - outputTrackerSimHitColls[icol]->setFlag(lcFlag_sim.getFlag()); - - // reco-sim relation output collections - outputTrackerHitRels[icol] = new LCCollectionVec(inputHitRels[icol]->getTypeName()); - LCFlagImpl lcFlag_rel(inputHitRels[icol]->getFlag()); - outputTrackerHitRels[icol]->setFlag(lcFlag_rel.getFlag()); - } + registerProcessorParameter("FillHistograms", "Flag to fill the diagnostic histograms", m_fillHistos, false); +} - // --- Loop over the tracker hits and select hits inside the chosen time window: - std::vector> hits_to_save(nTrackerHitCol); +void FilterTimeHits::init() { + streamlog_out(DEBUG) << " init called " << std::endl; - for (unsigned int icol = 0; icol < inputHitColls.size(); ++icol) - { + // --- Print the processor parameters: - LCCollection *hit_col = inputHitColls[icol]; - if (!hit_col) - continue; + printParameters(); - for (int ihit = 0; ihit < hit_col->getNumberOfElements(); ++ihit) - { + // --- Initialize the run and event counters: - if (TrackerHitPlane *hit = dynamic_cast(hit_col->getElementAt(ihit))){ - // Skipping the hit if its time is outside the acceptance time window - double hitT = hit->getTime(); + _nRun = 0; + _nEvt = 0; - dd4hep::rec::Vector3D pos = hit->getPosition(); - double hitR = pos.r(); + // --- Initialize the AIDAProcessor and book the diagnostic histograms: - // Correcting for the propagation time - double dt = hitR / (TMath::C() * m_beta / 1e6); - hitT -= dt; - streamlog_out(DEBUG3) << "corrected hit at R: " << hitR << " mm by propagation time: " << dt << " ns to T: " << hitT << " ns" << std::endl; + AIDAProcessor::histogramFactory(this); - //Apply time window selection - if (hitT < m_time_min || hitT > m_time_max) - { - streamlog_out(DEBUG4) << "hit at T: " << hitT << " ns is rejected by timing cuts" << std::endl; - continue; - } + m_corrected_time = new TH1F("m_corrected_time", "Corrected time of the hit [ps]", 1000, -250., 250.); +} - hits_to_save[icol].insert(ihit); +void FilterTimeHits::processRunHeader(LCRunHeader*) { _nRun++; } - if (m_fillHistos) - m_corrected_time->Fill(hitT); - } +void FilterTimeHits::processEvent(LCEvent* evt) { + streamlog_out(DEBUG) << " processing event: " << evt->getEventNumber() << " in run: " << evt->getRunNumber() + << std::endl; - } // ihit loop + // --- Check whether the number of input and output collections match - } // icol loop + if (m_inputTrackerHitsCollNames.size() != m_inputTrackerSimHitsCollNames.size() || + m_inputTrackerHitsCollNames.size() != m_inputTrackerHitRelNames.size()) { + std::stringstream err_msg; + err_msg << "Mismatch between the reco and sim hits input collections" << std::endl; - // --- Add the filtered hits to the output collections: + throw EVENT::Exception(err_msg.str()); + } - for (unsigned int icol = 0; icol < inputHitColls.size(); ++icol) - { + if (m_outputTrackerHitsCollNames.size() != m_outputTrackerSimHitsCollNames.size() || + m_outputTrackerHitsCollNames.size() != m_outputTrackerHitRelNames.size()) { + std::stringstream err_msg; + err_msg << "Mismatch between the reco and sim hits output collections" << std::endl; - for (auto &ihit : hits_to_save[icol]) - { + throw EVENT::Exception(err_msg.str()); + } - TrackerHitPlane *hit = dynamic_cast(inputHitColls[icol]->getElementAt(ihit)); - TrackerHitPlaneImpl *hit_new = new TrackerHitPlaneImpl(); + streamlog_out(DEBUG) << " passed container size checks" << std::endl; - hit_new->setCellID0(hit->getCellID0()); - hit_new->setCellID1(hit->getCellID1()); - hit_new->setType(hit->getType()); - hit_new->setPosition(hit->getPosition()); - hit_new->setU(hit->getU()); - hit_new->setV(hit->getV()); - hit_new->setdU(hit->getdU()); - hit_new->setdV(hit->getdV()); - hit_new->setEDep(hit->getEDep()); - hit_new->setEDepError(hit->getEDepError()); - hit_new->setTime(hit->getTime()); - hit_new->setQuality(hit->getQuality()); + // --- Get the input hit collections and create the corresponding output collections: - outputTrackerHitColls[icol]->addElement(hit_new); + const unsigned int nTrackerHitCol = m_inputTrackerHitsCollNames.size(); + std::vector inputHitColls(nTrackerHitCol); + std::vector inputSimHitColls(nTrackerHitCol); + std::vector inputHitRels(nTrackerHitCol); - LCRelation *rel = dynamic_cast(inputHitRels[icol]->getElementAt(ihit)); + LCCollectionVec* outputTrackerHitCol = 0; + LCCollectionVec* outputTrackerSimHitCol = 0; + LCCollection* outputTrackerHitRel = 0; - SimTrackerHit *simhit = dynamic_cast(rel->getTo()); - SimTrackerHitImpl *simhit_new = new SimTrackerHitImpl(); + for (unsigned int icol = 0; icol < nTrackerHitCol; ++icol) { + // get the reco hits + try { + inputHitColls[icol] = evt->getCollection(m_inputTrackerHitsCollNames[icol]); + } catch (lcio::DataNotAvailableException& e) { + streamlog_out(WARNING) << m_inputTrackerHitsCollNames[icol] << " collection not available" << std::endl; + continue; + } - simhit_new->setCellID0(simhit->getCellID0()); - simhit_new->setCellID1(simhit->getCellID1()); - simhit_new->setPosition(simhit->getPosition()); - simhit_new->setEDep(simhit->getEDep()); - simhit_new->setTime(simhit->getTime()); - simhit_new->setMCParticle(simhit->getMCParticle()); - simhit_new->setMomentum(simhit->getMomentum()); - simhit_new->setPathLength(simhit->getPathLength()); - simhit_new->setQuality(simhit->getQuality()); - simhit_new->setOverlay(simhit->isOverlay()); - simhit_new->setProducedBySecondary(simhit->isProducedBySecondary()); + // get the sim hits + try { + inputSimHitColls[icol] = evt->getCollection(m_inputTrackerSimHitsCollNames[icol]); + } catch (lcio::DataNotAvailableException& e) { + streamlog_out(WARNING) << m_inputTrackerSimHitsCollNames[icol] << " collection not available" << std::endl; + continue; + } - outputTrackerSimHitColls[icol]->addElement(simhit_new); + // get the reco-sim relations + try { + inputHitRels[icol] = evt->getCollection(m_inputTrackerHitRelNames[icol]); + } catch (lcio::DataNotAvailableException& e) { + streamlog_out(WARNING) << m_inputTrackerHitRelNames[icol] << " collection not available" << std::endl; + continue; + } + } + + // --- Loop over the tracker hits and select hits inside the chosen time window: + // std::vector> hits_to_save(nTrackerHitCol); + + for (unsigned int icol = 0; icol < inputHitColls.size(); ++icol) { + LCCollection* hit_col = inputHitColls[icol]; + if (!hit_col) + continue; + + // reco hit output collections + std::string encoderString = inputHitColls[icol]->getParameters().getStringVal("CellIDEncoding"); + outputTrackerHitCol = new LCCollectionVec(inputHitColls[icol]->getTypeName()); + outputTrackerHitCol->setSubset(true); + outputTrackerHitCol->parameters().setValue("CellIDEncoding", encoderString); + + // sim hit output collections + outputTrackerSimHitCol = new LCCollectionVec(inputSimHitColls[icol]->getTypeName()); + outputTrackerSimHitCol->parameters().setValue("CellIDEncoding", encoderString); + outputTrackerSimHitCol->setSubset(true); + + // reco-sim relation output collections + UTIL::LCRelationNavigator thitNav = UTIL::LCRelationNavigator(LCIO::TRACKERHITPLANE, LCIO::SIMTRACKERHIT); + + for (int ihit = 0; ihit < hit_col->getNumberOfElements(); ++ihit) { + if (TrackerHitPlane* hit = dynamic_cast(hit_col->getElementAt(ihit))) { + // Skipping the hit if its time is outside the acceptance time window + double hitT = hit->getTime(); + + dd4hep::rec::Vector3D pos = hit->getPosition(); + double hitR = pos.r(); + + // Correcting for the propagation time + double dt = hitR / (TMath::C() * m_beta / 1e6); + hitT -= dt; + streamlog_out(DEBUG3) << "corrected hit at R: " << hitR << " mm by propagation time: " << dt + << " ns to T: " << hitT << " ns" << std::endl; + + // Apply time window selection + if (hitT < m_time_min || hitT > m_time_max) { + streamlog_out(DEBUG4) << "hit at T: " << hitT << " ns is rejected by timing cuts" << std::endl; + continue; + } - LCRelationImpl *rel_new = new LCRelationImpl(); + // hits_to_save[icol].insert(ihit); + outputTrackerHitCol->addElement(hit); - rel_new->setFrom(hit_new); - rel_new->setTo(simhit_new); - rel_new->setWeight(rel->getWeight()); + LCRelation* rel = static_cast(inputHitRels[icol]->getElementAt(ihit)); + SimTrackerHit* simhit = static_cast(rel->getTo()); + outputTrackerSimHitCol->addElement(simhit); - outputTrackerHitRels[icol]->addElement(rel_new); + thitNav.addRelation(hit, simhit); - } // ihit loop + if (m_fillHistos) + m_corrected_time->Fill(hitT); + } - streamlog_out(MESSAGE) << " " << hits_to_save[icol].size() << " hits added to the collections: " - << m_outputTrackerHitsCollNames[icol] << ", " - << m_outputTrackerSimHitsCollNames[icol] << ", " - << m_outputTrackerHitRelNames[icol] << std::endl; + } // ihit loop - evt->addCollection(outputTrackerHitColls[icol], m_outputTrackerHitsCollNames[icol]); - evt->addCollection(outputTrackerSimHitColls[icol], m_outputTrackerSimHitsCollNames[icol]); - evt->addCollection(outputTrackerHitRels[icol], m_outputTrackerHitRelNames[icol]); + streamlog_out(MESSAGE) << " " << outputTrackerHitCol->getNumberOfElements() + << " hits added to the collections: " << m_outputTrackerHitsCollNames[icol] << ", " + << m_outputTrackerSimHitsCollNames[icol] << ", " << m_outputTrackerHitRelNames[icol] + << std::endl; - streamlog_out(DEBUG5) << " output collection " << m_outputTrackerHitsCollNames[icol] << " of type " - << outputTrackerHitColls[icol]->getTypeName() << " added to the event \n" - << " output collection " << m_outputTrackerSimHitsCollNames[icol] << " of type " - << outputTrackerSimHitColls[icol]->getTypeName() << " added to the event \n" - << " output collection " << m_outputTrackerHitRelNames[icol] << " of type " - << outputTrackerHitRels[icol]->getTypeName() << " added to the event " - << std::endl; + evt->addCollection(outputTrackerHitCol, m_outputTrackerHitsCollNames[icol]); + evt->addCollection(outputTrackerSimHitCol, m_outputTrackerSimHitsCollNames[icol]); + outputTrackerHitRel = thitNav.createLCCollection(); + evt->addCollection(outputTrackerHitRel, m_outputTrackerHitRelNames[icol]); - } // icol loop + streamlog_out(DEBUG) << " output collection " << m_outputTrackerHitsCollNames[icol] << " of type " + << outputTrackerHitCol->getTypeName() << " added to the event \n" + << " output collection " << m_outputTrackerSimHitsCollNames[icol] << " of type " + << outputTrackerSimHitCol->getTypeName() << " added to the event \n" + << " output collection " << m_outputTrackerHitRelNames[icol] << " of type " + << outputTrackerHitRel->getTypeName() << " added to the event " << std::endl; - streamlog_out(DEBUG) << " Event processed " << std::endl; + } // icol loop - _nEvt++; -} + streamlog_out(DEBUG) << " Event processed " << std::endl; -void FilterTimeHits::check(LCEvent *) -{ + _nEvt++; } -void FilterTimeHits::end() -{ +void FilterTimeHits::check(LCEvent*) {} - std::cout << "FilterTimeHits::end() " << name() - << " processed " << _nEvt << " events in " << _nRun << " runs " - << std::endl; +void FilterTimeHits::end() { + std::cout << "FilterTimeHits::end() " << name() << " processed " << _nEvt << " events in " << _nRun << " runs " + << std::endl; } diff --git a/source/Utils/src/FilterTracks.cc b/source/Utils/src/FilterTracks.cc index 4074751..5dda2f9 100644 --- a/source/Utils/src/FilterTracks.cc +++ b/source/Utils/src/FilterTracks.cc @@ -1,4 +1,4 @@ -# include "FilterTracks.h" +#include "FilterTracks.h" #include @@ -10,137 +10,86 @@ #include #include -FilterTracks aFilterTracks ; +FilterTracks aFilterTracks; -FilterTracks::FilterTracks() - : Processor("FilterTracks") -{ +FilterTracks::FilterTracks() : Processor("FilterTracks") { // modify processor description - _description = "FilterTracks processor filters a collection of tracks based on NHits and MinPt and outputs a filtered collection"; + _description = "FilterTracks processor filters a collection of tracks based on NHits and MinPt and outputs a " + "filtered collection"; // register steering parameters: name, description, class-variable, default value - registerProcessorParameter("BarrelOnly", - "If true, just keep tracks with only barrel hits", - _BarrelOnly, - _BarrelOnly - ); - + registerProcessorParameter("BarrelOnly", "If true, just keep tracks with only barrel hits", _BarrelOnly, _BarrelOnly); + registerProcessorParameter("HasCaloState", - "If true, just keep tracks that have a TrackState at the Calorimeter surface", - _HasCaloState, - _HasCaloState - ); - - registerProcessorParameter("NHitsTotal", - "Minimum number of hits on track", - _NHitsTotal, - _NHitsTotal - ); - - registerProcessorParameter("NHitsVertex", - "Minimum number of hits on vertex detector", - _NHitsVertex, - _NHitsVertex - ); - - registerProcessorParameter("NHitsInner", - "Minimum number of hits on inner tracker", - _NHitsInner, - _NHitsInner - ); - - registerProcessorParameter("NHitsOuter", - "Minimum number of hits on outer tracker", - _NHitsOuter, - _NHitsOuter - ); - - registerProcessorParameter("MaxHoles", - "Maximum number of holes on track", - _MaxHoles, - _MaxHoles - ); - - registerProcessorParameter("MinPt", - "Minimum transverse momentum", - _MinPt, - _MinPt - ); - - registerProcessorParameter("Chi2Spatial", - "Spatial chi squared", - _Chi2Spatial, - _Chi2Spatial - ); - - registerInputCollection( LCIO::TRACK, - "InputTrackCollectionName" , - "Name of the input collection", - _InputTrackCollection, - _InputTrackCollection - ); - - registerOutputCollection( LCIO::TRACK, - "OutputTrackCollectionName" , - "Name of output collection", - _OutputTrackCollection, - std::string("FilteredTracks") - ); + "If true, just keep tracks that have a TrackState at the Calorimeter surface", + _HasCaloState, _HasCaloState); + + registerProcessorParameter("NHitsTotal", "Minimum number of hits on track", _NHitsTotal, _NHitsTotal); + + registerProcessorParameter("NHitsVertex", "Minimum number of hits on vertex detector", _NHitsVertex, _NHitsVertex); + + registerProcessorParameter("NHitsInner", "Minimum number of hits on inner tracker", _NHitsInner, _NHitsInner); + + registerProcessorParameter("NHitsOuter", "Minimum number of hits on outer tracker", _NHitsOuter, _NHitsOuter); + + registerProcessorParameter("MaxHoles", "Maximum number of holes on track", _MaxHoles, _MaxHoles); + registerProcessorParameter("MinPt", "Minimum transverse momentum", _MinPt, _MinPt); + + registerProcessorParameter("Chi2Spatial", "Spatial chi squared", _Chi2Spatial, _Chi2Spatial); + + registerInputCollection(LCIO::TRACK, "InputTrackCollectionName", "Name of the input collection", + _InputTrackCollection, _InputTrackCollection); + + registerOutputCollection(LCIO::TRACK, "OutputTrackCollectionName", "Name of output collection", + _OutputTrackCollection, std::string("FilteredTracks")); } -void FilterTracks::init() -{ +void FilterTracks::init() { // Print the initial parameters - printParameters() ; - buildBfield() ; + printParameters(); + buildBfield(); } -void FilterTracks::processRunHeader( LCRunHeader* /*run*/) -{ } +void FilterTracks::processRunHeader(LCRunHeader* /*run*/) {} -void FilterTracks::buildBfield() -{ +void FilterTracks::buildBfield() { // Get the magnetic field dd4hep::Detector& lcdd = dd4hep::Detector::getInstance(); - const double position[3] = { - 0, 0, - 0}; // position to calculate magnetic field at (the origin in this case) - double magneticFieldVector[3] = { - 0, 0, 0}; // initialise object to hold magnetic field - lcdd.field().magneticField( - position, - magneticFieldVector); // get the magnetic field vector from DD4hep - _Bz = magneticFieldVector[2]/dd4hep::tesla; + const double position[3] = {0, 0, 0}; // position to calculate magnetic field at (the origin in this case) + double magneticFieldVector[3] = {0, 0, 0}; // initialise object to hold magnetic field + lcdd.field().magneticField(position, + magneticFieldVector); // get the magnetic field vector from DD4hep + _Bz = magneticFieldVector[2] / dd4hep::tesla; } -void FilterTracks::processEvent( LCEvent * evt ) -{ +void FilterTracks::processEvent(LCEvent* evt) { // Make the output track collection - LCCollectionVec *OutputTrackCollection = new LCCollectionVec(LCIO::TRACK); + LCCollectionVec* OutputTrackCollection = new LCCollectionVec(LCIO::TRACK); OutputTrackCollection->setSubset(true); // Get input collection - LCCollection* InputTrackCollection =evt->getCollection(_InputTrackCollection); + LCCollection* InputTrackCollection = evt->getCollection(_InputTrackCollection); - if( InputTrackCollection->getTypeName() != lcio::LCIO::TRACK ) - { throw EVENT::Exception( "Invalid collection type: " + InputTrackCollection->getTypeName() ) ; } + if (InputTrackCollection->getTypeName() != lcio::LCIO::TRACK) { + throw EVENT::Exception("Invalid collection type: " + InputTrackCollection->getTypeName()); + } // Filter std::string encoderString = lcio::LCTrackerCellID::encoding_string(); UTIL::CellIDDecoder decoder(encoderString); - for(int i=0; igetNumberOfElements(); i++) { - EVENT::Track *trk=dynamic_cast(InputTrackCollection->getElementAt(i)); + for (int i = 0; i < InputTrackCollection->getNumberOfElements(); i++) { + EVENT::Track* trk = dynamic_cast(InputTrackCollection->getElementAt(i)); - int nhittotal = trk->getTrackerHits().size(); + int nhittotal = trk->getTrackerHits().size(); const EVENT::IntVec& subdetectorHits = trk->getSubdetectorHitNumbers(); - int nhitvertex = subdetectorHits[1]+subdetectorHits[2]; - int nhitinner = subdetectorHits[3]+subdetectorHits[4]; - int nhitouter = subdetectorHits[5]+subdetectorHits[6]; + int nhitvertex = subdetectorHits[1] + subdetectorHits[2]; + int nhitinner = subdetectorHits[3] + subdetectorHits[4]; + int nhitouter = subdetectorHits[5] + subdetectorHits[6]; - float pt = fabs(0.3*_Bz/trk->getOmega()/1000); + float pt = fabs(0.3 * _Bz / trk->getOmega() / 1000); float chi2spatial = trk->getChi2(); @@ -148,38 +97,37 @@ void FilterTracks::processEvent( LCEvent * evt ) // Check if a TrackState at the calo surface exists const std::vector& trackStates = trk->getTrackStates(); - const auto foundCaloState = std::find_if(trackStates.begin(), trackStates.end(), - [](const auto ts) { return ts->getLocation() == EVENT::TrackState::AtCalorimeter; }) != trackStates.end(); - if (_HasCaloState && !foundCaloState) { continue; } + const auto foundCaloState = std::find_if(trackStates.begin(), trackStates.end(), [](const auto ts) { + return ts->getLocation() == EVENT::TrackState::AtCalorimeter; + }) != trackStates.end(); + if (_HasCaloState && !foundCaloState) { + streamlog_out(DEBUG) << "No calo state, skipping track!" << std::endl; + continue; + } - if(_BarrelOnly == true) { + if (_BarrelOnly == true) { bool endcaphits = false; - for(int j=0; jgetTrackerHits()[j])["system"]; - if(systemID == 2 || systemID == 4 || systemID == 6) { + if (systemID == 2 || systemID == 4 || systemID == 6) { endcaphits = true; break; } } - if(endcaphits == false) { OutputTrackCollection->addElement(trk); } + if (endcaphits == false) { + OutputTrackCollection->addElement(trk); + } } else { // track property cuts - if(nhittotal > _NHitsTotal && - nhitvertex > _NHitsVertex && - nhitinner > _NHitsInner && - nhitouter > _NHitsOuter && - pt > _MinPt && - chi2spatial > _Chi2Spatial && - nholes <= _MaxHoles) - { - OutputTrackCollection->addElement(trk); + if (nhittotal > _NHitsTotal && nhitvertex > _NHitsVertex && nhitinner > _NHitsInner && nhitouter > _NHitsOuter && + pt > _MinPt && chi2spatial > _Chi2Spatial && nholes <= _MaxHoles) { + OutputTrackCollection->addElement(trk); } } } // Save output track collection - evt->addCollection(OutputTrackCollection, _OutputTrackCollection); + evt->addCollection(OutputTrackCollection, _OutputTrackCollection); } -void FilterTracks::end() -{ } +void FilterTracks::end() {} diff --git a/source/Utils/src/SplitCollectionByLayer.cc b/source/Utils/src/SplitCollectionByLayer.cc index 6ad2133..57f2b12 100644 --- a/source/Utils/src/SplitCollectionByLayer.cc +++ b/source/Utils/src/SplitCollectionByLayer.cc @@ -1,90 +1,68 @@ #include "SplitCollectionByLayer.h" -#include -#include #include +#include +#include -#include +#include +#include #include #include #include -#include -#include +#include #include // #include "UTIL/LCTrackerConf.h" // #include +using namespace lcio; +using namespace marlin; -using namespace lcio ; -using namespace marlin ; - - -SplitCollectionByLayer aSplitCollectionByLayer ; - - - -SplitCollectionByLayer::SplitCollectionByLayer() : Processor("SplitCollectionByLayer") , - _nRun(0), _nEvt(0) { +SplitCollectionByLayer aSplitCollectionByLayer; +SplitCollectionByLayer::SplitCollectionByLayer() : Processor("SplitCollectionByLayer"), _nRun(0), _nEvt(0) { // modify processor description - _description = "split a hit collection based on the layer number of the hits " ; - + _description = "split a hit collection based on the layer number of the hits "; // register steering parameters: name, description, class-variable, default value - registerProcessorParameter("InputCollection" , - "Name of the input collection with hits" , - _colName , - std::string("FTDCollection") - ); - - StringVec outColEx ; - outColEx.push_back("FTD_PIXELCollection") ; - outColEx.push_back("0") ; - outColEx.push_back("1") ; - outColEx.push_back("FTD_STRIPCollection") ; - outColEx.push_back("2") ; - outColEx.push_back("6") ; - - registerProcessorParameter("OutputCollections" , - "Name of the output collection with start and end layer number" , - _outColAndLayers , - outColEx - ); - - registerProcessorParameter("KeepEmptyCollections" , - "Whether collections should be added to the event even if they are empty" , - _addEmptyCollections , - false - ); - - + registerProcessorParameter("InputCollection", "Name of the input collection with hits", _colName, + std::string("FTDCollection")); + + StringVec outColEx; + outColEx.push_back("FTD_PIXELCollection"); + outColEx.push_back("0"); + outColEx.push_back("1"); + outColEx.push_back("FTD_STRIPCollection"); + outColEx.push_back("2"); + outColEx.push_back("6"); + + registerProcessorParameter("OutputCollections", "Name of the output collection with start and end layer number", + _outColAndLayers, outColEx); + + registerProcessorParameter("KeepEmptyCollections", + "Whether collections should be added to the event even if they are empty", + _addEmptyCollections, false); } - template -long cellIDFromHit( const LCObject* o){ - long id = -1 ; - const T* h = dynamic_cast( o ) ; - if( h ) id = ( h->getCellID0() & 0xffffffff ) | ( ( long(h->getCellID1()) << 32 ) & 0xffffffff00000000 ) ; - return id ; +long cellIDFromHit(const LCObject* o) { + long id = -1; + const T* h = dynamic_cast(o); + if (h) + id = (h->getCellID0() & 0xffffffff) | ((long(h->getCellID1()) << 32) & 0xffffffff00000000); + return id; } - - void SplitCollectionByLayer::init() { - - streamlog_out(DEBUG) << " init called " << std::endl ; + streamlog_out(DEBUG) << " init called " << std::endl; // usually a good idea to - printParameters() ; - + printParameters(); - unsigned i=0; - while( i < _outColAndLayers.size() ){ - - std::string name = _outColAndLayers[ i++ ] ; - size_t layerStart = std::atoi( _outColAndLayers[ i++ ].c_str() ) ; - size_t layerEnd = std::atoi( _outColAndLayers[ i++ ].c_str() ) ; + unsigned i = 0; + while (i < _outColAndLayers.size()) { + std::string name = _outColAndLayers[i++]; + size_t layerStart = std::atoi(_outColAndLayers[i++].c_str()); + size_t layerEnd = std::atoi(_outColAndLayers[i++].c_str()); // Ensuring that a collection setup is created for this name if (_outCols.find(name) == _outCols.end()) { @@ -97,158 +75,132 @@ void SplitCollectionByLayer::init() { } } - streamlog_out(MESSAGE) << "Will split collection " << _colName << " by layers into the following collections:" << std::endl; - for(const auto &outCol : _outCols){ + streamlog_out(MESSAGE) << "Will split collection " << _colName + << " by layers into the following collections:" << std::endl; + for (const auto& outCol : _outCols) { streamlog_out(MESSAGE) << " " << outCol.first << ": "; - for(const auto &layer : outCol.second.layers) { + for (const auto& layer : outCol.second.layers) { streamlog_out(MESSAGE) << layer << " "; } streamlog_out(MESSAGE) << std::endl; } - _nRun = 0 ; - _nEvt = 0 ; - + _nRun = 0; + _nEvt = 0; } +void SplitCollectionByLayer::processRunHeader(LCRunHeader*) { _nRun++; } -void SplitCollectionByLayer::processRunHeader( LCRunHeader* ) { - - _nRun++ ; -} - +void SplitCollectionByLayer::processEvent(LCEvent* evt) { + LCCollection* col = 0; + try { + col = evt->getCollection(_colName); - -void SplitCollectionByLayer::processEvent( LCEvent * evt ) { - - - LCCollection* col = 0 ; - - try{ col = evt->getCollection( _colName ) ; - - } catch( lcio::DataNotAvailableException& e) { - - streamlog_out( DEBUG5 ) << " input collection not in event : " << _colName << " - nothing to do !!! " << std::endl ; - return ; + } catch (lcio::DataNotAvailableException& e) { + streamlog_out(DEBUG5) << " input collection not in event : " << _colName << " - nothing to do !!! " + << std::endl; + return; } // remember the type of the hit collection - if( col->getTypeName() == lcio::LCIO::SIMTRACKERHIT ) - _type = SimTrackerHitType ; - else if( col->getTypeName() == lcio::LCIO::TRACKERHIT ) - _type = TrackerHitType ; - else if( col->getTypeName() == lcio::LCIO::TRACKERHITPLANE ) - _type = TrackerHitPlaneType ; - else if( col->getTypeName() == lcio::LCIO::SIMCALORIMETERHIT ) - _type = SimCalorimeterHitType ; - else if( col->getTypeName() == lcio::LCIO::CALORIMETERHIT ) - _type = CalorimeterHitType ; + if (col->getTypeName() == lcio::LCIO::SIMTRACKERHIT) + _type = SimTrackerHitType; + else if (col->getTypeName() == lcio::LCIO::TRACKERHIT) + _type = TrackerHitType; + else if (col->getTypeName() == lcio::LCIO::TRACKERHITPLANE) + _type = TrackerHitPlaneType; + else if (col->getTypeName() == lcio::LCIO::SIMCALORIMETERHIT) + _type = SimCalorimeterHitType; + else if (col->getTypeName() == lcio::LCIO::CALORIMETERHIT) + _type = CalorimeterHitType; else - _type = UnkownType ; - + _type = UnkownType; - std::string encoderString = col->getParameters().getStringVal( "CellIDEncoding" ) ; + std::string encoderString = col->getParameters().getStringVal("CellIDEncoding"); - UTIL::BitField64 encoder( encoderString ) ; - - unsigned layerIndex = encoder.index("layer") ; + UTIL::BitField64 encoder(encoderString); + unsigned layerIndex = encoder.index("layer"); //---- create output collections - for(auto &outCol : _outCols){ - - LCCollectionVec* newCol = new LCCollectionVec( col->getTypeName() ) ; + for (auto& outCol : _outCols) { + LCCollectionVec* newCol = new LCCollectionVec(col->getTypeName()); - newCol->setSubset( true ) ; + newCol->setSubset(true); - newCol->parameters().setValue( "CellIDEncoding", encoderString ) ; + newCol->parameters().setValue("CellIDEncoding", encoderString); - outCol.second.collection = newCol ; - - streamlog_out( DEBUG5 ) << " created new output collection " << outCol.first << " of type " << col->getTypeName() << std::endl ; - } + outCol.second.collection = newCol; + streamlog_out(DEBUG5) << " created new output collection " << outCol.first << " of type " << col->getTypeName() + << std::endl; + } //---- loop over hits - int nHit = col->getNumberOfElements() ; + int nHit = col->getNumberOfElements(); - for(int iHit=0; iHit< nHit ; iHit++){ + for (int iHit = 0; iHit < nHit; iHit++) { + lcio::LCObject* h = col->getElementAt(iHit); - lcio::LCObject* h = col->getElementAt( iHit ) ; - - long id = -1 ; - - switch( _type ){ + long id = -1; + switch (_type) { case SimTrackerHitType: - id = cellIDFromHit( h ) ; - break ; + id = cellIDFromHit(h); + break; case TrackerHitType: - id = cellIDFromHit( h ) ; - break ; + id = cellIDFromHit(h); + break; case TrackerHitPlaneType: - id = cellIDFromHit( h ) ; - break ; + id = cellIDFromHit(h); + break; case SimCalorimeterHitType: - id = cellIDFromHit( h ) ; - break ; + id = cellIDFromHit(h); + break; case CalorimeterHitType: - id = cellIDFromHit( h ) ; - break ; + id = cellIDFromHit(h); + break; case UnkownType: - continue ; + continue; } - encoder.setValue( id ) ; - - size_t layerID = encoder[ layerIndex ] ; + encoder.setValue(id); + size_t layerID = encoder[layerIndex]; // check if we have an output collection for this layer - for(auto &outCol : _outCols){ - if (std::find(outCol.second.layers.begin(), outCol.second.layers.end(), layerID) == outCol.second.layers.end()) continue; - outCol.second.collection->addElement( h ) ; - streamlog_out( DEBUG0 ) << " adding hit for layer " << layerID << " to collection : " << outCol.first << std::endl ; + for (auto& outCol : _outCols) { + if (std::find(outCol.second.layers.begin(), outCol.second.layers.end(), layerID) == outCol.second.layers.end()) + continue; + outCol.second.collection->addElement(h); + streamlog_out(DEBUG0) << " adding hit for layer " << layerID << " to collection : " << outCol.first << std::endl; } } - // add non empty collections to the event - for(const auto &outCol : _outCols){ - - LCCollection* newCol = outCol.second.collection ; + for (const auto& outCol : _outCols) { + LCCollection* newCol = outCol.second.collection; - if( _addEmptyCollections || newCol->getNumberOfElements() > 0 ) { - evt->addCollection( newCol , outCol.first ) ; - streamlog_out( DEBUG5 ) << " Output collection " << outCol.first << " of type " << col->getTypeName() << " added to the event" << std::endl ; + if (_addEmptyCollections || newCol->getNumberOfElements() > 0) { + evt->addCollection(newCol, outCol.first); + streamlog_out(DEBUG5) << " Output collection " << outCol.first << " of type " << col->getTypeName() + << " added to the event" << std::endl; } } + streamlog_out(DEBUG3) << " processing event: " << evt->getEventNumber() << " in run: " << evt->getRunNumber() + << std::endl; - - - - streamlog_out( DEBUG3 ) << " processing event: " << evt->getEventNumber() - << " in run: " << evt->getRunNumber() << std::endl ; - - - _nEvt ++ ; + _nEvt++; } - - -void SplitCollectionByLayer::check( LCEvent* ) { +void SplitCollectionByLayer::check(LCEvent*) { // nothing to check here - could be used to fill checkplots in reconstruction processor } - -void SplitCollectionByLayer::end(){ - - streamlog_out( MESSAGE ) << "SplitCollectionByLayer::end() " << name() - << " processed " << _nEvt << " events in " << _nRun << " runs " - << std::endl ; - +void SplitCollectionByLayer::end() { + streamlog_out(MESSAGE) << "SplitCollectionByLayer::end() " << name() << " processed " << _nEvt << " events in " + << _nRun << " runs " << std::endl; } - diff --git a/source/Utils/src/SplitCollectionByPolarAngle.cc b/source/Utils/src/SplitCollectionByPolarAngle.cc index 04c7a9d..7497bc5 100644 --- a/source/Utils/src/SplitCollectionByPolarAngle.cc +++ b/source/Utils/src/SplitCollectionByPolarAngle.cc @@ -1,148 +1,98 @@ #include "SplitCollectionByPolarAngle.h" -#include #include +#include #include - #include -#include -#include -#include #include +#include +#include +#include #include #include #include "HelixClass_double.h" -using namespace lcio ; -using namespace marlin ; - - -SplitCollectionByPolarAngle aSplitCollectionByPolarAngle ; +using namespace lcio; +using namespace marlin; +SplitCollectionByPolarAngle aSplitCollectionByPolarAngle; SplitCollectionByPolarAngle::SplitCollectionByPolarAngle() : Processor("SplitCollectionByPolarAngle") { - // --- Processor description: _description = "SplitCollectionByPolarAngle selects tracker hits based on their polar angle"; - // --- Processor parameters: + registerProcessorParameter("TrackerHitInputCollections", "Name of the tracker hit input collections", + m_inputTrackerHitsCollNames, {}); - registerProcessorParameter("TrackerHitInputCollections", - "Name of the tracker hit input collections", - m_inputTrackerHitsCollNames, - {} ); - - registerProcessorParameter("TrackerSimHitInputCollections", - "Name of the tracker simhit input collections", - m_inputTrackerSimHitsCollNames, - {} ); - - registerProcessorParameter("TrackerHitInputRelations", - "Name of the tracker hit relation collections", - m_inputTrackerHitRelNames, - {} ); - - registerProcessorParameter("TrackerHitOutputCollections", - "Name of the tracker hit output collections", - m_outputTrackerHitsCollNames, - {} ); - - registerProcessorParameter("TrackerSimHitOutputCollections", - "Name of the tracker simhit output collections", - m_outputTrackerSimHitsCollNames, - {} ); + registerProcessorParameter("TrackerSimHitInputCollections", "Name of the tracker simhit input collections", + m_inputTrackerSimHitsCollNames, {}); - registerProcessorParameter("TrackerHitOutputRelations", - "Name of the tracker hit relation collections", - m_outputTrackerHitRelNames, - {} ); + registerProcessorParameter("TrackerHitInputRelations", "Name of the tracker hit relation collections", + m_inputTrackerHitRelNames, {}); + registerProcessorParameter("TrackerHitOutputCollections", "Name of the tracker hit output collections", + m_outputTrackerHitsCollNames, {}); - registerProcessorParameter( "FillHistograms", - "Flag to fill the diagnostic histograms", - m_fillHistos, - false ); + registerProcessorParameter("TrackerSimHitOutputCollections", "Name of the tracker simhit output collections", + m_outputTrackerSimHitsCollNames, {}); - registerProcessorParameter( "PolarAngleLowerLimit", - "Lower limit on the hit polar angle in degrees", - m_theta_min, - double(50.0) ); + registerProcessorParameter("TrackerHitOutputRelations", "Name of the tracker hit relation collections", + m_outputTrackerHitRelNames, {}); - registerProcessorParameter( "PolarAngleUpperLimit", - "Upper limit on the hit polar angle in degrees", - m_theta_max, - double(130.0) ); + registerProcessorParameter("FillHistograms", "Flag to fill the diagnostic histograms", m_fillHistos, false); + registerProcessorParameter("PolarAngleLowerLimit", "Lower limit on the hit polar angle in degrees", m_theta_min, + double(50.0)); + registerProcessorParameter("PolarAngleUpperLimit", "Upper limit on the hit polar angle in degrees", m_theta_max, + double(130.0)); } - - void SplitCollectionByPolarAngle::init() { - - streamlog_out(DEBUG) << " init called " << std::endl ; + streamlog_out(DEBUG) << " init called " << std::endl; // --- Print the processor parameters: - printParameters() ; - - + printParameters(); // --- Initialize the run and event counters: - _nRun = 0 ; - _nEvt = 0 ; - + _nRun = 0; + _nEvt = 0; // --- Initialize the AIDAProcessor and book the diagnostic histograms: AIDAProcessor::histogramFactory(this); m_theta = new TH1F("m_theta", "polar angle of the hit [rad]", 1000, 0., M_PI); - -} - - -void SplitCollectionByPolarAngle::processRunHeader( LCRunHeader* ) { - - _nRun++ ; - } +void SplitCollectionByPolarAngle::processRunHeader(LCRunHeader*) { _nRun++; } - -void SplitCollectionByPolarAngle::processEvent( LCEvent * evt) { - +void SplitCollectionByPolarAngle::processEvent(LCEvent* evt) { // --- Check whether the number of input and output collections match - if ( m_inputTrackerHitsCollNames.size() != m_inputTrackerSimHitsCollNames.size() || - m_inputTrackerHitsCollNames.size() != m_inputTrackerHitRelNames.size() ){ - + if (m_inputTrackerHitsCollNames.size() != m_inputTrackerSimHitsCollNames.size() || + m_inputTrackerHitsCollNames.size() != m_inputTrackerHitRelNames.size()) { std::stringstream err_msg; - err_msg << "Mismatch between the reco and sim hits input collections" - << std::endl ; - - throw EVENT::Exception( err_msg.str() ) ; + err_msg << "Mismatch between the reco and sim hits input collections" << std::endl; + throw EVENT::Exception(err_msg.str()); } - if ( m_outputTrackerHitsCollNames.size() != m_outputTrackerSimHitsCollNames.size() || - m_outputTrackerHitsCollNames.size() != m_outputTrackerHitRelNames.size() ){ - + if (m_outputTrackerHitsCollNames.size() != m_outputTrackerSimHitsCollNames.size() || + m_outputTrackerHitsCollNames.size() != m_outputTrackerHitRelNames.size()) { std::stringstream err_msg; - err_msg << "Mismatch between the reco and sim hits output collections" - << std::endl ; - - throw EVENT::Exception( err_msg.str() ) ; + err_msg << "Mismatch between the reco and sim hits output collections" << std::endl; + throw EVENT::Exception(err_msg.str()); } - // --- Get the input hit collections and create the corresponding output collections: const unsigned int nTrackerHitCol = m_inputTrackerHitsCollNames.size(); @@ -154,106 +104,93 @@ void SplitCollectionByPolarAngle::processEvent( LCEvent * evt) { std::vector outputTrackerSimHitColls(nTrackerHitCol); std::vector outputTrackerHitRels(nTrackerHitCol); - for (unsigned int icol=0; icolgetCollection(m_inputTrackerHitsCollNames[icol]); - } - catch( lcio::DataNotAvailableException& e ) { - streamlog_out(WARNING) << m_inputTrackerHitsCollNames[icol] - << " collection not available" << std::endl; + } catch (lcio::DataNotAvailableException& e) { + streamlog_out(WARNING) << m_inputTrackerHitsCollNames[icol] << " collection not available" << std::endl; continue; } // get the sim hits try { inputSimHitColls[icol] = evt->getCollection(m_inputTrackerSimHitsCollNames[icol]); - } - catch( lcio::DataNotAvailableException& e ) { - streamlog_out(WARNING) << m_inputTrackerSimHitsCollNames[icol] - << " collection not available" << std::endl; + } catch (lcio::DataNotAvailableException& e) { + streamlog_out(WARNING) << m_inputTrackerSimHitsCollNames[icol] << " collection not available" << std::endl; continue; } // get the reco-sim relations try { inputHitRels[icol] = evt->getCollection(m_inputTrackerHitRelNames[icol]); - } - catch( lcio::DataNotAvailableException& e ) { - streamlog_out(WARNING) << m_inputTrackerHitRelNames[icol] - << " collection not available" << std::endl; + } catch (lcio::DataNotAvailableException& e) { + streamlog_out(WARNING) << m_inputTrackerHitRelNames[icol] << " collection not available" << std::endl; continue; } // reco hit output collections - std::string encoderString = inputHitColls[icol]->getParameters().getStringVal( "CellIDEncoding" ); - outputTrackerHitColls[icol] = new LCCollectionVec( inputHitColls[icol]->getTypeName() ); - outputTrackerHitColls[icol]->parameters().setValue( "CellIDEncoding", encoderString ); + std::string encoderString = inputHitColls[icol]->getParameters().getStringVal("CellIDEncoding"); + outputTrackerHitColls[icol] = new LCCollectionVec(inputHitColls[icol]->getTypeName()); + outputTrackerHitColls[icol]->parameters().setValue("CellIDEncoding", encoderString); LCFlagImpl lcFlag(inputHitColls[icol]->getFlag()); outputTrackerHitColls[icol]->setFlag(lcFlag.getFlag()); - + // sim hit output collections - outputTrackerSimHitColls[icol] = new LCCollectionVec( inputSimHitColls[icol]->getTypeName() ); - outputTrackerSimHitColls[icol]->parameters().setValue( "CellIDEncoding", encoderString ); + outputTrackerSimHitColls[icol] = new LCCollectionVec(inputSimHitColls[icol]->getTypeName()); + outputTrackerSimHitColls[icol]->parameters().setValue("CellIDEncoding", encoderString); LCFlagImpl lcFlag_sim(inputSimHitColls[icol]->getFlag()); outputTrackerSimHitColls[icol]->setFlag(lcFlag_sim.getFlag()); // reco-sim relation output collections - outputTrackerHitRels[icol] = new LCCollectionVec( inputHitRels[icol]->getTypeName() ); - LCFlagImpl lcFlag_rel(inputHitRels[icol]->getFlag()) ; - outputTrackerHitRels[icol]->setFlag( lcFlag_rel.getFlag() ) ; - + outputTrackerHitRels[icol] = new LCCollectionVec(inputHitRels[icol]->getTypeName()); + LCFlagImpl lcFlag_rel(inputHitRels[icol]->getFlag()); + outputTrackerHitRels[icol]->setFlag(lcFlag_rel.getFlag()); } - // --- Loop over the tracker hits and select hits inside a cone around the particle trajectory: - std::vector > hits_to_save(nTrackerHitCol); - - - for (unsigned int icol=0; icol> hits_to_save(nTrackerHitCol); - LCCollection* hit_col = inputHitColls[icol]; - if( !hit_col ) continue ; - - for (int ihit=0; ihitgetNumberOfElements(); ++ihit){ + for (unsigned int icol = 0; icol < inputHitColls.size(); ++icol) { + LCCollection* hit_col = inputHitColls[icol]; + if (!hit_col) + continue; + for (int ihit = 0; ihit < hit_col->getNumberOfElements(); ++ihit) { TrackerHitPlane* hit = dynamic_cast(hit_col->getElementAt(ihit)); // --- Get the modulus of the position-vector of the hit : - double hit_position_mod = sqrt( hit->getPosition()[0]*hit->getPosition()[0] + hit->getPosition()[1]*hit->getPosition()[1] + hit->getPosition()[2]*hit->getPosition()[2] ); + double hit_position_mod = + sqrt(hit->getPosition()[0] * hit->getPosition()[0] + hit->getPosition()[1] * hit->getPosition()[1] + + hit->getPosition()[2] * hit->getPosition()[2]); // --- Get the polar angle theta of the hit : - double hit_theta = acos( hit->getPosition()[2] / hit_position_mod ); + double hit_theta = acos(hit->getPosition()[2] / hit_position_mod); // std::cout << "HIT POSITION AND ANGLE #################" << std::endl; - // std::cout << "(x,y,z) = " << "(" << hit->getPosition()[0] << "," << hit->getPosition()[1] << "," << hit->getPosition()[2] << ")" << std::endl; - // std::cout << "theta [rad] = " << hit_theta << " ; theta [deg] = " << hit_theta*180/M_PI << std::endl; - - + // std::cout << "(x,y,z) = " << "(" << hit->getPosition()[0] << "," << hit->getPosition()[1] << "," << + // hit->getPosition()[2] << ")" << std::endl; std::cout << "theta [rad] = " << hit_theta << " ; theta [deg] = " << + // hit_theta*180/M_PI << std::endl; // --- Add the filtered hits to the output collections : - if (hit_theta*180/M_PI < m_theta_min || hit_theta*180/M_PI > m_theta_max) continue; + if (hit_theta * 180 / M_PI < m_theta_min || hit_theta * 180 / M_PI > m_theta_max) + continue; hits_to_save[icol].insert(ihit); - if ( m_fillHistos ) m_theta->Fill(hit_theta); - - + if (m_fillHistos) + m_theta->Fill(hit_theta); } // ihit loop } // icol loop - // --- Add the filtered hits to the output collections: - for (unsigned int icol=0; icol(inputHitColls[icol]->getElementAt(ihit)); TrackerHitPlaneImpl* hit_new = new TrackerHitPlaneImpl(); @@ -270,12 +207,10 @@ void SplitCollectionByPolarAngle::processEvent( LCEvent * evt) { hit_new->setTime(hit->getTime()); hit_new->setQuality(hit->getQuality()); - outputTrackerHitColls[icol]->addElement( hit_new ); - + outputTrackerHitColls[icol]->addElement(hit_new); LCRelation* rel = dynamic_cast(inputHitRels[icol]->getElementAt(ihit)); - SimTrackerHit* simhit = dynamic_cast(rel->getTo()); SimTrackerHitImpl* simhit_new = new SimTrackerHitImpl(); @@ -291,8 +226,7 @@ void SplitCollectionByPolarAngle::processEvent( LCEvent * evt) { simhit_new->setOverlay(simhit->isOverlay()); simhit_new->setProducedBySecondary(simhit->isProducedBySecondary()); - outputTrackerSimHitColls[icol]->addElement( simhit_new ); - + outputTrackerSimHitColls[icol]->addElement(simhit_new); LCRelationImpl* rel_new = new LCRelationImpl(); @@ -300,47 +234,37 @@ void SplitCollectionByPolarAngle::processEvent( LCEvent * evt) { rel_new->setTo(simhit_new); rel_new->setWeight(rel->getWeight()); - outputTrackerHitRels[icol]->addElement( rel_new ); + outputTrackerHitRels[icol]->addElement(rel_new); } // ihit loop - streamlog_out( MESSAGE ) << " " << hits_to_save[icol].size() << " hits added to the collections: " - << m_outputTrackerHitsCollNames[icol] << ", " - << m_outputTrackerSimHitsCollNames[icol] << ", " - << m_outputTrackerHitRelNames[icol] << std::endl; + streamlog_out(MESSAGE) << " " << hits_to_save[icol].size() + << " hits added to the collections: " << m_outputTrackerHitsCollNames[icol] << ", " + << m_outputTrackerSimHitsCollNames[icol] << ", " << m_outputTrackerHitRelNames[icol] + << std::endl; - evt->addCollection( outputTrackerHitColls[icol], m_outputTrackerHitsCollNames[icol] ) ; - evt->addCollection( outputTrackerSimHitColls[icol], m_outputTrackerSimHitsCollNames[icol] ) ; - evt->addCollection( outputTrackerHitRels[icol], m_outputTrackerHitRelNames[icol] ) ; + evt->addCollection(outputTrackerHitColls[icol], m_outputTrackerHitsCollNames[icol]); + evt->addCollection(outputTrackerSimHitColls[icol], m_outputTrackerSimHitsCollNames[icol]); + evt->addCollection(outputTrackerHitRels[icol], m_outputTrackerHitRelNames[icol]); - streamlog_out( DEBUG5 ) << " output collection " << m_outputTrackerHitsCollNames[icol] << " of type " - << outputTrackerHitColls[icol]->getTypeName() << " added to the event \n" - << " output collection " << m_outputTrackerSimHitsCollNames[icol] << " of type " - << outputTrackerSimHitColls[icol]->getTypeName() << " added to the event \n" - << " output collection " << m_outputTrackerHitRelNames[icol] << " of type " - << outputTrackerHitRels[icol]->getTypeName() << " added to the event " - << std::endl ; + streamlog_out(DEBUG5) << " output collection " << m_outputTrackerHitsCollNames[icol] << " of type " + << outputTrackerHitColls[icol]->getTypeName() << " added to the event \n" + << " output collection " << m_outputTrackerSimHitsCollNames[icol] << " of type " + << outputTrackerSimHitColls[icol]->getTypeName() << " added to the event \n" + << " output collection " << m_outputTrackerHitRelNames[icol] << " of type " + << outputTrackerHitRels[icol]->getTypeName() << " added to the event " << std::endl; } // icol loop - streamlog_out(DEBUG) << " processing event: " << evt->getEventNumber() - << " in run: " << evt->getRunNumber() << std::endl ; + streamlog_out(DEBUG) << " processing event: " << evt->getEventNumber() << " in run: " << evt->getRunNumber() + << std::endl; - _nEvt ++ ; - -} - - - -void SplitCollectionByPolarAngle::check( LCEvent * ) { + _nEvt++; } +void SplitCollectionByPolarAngle::check(LCEvent*) {} - -void SplitCollectionByPolarAngle::end(){ - - std::cout << "SplitCollectionByPolarAngle::end() " << name() - << " processed " << _nEvt << " events in " << _nRun << " runs " - << std::endl ; - +void SplitCollectionByPolarAngle::end() { + std::cout << "SplitCollectionByPolarAngle::end() " << name() << " processed " << _nEvt << " events in " << _nRun + << " runs " << std::endl; }