diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/CMakeLists.txt b/InnerDetector/InDetRecTools/SiClusterizationTool/CMakeLists.txt index f50732775e67871e0a964aa35e0d74ca0ec8890d..3302f8da1c3a3f66ff74b6c0ab44fd721eef3a47 100644 --- a/InnerDetector/InDetRecTools/SiClusterizationTool/CMakeLists.txt +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/CMakeLists.txt @@ -13,6 +13,7 @@ atlas_depends_on_subdirs( PUBLIC DetectorDescription/Identifier Event/EventPrimitives GaudiKernel + InnerDetector/InDetRawEvent/InDetSimData InnerDetector/InDetConditions/InDetBeamSpotService InnerDetector/InDetConditions/InDetConditionsSummaryService InnerDetector/InDetDetDescr/InDetIdentifier @@ -44,13 +45,13 @@ atlas_add_library( SiClusterizationToolLib INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${EIGEN_INCLUDE_DIRS} PRIVATE_INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} PRIVATE_DEFINITIONS ${CLHEP_DEFINITIONS} - LINK_LIBRARIES ${ROOT_LIBRARIES} ${EIGEN_LIBRARIES} AthenaBaseComps AthenaKernel GeoPrimitives Identifier EventPrimitives GaudiKernel InDetIdentifier InDetReadoutGeometry InDetRawData InDetPrepRawData InDetRecToolInterfaces TrkParameters TrkNeuralNetworkUtilsLib + LINK_LIBRARIES ${ROOT_LIBRARIES} ${EIGEN_LIBRARIES} AthenaBaseComps AthenaKernel GeoPrimitives Identifier EventPrimitives GaudiKernel InDetSimData InDetIdentifier InDetReadoutGeometry InDetRawData InDetPrepRawData InDetRecToolInterfaces TrkParameters TrkNeuralNetworkUtilsLib PRIVATE_LINK_LIBRARIES ${CLHEP_LIBRARIES} AthenaPoolUtilities AtlasDetDescr TrkSurfaces TrkEventPrimitives VxVertex ) atlas_add_component( SiClusterizationTool src/components/*.cxx INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${EIGEN_INCLUDE_DIRS} - LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} ${EIGEN_LIBRARIES} AthenaBaseComps AthenaKernel GeoPrimitives Identifier EventPrimitives GaudiKernel InDetIdentifier InDetReadoutGeometry InDetRawData InDetPrepRawData InDetRecToolInterfaces TrkParameters AthenaPoolUtilities AtlasDetDescr TrkSurfaces TrkEventPrimitives VxVertex TrkNeuralNetworkUtilsLib SiClusterizationToolLib ) + LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} ${EIGEN_LIBRARIES} AthenaBaseComps AthenaKernel GeoPrimitives Identifier EventPrimitives GaudiKernel InDetSimData InDetIdentifier InDetReadoutGeometry InDetRawData InDetPrepRawData InDetRecToolInterfaces TrkParameters AthenaPoolUtilities AtlasDetDescr TrkSurfaces TrkEventPrimitives VxVertex TrkNeuralNetworkUtilsLib SiClusterizationToolLib ) # Install files from the package: atlas_install_joboptions( share/*.py ) diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/MergedPixelsTool.h b/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/MergedPixelsTool.h index e65c1586ea479286d20713d8b0267b77f300231b..89a30308d80a335030acf07aff3832f294d5d428 100755 --- a/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/MergedPixelsTool.h +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/MergedPixelsTool.h @@ -24,7 +24,6 @@ #include <vector> #include "GaudiKernel/ServiceHandle.h" -#include "GaudiKernel/IIncidentListener.h" // forward declare not possible (typedef) #include "InDetPrepRawData/PixelClusterCollection.h" @@ -47,7 +46,7 @@ namespace InDet { class IPixelClusterSplitProbTool; - class MergedPixelsTool : public PixelClusteringToolBase, virtual public IIncidentListener + class MergedPixelsTool : public PixelClusteringToolBase { public: @@ -96,8 +95,6 @@ namespace InDet { ///Statistics output StatusCode finalize(); - /** handle for incident service */ - void handle(const Incident& inc) ; private: typedef std::vector<Identifier> RDO_Vector; @@ -132,8 +129,33 @@ namespace InDet { InDetDD::SiDetectorElement* element, const PixelID& pixelID) const; + // Functions for merging broken cluster segments (to be used for upgrade studies) + // ITk: This function checks if two barrel clusters are potentially a single cluster which is broken into two pieces + bool mergeTwoBrokenClusters(const std::vector<Identifier>& group1, + const std::vector<Identifier>& group2, + InDetDD::SiDetectorElement* element, + const PixelID& pixelID) const; + // ITk: this function checks if two to-be-merged barrel proto-clusters have sizeZ consistent with cluster positions + bool mergeTwoClusters(const std::vector<Identifier>& group1, + const std::vector<Identifier>& group2, + InDetDD::SiDetectorElement* element, + const PixelID& pixelID) const; + + // ITk: checkSizeZ compares cluster sizeZ with expected cluster size for this cluster position (+/-200 mm for beam spread) + // checkSizeZ()=-1 if cluster is too small + // checkSizeZ()=0 if cluster sizeZ is within allowed range + // checkSizeZ()=1 if cluster is too large + // in the future, it may be changed to return deltaSizeZ + int checkSizeZ(int colmin, int colmax, int row, InDetDD::SiDetectorElement* element) const; + // this function returns expected sizeZ + int expectedSizeZ(int colmin, int colmax, int row, InDetDD::SiDetectorElement* element) const; + // this function returns size of the maximum gap between two cluster fragments + int maxGap(int colmin, int colmax, int row, InDetDD::SiDetectorElement* element) const; + //------- end of declaration of new functions + + ServiceHandle<IIncidentSvc> m_incidentSvc; //!< IncidentSvc to catch begin of event and end of envent - ServiceHandle<IBLParameterSvc> m_IBLParameterSvc; + ServiceHandle<IBLParameterSvc> m_IBLParameterSvc; bool m_emulateSplitter; //!< don't split - only emulate the split unsigned int m_minSplitSize; //!< minimum split size, regulates also the cluster splitting unsigned int m_maxSplitSize; //!< minimum split size, regulates also the cluster splitting @@ -143,15 +165,22 @@ namespace InDet { bool m_doIBLSplitting; bool m_IBLAbsent; - mutable InDet::PixelGangedClusterAmbiguities* m_splitClusterMap; //!< the actual split map - std::string m_splitClusterMapName; //!< split cluster ambiguity map + bool m_doMergeBrokenClusters; // ITk: switch to turn ON/OFF merging of broken clusters + bool m_doRemoveClustersWithToTequalSize; // ITk: switch to remove clusters with ToT=size + bool m_doCheckSizeBeforeMerging; // ITk: switch to check size of to-be-merged clusters + double m_beam_spread; // ITK: size of the luminous region, need to check cluster size + double m_lossThreshold; // maximum probability to loose N_mis consequitive pixels in a cluster + double m_pixelEff; // pixel efficiency (it depends on cluster eta; use smaller pixel efficiency) + + std::string m_splitClusterMapName;//No longer used Remove later //!< split cluster ambiguity map mutable unsigned int m_processedClusters; //!< statistics output mutable unsigned int m_modifiedOrigClusters; //!< statistics output mutable unsigned int m_splitOrigClusters; //!< statistics output mutable unsigned int m_splitProdClusters; //!< statistics output mutable unsigned int m_largeClusters; //!< statistics output - mutable int m_overflowIBLToT; + mutable int m_overflowIBLToT; + std::vector<int> m_minToT; /**< ToT cut */ ServiceHandle<IPixelOfflineCalibSvc> m_pixofflinecalibSvc; //ServiceHandle< StoreGateSvc > m_detStore; //const PixelID* m_idHelper; diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/NnClusterizationFactory.h b/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/NnClusterizationFactory.h index 9be495d60b8823216521317d5480c30ddea4f636..f27114f9939748a4e7d6ecdae8f89f8818ca1e37 100644 --- a/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/NnClusterizationFactory.h +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/NnClusterizationFactory.h @@ -215,6 +215,7 @@ namespace InDet { bool m_useRecenteringNNWithTracks; double m_correctLorShiftBarrelWithoutTracks; double m_correctLorShiftBarrelWithTracks; + }; }//end InDet namespace diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/SCT_ClusteringTool.h b/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/SCT_ClusteringTool.h index ba771504b193e6d47e35ab5cd1fae85e79766ff7..adf7d8b07f3965841794f83ad260309cd96919e4 100755 --- a/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/SCT_ClusteringTool.h +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/SCT_ClusteringTool.h @@ -28,6 +28,10 @@ class SCT_ChannelStatusAlg; class StatusCode; class IInDetConditionsSvc; +namespace InDetDD{ + class SCT_ModuleSideDesign; +} + namespace InDet { /** AlgTool for SCT_Clusterization. * Input is from RDOs, assumed to be sorted. They are then scanned @@ -64,7 +68,11 @@ namespace InDet { ToolHandle< ClusterMakerTool > m_clusterMaker; typedef std::vector<Identifier> IdVec_t; std::string m_timeBinStr; - int m_timeBinBits[3]; + int m_timeBinBits[3]; + bool m_innermostBarrelX1X; + bool m_innertwoBarrelX1X; + bool m_majority01X; + bool m_useRowInformation; ///Add strips to a cluster vector without checking for bad strips void addStripsToCluster(const Identifier & firstStripId, const unsigned int nStrips,IdVec_t & clusterVector,const SCT_ID& idHelper) const; @@ -73,6 +81,8 @@ namespace InDet { void addStripsToClusterWithChecks(const Identifier & firstStripId, const unsigned int nStrips, IdVec_t & clusterVector, std::vector<IdVec_t > & idGroups, const SCT_ID& idHelper) const; + void addStripsToClusterInclRows(const Identifier & firstStripId, const unsigned int nStrips, IdVec_t & clusterVector,std::vector<IdVec_t > & idGroups, const SCT_ID& idHelper) const; + /** Recluster the current vector, splitting on bad strips, and insert those new groups to the idGroups vector. * The cluster vector referenced will be changed by this, as well as the idGroups **/ @@ -87,7 +97,10 @@ namespace InDet { ///Calculate the cluster position and width given the first and last strip numbers for this element DimensionAndPosition clusterDimensions(const int firstStrip, const int lastStrip, const InDetDD::SiDetectorElement* element, - const SCT_ID& idHelper) const; + const SCT_ID& idHelper) const; + + DimensionAndPosition clusterDimensionsInclRow(const int firstStrip, const int lastStrip, const int row, const InDetDD::SiDetectorElement* element, const InDetDD::SCT_ModuleSideDesign* design) const; + /// In-class facade on the 'isGood' method for a strip identifier bool isBad(const Identifier & stripId) const; @@ -97,6 +110,8 @@ namespace InDet { StatusCode decodeTimeBin(const char& timeBin, int& bit); // Test the clusters time bin to see if matches pattern bool testTimeBins(int timeBin) const; + bool testTimeBins01X(int timeBin) const; + bool testTimeBinsX1X(int timeBin) const; };//end of class }//end of namespace diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/TruthClusterizationFactory.h b/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/TruthClusterizationFactory.h new file mode 100644 index 0000000000000000000000000000000000000000..f69758d9c7c545a66283cf1f768d5a698cba2a37 --- /dev/null +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/TruthClusterizationFactory.h @@ -0,0 +1,74 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + + #ifndef SICLUSTERIZATIONTOOL_TruthClusterizationFactory_C + #define SICLUSTERIZATIONTOOL_TruthClusterizationFactory_C + + /****************************************************** + @class TruthClusterizationFactory + @author Roland Jansky + Package : SiClusterizationTool + Created : April 2016 + DESCRIPTION: Emulates NN evaluation from truth (for ITK studies) + ********************************************************/ + + + #include "AthenaBaseComps/AthAlgTool.h" + #include "GaudiKernel/ToolHandle.h" + #include "GaudiKernel/IIncidentSvc.h" + + #include <vector> + #include <string> + #include <map> + + #include <TString.h> + #include "AthenaKernel/IOVSvcDefs.h" + + +//this is a typedef: no forward decl possible +#include "TrkParameters/TrackParameters.h" +#include "GeoPrimitives/GeoPrimitives.h" +#include "EventPrimitives/EventPrimitives.h" + + class InDetSimDataCollection; + +namespace InDet { + + class PixelCluster; + + static const InterfaceID IID_TruthClusterizationFactory("InDet::NnClusterizationFactory", 1, 0); + + class TruthClusterizationFactory : virtual public IIncidentListener, + public AthAlgTool { + + public: + + /** AlgTool interface methods */ + static const InterfaceID& interfaceID() { return IID_TruthClusterizationFactory; }; + + TruthClusterizationFactory(const std::string& name, + const std::string& n, const IInterface* p); + ~TruthClusterizationFactory(); + + virtual StatusCode initialize(); + virtual StatusCode finalize() { return StatusCode::SUCCESS; }; + + /** handle for incident service */ + virtual void handle(const Incident& inc); + + std::vector<double> estimateNumberOfParticles(const InDet::PixelCluster& pCluster); + + std::vector<Amg::Vector2D> estimatePositions(const InDet::PixelCluster&); + + private: + /** IncidentSvc to catch begining of event and end of event */ + ServiceHandle<IIncidentSvc> m_incidentSvc; + + std::string m_simDataCollectionName; //!< sim data collection name + mutable const InDetSimDataCollection* m_simDataCollection; //!< sim data collection - refreshed at BeginEvent incident + }; + + }//end InDet namespace + + #endif diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/TruthPixelClusterSplitProbTool.h b/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/TruthPixelClusterSplitProbTool.h new file mode 100644 index 0000000000000000000000000000000000000000..1080632e2926f9679292352fb3bbd6b1b2a0a333 --- /dev/null +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/TruthPixelClusterSplitProbTool.h @@ -0,0 +1,67 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +/////////////////////////////////////////////////////////////////// +//TruthPixelClusterSplitProbTool.h +/////////////////////////////////////////////////////////////////// +// (c) ATLAS Detector software +/////////////////////////////////////////////////////////////////// +// Estimate cluster split probability by using multivariate techniquew +/////////////////////////////////////////////////////////////////// +// +///////////////////////////////////////////////////////////////////////// +// @author Roland Jansky & Felix Cormier +// +/////////////////////////////////////////////////////////////////// + +#ifndef SiClusterizationTool_TruthPixelClusterSplitProbTool_H +#define SiClusterizationTool_TruthPixelClusterSplitProbTool_H + + +#include "GaudiKernel/ToolHandle.h" +#include "AthenaBaseComps/AthAlgTool.h" +#include "InDetRecToolInterfaces/IPixelClusterSplitProbTool.h" +//Beam Spot Condition +#include "InDetBeamSpotService/IBeamCondSvc.h" +#include "TrkParameters/TrackParameters.h" + +class IBeamCondSvc; + +namespace InDet { + + class TruthClusterizationFactory; + + class TruthPixelClusterSplitProbTool : public AthAlgTool, virtual public IPixelClusterSplitProbTool + { + public: + + TruthPixelClusterSplitProbTool(const std::string& t, const std::string& n, const IInterface* p); + + virtual ~TruthPixelClusterSplitProbTool() {}; + + StatusCode initialize(); + + virtual InDet::PixelClusterSplitProb splitProbability(const InDet::PixelCluster& origCluster ) const; + + virtual InDet::PixelClusterSplitProb splitProbability(const InDet::PixelCluster& origCluster, const Trk::TrackParameters& trackParameters ) const; + + + private: + + InDet::PixelClusterSplitProb compileSplitProbability(std::vector<double>& vectorOfProbs ) const; + + ToolHandle<TruthClusterizationFactory> m_truthClusterizationFactory; + ServiceHandle<IBeamCondSvc> m_iBeamCondSvc; + + std::vector<double> m_priorMultiplicityContent; + + bool m_useBeamSpotInfo; + + }; + +} + + +#endif + diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/TruthPixelClusterSplitter.h b/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/TruthPixelClusterSplitter.h new file mode 100644 index 0000000000000000000000000000000000000000..876163571280a3c441080c9cf9e56cbffdcb5762 --- /dev/null +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/SiClusterizationTool/TruthPixelClusterSplitter.h @@ -0,0 +1,65 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +/////////////////////////////////////////////////////////////////// +// TruthPixelClusterSplitter.h, (c) ATLAS Detector software +/////////////////////////////////////////////////////////////////// + +#ifndef INDETRECTOOLS_TruthPixelClusterSplitter_H +#define INDETRECTOOLS_TruthPixelClusterSplitter_H + + +#include "AthenaBaseComps/AthAlgTool.h" +#include "GaudiKernel/ToolHandle.h" + +#include "InDetRecToolInterfaces/IPixelClusterSplitter.h" +#include "InDetPrepRawData/PixelClusterParts.h" +#include "InDetPrepRawData/PixelClusterSplitProb.h" + +class IBeamCondSvc; + +namespace InDet +{ + + class TruthClusterizationFactory; + class PixelCluster; + + /** @class TruthPixelClusterSplitter + @author Roland.Jansky@cern.ch + */ + class TruthPixelClusterSplitter : public AthAlgTool, virtual public IPixelClusterSplitter { + public : + /** Constructor*/ + TruthPixelClusterSplitter(const std::string &type, + const std::string &name, + const IInterface *parent); + + /** Destructor*/ + ~TruthPixelClusterSplitter(); + + /** AthAlgTool interface methods */ + StatusCode initialize(); + StatusCode finalize(); + + /** take one, give zero or many */ + std::vector<InDet::PixelClusterParts> splitCluster(const InDet::PixelCluster& origCluster ) const; + + /** take one, give zero or many - with split probability object */ + std::vector<InDet::PixelClusterParts> splitCluster(const InDet::PixelCluster& origCluster, + const InDet::PixelClusterSplitProb& spo) const; + + private: + + ToolHandle<TruthClusterizationFactory> m_truthClusterizationFactory; + ServiceHandle<IBeamCondSvc> m_iBeamCondSvc; + + double m_thresholdSplittingIntoTwoClusters; + double m_thresholdSplittingIntoThreeClusters; + bool m_splitOnlyOnBLayer; + + bool m_useBeamSpotInfo; + + }; +} +#endif diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/cmt/requirements b/InnerDetector/InDetRecTools/SiClusterizationTool/cmt/requirements index c1182e4a0d61cde660564906858865ee52f18a0a..03f6fdc923d150a6c75eb4a31980faa7f63a7d03 100755 --- a/InnerDetector/InDetRecTools/SiClusterizationTool/cmt/requirements +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/cmt/requirements @@ -22,8 +22,8 @@ use InDetReadoutGeometry InDetReadoutGeometry-* InnerDetector/InDetDetDescr use InDetConditionsSummaryService InDetConditionsSummaryService-* InnerDetector/InDetConditions use TrkParameters TrkParameters-* Tracking/TrkEvent -use GeoPrimitives GeoPrimitives-* DetectorDescription -use EventPrimitives EventPrimitives-* Event +use GeoPrimitives GeoPrimitives-* DetectorDescription +use EventPrimitives EventPrimitives-* Event private use TrkEventPrimitives TrkEventPrimitives-* Tracking/TrkEvent @@ -39,6 +39,8 @@ use TrkNeuralNetworkUtils TrkNeuralNetworkUtils-* Tracking use VxVertex VxVertex-* Tracking/TrkEvent use TrkSurfaces TrkSurfaces-* Tracking/TrkDetDescr use PixelGeoModel PixelGeoModel-* InnerDetector/InDetDetDescr +use InDetSimData InDetSimData-* InnerDetector/InDetRawEvent +use AtlasHepMC AtlasHepMC-* External public diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/src/ClusterMakerTool.cxx b/InnerDetector/InDetRecTools/SiClusterizationTool/src/ClusterMakerTool.cxx index 5fe0882002d538943e990d94416dc84977b89420..de90aad3efa9518d7a8a3c08e67d73bacd891b38 100755 --- a/InnerDetector/InDetRecTools/SiClusterizationTool/src/ClusterMakerTool.cxx +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/src/ClusterMakerTool.cxx @@ -67,18 +67,7 @@ StatusCode ClusterMakerTool::initialize(){ // Code entered here will be executed once at program start. ATH_MSG_INFO ( name() << " initialize()" ); - -/* - if (m_IBLParameterSvc.retrieve().isFailure()) { - ATH_MSG_WARNING( "Could not retrieve IBLParameterSvc"); - } - else { - m_offlineCalibSvc.setTypeAndName(m_IBLParameterSvc->setStringParameters(m_offlineCalibSvc.typeAndName(),"PixelOfflineCalibSvc")); - m_calibSvc.setTypeAndName(m_IBLParameterSvc->setStringParameters(m_calibSvc.typeAndName(),"PixelCalibSvc")); - m_IBLParameterSvc->setBoolParameters(m_calibrateCharge,"UsePixelCalibCondDB"); - } -*/ - + // Protect from the situation in which the PixelOfflineCalibSvc is not // configured: that should be the case if no PixelRDO are read in. // AA 01/10/2009 @@ -195,7 +184,9 @@ PixelCluster* ClusterMakerTool::pixelCluster( for (int i=0; i<nRDO; i++) { Identifier pixid=rdoList[i]; int ToT=totList[i]; + float charge = m_calibSvc->getCharge(pixid,ToT); + chargeList.push_back(charge); } } @@ -313,7 +304,7 @@ PixelCluster* ClusterMakerTool::pixelCluster( double splitProb2) const{ - if (msgLvl(MSG::VERBOSE)) msg() << "ClusterMakerTool called, number " << endreq; + if (msgLvl(MSG::VERBOSE)) msg() << "ClusterMakerTool called, number " << endmsg; // Add protection in case m_offlineCalibSvc is not configured // but errorStrategy==2 is requested. // That should never happen, since the m_offlineCalibSvc should be switched @@ -350,13 +341,14 @@ PixelCluster* ClusterMakerTool::pixelCluster( for (int i=0; i<nRDO; i++) { Identifier pixid=rdoList[i]; int ToT=totList[i]; - + float charge = ToT; if (m_calibrateCharge){ + charge = m_calibSvc->getCharge(pixid,ToT); + chargeList.push_back(charge); } - // std::cout << "tot, charge = " << ToT << " " << charge << std::endl; int row = pixelID.phi_index(pixid); int col = pixelID.eta_index(pixid); @@ -391,7 +383,7 @@ PixelCluster* ClusterMakerTool::pixelCluster( if(qRowMin+qRowMax > 0) omegax = qRowMax/float(qRowMin+qRowMax); if(qColMin+qColMax > 0) omegay = qColMax/float(qColMin+qColMax); - if (msgLvl(MSG::VERBOSE)) msg() << "omega = " << omegax << " " << omegay << endreq; + if (msgLvl(MSG::VERBOSE)) msg() << "omega = " << omegax << " " << omegay << endmsg; // ask for Lorentz correction, get global position double shift = element->getLorentzCorrection(); @@ -558,7 +550,7 @@ SCT_Cluster* ClusterMakerTool::sctCluster( } // rotation for endcap SCT - if(element->design().shape() == InDetDD::Trapezoid) { + if(element->design().shape() == InDetDD::Trapezoid || element->design().shape() == InDetDD::Annulus) { double sn = element->sinStereoLocal(localPos); double sn2 = sn*sn; double cs2 = 1.-sn2; @@ -603,7 +595,7 @@ double ClusterMakerTool::getPixelCTBPhiError(int layer, int phi, // shouldn't really happen... if(msgLvl(MSG::WARNING)) msg() << "Unexpected layer and phi numbers: layer = " - << layer << " and phi = " << phi << endreq; + << layer << " and phi = " << phi << endmsg; return 14.6*micrometer; } diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/src/MergedPixelsTool.cxx b/InnerDetector/InDetRecTools/SiClusterizationTool/src/MergedPixelsTool.cxx index f34c8635c040575d5c21ff336094e6824b166072..b888b8a1a7523ed1732e46957a2cedd4f0c81e6a 100755 --- a/InnerDetector/InDetRecTools/SiClusterizationTool/src/MergedPixelsTool.cxx +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/src/MergedPixelsTool.cxx @@ -8,6 +8,8 @@ // input from RDOs. // // Changes +// 18/03/2016 Sasha Pranko +// Added methods to merge segments of broken clusters: need for ITK studies // 29/1/2011 Andi Salzburger & Giacinto Piacquadio (Pixel Clusterization TF) // Added support for splitting & use of more refined clustering methods /////////////////////////////////////////////////////////////////// @@ -58,11 +60,16 @@ namespace InDet { m_minSplitSize(1), m_maxSplitSize(1000), m_minSplitProbability(0.), - m_splitProbTool(""), - m_clusterSplitter(""), + m_splitProbTool("", this), + m_clusterSplitter("", this), m_doIBLSplitting(false), m_IBLAbsent(true), - m_splitClusterMap(0), + m_doMergeBrokenClusters(false), /* ITk: switch to turn ON/OFF merging of broken clusters */ + m_doRemoveClustersWithToTequalSize(false), /* ITk: switch to remove clusters with ToT=size */ + m_doCheckSizeBeforeMerging(false), /* ITk: switch to check size to-be-merged clusters */ + m_beam_spread(200.0), /* ITk: size of luminous region, needed for cluster size check */ + m_lossThreshold(0.001), /* ITk: maximum probability to loose N_mis consequitive pixels in a cluster */ + m_pixelEff(0.90), /* ITk: pixel efficiency (it depends on cluster eta; use smaller pixel efficiency) */ m_splitClusterMapName("SplitClusterAmbiguityMap"), m_processedClusters(0), m_modifiedOrigClusters(0), @@ -70,6 +77,7 @@ namespace InDet { m_splitProdClusters(0), m_largeClusters(0), m_overflowIBLToT(0), + m_minToT({0,0,0,0,0,0,0}), m_pixofflinecalibSvc("PixelOfflineCalibSvc", name) //m_detStore("DetectorStore", name), //m_idHelper(0) @@ -83,7 +91,14 @@ namespace InDet { declareProperty("SplitProbTool", m_splitProbTool); declareProperty("ClusterSplitter", m_clusterSplitter); declareProperty("DoIBLSplitting", m_doIBLSplitting); - declareProperty("SplitClusterAmbiguityMap", m_splitClusterMapName); + declareProperty("SplitClusterAmbiguityMap", m_splitClusterMapName); //No longer used Remove later + declareProperty("DoMergeBrokenClusters", m_doMergeBrokenClusters); // ITk: switch to turn ON/OFF merging of broken clusters + declareProperty("DoRemoveClustersWithToTequalSize",m_doRemoveClustersWithToTequalSize); // ITk: switch to remove clusters with ToT=size + declareProperty("DoCheckSizeBeforeMerging", m_doCheckSizeBeforeMerging); // ITk: switch to check size to-be-merged clusters + declareProperty("BeamSpread", m_beam_spread); // ITk: size of luminous region, needed for cluster size check + declareProperty("LossProbability", m_lossThreshold); // ITk: maximum probability to loose N_mis consequitive pixels in a cluster + declareProperty("MinPixelEfficiency", m_pixelEff); // ITk: pixel efficiency (it depends on cluster eta; use smaller pixel efficiency) + declareProperty("ToTMinCut", m_minToT, "Minimum ToT cut [IBL, b-layer, L1, L2, Endcap, DBM, ITk extra"); } //--------------------------------------------------------------------------- @@ -109,12 +124,12 @@ namespace InDet { //return StatusCode::FAILURE; //} //else { - // msg(MSG::DEBUG) << "DetectorStore service found" << endreq; + // msg(MSG::DEBUG) << "DetectorStore service found" << endmsg; //} //// Get the PixelID Helper //if (m_detStore->retrieve(m_idHelper, "PixelID").isFailure()) { - // msg(MSG::FATAL) << "Could not get Pixel ID helper" << endreq; + // msg(MSG::FATAL) << "Could not get Pixel ID helper" << endmsg; //return StatusCode::FAILURE; //} @@ -134,14 +149,19 @@ namespace InDet { } // register to the incident service: // (a) BeginEvent needed for the recording of the Split ambiguity map (always need to be there for merging) - m_incidentSvc->addListener( this, IncidentType::BeginEvent); +// m_incidentSvc->addListener( this, IncidentType::BeginEvent); // (b) EndEvent needed for memory cleanup - m_incidentSvc->addListener( this, IncidentType::EndEvent); +// m_incidentSvc->addListener( this, IncidentType::EndEvent); if (m_pixofflinecalibSvc.retrieve().isFailure()){ ATH_MSG_ERROR("Could not retrieve " << m_pixofflinecalibSvc); return StatusCode::FAILURE; } + + if (m_minToT.size() != 7){ + ATH_MSG_ERROR("Number of entries for ToT Cut is:" << m_minToT.size() << " . 7 Values are needed, so fix jO."); + return StatusCode::FAILURE; + } return PixelClusteringToolBase::initialize(); } @@ -166,205 +186,211 @@ namespace InDet { // module (with non-empty RDO collection...). // It clusters together the RDOs with a pixell cell side in common. - PixelClusterCollection* MergedPixelsTool::clusterize( - const InDetRawDataCollection<PixelRDORawData> &collection, - const InDetDD::SiDetectorManager& manager, - const PixelID& pixelID) const - { + PixelClusterCollection* MergedPixelsTool::clusterize( + const InDetRawDataCollection<PixelRDORawData> &collection, + const InDetDD::SiDetectorManager& manager, + const PixelID& pixelID) const + { // Get the messaging service, print where you are - - Identifier elementID = collection.identify(); - + + Identifier elementID = collection.identify(); + // Get hash Identifier for these RDOs - IdentifierHash idHash = collection.identifyHash(); - + IdentifierHash idHash = collection.identifyHash(); + // Size of RDO's collection: - unsigned int RDO_size = collection.size(); - if ( RDO_size==0) - { - // Empty RDO collection - ATH_MSG_DEBUG (" areNeighbours - problems "); - return 0; - } - - int clusterNumber = 0; - + unsigned int RDO_size = collection.size(); + if ( RDO_size==0) + { + // Empty RDO collection + ATH_MSG_DEBUG (" areNeighbours - problems "); + return 0; + } + + int clusterNumber = 0; + // If module is bad, do not create a cluster collection - if (m_useModuleMap && !(m_summarySvc->isGood(idHash))) - return 0; - - - typedef InDetRawDataCollection<PixelRDORawData> RDO_Collection; - - RDO_Collection::const_iterator firstRDO(collection.begin()); - RDO_Collection::const_iterator lastRDO(collection.end()); - - // Get detector info. + if (m_useModuleMap && !(m_summarySvc->isGood(idHash))) + return 0; + + + typedef InDetRawDataCollection<PixelRDORawData> RDO_Collection; + + RDO_Collection::const_iterator firstRDO(collection.begin()); + RDO_Collection::const_iterator lastRDO(collection.end()); + + // Get detector info. // Find detector element for these RDOs - - InDetDD::SiDetectorElement* element = manager.getDetectorElement(elementID); - - const Trk::RectangleBounds *mybounds=dynamic_cast<const Trk::RectangleBounds *>(&element->surface().bounds()); - if (not mybounds){ - ATH_MSG_ERROR("Dynamic cast failed at "<<__LINE__<<" of MergedPixelsTool.cxx."); - return nullptr; - } + + InDetDD::SiDetectorElement* element = manager.getDetectorElement(elementID); + + const Trk::RectangleBounds *mybounds=dynamic_cast<const Trk::RectangleBounds *>(&element->surface().bounds()); + if (not mybounds){ + ATH_MSG_ERROR("Dynamic cast failed at "<<__LINE__<<" of MergedPixelsTool.cxx."); + return nullptr; + } // RDOs will be accumulated as a vector of Identifier, which will // be added to a vector of groups. - - RDO_GroupVector rdoGroups; - rdoGroups.reserve(RDO_size); - + + RDO_GroupVector rdoGroups; + rdoGroups.reserve(RDO_size); + // A similar structure is created for TOT information - TOT_GroupVector totGroups; - totGroups.reserve(RDO_size); - TOT_GroupVector lvl1Groups; - lvl1Groups.reserve(RDO_size); - + TOT_GroupVector totGroups; + totGroups.reserve(RDO_size); + TOT_GroupVector lvl1Groups; + lvl1Groups.reserve(RDO_size); + // loop on all RDOs - for(RDO_Collection::const_iterator nextRDO=collection.begin() ; + for(RDO_Collection::const_iterator nextRDO=collection.begin() ; nextRDO!=collection.end() ; ++nextRDO) - { - Identifier rdoID= (*nextRDO)->identify(); - // if a pixel is not skip it in clusterization - if ( m_usePixelMap && !(m_summarySvc->isGood(idHash,rdoID)) ) continue; - int tot = (*nextRDO)->getToT(); - int lvl1= (*nextRDO)->getLVL1A(); - // check if this is a ganged pixel - Identifier gangedID; - // bool ganged= isGanged(rdoID, elementID, pixelID, gangedID); - bool ganged= isGanged(rdoID, element, gangedID); - if(ganged){ - ATH_MSG_VERBOSE("Ganged Pixel, row = " << pixelID.phi_index(rdoID) - << "Ganged row = " << pixelID.phi_index(gangedID)); - } - else - ATH_MSG_VERBOSE("Not ganged Pixel, row = " << pixelID.phi_index(rdoID)); - - // loop on all existing RDO groups, until you find that the RDO - // is a neighbour of the group. - bool found= false; - - RDO_GroupVector::iterator firstGroup = rdoGroups.begin(); - RDO_GroupVector::iterator lastGroup = rdoGroups.end(); - TOT_GroupVector::iterator totGroup = totGroups.begin(); - TOT_GroupVector::iterator lvl1Group = lvl1Groups.begin(); - - while( !found && firstGroup!= lastGroup) - { - // Check if RDO is neighbour of the cluster - if (areNeighbours(**firstGroup, rdoID, element, pixelID)) - { - - // if RDO is a duplicate of one in the cluster, do not add it. - // Instead the method isDuplicated check wether the new - // one has a larger LVL1 - if so it does replace the old - // lvl1 with the new one. - if(!isDuplicated(**firstGroup, **lvl1Group, rdoID, lvl1, pixelID)){ - (*firstGroup)->push_back(rdoID); - (*totGroup)->push_back(tot); - (*lvl1Group)->push_back(lvl1); + { + Identifier rdoID= (*nextRDO)->identify(); + // if a pixel is not skip it in clusterization + if ( m_usePixelMap && !(m_summarySvc->isGood(idHash,rdoID)) ) continue; + int tot = (*nextRDO)->getToT(); + int layerIndex = pixelID.layer_disk(rdoID); + if (layerIndex>=4) { layerIndex=6; } // ITk 5th layer + if (abs(pixelID.barrel_ec(rdoID))==2) { layerIndex=4; } // disks + if (abs(pixelID.barrel_ec(rdoID))==4) { layerIndex=5; } // DBM + // cut on minimum ToT + if (tot<m_minToT.at(layerIndex)) { continue; } // skip hits with ToT less than ToT cut + int lvl1= (*nextRDO)->getLVL1A(); + // check if this is a ganged pixel + Identifier gangedID; + // bool ganged= isGanged(rdoID, elementID, pixelID, gangedID); + bool ganged= isGanged(rdoID, element, gangedID); + if(ganged){ + ATH_MSG_VERBOSE("Ganged Pixel, row = " << pixelID.phi_index(rdoID) + << "Ganged row = " << pixelID.phi_index(gangedID)); + } + else + ATH_MSG_VERBOSE("Not ganged Pixel, row = " << pixelID.phi_index(rdoID)); + + // loop on all existing RDO groups, until you find that the RDO + // is a neighbour of the group. + bool found= false; + + RDO_GroupVector::iterator firstGroup = rdoGroups.begin(); + RDO_GroupVector::iterator lastGroup = rdoGroups.end(); + TOT_GroupVector::iterator totGroup = totGroups.begin(); + TOT_GroupVector::iterator lvl1Group = lvl1Groups.begin(); + + while( !found && firstGroup!= lastGroup) + { + // Check if RDO is neighbour of the cluster + if (areNeighbours(**firstGroup, rdoID, element, pixelID)) + { + + // if RDO is a duplicate of one in the cluster, do not add it. + // Instead the method isDuplicated check wether the new + // one has a larger LVL1 - if so it does replace the old + // lvl1 with the new one. + if(!isDuplicated(**firstGroup, **lvl1Group, rdoID, lvl1, pixelID)){ + (*firstGroup)->push_back(rdoID); + (*totGroup)->push_back(tot); + (*lvl1Group)->push_back(lvl1); // see if it is a neighbour to any other groups - checkForMerge(rdoID, firstGroup, lastGroup, - totGroup, lvl1Group, element, pixelID); - } - else{ - ATH_MSG_VERBOSE("duplicate found"); - } - found = true; - } - ++firstGroup; - ++totGroup; - ++lvl1Group; - } - - // if RDO is isolated, create new cluster. - - if(!found) - { - RDO_Vector* newGroup= new RDO_Vector; - rdoGroups.push_back(newGroup); - newGroup->push_back(rdoID); - TOT_Vector* newtotGroup = new TOT_Vector; - totGroups.push_back(newtotGroup); - newtotGroup->push_back(tot); - TOT_Vector* newlvl1Group = new TOT_Vector; - lvl1Groups.push_back(newlvl1Group); - newlvl1Group->push_back(lvl1); - } - - // Repeat for ganged pixel if necessary - if (! ganged) continue; - - ATH_MSG_VERBOSE("Ganged pixel, row = " << pixelID.phi_index(gangedID)); - found= false; - - firstGroup = rdoGroups.begin(); - lastGroup = rdoGroups.end(); - totGroup = totGroups.begin(); - lvl1Group = lvl1Groups.begin(); - - while( !found && firstGroup!= lastGroup) - { - // if neighbour of the cluster, add it to the cluster - - if ( areNeighbours(**firstGroup, gangedID, element, pixelID) ) - { - ATH_MSG_VERBOSE("Ganged pixel is neighbour of a cluster"); - - // if RDO is a duplicate of one in the cluster, do not add it. - // Instead the method isDuplicated check wether the new - // one has a larger LVL1 - if so it does replace the old - // lvl1 with the new one. - if(!isDuplicated(**firstGroup, **lvl1Group, gangedID, lvl1, pixelID)){ - - (*firstGroup)->push_back(gangedID); - (*totGroup)->push_back(tot); - (*lvl1Group)->push_back(lvl1); - checkForMerge(gangedID, firstGroup, lastGroup, - totGroup, lvl1Group, element, pixelID); - } - else{ - ATH_MSG_VERBOSE("duplicate found"); - } - found = true; - - } - ++firstGroup; - ++totGroup; - ++lvl1Group; - } - - // if isolated, create new cluster. - - if(!found) - { - ATH_MSG_VERBOSE("New cluster with ganged pixel"); - RDO_Vector* newGroup= new RDO_Vector; - rdoGroups.push_back(newGroup); - newGroup->push_back(gangedID); - TOT_Vector* newtotGroup = new TOT_Vector; - totGroups.push_back(newtotGroup); - newtotGroup->push_back(tot); - TOT_Vector* newlvl1Group = new TOT_Vector; - lvl1Groups.push_back(newlvl1Group); - newlvl1Group->push_back(lvl1); - } - - } - - if(totGroups.size() != rdoGroups.size()) - ATH_MSG_ERROR("Mismatch between RDO identifier and TOT info!"); + checkForMerge(rdoID, firstGroup, lastGroup, + totGroup, lvl1Group, element, pixelID); + } + else{ + ATH_MSG_VERBOSE("duplicate found"); + } + found = true; + } + ++firstGroup; + ++totGroup; + ++lvl1Group; + } + + // if RDO is isolated, create new cluster. + + if(!found) + { + RDO_Vector* newGroup= new RDO_Vector; + rdoGroups.push_back(newGroup); + newGroup->push_back(rdoID); + TOT_Vector* newtotGroup = new TOT_Vector; + totGroups.push_back(newtotGroup); + newtotGroup->push_back(tot); + TOT_Vector* newlvl1Group = new TOT_Vector; + lvl1Groups.push_back(newlvl1Group); + newlvl1Group->push_back(lvl1); + } + + // Repeat for ganged pixel if necessary + if (! ganged) continue; + + ATH_MSG_VERBOSE("Ganged pixel, row = " << pixelID.phi_index(gangedID)); + found= false; + + firstGroup = rdoGroups.begin(); + lastGroup = rdoGroups.end(); + totGroup = totGroups.begin(); + lvl1Group = lvl1Groups.begin(); + + while( !found && firstGroup!= lastGroup) + { + // if neighbour of the cluster, add it to the cluster + + if ( areNeighbours(**firstGroup, gangedID, element, pixelID) ) + { + ATH_MSG_VERBOSE("Ganged pixel is neighbour of a cluster"); + + // if RDO is a duplicate of one in the cluster, do not add it. + // Instead the method isDuplicated check wether the new + // one has a larger LVL1 - if so it does replace the old + // lvl1 with the new one. + if(!isDuplicated(**firstGroup, **lvl1Group, gangedID, lvl1, pixelID)){ + + (*firstGroup)->push_back(gangedID); + (*totGroup)->push_back(tot); + (*lvl1Group)->push_back(lvl1); + checkForMerge(gangedID, firstGroup, lastGroup, + totGroup, lvl1Group, element, pixelID); + } + else{ + ATH_MSG_VERBOSE("duplicate found"); + } + found = true; + + } + ++firstGroup; + ++totGroup; + ++lvl1Group; + } + + // if isolated, create new cluster. + + if(!found) + { + ATH_MSG_VERBOSE("New cluster with ganged pixel"); + RDO_Vector* newGroup= new RDO_Vector; + rdoGroups.push_back(newGroup); + newGroup->push_back(gangedID); + TOT_Vector* newtotGroup = new TOT_Vector; + totGroups.push_back(newtotGroup); + newtotGroup->push_back(tot); + TOT_Vector* newlvl1Group = new TOT_Vector; + lvl1Groups.push_back(newlvl1Group); + newlvl1Group->push_back(lvl1); + } + + } + + if(totGroups.size() != rdoGroups.size()) + ATH_MSG_ERROR("Mismatch between RDO identifier and TOT info!"); // We now have groups of contiguous RDOs. Make clusters // Make a new empty cluster collection - PixelClusterCollection *clusterCollection = - new PixelClusterCollection(idHash); - clusterCollection->setIdentifier(elementID); - clusterCollection->reserve(rdoGroups.size()); - - std::vector<TOT_Vector *>::iterator totgroup = totGroups.begin(); - std::vector<TOT_Vector *>::iterator lvl1group= lvl1Groups.begin(); - + PixelClusterCollection *clusterCollection = + new PixelClusterCollection(idHash); + clusterCollection->setIdentifier(elementID); + clusterCollection->reserve(rdoGroups.size()); + + std::vector<TOT_Vector *>::iterator totgroup = totGroups.begin(); + std::vector<TOT_Vector *>::iterator lvl1group= lvl1Groups.begin(); + // LOOP over the RDO-groups to be clustered ---------------------------------------------------------------- // Logics in splitting is: the splitter is only activated clusters within min/max split size, // single pixel clusters are not split into more clusters, but can be modified if included by split size req. @@ -378,267 +404,274 @@ namespace InDet { // // (MAIN CLUSTERIZATION LOOP after connected component finding) - for(std::vector<RDO_Vector *>::iterator group = rdoGroups.begin() ; group!= rdoGroups.end() ; ++group) - { - - // the split probabilities - // writing the split boolean is done the following way - // - if in emulation mode: always write the output of the splitter for validation - // - if in pseudo-emulation mode (1-pixel clusters): set boolean to false, but keep split probs - bool clusterModified = false; - bool clusterSplit = false; - double clusterSplitP1 = 0.; - double clusterSplitP2 = 0.; - bool singlePixelSplitCase = 0; - // the size of the inition rdo group - size_t groupSize = (**group).size(); - // If cluster is empty, i.e. it has been merged with another, - // do not attempt to make cluster. - if ( groupSize > 0 ) - { - ++m_processedClusters; - // create the original cluster - split & split probs are by default false and 0 - PixelCluster* cluster = makeCluster(**group, - **totgroup, - **lvl1group, - element, - pixelID, - ++clusterNumber); - // check for splitting - if ( groupSize >= m_minSplitSize && groupSize <= m_maxSplitSize ) { - - // prepare for the return value of the pixel cluster - std::vector<InDet::PixelClusterParts> splitClusterParts; - if ( !m_splitProbTool.empty() && (m_doIBLSplitting || m_IBLAbsent || !element->isBlayer())){ - InDet::PixelClusterSplitProb splitProbObj = m_splitProbTool->splitProbability(*cluster); - clusterSplitP1 = splitProbObj.splitProbability(2); - clusterSplitP2 = splitProbObj.splitProbability(3); - ATH_MSG_VERBOSE( "Obtained split prob object with split prob: " << splitProbObj.splitProbability()); - if ( splitProbObj.splitProbability() > m_minSplitProbability ) { - ATH_MSG_VERBOSE( "Trying to split cluster ... "); - splitClusterParts = m_clusterSplitter->splitCluster(*cluster,splitProbObj); - } - } else if ( !m_clusterSplitter.empty() && (m_doIBLSplitting || m_IBLAbsent || !element->isBlayer())) - splitClusterParts = m_clusterSplitter->splitCluster(*cluster); - // check if splitting worked - clusterModified = !m_emulateSplitter && splitClusterParts.size() > 0; - clusterSplit = (splitClusterParts.size() > 1); - // exclusion: do not allow - singlePixelSplitCase = (groupSize==1) && clusterSplit; - - // CASE A: perform the actual split & create new clusters - if (clusterModified && !singlePixelSplitCase){ - if ( splitClusterParts.size() > 1) - ATH_MSG_VERBOSE( "--> Cluster with " << groupSize << " pixels is split into " << splitClusterParts.size() << " parts."); - else - ATH_MSG_VERBOSE( "--> Cluster is not actually split, but eventually modified, filling isSplit as: " << (clusterSplit || groupSize == 1)); - - std::vector<InDet::PixelCluster*> splitClusters; - // statistics output - if (clusterSplit){ - ++m_splitOrigClusters; - m_splitProdClusters += splitClusterParts.size(); - } else - ++m_modifiedOrigClusters; - - ATH_MSG_VERBOSE( "--> Non-zero cluster split size. Try to use new clusterization..."); - // iterate and make clusters - std::vector<InDet::PixelClusterParts>::iterator splitClusterPartsIter = splitClusterParts.begin(); - std::vector<InDet::PixelClusterParts>::iterator splitClusterPartsIterEnd = splitClusterParts.end(); - // use internal clustering if no position is estimated by clustersplitter - if ( !(*splitClusterPartsIter).localPosition()){ - - ATH_MSG_VERBOSE( "--> Position estimate from new clusterization not available... Use old clusterization"); - for ( ; splitClusterPartsIter != splitClusterPartsIterEnd; ++splitClusterPartsIter ){ - // make a new cluster, use standard clusterization to have a consistent clustering used - PixelCluster* splitCluster = makeCluster((*splitClusterPartsIter).identifierGroup(), - (*splitClusterPartsIter).totGroup(), - (*splitClusterPartsIter).lvl1Group(), - element, - pixelID, - ++clusterNumber, - (clusterSplit || groupSize == 1), - clusterSplitP1, - clusterSplitP2); - splitCluster->setHashAndIndex(clusterCollection->identifyHash(), clusterCollection->size()); - /** end new stuff */ - clusterCollection->push_back(splitCluster); - // @TODO: fill these clusters into an ambiguity map - splitClusters.push_back(splitCluster); - } - } else { - - std::vector<InDet::PixelClusterParts>::iterator splitClusterPartsIter = splitClusterParts.begin(); - std::vector<InDet::PixelClusterParts>::iterator splitClusterPartsIterEnd = splitClusterParts.end(); - - ATH_MSG_VERBOSE( "--> Processing new splitCluster with n. " << splitClusterParts.size() << " subClusters. "); - - for ( size_t iclus = 0 ; splitClusterPartsIter != splitClusterPartsIterEnd; ++splitClusterPartsIter, ++iclus ){ - - const Amg::Vector2D& position=(*(*splitClusterPartsIter).localPosition()); - const Amg::MatrixX& error=(*(*splitClusterPartsIter).errorMatrix()); - const std::vector<int>& totGroup=(*splitClusterPartsIter).totGroup(); - const std::vector<int>& lvl1Group=(*splitClusterPartsIter).lvl1Group(); - const std::vector<Identifier>& identifierGroup=(*splitClusterPartsIter).identifierGroup(); - - ATH_MSG_VERBOSE( "--> New Cluster :" << iclus << " - Position: " << position << " error: " << error ); - - InDetDD::SiLocalPosition subPos(position); - if (std::abs(subPos.xPhi())>mybounds->halflengthPhi()) { - double newxphi= (subPos.xPhi()>0) ? mybounds->halflengthPhi()-0.001 : -mybounds->halflengthPhi()+0.001; - subPos.xPhi(newxphi); - } - if (std::abs(subPos.xEta())>mybounds->halflengthEta()) { - double newxeta= (subPos.xEta()>0) ? mybounds->halflengthEta()-0.001 : -mybounds->halflengthEta()+0.001; - subPos.xEta(newxeta); - } - - const Identifier idSubCluster = element->identifierOfPosition(subPos); - - - std::vector<Identifier>::const_iterator rdosBegin = identifierGroup.begin(); - std::vector<Identifier>::const_iterator rdosEnd = identifierGroup.end(); - - int rowMin = int(2*(element->width()/element->phiPitch()))+1; - int rowMax = 0; - int colMin = int(2*(element->length()/element->etaPitch()))+1;; - int colMax = 0; - - for (; rdosBegin!= rdosEnd; ++rdosBegin) - { - Identifier rId = *rdosBegin; - int row = pixelID.phi_index(rId); - int col = pixelID.eta_index(rId); - - if (row < rowMin){ - rowMin = row; - } - if (row > rowMax){ - rowMax = row; - } - if (col < colMin){ - colMin = col; - } - if (col > colMax){ - colMax = col; - } - } - const InDetDD::PixelModuleDesign* design - (dynamic_cast<const InDetDD::PixelModuleDesign*>(&element->design())); - if (not design){ - ATH_MSG_ERROR("Dynamic cast failed at "<<__LINE__<<" of MergedPixelsTool.cxx."); - return nullptr; - } - int colWidth = colMax-colMin+1; - int rowWidth = rowMax-rowMin+1; - double etaWidth = design->widthFromColumnRange(colMin, colMax); - double phiWidth = design->widthFromRowRange(rowMin, rowMax); - SiWidth siWidth(Amg::Vector2D(rowWidth,colWidth), Amg::Vector2D(phiWidth,etaWidth) ); - - // create the new cluster - PixelCluster* splitCluster = new PixelCluster( - idSubCluster, - position, - identifierGroup, - lvl1Group[0], - totGroup, - cluster->chargeList(), - siWidth, - element, - new Amg::MatrixX(error), - cluster->omegax(),cluster->omegay(), - (clusterSplit || groupSize == 1), - clusterSplitP1, - clusterSplitP2 - ); - - splitCluster->setHashAndIndex(clusterCollection->identifyHash(), clusterCollection->size()); - clusterCollection->push_back(splitCluster); - splitClusters.push_back(splitCluster); - - }//end iteration on split clusters - } //end if no content in split clusters - - // fill the clusters into the ambiguity map - std::vector<InDet::PixelCluster*>::iterator iterOne = splitClusters.begin(); - std::vector<InDet::PixelCluster*>::iterator endOne = splitClusters.end(); - for ( ; iterOne != endOne; ++iterOne){ - std::vector<InDet::PixelCluster*>::iterator iterTwo = splitClusters.begin(); - std::vector<InDet::PixelCluster*>::iterator endTwo = splitClusters.end(); - for ( ; iterTwo != endTwo; ++iterTwo){ - if (iterTwo != iterOne) - (m_splitClusterMap)->insert(std::make_pair((*iterOne),(*iterTwo))); - } - } - // delete the original cluster - delete cluster; cluster = 0; - }//end if split clusters size is 0 - else if (m_emulateSplitter || singlePixelSplitCase){ - // create a new cluser with updated split information - if (singlePixelSplitCase) - ATH_MSG_VERBOSE( "--> Cluster is a single pixel cluser - no split performed, only fill split information - do not flag as split."); - else - ATH_MSG_VERBOSE( "--> Emulation mode: no split performed, only fill split information - flagged split."); - - PixelCluster* emulatedCluster = new PixelCluster( - cluster->identify(), - cluster->localPosition(), - cluster->rdoList(), - cluster->LVL1A(), - cluster->totList(), - cluster->chargeList(), - cluster->width(), - element, - new Amg::MatrixX(cluster->localCovariance()), - cluster->omegax(),cluster->omegay(), - (m_emulateSplitter ? clusterSplit : false), - clusterSplitP1, - clusterSplitP2 - ); - - emulatedCluster->setHashAndIndex(clusterCollection->identifyHash(), clusterCollection->size()); - clusterCollection->push_back(emulatedCluster); - //splitClusters.push_back(splitCluster); // check if we want it there or not - - // delete the original cluster - delete cluster; cluster = 0; - } else { - ATH_MSG_VERBOSE( "ZERO cluster split size is. Not replacing old cluster..."); - } - - - } - // no merging has been done; - if (cluster){ - // statistics output - if ((**group).size() >= m_maxSplitSize ) ++m_largeClusters; - - /** new: store hash id and index in the cluster: - * hash can be obtained from collection - * index is just the size of the coll before the push back - this is needed for later to make the EL to IDC valid in the - RIO_OnTrack objects set method might be temporary, this tool - (MergedPixelsTool could be friend of the cluster objects) - */ - cluster->setHashAndIndex(clusterCollection->identifyHash(), - clusterCollection->size()); - /** end new stuff */ - clusterCollection->push_back(cluster); - } - (**group).clear(); - } - delete *group; // now copied into cluster - *group = 0; - delete *totgroup; // won't be used any more (cluster class has not - // TOT info at the moment). - *totgroup=0; - totgroup++; - delete *lvl1group; // now used into cluster - *lvl1group=0; - lvl1group++; -} -return clusterCollection; -} + for(std::vector<RDO_Vector *>::iterator group = rdoGroups.begin() ; group!= rdoGroups.end() ; ++group) + { + + // the split probabilities + // writing the split boolean is done the following way + // - if in emulation mode: always write the output of the splitter for validation + // - if in pseudo-emulation mode (1-pixel clusters): set boolean to false, but keep split probs + bool clusterModified = false; + bool clusterSplit = false; + double clusterSplitP1 = 0.; + double clusterSplitP2 = 0.; + bool singlePixelSplitCase = 0; + // the size of the inition rdo group + size_t groupSize = (**group).size(); + + // ITk: check for clusters with ToT=size, this might be needed for upgrade simulation + bool isToTequalSize=false; + if(m_doRemoveClustersWithToTequalSize) + { + std::vector<int>::const_iterator tot_begin = (**totgroup).begin(); + std::vector<int>::const_iterator tot_end = (**totgroup).end(); + unsigned int cluster_total_tot=0; + for (; tot_begin!= tot_end; ++tot_begin) + { + cluster_total_tot += *tot_begin; + } + isToTequalSize = (groupSize==cluster_total_tot); + } + // -------- end of checks for clusters with ToT=size + + // If cluster is empty, i.e. it has been merged with another, + // do not attempt to make cluster. + if ( groupSize > 0 && isToTequalSize==false) // ITk: modification to remove clusters with ToT=size, may be needed for upgrade + // if ( groupSize > 0 ) + { + ++m_processedClusters; + // create the original cluster - split & split probs are by default false and 0 + PixelCluster* cluster = makeCluster(**group, + **totgroup, + **lvl1group, + element, + pixelID, + ++clusterNumber); + // check for splitting + if ( groupSize >= m_minSplitSize && groupSize <= m_maxSplitSize ) { + + // prepare for the return value of the pixel cluster + std::vector<InDet::PixelClusterParts> splitClusterParts; + if ( !m_splitProbTool.empty() && (m_doIBLSplitting || m_IBLAbsent || !element->isBlayer())){ + InDet::PixelClusterSplitProb splitProbObj = m_splitProbTool->splitProbability(*cluster); + clusterSplitP1 = splitProbObj.splitProbability(2); + clusterSplitP2 = splitProbObj.splitProbability(3); + ATH_MSG_VERBOSE( "Obtained split prob object with split prob: " << splitProbObj.splitProbability()); + if ( splitProbObj.splitProbability() > m_minSplitProbability ) { + ATH_MSG_VERBOSE( "Trying to split cluster ... "); + splitClusterParts = m_clusterSplitter->splitCluster(*cluster,splitProbObj); + } + } else if ( !m_clusterSplitter.empty() && (m_doIBLSplitting || m_IBLAbsent || !element->isBlayer())) + splitClusterParts = m_clusterSplitter->splitCluster(*cluster); + // check if splitting worked + clusterModified = !m_emulateSplitter && splitClusterParts.size() > 0; + clusterSplit = (splitClusterParts.size() > 1); + // exclusion: do not allow + singlePixelSplitCase = (groupSize==1) && clusterSplit; + + // CASE A: perform the actual split & create new clusters + if (clusterModified && !singlePixelSplitCase){ + if ( splitClusterParts.size() > 1) + ATH_MSG_VERBOSE( "--> Cluster with " << groupSize << " pixels is split into " << splitClusterParts.size() << " parts."); + else + ATH_MSG_VERBOSE( "--> Cluster is not actually split, but eventually modified, filling isSplit as: " << (clusterSplit || groupSize == 1)); + + std::vector<InDet::PixelCluster*> splitClusters; + // statistics output + if (clusterSplit){ + ++m_splitOrigClusters; + m_splitProdClusters += splitClusterParts.size(); + } else + ++m_modifiedOrigClusters; + + ATH_MSG_VERBOSE( "--> Non-zero cluster split size. Try to use new clusterization..."); + // iterate and make clusters + std::vector<InDet::PixelClusterParts>::iterator splitClusterPartsIter = splitClusterParts.begin(); + std::vector<InDet::PixelClusterParts>::iterator splitClusterPartsIterEnd = splitClusterParts.end(); + // use internal clustering if no position is estimated by clustersplitter + if ( !(*splitClusterPartsIter).localPosition()){ + + ATH_MSG_VERBOSE( "--> Position estimate from new clusterization not available... Use old clusterization"); + for ( ; splitClusterPartsIter != splitClusterPartsIterEnd; ++splitClusterPartsIter ){ + // make a new cluster, use standard clusterization to have a consistent clustering used + PixelCluster* splitCluster = makeCluster((*splitClusterPartsIter).identifierGroup(), + (*splitClusterPartsIter).totGroup(), + (*splitClusterPartsIter).lvl1Group(), + element, + pixelID, + ++clusterNumber, + (clusterSplit || groupSize == 1), + clusterSplitP1, + clusterSplitP2); + splitCluster->setHashAndIndex(clusterCollection->identifyHash(), clusterCollection->size()); + /** end new stuff */ + clusterCollection->push_back(splitCluster); + // @TODO: fill these clusters into an ambiguity map + splitClusters.push_back(splitCluster); + } + } else { + + std::vector<InDet::PixelClusterParts>::iterator splitClusterPartsIter = splitClusterParts.begin(); + std::vector<InDet::PixelClusterParts>::iterator splitClusterPartsIterEnd = splitClusterParts.end(); + + ATH_MSG_VERBOSE( "--> Processing new splitCluster with n. " << splitClusterParts.size() << " subClusters. "); + + for ( size_t iclus = 0 ; splitClusterPartsIter != splitClusterPartsIterEnd; ++splitClusterPartsIter, ++iclus ){ + + const Amg::Vector2D& position=(*(*splitClusterPartsIter).localPosition()); + const Amg::MatrixX& error=(*(*splitClusterPartsIter).errorMatrix()); + const std::vector<int>& totGroup=(*splitClusterPartsIter).totGroup(); + const std::vector<int>& lvl1Group=(*splitClusterPartsIter).lvl1Group(); + const std::vector<Identifier>& identifierGroup=(*splitClusterPartsIter).identifierGroup(); + + ATH_MSG_VERBOSE( "--> New Cluster :" << iclus << " - Position: " << position << " error: " << error ); + + InDetDD::SiLocalPosition subPos(position); + if (std::abs(subPos.xPhi())>mybounds->halflengthPhi()) { + double newxphi= (subPos.xPhi()>0) ? mybounds->halflengthPhi()-0.001 : -mybounds->halflengthPhi()+0.001; + subPos.xPhi(newxphi); + } + if (std::abs(subPos.xEta())>mybounds->halflengthEta()) { + double newxeta= (subPos.xEta()>0) ? mybounds->halflengthEta()-0.001 : -mybounds->halflengthEta()+0.001; + subPos.xEta(newxeta); + } + + const Identifier idSubCluster = element->identifierOfPosition(subPos); + + + std::vector<Identifier>::const_iterator rdosBegin = identifierGroup.begin(); + std::vector<Identifier>::const_iterator rdosEnd = identifierGroup.end(); + + int rowMin = int(2*(element->width()/element->phiPitch()))+1; + int rowMax = 0; + int colMin = int(2*(element->length()/element->etaPitch()))+1;; + int colMax = 0; + + for (; rdosBegin!= rdosEnd; ++rdosBegin) + { + Identifier rId = *rdosBegin; + int row = pixelID.phi_index(rId); + int col = pixelID.eta_index(rId); + + if (row < rowMin){ + rowMin = row; + } + if (row > rowMax){ + rowMax = row; + } + if (col < colMin){ + colMin = col; + } + if (col > colMax){ + colMax = col; + } + } + const InDetDD::PixelModuleDesign* design + (dynamic_cast<const InDetDD::PixelModuleDesign*>(&element->design())); + if (not design){ + ATH_MSG_ERROR("Dynamic cast failed at "<<__LINE__<<" of MergedPixelsTool.cxx."); + return nullptr; + } + int colWidth = colMax-colMin+1; + int rowWidth = rowMax-rowMin+1; + double etaWidth = design->widthFromColumnRange(colMin, colMax); + double phiWidth = design->widthFromRowRange(rowMin, rowMax); + SiWidth siWidth(Amg::Vector2D(rowWidth,colWidth), Amg::Vector2D(phiWidth,etaWidth) ); + + // create the new cluster + PixelCluster* splitCluster = new PixelCluster( + idSubCluster, + position, + identifierGroup, + lvl1Group[0], + totGroup, + cluster->chargeList(), + siWidth, + element, + new Amg::MatrixX(error), + cluster->omegax(),cluster->omegay(), + (clusterSplit || groupSize == 1), + clusterSplitP1, + clusterSplitP2 + ); + + splitCluster->setHashAndIndex(clusterCollection->identifyHash(), clusterCollection->size()); + clusterCollection->push_back(splitCluster); + splitClusters.push_back(splitCluster); + + }//end iteration on split clusters + } //end if no content in split clusters + + + // delete the original cluster + delete cluster; cluster = 0; + }//end if split clusters size is 0 + else if (m_emulateSplitter || singlePixelSplitCase){ + // create a new cluser with updated split information + if (singlePixelSplitCase) + ATH_MSG_VERBOSE( "--> Cluster is a single pixel cluser - no split performed, only fill split information - do not flag as split."); + else + ATH_MSG_VERBOSE( "--> Emulation mode: no split performed, only fill split information - flagged split."); + + PixelCluster* emulatedCluster = new PixelCluster( + cluster->identify(), + cluster->localPosition(), + cluster->rdoList(), + cluster->LVL1A(), + cluster->totList(), + cluster->chargeList(), + cluster->width(), + element, + new Amg::MatrixX(cluster->localCovariance()), + cluster->omegax(),cluster->omegay(), + (m_emulateSplitter ? clusterSplit : false), + clusterSplitP1, + clusterSplitP2 + ); + + emulatedCluster->setHashAndIndex(clusterCollection->identifyHash(), clusterCollection->size()); + clusterCollection->push_back(emulatedCluster); + //splitClusters.push_back(splitCluster); // check if we want it there or not + + // delete the original cluster + delete cluster; cluster = 0; + } else { + ATH_MSG_VERBOSE( "ZERO cluster split size is. Not replacing old cluster..."); + } + + + } + // no merging has been done; + if (cluster){ + // statistics output + if ((**group).size() >= m_maxSplitSize ) ++m_largeClusters; + + /** new: store hash id and index in the cluster: + * hash can be obtained from collection + * index is just the size of the coll before the push back + this is needed for later to make the EL to IDC valid in the + RIO_OnTrack objects set method might be temporary, this tool + (MergedPixelsTool could be friend of the cluster objects) + */ + cluster->setHashAndIndex(clusterCollection->identifyHash(), + clusterCollection->size()); + /** end new stuff */ + clusterCollection->push_back(cluster); + } + (**group).clear(); + } + delete *group; // now copied into cluster + *group = 0; + delete *totgroup; // won't be used any more (cluster class has not + // TOT info at the moment). + *totgroup=0; + totgroup++; + delete *lvl1group; // now used into cluster + *lvl1group=0; + lvl1group++; + } + return clusterCollection; + } //----------------------------------------------------------------------- // Once the lists of RDOs which makes up the clusters have been found by the @@ -664,7 +697,7 @@ PixelCluster* MergedPixelsTool::makeCluster double splitProb1, double splitProb2) const { - if ( outputLevel() < MSG::INFO ) + if ( msgLevel() < MSG::INFO ) ATH_MSG_VERBOSE("makeCluster called, number " << clusterNumber); Identifier gangedID; @@ -717,7 +750,7 @@ PixelCluster* MergedPixelsTool::makeCluster int realtot = *tot; if( m_IBLParameterSvc->containsIBL() && pixelID.barrel_ec(rId) == 0 && pixelID.layer_disk(rId) == 0 ) { if (*tot >= m_overflowIBLToT ) realtot = m_overflowIBLToT; - msg(MSG::DEBUG) << "barrel_ec = " << pixelID.barrel_ec(rId) << " layer_disque = " << pixelID.layer_disk(rId) << " ToT = " << *tot << " Real ToT = " << realtot << endreq; + msg(MSG::DEBUG) << "barrel_ec = " << pixelID.barrel_ec(rId) << " layer_disque = " << pixelID.layer_disk(rId) << " ToT = " << *tot << " Real ToT = " << realtot << endmsg; } if (row == rowMin) qRowMin += realtot; @@ -833,14 +866,14 @@ PixelCluster* MergedPixelsTool::makeCluster Amg::Vector2D position(centroid.xPhi(),centroid.xEta()); if (msgLvl(MSG::VERBOSE)){ - msg() << "Cluster ID =" << id << endreq; + msg() << "Cluster ID =" << id << endmsg; msg() << "Cluster width (eta x phi) = " << colWidth - << " x " << rowWidth << endreq; + << " x " << rowWidth << endmsg; msg() << "Cluster width (eta x phi) = " << etaWidth - << " x " << phiWidth << endreq; + << " x " << phiWidth << endmsg; msg() << "Cluster local position (eta,phi) = " << (position)[0] << " " - << (position)[1] << endreq; + << (position)[1] << endmsg; } if(!m_clusterMaker){ @@ -893,30 +926,6 @@ bool MergedPixelsTool::isGanged(const Identifier& rdoID, } -void MergedPixelsTool::handle(const Incident& inc) { - - - - // the cluster ambiguity map - if ( inc.type() == IncidentType::BeginEvent ){ - // record the Split ambiguity map - if (!m_splitClusterMap && !m_clusterSplitter.empty() && !m_emulateSplitter){ - m_splitClusterMap = new InDet::PixelGangedClusterAmbiguities; - if ( evtStore()->record(m_splitClusterMap,m_splitClusterMapName).isFailure()){ - ATH_MSG_WARNING("Could not record split cluster ambiguity map."); - delete m_splitClusterMap; m_splitClusterMap = 0; - } else - ATH_MSG_VERBOSE("Cluster split ambiguity map recorded as '" << m_splitClusterMapName <<"'."); - } - } - - - if ( inc.type() == IncidentType::EndEvent ){ - ATH_MSG_VERBOSE("'EndEvent' incident caught. Refreshing Cache."); - m_splitClusterMap = 0; - } -} - //----------------------------------------------------------------------- // Checks if two RDO lists (would be clusters) whould be merged, @@ -943,7 +952,13 @@ void MergedPixelsTool::checkForMerge(const Identifier& id, for (; nextGroup!= lastGroup; ++nextGroup) { - if (areNeighbours(**nextGroup, id, element, pixelID)) + // modified for ITk samples: + // mergeTwoClusters() checks if neighboring clusters are not too large to be merged + // mergeTwoBrokenClusters() checks if two separated clusters are from the same broken cluster + if ((areNeighbours(**nextGroup, id, element, pixelID) + && mergeTwoClusters(**baseGroup,**nextGroup,element,pixelID)) + || mergeTwoBrokenClusters(**baseGroup,**nextGroup,element,pixelID)) + // if (areNeighbours(**nextGroup, id, element, pixelID)) { // merge the RDO identifier groups idIterator firstID = (*nextGroup)->begin(); @@ -976,6 +991,246 @@ void MergedPixelsTool::checkForMerge(const Identifier& id, } } + +// ------------ New functions to merge broken clusters for ITk upgrade studies +bool MergedPixelsTool::mergeTwoBrokenClusters(const std::vector<Identifier>& group1, + const std::vector<Identifier>& group2, + InDetDD::SiDetectorElement* element, + const PixelID& pixelID) const +{ + bool mergeClusters=false; + if(m_doMergeBrokenClusters) + { + std::vector<Identifier>::const_iterator group1Begin = group1.begin(); // first pixel in cluster-1 + std::vector<Identifier>::const_iterator group1End = group1.end(); // last pixel in cluster-1 + std::vector<Identifier>::const_iterator group2Begin = group2.begin(); // first pixel in cluster-2 + std::vector<Identifier>::const_iterator group2End = group2.end(); // last pixel in cluster-2 + Identifier id1b= *group1Begin; // first pixel in cluster1 + Identifier id2b= *group2Begin; // first pixel in cluster2 + if(pixelID.is_barrel(id1b) && pixelID.is_barrel(id2b)) // make sure both clusters are in barrel + { + int cl1_rowMin= 1000; + int cl1_rowMax= 0; + int cl2_rowMin= 1000; + int cl2_rowMax= 0; + int cl1_colMin= 1000; + int cl1_colMax= 0; + int cl2_colMin= 1000; + int cl2_colMax= 0; + + float cl1_ave_row=0.0; // average row # in cluster1 + float cl2_ave_row=0.0; // average row # in cluster2 + unsigned int cl1_size=group1.size(); // number of pixels in cluster-1 + unsigned int cl2_size=group2.size(); // number of pixels in cluster-2 + + // determin cluster-1 boundaries + for(; group1Begin!=group1End; ++group1Begin) + { + Identifier id=*group1Begin; + int row=pixelID.phi_index(id); + int col=pixelID.eta_index(id); + + cl1_ave_row=cl1_ave_row + row/(1.0*cl1_size); + + if(col<cl1_colMin) cl1_colMin=col; + if(col>cl1_colMax) cl1_colMax=col; + if(row<cl1_rowMin) cl1_rowMin=row; + if(row>cl1_rowMax) cl1_rowMax=row; + } + + // determin cluster-2 boundaries + for(; group2Begin!=group2End; ++group2Begin) + { + Identifier id=*group2Begin; + int row=pixelID.phi_index(id); + int col=pixelID.eta_index(id); + + cl2_ave_row=cl2_ave_row + row/(1.0*cl2_size); + + if(col<cl2_colMin) cl2_colMin=col; + if(col>cl2_colMax) cl2_colMax=col; + if(row<cl2_rowMin) cl2_rowMin=row; + if(row>cl2_rowMax) cl2_rowMax=row; + } + + int cl1_sizeZ=cl1_colMax - cl1_colMin + 1; + int cl2_sizeZ=cl2_colMax - cl2_colMin + 1; + int cl1_sizePhi=cl1_rowMax - cl1_rowMin + 1; + int cl2_sizePhi=cl2_rowMax - cl2_rowMin + 1; + + int col_last= 0; + int col_first= 0; + if(cl1_colMax<cl2_colMax) + { + col_last=cl1_colMax; + col_first=cl2_colMin; + } + else + { + col_last=cl2_colMax; + col_first=cl1_colMin; + } + + int clMerged_sizePhi= std::max(cl1_rowMax,cl2_rowMax) - std::min(cl1_rowMin,cl2_rowMin) + 1; + + // these are intial values for merging, to be optimized in studies + if( (cl1_colMin<=cl2_colMin && cl2_colMin<=cl1_colMax) + || (cl2_colMin<=cl1_colMin && cl1_colMin<=cl2_colMax) ) return mergeClusters; // don't merge, clusters overlap in column number; original default + if(cl1_sizePhi>cl1_sizeZ || cl2_sizePhi>cl2_sizeZ) return mergeClusters; // for now, don't merge clusters if sizePhi>sizeZ; original default + + int maxNmis=maxGap(std::min(cl1_colMin,cl2_colMin),std::max(cl1_colMax,cl2_colMax),cl1_rowMin,element); + if(col_first - col_last > maxNmis) return mergeClusters; // don't merge clusters because gap is too large + if(col_first - col_last > std::max(cl1_sizeZ,cl2_sizeZ)) return mergeClusters; // don't merge clusters where gap in eta is larger than max size + + if(checkSizeZ(cl1_colMin,cl1_colMax,cl1_rowMin,element)>0) return mergeClusters; // don't merge, cluster1 is too large; original default + if(checkSizeZ(cl2_colMin,cl2_colMax,cl2_rowMin,element)>0) return mergeClusters; // don't merge, cluster2 is too large; original default + if(checkSizeZ(std::min(cl1_colMin,cl2_colMin),std::max(cl1_colMax,cl2_colMax),cl2_rowMin,element)>0) return mergeClusters; // don't merge, new cluster is too large; original default + + if(clMerged_sizePhi>(cl1_sizePhi+cl2_sizePhi-1)) return mergeClusters; // don't merge, cluster1 is too large in sizePhi; to be replaced by eta-dependent cut in the future; original default + // Note that the last condition (if satisfied) should always be mergeClusters=true + if(fabs(cl1_ave_row - cl2_ave_row)<1.0) mergeClusters=true; // merge clusters only if ave_row(CL1) and ave_row(CL2) are less than one pixel apart (test-2) + } + } + return mergeClusters; +} + + +bool MergedPixelsTool::mergeTwoClusters(const std::vector<Identifier>& group1, + const std::vector<Identifier>& group2, + InDetDD::SiDetectorElement* element, + const PixelID& pixelID) const +{ + bool mergeClusters=true; + if(m_doCheckSizeBeforeMerging) + { + std::vector<Identifier>::const_iterator group1Begin = group1.begin(); // first pixel in cluster-1 + std::vector<Identifier>::const_iterator group1End = group1.end(); // last pixel in cluster-1 + std::vector<Identifier>::const_iterator group2Begin = group2.begin(); // first pixel in cluster-2 + std::vector<Identifier>::const_iterator group2End = group2.end(); // last pixel in cluster-2 + Identifier id1b= *group1Begin; // first pixel in cluster1 + Identifier id2b= *group2Begin; // first pixel in cluster2 + if(pixelID.is_barrel(id1b) && pixelID.is_barrel(id2b)) // make sure both clusters are in barrel + { + int cl1_rowMin= 1000; + int cl1_rowMax= 0; + int cl2_rowMin= 1000; + int cl2_rowMax= 0; + int cl1_colMin= 1000; + int cl1_colMax= 0; + int cl2_colMin= 1000; + int cl2_colMax= 0; + + // determin cluster-1 boundaries + for(; group1Begin!=group1End; ++group1Begin) + { + Identifier id=*group1Begin; + int row=pixelID.phi_index(id); + int col=pixelID.eta_index(id); + if(col<cl1_colMin) cl1_colMin=col; + if(col>cl1_colMax) cl1_colMax=col; + if(row<cl1_rowMin) cl1_rowMin=row; + if(row>cl1_rowMax) cl1_rowMax=row; + } + + // determin cluster-2 boundaries + for(; group2Begin!=group2End; ++group2Begin) + { + Identifier id=*group2Begin; + int row=pixelID.phi_index(id); + int col=pixelID.eta_index(id); + if(col<cl2_colMin) cl2_colMin=col; + if(col>cl2_colMax) cl2_colMax=col; + if(row<cl2_rowMin) cl2_rowMin=row; + if(row>cl2_rowMax) cl2_rowMax=row; + } + + // it doesn't matter that much which row number is used as input to checkSizeZ() + if(checkSizeZ(cl1_colMin,cl1_colMax,cl1_rowMin,element)>0) return false; // don't merge, cluster1 is too large (original default) + if(checkSizeZ(cl2_colMin,cl2_colMax,cl2_rowMin,element)>0) return false; // don't merge, cluster2 is too large (original default) + if(checkSizeZ(std::min(cl1_colMin,cl2_colMin),std::max(cl1_colMax,cl2_colMax),cl2_rowMin,element)>0) return false; // don't merge, new cluster is too large + } + } + return mergeClusters; +} + +// checkSizeZ compares cluster sizeZ with expected cluster size for this cluster position (+/-200 mm for beam spread) +// checkSizeZ()=-1 if cluster is too small +// checkSizeZ()=0 if cluster sizeZ is within allowed range +// checkSizeZ()=1 if cluster is too large +// in the future, it may be changed to return deltaSizeZ +int MergedPixelsTool::checkSizeZ(int colmin, int colmax, int row, InDetDD::SiDetectorElement* element) const +{ + int pass_code=0; + + // const float m_beam_spread=200.0; + const InDetDD::PixelModuleDesign* design(dynamic_cast<const InDetDD::PixelModuleDesign*>(&element->design())); + if (not design) + { + ATH_MSG_ERROR("Dynamic cast failed at "<<__LINE__<<" of MergedPixelsTool. This cluster will not be considered for merging"); + return 1; // don't merge + } + + // calculating cluster position in the global coordinate system + InDetDD::SiLocalPosition sumOfPositions(0,0,0); + InDetDD::SiLocalPosition siLocalPosition1(design->positionFromColumnRow(colmin,row)); + InDetDD::SiLocalPosition siLocalPosition2(design->positionFromColumnRow(colmax,row)); + sumOfPositions = siLocalPosition1+siLocalPosition2; + InDetDD::SiLocalPosition centroid(sumOfPositions/2); + Amg::Vector3D globalPos = element->globalPosition(centroid); + // calculating min and max size of a cluster assuming its global position and +/-20 cm luminous region + float sensorThickness = element->thickness(); + float pitch = design->etaPitch(); + + if(globalPos.perp()>0.0) + { + float length1= fabs(globalPos.z()+m_beam_spread)*sensorThickness/globalPos.perp(); + float length2= fabs(globalPos.z()-m_beam_spread)*sensorThickness/globalPos.perp(); + if((colmax-colmin+1)*pitch > (std::max(length1,length2)+pitch)) pass_code=1; + if((colmax-colmin+1)*pitch < (std::max(length1,length2)-pitch)) pass_code=-1; + } + else return 1; // don't merge + + return pass_code; +} + +// this function returns expected sizeZ +int MergedPixelsTool::expectedSizeZ(int colmin, int colmax, int row, InDetDD::SiDetectorElement* element) const { + int exp_sizeZ=1; + const InDetDD::PixelModuleDesign* design(dynamic_cast<const InDetDD::PixelModuleDesign*>(&element->design())); + if (not design) + { + ATH_MSG_ERROR("Dynamic cast failed at "<<__LINE__<<" of MergedPixelsTool. This cluster will not be considered for merging"); + return 1; + } + // calculating cluster position in the global coordinate system + InDetDD::SiLocalPosition sumOfPositions(0,0,0); + InDetDD::SiLocalPosition siLocalPosition1(design->positionFromColumnRow(colmin,row)); + InDetDD::SiLocalPosition siLocalPosition2(design->positionFromColumnRow(colmax,row)); + sumOfPositions = siLocalPosition1+siLocalPosition2; + InDetDD::SiLocalPosition centroid(sumOfPositions/2); + Amg::Vector3D globalPos = element->globalPosition(centroid); + // calculating min and max size of a cluster assuming its global position and +/-20 cm luminous region + float sensorThickness = element->thickness(); + float pitch = design->etaPitch(); + float length1= fabs(globalPos.z()+m_beam_spread)*sensorThickness/globalPos.perp(); + float length2= fabs(globalPos.z()-m_beam_spread)*sensorThickness/globalPos.perp(); + float length=std::max(length1,length2); + exp_sizeZ=int(length/pitch+0.5); + return exp_sizeZ; +} + +// this function returns size of the maximum gap between two cluster fragments +int MergedPixelsTool::maxGap(int colmin, int colmax, int row, InDetDD::SiDetectorElement* element) const { + int Nmis=1; + int Nexp=expectedSizeZ(colmin,colmax,row,element); + if(Nexp<=3) return 0; // there should not be any gap for very small clusters + int gap=(int)((log(m_lossThreshold)-log(1.0*(Nexp-2)))/log(1.0-m_pixelEff)+0.5); + Nmis=std::min(gap,Nexp-2); + return Nmis; +} + + + } //---------------------------------------------------------------------------- diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/src/NnClusterizationFactory.cxx b/InnerDetector/InDetRecTools/SiClusterizationTool/src/NnClusterizationFactory.cxx index 8b16b751ee70949966616731d5d261aef2338e90..4b45318efb2ce5d5515053c67c5b9843d7bf9787 100644 --- a/InnerDetector/InDetRecTools/SiClusterizationTool/src/NnClusterizationFactory.cxx +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/src/NnClusterizationFactory.cxx @@ -63,7 +63,7 @@ namespace InDet { m_layerPrefix("Layer"), m_weightIndicator("_weights"), m_thresholdIndicator("_thresholds"), - m_networkToHistoTool("Trk::NeuralNetworkToHistoTool/NeuralNetworkToHistoTool"), + m_networkToHistoTool("Trk::NeuralNetworkToHistoTool/NeuralNetworkToHistoTool", this), m_calibSvc("PixelCalibSvc", name), m_useToT(true), m_addIBL(false), @@ -1002,7 +1002,9 @@ if(m_doRunI){ return assembleInputRunI( input, sizeX, sizeY ); }els for ( ; rdosBegin!= rdosEnd && tot != totList.end(); ++tot, ++rdosBegin, ++totRecreated ){ // recreate the charge: should be a method of the calibSvc int tot0 = *tot; + float ch = m_calibSvc->getCharge(*rdosBegin,tot0); + chListRecreated.push_back(ch); totListRecreated.push_back(tot0); } diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/src/NnPixelClusterSplitProbTool.cxx b/InnerDetector/InDetRecTools/SiClusterizationTool/src/NnPixelClusterSplitProbTool.cxx index 68b0ad433e873dc66601526778017ee34b7c3969..ae145d6412162580a9c8fc2adf58042535911966 100644 --- a/InnerDetector/InDetRecTools/SiClusterizationTool/src/NnPixelClusterSplitProbTool.cxx +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/src/NnPixelClusterSplitProbTool.cxx @@ -28,7 +28,7 @@ namespace InDet NnPixelClusterSplitProbTool::NnPixelClusterSplitProbTool(const std::string& t, const std::string& n, const IInterface* p) :AthAlgTool(t,n,p), - m_NnClusterizationFactory("InDet::NnClusterizationFactory/NnClusterizationFactory"), + m_NnClusterizationFactory("InDet::NnClusterizationFactory/NnClusterizationFactory", this), m_iBeamCondSvc("BeamCondSvc",n), m_useBeamSpotInfo(true) { diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/src/NnPixelClusterSplitter.cxx b/InnerDetector/InDetRecTools/SiClusterizationTool/src/NnPixelClusterSplitter.cxx index 1258ae075cdddb64b00ff4d6fbb183f6c1121084..651263b4c6e106921bb7b0f70be1f4ed110e5e75 100644 --- a/InnerDetector/InDetRecTools/SiClusterizationTool/src/NnPixelClusterSplitter.cxx +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/src/NnPixelClusterSplitter.cxx @@ -34,7 +34,7 @@ InDet::NnPixelClusterSplitter::NnPixelClusterSplitter(const std::string &type, const std::string &name, const IInterface *parent) : AthAlgTool(type,name,parent), - m_NnClusterizationFactory("InDet::NnClusterizationFactory/NnClusterizationFactory"), + m_NnClusterizationFactory("InDet::NnClusterizationFactory/NnClusterizationFactory", this), m_iBeamCondSvc("BeamCondSvc",name), m_useBeamSpotInfo(true) { diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/src/PixelClusteringToolBase.cxx b/InnerDetector/InDetRecTools/SiClusterizationTool/src/PixelClusteringToolBase.cxx index 2f23552e86707663c5b9d3dae795f88fd747af7d..786a00bc0411e948d25c9aed01878c39f95a349c 100755 --- a/InnerDetector/InDetRecTools/SiClusterizationTool/src/PixelClusteringToolBase.cxx +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/src/PixelClusteringToolBase.cxx @@ -20,7 +20,7 @@ namespace InDet const std::string &name, const IInterface *parent) : AthAlgTool(type,name,parent), - m_clusterMaker("InDet::ClusterMakerTool"), + m_clusterMaker("InDet::ClusterMakerTool", this), m_posStrategy(0), m_errorStrategy(1), m_acceptDiagonalClusters(1), @@ -41,27 +41,27 @@ namespace InDet StatusCode PixelClusteringToolBase::initialize() { - msg(MSG::INFO) << "initialize()" << endreq; + msg(MSG::INFO) << "initialize()" << endmsg; if ( m_clusterMaker.retrieve().isFailure() ) { - msg(MSG::FATAL) << m_clusterMaker.propertyName() << ": Failed to retrieve tool " << m_clusterMaker.type() << endreq; + msg(MSG::FATAL) << m_clusterMaker.propertyName() << ": Failed to retrieve tool " << m_clusterMaker.type() << endmsg; return StatusCode::FAILURE; } else { - msg(MSG::INFO) << m_clusterMaker.propertyName() << ": Retrieved tool " << m_clusterMaker.type() << endreq; + msg(MSG::INFO) << m_clusterMaker.propertyName() << ": Retrieved tool " << m_clusterMaker.type() << endmsg; } StatusCode sc = m_summarySvc.retrieve(); if (sc.isFailure() || !m_summarySvc) { - msg(MSG::WARNING) << m_summarySvc.type() << " not found! "<< endreq; + msg(MSG::WARNING) << m_summarySvc.type() << " not found! "<< endmsg; if ( m_usePixelMap || m_useModuleMap ) { - msg(MSG::FATAL) << m_summarySvc.type() << " is compulsory with this tool configuration" << endreq; + msg(MSG::FATAL) << m_summarySvc.type() << " is compulsory with this tool configuration" << endmsg; return StatusCode::FAILURE; } } else{ - msg(MSG::INFO) << "Retrieved service " << m_summarySvc.type() << endreq; + msg(MSG::INFO) << "Retrieved service " << m_summarySvc.type() << endmsg; } return StatusCode::SUCCESS; diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/src/PixelGangedAmbiguitiesFinder.cxx b/InnerDetector/InDetRecTools/SiClusterizationTool/src/PixelGangedAmbiguitiesFinder.cxx index 9f14c4487d721f57b5b78a7b19d0d63de5fa6d24..1546b9a4a416edf3353965781d188b1099f639aa 100755 --- a/InnerDetector/InDetRecTools/SiClusterizationTool/src/PixelGangedAmbiguitiesFinder.cxx +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/src/PixelGangedAmbiguitiesFinder.cxx @@ -61,8 +61,8 @@ void PixelGangedAmbiguitiesFinder::execute( if (collection->size()<2) return; if (msgLvl(MSG::DEBUG)){ - msg() << collection->size() << " clusters" << endreq; - msg() << "The map has " << theMap.size() << " entries already" << endreq; + msg() << collection->size() << " clusters" << endmsg; + msg() << "The map has " << theMap.size() << " entries already" << endmsg; } Identifier elementID = collection->identify(); @@ -120,7 +120,7 @@ void PixelGangedAmbiguitiesFinder::execute( if(hasGanged){ - if(msgLvl(MSG::DEBUG)) msg() << "Ganged pixel, find combi..." << endreq; + if(msgLvl(MSG::DEBUG)) msg() << "Ganged pixel, find combi..." << endmsg; std::vector<Identifier>::const_iterator gangedPixelsBegin = gangedPixels.begin(); std::vector<Identifier>::const_iterator gangedPixelsEnd = gangedPixels.end(); PixelClusterCollection::iterator cluster2=cluster+1; @@ -128,7 +128,7 @@ void PixelGangedAmbiguitiesFinder::execute( if(msgLvl(MSG::DEBUG)) msg() << "Comparing " << std::hex << (*cluster)->identify() << " and " << (*cluster2)->identify() - << std::dec << endreq; + << std::dec << endmsg; bool sharedGanged = false; int rmin2=999, cmin2=999; // bottom left corner of cluster2 const std::vector<Identifier>& rdos2 = (*cluster2)->rdoList(); @@ -167,7 +167,7 @@ void PixelGangedAmbiguitiesFinder::execute( rmList.push_back(cluster); if(msgLvl(MSG::DEBUG)) msg() << std::hex << ": deleted " << (*cluster)->identify() - << std::dec << endreq; + << std::dec << endmsg; break; } else (*cluster)->setFake(true); @@ -178,7 +178,7 @@ void PixelGangedAmbiguitiesFinder::execute( rmList.push_back(cluster2); if(msgLvl(MSG::DEBUG)) msg() << std::hex << ": deleted " << (*cluster2)->identify() - << std::dec << endreq; + << std::dec << endmsg; continue; } else (*cluster2)->setFake(true); @@ -243,7 +243,7 @@ void PixelGangedAmbiguitiesFinder::execute( rmList.push_back(cluster2); if(msgLvl(MSG::DEBUG)) msg() << std::hex << ": deleted " << (*cluster2)->identify() - << std::dec << endreq; + << std::dec << endmsg; continue; } else { @@ -257,7 +257,7 @@ void PixelGangedAmbiguitiesFinder::execute( rmList.push_back(cluster); if(msgLvl(MSG::DEBUG)) msg() << std::hex << ": deleted " << (*cluster)->identify() - << std::dec << endreq; + << std::dec << endmsg; continue; } else { @@ -287,7 +287,7 @@ void PixelGangedAmbiguitiesFinder::execute( theMap.insert(std::make_pair(*cluster2,*cluster)); if(msgLvl(MSG::DEBUG)) msg() << std::hex << ": added ambiguity entry" - << std::dec << endreq; + << std::dec << endmsg; } } } @@ -296,7 +296,7 @@ void PixelGangedAmbiguitiesFinder::execute( if (msgLvl(MSG::DEBUG)) { msg(MSG::DEBUG) << "The map has " << theMap.size() << " entries " - << endreq; + << endmsg; } @@ -311,11 +311,11 @@ void PixelGangedAmbiguitiesFinder::execute( Identifier gangedID; if(msgLvl(MSG::DEBUG)) msg() << "Removed " << rmNumber+1 << " cluster: " << std::hex << (*(*rmit-rmNumber))->identify() << std::dec - << endreq ; + << endmsg ; collection->erase(*rmit-rmNumber); // The position of the iterator rmNumber++; } - if(msgLvl(MSG::DEBUG)) msg() << rmNumber << " fake clusters from ganged pixel have been removed" << endreq; + if(msgLvl(MSG::DEBUG)) msg() << rmNumber << " fake clusters from ganged pixel have been removed" << endmsg; } diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/src/SCT_ClusteringTool.cxx b/InnerDetector/InDetRecTools/SiClusterizationTool/src/SCT_ClusteringTool.cxx index e2b5772ab5dc92382fe8225991a28a35a66351d3..e69a506b6cdcbe112dea36a9eeaed732fc237495 100755 --- a/InnerDetector/InDetRecTools/SiClusterizationTool/src/SCT_ClusteringTool.cxx +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/src/SCT_ClusteringTool.cxx @@ -56,8 +56,12 @@ namespace InDet{ m_errorStrategy(1), m_checkBadChannels(true), m_conditionsSvc("SCT_ConditionsSummarySvc",name), - m_clusterMaker("InDet::ClusterMakerTool"), - m_timeBinStr("") + m_clusterMaker("InDet::ClusterMakerTool", this), + m_timeBinStr(""), + m_innermostBarrelX1X(false), + m_innertwoBarrelX1X(false), + m_majority01X(false), + m_useRowInformation(false) { declareInterface<ISCT_ClusteringTool>(this); declareProperty("globalPosAlg" , m_clusterMaker); @@ -65,6 +69,10 @@ namespace InDet{ declareProperty("checkBadChannels",m_checkBadChannels); declareProperty("conditionsService" , m_conditionsSvc); declareProperty("timeBins" , m_timeBinStr); + declareProperty("majority01X" , m_majority01X); + declareProperty("innermostBarrelX1X" , m_innermostBarrelX1X); + declareProperty("innertwoBarrelX1X" , m_innertwoBarrelX1X); + declareProperty("useRowInformation", m_useRowInformation); } StatusCode SCT_ClusteringTool::decodeTimeBins() { @@ -75,7 +83,7 @@ namespace InDet{ if (m_timeBinStr.empty()) return StatusCode::SUCCESS; if (m_timeBinStr.size() != 3) { - msg(MSG::ERROR) << "Time bine string must only contain 3 bins" << endreq; + msg(MSG::ERROR) << "Time bine string must only contain 3 bins" << endmsg; return StatusCode::FAILURE; } std::transform(m_timeBinStr.begin(), m_timeBinStr.end(), m_timeBinStr.begin(), ::toupper); @@ -95,7 +103,7 @@ namespace InDet{ if (timeBin == '0') {bit = 0; return StatusCode::SUCCESS;} if (timeBin == '1') {bit = 1; return StatusCode::SUCCESS;} - msg(MSG::ERROR) << "Invalid time bin string " << timeBin << endreq; + msg(MSG::ERROR) << "Invalid time bin string " << timeBin << endmsg; return StatusCode::FAILURE; } @@ -112,8 +120,31 @@ namespace InDet{ return pass; } + bool SCT_ClusteringTool::testTimeBins01X(int timeBin) const { + // Convert the given timebin to a bit set and test each bit + // if bit is -1 (i.e. X) it always passes, otherwise require exact match of 0/1 + // N.B bitset has opposite order to the bit pattern we define + + bool pass(true); + std::bitset<3> timePattern(static_cast<unsigned long>(timeBin)); + if (timePattern.test(2) != bool(0)) pass=false; + if (timePattern.test(1) != bool(1)) pass=false; + return pass; + } + + bool SCT_ClusteringTool::testTimeBinsX1X(int timeBin) const { + // Convert the given timebin to a bit set and test each bit + // if bit is -1 (i.e. X) it always passes, otherwise require exact match of 0/1 + // N.B bitset has opposite order to the bit pattern we define + + bool pass(true); + std::bitset<3> timePattern(static_cast<unsigned long>(timeBin)); + if (timePattern.test(1) != bool(1)) pass=false; + return pass; + } + StatusCode SCT_ClusteringTool::initialize(){ - msg(MSG::INFO) << "Initialize clustering tool" << endreq; + msg(MSG::INFO) << "Initialize clustering tool" << endmsg; static const StatusCode fail(StatusCode::FAILURE); @@ -123,6 +154,33 @@ namespace InDet{ } if (decodeTimeBins().isFailure()) return StatusCode::FAILURE; + + /// only one of m_majority01X, m_innermostBarrelX1X and m_innertwoBarrelX1X + /// can be true - check if this is the case + /// If none of them is true, m_timeBinStr is used. + int countTrueSettings=0; + if (m_majority01X) countTrueSettings++; + if (m_innermostBarrelX1X) countTrueSettings++; + if (m_innertwoBarrelX1X) countTrueSettings++; + if (countTrueSettings !=1) { + if(!m_timeBinStr.empty()) { + ATH_MSG_INFO("Timing requirement: m_timeBinStr " << m_timeBinStr << " is used for clustering"); + } else { + if(countTrueSettings==0) { + ATH_MSG_INFO("Timing requirement is not used for clustering"); + } else { + ATH_MSG_ERROR("One and only one of m_majority01X, m_innermostBarrelX1X and m_innertwoBarrelX1X should be set to True!"); + return StatusCode::FAILURE; + } + } + } else { + ATH_MSG_INFO("Timing requirement: " << + (m_majority01X ? "m_majority01X" : "") << + (m_innermostBarrelX1X ? "m_innermostBarrelX1X" : "") << + (m_innertwoBarrelX1X ? "m_innertwoBarrelX1X" : "") << + " is true and used for clustering"); + } + return StatusCode::SUCCESS; } @@ -189,6 +247,42 @@ namespace InDet{ // After this, the cluster vector is either empty or has the last good cluster } } + + void SCT_ClusteringTool::addStripsToClusterInclRows(const Identifier & firstStripId, const unsigned int nStrips, std::vector<Identifier> & clusterVector, + std::vector<std::vector<Identifier> > & idGroups, const SCT_ID& idHelper) const{ + + unsigned int firstStripNumber(idHelper.strip(firstStripId)); + unsigned int firstRowNumber(idHelper.row(firstStripId)); + unsigned int endStripNumber(firstStripNumber + nStrips); // one-past-the-end + Identifier waferId(idHelper.wafer_id(firstStripId)); + unsigned int stripNumber(firstStripNumber); + unsigned int nBadStrips(0); + clusterVector.reserve(clusterVector.size() + nStrips); + const Identifier badId; + for (; stripNumber not_eq endStripNumber; ++stripNumber){ + Identifier stripId(idHelper.strip_id(waferId,firstRowNumber,stripNumber)); + if (isBad(stripId)) { + ++nBadStrips; + stripId = badId; + } + clusterVector.push_back(stripId); + } + + // Maybe all the strips are bad, clear the cluster vector + if (clusterVector.size() == nBadStrips){ + clusterVector.clear(); + // No need to recluster if vector is empty (same true if empty for other reasons) + return; + } + + // Now we have one vector of stripIds, some of which may be 'bad' + if (nBadStrips != 0) { + + clusterVector=recluster(clusterVector, idGroups); + // After this, the cluster vector is either empty or has the last good cluster + } + } + /** * recluster starts with a vector of Ids, some of which may be invalid due to them being bad strips, @@ -228,7 +322,7 @@ namespace InDet{ ATH_MSG_VERBOSE ("SCT_ClusteringTool::clusterize()"); SCT_ClusterCollection* nullResult(0); - if (collection.empty()) return (msg(MSG::DEBUG) << "Empty RDO collection"<< endreq), nullResult; + if (collection.empty()) return (msg(MSG::DEBUG) << "Empty RDO collection"<< endmsg), nullResult; // Make a copy of the collection for sorting (no need to sort if theres only one RDO) std::vector<const SCT_RDORawData*> collectionCopy(collection.begin(),collection.end()); @@ -241,6 +335,8 @@ namespace InDet{ // Vector of clusters to make the cluster collection (most likely equal to collection size) std::vector<IdVec_t> idGroups; idGroups.reserve(collection.size()); + int n01X=0; + int n11X=0; std::vector<uint16_t> tbinGroups; tbinGroups.reserve(collection.size()); @@ -253,19 +349,31 @@ namespace InDet{ Identifier firstStripId(pRawData->identify()); unsigned int nStrips(pRawData->getGroupSize()); int thisStrip(idHelper.strip(firstStripId)); + int BEC(idHelper.barrel_ec(firstStripId)); + int layer(idHelper.layer_disk(firstStripId)); // Flushes the vector every time a non-adjacent strip is found if (not adjacent(thisStrip, previousStrip) and not(currentVector.empty())){ - // Add this group to existing groups (and flush) - idGroups.push_back(currentVector); + if (m_majority01X) { + if (n01X >= n11X){ + idGroups.push_back(currentVector); + } + } else { + // Add this group to existing groups (and flush) + idGroups.push_back(currentVector); + } currentVector.clear(); - tbinGroups.push_back(hitsInThirdTimeBin); - hitsInThirdTimeBin =0; - stripCount = 0; + n01X=0; + n11X=0; + tbinGroups.push_back(hitsInThirdTimeBin); + hitsInThirdTimeBin =0; + stripCount = 0; } // Only use clusters with certain time bit patterns if m_timeBinStr set bool passTiming(true); + bool pass01X(true); + bool passX1X(true); const SCT3_RawData* pRawData3 = dynamic_cast<const SCT3_RawData*>(pRawData); //sroe: coverity 31562 if (!pRawData3) { @@ -276,36 +384,56 @@ namespace InDet{ std::bitset<3> timePattern(static_cast<unsigned long>(timeBin)); if (!m_timeBinStr.empty()) passTiming = testTimeBins(timeBin); + passX1X = testTimeBinsX1X(pRawData3->getTimeBin()); + if(passX1X) pass01X = testTimeBins01X(pRawData3->getTimeBin()); + if (pass01X) n01X++; + if (passX1X && (!pass01X)) n11X++; + if (m_innermostBarrelX1X) { + if ((BEC==0) && (layer==0) && passX1X) passTiming=true; + else passTiming = pass01X; + } else if (m_innertwoBarrelX1X) { + if ((BEC==0) && (layer==0 || layer==1) && passX1X) passTiming=true; + else passTiming = pass01X; + } + // Now we are either (a) pushing more contiguous strips onto an existing vector // or (b) pushing a new set of ids onto an empty vector - if (passTiming) { - if (not m_checkBadChannels){ - addStripsToCluster(firstStripId, nStrips, currentVector, idHelper); // Note this takes the current vector only - if (stripCount < 16) hitsInThirdTimeBin = hitsInThirdTimeBin | (timePattern.test(0) << stripCount); - stripCount++; - } else { - addStripsToClusterWithChecks(firstStripId, nStrips, currentVector,idGroups, idHelper); // This one includes the groups of vectors as well - if (stripCount < 16) hitsInThirdTimeBin = hitsInThirdTimeBin | (timePattern.test(0) << stripCount); - stripCount++; - } + if (passTiming || m_majority01X) { + if(m_useRowInformation){ + addStripsToClusterInclRows(firstStripId, nStrips, currentVector,idGroups, idHelper); // Note this takes the current vector only + if (stripCount < 16) hitsInThirdTimeBin = hitsInThirdTimeBin | (timePattern.test(0) << stripCount); + stripCount++; + } + else if (not m_checkBadChannels){ + addStripsToCluster(firstStripId, nStrips, currentVector, idHelper); // Note this takes the current vector only + if (stripCount < 16) hitsInThirdTimeBin = hitsInThirdTimeBin | (timePattern.test(0) << stripCount); + stripCount++; + } + else { + addStripsToClusterWithChecks(firstStripId, nStrips, currentVector,idGroups, idHelper); // This one includes the groups of vectors as well + if (stripCount < 16) hitsInThirdTimeBin = hitsInThirdTimeBin | (timePattern.test(0) << stripCount); + stripCount++; + } } if (not currentVector.empty()) { - // Gives the last strip number in the cluster - previousStrip = idHelper.strip(currentVector.back()); + // Gives the last strip number in the cluster + previousStrip = idHelper.strip(currentVector.back()); } } // Still need to add this last vector if (not currentVector.empty()) { - idGroups.push_back(currentVector); - tbinGroups.push_back(hitsInThirdTimeBin); - hitsInThirdTimeBin=0; + if ((!m_majority01X) || (n01X >= n11X)){ + idGroups.push_back(currentVector); + tbinGroups.push_back(hitsInThirdTimeBin); + hitsInThirdTimeBin=0; + } } // Find detector element for these digits Identifier elementID(collection.identify()); const InDetDD::SiDetectorElement* element = manager.getDetectorElement(elementID); - if (!element) return msg(MSG::WARNING) << "Element not in the element map, ID = "<< elementID << endreq, nullResult; + if (!element) return msg(MSG::WARNING) << "Element not in the element map, ID = "<< elementID << endmsg, nullResult; const InDetDD::SCT_ModuleSideDesign* design; if (idHelper.is_barrel(elementID)){ @@ -335,8 +463,11 @@ namespace InDet{ const int nStrips(pGroup->size()); if (nStrips == 0) continue; // - DimensionAndPosition clusterDim = clusterDimensions(idHelper.strip(pGroup->front()), idHelper.strip(pGroup->back()), element, idHelper); - Amg::Vector2D localPos(clusterDim.centre.xPhi(),clusterDim.centre.xEta()); + const InDetDD::SiLocalPosition dummyPos(1,0); + DimensionAndPosition clusterDim(dummyPos,1.0);//just initialize with arbitrary values, will be set properly in the next two lines... + if(m_useRowInformation) clusterDim = clusterDimensionsInclRow(idHelper.strip(pGroup->front()), idHelper.strip(pGroup->back()), idHelper.row(pGroup->front()), element, design); + else clusterDim = clusterDimensions(idHelper.strip(pGroup->front()), idHelper.strip(pGroup->back()), element, idHelper); + Amg::Vector2D localPos(clusterDim.centre.xPhi(),clusterDim.centre.xEta()); // Since clusterId is arbitary (it only needs to be unique) just use ID of first strip //const Identifier clusterId = element->identifierOfPosition(clusterDim.centre); const Identifier clusterId(pGroup->front()); @@ -351,6 +482,7 @@ namespace InDet{ SiWidth siWidth(Amg::Vector2D(nStrips,1), Amg::Vector2D(clusterDim.width,stripLength) ); // //sroe: coverity 31564 + SCT_Cluster* cluster = (m_clusterMaker)? (m_clusterMaker->sctCluster(clusterId,localPos,stripGroup,siWidth,element,m_errorStrategy)) : (new SCT_Cluster(clusterId,localPos,stripGroup,siWidth,element,0)); cluster->setHashAndIndex(clusterCollection->identifyHash(), clusterCollection->size()); @@ -361,7 +493,7 @@ namespace InDet{ if (badStripInClusterOnThisModuleSide) cluster->setHitsInThirdTimeBin(0); /// clusters had been split - recalculating HitsInThirdTimeBin too difficult to be worthwhile for this rare corner case.. clusterCollection->push_back(cluster); - + } return clusterCollection; @@ -382,6 +514,23 @@ namespace InDet{ const InDetDD::SiLocalPosition centre((firstStripPos+lastStripPos)/2.0); return SCT_ClusteringTool::DimensionAndPosition(centre, width); } + + SCT_ClusteringTool::DimensionAndPosition + SCT_ClusteringTool::clusterDimensionsInclRow(const int firstStrip, const int lastStrip, const int row, + const InDetDD::SiDetectorElement* pElement, const InDetDD::SCT_ModuleSideDesign* design) const{ + const int nStrips(lastStrip - firstStrip + 1); + const int firstStrip1D = design->strip1Dim(firstStrip,row); + const int lastStrip1D = design->strip1Dim(lastStrip, row); + const double stripPitch = design->stripPitch(); + const InDetDD::SiCellId cell1(firstStrip1D); + const InDetDD::SiCellId cell2(lastStrip1D); + const InDetDD::SiLocalPosition firstStripPos(pElement->rawLocalPositionOfCell(cell1)); + const InDetDD::SiLocalPosition lastStripPos(pElement->rawLocalPositionOfCell(cell2)); + const double width(nStrips*stripPitch); + const InDetDD::SiLocalPosition centre((firstStripPos+lastStripPos)/2.0); + return SCT_ClusteringTool::DimensionAndPosition(centre, width); + } + bool SCT_ClusteringTool::isBad(const Identifier & stripId) const{ return (not m_conditionsSvc->isGood(stripId, InDetConditions::SCT_STRIP)); diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/src/TotPixelClusterSplitter.cxx b/InnerDetector/InDetRecTools/SiClusterizationTool/src/TotPixelClusterSplitter.cxx index f6e91a74263500338d55f3425f277862d3ea06a4..a183e8631a027243dbccb6ba70958ea64bdaf800 100644 --- a/InnerDetector/InDetRecTools/SiClusterizationTool/src/TotPixelClusterSplitter.cxx +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/src/TotPixelClusterSplitter.cxx @@ -14,9 +14,7 @@ #include "PixelConditionsServices/IPixelCalibSvc.h" -InDet::TotPixelClusterSplitter::TotPixelClusterSplitter(const std::string & type, - const std::string & name, - const IInterface * parent) : +InDet::TotPixelClusterSplitter::TotPixelClusterSplitter(const std::string & type, const std::string & name, const IInterface * parent) : AthAlgTool(type, name, parent), m_calibSvc("PixelCalibSvc", name), m_minPixels(3), @@ -30,22 +28,19 @@ InDet::TotPixelClusterSplitter::TotPixelClusterSplitter(const std::string & type InDet::TotPixelClusterSplitter::~TotPixelClusterSplitter() {} -StatusCode InDet::TotPixelClusterSplitter::initialize() -{ +StatusCode InDet::TotPixelClusterSplitter::initialize() { + CHECK(m_calibSvc.retrieve()); return StatusCode::SUCCESS; } -StatusCode InDet::TotPixelClusterSplitter::finalize() -{ -// ATH_MSG_INFO("finalize()"); +StatusCode InDet::TotPixelClusterSplitter::finalize() { return StatusCode::SUCCESS; } -std::vector<InDet::PixelClusterParts> InDet::TotPixelClusterSplitter::splitCluster( - const InDet::PixelCluster & OrigCluster) const -{ +std::vector<InDet::PixelClusterParts> InDet::TotPixelClusterSplitter::splitCluster( const InDet::PixelCluster & OrigCluster) const { + std::vector<InDet::PixelClusterParts> Parts; // ATH_MSG_INFO("splitCluster() called!"); @@ -218,7 +213,7 @@ std::vector<InDet::PixelClusterParts> InDet::TotPixelClusterSplitter::splitClust if (!pixelIDp){ ATH_MSG_ERROR("Could not get PixelID pointer"); } - const PixelID& pixelID = *pixelIDp; + // const PixelID& pixelID = *pixelIDp; for (unsigned int i = 0; i < NumPixels; i++) diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/src/TruthClusterizationFactory.cxx b/InnerDetector/InDetRecTools/SiClusterizationTool/src/TruthClusterizationFactory.cxx new file mode 100644 index 0000000000000000000000000000000000000000..9c40cfee851b8a81a9a9d0b7caa27061465f522f --- /dev/null +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/src/TruthClusterizationFactory.cxx @@ -0,0 +1,161 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +///////////////////////////////////////////////////////////////////////////////////////////////////// +/// Name : TruthClusterizationFactory.cxx +/// Package : SiClusterizationTool +/// Author : Roland Jansky & Felix Cormier +/// Created : April 2016 +/// +/// DESCRIPTION: Emulates NN evaluation from truth (for ITK studies) +/////////////////////////////////////////////////////////////////////////////////////////////////////// + + + +#include <TMath.h> +#include "SiClusterizationTool/TruthClusterizationFactory.h" + +//for position estimate and clustering +#include "InDetIdentifier/PixelID.h" +#include "InDetPrepRawData/PixelCluster.h" + +#include "GaudiKernel/ServiceHandle.h" +#include "GaudiKernel/Incident.h" +#include "InDetSimData/InDetSimDataCollection.h" +#include "InDetSimData/InDetSimData.h" +#include "HepMC/GenParticle.h" + +#include "TrkEventPrimitives/ParamDefs.h" + +#include "TRandom.h" + +namespace InDet { + + TruthClusterizationFactory::TruthClusterizationFactory(const std::string& name, + const std::string& n, const IInterface* p): + AthAlgTool(name, n,p), + m_incidentSvc("IncidentSvc", n), + m_simDataCollectionName("PixelSDO_Map"), + m_simDataCollection(0) + { + // further properties + declareProperty("IncidentService" , m_incidentSvc ); + + declareInterface<TruthClusterizationFactory>(this); + } + +///////////////////////////////////////////////////////////////////////////////////// +/// Destructor - check up memory allocation +/// delete any memory allocation on the heap + + TruthClusterizationFactory::~TruthClusterizationFactory() {} + + StatusCode TruthClusterizationFactory::initialize() { + if (m_incidentSvc.retrieve().isFailure()){ + ATH_MSG_WARNING("Can not retrieve " << m_incidentSvc << ". Exiting."); + return StatusCode::FAILURE; + } + + // register to the incident service + m_incidentSvc->addListener( this, "BeginEvent"); + + msg(MSG::INFO) << "initialize() successful in " << name() << endmsg; + return StatusCode::SUCCESS; + } + + void TruthClusterizationFactory::handle(const Incident& inc) + { + if ( inc.type() == IncidentType::BeginEvent ){ + // record the SDO collection + if (evtStore()->retrieve(m_simDataCollection,m_simDataCollectionName).isFailure()){ + ATH_MSG_WARNING("Could not retrieve the InDetSimDataCollection with name " << m_simDataCollectionName); + m_simDataCollection = 0; + } else { + ATH_MSG_VERBOSE("Successfully retrieved the InDetSimDataCollection with name " << m_simDataCollectionName); + } + } + } + + std::vector<double> TruthClusterizationFactory::estimateNumberOfParticles(const InDet::PixelCluster& pCluster) + { + std::vector<double> probabilities(3,0.); + auto rdos = pCluster.rdoList(); + bool crazycluster(true); + unsigned int nPartContributing = 0; + //Initialize vector for a list of UNIQUE barcodes for the cluster + std::vector<int> barcodes; + //Loop over all elements (pixels/strips) in the cluster + for (auto rdoIter : rdos){ + if (m_simDataCollection){ + auto simDataIter = m_simDataCollection->find(rdoIter); + if (simDataIter != m_simDataCollection->end()){ + crazycluster = false; + // get the SimData and count the individual contributions + auto simData = (simDataIter->second); + //auto simDataDeposits = simData.getdeposits(); + for( auto deposit : simData.getdeposits() ){ + //If deposit exists + if (deposit.first){ + //Now iterate over all barcodes + std::vector<int>::iterator barcodeIterator; + //Look for the barcode of the specific particle depositing energy in the barcodes vector + barcodeIterator = find(barcodes.begin(), barcodes.end(), deposit.first->barcode()); + //If this barcode is not found + if (!(barcodeIterator != barcodes.end())){ + //Add the barcode to the barcodes vector + barcodes.push_back(deposit.first->barcode()); + } + } + else ATH_MSG_WARNING("No deposits found"); + } + //nPartContributing = simDataDeposits.size() > nPartContributing ? simDataDeposits.size() : nPartContributing; + } + } + } + //Barcodes vector is then a list of the total number of UNIQUE + //barcodes in the cluster, each corresponding to a truth particle + nPartContributing = barcodes.size(); + ATH_MSG_VERBOSE("n Part Contributing: " << nPartContributing); + //Set random seed for smearing NN efficiency + gRandom->SetSeed(0); + ATH_MSG_VERBOSE("Smearing TruthClusterizationFactory probability output for TIDE studies"); + //If only 1 truth particles found + if (nPartContributing==1) { + //NN will always return 100% chance of there being only 1 particle + probabilities[0] = 1.0; + } + //If two unique truth particles found in cluster + else if (nPartContributing==2) { + //90% chance NN returns high probability of there being 2 particles + if (gRandom->Uniform(0,1) < 0.9) probabilities[1] = 1.0; + //Other 10% NN returns high probability of there being 1 particle + else probabilities[0] = 1.0; + } + //If greater than 2 unique truth particles in cluster + else if (nPartContributing>2) { + //90% chance NN returns high probability of there being >2 particles + if (gRandom->Uniform(0,1) < 0.9) probabilities[2] = 1.0; + //Other 10% NN returns high probability of there being 1 particle + else probabilities[0] = 1.0; + } + + //if truth collection not found + if(crazycluster) { + std::vector<double> noprobabilities; + return noprobabilities; + } + //Else return probabilities calculated above + + return probabilities; + + } + + std::vector<Amg::Vector2D> TruthClusterizationFactory::estimatePositions(const InDet::PixelCluster& ) + { + ATH_MSG_ERROR("TruthClusterizationFactory::estimatePositions called for ITk ambiguity setup, should never happen! Digital clustering should be run for positions & errors."); + return std::vector<Amg::Vector2D>(2,Amg::Vector2D (2,0.)); + } + +}//end InDet namespace + diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/src/TruthPixelClusterSplitProbTool.cxx b/InnerDetector/InDetRecTools/SiClusterizationTool/src/TruthPixelClusterSplitProbTool.cxx new file mode 100644 index 0000000000000000000000000000000000000000..1e023a9b6c8f3f599aa9098bdad760c431fdd741 --- /dev/null +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/src/TruthPixelClusterSplitProbTool.cxx @@ -0,0 +1,204 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +///////////////////////////////////////////////////////////////////////////////////////////////////// +/// Name : TruthPixelClusterSplitProbTool.cxx +/// Package : SiClusterizationTool +/// Author : Roland Jansky & Felix Cormier +/// Created : April 2016 +/// +/// DESCRIPTION: Compute cluster splitting probabilities +/// (for splitting a cluster into 2 .... N subclusters) +/////////////////////////////////////////////////////////////////////////////////////////////////////// + +#include "SiClusterizationTool/TruthPixelClusterSplitProbTool.h" +#include "InDetRecToolInterfaces/IPixelClusterSplitProbTool.h" +#include "SiClusterizationTool/TruthClusterizationFactory.h" +#include "InDetIdentifier/PixelID.h" +#include "InDetPrepRawData/PixelClusterSplitProb.h" +#include "VxVertex/RecVertex.h" +#include "InDetBeamSpotService/IBeamCondSvc.h" + + + +namespace InDet +{ + + + TruthPixelClusterSplitProbTool::TruthPixelClusterSplitProbTool(const std::string& t, const std::string& n, const IInterface* p) + :AthAlgTool(t,n,p), + m_truthClusterizationFactory("InDet::NnClusterizationFactory/TruthClusterizationFactory", this), + m_iBeamCondSvc("BeamCondSvc",n), + m_useBeamSpotInfo(true) + { + + m_priorMultiplicityContent.push_back(2793337); + m_priorMultiplicityContent.push_back(82056); + m_priorMultiplicityContent.push_back(19944); + + + declareInterface<IPixelClusterSplitProbTool>(this); + + declareProperty("NnClusterizationFactory",m_truthClusterizationFactory); + declareProperty("BeamCondSv",m_iBeamCondSvc); + declareProperty("PriorMultiplicityContent",m_priorMultiplicityContent); + declareProperty("useBeamSpotInfo",m_useBeamSpotInfo); + + + } + + + StatusCode TruthPixelClusterSplitProbTool::initialize() + { + + + if (m_truthClusterizationFactory.retrieve().isFailure()) + { + ATH_MSG_ERROR(" Unable to retrieve "<< m_truthClusterizationFactory ); + return StatusCode::FAILURE; + } + + if (m_iBeamCondSvc.retrieve().isFailure()) + { + ATH_MSG_ERROR( "Could not find BeamCondSvc." ); + return StatusCode::FAILURE; + } + + ATH_MSG_INFO(" Cluster split prob tool initialized successfully "<< m_truthClusterizationFactory ); + return StatusCode::SUCCESS; + } + + + InDet::PixelClusterSplitProb TruthPixelClusterSplitProbTool::splitProbability(const InDet::PixelCluster& origCluster ) const + { + std::vector<double> vectorOfProbs=m_truthClusterizationFactory->estimateNumberOfParticles(origCluster); + + ATH_MSG_VERBOSE(" Got splitProbability, size of vector: " << vectorOfProbs.size() ); + + if (vectorOfProbs.size()==0) + { + std::vector<double> vectorOfSplitProbs; + vectorOfSplitProbs.push_back(-100); + PixelClusterSplitProb clusterSplitProb(vectorOfSplitProbs); + ATH_MSG_VERBOSE(" Returning single split prob equal to -100 " ); + return clusterSplitProb; + } + + + ATH_MSG_VERBOSE( + " P(1): " << vectorOfProbs[0] << + " P(2): " << vectorOfProbs[1] << + " P(>=3): " << vectorOfProbs[2] ); + + + return compileSplitProbability(vectorOfProbs); + } + + InDet::PixelClusterSplitProb TruthPixelClusterSplitProbTool::splitProbability(const InDet::PixelCluster& origCluster, const Trk::TrackParameters& trackParameters ) const + { + //ugly thing to stop unused variable warning + if(&trackParameters) {}; + std::vector<double> vectorOfProbs=m_truthClusterizationFactory->estimateNumberOfParticles(origCluster); + + ATH_MSG_VERBOSE(" Got splitProbability, size of vector: " << vectorOfProbs.size() ); + + if (vectorOfProbs.size()==0) + { + std::vector<double> vectorOfSplitProbs; + vectorOfSplitProbs.push_back(-100); + PixelClusterSplitProb clusterSplitProb(vectorOfSplitProbs); + ATH_MSG_VERBOSE(" Returning single split prob equal to -100 " ); + return clusterSplitProb; + } + + + ATH_MSG_VERBOSE( + " P(1): " << vectorOfProbs[0] << + " P(2): " << vectorOfProbs[1] << + " P(>=3): " << vectorOfProbs[2] ); + + + return compileSplitProbability(vectorOfProbs); + } + + + + InDet::PixelClusterSplitProb TruthPixelClusterSplitProbTool::compileSplitProbability(std::vector<double>& vectorOfProbs ) const + { + + + double sum=0; + + std::vector<double>::iterator begin=vectorOfProbs.begin(); + std::vector<double>::iterator end=vectorOfProbs.end(); + + for (std::vector<double>::iterator iter=begin;iter!=end;++iter) + { + sum+=*iter; + } + + + ATH_MSG_VERBOSE(" Sum of cluster probabilities is: "<<sum); + + std::vector<double> vectorOfSplitProbs; + + for (std::vector<double>::iterator iter=begin;iter!=end;++iter) + { + (*iter)/=sum; + } + + if (m_priorMultiplicityContent.size()<vectorOfProbs.size()) + { + ATH_MSG_ERROR("Prior compatibilities count " << m_priorMultiplicityContent.size() << " is too small: please correct through job properties."); + return InDet::PixelClusterSplitProb(std::vector<double>()); + } + + double psum=0; + int count=0; + for (std::vector<double>::iterator iter=begin;iter!=end;++iter,++count) + { + psum+=(*iter)/m_priorMultiplicityContent[count]; + } + + count=0; + for (std::vector<double>::iterator iter=begin;iter!=end;++iter,++count) + { + (*iter)/=m_priorMultiplicityContent[count]; + (*iter)/=psum; + } + + double sumTest=0; + + for (std::vector<double>::iterator iter=begin;iter!=end;++iter) + { + ATH_MSG_VERBOSE("After update prob is: " << *iter); + sumTest+=*iter; + } + + ATH_MSG_VERBOSE(" Sum of cluster probabilities is: "<<sumTest); + for (std::vector<double>::iterator iter=begin;iter!=end;++iter) + { + if (iter!=begin) + { + vectorOfSplitProbs.push_back(*iter); + } + } + + ATH_MSG_VERBOSE(" normalized P(1->2): " << vectorOfSplitProbs[0] << " P(2->3): " << vectorOfSplitProbs[1] ); + + PixelClusterSplitProb clusterSplitProb(vectorOfSplitProbs); + + ATH_MSG_VERBOSE("SplitProb: " << clusterSplitProb.splitProbability(2) << " -->3 " << clusterSplitProb.splitProbability(3) ); + + return clusterSplitProb; + } + + + + + +}//end namespace + + + diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/src/TruthPixelClusterSplitter.cxx b/InnerDetector/InDetRecTools/SiClusterizationTool/src/TruthPixelClusterSplitter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..5330993c7c0b154c826ad11c3483bc1961dc759c --- /dev/null +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/src/TruthPixelClusterSplitter.cxx @@ -0,0 +1,259 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +///////////////////////////////////////////////////////////////////////////////////////////////////// +/// Name : TruthPixelClusterSplitter.cxx +/// Package : SiClusterizationTool +/// Author : Roland Jansky & Felix Cormier +/// Created : April 2016 +/// +/// DESCRIPTION: Split cluster given the truth information +/////////////////////////////////////////////////////////////////////////////////////////////////////// + +#include "SiClusterizationTool/TruthPixelClusterSplitter.h" +#include "VxVertex/RecVertex.h" +#include "InDetPrepRawData/PixelCluster.h" +#include "InDetBeamSpotService/IBeamCondSvc.h" +#include "SiClusterizationTool/TruthClusterizationFactory.h" +#include "InDetIdentifier/PixelID.h" +#include "InDetReadoutGeometry/SiDetectorElement.h" +#include "GeoPrimitives/GeoPrimitives.h" +#include "EventPrimitives/EventPrimitives.h" + +InDet::TruthPixelClusterSplitter::TruthPixelClusterSplitter(const std::string &type, + const std::string &name, + const IInterface *parent) : + AthAlgTool(type,name,parent), + m_truthClusterizationFactory("InDet::NnClusterizationFactory/TruthClusterizationFactory", this), + m_iBeamCondSvc("BeamCondSvc",name), + m_useBeamSpotInfo(true) +{ + declareInterface<IPixelClusterSplitter>(this); + + declareProperty("NnClusterizationFactory",m_truthClusterizationFactory); + declareProperty("BeamCondSv",m_iBeamCondSvc); + declareProperty("ThresholdSplittingIntoTwoClusters",m_thresholdSplittingIntoTwoClusters=0.95); + declareProperty("ThresholdSplittingIntoThreeClusters",m_thresholdSplittingIntoThreeClusters=0.90); + declareProperty("SplitOnlyOnBLayer",m_splitOnlyOnBLayer=true); + declareProperty("useBeamSpotInfo",m_useBeamSpotInfo); + +} + +InDet::TruthPixelClusterSplitter::~TruthPixelClusterSplitter() +{} + +StatusCode InDet::TruthPixelClusterSplitter::initialize() { + + if (m_truthClusterizationFactory.retrieve().isFailure()) + { + ATH_MSG_ERROR(" Unable to retrieve "<< m_truthClusterizationFactory ); + return StatusCode::FAILURE; + } + + if (m_iBeamCondSvc.retrieve().isFailure()) + { + ATH_MSG_ERROR( "Could not find BeamCondSvc." ); + return StatusCode::FAILURE; + } + + ATH_MSG_INFO(" Cluster splitter initialized successfully "<< m_truthClusterizationFactory ); + return StatusCode::SUCCESS; +} + +StatusCode InDet::TruthPixelClusterSplitter::finalize() { + return StatusCode::SUCCESS; +} + +/* default method which simply splits cluster into 2 */ +std::vector<InDet::PixelClusterParts> InDet::TruthPixelClusterSplitter::splitCluster(const InDet::PixelCluster& origCluster ) const +{ + + //add treatment for b-layer only HERE + + Trk::RecVertex beamposition(m_iBeamCondSvc->beamVtx()); + + const std::vector<Identifier>& rdos = origCluster.rdoList(); + const std::vector<int>& totList = origCluster.totList(); + + //fill lvl1group all with the same value... (not best way but ...) + std::vector<int> lvl1group; + lvl1group.reserve(rdos.size()); + std::vector<Identifier>::const_iterator rdoBegin=rdos.begin(); + std::vector<Identifier>::const_iterator rdoEnd=rdos.end(); + for (std::vector<Identifier>::const_iterator rdoIter=rdoBegin;rdoIter!=rdoEnd;rdoIter++) + { + lvl1group.push_back(origCluster.LVL1A()); + } + + + std::vector<Amg::Vector2D> allLocalPositions; + std::vector<Amg::MatrixX> allErrorMatrix; + + + std::vector<Amg::MatrixX> errorMatrix; + + Amg::Vector3D beamSpotPosition=Amg::Vector3D( + beamposition.position()[0], + beamposition.position()[1], + beamposition.position()[2]); + + if (!m_useBeamSpotInfo) beamSpotPosition=Amg::Vector3D(0,0,0); + + std::vector<Amg::Vector2D> localPosition=m_truthClusterizationFactory->estimatePositions(origCluster); + + + + if (errorMatrix.size()!=2 || localPosition.size()!=2) + { + ATH_MSG_WARNING("Error matrix or local position vector size is not 2, it is:" << errorMatrix.size() << " or " << localPosition.size() << "."); + } + + std::vector<InDet::PixelClusterParts> allMultiPClusters; + + + allMultiPClusters.push_back(PixelClusterParts(rdos,totList,lvl1group,localPosition[0],errorMatrix[0])); + allMultiPClusters.push_back(PixelClusterParts(rdos,totList,lvl1group,localPosition[1],errorMatrix[1])); + + + return allMultiPClusters; + +} + +std::vector<InDet::PixelClusterParts> InDet::TruthPixelClusterSplitter::splitCluster(const InDet::PixelCluster& origCluster, + const InDet::PixelClusterSplitProb& splitProb) const +{ + + + if (m_splitOnlyOnBLayer) + { + const InDetDD::SiDetectorElement* element=origCluster.detectorElement(); + if (element==0) { + ATH_MSG_WARNING("Could not get detector element"); + return std::vector<InDet::PixelClusterParts>(); + } + const AtlasDetectorID* aid = element->getIdHelper(); + if (aid==0) + { + ATH_MSG_WARNING("Could not get ATLASDetectorID"); + return std::vector<InDet::PixelClusterParts>(); + } + const PixelID* pixelIDp=dynamic_cast<const PixelID*>(aid); + if (!pixelIDp) + { + ATH_MSG_WARNING("Could not get PixelID pointer"); + return std::vector<InDet::PixelClusterParts>(); + } + //check if original pixel is on b-layer and if yes continue, otherwise interrupt... + Identifier pixelId = origCluster.identify(); + if (!pixelIDp->is_blayer(pixelId)) + { + //return empty object... + ATH_MSG_VERBOSE(" Cluster not on b-layer. Return empty object-->back to default clustering." );; + + return std::vector<InDet::PixelClusterParts>(); + } + } + + //add treatment for b-layer only HERE + + Trk::RecVertex beamposition(m_iBeamCondSvc->beamVtx()); + + const std::vector<Identifier>& rdos = origCluster.rdoList(); + const std::vector<int>& totList = origCluster.totList(); + + //fill lvl1group all with the same value... (not best way but ...) + std::vector<int> lvl1group; + lvl1group.reserve(rdos.size()); + std::vector<Identifier>::const_iterator rdoBegin=rdos.begin(); + std::vector<Identifier>::const_iterator rdoEnd=rdos.end(); + for (std::vector<Identifier>::const_iterator rdoIter=rdoBegin;rdoIter!=rdoEnd;rdoIter++) + { + lvl1group.push_back(origCluster.LVL1A()); + } + + + if (splitProb.getHighestSplitMultiplicityStored()<3) return std::vector<InDet::PixelClusterParts>(); + + double splitProb2=splitProb.splitProbability(2); + double splitProb3rel=splitProb.splitProbability(3); + + double splitProb3=splitProb3rel/(splitProb3rel+splitProb2); + + ATH_MSG_VERBOSE( " SplitProb -->2 " << splitProb2 << " SplitProb -->3 " << splitProb3 );; + + int nParticles=1; + + if (splitProb2>m_thresholdSplittingIntoTwoClusters) + { + if (splitProb3>m_thresholdSplittingIntoThreeClusters) + { + nParticles=3; + } + else + { + nParticles=2; + } + } + + + ATH_MSG_VERBOSE( " Decided for n. particles: " << nParticles << "." );; + + std::vector<Amg::Vector2D> allLocalPositions; + std::vector<Amg::MatrixX> allErrorMatrix; + + Amg::Vector3D beamSpotPosition=Amg::Vector3D( + beamposition.position()[0], + beamposition.position()[1], + beamposition.position()[2]); + + if (!m_useBeamSpotInfo) beamSpotPosition=Amg::Vector3D(0,0,0); + + std::vector<InDet::PixelClusterParts> allMultiPClusters; + + if (nParticles==1) + { + std::vector<Amg::MatrixX> errorMatrix; + std::vector<Amg::Vector2D> localPosition=m_truthClusterizationFactory->estimatePositions(origCluster); + + if (errorMatrix.size()!=1 || localPosition.size()!=1) + { + ATH_MSG_ERROR("Error matrix or local position vector size is not 1, it is:" << errorMatrix.size() << " or " << localPosition.size() << "."); + } + + allMultiPClusters.push_back(PixelClusterParts(rdos,totList,lvl1group,localPosition[0],errorMatrix[0])); + } + else if (nParticles==2) + { + + std::vector<Amg::MatrixX> errorMatrix; + std::vector<Amg::Vector2D> localPosition=m_truthClusterizationFactory->estimatePositions(origCluster); + + if (errorMatrix.size()!=2 || localPosition.size()!=2) + { + ATH_MSG_ERROR("Error matrix or local position vector size is not 2, it is:" << errorMatrix.size() << " or " << localPosition.size() << "."); + } + + allMultiPClusters.push_back(PixelClusterParts(rdos,totList,lvl1group,localPosition[0],errorMatrix[0])); + allMultiPClusters.push_back(PixelClusterParts(rdos,totList,lvl1group,localPosition[1],errorMatrix[1])); + } + else if (nParticles==3) + { + + std::vector<Amg::MatrixX> errorMatrix; + std::vector<Amg::Vector2D> localPosition=m_truthClusterizationFactory->estimatePositions(origCluster); + + if (errorMatrix.size()!=3 || localPosition.size()!=3) + { + ATH_MSG_ERROR("Error matrix or local position vector size is not 2, it is:" << errorMatrix.size() << " or " << localPosition.size() << "."); + } + + + allMultiPClusters.push_back(PixelClusterParts(rdos,totList,lvl1group,localPosition[0],errorMatrix[0])); + allMultiPClusters.push_back(PixelClusterParts(rdos,totList,lvl1group,localPosition[1],errorMatrix[1])); + allMultiPClusters.push_back(PixelClusterParts(rdos,totList,lvl1group,localPosition[2],errorMatrix[2])); + } + + return allMultiPClusters; + +} diff --git a/InnerDetector/InDetRecTools/SiClusterizationTool/src/components/SiClusterizationTool_entries.cxx b/InnerDetector/InDetRecTools/SiClusterizationTool/src/components/SiClusterizationTool_entries.cxx index 9ee549093c806ffaf0a5a8be1cc7997b7c05bd82..91db630a96ff91159cdbd4fcef5bb04d5c986884 100755 --- a/InnerDetector/InDetRecTools/SiClusterizationTool/src/components/SiClusterizationTool_entries.cxx +++ b/InnerDetector/InDetRecTools/SiClusterizationTool/src/components/SiClusterizationTool_entries.cxx @@ -12,6 +12,9 @@ #include "SiClusterizationTool/NnPixelClusterSplitter.h" #include "SiClusterizationTool/NnClusterizationFactory.h" #include "SiClusterizationTool/NnPixelClusterSplitProbTool.h" +#include "SiClusterizationTool/TruthPixelClusterSplitter.h" +#include "SiClusterizationTool/TruthClusterizationFactory.h" +#include "SiClusterizationTool/TruthPixelClusterSplitProbTool.h" using namespace InDet; @@ -23,6 +26,9 @@ DECLARE_NAMESPACE_TOOL_FACTORY( InDet, TotPixelClusterSplitter ) DECLARE_NAMESPACE_TOOL_FACTORY( InDet, NnPixelClusterSplitter ) DECLARE_NAMESPACE_TOOL_FACTORY( InDet, NnClusterizationFactory ) DECLARE_NAMESPACE_TOOL_FACTORY( InDet, NnPixelClusterSplitProbTool ) +DECLARE_NAMESPACE_TOOL_FACTORY( InDet, TruthPixelClusterSplitter ) +DECLARE_NAMESPACE_TOOL_FACTORY( InDet, TruthClusterizationFactory ) +DECLARE_NAMESPACE_TOOL_FACTORY( InDet, TruthPixelClusterSplitProbTool ) DECLARE_FACTORY_ENTRIES( SiClusterizationTool ){ @@ -34,6 +40,9 @@ DECLARE_FACTORY_ENTRIES( SiClusterizationTool ){ DECLARE_NAMESPACE_TOOL( InDet, NnPixelClusterSplitter ) DECLARE_NAMESPACE_TOOL( InDet, NnClusterizationFactory ) DECLARE_NAMESPACE_TOOL( InDet, NnPixelClusterSplitProbTool ) + DECLARE_NAMESPACE_TOOL( InDet, TruthPixelClusterSplitter ) + DECLARE_NAMESPACE_TOOL( InDet, TruthClusterizationFactory ) + DECLARE_NAMESPACE_TOOL( InDet, TruthPixelClusterSplitProbTool ) }