diff --git a/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellHelpers.h b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellHelpers.h
new file mode 100644
index 0000000000000000000000000000000000000000..24f6fb5040b585885c25ed9499ee9456dd61f749
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellHelpers.h
@@ -0,0 +1,99 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef REC_CALOCELLHELPERS_H
+#define REC_CALOCELLHELPERS_H
+
+#include <cmath>
+#include "CaloEvent/CaloCell.h"
+#include "CaloGeoHelpers/CaloPhiRange.h"
+
+inline double phiMean(double a, double b) { return 0.5*(a+b) + (a*b < 0)*M_PI; }
+
+/** Return true if the cell crossed was crossed by the track in phi **/
+inline bool crossedPhi(const CaloCell& cell, double phi_entrance, double phi_exit) {
+ double mean_phi = phiMean(phi_entrance, phi_exit);
+ double dphi = fabs( CaloPhiRange::diff(phi_entrance, phi_exit) );
+ double phi_min = mean_phi - dphi, phi_max = mean_phi + dphi;
+
+ return (CaloPhiRange::diff(cell.phi() + cell.caloDDE()->dphi()/2., phi_min) > 0 &&
+ CaloPhiRange::diff(phi_max, cell.phi() - cell.caloDDE()->dphi()/2.) > 0);
+}
+
+/** Return the % of path length crossed by the track inside a cell in eta **/
+inline double getPathLengthInEta(const CaloCell& cell, double eta_entrance, double eta_exit) {
+ double etaMin = cell.eta() - 0.5*cell.caloDDE()->deta();
+ double etaMax = cell.eta() + 0.5*cell.caloDDE()->deta();
+ if ( fabs(eta_entrance - eta_exit) < 1e-6 ) // to avoid FPE
+ return eta_entrance > etaMin && eta_entrance < etaMax;
+
+ double etaMinTrack = std::min(eta_entrance, eta_exit);
+ double etaMaxTrack = std::max(eta_entrance, eta_exit);
+ return (std::min(etaMax, etaMaxTrack) - std::max(etaMin, etaMinTrack))/(etaMaxTrack - etaMinTrack);
+}
+
+/** Return the % of path length crossed by the track inside a cell in Z **/
+inline double getPathLengthInZ(double zMin, double zMax, double z_entrance, double z_exit) {
+ if ( fabs(z_entrance - z_exit) < 1e-6 ) // to avoid FPE
+ return z_entrance > zMin && z_entrance < zMax;
+
+ double zMinTrack = std::min(z_entrance, z_exit);
+ double zMaxTrack = std::max(z_entrance, z_exit);
+ return (std::min(zMax, zMaxTrack) - std::max(zMin, zMinTrack))/(zMaxTrack - zMinTrack);
+}
+
+/** Return the % of path length crossed by the track inside a cell in Z **/
+inline double getPathLengthInZ(const CaloCell& cell, double z_entrance, double z_exit) {
+ return getPathLengthInZ(cell.z() - 0.5*cell.caloDDE()->dz(), cell.z() + 0.5*cell.caloDDE()->dz(), z_entrance, z_exit);
+}
+
+// /** Return the % of path length crossed by the track inside a cell in Z for a ladder shaped cell **/
+// inline double getPathLengthInZTileBar1(const CaloCell& cell, double z_entrance, double z_exit) {
+// // Divide the problem in 2 rectangles: the cell contains 6 tile rows, 3 in the bottom part, 3 in the top
+// // Calculate the point where the track crossed the boundary between the bottom and top parts
+// // and determine the path lenght inside each part
+// TileCellDim *cell_dim = m_tileDDM->get_cell_dim(cell.caloDDE()->identify());
+// if (!cell_dim || cell_dim->getNRows() != 6) // should always be the case for this sampling
+// return 0;
+
+// // Get the percentage in depth covered by the bottom part of the cell
+// double R = (cell_dim->getRMax(2) - cell_dim->getRMin(0))/(cell_dim->getRMax(5) - cell_dim->getRMin(0));
+
+// // The point where the track crossed the boundary
+// double z_middle = z_entrance + R*(z_exit - z_entrance);
+
+// // Get the boundaries of the 2 half-cells
+// double zBottom_min = cell_dim->getZMin(0), zBottom_max = cell_dim->getZMax(0);
+// double zTop_min = cell_dim->getZMin(3), zTop_max = cell_dim->getZMax(3);
+// double pathBottom = getPathLengthInZ(zBottom_min, zBottom_max, z_entrance, z_middle);
+
+// // Calculate the path traversed in the bottom and top parts
+// if (zTop_min == zTop_max) { // top part of B9 cell has 0 width. Don't take the mean in this case
+// if ( fabs(z_middle) > fabs(zBottom_max) && pathBottom > 0)
+// pathBottom = 1;
+// return pathBottom;
+// }
+// double pathTop = getPathLengthInZ(zTop_min, zTop_max, z_middle, z_exit);
+// // Take the weighted mean using the depth ratio since the method gives the %
+// return R*pathBottom + (1.-R)*pathTop;
+// }
+
+/** Return the % of the path crossed inside the cell, given the parameters for the extrapolation at entrance and exit of the layer **/
+inline double pathInsideCell(const CaloCell& cell, const Amg::Vector3D& entry, const Amg::Vector3D& exit ) {
+ if (!crossedPhi(cell, entry.phi(), exit.phi()))
+ return 0;
+ double pathCrossed = 0;
+ if (cell.caloDDE()->getSubCalo() != CaloCell_ID::TILE)
+ pathCrossed = getPathLengthInEta(cell, entry.eta(), exit.eta());
+ else if (cell.caloDDE()->getSampling() == CaloCell_ID::TileBar1) // ladder shape cells, divide the problem in 2
+ pathCrossed = 0; //getPathLengthInZTileBar1(cell, entry.z(), exit.z());
+ else
+ pathCrossed = getPathLengthInZ(cell, entry.z(), exit.z());
+ if (pathCrossed <= 0)
+ return 0;
+ return pathCrossed;
+}
+
+
+#endif
diff --git a/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellSelectorLayerdR.h b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellSelectorLayerdR.h
new file mode 100644
index 0000000000000000000000000000000000000000..6fef33a72e79fc729a0b4c9ea5e6aed27fd1c7b5
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellSelectorLayerdR.h
@@ -0,0 +1,43 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// CaloCellSelectorLayerdR.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+#ifndef CALOCELLSELECTORLAYERDR_H
+#define CALOCELLSELECTORLAYERDR_H
+
+#include "RecoToolInterfaces/ICaloCellSelector.h"
+#include "TrkParametersIdentificationHelpers/TrackParametersIdHelper.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+#include "TrkCaloExtension/CaloExtensionHelpers.h"
+
+namespace Trk
+{
+
+ //Recode of original method by TrackInCaloTools to select calo cells around track
+ //get mid point for each calo layer from the entry & exit point
+ //for each layer's mid point select calo cell within dR in the SAME layer
+ //Fast (?) but
+ //1. may miss calo cell at another layer which nevertheless is close to the crossing
+ //2. current implementation will fail in case where there's multiple entry to a calo layer
+
+ class CaloCellSelectorLayerdR : public ICaloCellSelector {
+ public:
+ CaloCellSelectorLayerdR(double coneSize);
+ ~CaloCellSelectorLayerdR ();
+
+ void setConeSize( double coneSize ) { m_coneSize2 = coneSize*coneSize; }
+
+ bool preSelectAction( const Trk::CaloExtension& caloExtension );
+ bool select( const CaloCell& cell ) const; // select or reject the cell
+
+ private:
+ CaloExtensionHelpers::EtaPhiHashLookupVector m_midPoints;
+ double m_coneSize2;
+ TrackParametersIdHelper parsIdHelper;
+ };
+} // end of namespace
+
+#endif
diff --git a/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellSelectorMinPerp.h b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellSelectorMinPerp.h
new file mode 100644
index 0000000000000000000000000000000000000000..74b0ad91d00ab759ddcb81cf311c8c23a45be080
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellSelectorMinPerp.h
@@ -0,0 +1,41 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// CaloCellSelectorMinPerp.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+#ifndef CALOCELLSELECTORLAYERDR_H
+#define CALOCELLSELECTORLAYERDR_H
+
+#include "RecoToolInterfaces/ICaloCellSelector.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+
+
+namespace Trk
+{
+
+ //Using min perperdicular distance from the nearest interpolated track segment as selection criteria
+ //the nearest interpolation point need not be at the same layer as the cell
+ //so the "missing calo cell at another layer" problem in the "classic" method
+ //is avoided.
+ //Slower as need to find nearest interpoaltion point for each calo cell
+
+ class CaloCellSelectorMinPerp : public ICaloCellSelector {
+ public:
+ CaloCellSelectorMinPerp(double coneSize);
+ ~CaloCellSelectorMinPerp () = default;
+
+ bool preSelectAction( const Trk::CaloExtension& caloExtension );
+
+ bool select( const CaloCell& cell )const; // select or reject the cell
+
+ private:
+ const Trk::CaloExtension* m_caloExtension;
+ double m_coneSize;
+ Amg::Vector3D m_meanPos;
+ double m_perp2cut;
+ };
+} // end of namespace
+
+#endif
diff --git a/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellSelectorNearestdR.h b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellSelectorNearestdR.h
new file mode 100644
index 0000000000000000000000000000000000000000..4eee779133cedb022dadf690ad9c7830601daf65
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellSelectorNearestdR.h
@@ -0,0 +1,41 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// CaloCellSelectorNearestdR.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+#ifndef CALOCELLSELECTORLAYERDR_H
+#define CALOCELLSELECTORLAYERDR_H
+
+#include "RecoToolInterfaces/ICaloCellSelector.h"
+
+
+namespace Trk
+{
+ class CaloExtension;
+
+ //Using dR from the nearest interpolation point as selection criteria
+ //the nearest interpolation point need not be at the same layer as the cell
+ //so the "missing calo cell at another layer" problem in the "classic" method
+ //is avoided.
+ //Slower as need to find nearest interpoaltion point for each calo cell
+
+ class CaloCellSelectorNearestdR : public ICaloCellSelector {
+ public:
+ CaloCellSelectorNearestdR(double coneSize);
+ ~CaloCellSelectorNearestdR ();
+
+ bool preSelectAction( const Trk::CaloExtension& caloExtension );
+
+ bool select( const CaloCell& cell )const; // select or reject the cell
+
+ private:
+ const Trk::CaloExtension* m_caloExtension;
+ double m_coneSize2;
+ };
+
+
+} // end of namespace
+
+#endif
diff --git a/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellSelectorRoughdR.h b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellSelectorRoughdR.h
new file mode 100644
index 0000000000000000000000000000000000000000..a6d4528bf45b6dcb1afcd09d68abf861314979c9
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellSelectorRoughdR.h
@@ -0,0 +1,45 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// CaloCellSelectorRoughdR.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+#ifndef CALOCELLSELECTORLAYERDR_H
+#define CALOCELLSELECTORLAYERDR_H
+
+#include "RecoToolInterfaces/ICaloCellSelector.h"
+#include "TrkParametersIdentificationHelpers/TrackParametersIdHelper.h"
+
+#include "TrkCaloExtension/CaloExtensionHelpers.h"
+
+namespace Trk
+{
+ //Rough method to select calo cells around track
+ //mean eta phi was found for all the calorimeter layer crossing points
+ //cells within dR+tolerance around the mean are selected
+ //tolerance is the max deviation of individual cross points from the mean
+ //Fast and easy, but not very valid when the track is highly curved
+
+ class CaloCellSelectorRoughdR : public ICaloCellSelector {
+ public:
+ CaloCellSelectorRoughdR(double coneSize);
+ ~CaloCellSelectorRoughdR ();
+
+ bool preSelectAction( const Trk::CaloExtension& caloExtension );
+
+ bool select( const CaloCell& cell )const; // select or reject the cell
+
+ private:
+ double m_coneSize, m_coneSize2;
+ double m_midEta;
+ double m_midPhi;
+ double m_maxDiff;
+
+ TrackParametersIdHelper parsIdHelper;
+ CaloExtensionHelpers::EntryExitPerLayerVector m_crossPoints;
+ };
+
+} // end of namespace
+
+#endif
diff --git a/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellSelectorUtils.h b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellSelectorUtils.h
new file mode 100644
index 0000000000000000000000000000000000000000..7547e9b5b29464f0f98fc94839d75bee88d676b3
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CaloCellSelectorUtils.h
@@ -0,0 +1,24 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef CALOCELLSELECTORUTILS_H
+#define CALOCELLSELECTORUTILS_H
+
+#include "TrkCaloExtension/CaloExtension.h"
+
+namespace Utils
+{
+ double deltaPhi(double phi1, double phi2);
+
+ double deltaR2(double eta1, double eta2, double phi1, double phi2);
+
+ double deltaR(double eta1, double eta2, double phi1, double phi2);
+
+ void findNearestPoint( const Amg::Vector3D& inputPos, const Trk::CaloExtension* caloExtension,
+ int& nearestIdx,
+ Amg::Vector3D& nearestPos,
+ Amg::Vector3D& nearestMom );
+}
+
+#endif
diff --git a/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CrossedCaloCellHelper.h b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CrossedCaloCellHelper.h
new file mode 100644
index 0000000000000000000000000000000000000000..67c396a3e6a837fbefd8271472b476c333ece863
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/CrossedCaloCellHelper.h
@@ -0,0 +1,42 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef CROSSEDCALOCELLHELPER_H
+#define CROSSEDCALOCELLHELPER_H
+
+#include "CaloUtils/CaloClusterStoreHelper.h"
+#include "CaloEvent/CaloCellContainer.h"
+#include "xAODCaloEvent/CaloClusterContainer.h"
+#include "xAODCaloEvent/CaloCluster.h"
+#include "ParticleCaloExtension/ParticleCellAssociation.h"
+
+namespace Rec {
+
+ class CrossedCaloCellHelper {
+ public:
+
+ static xAOD::CaloCluster* crossedCells( const Rec::ParticleCellAssociation& association,
+ const CaloCellContainer& cellContainer,
+ xAOD::CaloClusterContainer& clusterContainer ) {
+
+ // create cluster
+ xAOD::CaloCluster* cluster = CaloClusterStoreHelper::makeCluster(&clusterContainer,&cellContainer);
+ if( !cluster ){
+ return nullptr;
+ }
+
+ // loop over intersections and add cells to cluster
+ for( auto entry : association.cellIntersections() ){
+ if( !entry.first || !entry.first->caloDDE() ) continue;
+ int index = cellContainer.findIndex(entry.first->caloDDE()->calo_hash());
+ if( index == -1 ) continue;
+ cluster->addCell(index,1.);
+ }
+ return cluster;
+ }
+ };
+
+}
+
+#endif
diff --git a/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/ExtrapolateToCaloTool.h b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/ExtrapolateToCaloTool.h
deleted file mode 100644
index ee4b80bb19b707ac8989c4cb8d68e62fd1faa2b8..0000000000000000000000000000000000000000
--- a/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/ExtrapolateToCaloTool.h
+++ /dev/null
@@ -1,235 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-// ***************************************************************************
-// ExtrapolateToCaloTool
-// ***************************************************************************
-
-#ifndef TRACKTOCALO_EXTRAPOLATETOCALOTOOL_H
-#define TRACKTOCALO_EXTRAPOLATETOCALOTOOL_H
-
-#include "RecoToolInterfaces/IExtrapolateToCaloTool.h"
-
-#include "AthenaBaseComps/AthAlgTool.h"
-#include "GaudiKernel/ToolHandle.h"
-#include "TrkParameters/TrackParameters.h" // typedef
-
-// xAOD
-#include "xAODTracking/NeutralParticle.h" // typedef
-#include "xAODTracking/TrackParticle.h" // typedef
-
-
-class CaloDepthTool;
-class ICaloSurfaceBuilder;
-
-class ImpactInCalo;
-
-namespace Trk {
- class Track;
- class IExtrapolator;
- class Surface;
-}
-
-//<<<<<< INCLUDES >>>>>>
-
-/** @class ExtrapolateToCaloTool
-This class takes a given Trk::Track, extrapolates it to the Calo
- and returns the variables needed to compare with the cluster
- position.
-
- It heavily uses the Tracking Extrapolator class.
-
- The client is asked to select explicitely the CaloCell_ID::CaloSample where he wants to extrapolate.
- The 3 steps are combined :
- - first guess of eta : take it from last point measured for Trk::Track or the input parameters
- - CreateUserSurface(sample,offset,trketa)
- - extrapolate
-
- Other methods exist which return a vector of Trk::TrackParameters at all calorimeter layers.
- The CaloDepthTool is used to decide which sample is hit at a layer.
-
-
- INTERFACE :
-
- Calorimeter : use CaloCell_ID::CaloSample to choose the sample to extrapolate to
-
- Trk::Trk is recommended while running the full reconstruction or on ESD. If so, the
- extrapolation starts from the last measured point.
-
- Trk::TrackParameters is mandatory while running on AOD's and can be used if you want to start
- the extrapolation from the perigee.
-
- @author Sebastian.Fleischmann@cern.ch
- */
-
-class ExtrapolateToCaloTool : public AthAlgTool,
- virtual public IExtrapolateToCaloTool
-{
-public:
-
- // Constructor:
-
- ExtrapolateToCaloTool(const std::string& type,
- const std::string& name,
- const IInterface* parent);
-
- ~ExtrapolateToCaloTool();
- StatusCode initialize();
- StatusCode finalize();
-
-
- /** get Trk::TrackParameters at the calorimeter sample from Trk::Track
- Particle hypothesis of the track is used, if given hypothesis is set to Trk::undefined */
- virtual const Trk::TrackParameters* extrapolate (const Trk::Track& trk,
- const CaloCell_ID::CaloSample sample,
- const double offset,
- const Trk::PropDirection dir,
- const Trk::ParticleHypothesis particle) const;
-
-
- /** get Trk::TrackParameters at the calorimeter sample from Trk::TrackParticleBase.
- If the Trk::ParticleHypothesis is undefined the tool tries to extract the hypothesis
- from the given Trk::TrackParticleBase (if this does not succeed Trk::nonInteracting is used) */
- virtual const Trk::TrackParameters* extrapolate (const Trk::TrackParticleBase& trk,
- const CaloCell_ID::CaloSample sample,
- const double offset,
- const Trk::PropDirection dir,
- const Trk::ParticleHypothesis particle) const;
-
- /** get Trk::TrackParameters at the calorimeter sample from Trk::TrackParameters */
- virtual const Trk::TrackParameters* extrapolate (const Trk::TrackParameters& trkPar,
- const CaloCell_ID::CaloSample sample,
- const double offset,
- Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir) const;
-
- /** get Trk::NeutralParameters at the calorimeter sample from Trk::NeutralParameters */
- virtual const Trk::NeutralParameters* extrapolate (const Trk::NeutralParameters& trkPar,
- const CaloCell_ID::CaloSample sample,
- const double offset,
- Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir) const;
-
- /** get Trk::TrackParameters at the calorimeter sample from xAOD::TrackParticle */
- virtual const Trk::TrackParameters* extrapolate (const xAOD::TrackParticle& trk,
- const CaloCell_ID::CaloSample sample,
- const double offset,
- const Trk::PropDirection dir=Trk::alongMomentum,
- const Trk::ParticleHypothesis particle=Trk::undefined) const;
-
- /** get Trk::NeutralParameters at the calorimeter sample from xAOD::NeutralParticle */
- virtual const Trk::NeutralParameters* extrapolate (const xAOD::NeutralParticle& trk,
- const CaloCell_ID::CaloSample sample,
- const double offset,
- const Trk::PropDirection dir=Trk::alongMomentum ) const;
-
- /** get vector of Trk::TrackParameters
- at the different calorimeter samples from Trk::TrackParameters */
- virtual const DataVector< const Trk::TrackParameters >* getParametersInCalo (const Trk::TrackParameters& parameters,
- Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir) const;
-
-
- /** get vector of Trk::NeutralParameters
- at the different calorimeter samples from Trk::NeutralParameters */
- virtual const DataVector< const Trk::NeutralParameters >* getParametersInCalo (const Trk::NeutralParameters& parameters,
- Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir) const;
-
-
- /** get DataVector of Trk::TrackParameters
- at the different calorimeter samples from Trk::TrackParticleBase. The DataVector owns the Trk::TrackParameters,
- but the DataVector itself needs to be deleted by the client.
- If the Trk::ParticleHypothesis is undefined the tool tries to extract the hypothesis
- from the given Trk::TrackParticleBase (if this does not succeed Trk::nonInteracting is used)*/
- virtual const DataVector< const Trk::TrackParameters >* getParametersInCalo (const Trk::TrackParticleBase& trk,
- const Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir) const;
-
- virtual const DataVector< const Trk::TrackParameters >* getParametersInCalo (const xAOD::TrackParticle& trk,
- const Trk::ParticleHypothesis particle=Trk::undefined,
- const Trk::PropDirection dir=Trk::alongMomentum) const;
-
- /** get DataVector of Trk::TrackParameters
- at the different calorimeter samples from Trk::Track. The DataVector owns the Trk::TrackParameters,
- but the DataVector itself needs to be deleted by the client.
- Particle hypothesis of the track is used, if given hypothesis is set to Trk::undefined */
- virtual const DataVector< const Trk::TrackParameters >* getParametersInCalo (const Trk::Track& trk,
- const Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir) const;
-
- /** get ImpactInCalo from Trk::Track
- Particle hypothesis of the track is used, if given hypothesis is set to Trk::undefined */
- virtual ImpactInCalo* getImpactInCalo ( const Trk::Track& trk,
- const Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir) const;
-
- /** get ImpactInCalo from Trk::TrackParticleBase
- If the Trk::ParticleHypothesis is undefined the tool tries to extract the hypothesis
- from the given Trk::TrackParticleBase (if this does not succeed Trk::nonInteracting is used)*/
- virtual ImpactInCalo* getImpactInCalo ( const Trk::TrackParticleBase& trk,
- const Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir) const;
-
- /** get ImpactInCalo from Trk::TrackParameters*/
- virtual ImpactInCalo* getImpactInCalo (const Trk::TrackParameters& parameters,
- Trk::ParticleHypothesis particleHypo,
- const Trk::PropDirection dir) const;
-
- /** get ImpactInCalo from Trk::NeutralParameters*/
- virtual ImpactInCalo* getImpactInCalo (const Trk::NeutralParameters& parameters,
- Trk::ParticleHypothesis particleHypo,
- const Trk::PropDirection dir) const;
-
- /** get sum of the momenta at the vertex (designed for conversions)
- If reuse=true, try to get it from auxdata first **/
- virtual Amg::Vector3D getMomentumAtVertex(const xAOD::Vertex&, bool reuse=true) const;
-
- /** get straight line intersection in calo given position and momentum vectors **/
- virtual Trk::SurfaceIntersection getIntersectionInCalo(const Amg::Vector3D& position,
- const Amg::Vector3D& momentum, const CaloCell_ID::CaloSample) const;
-
- /** access to the private tool used to define the extrapolation depth, needed to play with delta-eta */
- virtual ToolHandle<CaloDepthTool> getCaloDepth();
-
-private:
-
- ImpactInCalo* getImpactInCaloFromParametersInCalo(
- const DataVector< const Trk::TrackParameters >& parametersInCalo) const;
-
- ImpactInCalo* getImpactInCaloFromParametersInCalo(
- const DataVector< const Trk::NeutralParameters >& parametersInCalo) const;
-
- /** helper function to take last known parameter in calo and extrapolate it stepwise to next samplings, adds parameters to the DataVector */
- void extrapolateStepwiseInCalo ( const Trk::TrackParameters& parameters,
- Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir,
- DataVector< const Trk::TrackParameters >& parameterVector) const;
-
- /** helper function to take last known parameter in calo and extrapolate it stepwise to next samplings, adds parameters to the DataVector */
- void extrapolateStepwiseInCalo ( const Trk::NeutralParameters& parameters,
- Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir,
- DataVector< const Trk::NeutralParameters >& parameterVector) const;
-
- // Defines the surfaces for extrapolation :
- ToolHandle<ICaloSurfaceBuilder> m_calosurf; //!< tool handle to the calo surface builder
- // gives the distance to calo border to decide if barrel or endcap is hit
- ToolHandle<CaloDepthTool> m_calodepth; //!< tool handle to the calo depth tool
-
- // Pre-configured extrapolator
- ToolHandle< Trk::IExtrapolator > m_extrapolator; //!< tool handle to the extrapolator
-
- double m_offsetPresampler; //!< jobOption: offset for presampler
- double m_offsetStrip; //!< jobOption: offset for strip
- double m_offsetMiddle; //!< jobOption: offset for em middle
- double m_offsetBack; //!< jobOption: offset for em back
- double m_offsetTile; //!< jobOption: offset for tile/hec
- bool m_assumeOrderedTracks; //!< jobOption: assume ordered tracks
- bool m_ignoreMaterialEffectsInsideCalo; //!< jobOption: ignore material effects inside calo, if step-wise extrapolation to all layers of calo samplings is performed (used for getParametersInCalo() and getImpactInCalo() )
-
-};
-
-
-#endif // TRACKTOCALO_EXTRAPOLATETOCALOTOOL_H
diff --git a/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/ImpactInCalo.h b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/ImpactInCalo.h
deleted file mode 100755
index 8a31996edf30736581835f2ed1badff04cb5773a..0000000000000000000000000000000000000000
--- a/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/ImpactInCalo.h
+++ /dev/null
@@ -1,175 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#ifndef IMPACTINCALO_H
-#define IMPACTINCALO_H
-/*************************************************************************
- Package:
- File: ImpactInCalo.h
- Description:
-
-*************************************************************************/
-
-//<<<<<< INCLUDES
-
-/**
-
-The ImpactInCalo collection is created by the TrackToCaloAlg algorithm,
-which loops on all Tracks of a given collection and extrapolates them
-to the LAr presampler (_0), strips (_1), middle (_2), back (_3) and to the Tiles.
-
-For each impact point, the cartesian coordinate x,y,z are in the
-official coordinate system, whereas etaCaloLocal and phiCaloLocal are the variables
-to be compared with the CaloCluster's.
-
-The Track DIRECTION at the impact point is also given, to avoid confusion with the
-POSITION of the impact point.
-
-This collection is used by the CBNT_TrackToCalo class.
-
- */
-
-class ImpactInCalo {
-
- public:
- // constructors
- ImpactInCalo() {};
-
- ImpactInCalo( double x_0, double y_0, double z_0,
- double etaCaloLocal_0, double phiCaloLocal_0,
- double trketa_at_0, double trkphi_at_0,
- double x_1, double y_1, double z_1,
- double etaCaloLocal_1, double phiCaloLocal_1,
- double trketa_at_1, double trkphi_at_1,
- double x_2, double y_2, double z_2,
- double etaCaloLocal_2, double phiCaloLocal_2,
- double trketa_at_2, double trkphi_at_2,
- double x_3, double y_3, double z_3,
- double etaCaloLocal_3, double phiCaloLocal_3,
- double trketa_at_3, double trkphi_at_3,
- double x_tile, double y_tile, double z_tile,
- double etaCaloLocal_tile, double phiCaloLocal_tile,
- double trketa_at_tile, double trkphi_at_tile )
- :
- m_x_0(x_0), m_y_0(y_0), m_z_0(z_0),
- m_etaCaloLocal_0(etaCaloLocal_0), m_phiCaloLocal_0(phiCaloLocal_0),
- m_trketa_at_0(trketa_at_0),m_trkphi_at_0(trkphi_at_0),
- m_x_1(x_1), m_y_1(y_1), m_z_1(z_1),
- m_etaCaloLocal_1(etaCaloLocal_1), m_phiCaloLocal_1(phiCaloLocal_1),
- m_trketa_at_1(trketa_at_1),m_trkphi_at_1(trkphi_at_1),
- m_x_2(x_2), m_y_2(y_2), m_z_2(z_2),
- m_etaCaloLocal_2(etaCaloLocal_2), m_phiCaloLocal_2(phiCaloLocal_2),
- m_trketa_at_2(trketa_at_2),m_trkphi_at_2(trkphi_at_2),
- m_x_3(x_3), m_y_3(y_3), m_z_3(z_3),
- m_etaCaloLocal_3(etaCaloLocal_3), m_phiCaloLocal_3(phiCaloLocal_3),
- m_trketa_at_3(trketa_at_3),m_trkphi_at_3(trkphi_at_3),
- m_x_tile(x_tile), m_y_tile(y_tile), m_z_tile(z_tile),
- m_etaCaloLocal_tile(etaCaloLocal_tile), m_phiCaloLocal_tile(phiCaloLocal_tile),
- m_trketa_at_tile(trketa_at_tile),m_trkphi_at_tile(trkphi_at_tile)
- {};
-
- // destructor
- virtual ~ImpactInCalo(){};
- // virtual ~ImpactInCalo(){};
-
- void print() const ;
-
- // gets
-
- // Presampler
- inline double x_0() const { return m_x_0; }
- inline double y_0() const { return m_y_0; }
- inline double z_0() const { return m_z_0; }
- inline double etaCaloLocal_0() const { return m_etaCaloLocal_0; }
- inline double phiCaloLocal_0() const { return m_phiCaloLocal_0; }
- inline double trketa_at_0() const { return m_trketa_at_0; };
- inline double trkphi_at_0() const { return m_trkphi_at_0; };
-
- // Strip
- inline double x_1() const { return m_x_1; }
- inline double y_1() const { return m_y_1; }
- inline double z_1() const { return m_z_1; }
- inline double etaCaloLocal_1() const { return m_etaCaloLocal_1; }
- inline double phiCaloLocal_1() const { return m_phiCaloLocal_1; }
- inline double trketa_at_1() const { return m_trketa_at_1; };
- inline double trkphi_at_1() const { return m_trkphi_at_1; };
-
- // Middle
- inline double x_2() const { return m_x_2; }
- inline double y_2() const { return m_y_2; }
- inline double z_2() const { return m_z_2; }
- inline double etaCaloLocal_2() const { return m_etaCaloLocal_2; }
- inline double phiCaloLocal_2() const { return m_phiCaloLocal_2; }
- inline double trketa_at_2() const { return m_trketa_at_2; };
- inline double trkphi_at_2() const { return m_trkphi_at_2; };
-
- // Back
- inline double x_3() const { return m_x_3; }
- inline double y_3() const { return m_y_3; }
- inline double z_3() const { return m_z_3; }
- inline double etaCaloLocal_3() const { return m_etaCaloLocal_3; }
- inline double phiCaloLocal_3() const { return m_phiCaloLocal_3; }
- inline double trketa_at_3() const { return m_trketa_at_3; };
- inline double trkphi_at_3() const { return m_trkphi_at_3; };
-
- // Tile
- inline double x_tile() const { return m_x_tile; }
- inline double y_tile() const { return m_y_tile; }
- inline double z_tile() const { return m_z_tile; }
- inline double etaCaloLocal_tile() const { return m_etaCaloLocal_tile; }
- inline double phiCaloLocal_tile() const { return m_phiCaloLocal_tile; }
- inline double trketa_at_tile() const { return m_trketa_at_tile; };
- inline double trkphi_at_tile() const { return m_trkphi_at_tile; };
-
- private:
-
- // Presampler
- double m_x_0;
- double m_y_0;
- double m_z_0;
- double m_etaCaloLocal_0;
- double m_phiCaloLocal_0;
- double m_trketa_at_0;
- double m_trkphi_at_0;
-
- // Strip
- double m_x_1;
- double m_y_1;
- double m_z_1;
- double m_etaCaloLocal_1;
- double m_phiCaloLocal_1;
- double m_trketa_at_1;
- double m_trkphi_at_1;
-
- // Middle
- double m_x_2;
- double m_y_2;
- double m_z_2;
- double m_etaCaloLocal_2;
- double m_phiCaloLocal_2;
- double m_trketa_at_2;
- double m_trkphi_at_2;
-
- // Back
- double m_x_3;
- double m_y_3;
- double m_z_3;
- double m_etaCaloLocal_3;
- double m_phiCaloLocal_3;
- double m_trketa_at_3;
- double m_trkphi_at_3;
-
- // Tile
- double m_x_tile;
- double m_y_tile;
- double m_z_tile;
- double m_etaCaloLocal_tile;
- double m_phiCaloLocal_tile;
- double m_trketa_at_tile;
- double m_trkphi_at_tile;
-
-};
-
-
-#endif // IMPACTINCALO_H
diff --git a/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/ImpactInCaloAlg.h b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/ImpactInCaloAlg.h
deleted file mode 100755
index de50a6cd1167f38160258a938ac4de980ad15b65..0000000000000000000000000000000000000000
--- a/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/ImpactInCaloAlg.h
+++ /dev/null
@@ -1,136 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-///////////////////////////////////////////////////////////////////
-// ImpactInCaloAlg.h
-// Header file for class ImpactInCaloAlg
-///////////////////////////////////////////////////////////////////
-// Top algorithm class for extrapolating tracks to Calo
-///////////////////////////////////////////////////////////////////
-
-
-#ifndef IMPACTINCALOALG_H
-#define IMPACTINCALOALG_H
-
-// Base class
-#include "AthenaBaseComps/AthAlgorithm.h"
-#include "GaudiKernel/ToolHandle.h"
-//#include "StoreGate/DataHandle.h"
-
-//#include "TrkTrack/Track.h"
-#include "TrkTrack/TrackCollection.h" // typedef
-
-// #include "Particle/TrackParticleContainer.h"
-// #include "Particle/TrackParticle.h"
-
-//#include "TrackToCalo/ImpactInCalo.h"
-
-//#include "CaloDetDescr/CaloPhiRange.h"
-#include "CaloIdentifier/CaloCell_ID.h"
-
-
-class StoreGateSvc;
-class IExtrapolateToCaloTool;
-class ICaloCoordinateTool;
-class CaloDetDescrManager;
-class CaloCellList;
-class CaloDepthTool;
-class ImpactInCalo;
-
-namespace Trk {
- class Track;
- class ITrackSelectorTool;
-}
-namespace Rec {
- class TrackParticle;
- class TrackParticleContainer;
-}
-
-/**
-@class ImpactInCaloAlg
-This Algorithm is meant to be an example of use of the TrackToCalo tools.
-It is based on the TrackToCaloAlg, but uses the IExtrapolateToCaloTool:
-
- - It loops on a given Track collection ( choosen by jobOption ) and extrapolates all tracks
- to the layers of the calorimeter using the IExtrapolateToCaloTool.
-
- - Impacts are stored in an ImpactInCaloCollection, which can then be used to fill CBNT
-
- - An example of primitive Track to Cluster matching is provided, as well as an example of how
- one can retreive all the CaloCell's around a track extrapolation, using the CaloCellList
- class ( these 2 are off by default )
-
- @author Sebastian.Fleischmann@cern.ch
-*/
-
-class ImpactInCaloAlg : public AthAlgorithm {
-
-public:
-
- // Constructor with parameters:
- ImpactInCaloAlg(const std::string &name,ISvcLocator *pSvcLocator);
-
- ///////////////////////////////////////////////////////////////////
- // Non-const methods:
- ///////////////////////////////////////////////////////////////////
-
- // Basic algorithm methods:
- virtual StatusCode initialize();
- virtual StatusCode execute();
- virtual StatusCode finalize();
-
- /** Loop on Trk::Tracks and create ImpactInCaloCollection */
- StatusCode CreateTrkImpactInCalo();
-
- /** Retreive ImpactInCaloCollection and compare to CaloClusters */
- StatusCode CompareImpactWithCluster();
- /** Debug output of impacts */
- StatusCode PrintImpact();
-
- /** List of cells crossed by Trk::Tracks
- ( neta and nphi are the total window width, e.g .5x5 ) */
- CaloCellList* CellsCrossedByTrack(const Trk::Track& trk ,
- const CaloCell_ID::CaloSample sam,
- int neta, int nphi);
- /** Debug output of cells crossed*/
- StatusCode PrintCellsCrossed();
-
- ///////////////////////////////////////////////////////////////////
- // Private methods:
- ///////////////////////////////////////////////////////////////////
-//private:
-// ImpactInCaloAlg();
-// ImpactInCaloAlg(const ImpactInCaloAlg&);
-// ImpactInCaloAlg &operator=(const ImpactInCaloAlg&);
-
- ///////////////////////////////////////////////////////////////////
- // Private data:
- ///////////////////////////////////////////////////////////////////
-private:
-
- std::string m_TrackName; //!< jobOption: Track collection name
- const TrackCollection *m_tracks; //!< pointer to TrackCollection
- std::string m_TrackParticleName; //!< jobOption: TrackParticle collection name
- const Rec::TrackParticleContainer *m_particle;
- std::string m_trkinput; //!< jobOption: which input type to use
-
- std::string m_cluster_container; //!< jobOption: Cluster container name
- std::string m_cell_container; //!< jobOption: Cell container name
- const CaloCell_ID* m_calo_id;
- const CaloDetDescrManager* m_calo_dd;
-
- std::string m_ImpactInCalosOutputName; //!< jobOption: Output name of ImpactInCaloCollection
- ToolHandle < IExtrapolateToCaloTool > m_toCalo; //!< the ExtrapolateToCaloTool
- ToolHandle < Trk::ITrackSelectorTool> m_trackSelectorTool;//!< handle to ITrackSelectorTool
- //CaloDepthTool* m_calodepth;
- //CaloPhiRange m_phiRange;
-
- bool m_extrapolateInBothDirections; //!< jobOption: extrapolate each track in both directions (along and opposite momentum)?
- bool m_compareImpWithCluster; //!< jobOption: compare impacts with clusters?
- bool m_printCellsCrossed; //!< jobOption: print cells crossed by track?
-
-};
-
-#endif //IMPACTINCALOALG_H
-
diff --git a/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/ImpactInCaloCollection.h b/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/ImpactInCaloCollection.h
deleted file mode 100755
index 2eed0a6c82f68695b622dc8531cae5b2f22b54d0..0000000000000000000000000000000000000000
--- a/Reconstruction/RecoTools/TrackToCalo/TrackToCalo/ImpactInCaloCollection.h
+++ /dev/null
@@ -1,27 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-#ifndef IMPACTINCALOCOLLECTION_H
-#define IMPACTINCALOCOLLECTION_H
-
-#include "DataModel/DataVector.h"
-#include "CLIDSvc/CLASS_DEF.h"
-#include "TrackToCalo/ImpactInCalo.h"
-
-//class ImpactInCalo;
-class ImpactInCaloCollection : public DataVector<ImpactInCalo> {
-public:
- virtual ~ImpactInCaloCollection() {};
-
-
-
- void print() const;
-};
-
-
-
-CLASS_DEF(ImpactInCaloCollection,1237752222,1)
-
-
-#endif // IMPACTINCALOCOLLECTION_H
diff --git a/Reconstruction/RecoTools/TrackToCalo/cmt/requirements b/Reconstruction/RecoTools/TrackToCalo/cmt/requirements
index bfc823f2260841b6a8dbb0f65bd72c1e7651da2e..7e6f900d9dd8b338f599584d4e1a35868716a3ec 100755
--- a/Reconstruction/RecoTools/TrackToCalo/cmt/requirements
+++ b/Reconstruction/RecoTools/TrackToCalo/cmt/requirements
@@ -1,46 +1,43 @@
package TrackToCalo
-author Claire Adam <claire.bourdarios@cern.ch>
+author Niels van Eldik
+author Yat Long Chan
+author Marco van Woerden
use AtlasPolicy AtlasPolicy-*
use GaudiInterface GaudiInterface-* External
-use AthenaBaseComps AthenaBaseComps-* Control
-use CLIDSvc CLIDSvc-* Control
use RecoToolInterfaces RecoToolInterfaces-* Reconstruction/RecoTools
-use CaloIdentifier CaloIdentifier-* Calorimeter
-use DataModel DataModel-* Control
-use TrkTrack TrkTrack-* Tracking/TrkEvent
-use TrkParameters TrkParameters-* Tracking/TrkEvent
-use xAODTracking xAODTracking-* Event/xAOD
-
+use TrkCaloExtension TrkCaloExtension-* Tracking/TrkEvent
+use GeoPrimitives GeoPrimitives-* DetectorDescription
+use TrkParametersIdentificationHelpers TrkParametersIdentificationHelpers-* Tracking/TrkEvent
+use CaloEvent CaloEvent-* Calorimeter
+use CaloUtils CaloUtils-* Calorimeter
+use CaloGeoHelpers CaloGeoHelpers-* Calorimeter
+use ParticleCaloExtension ParticleCaloExtension-* Reconstruction/RecoEvent
+use xAODCaloEvent xAODCaloEvent-* Event/xAOD
private
-use StoreGate StoreGate-* Control
-use TrkParticleBase TrkParticleBase-* Tracking/TrkEvent
+use AthenaBaseComps AthenaBaseComps-* Control
+use xAODTracking xAODTracking-* Event/xAOD
+use xAODTruth xAODTruth-* Event/xAOD
+use xAODMuon xAODMuon-* Event/xAOD
+use TrkTrack TrkTrack-* Tracking/TrkEvent
+use TrkParameters TrkParameters-* Tracking/TrkEvent
use TrkEventPrimitives TrkEventPrimitives-* Tracking/TrkEvent
-use EventPrimitives EventPrimitives-* Event
-use CaloEvent CaloEvent-* Calorimeter
use CaloUtils CaloUtils-* Calorimeter
-use CaloDetDescr CaloDetDescr-* Calorimeter
-use CaloGeoHelpers CaloGeoHelpers-* Calorimeter
-use CaloTrackingGeometry CaloTrackingGeometry-* Calorimeter
+use CaloInterface CaloInterface-* Calorimeter
use TrkExInterfaces TrkExInterfaces-* Tracking/TrkExtrapolation
use TrkToolInterfaces TrkToolInterfaces-* Tracking/TrkTools
use TrkSurfaces TrkSurfaces-* Tracking/TrkDetDescr
-use Particle Particle-* Reconstruction
use AtlasDetDescr AtlasDetDescr-* DetectorDescription
-use MagFieldInterfaces MagFieldInterfaces-* MagneticField
-
+use ParticlesInConeTools ParticlesInConeTools-* Reconstruction/RecoTools
+use FourMomUtils FourMomUtils-* Event
public
-apply_pattern declare_joboptions files="*.py"
-
-apply_pattern declare_python_modules files="*.py"
-
-#apply_pattern dual_use_library files=*.cxx
+apply_pattern dual_use_library files=*.cxx
-library TrackToCalo *.cxx components/*.cxx
-apply_pattern component_library
+#library TrackToCalo *.cxx components/*.cxx
+#apply_pattern component_library
#private
#macro cppdebugflags '$(cppdebugflags_s)'
diff --git a/Reconstruction/RecoTools/TrackToCalo/doc/mainpage.h b/Reconstruction/RecoTools/TrackToCalo/doc/mainpage.h
index 5e04e6c0c4db5db0f2b0fede2ec8a149333949ee..5f70ef7498b78ee990810c4e074319ff3ace49eb 100755
--- a/Reconstruction/RecoTools/TrackToCalo/doc/mainpage.h
+++ b/Reconstruction/RecoTools/TrackToCalo/doc/mainpage.h
@@ -15,14 +15,7 @@ The following classes are defined in this package:
@htmlinclude annotated.html
- - ExtrapolateToCaloTool : Extrapolates Trk::Track, Trk::TrackParameters to the calorimeter samplings using the Trk::IExtrapolator (extended and revised version of ExtrapolTrackToCaloTool)
- - ExtrapolTrackToCaloTool : Extrapolates Trk::Track to the calorimeter samplings
- - TrackToCaloAlg : Example algorithm which uses ExtrapolTrackToCaloTool to create a ImpactInCaloCollection
- - ImpactInCaloAlg : Example algorithm which uses ExtrapolateToCaloTool to create a ImpactInCaloCollection
- - CBNT_TrackToCalo : Writes contents of ImpactInCaloCollection into CBNT
- - CBNTAA_TrackToCalo : Writes contents of ImpactInCaloCollection into CBNTAA
- - ImpactInCalo : data class containing impact positions and directions at the various calorimter samplings
- - ImpactInCaloCollection : collection of ImpactInCalo
+
@section ExtrasTrackToCalo Extra Pages
diff --git a/Reconstruction/RecoTools/TrackToCalo/python/ExtrapolateToCaloToolBase.py b/Reconstruction/RecoTools/TrackToCalo/python/ExtrapolateToCaloToolBase.py
deleted file mode 100644
index 476e780459201f84e986d5205bebb149140255a1..0000000000000000000000000000000000000000
--- a/Reconstruction/RecoTools/TrackToCalo/python/ExtrapolateToCaloToolBase.py
+++ /dev/null
@@ -1,98 +0,0 @@
-# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-
-
-from TrackToCalo.TrackToCaloConf import ExtrapolateToCaloTool
-
-class ExtrapolateToCaloToolBase ( ExtrapolateToCaloTool ) :
-
- def __init__(self):
- mlog = logging.getLogger( name+'(Base)::__init__ dummy for instance copy' )
- mlog.info('entering')
- ExtrapolateToCaloTool.__init__( self) # call base class constructor
- pass
-
- def __init__(self, name="ExtrapolateToCaloToolBaseUNSET", depth="showerdefault", straightLine=False):
-
- ExtrapolateToCaloTool.__init__( self,name) # call base class constructor
-
- from AthenaCommon.Logging import logging
- mlog = logging.getLogger( name+'(Base)::__init__ ' )
- mlog.info("entering")
-
-
- if name=="ExtrapolateToCaloToolBaseUNSET":
- mlog.info ("name should be explicitly set")
- return
-
- # this creates CaloDepthTool and a CaloSurfaceBuilder where the depth is set to the given value:
- from CaloTrackingGeometry.CaloSurfaceBuilderBase import CaloSurfaceBuilderFactory
- theCaloSurfaceBuilder = CaloSurfaceBuilderFactory(depth)
- # set as public tool
- from AthenaCommon.AppMgr import ToolSvc
- ToolSvc += theCaloSurfaceBuilder
- # use the CaloSurfaceBuilder and the CaloDepthTool
- self.CaloSurfaceBuilder = theCaloSurfaceBuilder
- self.CaloDepthTool = theCaloSurfaceBuilder.CaloDepthTool
-
- #if not straightLine:
- #theExtrapolator=AtlasExtrapolator(name = 'TrackToCaloExtrapolator')
- #theExtrapolator.DoCaloDynamic = False # this turns off dynamic calculation of eloss in calorimeters
- ## all left to MaterialEffects/EnergyLossUpdators
-
- #from TrkExTools.TrkExToolsConf import Trk__MaterialEffectsUpdator as MaterialEffectsUpdator
- #AtlasMaterialEffectsUpdator = MaterialEffectsUpdator(name = 'AtlasMaterialEffectsUpdator')
- #ToolSvc += AtlasMaterialEffectsUpdator #default material effects updator
- #NoElossMaterialEffectsUpdator = MaterialEffectsUpdator(name = 'NoElossMaterialEffectsUpdator')
- #NoElossMaterialEffectsUpdator.EnergyLoss = False
- #ToolSvc += NoElossMaterialEffectsUpdator
-
- ##setup MaterialEffectsUpdator arrays
- #MyUpdators = []
- #MyUpdators += [AtlasMaterialEffectsUpdator] # for ID
- #MyUpdators += [NoElossMaterialEffectsUpdator] # for Calo
- ##MyUpdators += [NoElossMaterialEffectsUpdator] # for muon
-
- #MySubUpdators = []
- #MySubUpdators += [AtlasMaterialEffectsUpdator.name()] # for ID
- #MySubUpdators += [NoElossMaterialEffectsUpdator.name()] # for Calo
- #MySubUpdators += [NoElossMaterialEffectsUpdator.name()] # for muon
-
- #theExtrapolator.MaterialEffectsUpdators = MyUpdators
- #theExtrapolator.SubMEUpdators = MySubUpdators
- #else:
- ##theExtrapolator=AtlasExtrapolator("ATLASExtrapolatorStraightLine")
- ##from TrkExSlPropagator.TrkExSlPropagatorConf import Trk__StraightLinePropagator as Propagator
- ##SlPropagator = Propagator(name='TTC_SlPropagator')
- ##ToolSvc += SlPropagator
- ##theExtrapolator.Propagators = [ SlPropagator ]
-
- #theExtrapolator=AtlasExtrapolator("TrackToCaloExtrapolatorStraightLine")
- #from TrkExSlPropagator.TrkExSlPropagatorConf import Trk__StraightLinePropagator as Propagator
- #SlPropagator = Propagator(name='TTC_SlPropagator')
- #ToolSvc += SlPropagator
- #theExtrapolator.Propagators = [ SlPropagator ]
-
- #Get the extrapolator
- from TrkExTools.AtlasExtrapolator import AtlasExtrapolator
-
- theExtrapolator=AtlasExtrapolator()
-
- ## need to add to ToolSvc before putting in ToolHandle
- ToolSvc+=theExtrapolator
- self.Extrapolator=theExtrapolator
-
-
-
-# factory : autogenerates the name
-# usage :
-# from TrackToCalo.ExtrapolateToCaloToolBase import ExtrapolateToCaloToolFactory
-# exToCalo = ExtrapolateToCaloToolFactory(depth="showerdefault")
-# ToolSvc+=exToCalo
-def ExtrapolateToCaloToolFactory( depth="showerdefault", straightLine=False ):
- # build the tool name by appending depth valure
- #if straightLine:
- #straightName="Straight"
- #else:
- #straightName=""
- #return ExtrapolateToCaloToolBase( "ExtrapolateToCaloTool"+depth+straightName, depth=depth,straightLine=straightLine )
- return ExtrapolateToCaloToolBase( "ExtrapolateToCaloTool"+depth, depth=depth )
diff --git a/Reconstruction/RecoTools/TrackToCalo/share/ExtrapolToCaloTool_joboptions.py b/Reconstruction/RecoTools/TrackToCalo/share/ExtrapolToCaloTool_joboptions.py
deleted file mode 100755
index 4c3916e8d94aed643dd3d7ca5c1d5ef1413434bd..0000000000000000000000000000000000000000
--- a/Reconstruction/RecoTools/TrackToCalo/share/ExtrapolToCaloTool_joboptions.py
+++ /dev/null
@@ -1,16 +0,0 @@
-print "WARNING obsolete joboption. Please Reconstruction/RecoTools/TrackToCalo/python/ExtrapolTrackToCalo configurables instead "
-
-
-if not 'TTC_ExtrapolatorInstance' in dir():
- TTC_ExtrapolatorInstance = 'TTC_Extrapolator'
-# from AW
-from TrkConfigurableWrapper.ConfigurableWrapper import ConfigurableWrapper
-
-
-from TrkExTools.AtlasExtrapolator import AtlasExtrapolator
-
-theAtlasExtrapolator=AtlasExtrapolator()
-ToolSvc+=theAtlasExtrapolator
-TTC_Extrapolator = ConfigurableWrapper(theAtlasExtrapolator)
-# default name Trk::Extrapolator/AtlasExtrapolor used by default ( to be changed)
-
diff --git a/Reconstruction/RecoTools/TrackToCalo/share/MuToCaloAlg_jobOptions.py b/Reconstruction/RecoTools/TrackToCalo/share/MuToCaloAlg_jobOptions.py
deleted file mode 100755
index a2ab9b54976b87c0113cf140ceb49827864376a5..0000000000000000000000000000000000000000
--- a/Reconstruction/RecoTools/TrackToCalo/share/MuToCaloAlg_jobOptions.py
+++ /dev/null
@@ -1,31 +0,0 @@
-# JobOptions to extrapolate a track to Calorimeters
-
-# ------ Specific to this alg :
-
-theApp.Dlls += [ "TrackToCalo" ]
-theApp.TopAlg += ["TrackToCaloAlg/MuToCaloAlg"]
-MuToCaloAlg = Algorithm( "MuToCaloAlg" )
-
-MuToCaloAlg.ClusterContainerName = "CaloTopoClusterEM"
-MuToCaloAlg.CaloCellContainerName = "AllCalo"
-MuToCaloAlg.ImpactInCaloContainerName = "MuImpactInCalo"
-
-# If type is anything else than TrackParticleCandidates, will take Tracks
-# MuToCaloAlg.TrackInputType = "TrackParticleCandidates"
-# MuToCaloAlg.TrackParticleName = "TrackParticleCandidate"
-MuToCaloAlg.TrackInputType = "Tracks"
-MuToCaloAlg.TrackName = "ConvertedMBoyTracks"
-
-# ------ More general part to be used as an example : setup Extrapolator instance
-
-TTC_ExtrapolatorInstance = 'MuToCaloExtrapolator'
-include( "TrackToCalo/ExtrapolToCaloTool_joboptions.py" )
-
-# configure private ExtrapolTrackToCaloTool tool
-MuToCaloAlg.ExtrapolTrackToCaloTool.ExtrapolatorName = TTC_Extrapolator.name()
-MuToCaloAlg.ExtrapolTrackToCaloTool.ExtrapolatorInstanceName = TTC_Extrapolator.instance()
-MuToCaloAlg.ExtrapolTrackToCaloTool.CaloDepthTool.DepthChoice= "entrance"
-
-# the type of extrapolation is set automatically. To force it add here :
-# doStraightToCalo=False will use RungeKutta (default for Atlas)
-# doStraightToCalo=True will use StraightLine (default for ctb)
diff --git a/Reconstruction/RecoTools/TrackToCalo/share/TrackToCaloAlg_jobOptions.py b/Reconstruction/RecoTools/TrackToCalo/share/TrackToCaloAlg_jobOptions.py
deleted file mode 100755
index a8cc6b6aa1e8a02e9a6c7f519287be012301b966..0000000000000000000000000000000000000000
--- a/Reconstruction/RecoTools/TrackToCalo/share/TrackToCaloAlg_jobOptions.py
+++ /dev/null
@@ -1,33 +0,0 @@
-# JobOptions to extrapolate a track to Calorimeters
-
-# ------ Specific to this alg :
-
-theApp.Dlls += [ "TrackToCalo" ]
-theApp.TopAlg += ["TrackToCaloAlg"]
-TrackToCaloAlg = Algorithm( "TrackToCaloAlg" )
-TrackToCaloAlg.ClusterContainerName = "CaloTopoClusterEM"
-TrackToCaloAlg.CaloCellContainerName = "AllCalo"
-TrackToCaloAlg.ImpactInCaloContainerName = "ImpactInCalo"
-
-# If type is anything else than TrackParticleCandidates, will take Tracks
-#TrackToCaloAlg.TrackInputType = "Tracks"
-TrackToCaloAlg.TrackInputType = "TrackParticleCandidates"
-TrackToCaloAlg.TrackName = "Tracks"
-TrackToCaloAlg.TrackParticleName = "TrackParticleCandidate"
-
-# ------ More general part to be used as an example : setup Extrapolator instance
-
-TTC_ExtrapolatorInstance = 'TrackToCaloExtrapolator'
-include( "TrackToCalo/ExtrapolToCaloTool_joboptions.py" )
-
-# configure private ExtrapolTrackToCaloTool tool
-TrackToCaloAlg.ExtrapolTrackToCaloTool.ExtrapolatorName = TTC_Extrapolator.name()
-TrackToCaloAlg.ExtrapolTrackToCaloTool.ExtrapolatorInstanceName = TTC_Extrapolator.instance()
-TrackToCaloAlg.ExtrapolTrackToCaloTool.CaloDepthTool.DepthChoice= "entrance"
-
-# the type of extrapolation is set automatically. To force it add here :
-# doStraightToCalo=False will use RungeKutta (default for Atlas)
-# doStraightToCalo=True will use StraightLine (default for ctb)
-
-
-
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/CaloCellSelectorLayerdR.cxx b/Reconstruction/RecoTools/TrackToCalo/src/CaloCellSelectorLayerdR.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0eb73e0173b370f7a0740a659316ba625442fe56
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/src/CaloCellSelectorLayerdR.cxx
@@ -0,0 +1,54 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// CaloCellSelectorLayerdR.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#include "TrackToCalo/CaloCellSelectorUtils.h"
+#include "TrackToCalo/CaloCellSelectorLayerdR.h"
+#include "CaloEvent/CaloCell.h"
+#include "TrkCaloExtension/CaloExtension.h"
+namespace Trk
+{
+
+ CaloCellSelectorLayerdR::CaloCellSelectorLayerdR(double coneSize):
+ m_coneSize2( coneSize*coneSize ){
+ }
+
+ CaloCellSelectorLayerdR::~CaloCellSelectorLayerdR(){}
+
+ bool CaloCellSelectorLayerdR::preSelectAction( const Trk::CaloExtension& caloExtension ){
+
+ //std::cout << "CaloCellSelectorLayerdR p00" << std::endl;
+ if( caloExtension.caloLayerIntersections().empty() ) return false;
+ //std::cout << "CaloCellSelectorLayerdR p01" << std::endl;
+ CaloExtensionHelpers::midPointEtaPhiHashLookupVector( caloExtension, m_midPoints );
+ return true;
+ }
+
+ bool CaloCellSelectorLayerdR::select( const CaloCell& cell )const {
+ // select cell within dR from the midPoint of the same calo layer
+ const CaloDetDescrElement* dde = cell.caloDDE();
+ //std::cout << "CaloCellSelectorLayerdR p1" << std::endl;
+ if(!dde) return false;
+
+ int samplingID = dde->getSampling();
+ //std::cout << "CaloCellSelectorLayerdR p2 " << samplingID << std::endl;
+ //for (auto entry : m_midPoints){
+ // std::cout << std::get<0>(entry) << std::endl;
+ //}
+ if ( !std::get<0>(m_midPoints[samplingID]) ) return false;
+ double dr = Utils::deltaR2( std::get<1>(m_midPoints[samplingID]),dde->eta(),std::get<2>(m_midPoints[samplingID]),dde->phi());
+ if( dr < m_coneSize2){
+ //std::cout << "CaloCellSelectorLayerdR p3" << std::endl;
+ return true;
+ }
+ //std::cout << "CaloCellSelectorLayerdR p4" << std::endl;
+ return false;
+ }
+
+
+} // end of namespace
+
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/CaloCellSelectorMinPerp.cxx b/Reconstruction/RecoTools/TrackToCalo/src/CaloCellSelectorMinPerp.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..eba7154f683f48e01ac9fbbf8e403e976d2b0640
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/src/CaloCellSelectorMinPerp.cxx
@@ -0,0 +1,93 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// CaloCellSelectorMinPerp.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#include "TrackToCalo/CaloCellSelectorUtils.h"
+#include "TrackToCalo/CaloCellSelectorMinPerp.h"
+#include "CaloEvent/CaloCell.h"
+#include "TrkCaloExtension/CaloExtension.h"
+
+
+namespace Trk
+{
+
+ CaloCellSelectorMinPerp::CaloCellSelectorMinPerp(double coneSize ):
+ m_caloExtension(0),
+ m_coneSize( coneSize ),
+ m_perp2cut(0) {
+ }
+
+
+ bool CaloCellSelectorMinPerp::preSelectAction( const Trk::CaloExtension& caloExtension ){
+
+ m_caloExtension = &caloExtension;
+
+ //find a circular cone (defined by cut in perpendicular distance form track i.e. perp2)
+ //with approx area as the dR cone
+ //the dR cone is roughly thought as an elipse
+ if( m_caloExtension->caloLayerIntersections().empty() ) return false;
+
+ Amg::Vector3D pos = m_caloExtension->caloLayerIntersections()[0]->position();
+
+ double a = fabs(pos.eta());
+ double b = exp(-a); // = tan(theta/2)
+ double c = 2.*b/(1.-b*b); // = tan(theta)
+ double r = c/sqrt(1.+c*c); // distance from beam axis
+ double t1 = 2.*atan(exp(-(a-m_coneSize)));
+ double t2 = 2.*atan(exp(-(a+m_coneSize)));
+ double dRad1 = fabs(t1-t2)/2.;
+ //double dRad1 = 2.*b/(1.+b*b) * m_coneSize; //semi-major axis of the dR cone elipse
+ double dRad2 = r * m_coneSize; //semi minor axis of the dR cone elipse
+ m_perp2cut = dRad1*dRad2; // force area of circle = area of elipse, i.e. pi*r*r = pi*dRad1*dRad2
+ return true;
+ }
+
+ bool CaloCellSelectorMinPerp::select( const CaloCell& cell )const {
+
+ if( !m_caloExtension || m_caloExtension->caloLayerIntersections().empty() ) return false;
+
+ const CaloDetDescrElement* dde = cell.caloDDE();
+ if(!dde) return false;
+
+ Amg::Vector3D cellPos(dde->x(), dde->y(), dde->z());
+
+ int nearestIdx;
+ Amg::Vector3D nearestPos,nearestMom;
+ Utils::findNearestPoint( cellPos, m_caloExtension, nearestIdx, nearestPos, nearestMom);
+
+
+ //get the perp2 from track
+ Amg::Vector3D dPos = cellPos-nearestPos;
+ float perp2 = dPos.perp2(nearestMom);
+
+ //get the total track length from IP to point of cloest approach to the cell
+ //scale the perp2 cut with this length
+ float totTrkLen = sqrt(dPos.mag2() - perp2);
+ if (dPos.dot(nearestMom)<0) {totTrkLen = -totTrkLen;}
+
+ Amg::Vector3D pos;
+ Amg::Vector3D oldPos(0.,0.,0.);
+ const std::vector<const Trk::CurvilinearParameters*>& intersections = m_caloExtension->caloLayerIntersections();
+ for (int i=0;i<=nearestIdx;++i){
+ pos = intersections[i]->position();
+ totTrkLen += (pos-oldPos).mag();
+ std::swap(oldPos, pos);
+ }
+
+ //prevent cell at exact opposite of the track from being selected..
+ if (totTrkLen<0) return false;
+
+ if( perp2 < (m_perp2cut*totTrkLen*totTrkLen)){
+ return true;
+ } // IF
+
+ return false;
+
+ }
+
+} // end of namespace
+
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/CaloCellSelectorNearestdR.cxx b/Reconstruction/RecoTools/TrackToCalo/src/CaloCellSelectorNearestdR.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..ee413f527c6443301a0cd7638bf8d1d711e1bb3c
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/src/CaloCellSelectorNearestdR.cxx
@@ -0,0 +1,49 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// CaloCellSelectorNearestdR.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#include "TrackToCalo/CaloCellSelectorUtils.h"
+#include "TrackToCalo/CaloCellSelectorNearestdR.h"
+#include "CaloEvent/CaloCell.h"
+#include "TrkCaloExtension/CaloExtension.h"
+
+namespace Trk
+{
+
+ CaloCellSelectorNearestdR::CaloCellSelectorNearestdR( double coneSize ):
+ m_coneSize2( coneSize*coneSize ){
+ }
+
+ CaloCellSelectorNearestdR::~CaloCellSelectorNearestdR(){}
+
+ bool CaloCellSelectorNearestdR::preSelectAction( const Trk::CaloExtension& caloExtension ) {
+ m_caloExtension = &caloExtension;
+ return true;
+ }
+
+
+ bool CaloCellSelectorNearestdR::select( const CaloCell& cell )const {
+ if( !m_caloExtension ) return false;
+
+ const CaloDetDescrElement* dde = cell.caloDDE();
+ if(!dde) return false;
+
+ Amg::Vector3D cellPos(dde->x(), dde->y(), dde->z());
+
+ int nearestIdx;
+ Amg::Vector3D nearestPos,nearestMom;
+ Utils::findNearestPoint( cellPos, m_caloExtension, nearestIdx, nearestPos, nearestMom);
+
+ if( Utils::deltaR2( nearestPos.eta(),dde->eta(),nearestPos.phi(),dde->phi()) < m_coneSize2){
+ return true;
+ }
+
+ return false;
+ }
+
+} // end of namespace
+
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/CaloCellSelectorRoughdR.cxx b/Reconstruction/RecoTools/TrackToCalo/src/CaloCellSelectorRoughdR.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..e530ddff2f5ea22f1ab538b3af322defb3c56ba7
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/src/CaloCellSelectorRoughdR.cxx
@@ -0,0 +1,106 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// CaloCellSelectorRoughdR.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#include "TrackToCalo/CaloCellSelectorUtils.h"
+#include "TrackToCalo/CaloCellSelectorRoughdR.h"
+#include "CaloEvent/CaloCell.h"
+#include "TrkCaloExtension/CaloExtension.h"
+
+
+namespace Trk
+{
+
+ CaloCellSelectorRoughdR::CaloCellSelectorRoughdR(double coneSize):
+ m_coneSize( coneSize ), m_midEta(0),m_midPhi(0),m_maxDiff(0)
+ {
+ }
+
+ CaloCellSelectorRoughdR::~CaloCellSelectorRoughdR(){}
+
+ bool CaloCellSelectorRoughdR::preSelectAction( const Trk::CaloExtension& caloExtension ){
+
+ if( caloExtension.caloLayerIntersections().empty() ) return false;
+
+ //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ CaloExtensionHelpers::entryExitPerLayerVector( caloExtension, m_crossPoints );
+ // get mean direction of the calo crossings
+ Amg::Vector3D meanPos(0., 0., 0.);
+ for (auto entry : m_crossPoints){
+ int code = std::get<0>(entry);
+ if ( !(code >= 0 && code < 24) ) { continue; } // not a intersection with a calo layer
+ meanPos += std::get<1>(entry).unit();
+ }
+
+ m_midEta = meanPos.eta();
+ m_midPhi = meanPos.phi();
+
+ //get individual crossings max deviation from the mean direction as tolerance
+ m_maxDiff = 0.;
+ for (auto entry : m_crossPoints){
+ int code = std::get<0>(entry);
+ if ( !(code >= 0 && code < 24) ) { continue; } // not a intersection with a calo layer
+ Amg::Vector3D pos = std::get<1>(entry);
+ double rDiff = Utils::deltaR(pos.eta(), m_midEta, pos.phi(), m_midPhi);
+ if (rDiff>m_maxDiff) m_maxDiff = rDiff;
+ }
+ m_maxDiff += m_coneSize;
+ m_maxDiff *= m_maxDiff; //do comparison in dR*dR
+ //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+
+
+ //const std::vector<const CurvilinearParameters*>& intersections = caloExtension.caloLayerIntersections();
+ //int nPts = intersections.size();
+
+ //Amg::Vector3D meanPos(0., 0., 0.);
+ //m_maxDiff = 0.;
+
+ //// get mean direction of the calo crossings
+ //for (int i=0;i<nPts;++i){
+ // //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ // int code = parsIdHelper.caloSample(intersections[i]->cIdentifier());
+ // if (!parsIdHelper.isEntryToVolume(intersections[i]->cIdentifier())) code = -code;
+ // //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ // if ( !(code >= 0 && code < 24) ) { continue; } // not a intersection with a calo layer
+ // meanPos += intersections[i]->position().unit();
+ //}
+
+ //m_midEta = meanPos.eta();
+ //m_midPhi = meanPos.phi();
+
+ ////get individual crossings max deviation from the mean direction as tolerance
+ //for (int i=0;i<nPts;++i){
+ // //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ // int code = parsIdHelper.caloSample(intersections[i]->cIdentifier());
+ // if (!parsIdHelper.isEntryToVolume(intersections[i]->cIdentifier())) code = -code;
+ // //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ // if ( !(code >= 0 && code < 24) ) { continue; } // not a intersection with a calo layer
+ // Amg::Vector3D pos = intersections[i]->position();
+ // double rDiff = Utils::deltaR(pos.eta(), m_midEta, pos.phi(), m_midPhi);
+ // if (rDiff>m_maxDiff) m_maxDiff = rDiff;
+ //}
+
+ ////select cell within requested cone size+tolerence
+ //m_maxDiff += m_coneSize;
+ //m_maxDiff *= m_maxDiff; //do comparison in dR*dR
+
+ return true;
+ }
+
+ bool CaloCellSelectorRoughdR::select( const CaloCell& cell )const {
+ const CaloDetDescrElement* dde = cell.caloDDE();
+ if(!dde) return false;
+ if( Utils::deltaR2( m_midEta,dde->eta(),m_midPhi,dde->phi()) < m_maxDiff){
+ return true;
+ }
+ return false;
+ }
+
+} // end of namespace
+
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/CaloCellSelectorUtils.cxx b/Reconstruction/RecoTools/TrackToCalo/src/CaloCellSelectorUtils.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..56505c177f8234de139f9d4eda3bb66d130af0d7
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/src/CaloCellSelectorUtils.cxx
@@ -0,0 +1,76 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include <cmath>
+#include <algorithm>
+#include "TrackToCalo/CaloCellSelectorUtils.h"
+
+
+namespace Utils
+{
+ double deltaPhi(double phi1, double phi2) {
+ double dPhi=std::fabs(phi1-phi2);
+ if (dPhi>M_PI) dPhi=2*M_PI-dPhi;
+ return dPhi;
+ }
+
+ double deltaR2(double eta1, double eta2, double phi1, double phi2) {
+ double dPhi=deltaPhi(phi1,phi2);
+ double dEta=std::fabs(eta1-eta2);
+ return dEta*dEta + dPhi*dPhi;
+ }
+
+ double deltaR(double eta1, double eta2, double phi1, double phi2) {
+ return std::sqrt(deltaR2(eta1, eta2, phi1, phi2));
+ }
+
+ void findNearestPoint( const Amg::Vector3D& inputPos, const Trk::CaloExtension* caloExtension,
+ int& nearestIdx,
+ Amg::Vector3D& nearestPos,
+ Amg::Vector3D& nearestMom ){
+
+ const std::vector<const Trk::CurvilinearParameters*>& intersections = caloExtension->caloLayerIntersections();
+ int nPts = intersections.size();
+
+ int idxL, idxR, idxMid;
+ Amg::Vector3D pos (0,0,0), mom (0,0,0); //initialization should be(?) unnecessary, just to suppress scary warning msg...
+ Amg::Vector3D posL(0,0,0), momL(0,0,0);
+ Amg::Vector3D posR(0,0,0), momR(0,0,0);
+
+ // find nearest crossing point
+ idxL = 0;
+ idxR = nPts-1;
+ while ((idxR-idxL)>1){
+ idxMid = (idxL+idxR)/2;
+ pos = intersections[idxMid]->position();
+ mom = intersections[idxMid]->momentum();
+
+ if ( (inputPos-pos).dot(mom)>0) {
+ idxL = idxMid;
+ std::swap(posL, pos);
+ std::swap(momL, mom);
+ }
+ else {
+ idxR = idxMid;
+ std::swap(posR, pos);
+ std::swap(momR, mom);
+ }
+ }
+
+ if (idxL==0){
+ posL = intersections[0]->position();
+ momL = intersections[0]->momentum();
+ }
+ if (idxR==(nPts-1)){
+ posR = intersections[nPts-1]->position();
+ momR = intersections[nPts-1]->momentum();
+ }
+
+ float mag2L = (inputPos-posL).mag2();
+ float mag2R = (inputPos-posR).mag2();
+ nearestIdx = (mag2L<mag2R)? idxL : idxR;
+ nearestPos = (mag2L<mag2R)? posL : posR;
+ nearestMom = (mag2L<mag2R)? momL : momR;
+ }
+}
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/ExtrapolateToCaloTool.cxx b/Reconstruction/RecoTools/TrackToCalo/src/ExtrapolateToCaloTool.cxx
deleted file mode 100644
index fe5b7b8693c56b138e4d12ade0ba0767b26c2fb9..0000000000000000000000000000000000000000
--- a/Reconstruction/RecoTools/TrackToCalo/src/ExtrapolateToCaloTool.cxx
+++ /dev/null
@@ -1,969 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-// ***************************************************************************
-// ExtrapolateToCaloTool
-//****************************************************************************
-
-#include "TrackToCalo/ExtrapolateToCaloTool.h"
-#include "TrackToCalo/ImpactInCalo.h"
-
-#include "GaudiKernel/MsgStream.h"
-/*#include "GaudiKernel/Property.h"
-#include "GaudiKernel/IService.h"
-#include "GaudiKernel/IToolSvc.h"
-#include "GaudiKernel/IMessageSvc.h"
-#include "StoreGate/StoreGateSvc.h"*/
-#include "GaudiKernel/ISvcLocator.h"
-#include <vector>
-
-// Stuff needed for the extrapolation :
-//#include "TrkSurfaces/CylinderSurface.h"
-#include "TrkSurfaces/PerigeeSurface.h"
-#include "TrkParameters/TrackParameters.h"
-// #include "TrkNeutralParameters/NeutralParameters.h"
-// #include "TrkNeutralParameters/NeutralPerigee.h"
-#include "TrkTrack/Track.h"
-#include "TrkParticleBase/TrackParticleBase.h"
-#include "TrkEventPrimitives/PropDirection.h"
-#include "TrkEventPrimitives/ParticleHypothesis.h"
-#include "TrkExInterfaces/IExtrapolator.h"
-// CLHEP
-// #include "CLHEP/Units/SystemOfUnits.h"
-// #include "CLHEP/Vector/Rotation.h"
-// #include "CLHEP/Vector/ThreeVector.h"
-
-// Calo specific stuff :
-/*#include "CaloDetDescr/CaloDetDescrManager.h"
-#include "CaloDetDescr/ICaloCoordinateTool.h"*/
-#include "CaloTrackingGeometry/ICaloSurfaceBuilder.h"
-#include "CaloDetDescr/CaloDepthTool.h"
-
-#include "xAODTracking/Vertex.h"
-//#include <cmath>
-
-
-ExtrapolateToCaloTool::ExtrapolateToCaloTool(const std::string& type,
- const std::string& name,
- const IInterface* parent) :
- AthAlgTool(type, name, parent),
- m_calosurf("CaloSurfaceBuilder"),
- m_calodepth("CaloDepthTool"),
- m_extrapolator("Trk::Extrapolator/AtlasExtrapolator"),
- m_offsetPresampler(0.),
- m_offsetStrip(0.),
- m_offsetMiddle(0.),
- m_offsetBack(0.),
- m_offsetTile(0.),
- m_assumeOrderedTracks(true),
- m_ignoreMaterialEffectsInsideCalo(true) {
- declareInterface<IExtrapolateToCaloTool>( this );
-
- declareProperty ("Extrapolator", m_extrapolator, "Extrapolation tool");
- declareProperty ("CaloSurfaceBuilder", m_calosurf, "Calo Surface builder tool");
- declareProperty ("CaloDepthTool", m_calodepth, "Calo Depth tool");
- declareProperty ("OffsetPresampler", m_offsetPresampler, "extrapolation depth (offset) for presampler");
- declareProperty ("OffsetStrip", m_offsetStrip, "extrapolation depth (offset) for strip");
- declareProperty ("OffsetMiddle", m_offsetMiddle, "extrapolation depth (offset) for em middle");
- declareProperty ("OffsetBack", m_offsetBack, "extrapolation depth (offset) for em back");
- declareProperty ("OffsetTile", m_offsetTile, "extrapolation depth (offset) for em tile");
- declareProperty ("OrderedTracks", m_assumeOrderedTracks, "assume ordered tracks");
- declareProperty ("IgnoreMaterialEffectsInsideCalo", m_ignoreMaterialEffectsInsideCalo, "ignore material effects inside calo, if step-wise extrapolation to all layers of calo samplings is performed (used for getParametersInCalo() and getImpactInCalo() )");
-}
-
-ExtrapolateToCaloTool::~ExtrapolateToCaloTool() {}
-
-StatusCode ExtrapolateToCaloTool::initialize() {
- ATH_MSG_INFO("ExtrapolateToCaloTool initialize()" );
-
- // Get Extrapolator from ToolService
- if (m_extrapolator.retrieve().isFailure()) {
- ATH_MSG_FATAL( "Could not retrieve " << m_extrapolator );
- return StatusCode::FAILURE;
- }
- // Get CaloSurfaceBuilder from ToolService
- if (m_calosurf.retrieve().isFailure()) {
- ATH_MSG_FATAL( "Could not retrieve " << m_calosurf );
- return StatusCode::FAILURE;
- }
- // Get CaloDepthTool from ToolService
- if (m_calodepth.retrieve().isFailure()) {
- ATH_MSG_FATAL( "Could not retrieve " << m_calodepth );
- return StatusCode::FAILURE;
- }
- return StatusCode::SUCCESS;
-}
-
-StatusCode ExtrapolateToCaloTool::finalize() {
- return StatusCode::SUCCESS;
-}
-
-const Trk::TrackParameters* ExtrapolateToCaloTool::extrapolate (const Trk::Track& trk,
- const CaloCell_ID::CaloSample sample,
- const double offset,
- const Trk::PropDirection dir,
- const Trk::ParticleHypothesis particle) const {
-
- const Trk::TrackParameters* parametersAtCalo = 0;
- // get the DataVector of TrackParameters
- const DataVector<const Trk::TrackParameters>* parameterVector = trk.trackParameters();
- if (parameterVector) {
- // use last track parameters on track
- DataVector<const Trk::TrackParameters>::const_iterator curParIter = parameterVector->end();
- --curParIter;
- const Trk::TrackParameters* startParameters = (*curParIter);
- // get the destination Surface
- const double etaDir = (dir==Trk::oppositeMomentum ? -1.*(startParameters->eta()) : startParameters->eta());
- const Trk::Surface* destinationSurface = m_calosurf->CreateUserSurface (sample, offset, etaDir);
- if (!destinationSurface)
- return 0;
- // use particle hypothesis of the track, if given hypothesis is set to undefined
- Trk::ParticleHypothesis particleHypo = (particle == Trk::undefined ? trk.info().particleHypothesis() : particle);
- // extrapolate to calorimeter surface
- if (m_assumeOrderedTracks)
- parametersAtCalo = m_extrapolator->extrapolate(*startParameters,
- *destinationSurface,
- dir,
- false,
- particleHypo);
- else
- parametersAtCalo = m_extrapolator->extrapolate(trk,
- *destinationSurface,
- dir,
- false,
- particleHypo);
- delete destinationSurface;
- destinationSurface = 0;
- }
- return parametersAtCalo;
-}
-
-const Trk::TrackParameters* ExtrapolateToCaloTool::extrapolate (const Trk::TrackParticleBase& trk,
- const CaloCell_ID::CaloSample sample,
- const double offset,
- const Trk::PropDirection dir,
- const Trk::ParticleHypothesis particle) const {
-
- // if the TrackParticle contains the pointer to the original track, use this one:
- if (trk.originalTrack()) {
- return extrapolate(*(trk.originalTrack()), sample, offset, dir, particle);
- }
-
- // get the destination Surface
- const double etaDir = (dir==Trk::oppositeMomentum ? -1.*(trk.definingParameters().momentum().eta()) : trk.definingParameters().momentum().eta());
- const Trk::Surface* destinationSurface = m_calosurf->CreateUserSurface (sample, offset, etaDir);
- if (!destinationSurface) return 0;
- // Without the original track we cannot get the particle hypothesis used during the track fit
- // TODO: (will be added to TrackParticle at some point)!
- // use nonInteracting, if no hypothesis given:
- Trk::ParticleHypothesis hypo = (particle == Trk::undefined ? Trk::nonInteracting : particle);
-
- const Trk::TrackParameters* parametersAtCalo = m_extrapolator->extrapolate(trk.definingParameters(),
- *destinationSurface,
- dir,
- false,
- hypo);
- delete destinationSurface;
- return parametersAtCalo;
-}
-
-
-const Trk::TrackParameters* ExtrapolateToCaloTool::extrapolate (const Trk::TrackParameters& trkPar,
- const CaloCell_ID::CaloSample sample,
- const double offset, Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir) const {
- const Trk::TrackParameters* parametersAtCalo = 0;
- // get the destination Surface
- const double etaDir = (dir==Trk::oppositeMomentum ? -1.*(trkPar.eta()) : trkPar.eta());
- const Trk::Surface* destinationSurface = m_calosurf->CreateUserSurface (sample, offset, etaDir);
- if (!destinationSurface)
- return 0;
- // extrapolate to calorimeter surface
- parametersAtCalo = m_extrapolator->extrapolate(trkPar,
- *destinationSurface,
- dir,
- false,
- particle);
- delete destinationSurface;
- destinationSurface = 0;
- return parametersAtCalo;
-}
-
-const Trk::NeutralParameters* ExtrapolateToCaloTool::extrapolate (const Trk::NeutralParameters& parameters,
- const CaloCell_ID::CaloSample sample,
- const double offset,
- Trk::ParticleHypothesis, // TODO: drop!
- const Trk::PropDirection dir) const {
- const Trk::NeutralParameters* parametersAtCalo = 0;
- // get the destination Surface
- const double etaDir = (dir==Trk::oppositeMomentum ? -1.*(parameters.momentum().eta()) : parameters.momentum().eta());
- const Trk::Surface* destinationSurface = m_calosurf->CreateUserSurface (sample, offset, etaDir);
- if (!destinationSurface) return 0;
- // extrapolate to calorimeter surface
- parametersAtCalo = m_extrapolator->extrapolate(parameters,
- *destinationSurface,
- dir,
- false);
- delete destinationSurface;
- destinationSurface = 0;
- return parametersAtCalo;
-}
-
-const Trk::TrackParameters* ExtrapolateToCaloTool::extrapolate (const xAOD::TrackParticle& trk,
- const CaloCell_ID::CaloSample sample,
- const double offset,
- const Trk::PropDirection dir,
- const Trk::ParticleHypothesis particle) const {
- const Trk::TrackParameters* parametersAtCalo = 0;
- // get the destination Surface
- const double etaDir = (dir==Trk::oppositeMomentum ? -1.*(trk.eta()) : trk.eta());
- const Trk::Surface* destinationSurface = m_calosurf->CreateUserSurface (sample, offset, etaDir);
- if (!destinationSurface)
- return 0;
-
- Trk::ParticleHypothesis hypo = (particle == Trk::undefined ? Trk::nonInteracting : particle);
-
- // extrapolate to calorimeter surface
- parametersAtCalo = m_extrapolator->extrapolate(trk,
- *destinationSurface,
- dir,
- false,
- hypo);
- delete destinationSurface;
- destinationSurface = 0;
- return parametersAtCalo;
-}
-
-const Trk::NeutralParameters* ExtrapolateToCaloTool::extrapolate (const xAOD::NeutralParticle& trk,
- const CaloCell_ID::CaloSample sample,
- const double offset,
- const Trk::PropDirection dir ) const {
- const Trk::NeutralParameters* parametersAtCalo = 0;
- // get the destination Surface
- const double etaDir = (dir==Trk::oppositeMomentum ? -1.*(trk.eta()) : trk.eta());
- const Trk::Surface* destinationSurface = m_calosurf->CreateUserSurface (sample, offset, etaDir);
- if (!destinationSurface) return 0;
- // extrapolate to calorimeter surface
- parametersAtCalo = m_extrapolator->extrapolate(trk,
- *destinationSurface,
- dir,
- false);
- delete destinationSurface;
- destinationSurface = 0;
- return parametersAtCalo;
-}
-
-const DataVector< const Trk::TrackParameters >* ExtrapolateToCaloTool::getParametersInCalo (
- const Trk::TrackParameters& parameters,
- Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir) const {
-
- // the return vector
- DataVector< const Trk::TrackParameters >* parameterVector = new DataVector< const Trk::TrackParameters >(IExtrapolateToCaloTool::NCaloLayers, SG::OWN_ELEMENTS);
- // variables to select sampler (barrel or endcap)
- double distbar = 0.; // distance to barrel
- double distec = 0.; // distance to end caps
- CaloCell_ID::CaloSample sample;
- // get eta direction of particle:
- double trketa = parameters.momentum().eta();
-
- // PreSampler :
- // figure out which sampling will be hit
- distbar = m_calodepth->deta(CaloCell_ID::PreSamplerB,trketa);
- distec = m_calodepth->deta(CaloCell_ID::PreSamplerE,trketa);
-
- ATH_MSG_DEBUG(" TrackTo ...PS : for eta= " << trketa << " dist to Barrel =" << distbar
- << " to endcap =" << distec);
- if (distbar < 0 ) sample = CaloCell_ID::PreSamplerB;
- else if (distec < 0 ) sample = CaloCell_ID::PreSamplerE;
- else if ( distbar < distec) sample = CaloCell_ID::PreSamplerB;
- else sample = CaloCell_ID::PreSamplerE;
- ATH_MSG_DEBUG( " => will shoot in sample " << sample );
-
- (*parameterVector)[IExtrapolateToCaloTool::PreSampler] = extrapolate(parameters, sample, m_offsetPresampler, particle, dir);
-
- // do rest of the extrapolation step-wise:
- extrapolateStepwiseInCalo ( parameters,
- particle,
- dir,
- *parameterVector);
-
- return parameterVector;
-}
-
-const DataVector< const Trk::NeutralParameters >* ExtrapolateToCaloTool::getParametersInCalo (
- const Trk::NeutralParameters& parameters,
- Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir) const {
-
- // the return vector
- DataVector< const Trk::NeutralParameters >* parameterVector = new DataVector< const Trk::NeutralParameters >(IExtrapolateToCaloTool::NCaloLayers, SG::OWN_ELEMENTS);
- // variables to select sampler (barrel or endcap)
- double distbar = 0.; // distance to barrel
- double distec = 0.; // distance to end caps
- CaloCell_ID::CaloSample sample;
- // get eta direction of particle:
- double trketa = parameters.momentum().eta();
-
- // PreSampler :
- // figure out which sampling will be hit
- distbar = m_calodepth->deta(CaloCell_ID::PreSamplerB,trketa);
- distec = m_calodepth->deta(CaloCell_ID::PreSamplerE,trketa);
-
- ATH_MSG_DEBUG(" TrackTo ...PS : for eta= " << trketa << " dist to Barrel =" << distbar
- << " to endcap =" << distec);
- if (distbar < 0 ) sample = CaloCell_ID::PreSamplerB;
- else if (distec < 0 ) sample = CaloCell_ID::PreSamplerE;
- else if ( distbar < distec) sample = CaloCell_ID::PreSamplerB;
- else sample = CaloCell_ID::PreSamplerE;
- ATH_MSG_DEBUG( " => will shoot in sample " << sample );
-
- (*parameterVector)[IExtrapolateToCaloTool::PreSampler] = extrapolate(parameters, sample, m_offsetPresampler, particle, dir);
-
- // do rest of the extrapolation step-wise:
- extrapolateStepwiseInCalo ( parameters,
- particle,
- dir,
- *parameterVector);
-
- return parameterVector;
-}
-
-
-ImpactInCalo* ExtrapolateToCaloTool::getImpactInCalo (const Trk::Track& trk,
- const Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir) const {
- const DataVector< const Trk::TrackParameters >* parametersInCalo = getParametersInCalo (trk, particle, dir);
- if (!parametersInCalo) return 0;
- ImpactInCalo* imp = getImpactInCaloFromParametersInCalo(*parametersInCalo);
- // delete extrapolated parameters
- delete parametersInCalo;
- return imp;
-}
-
-ImpactInCalo* ExtrapolateToCaloTool::getImpactInCalo (const Trk::TrackParticleBase& trk,
- const Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir) const {
- const DataVector< const Trk::TrackParameters >* parametersInCalo = getParametersInCalo (trk, particle, dir);
- if (!parametersInCalo) return 0;
- ImpactInCalo* imp = getImpactInCaloFromParametersInCalo(*parametersInCalo);
- // delete extrapolated parameters
- delete parametersInCalo;
- return imp;
-}
-
-ImpactInCalo* ExtrapolateToCaloTool::getImpactInCalo (const Trk::TrackParameters& parameters,
- Trk::ParticleHypothesis particleHypo,
- const Trk::PropDirection dir) const {
-
- const DataVector< const Trk::TrackParameters >* parametersInCalo = getParametersInCalo (parameters, particleHypo, dir);
- if (!parametersInCalo) return 0;
- ImpactInCalo* imp = getImpactInCaloFromParametersInCalo(*parametersInCalo);
- // delete extrapolated parameters
- delete parametersInCalo;
- return imp;
-}
-
-ImpactInCalo* ExtrapolateToCaloTool::getImpactInCalo (const Trk::NeutralParameters& parameters,
- Trk::ParticleHypothesis particleHypo,
- const Trk::PropDirection dir) const {
-
- const DataVector< const Trk::NeutralParameters >* parametersInCalo = getParametersInCalo (parameters, particleHypo, dir);
- if (!parametersInCalo) return 0;
- ImpactInCalo* imp = getImpactInCaloFromParametersInCalo(*parametersInCalo);
- // delete extrapolated parameters
- delete parametersInCalo;
- return imp;
-}
-
-ImpactInCalo* ExtrapolateToCaloTool::getImpactInCaloFromParametersInCalo(const DataVector< const Trk::TrackParameters >& parametersInCalo) const {
- double x_0 = 0.; double y_0 = 0.; double z_0 = 0.;
- double etaCaloLocal_0 = 0.; double phiCaloLocal_0 = 0.;
- double trketa_at_0 = 0.; double trkphi_at_0 = 0.;
- double x_1 = 0.; double y_1 = 0.; double z_1 = 0.;
- double etaCaloLocal_1 = 0.; double phiCaloLocal_1 = 0.;
- double trketa_at_1 = 0.; double trkphi_at_1 = 0.;
- double x_2 = 0.; double y_2 = 0.; double z_2 = 0.;
- double etaCaloLocal_2 = 0.; double phiCaloLocal_2 = 0.;
- double trketa_at_2 = 0.; double trkphi_at_2 = 0.;
- double x_3 = 0.; double y_3 = 0.; double z_3 = 0.;
- double etaCaloLocal_3 = 0.; double phiCaloLocal_3 = 0.;
- double trketa_at_3 = 0.; double trkphi_at_3 = 0.;
- double x_tile = 0.; double y_tile = 0.; double z_tile = 0.;
- double etaCaloLocal_tile = 0.; double phiCaloLocal_tile = 0.;
- double trketa_at_tile = 0.; double trkphi_at_tile = 0.;
-
- if (parametersInCalo[IExtrapolateToCaloTool::PreSampler]) {
- x_0 = parametersInCalo[IExtrapolateToCaloTool::PreSampler]->position().x();
- y_0 = parametersInCalo[IExtrapolateToCaloTool::PreSampler]->position().y();
- z_0 = parametersInCalo[IExtrapolateToCaloTool::PreSampler]->position().z();
- etaCaloLocal_0 = parametersInCalo[IExtrapolateToCaloTool::PreSampler]->position().eta();
- phiCaloLocal_0 = parametersInCalo[IExtrapolateToCaloTool::PreSampler]->position().phi();
- trketa_at_0 = parametersInCalo[IExtrapolateToCaloTool::PreSampler]->momentum().eta();
- trkphi_at_0 = parametersInCalo[IExtrapolateToCaloTool::PreSampler]->momentum().phi();
- }
- if (parametersInCalo[IExtrapolateToCaloTool::Strips]) {
- x_1 = parametersInCalo[IExtrapolateToCaloTool::Strips]->position().x();
- y_1 = parametersInCalo[IExtrapolateToCaloTool::Strips]->position().y();
- z_1 = parametersInCalo[IExtrapolateToCaloTool::Strips]->position().z();
- etaCaloLocal_1 = parametersInCalo[IExtrapolateToCaloTool::Strips]->position().eta();
- phiCaloLocal_1 = parametersInCalo[IExtrapolateToCaloTool::Strips]->position().phi();
- trketa_at_1 = parametersInCalo[IExtrapolateToCaloTool::Strips]->momentum().eta();
- trkphi_at_1 = parametersInCalo[IExtrapolateToCaloTool::Strips]->momentum().phi();
- }
- if (parametersInCalo[IExtrapolateToCaloTool::Middle]) {
- x_2 = parametersInCalo[IExtrapolateToCaloTool::Middle]->position().x();
- y_2 = parametersInCalo[IExtrapolateToCaloTool::Middle]->position().y();
- z_2 = parametersInCalo[IExtrapolateToCaloTool::Middle]->position().z();
- etaCaloLocal_2 = parametersInCalo[IExtrapolateToCaloTool::Middle]->position().eta();
- phiCaloLocal_2 = parametersInCalo[IExtrapolateToCaloTool::Middle]->position().phi();
- trketa_at_2 = parametersInCalo[IExtrapolateToCaloTool::Middle]->momentum().eta();
- trkphi_at_2 = parametersInCalo[IExtrapolateToCaloTool::Middle]->momentum().phi();
- }
- if (parametersInCalo[IExtrapolateToCaloTool::Back]) {
- x_3 = parametersInCalo[IExtrapolateToCaloTool::Back]->position().x();
- y_3 = parametersInCalo[IExtrapolateToCaloTool::Back]->position().y();
- z_3 = parametersInCalo[IExtrapolateToCaloTool::Back]->position().z();
- etaCaloLocal_3 = parametersInCalo[IExtrapolateToCaloTool::Back]->position().eta();
- phiCaloLocal_3 = parametersInCalo[IExtrapolateToCaloTool::Back]->position().phi();
- trketa_at_3 = parametersInCalo[IExtrapolateToCaloTool::Back]->momentum().eta();
- trkphi_at_3 = parametersInCalo[IExtrapolateToCaloTool::Back]->momentum().phi();
- }
- if (parametersInCalo[IExtrapolateToCaloTool::Tile]) {
- x_tile = parametersInCalo[IExtrapolateToCaloTool::Tile]->position().x();
- y_tile = parametersInCalo[IExtrapolateToCaloTool::Tile]->position().y();
- z_tile = parametersInCalo[IExtrapolateToCaloTool::Tile]->position().z();
- etaCaloLocal_tile = parametersInCalo[IExtrapolateToCaloTool::Tile]->position().eta();
- phiCaloLocal_tile = parametersInCalo[IExtrapolateToCaloTool::Tile]->position().phi();
- trketa_at_tile = parametersInCalo[IExtrapolateToCaloTool::Tile]->momentum().eta();
- trkphi_at_tile = parametersInCalo[IExtrapolateToCaloTool::Tile]->momentum().phi();
- }
- return new ImpactInCalo(x_0, y_0, z_0,
- etaCaloLocal_0, phiCaloLocal_0,
- trketa_at_0, trkphi_at_0,
- x_1, y_1, z_1,
- etaCaloLocal_1, phiCaloLocal_1,
- trketa_at_1, trkphi_at_1,
- x_2, y_2, z_2,
- etaCaloLocal_2, phiCaloLocal_2,
- trketa_at_2, trkphi_at_2,
- x_3, y_3, z_3,
- etaCaloLocal_3, phiCaloLocal_3,
- trketa_at_3, trkphi_at_3,
- x_tile, y_tile, z_tile,
- etaCaloLocal_tile, phiCaloLocal_tile,
- trketa_at_tile, trkphi_at_tile );
-
-}
-
-ImpactInCalo* ExtrapolateToCaloTool::getImpactInCaloFromParametersInCalo(const DataVector< const Trk::NeutralParameters >& parametersInCalo) const {
- double x_0 = 0.; double y_0 = 0.; double z_0 = 0.;
- double etaCaloLocal_0 = 0.; double phiCaloLocal_0 = 0.;
- double trketa_at_0 = 0.; double trkphi_at_0 = 0.;
- double x_1 = 0.; double y_1 = 0.; double z_1 = 0.;
- double etaCaloLocal_1 = 0.; double phiCaloLocal_1 = 0.;
- double trketa_at_1 = 0.; double trkphi_at_1 = 0.;
- double x_2 = 0.; double y_2 = 0.; double z_2 = 0.;
- double etaCaloLocal_2 = 0.; double phiCaloLocal_2 = 0.;
- double trketa_at_2 = 0.; double trkphi_at_2 = 0.;
- double x_3 = 0.; double y_3 = 0.; double z_3 = 0.;
- double etaCaloLocal_3 = 0.; double phiCaloLocal_3 = 0.;
- double trketa_at_3 = 0.; double trkphi_at_3 = 0.;
- double x_tile = 0.; double y_tile = 0.; double z_tile = 0.;
- double etaCaloLocal_tile = 0.; double phiCaloLocal_tile = 0.;
- double trketa_at_tile = 0.; double trkphi_at_tile = 0.;
-
- if (parametersInCalo[IExtrapolateToCaloTool::PreSampler]) {
- x_0 = parametersInCalo[IExtrapolateToCaloTool::PreSampler]->position().x();
- y_0 = parametersInCalo[IExtrapolateToCaloTool::PreSampler]->position().y();
- z_0 = parametersInCalo[IExtrapolateToCaloTool::PreSampler]->position().z();
- etaCaloLocal_0 = parametersInCalo[IExtrapolateToCaloTool::PreSampler]->position().eta();
- phiCaloLocal_0 = parametersInCalo[IExtrapolateToCaloTool::PreSampler]->position().phi();
- trketa_at_0 = parametersInCalo[IExtrapolateToCaloTool::PreSampler]->momentum().eta();
- trkphi_at_0 = parametersInCalo[IExtrapolateToCaloTool::PreSampler]->momentum().phi();
- }
- if (parametersInCalo[IExtrapolateToCaloTool::Strips]) {
- x_1 = parametersInCalo[IExtrapolateToCaloTool::Strips]->position().x();
- y_1 = parametersInCalo[IExtrapolateToCaloTool::Strips]->position().y();
- z_1 = parametersInCalo[IExtrapolateToCaloTool::Strips]->position().z();
- etaCaloLocal_1 = parametersInCalo[IExtrapolateToCaloTool::Strips]->position().eta();
- phiCaloLocal_1 = parametersInCalo[IExtrapolateToCaloTool::Strips]->position().phi();
- trketa_at_1 = parametersInCalo[IExtrapolateToCaloTool::Strips]->momentum().eta();
- trkphi_at_1 = parametersInCalo[IExtrapolateToCaloTool::Strips]->momentum().phi();
- }
- if (parametersInCalo[IExtrapolateToCaloTool::Middle]) {
- x_2 = parametersInCalo[IExtrapolateToCaloTool::Middle]->position().x();
- y_2 = parametersInCalo[IExtrapolateToCaloTool::Middle]->position().y();
- z_2 = parametersInCalo[IExtrapolateToCaloTool::Middle]->position().z();
- etaCaloLocal_2 = parametersInCalo[IExtrapolateToCaloTool::Middle]->position().eta();
- phiCaloLocal_2 = parametersInCalo[IExtrapolateToCaloTool::Middle]->position().phi();
- trketa_at_2 = parametersInCalo[IExtrapolateToCaloTool::Middle]->momentum().eta();
- trkphi_at_2 = parametersInCalo[IExtrapolateToCaloTool::Middle]->momentum().phi();
- }
- if (parametersInCalo[IExtrapolateToCaloTool::Back]) {
- x_3 = parametersInCalo[IExtrapolateToCaloTool::Back]->position().x();
- y_3 = parametersInCalo[IExtrapolateToCaloTool::Back]->position().y();
- z_3 = parametersInCalo[IExtrapolateToCaloTool::Back]->position().z();
- etaCaloLocal_3 = parametersInCalo[IExtrapolateToCaloTool::Back]->position().eta();
- phiCaloLocal_3 = parametersInCalo[IExtrapolateToCaloTool::Back]->position().phi();
- trketa_at_3 = parametersInCalo[IExtrapolateToCaloTool::Back]->momentum().eta();
- trkphi_at_3 = parametersInCalo[IExtrapolateToCaloTool::Back]->momentum().phi();
- }
- if (parametersInCalo[IExtrapolateToCaloTool::Tile]) {
- x_tile = parametersInCalo[IExtrapolateToCaloTool::Tile]->position().x();
- y_tile = parametersInCalo[IExtrapolateToCaloTool::Tile]->position().y();
- z_tile = parametersInCalo[IExtrapolateToCaloTool::Tile]->position().z();
- etaCaloLocal_tile = parametersInCalo[IExtrapolateToCaloTool::Tile]->position().eta();
- phiCaloLocal_tile = parametersInCalo[IExtrapolateToCaloTool::Tile]->position().phi();
- trketa_at_tile = parametersInCalo[IExtrapolateToCaloTool::Tile]->momentum().eta();
- trkphi_at_tile = parametersInCalo[IExtrapolateToCaloTool::Tile]->momentum().phi();
- }
- return new ImpactInCalo(x_0, y_0, z_0,
- etaCaloLocal_0, phiCaloLocal_0,
- trketa_at_0, trkphi_at_0,
- x_1, y_1, z_1,
- etaCaloLocal_1, phiCaloLocal_1,
- trketa_at_1, trkphi_at_1,
- x_2, y_2, z_2,
- etaCaloLocal_2, phiCaloLocal_2,
- trketa_at_2, trkphi_at_2,
- x_3, y_3, z_3,
- etaCaloLocal_3, phiCaloLocal_3,
- trketa_at_3, trkphi_at_3,
- x_tile, y_tile, z_tile,
- etaCaloLocal_tile, phiCaloLocal_tile,
- trketa_at_tile, trkphi_at_tile );
-
-}
-
-// ==================================================
-ToolHandle<CaloDepthTool> ExtrapolateToCaloTool::getCaloDepth()
-{
- //std::cout << " m_calodepth = " << m_calodepth << std::endl;
- return m_calodepth;
-}
-
-const DataVector< const Trk::TrackParameters >* ExtrapolateToCaloTool::getParametersInCalo (const Trk::TrackParticleBase& trk,
- const Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir) const {
-
- DataVector< const Trk::TrackParameters >* parametersInCalo = new DataVector< const Trk::TrackParameters >(IExtrapolateToCaloTool::NCaloLayers, SG::OWN_ELEMENTS);
- // variables to select sampler (barrel or endcap)
- double distbar = 0.; // distance to barrel
- double distec = 0.; // distance to end caps
- CaloCell_ID::CaloSample sample;
- // get eta direction of particle:
- double trketa = trk.definingParameters().momentum().eta();
-
- // PreSampler :
- // figure out which sampling will be hit
- distbar = m_calodepth->deta(CaloCell_ID::PreSamplerB,trketa);
- distec = m_calodepth->deta(CaloCell_ID::PreSamplerE,trketa);
-
- ATH_MSG_DEBUG(" TrackTo ...PS : for eta= " << trketa << " dist to Barrel =" << distbar
- << " to endcap =" << distec);
- if (distbar < 0 ) sample = CaloCell_ID::PreSamplerB;
- else if (distec < 0 ) sample = CaloCell_ID::PreSamplerE;
- else if ( distbar < distec) sample = CaloCell_ID::PreSamplerB;
- else sample = CaloCell_ID::PreSamplerE;
- ATH_MSG_DEBUG( " => will shoot in sample " << sample );
-
- (*parametersInCalo)[IExtrapolateToCaloTool::PreSampler] = extrapolate(trk, sample, m_offsetPresampler, dir, particle);
-
- Trk::ParticleHypothesis hypo = particle;
- if (particle == Trk::undefined) {
- if (trk.originalTrack()) {
- hypo = trk.originalTrack()->info().particleHypothesis();
- } else {
- // Without the original track we cannot get the particle hypothesis used during the track fit
- // TODO: (will be added to TrackParticle at some point)!
- // use nonInteracting, if no hypothesis given:
- hypo = Trk::nonInteracting;
- }
- }
- // do rest of the extrapolation step-wise:
- extrapolateStepwiseInCalo ( trk.definingParameters(),
- hypo,
- dir,
- *parametersInCalo);
-
- return parametersInCalo;
-}
-
-const DataVector< const Trk::TrackParameters >* ExtrapolateToCaloTool::getParametersInCalo (const xAOD::TrackParticle& trk,
- const Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir) const {
-
- DataVector< const Trk::TrackParameters >* parametersInCalo = new DataVector< const Trk::TrackParameters >(IExtrapolateToCaloTool::NCaloLayers, SG::OWN_ELEMENTS);
- // variables to select sampler (barrel or endcap)
- double distbar = 0.; // distance to barrel
- double distec = 0.; // distance to end caps
- CaloCell_ID::CaloSample sample;
- // get eta direction of particle:
- double trketa = trk.eta();
-
- // PreSampler :
- // figure out which sampling will be hit
- distbar = m_calodepth->deta(CaloCell_ID::PreSamplerB,trketa);
- distec = m_calodepth->deta(CaloCell_ID::PreSamplerE,trketa);
-
- ATH_MSG_DEBUG(" TrackTo ...PS : for eta= " << trketa << " dist to Barrel =" << distbar
- << " to endcap =" << distec);
- if (distbar < 0 ) sample = CaloCell_ID::PreSamplerB;
- else if (distec < 0 ) sample = CaloCell_ID::PreSamplerE;
- else if ( distbar < distec) sample = CaloCell_ID::PreSamplerB;
- else sample = CaloCell_ID::PreSamplerE;
- ATH_MSG_DEBUG( " => will shoot in sample " << sample );
-
- (*parametersInCalo)[IExtrapolateToCaloTool::PreSampler] = extrapolate(trk, sample, m_offsetPresampler, dir, particle);
-
- Trk::ParticleHypothesis hypo = particle;
- if (particle == Trk::undefined) {
- // not very nice, will hopefully have only a single (forwarded) enum at some point,
- // instead of having xAOD::ParticleHypothesis and Trk::ParticleHypothesis
- hypo = static_cast<Trk::ParticleHypothesis>(trk.particleHypothesis());
- }
- // do rest of the extrapolation step-wise:
- extrapolateStepwiseInCalo ( trk.perigeeParameters(),
- hypo,
- dir,
- *parametersInCalo);
-
- return parametersInCalo;
-}
-
-const DataVector< const Trk::TrackParameters >* ExtrapolateToCaloTool::getParametersInCalo (const Trk::Track& trk,
- const Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir) const {
-
- DataVector< const Trk::TrackParameters >* parametersInCalo = new DataVector< const Trk::TrackParameters >(IExtrapolateToCaloTool::NCaloLayers, SG::OWN_ELEMENTS);
- // variables to select sampler (barrel or endcap)
- double distbar = 0.; // distance to barrel
- double distec = 0.; // distance to end caps
- CaloCell_ID::CaloSample sample;
-
- double trketa = 0;
- const Trk::TrackParameters* firstParam = 0;
- if (trk.perigeeParameters()) {
- trketa = trk.perigeeParameters()->eta();
- firstParam = trk.perigeeParameters();
- } else {
- if (trk.trackParameters()) {
- firstParam = *(trk.trackParameters()->begin());
- trketa = firstParam->momentum().eta();
- } else {
- return parametersInCalo;
- }
- }
-
- // PreSampler :
- // figure out which sampling will be hit
- distbar = m_calodepth->deta(CaloCell_ID::PreSamplerB,trketa);
- distec = m_calodepth->deta(CaloCell_ID::PreSamplerE,trketa);
-
- ATH_MSG_DEBUG(" TrackTo ...PS : for eta= " << trketa << " dist to Barrel =" << distbar
- << " to endcap =" << distec);
- if (distbar < 0 ) sample = CaloCell_ID::PreSamplerB;
- else if (distec < 0 ) sample = CaloCell_ID::PreSamplerE;
- else if ( distbar < distec) sample = CaloCell_ID::PreSamplerB;
- else sample = CaloCell_ID::PreSamplerE;
- ATH_MSG_DEBUG( " => will shoot in sample " << sample );
-
- (*parametersInCalo)[IExtrapolateToCaloTool::PreSampler] = extrapolate(trk, sample, m_offsetPresampler, dir, particle);
-
- Trk::ParticleHypothesis hypo = (particle == Trk::undefined ? trk.info().particleHypothesis() : particle);
- // do rest of the extrapolation step-wise:
- extrapolateStepwiseInCalo ( *firstParam,
- hypo,
- dir,
- *parametersInCalo);
-
- return parametersInCalo;
-}
-
-
-void ExtrapolateToCaloTool::extrapolateStepwiseInCalo (
- const Trk::TrackParameters& parameters,
- Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir,
- DataVector< const Trk::TrackParameters >& parameterVector) const {
- double distbar = 0.; // distance to barrel
- double distec = 0.; // distance to end caps
- CaloCell_ID::CaloSample sample;
-
- // check whether we got the parameters at the first sampling:
- if (parameterVector[IExtrapolateToCaloTool::PreSampler] == 0) return;
- // the last successfully extrapolated parameter in the calo:
- const Trk::TrackParameters* lastParamInCalo = parameterVector[IExtrapolateToCaloTool::PreSampler];
- // take position on first sampling to decide which sampling will be hit next:
- double trketa = lastParamInCalo->position().eta();
- // particle hypo for extrapolation inside calo
- Trk::ParticleHypothesis particleHypoStepWiseInCalo = (m_ignoreMaterialEffectsInsideCalo ? Trk::nonInteracting : particle);
-
- // strip :
- distbar = m_calodepth->deta(CaloCell_ID::EMB1,trketa);
- distec = m_calodepth->deta(CaloCell_ID::EME1,trketa);
- ATH_MSG_DEBUG(" TrackTo ...Strip : for eta= " << trketa << " dist to Barrel =" << distbar
- << " to endcap =" << distec);
- if (distbar < 0 ) sample = CaloCell_ID::EMB1;
- else if (distec < 0 ) sample = CaloCell_ID::EME1;
- else if ( distbar < distec) sample = CaloCell_ID::EMB1;
- else sample = CaloCell_ID::EME1;
- ATH_MSG_DEBUG( " => will shoot in sample " << sample );
-
- parameterVector[IExtrapolateToCaloTool::Strips] = extrapolate(*lastParamInCalo, sample, m_offsetStrip, particleHypoStepWiseInCalo, dir);
-
- // fall back to original parameters if extrapolation did not succeed:
- if (!parameterVector[IExtrapolateToCaloTool::Strips]) {
- parameterVector[IExtrapolateToCaloTool::Strips] = extrapolate(parameters, sample, m_offsetStrip, particle, dir);
- }
- if (parameterVector[IExtrapolateToCaloTool::Strips]) {
- lastParamInCalo = parameterVector[IExtrapolateToCaloTool::Strips];
- }
-
- // EM middle :
- distbar = m_calodepth->deta(CaloCell_ID::EMB2,trketa);
- distec = m_calodepth->deta(CaloCell_ID::EME2,trketa);
- ATH_MSG_DEBUG(" TrackTo ...Middle : for eta= " << trketa << " dist to Barrel =" << distbar
- << " to endcap =" << distec);
- if (distbar < 0 ) sample = CaloCell_ID::EMB2;
- else if (distec < 0 ) sample = CaloCell_ID::EME2;
- else if ( distbar < distec) sample = CaloCell_ID::EMB2;
- else sample = CaloCell_ID::EME2;
- ATH_MSG_DEBUG( " => will shoot in sample " << sample );
- parameterVector[IExtrapolateToCaloTool::Middle] = extrapolate(*lastParamInCalo, sample, m_offsetMiddle, particleHypoStepWiseInCalo, dir);
- if (!parameterVector[IExtrapolateToCaloTool::Middle]) {
- parameterVector[IExtrapolateToCaloTool::Middle] = extrapolate(parameters, sample, m_offsetMiddle, particle, dir);
- }
- if (parameterVector[IExtrapolateToCaloTool::Middle]) {
- lastParamInCalo = parameterVector[IExtrapolateToCaloTool::Middle];
- }
-
- // EM back :
- distbar = m_calodepth->deta(CaloCell_ID::EMB3,trketa);
- distec = m_calodepth->deta(CaloCell_ID::EME3,trketa);
- ATH_MSG_DEBUG(" TrackTo ...Back : for eta= " << trketa << " dist to Barrel =" << distbar
- << " to endcap =" << distec);
- if (distbar < 0 ) sample = CaloCell_ID::EMB3;
- else if (distec < 0 ) sample = CaloCell_ID::EME3;
- else if ( distbar < distec) sample = CaloCell_ID::EMB3;
- else sample = CaloCell_ID::EME3;
- ATH_MSG_DEBUG( " => will shoot in sample " << sample );
- parameterVector[IExtrapolateToCaloTool::Back] = extrapolate(*lastParamInCalo, sample, m_offsetBack, particleHypoStepWiseInCalo, dir);
- if (!parameterVector[IExtrapolateToCaloTool::Back]) {
- parameterVector[IExtrapolateToCaloTool::Back] = extrapolate(parameters, sample, m_offsetBack, particle, dir);
- }
- if (parameterVector[IExtrapolateToCaloTool::Back]) {
- lastParamInCalo = parameterVector[IExtrapolateToCaloTool::Back];
- }
-
- // Tile or HEC0 :
- distbar = m_calodepth->deta(CaloCell_ID::TileBar0,trketa);
- distec = m_calodepth->deta(CaloCell_ID::HEC0,trketa);
- ATH_MSG_DEBUG(" TrackTo ...Tile : for eta= " << trketa << " dist to Barrel =" << distbar
- << " to endcap =" << distec);
- if (distbar < 0 ) sample = CaloCell_ID::TileBar0;
- else if (distec < 0 ) sample = CaloCell_ID::HEC0;
- else if ( distbar > distec && distec < 10. ) sample = CaloCell_ID::HEC0;
- else sample = CaloCell_ID::TileBar0;
- ATH_MSG_DEBUG( " => will shoot in sample " << sample );
- parameterVector[IExtrapolateToCaloTool::Tile] = extrapolate(*lastParamInCalo, sample, m_offsetTile, particleHypoStepWiseInCalo, dir);
- if (!parameterVector[IExtrapolateToCaloTool::Tile]) {
- parameterVector[IExtrapolateToCaloTool::Tile] = extrapolate(parameters, sample, m_offsetTile, particle, dir);
- }
- return;
-}
-
-
-void ExtrapolateToCaloTool::extrapolateStepwiseInCalo (
- const Trk::NeutralParameters& parameters,
- Trk::ParticleHypothesis particle,
- const Trk::PropDirection dir,
- DataVector< const Trk::NeutralParameters >& parameterVector) const {
- double distbar = 0.; // distance to barrel
- double distec = 0.; // distance to end caps
- CaloCell_ID::CaloSample sample;
-
- // check whether we got the parameters at the first sampling:
- if (parameterVector[IExtrapolateToCaloTool::PreSampler] == 0) return;
- // the last successfully extrapolated parameter in the calo:
- const Trk::NeutralParameters* lastParamInCalo = parameterVector[IExtrapolateToCaloTool::PreSampler];
- // take position on first sampling to decide which sampling will be hit next:
- double trketa = lastParamInCalo->position().eta();
- // particle hypo for extrapolation inside calo
- Trk::ParticleHypothesis particleHypoStepWiseInCalo = (m_ignoreMaterialEffectsInsideCalo ? Trk::nonInteracting : particle);
-
- // strip :
- distbar = m_calodepth->deta(CaloCell_ID::EMB1,trketa);
- distec = m_calodepth->deta(CaloCell_ID::EME1,trketa);
- ATH_MSG_DEBUG(" TrackTo ...Strip : for eta= " << trketa << " dist to Barrel =" << distbar
- << " to endcap =" << distec);
- if (distbar < 0 ) sample = CaloCell_ID::EMB1;
- else if (distec < 0 ) sample = CaloCell_ID::EME1;
- else if ( distbar < distec) sample = CaloCell_ID::EMB1;
- else sample = CaloCell_ID::EME1;
- ATH_MSG_DEBUG( " => will shoot in sample " << sample );
-
- parameterVector[IExtrapolateToCaloTool::Strips] = extrapolate(*lastParamInCalo, sample, m_offsetStrip, particleHypoStepWiseInCalo, dir);
-
- // fall back to original parameters if extrapolation did not succeed:
- if (!parameterVector[IExtrapolateToCaloTool::Strips]) {
- parameterVector[IExtrapolateToCaloTool::Strips] = extrapolate(parameters, sample, m_offsetStrip, particle, dir);
- }
- if (parameterVector[IExtrapolateToCaloTool::Strips]) {
- lastParamInCalo = parameterVector[IExtrapolateToCaloTool::Strips];
- }
-
- // EM middle :
- distbar = m_calodepth->deta(CaloCell_ID::EMB2,trketa);
- distec = m_calodepth->deta(CaloCell_ID::EME2,trketa);
- ATH_MSG_DEBUG(" TrackTo ...Middle : for eta= " << trketa << " dist to Barrel =" << distbar
- << " to endcap =" << distec);
- if (distbar < 0 ) sample = CaloCell_ID::EMB2;
- else if (distec < 0 ) sample = CaloCell_ID::EME2;
- else if ( distbar < distec) sample = CaloCell_ID::EMB2;
- else sample = CaloCell_ID::EME2;
- ATH_MSG_DEBUG( " => will shoot in sample " << sample );
- parameterVector[IExtrapolateToCaloTool::Middle] = extrapolate(*lastParamInCalo, sample, m_offsetMiddle, particleHypoStepWiseInCalo, dir);
- if (!parameterVector[IExtrapolateToCaloTool::Middle]) {
- parameterVector[IExtrapolateToCaloTool::Middle] = extrapolate(parameters, sample, m_offsetMiddle, particle, dir);
- }
- if (parameterVector[IExtrapolateToCaloTool::Middle]) {
- lastParamInCalo = parameterVector[IExtrapolateToCaloTool::Middle];
- }
-
- // EM back :
- distbar = m_calodepth->deta(CaloCell_ID::EMB3,trketa);
- distec = m_calodepth->deta(CaloCell_ID::EME3,trketa);
- ATH_MSG_DEBUG(" TrackTo ...Back : for eta= " << trketa << " dist to Barrel =" << distbar
- << " to endcap =" << distec);
- if (distbar < 0 ) sample = CaloCell_ID::EMB3;
- else if (distec < 0 ) sample = CaloCell_ID::EME3;
- else if ( distbar < distec) sample = CaloCell_ID::EMB3;
- else sample = CaloCell_ID::EME3;
- ATH_MSG_DEBUG( " => will shoot in sample " << sample );
- parameterVector[IExtrapolateToCaloTool::Back] = extrapolate(*lastParamInCalo, sample, m_offsetBack, particleHypoStepWiseInCalo, dir);
- if (!parameterVector[IExtrapolateToCaloTool::Back]) {
- parameterVector[IExtrapolateToCaloTool::Back] = extrapolate(parameters, sample, m_offsetBack, particle, dir);
- }
- if (parameterVector[IExtrapolateToCaloTool::Back]) {
- lastParamInCalo = parameterVector[IExtrapolateToCaloTool::Back];
- }
-
- // Tile or HEC0 :
- distbar = m_calodepth->deta(CaloCell_ID::TileBar0,trketa);
- distec = m_calodepth->deta(CaloCell_ID::HEC0,trketa);
- ATH_MSG_DEBUG(" TrackTo ...Tile : for eta= " << trketa << " dist to Barrel =" << distbar
- << " to endcap =" << distec);
- if (distbar < 0 ) sample = CaloCell_ID::TileBar0;
- else if (distec < 0 ) sample = CaloCell_ID::HEC0;
- else if ( distbar > distec && distec < 10. ) sample = CaloCell_ID::HEC0;
- else sample = CaloCell_ID::TileBar0;
- ATH_MSG_DEBUG( " => will shoot in sample " << sample );
- parameterVector[IExtrapolateToCaloTool::Tile] = extrapolate(*lastParamInCalo, sample, m_offsetTile, particleHypoStepWiseInCalo, dir);
- if (!parameterVector[IExtrapolateToCaloTool::Tile]) {
- parameterVector[IExtrapolateToCaloTool::Tile] = extrapolate(parameters, sample, m_offsetTile, particle, dir);
- }
- return;
-}
-
-Amg::Vector3D ExtrapolateToCaloTool::getMomentumAtVertex(const xAOD::Vertex& vertex,
- bool reuse /* = true */) const
-{
- Amg::Vector3D momentum(0., 0., 0.);
-
- // xAOD::Vertex does not have method isAvailable for the moment
- // create a lambda function
- auto isAvailable = [](const xAOD::Vertex& vertex, std::string name) {
- SG::AuxElement::Accessor<float> acc(name, "");
- return acc.isAvailable(vertex); } ;
-
- if (vertex.nTrackParticles() == 0)
- {
- ATH_MSG_WARNING("getMomentumAtVertex : vertex has no track particles!");
- return momentum;
- }
-
- if (reuse &&
- isAvailable(vertex, "px") &&
- isAvailable(vertex, "py") &&
- isAvailable(vertex, "pz") )
- {
- // Already decorated with parameters at vertex
- ATH_MSG_DEBUG("getMomentumAtVertex : getting from auxdata");
- return Amg::Vector3D( vertex.auxdata<float>("px"),
- vertex.auxdata<float>("py"),
- vertex.auxdata<float>("pz") );
- }
- else if (vertex.vxTrackAtVertexAvailable() && vertex.vxTrackAtVertex().size())
- {
- // Use the parameters at the vertex
- // (the tracks should be parallel but we will do the sum anyway)
- ATH_MSG_DEBUG("getMomentumAtVertex : getting from vxTrackAtVertex");
- for (const auto& trkAtVertex : vertex.vxTrackAtVertex())
- {
- const Trk::TrackParameters* paramAtVertex = trkAtVertex.perigeeAtVertex();
- if (!paramAtVertex)
- ATH_MSG_WARNING("VxTrackAtVertex does not have perigee at vertex");
- else
- momentum += paramAtVertex->momentum();
- }
- }
- else if (vertex.nTrackParticles() == 1)
- {
- // Use the first measurement
- ATH_MSG_DEBUG("getMomentumAtVertex : 1 track only, getting from first measurement");
- const xAOD::TrackParticle *tp = vertex.trackParticle(0);
- unsigned int index(0);
- if (!tp || !tp->indexOfParameterAtPosition(index, xAOD::FirstMeasurement))
- {
- ATH_MSG_WARNING("No TrackParticle or no have first measurement");
- }
- else
- momentum += tp->curvilinearParameters(index).momentum();
- // OR last 3 values of trackParameters(index)
- }
- else
- {
- // Extrapolate track particles to vertex
- // (the tracks should be parallel but we will do the sum anyway)
- ATH_MSG_DEBUG("getMomentumAtVertex : extrapolating to perigee surface");
- const Trk::PerigeeSurface *surface = new Trk::PerigeeSurface(vertex.position());
- for (unsigned int i = 0; i < vertex.nTrackParticles(); ++i)
- {
- const xAOD::TrackParticle* tp = vertex.trackParticle( i );
- if (!tp)
- {
- ATH_MSG_WARNING("NULL pointer to TrackParticle in vertex");
- continue;
- }
- const Trk::TrackParameters* params = m_extrapolator->extrapolate(*tp, *surface, Trk::alongMomentum);
- if (!params)
- ATH_MSG_DEBUG("Extrapolation to vertex (perigee) failed");
- else
- {
- momentum += params->momentum();
- delete params;
- }
- }
- delete surface;
- }
-
- return momentum;
-}
-
-Trk::SurfaceIntersection ExtrapolateToCaloTool::getIntersectionInCalo(const Amg::Vector3D& position, const Amg::Vector3D& momentum, const CaloCell_ID::CaloSample sample) const
-{
- Trk::SurfaceIntersection result{Amg::Vector3D(0., 0., 0.), 0., false};
-
- // get the destination Surface
- const Trk::Surface* surface = m_calosurf->CreateUserSurface (sample, 0. /* offset */, momentum.eta());
- if (!surface)
- return result;
-
- // intersect with calorimeter surface
- result = surface->straightLineIntersection(position, momentum, true);
- delete surface;
- return result;
-}
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/ImpactInCalo.cxx b/Reconstruction/RecoTools/TrackToCalo/src/ImpactInCalo.cxx
deleted file mode 100755
index 8d7a1a4133bd322afb058c9a69035fbf2c9ee49c..0000000000000000000000000000000000000000
--- a/Reconstruction/RecoTools/TrackToCalo/src/ImpactInCalo.cxx
+++ /dev/null
@@ -1,48 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-
-#include "TrackToCalo/ImpactInCalo.h"
-
-#ifdef HAVE_NEW_IOSTREAMS
-#include <iostream>
-#else
-#include <iostream.h>
-#define std
-#endif
-
-
-
-void ImpactInCalo::print() const {
-
- std::cout << "Impact in PS: x,y,z=" << m_x_0 << " " << m_y_0
- << " " << m_z_0
- << " local eta,phi=" << m_etaCaloLocal_0 << " "
- << m_phiCaloLocal_0
- << std::endl;
-
- std::cout << "Impact in Strip: x,y,z=" << m_x_1 << " " << m_y_1
- << " " << m_z_1
- << " local eta,phi=" << m_etaCaloLocal_1 << " "
- << m_phiCaloLocal_1
- << std::endl;
-
- std::cout << "Impact in Middle: x,y,z=" << m_x_2 << " " << m_y_2
- << " " << m_z_2
- << " local eta,phi=" << m_etaCaloLocal_2
- << " " << m_phiCaloLocal_2
- << std::endl;
-
- std::cout << "Impact in Back: x,y,z=" << m_x_3 << " " << m_y_3
- << " " << m_z_3
- << " local eta,phi=" << m_etaCaloLocal_3 << " "
- << m_phiCaloLocal_3
- << std::endl;
-
- std::cout << "Impact in Tiles: x,y,z=" << m_x_tile << " "
- << m_y_tile << " " << m_z_tile
- << " local eta,phi=" << m_etaCaloLocal_tile
- << " " << m_phiCaloLocal_tile
- << std::endl;
-}
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/ImpactInCaloAlg.cxx b/Reconstruction/RecoTools/TrackToCalo/src/ImpactInCaloAlg.cxx
deleted file mode 100755
index 6090ce695ccc0c55e68e171165fd1f38d76605cf..0000000000000000000000000000000000000000
--- a/Reconstruction/RecoTools/TrackToCalo/src/ImpactInCaloAlg.cxx
+++ /dev/null
@@ -1,463 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-///////////////////////////////////////////////////////////////////
-// ImpactInCaloAlg.cxx
-// Implementation file for class ImpactInCaloAlg
-///////////////////////////////////////////////////////////////////
-
-#include "TrackToCalo/ImpactInCaloAlg.h"
-
-// Gaudi includes
-#include "GaudiKernel/MsgStream.h"
-#include "StoreGate/StoreGate.h"
-//#include "GaudiKernel/AlgFactory.h"
-#include "StoreGate/StoreGateSvc.h"
-
-// Tracking includes
-#include "TrkParameters/TrackParameters.h"
-#include <vector>
-#include "Particle/TrackParticleContainer.h"
-#include "Particle/TrackParticle.h"
-
-#include "TrackToCalo/ImpactInCalo.h"
-
-#include "CaloEvent/CaloCluster.h"
-#include "CaloEvent/CaloClusterContainer.h"
-#include "CaloDetDescr/CaloDetDescrManager.h"
-#include "CaloGeoHelpers/CaloPhiRange.h"
-#include "CaloUtils/CaloCellList.h"
-#include "CaloEvent/CaloCell.h"
-#include "CaloEvent/CaloCellContainer.h"
-#include "CaloDetDescr/ICaloCoordinateTool.h"
-#include "CaloDetDescr/CaloDepthTool.h"
-
-#include "RecoToolInterfaces/IExtrapolateToCaloTool.h"
-#include "TrackToCalo/ImpactInCaloCollection.h"
-#include "TrkToolInterfaces/ITrackSelectorTool.h"
-
-// Constructor with parameters:
-ImpactInCaloAlg::ImpactInCaloAlg(const std::string &name,
- ISvcLocator *pSvcLocator) :
- AthAlgorithm(name,pSvcLocator),
- m_TrackName("Tracks"),
- m_tracks(0),
- m_particle(0),
- m_calo_id(0),
- m_calo_dd(0),
- m_toCalo("ExtrapolateToCaloTool"),
- m_trackSelectorTool(""),
- m_extrapolateInBothDirections(false),
- m_compareImpWithCluster(false),
- m_printCellsCrossed(false) {
- // Get parameter values from jobOptions file
- declareProperty("TrackName", m_TrackName);
- declareProperty("TrackParticleName", m_TrackParticleName);
- declareProperty("ClusterContainerName", m_cluster_container);
- declareProperty("CaloCellContainerName", m_cell_container);
- declareProperty("ImpactInCaloContainerName", m_ImpactInCalosOutputName);
- declareProperty("TrackInputType", m_trkinput);
- declareProperty("ExtrapolateToCaloTool", m_toCalo, "the extrapolate to calo tool");
- declareProperty("TrackSelectorTool", m_trackSelectorTool, "Tool handle to TrackSelectorTool (if no tool given: use all tracks)");
- declareProperty("CompareImpactWithCluster", m_compareImpWithCluster, "compare impact with clusters?");
- declareProperty("PrintCellsCrossed", m_printCellsCrossed, "print all cells crossed by tracks?");
- declareProperty("ExtrapolateInBothDirections", m_extrapolateInBothDirections, "extrapolate each track in both directions (along and opposite momentum)?");
-}
-
-// Initialize method:
-StatusCode ImpactInCaloAlg::initialize() {
-
- // Get the messaging service, print where you are
- ATH_MSG_DEBUG( "ImpactInCaloAlg::initialize()" );
-
- //StatusCode sc;
-
- m_calo_dd = CaloDetDescrManager::instance();
- m_calo_id = m_calo_dd->getCaloCell_ID();
-
- if (m_toCalo.retrieve().isFailure()) {
- ATH_MSG_FATAL( "Could not retrieve " << m_toCalo );
- return StatusCode::FAILURE;
- }
-
- ATH_MSG_INFO( "ImpactInCaloAlg initialisation OK" );
- return StatusCode::SUCCESS;
-}
-
-// Execute method:
-StatusCode ImpactInCaloAlg::execute() {
- // Get the messaging service, print where you are
- //StatusCode sc;
- //ATH_MSG_DEBUG( "ImpactInCaloAlg::execute()" );
-
- // Create the impacts :
- StatusCode sc = CreateTrkImpactInCalo();
- if (sc.isFailure())
- return StatusCode::SUCCESS;
-
- // Example 1 : you want to read the impacts and compare to clusters
- if (m_compareImpWithCluster) {
- sc = PrintImpact();
- if (sc.isFailure())
- return StatusCode::SUCCESS;
- sc = CompareImpactWithCluster();
- if (sc.isFailure())
- return StatusCode::SUCCESS;
- }
-
- // Example 2 : you want to know the list of cells crossed by the track
- if (m_printCellsCrossed) {
- sc = PrintCellsCrossed();
- if (sc.isFailure())
- return StatusCode::SUCCESS;
- }
-
- return StatusCode::SUCCESS;
-}
-
-// Finalize method:
-StatusCode ImpactInCaloAlg::finalize() {
- ATH_MSG_DEBUG( "ImpactInCaloAlg::finalize()" );
-
- return StatusCode::SUCCESS;
-}
-
-
-// Loop on Tracks/TrackParticles and create ImpactInCaloCollection
-StatusCode ImpactInCaloAlg::CreateTrkImpactInCalo() {
- //StatusCode sc;
- ATH_MSG_DEBUG( "ImpactInCaloAlg::CreateTrkImpactInCalo()" );
-
- //create and record new ImpactInCalo container
- ImpactInCaloCollection* outputContainer = new ImpactInCaloCollection();
- StatusCode sc=evtStore()->record(outputContainer,m_ImpactInCalosOutputName,false); //"false" indicates it can't be modified by another algorithm
- if(sc != StatusCode::SUCCESS) {
- ATH_MSG_WARNING(" Could not record ImpactInCaloCollection"
- << m_ImpactInCalosOutputName);
- return StatusCode::FAILURE;
- }
-
- m_particle = 0;
- m_tracks = 0;
-
- if( m_trkinput == "TrackParticleCandidates") {
- // use track particle collection
- if (m_TrackParticleName == "") {
- ATH_MSG_WARNING("TrackParticleName not set" );
- return StatusCode::FAILURE;
- }
-
- sc = evtStore()->retrieve(m_particle, m_TrackParticleName);
-
- if (sc.isFailure()) {
- ATH_MSG_WARNING("TrackParticle not found: " << m_TrackParticleName );
- return StatusCode::FAILURE;
- }
- ATH_MSG_DEBUG("TrackParticle found in StoreGate" );
- // loop over track particles and get impacts in calo:
- for (Rec::TrackParticleContainer::const_iterator itr = (*m_particle).begin(); itr != (*m_particle).end(); itr++) {
- if (!(*itr)) continue;
- if (!m_trackSelectorTool.empty()) {
- if (!m_trackSelectorTool->decision(**itr)) {
- // track out of selection cuts
- continue;
- }
- }
- // get impacts from ExtrapolateToCaloTool
- ImpactInCalo* imp = m_toCalo->getImpactInCalo(**itr, Trk::undefined, Trk::alongMomentum);
- if(imp)
- outputContainer->push_back(imp);
- else
- ATH_MSG_DEBUG(" ImpactInCalo pointer not valid for this track");
- if (m_extrapolateInBothDirections) {
- ImpactInCalo* imp = m_toCalo->getImpactInCalo(**itr, Trk::undefined, Trk::oppositeMomentum);
- if (imp) outputContainer->push_back(imp);
- }
- } // end loop over track particles
- } else {
- // use track collection
- if (m_TrackName == "") {
- ATH_MSG_WARNING("TrackName not set" );
- return StatusCode::FAILURE;
- }
-
- sc = evtStore()->retrieve(m_tracks, m_TrackName);
-
- if (sc.isFailure()) {
- ATH_MSG_WARNING("Tracks not found: " << m_TrackName );
- return StatusCode::FAILURE;
- }
- ATH_MSG_DEBUG("Tracks found in StoreGate" );
- // loop over tracks
- for (TrackCollection::const_iterator itr = (*m_tracks).begin(); itr < (*m_tracks).end(); itr++) {
- if (!(*itr))
- continue;
- if (!m_trackSelectorTool.empty()) {
- if (!m_trackSelectorTool->decision(**itr)) {
- // track out of selection cuts
- continue;
- }
- }
- // get impacts from ExtrapolateToCaloTool
- ImpactInCalo* imp = m_toCalo->getImpactInCalo(*(*itr), Trk::undefined, Trk::alongMomentum);
- if(imp)
- outputContainer->push_back(imp);
- else
- ATH_MSG_DEBUG(" ImpactInCalo pointer not valid for this track");
- if (m_extrapolateInBothDirections) {
- ImpactInCalo* imp = m_toCalo->getImpactInCalo(*(*itr), Trk::undefined, Trk::oppositeMomentum);
- if (imp) outputContainer->push_back(imp);
- }
- } // end loop over tracks
- }
- return StatusCode::SUCCESS;
-}
-
-
-// Retreive ImpactInCaloCollection and compare to CaloClusters
-StatusCode ImpactInCaloAlg::CompareImpactWithCluster() {
- StatusCode sc;
-
- ATH_MSG_DEBUG( "ImpactInCaloAlg::CompareImpactWithCluster()" );
-
- // loop on clusters
-
- const CaloClusterContainer* cluster_container;
- sc=evtStore()->retrieve(cluster_container,m_cluster_container);
- if (sc.isFailure()) {
- ATH_MSG_WARNING( "Cluster Container could not be retrieved: " << m_cluster_container );
- return sc;
- }
-
- const ImpactInCaloCollection* impact_collection;
- sc=evtStore()->retrieve(impact_collection,m_ImpactInCalosOutputName);
- if (sc.isFailure()) {
- ATH_MSG_WARNING( "Impact collection could not be retrieved: " << m_ImpactInCalosOutputName );
- return sc;
- }
-
- typedef CaloClusterContainer::const_iterator cluster_iterator;
- cluster_iterator f_clu =cluster_container->begin();
- cluster_iterator l_clu = cluster_container->end();
-
- typedef ImpactInCaloCollection::const_iterator impact_iterator;
- impact_iterator f_imp = impact_collection->begin();
- impact_iterator l_imp = impact_collection->end();
-
- for ( ; f_clu!=l_clu; f_clu++) {
- const CaloCluster* cluster = (*f_clu);
- double hecluster = cluster->energy()/1000.;
- double heta = cluster->eta();
- double hphi = cluster->phi();
-
- ATH_MSG_INFO("Found a cluster : E= " << hecluster
- << "(GeV), etaCaloLocal=" << heta
- << ", phiCaloLocal=" << hphi);
-
- for ( ; f_imp!=l_imp; f_imp++) {
- const ImpactInCalo* impact = (*f_imp);
-
- ATH_MSG_INFO("==> Comparison between cluster and impact in Middle : deta="
- << heta - impact->etaCaloLocal_2()
- << " , dphi=" << hphi - impact->phiCaloLocal_2());
- }
- }
- return StatusCode::SUCCESS;
-}
-
-// Retreive ImpactInCaloCollection and print content
-StatusCode ImpactInCaloAlg::PrintImpact() {
- StatusCode sc;
-
-
- ATH_MSG_INFO( " Method PrintImpacts : " );
-
- ATH_MSG_INFO( " " );
- ATH_MSG_INFO( " Start with Impacts : " );
-
- const ImpactInCaloCollection* impact_collection;
- sc=evtStore()->retrieve(impact_collection,m_ImpactInCalosOutputName);
-
- if(sc.isFailure()) {
- ATH_MSG_WARNING( "Could not retrieve impact colection: " << m_ImpactInCalosOutputName );
- return sc;
- }
-
- typedef ImpactInCaloCollection::const_iterator impact_iterator;
- impact_iterator f_imp = impact_collection->begin();
- impact_iterator l_imp = impact_collection->end();
-
- for ( ; f_imp!=l_imp; f_imp++) {
- const ImpactInCalo* impact = (*f_imp);
-
- const double impPhi = impact->phiCaloLocal_1();
- const double impeta = impact->etaCaloLocal_1();
-
- ATH_MSG_INFO("Found an impact in strips : parameters are eta = " << impeta
- << " phi = " << impPhi);
-
- //impact->print();
- }
-
- if (m_tracks) {
-
- ATH_MSG_INFO( " " );
- ATH_MSG_INFO(" Now loop on Trk::Track collection " );
- ATH_MSG_INFO( " " );
-
- for (TrackCollection::const_iterator itr =
- (*m_tracks).begin(); itr < (*m_tracks).end(); itr++) {
- const Trk::Perigee *aPer=(*itr)->perigeeParameters();
- if (!aPer) {
- ATH_MSG_WARNING( "Could not get Trk::Perigee from Track" );
- } else {
-
- const double trkphi = aPer->parameters()[Trk::phi];
- const double trketa = aPer->eta();
-
- ATH_MSG_INFO("Found a Trk::Track : parameters are eta = " << trketa
- << " phi = " << trkphi);
- }
- }
- }
-
- if (m_particle) {
-
- ATH_MSG_INFO( " " );
- ATH_MSG_INFO(" Now loop on Trk::TrackParticle collection " );
- ATH_MSG_INFO( " " );
-
- for (Rec::TrackParticleContainer::const_iterator itr = m_particle->begin();
- itr != m_particle->end(); ++itr ) {
- if (!(*itr)->originalTrack()) {
- ATH_MSG_WARNING( "Could not get Track from TrackParticle" );
- continue;
- }
- const Trk::Perigee *aPer=((*itr)->originalTrack())->perigeeParameters();
- if (!aPer) {
- ATH_MSG_WARNING( "Could not get Perigee from TrackParticle" );
- } else {
-
- const double partphi = aPer->parameters()[Trk::phi];
- const double parteta = aPer->eta();
-
- ATH_MSG_INFO("Found a trackparticle : parameters are eta = " << parteta
- << " phi = " << partphi);
- }
- }
- }
- return StatusCode::SUCCESS;
-}
-
-// List of cells crossed by Trk::Tracks in CaloSample
-// ( neta and nphi are the total window width, e.g .5x5 )
-CaloCellList* ImpactInCaloAlg::CellsCrossedByTrack(const Trk::Track& trk,
- const CaloCell_ID::CaloSample sam,
- int neta, int nphi) {
- CaloCellList* my_list = 0;
-
- // Retreive CaloCell's from StoreGate :
-
- const CaloCellContainer* cell_container;
- StatusCode sc=evtStore()->retrieve(cell_container,m_cell_container);
- if ( sc.isFailure() )
- return 0;
-
- // Where is the track shooting ?
- double offset = 0.;
-
- const Trk::TrackParameters* parametersInCalo = m_toCalo->extrapolate (trk, sam, offset);
-
- if (!parametersInCalo)
- return 0;
-
- double eta = parametersInCalo->position().eta();
- double phi = parametersInCalo->position().phi();
-
- // CaloCellList needs both enums: subCalo and CaloSample
- CaloCell_ID::SUBCALO subcalo;
- bool barrel;
- int sampling_or_module;
- m_calo_dd->decode_sample (subcalo, barrel, sampling_or_module, sam);
-
- // Get the corresponding granularities : needs to know where you are
- // the easiest is to look for the CaloDetDescrElement
- const CaloDetDescrElement* dde =
- m_calo_dd->get_element(subcalo,sampling_or_module,barrel,eta,phi);
-
- double deta = int(neta/2)*dde->deta();
- double dphi = int(nphi/2)*dde->dphi();
-
- // Construct the list :
- my_list = new CaloCellList(cell_container,subcalo);
- my_list->select(eta,phi,deta,dphi, (int) sam);
-
- // cleanup
- delete parametersInCalo;
-
- return my_list ;
-}
-
-StatusCode ImpactInCaloAlg::PrintCellsCrossed() {
-
- // Get the messaging service, print where you are
- StatusCode sc;
- ATH_MSG_INFO( "ImpactInCaloAlg::PrintCellsCrossed()" );
-
- // Here we are :
- CaloCell_ID::CaloSample sam = CaloCell_ID::EMB2;
- int neta = 5;
- int nphi = 5;
-
- // get tracks from TDS
- if (m_TrackName == "") {
- ATH_MSG_WARNING("TrackName not set" );
- return StatusCode::FAILURE;
- }
-
- sc = evtStore()->retrieve(m_tracks, m_TrackName);
- if (sc.isFailure()) {
- ATH_MSG_WARNING("Tracks not found: " << m_TrackName );
- return sc;
- }
- ATH_MSG_DEBUG("Tracks found in StoreGate" );
-
- for (TrackCollection::const_iterator itr = (*m_tracks).begin(); itr != (*m_tracks).end(); itr++) {
- if (!(*itr)) {
- ATH_MSG_INFO( "Trackcollection contains empty entry" );
- continue;
- }
- const Trk::Perigee *aPer= (*itr)->perigeeParameters();
- if (!aPer) {
- ATH_MSG_INFO( "Could not get Trk::Perigee from Track" );
- } else {
- double d0 = aPer->parameters()[Trk::d0];
- double z0 = aPer->parameters()[Trk::z0];
- double phi0 = aPer->parameters()[Trk::phi0];
- double theta = aPer->parameters()[Trk::theta];
- double qOverP = aPer->parameters()[Trk::qOverP];
- ATH_MSG_INFO( " " );
- ATH_MSG_INFO("Found a track: parameters are (" << d0 << ", "
- << z0 << ", " << phi0 << ", " << theta << ", " << qOverP << ")");
- }
-
- // Warning : if anything fails, CellsCrossedByTrack will
- // simply return a null pointer
- // if it works, it does a new CaloCellList
- // ==> the client has to do the delete !!!!
-
- CaloCellList* my_list = CellsCrossedByTrack(*(*itr), sam, neta, nphi);
-
- if (my_list) {
-
- for ( CaloCellList::list_iterator itr = my_list->begin(); itr < my_list->end(); itr++) {
- ATH_MSG_INFO("found cell ! eta=" << (*itr)->eta()
- << " phi=" << (*itr)->phi() << " energy=" << (*itr)->energy());
- }
- delete my_list;
- } else
- ATH_MSG_INFO( "could not get cell list " );
- }
- return StatusCode::SUCCESS;
-}
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/ImpactInCaloCollection.cxx b/Reconstruction/RecoTools/TrackToCalo/src/ImpactInCaloCollection.cxx
deleted file mode 100755
index 70a1c72fbba557caaeb55707352a44605a5b320e..0000000000000000000000000000000000000000
--- a/Reconstruction/RecoTools/TrackToCalo/src/ImpactInCaloCollection.cxx
+++ /dev/null
@@ -1,13 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-
-#include "TrackToCalo/ImpactInCaloCollection.h"
-#include "TrackToCalo/ImpactInCalo.h"
-
-void ImpactInCaloCollection::print() const{
- ImpactInCaloCollection::const_iterator iter;
- for( iter=begin(); iter!=end(); ++iter ) (*iter)->print();
-}
-
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/MuonCaloEnergyTool.cxx b/Reconstruction/RecoTools/TrackToCalo/src/MuonCaloEnergyTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..1175e27a17ed57b815dc97ec4be27e0f1eaf2fde
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/src/MuonCaloEnergyTool.cxx
@@ -0,0 +1,428 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "MuonCaloEnergyTool.h"
+// forward declares
+#include "ParticleCaloExtension/ParticleCellAssociationCollection.h"
+
+#include "TrackToCalo/CaloCellHelpers.h"
+#include "ParticleCaloExtension/ParticleCellIntersection.h"
+#include "ParticleCaloExtension/ParticleCaloAssociation.h"
+
+#include "CaloUtils/CaloCellList.h"
+#include "CaloEvent/CaloCellContainer.h"
+#include "TrkCaloExtension/CaloExtension.h"
+#include "TrkCaloExtension/CaloExtensionHelpers.h"
+
+#include "xAODTracking/TrackParticle.h"
+#include "xAODTracking/TrackingPrimitives.h"
+#include <math.h>
+
+namespace Rec {
+
+ MuonCaloEnergyTool::MuonCaloEnergyTool(const std::string& t, const std::string& n, const IInterface* p )
+ : AthAlgTool(t,n,p),
+ m_caloExtensionTool("Trk::ParticleCaloExtensionTool/ParticleCaloExtensionTool"),
+ m_caloCellAssociationTool("Rec::ParticleCaloCellAssociationTool/ParticleCaloCellAssociationTool"),
+ m_particleCreator("Trk::TrackParticleCreatorTool/MuonCaloParticleCreator"),
+ m_caloNoiseTool("CaloNoiseToolDefault"),
+ m_sigmasAboveNoise(4.),
+ m_emEtCut(2.5*CLHEP::GeV),
+ m_emF1Cut(0.15),
+ m_emipEM(0.78),
+ m_emipTile(0.86),
+ m_emipHEC(0.94)
+ {
+ declareInterface<IMuonCaloEnergyTool>(this);
+ declareProperty("ParticleCaloExtensionTool", m_caloExtensionTool );
+ declareProperty("ParticleCaloCellAssociationTool",m_caloCellAssociationTool );
+ declareProperty("TrackParticleCreator", m_particleCreator );
+ declareProperty("CaloNoiseTool", m_caloNoiseTool);
+ declareProperty("TrackParticleLocation", m_indetTrackParticleLocation = "InDetTrackParticles");
+ declareProperty("MuonSpectrometerTrackPArticleLocation",m_muonTrackParticleLocation = "MuonSpectrometerTrackParticles");
+
+ //coneSize for including calo cells around track
+ declareProperty("SigmasAboveNoise", m_sigmasAboveNoise = 4.);
+ declareProperty("EmEtCut", m_emEtCut = 2.5*1000.); // in MeV
+ declareProperty("EmF1Cut", m_emF1Cut = 0.15);
+ }
+
+ MuonCaloEnergyTool::~MuonCaloEnergyTool() {}
+
+ StatusCode MuonCaloEnergyTool::initialize() {
+
+ // RETRIEVE TOOLS
+ ATH_CHECK(m_caloExtensionTool.retrieve());
+ ATH_CHECK(m_caloCellAssociationTool.retrieve());
+ ATH_CHECK(m_particleCreator.retrieve());
+ ATH_CHECK(m_caloNoiseTool.retrieve());
+
+ m_indetTrackParticles = 0;
+ if(evtStore()->contains<xAOD::TrackParticleContainer>(m_indetTrackParticleLocation)) {
+ if(evtStore()->retrieve(m_indetTrackParticles,m_indetTrackParticleLocation).isFailure()) {
+ ATH_MSG_FATAL( "Unable to retrieve " << m_indetTrackParticleLocation );
+ return StatusCode::FAILURE;
+ }
+ }
+
+ m_muonTrackParticles = 0;
+ if(evtStore()->contains<xAOD::TrackParticleContainer>(m_muonTrackParticleLocation)) {
+ if(evtStore()->retrieve(m_muonTrackParticles,m_muonTrackParticleLocation).isFailure()) {
+ ATH_MSG_FATAL( "Unable to retrieve " << m_muonTrackParticleLocation );
+ return StatusCode::FAILURE;
+ }
+ }
+
+
+ return StatusCode::SUCCESS;
+ }
+
+ StatusCode MuonCaloEnergyTool::finalize() {
+ return StatusCode::SUCCESS;
+ }
+
+ void MuonCaloEnergyTool::calculateMuonEnergies( const Trk::Track* trk,
+ double deltaE, double meanIoni, double sigmaIoni,
+ double& E, double& sigma, double& E_FSR , double& E_expected,
+ double &E_em_meas, double &E_em_exp, double &E_tile_meas, double &E_tile_exp,
+ double &E_HEC_meas, double &E_HEC_exp, double &E_dead_exp ) const {
+
+//
+// Input parameters trk: (muon) track pointer
+// deltaE: Mean Energy loss in Calorimeter
+// meanIoni: Mean ionization Energy loss in Calorimeter
+// sigmaIoni: sigma of the ionization Energy loss in Calorimeter
+//
+// Ouput parameters E: Calorimeter energy measured (corrected for dead material etc.)
+// sigma: Error on measured energy
+// E_FSR: Energy of Final State Radiation: e.m. = photon energy (if above ET and F1 cut)
+// if EFSR>0 the corrected hadronic energy is used for the measured Energy E
+// E_expected: expected meanIonization Energyloss from the TrackingGeometry
+
+ E = 0.;
+ sigma = 0.;
+ E_FSR = 0.;
+ E_expected = 0.;
+
+//
+// For the expected Eloss in the dead or not instrumented material we will use the meanIonization loss
+// this meanIoni is stored on the extended Eloss object
+// the sigmaIoni corresponds to the Landau width
+// a factor 3.59524 gives the sigma of the upper tail
+//
+// The Eloss derived from the momentum differences on the calo extension correponds to
+// deltaE on the Eloss object = mean Eloss including ionization and radiation
+//
+// To go to the meanIonization loss one needs a scale factor scale_Ionization
+//
+ double scale_Ionization = 0.9;
+ if(fabs(deltaE)>0&&fabs(meanIoni)>0) scale_Ionization = meanIoni / deltaE;
+ double error_ratio = 0.3;
+ if(fabs(meanIoni)>0&&sigmaIoni>0) error_ratio = 3.59524*fabs(sigmaIoni/meanIoni);
+
+ ATH_MSG_DEBUG( " Inputs deltaE " << deltaE << " meanIoni " << meanIoni << " sigmaIoni " << sigmaIoni );
+ if(deltaE==0||meanIoni==0||sigmaIoni<0) ATH_MSG_WARNING( " Strange Inputs deltaE " << deltaE << " meanIoni " << meanIoni << " sigmaIoni " << sigmaIoni );
+
+ // associate muon to calorimeter
+
+ const xAOD::TrackParticle* tp = 0;
+
+// check ID trackparticles
+
+ if(m_indetTrackParticles) {
+ for(auto it : *m_indetTrackParticles){
+ if((*it).track()==trk) {
+ tp = &(*it);
+ break;
+ }
+ }
+ if(tp) ATH_MSG_DEBUG( " ID xAOD::TrackParticle found " << tp );
+ }
+
+// look for Muon trackparticles
+
+ if(m_muonTrackParticles&&!tp) {
+ for(auto it : *m_muonTrackParticles){
+ if((*it).track()==trk) {
+ tp = &(*it);
+ break;
+ }
+ }
+ if(tp) ATH_MSG_DEBUG( " Muon xAOD::TrackParticle found " << tp );
+ }
+
+ if(!tp) {
+ tp = m_particleCreator->createParticle(*trk, NULL, NULL, xAOD::muon);
+ if(tp) ATH_MSG_DEBUG( " xAOD::TrackParticle created from scratch " << tp );
+ }
+
+ if(!tp) return;
+
+ const Rec::ParticleCellAssociation* association = 0;
+ m_caloCellAssociationTool->particleCellAssociation(*tp,association,0.1,NULL,true);
+
+ if(!association) return;
+ ATH_MSG_DEBUG( " particleCellAssociation done " << association );
+
+ std::vector< std::pair<const CaloCell*,Rec::ParticleCellIntersection*> > cellIntersections = association->cellIntersections();
+
+ const Trk::CaloExtension& caloExtension = association->caloExtension();
+
+ if(!(&caloExtension)) {
+ ATH_MSG_WARNING( " No caloExtension found ");
+ return;
+ }
+
+ if(!caloExtension.caloEntryLayerIntersection()) {
+ ATH_MSG_WARNING( " No caloEntryLayerIntersection found ");
+ return;
+ }
+
+ ATH_MSG_DEBUG( " nr of cell intersections " << cellIntersections.size() );
+ if( cellIntersections.size() < 3) ATH_MSG_WARNING( " Strange nr of calorimeter cell intersections " << cellIntersections.size() );
+
+ double theta = caloExtension.caloEntryLayerIntersection()->position().theta();
+
+ double Eloss = 0.;
+ if(!caloExtension.muonEntryLayerIntersection()) {
+ ATH_MSG_WARNING( " No muonEntryLayerIntersection found and therefore the expected Eloss is not calculated properly ");
+ } else {
+ Eloss = caloExtension.caloEntryLayerIntersection()->momentum().mag() - caloExtension.muonEntryLayerIntersection()->momentum().mag();
+ ATH_MSG_DEBUG( " Energy loss from CaloExtension " << Eloss << " R muon Entry " << caloExtension.muonEntryLayerIntersection()->position().perp() << " Z muon Entry " << caloExtension.muonEntryLayerIntersection()->position().z() << " R calo entry " << caloExtension.caloEntryLayerIntersection()->position().perp() << " Z calo entry " << caloExtension.caloEntryLayerIntersection()->position().z() );
+ }
+
+ //auto cellcoll = association->data();
+
+ // measured and expected energies
+
+ double E_em1 = 0.;
+ double E_em = 0.;
+ double E_em_expected = 0.;
+ double E_em_exptot = 0.;
+ double E_em_expall = 0.;
+ double E_em_threshold = 0.;
+ int nlay_em = 0;
+ double sigma_Noise_em = 0.;
+
+ double E_HEC = 0.;
+ double E_HEC_expected = 0.;
+ double E_HEC_threshold = 0.;
+ double E_HEC_exptot = 0.;
+ double E_HEC_expall = 0.;
+ int nlay_HEC = 0;
+ double sigma_Noise_HEC = 0.;
+
+ double E_tile = 0.;
+ double E_tile_expected = 0.;
+ double E_tile_threshold = 0.;
+ double E_tile_exptot = 0.;
+ double E_tile_expall = 0.;
+ int nlay_tile = 0;
+ double sigma_Noise_tile = 0.;
+
+// const char* sampname[24] = {
+// "PreSamplerB", "EMB1", "EMB2", "EMB3",
+// "PreSamplerE", "EME1", "EME2", "EME3",
+// "HEC0", "HEC1","HEC2", "HEC3",
+// "TileBar0", "TileBar1", "TileBar2",
+// "TileExt0", "TileExt1", "TileExt2",
+// "TileGap1", "TileGap2", "TileGap3",
+// "FCAL0", "FCAL1", "FCAL2"};
+// for(const auto& curr_cell : cellcoll ){
+
+ E_expected = 0.;
+ double phiPos = caloExtension.caloEntryLayerIntersection()->position().phi();
+ double thetaPos = caloExtension.caloEntryLayerIntersection()->position().theta();
+
+ for(auto it : cellIntersections){
+
+ const CaloCell* curr_cell = it.first;
+ int cellSampling = curr_cell->caloDDE()->getSampling();
+ bool badCell = curr_cell->badcell();
+ double cellEn = curr_cell->energy();
+
+ double f_exp = (it.second)->pathLength();
+ double E_exp = (it.second)->expectedEnergyLoss();
+
+// if(f_exp<0.1) f_exp = 0.1;
+
+// cellnoisedb = m_caloNoiseTool->getNoise(cell,ICalorimeterNoiseTool::ELECTRONICNOISE);
+ double sigma_Noise = m_caloNoiseTool->getEffectiveSigma(curr_cell,ICalorimeterNoiseTool::MAXSYMMETRYHANDLING,ICalorimeterNoiseTool::ELECTRONICNOISE);
+// double sigma_NoiseA = m_caloNoiseTool->totalNoiseRMS(curr_cell);
+ double thetaCell = atan2(sqrt(curr_cell->x()*curr_cell->x()+curr_cell->y()*curr_cell->y()),curr_cell->z());
+ double phiCell = atan2(curr_cell->y(),curr_cell->x());
+
+ double dtheta = thetaCell - thetaPos;
+ double dphi = acos(cos(phiPos)*cos(phiCell)+sin(phiPos)*sin(phiCell));
+ double celldr = curr_cell->caloDDE()->dr();
+ double celldz = curr_cell->caloDDE()->dz();
+ double celldtheta = celldr/sqrt(curr_cell->x()*curr_cell->x()+curr_cell->y()*curr_cell->y());
+ double celldphi = curr_cell->caloDDE()->dphi();
+
+
+ ATH_MSG_DEBUG( " cell sampling " << cellSampling << " energy " << cellEn << " position R " << sqrt(curr_cell->x()*curr_cell->x()+curr_cell->y()*curr_cell->y()) << " z " << curr_cell->z() << " f_exp " << f_exp << " E_exp " << E_exp << " dtheta " << dtheta << " dphi " << dphi << " cell dtheta " << celldtheta << " cell dr " << celldr << " cell dz " << celldz << " cell dphi " << celldphi );
+// ATH_MSG_DEBUG( " sigma_Noise totalNoiseRMS from Alan " << sigma_NoiseA << " sigma_Noise getEffectiveSigma " << sigma_Noise);
+
+// Use expected energy if measured energy below noise level
+// - this will bias the energy measurement towards too high values
+// - the bias is corrected in the thresholdCorrection function below
+//
+// sum measured, expected energies for crossed cells after noise cuts
+//
+ if(cellSampling == CaloSampling::PreSamplerB || cellSampling == CaloSampling::PreSamplerE) {
+ if(f_exp>0&&cellEn>m_sigmasAboveNoise*sigma_Noise&&!badCell) E_em1 += cellEn;
+ }
+ if(curr_cell->caloDDE()->getSubCalo() == CaloCell_ID::LAREM) {
+ E_em_exptot += scale_Ionization*f_exp*E_exp;
+ if(f_exp>0) E_em_expall += scale_Ionization*E_exp;
+ if(f_exp>0&&cellEn>m_sigmasAboveNoise*sigma_Noise&&!badCell) {
+ E_em += cellEn;
+ E_em_expected += scale_Ionization*f_exp*E_exp;
+ ATH_MSG_VERBOSE( " EM cell " << cellEn << " sigma_Noise " << sigma_Noise << " f_exp " << f_exp << " E_exp " << E_exp);
+ }
+ if(f_exp>0&&!badCell) nlay_em++;
+ if(f_exp>0&&!badCell) sigma_Noise_em += sigma_Noise;
+ } else if(curr_cell->caloDDE()->getSubCalo() == CaloCell_ID::TILE) {
+ E_tile_exptot += scale_Ionization*f_exp*E_exp;
+ if(f_exp>0) E_tile_expall += scale_Ionization*E_exp;
+ if(f_exp>0&&cellEn>m_sigmasAboveNoise*sigma_Noise&&!badCell) {
+// &&f_exp*E_exp>m_sigmasAboveNoise*sigma_Noise/4) {
+ E_tile += cellEn;
+ E_tile_expected += scale_Ionization*f_exp*E_exp;
+ ATH_MSG_VERBOSE( " Tile cell " << cellEn << " sigma_Noise " << sigma_Noise << " f_exp " << f_exp << " E_exp " << E_exp);
+ }
+ if(f_exp>0&&!badCell) nlay_tile++;
+ if(f_exp>0&&!badCell) sigma_Noise_tile += sigma_Noise;
+ } else if(curr_cell->caloDDE()->getSubCalo() == CaloCell_ID::LARHEC) {
+ E_HEC_exptot += scale_Ionization*f_exp*E_exp;
+ if(f_exp>0) E_HEC_expall += scale_Ionization*E_exp;
+// lower threshold for HEC
+ if(f_exp>0&&cellEn>m_sigmasAboveNoise*sigma_Noise/2.&&!badCell) {
+// &&f_exp*E_exp>m_sigmasAboveNoise*sigma_Noise/4) {
+ E_HEC += cellEn;
+ E_HEC_expected += scale_Ionization*f_exp*E_exp;
+ ATH_MSG_VERBOSE( " HEC cell " << cellEn << " sigma_Noise " << sigma_Noise << " f_exp " << f_exp << " E_exp " << E_exp);
+ }
+ if(f_exp>0&&!badCell) nlay_HEC++;
+ if(f_exp>0&&!badCell) sigma_Noise_HEC += sigma_Noise;
+ }
+ E_expected += scale_Ionization*f_exp*E_exp;
+ }
+ ATH_MSG_DEBUG( " After cuts measured Energies em " << E_em << " tile " << E_tile << " HEC " << E_HEC);
+ ATH_MSG_DEBUG( " After cuts expected Energies em " << E_em_expected << " tile " << E_tile_expected << " HEC " << E_HEC_expected);
+ ATH_MSG_DEBUG( " No cuts with length expected Energies em " << E_em_exptot << " tile " << E_tile_exptot << " HEC " << E_HEC_exptot);
+ ATH_MSG_DEBUG( " No cuts NO length expected Energies em " << E_em_expall << " tile " << E_tile_expall << " HEC " << E_HEC_expall);
+
+// E resolution of calorimeters
+
+ double sigma_em = sqrt( 500.*E_em);
+ double sigma_tile = sqrt(1000.*E_tile);
+ double sigma_HEC = sqrt(1000.*E_HEC);
+
+// go from e.m. scale to Muon Energy scale
+
+ E_em /= m_emipEM;
+ E_tile /= m_emipTile;
+ E_HEC /= m_emipHEC;
+ ATH_MSG_DEBUG( " e.m. to Muon scale measured Energies em " << E_em << " tile " << E_tile << " HEC " << E_HEC);
+
+// also for errors
+
+ sigma_em /= m_emipEM;
+ sigma_tile /= m_emipTile;
+ sigma_HEC /= m_emipHEC;
+
+// Average Noise per layer
+
+ if(nlay_em>0) sigma_Noise_em /= nlay_em;
+ if(nlay_tile>0) sigma_Noise_tile /= nlay_tile;
+ if(nlay_HEC>0) sigma_Noise_HEC /= nlay_HEC;
+ ATH_MSG_DEBUG( " Average Noise em " << sigma_Noise_em << " nlay " << nlay_em << " tile " << sigma_Noise_tile << " nlay " << nlay_tile << " HEC " << sigma_Noise_HEC << " nlay " << nlay_HEC );
+
+// apply energy correction using the expected Eloss for noise cut
+
+ E_em_threshold += thresholdCorrection(E_em,E_em_expected,m_sigmasAboveNoise*sigma_Noise_em/m_emipEM/4);
+ E_tile_threshold += thresholdCorrection(E_tile,E_tile_expected,m_sigmasAboveNoise*sigma_Noise_tile/m_emipTile/4);
+ E_HEC_threshold += thresholdCorrection(E_HEC,E_HEC_expected,m_sigmasAboveNoise*sigma_Noise_HEC/m_emipHEC/8);
+
+ ATH_MSG_DEBUG( " Threshold correction to Energies em " << E_em_threshold << " tile " << E_tile_threshold << " HEC " << E_HEC_threshold);
+
+ ATH_MSG_DEBUG( " total Energies em " << E_em+E_em_threshold << " tile " << E_tile+E_tile_threshold << " HEC " << E_HEC+E_HEC_threshold);
+
+// Eloss in dead material (so everything not accounted for in the measurement)
+
+ if(scale_Ionization*Eloss>E_expected) E_expected = scale_Ionization*Eloss;
+
+ double E_dead = E_expected - E_em_expected - E_tile_expected - E_HEC_expected;
+
+// treatment of FSR
+
+ E_FSR = 0.;
+ double E_measured = 0.;
+ double E_measured_expected = E_em_expected + E_tile_expected + E_HEC_expected;
+// if(E_em*cos(theta)>m_emEtCut&&E_em1>0.15*E_em) {
+ if(E_em*cos(theta)>m_emEtCut) {
+// large e.m. deposit starting in first e.m. layer
+ E_FSR = E_em;
+// do not use measured e.m. energy for muons and use expected (tile and HEC are fine)
+ E = E_em_expected + E_tile + E_tile_threshold + E_HEC + E_HEC_threshold + E_dead;
+ double sigma_expdead = error_ratio*(E_em_expected + E_dead);
+ sigma = sqrt(sigma_tile*sigma_tile + sigma_HEC*sigma_HEC + sigma_expdead*sigma_expdead);
+ E_measured = E_tile + E_tile_threshold + E_HEC + E_HEC_threshold;
+ } else {
+// no FSR
+ E = E_em + E_em_threshold + E_tile + E_tile_threshold + E_HEC + E_HEC_threshold + E_dead;
+ double sigma_threshold = error_ratio*E_dead;
+ sigma = sqrt(sigma_em*sigma_em + sigma_tile*sigma_tile + sigma_HEC*sigma_HEC+sigma_threshold*sigma_threshold);
+ E_measured = E_em + E_em_threshold + E_tile + E_tile_threshold + E_HEC + E_HEC_threshold;
+ }
+
+ ATH_MSG_DEBUG( " Total energy " << E << " sigma " << sigma << " E calo measured in cells " << E_measured << " E calo expected in cells " << E_measured_expected << " E_expected meanIoni from TG " << E_expected );
+//
+// add for validation (temporarily)
+//
+// E_em_meas,E_em_exp,E_tile_meas,E_tile_exp,E_HEC_meas,E_HEC_exp,E_dead_exp
+//
+ E_em_meas = E_em + E_em_threshold;
+ E_em_exp = E_em_expected;
+ E_tile_meas = E_tile + E_tile_threshold;
+ E_tile_exp = E_tile_expected;
+ E_HEC_meas = E_HEC + E_HEC_threshold;
+ E_HEC_exp = E_HEC_expected;
+ E_dead_exp = E_dead;
+
+ } // calculateMuonEnergies
+
+
+ double MuonCaloEnergyTool::thresholdCorrection(double E_observed,double E_expected,double sigma_Noise) const {
+
+//
+// Total energy observed and expected in cells of LAr, Tile, or HEC after 4*sigma_Noise cut
+//
+// returns a correction to the energy based on the expected energy
+//
+// one should so use: Etotal = E_observed + E_dead + ThresholdCorrection(E_observed,E_expected,sigma_Noise)
+
+ Double_t par[5] = {1.33548e+00,-1.00854e+01,5.38111e+01,-9.33379e+00,5.32913e+01};
+
+ double E = E_expected;
+
+ if(E==0) return 0.;
+
+ double x = sigma_Noise/E;
+ if(x>1.) x = 1.;
+
+ // formula is for a 4 sigma Noise cut and obtained by a small simulation programm
+
+ double fraction = (par[0]+par[1]*x+par[2]*x*x)/(1.+par[3]*x+par[4]*x*x);
+
+ // for low x values (low thresholds) the fraction is bigger than 1
+ // (1.33) because the observed energy overestimates
+
+ ATH_MSG_DEBUG( " ThresholdCorrection E " << E << " E_observed not used " << E_observed << " sigma_Noise " << sigma_Noise << " x = sigma_Noise/E " << x << " fraction " << fraction << " E correction " << (1-fraction)*E );
+
+ return (1.-fraction)*E;
+ } // thresholdCorrection
+
+} // end of namespace Trk
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/MuonCaloEnergyTool.h b/Reconstruction/RecoTools/TrackToCalo/src/MuonCaloEnergyTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..1903eaa3e148d834d94b6c98cc82f7c4afa2d323
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/src/MuonCaloEnergyTool.h
@@ -0,0 +1,89 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/***************************************************************************
+MuonCaloEnergyTool.h - Description
+-------------------
+begin : Summer 2014
+authors : Niels van Eldik (CERN PH-ATC)
+***************************************************************************/
+#ifndef MUONCALOENERGYTOOL_H
+#define MUONCALOENERGYTOOL_H
+
+#include "RecoToolInterfaces/IMuonCaloEnergyTool.h"
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+
+#include "ParticleCaloExtension/ParticleCellAssociation.h"
+#include "TrackToCalo/CaloCellSelectorLayerdR.h"
+
+#include "RecoToolInterfaces/IParticleCaloExtensionTool.h"
+#include "RecoToolInterfaces/IParticleCaloCellAssociationTool.h"
+#include "CaloInterface/ICaloNoiseTool.h"
+
+#include "TrkToolInterfaces/ITrackParticleCreatorTool.h"
+
+#include "PathLengthUtils.h"
+
+class ICaloNoiseTool;
+
+namespace Rec {
+ class IParticleCaloCellAssociationTool;
+}
+
+namespace Trk {
+ class IParticleCaloExtensionTool;
+ class ITrackParticleCreatorTool;
+ class CaloExtension;
+}
+
+namespace Rec {
+
+ class MuonCaloEnergyTool : virtual public IMuonCaloEnergyTool, public AthAlgTool {
+ public:
+
+ MuonCaloEnergyTool(const std::string&,const std::string&,const IInterface*);
+
+ virtual ~MuonCaloEnergyTool();
+
+ virtual StatusCode initialize();
+ virtual StatusCode finalize();
+
+ void calculateMuonEnergies( const Trk::Track* trk,
+ double deltaE, double meanIoni, double sigmaIoni,
+ double& E, double& sigma, double& E_FSR, double& E_expected,
+ double &E_em_meas, double &E_em_exp, double &E_tile_meas, double &E_tile_exp,
+ double &E_HEC_meas, double &E_HEC_exp, double &E_dead_exp ) const;
+
+ private:
+
+ StoreGateSvc* m_storeGate;
+ ToolHandle <Trk::IParticleCaloExtensionTool> m_caloExtensionTool; //!< Tool to make the step-wise extrapolation
+ ToolHandle <Rec::IParticleCaloCellAssociationTool> m_caloCellAssociationTool; //!< Tool to make the step-wise extrapolation
+ ToolHandle< Trk::ITrackParticleCreatorTool > m_particleCreator; /**< The CB Particle Creator Tool */
+
+ ToolHandle <ICaloNoiseTool> m_caloNoiseTool; //!< Tool to quantify electronic noise in calorimeter
+
+ // DATA MEMBERS
+ double m_sigmasAboveNoise; // 4.
+ double m_emEtCut; // 2.5 GeV
+ double m_emF1Cut; // 0.15
+ const double m_emipEM; // 0.78
+ const double m_emipTile; // 0.86
+ const double m_emipHEC; // 0.94
+
+ std::string m_indetTrackParticleLocation;
+ const xAOD::TrackParticleContainer* m_indetTrackParticles;
+ std::string m_muonTrackParticleLocation;
+ const xAOD::TrackParticleContainer* m_muonTrackParticles;
+
+ // PRIVATE METHODS
+ double thresholdCorrection(double E_observed,double E_expected,double sigma_Noise) const;
+
+ };
+
+
+}
+
+#endif
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloCellAssociationTool.cxx b/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloCellAssociationTool.cxx
index 0dc158941f1a6eee37c93848ef4a97c65ec4058a..136c0fec5005cf49bdf383d3eab806526d091ada 100644
--- a/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloCellAssociationTool.cxx
+++ b/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloCellAssociationTool.cxx
@@ -6,15 +6,37 @@
// forward declares
#include "RecoToolInterfaces/IParticleCaloExtensionTool.h"
+#include "ParticleCaloExtension/ParticleCellAssociationCollection.h"
-namespace Trk {
+#include "TrackToCalo/CaloCellHelpers.h"
+
+#include "CaloUtils/CaloCellList.h"
+#include "CaloEvent/CaloCellContainer.h"
+#include "TrkCaloExtension/CaloExtension.h"
+#include "TrkCaloExtension/CaloExtensionHelpers.h"
+
+#include "CaloUtils/CaloCellList.h"
+
+#include "xAODTracking/TrackingPrimitives.h"
+#include <math.h>
+
+namespace Rec {
ParticleCaloCellAssociationTool::ParticleCaloCellAssociationTool(const std::string& t, const std::string& n, const IInterface* p )
: AthAlgTool(t,n,p),
- m_caloExtensionTool("Trk::ParticleCaloExtensionTool/ParticleCaloExtensionTool") {
+ m_caloExtensionTool("Trk::ParticleCaloExtensionTool/ParticleCaloExtensionTool"),
+ m_defaultSelector(0.4)
+ {
declareInterface<IParticleCaloCellAssociationTool>(this);
declareProperty("ParticleCaloExtensionTool", m_caloExtensionTool );
+
+ //Default data source for the calocells
+ declareProperty("CaloCellContainer", m_cellContainerName="AllCalo");
+
+ //coneSize for including calo cells around track
+ declareProperty("ConeSize", m_coneSize = 0.1);
+
}
ParticleCaloCellAssociationTool::~ParticleCaloCellAssociationTool() {}
@@ -23,6 +45,8 @@ namespace Trk {
/* Retrieve track extrapolator from ToolService */
ATH_CHECK( m_caloExtensionTool.retrieve() );
+ m_defaultSelector.setConeSize(m_coneSize);
+
return StatusCode::SUCCESS;
}
@@ -30,15 +54,307 @@ namespace Trk {
return StatusCode::SUCCESS;
}
- const xAOD::ParticleCaloExtension* ParticleCaloCellAssociationTool::caloAssociation( xAOD::TrackParticle& trackParticle ) const {
-
+ const CaloCellContainer* ParticleCaloCellAssociationTool::getCellContainer() const {
+
+ const CaloCellContainer* container = 0;
+ //retrieve the cell container
+ if( evtStore()->retrieve(container, m_cellContainerName).isFailure() || !container ) {
+ ATH_MSG_WARNING( "Unable to retrieve the cell container " << m_cellContainerName << " container ptr " << container );
+ return 0;
+ }
+ if( container ) ATH_MSG_DEBUG("Retrieved cell container " << container->size());
+ return container;
+ }
+
+ void ParticleCaloCellAssociationTool::getCellIntersections( const Trk::CaloExtension& extension,
+ const std::vector<const CaloCell*>& cells,
+ ParticleCellAssociation::CellIntersections& cellIntersections ) const {
+
+// use 3D pathlength in cells
+
+ bool use3D = true;
+
+ cellIntersections.reserve(extension.caloLayerIntersections().size()*1.3);
+
+
+ CaloExtensionHelpers::EntryExitLayerMap entryExitLayerMap;
+ CaloExtensionHelpers::entryExitLayerMap( extension, entryExitLayerMap );
+ ATH_MSG_DEBUG("EntryExitLayerMap " << entryExitLayerMap.size() );
+
+ CaloExtensionHelpers::ScalarLayerMap eLossLayerMap, pathLenLayerMap;
+ CaloExtensionHelpers::eLossLayerMap( extension, eLossLayerMap );
+ CaloExtensionHelpers::pathLenLayerMap( extension, pathLenLayerMap );
+
+ ATH_MSG_DEBUG("Getting cells intersections using cells " << cells.size() );
+ for( auto cell : cells ){
+ // get sampling and look up entry/exit points
+ CaloSampling::CaloSample sample = cell->caloDDE()->getSampling();
+
+ auto pos = entryExitLayerMap.find(sample);
+ if( pos == entryExitLayerMap.end() ) continue;
+
+ //// calculate 3D path length
+ double path = 0.;
+
+ double drFix = cell->caloDDE()->dr();
+ double dzFix = cell->caloDDE()->dz();
+// double dphi = cell->caloDDE()->dphi();
+
+ int isample = cell->caloDDE()->getSampling();
+ bool barrel = false;
+ if(cell->caloDDE()->getSubCalo() == CaloCell_ID::TILE) barrel = true;
+ if(sample== CaloSampling::PreSamplerB || sample== CaloSampling::EMB1 || sample== CaloSampling::EMB2 || sample== CaloSampling::EMB3) barrel = true;
+
+
+ double drTG = fabs((pos->second.first-pos->second.second).perp());
+ double dzTG = fabs((pos->second.first-pos->second.second).z());
+
+ if(barrel) ATH_MSG_VERBOSE(" barrel cell sampling " << cell->caloDDE()->getSampling() << " dr " << cell->caloDDE()->dr() << " drTG " << drTG );
+ if(!barrel) ATH_MSG_VERBOSE(" endcap cell sampling " << cell->caloDDE()->getSampling() << " dz " << cell->caloDDE()->dz() << " dzTG " << dzTG );
+//
+// always use fixed values that correspond to the Calorimeter Tracking Geometry
+// these are different from the CaloCell values
+//
+ if(!barrel) dzFix = dzTG;
+ if(barrel) drFix = drTG;
+
+
+ if(drFix==0.) {
+// recalculate the r values from the other cells
+// BUG: extract dr from cell container for sampling 4 5 6 7 needed EME
+// BUG: extract dr from cell container for sampling 8 9 10 11 needed HEC
+ if(cell->caloDDE()->deta()>0) {
+ double dtheta = cell->caloDDE()->deta();
+ double theta = atan2(cell->caloDDE()->r(),cell->z());
+ if( theta+dtheta < M_PI ) {
+ double deta = -log(tan((theta+dtheta)/2.)) + log(tan((theta)/2.));
+ double scale = fabs(deta)/cell->caloDDE()->deta();
+ double dr = fabs(cell->z()*tan(theta+dtheta/scale) - cell->z()*tan(theta));
+ drFix = dr;
+ ATH_MSG_VERBOSE(" FIX cell sampling " << cell->caloDDE()->getSampling() << " deta " << cell->caloDDE()->deta() << " drFix " << drFix);
+ }else{
+ ATH_MSG_WARNING(" FIXR cell sampling failed: theta " << theta << " dtheta " << dtheta << " sum/pi " << (theta+dtheta)/M_PI
+ << " deta " << cell->caloDDE()->deta());
+ }
+// ATH_MSG_VERBOSE(" FIX cell sampling deta " << deta << " dtheta " << dtheta << " scale " << scale << " theta " << theta );
+ } else {
+ double drMin = 100000.;
+ int dscut = 1;
+ if(!barrel) dscut = 0;
+ const CaloCell* cellFound = 0;
+ for( auto celln : cells ){
+ if(cell==celln) continue;
+ if(cell->caloDDE()->getSubCalo() == celln->caloDDE()->getSubCalo()) {
+ int dsample = isample-celln->caloDDE()->getSampling();
+ if(abs(dsample)==dscut) {
+ double drNew = fabs(cell->caloDDE()->r()-celln->caloDDE()->r());
+ if(drNew<1) continue;
+ if(drNew<drMin) {
+ drMin = drNew;
+ cellFound = celln;
+ }
+ }
+ }
+ }
+ drFix = drMin;
+ ATH_MSG_VERBOSE(" Problem cell sampling " << cell->caloDDE()->getSampling() << " x " << cell->caloDDE()->x() << " y " << cell->caloDDE()->y() << " z " << cell->caloDDE()->z() << " dr " << cell->caloDDE()->dr() << " drFix " << drFix << " drTG " << drTG );
+ if(cellFound) ATH_MSG_VERBOSE(" cellFound sampling " << cellFound->caloDDE()->getSampling() << " x " << cellFound->caloDDE()->x() << " y " << cellFound->caloDDE()->y() << " z " << cellFound->caloDDE()->z() << " dr " << cellFound->caloDDE()->dr() << " dscut " << dscut << " drFix " << drFix );
+ }
+ }
+
+ if(dzFix==0.) {
+// recalculate z values from the other cells
+// BUG: extract dz from cell container for sampling 0 1 2 3 needed EMB
+ if(cell->caloDDE()->deta()>0) {
+ double dtheta = cell->caloDDE()->deta();
+ double theta = atan2(cell->caloDDE()->r(),cell->z());
+ if( theta+dtheta < M_PI ) {
+ double deta = -log(tan((theta+dtheta)/2.)) + log(tan((theta)/2.));
+ double scale = fabs(deta)/cell->caloDDE()->deta();
+ double dz = fabs(cell->caloDDE()->r()/tan(theta+dtheta/scale) - cell->caloDDE()->r()/tan(theta));
+ dzFix = dz;
+ }else{
+ ATH_MSG_WARNING(" FIXZ cell sampling failed: theta " << theta << " dtheta " << dtheta << " sum/pi " << (theta+dtheta)/M_PI
+ << " deta " << cell->caloDDE()->deta());
+ }
+ ATH_MSG_VERBOSE(" Fix cell sampling " << cell->caloDDE()->getSampling() << " deta " << cell->caloDDE()->deta() << " dzFix " << dzFix);
+ // ATH_MSG_VERBOSE(" FIX cell sampling deta " << deta << " dtheta " << dtheta << " scale " << scale << " theta " << theta );
+ } else {
+ double dzMin = 100000.;
+ int dscut = 1;
+ if(barrel) dscut = 0;
+ const CaloCell* cellFound = 0;
+ for( auto celln : cells ){
+ if(cell==celln) continue;
+ if(cell->caloDDE()->getSubCalo() == celln->caloDDE()->getSubCalo()) {
+ int isample2 = celln->caloDDE()->getSampling();
+ if(abs(isample-isample2)==dscut) {
+ double dzNew = fabs(cell->caloDDE()->z()-celln->caloDDE()->z());
+ if(dzNew<1) continue;
+ if(dzNew<dzMin) {
+ dzMin = dzNew;
+ cellFound = celln;
+ }
+ }
+ }
+ }
+ dzFix = dzMin;
+ ATH_MSG_VERBOSE(" Problem cell sampling " << cell->caloDDE()->getSampling() << " x " << cell->caloDDE()->x() << " y " << cell->caloDDE()->y() << " z " << cell->caloDDE()->z() << " dz " << cell->caloDDE()->dz() << " dzFix " << dzFix << " dzTG " << dzTG );
+ if(cellFound) ATH_MSG_VERBOSE(" cellFound sampling " << cellFound->caloDDE()->getSampling() << " x " << cellFound->caloDDE()->x() << " y " << cellFound->caloDDE()->y() << " z " << cellFound->caloDDE()->z() << " dz " << cellFound->caloDDE()->dz() << " dscut " << dscut << " dzFix " << dzFix );
+ }
+ }
+
+ if(use3D) {
+ //pathLenUtil.pathInsideCell( *cell, entryExitLayerMap);
+ double pathInMM = pathLenUtil.get3DPathLength(*cell, pos->second.first, pos->second.second,drFix,dzFix);
+ double totpath = (pos->second.first-pos->second.second).mag();
+ path = totpath!=0 ? pathInMM / totpath : 0.;
+ }
+
+ //// calculate 2D path length (method2)
+ double path2 = 0.;
+
+ if(!use3D) path2 = pathInsideCell(*cell, pos->second.first, pos->second.second );
+
+ if( path2 <= 0. && path <= 0. ) continue;
+
+ // auto entrancePair = entryExitLayerMap.find(entranceID);
+ auto eLossPair = eLossLayerMap.find(sample);
+ double eLoss = 0.;
+//
+// DO NOT scale eloss with tracklength in cell just store total expected eloss
+//
+ if( eLossPair != eLossLayerMap.end() ){
+ eLoss = eLossPair->second;
+ } // IF
+
+ ATH_MSG_DEBUG(" PATH3D = " << path << " PATH2D = " << path2 << " eLoss " << eLoss << " cell energy " << (cell)->energy() << " radius " << cell->caloDDE()->r() << " phi " << cell->caloDDE()->phi() << " dr " << cell->caloDDE()->dr() << " dphi " << cell->caloDDE()->dphi() << " x " << cell->caloDDE()->x() << " y " << cell->caloDDE()->y() << " z " << cell->caloDDE()->z() << " dx " << cell->caloDDE()->dx() << " dy " << cell->caloDDE()->dy() << " dz " << cell->caloDDE()->dz() << " volume " << cell->caloDDE()->volume());
+
+ cellIntersections.push_back( std::make_pair(cell, new ParticleCellIntersection( *cell, eLoss, use3D?path:path2) ) );
+ //cellIntersections.push_back( std::make_pair(cell, new ParticleCellIntersection( *cell, eLoss, path/pathLenLayerMap[sample]) ) );
+ //cellIntersections.push_back( std::make_pair(cell, new ParticleCellIntersection( *cell, path2, path/pathLenLayerMap[sample]) ) ); //tmp hack
+ }
+ ATH_MSG_DEBUG(" added cell intersections " << cellIntersections.size() );
+ }
+
+
+ bool ParticleCaloCellAssociationTool::particleCellAssociation( const xAOD::IParticle& particle,
+ const ParticleCellAssociation*& association, float dr,
+ const CaloCellContainer* container, bool useCaching ) const {
+
+
+ ATH_MSG_DEBUG(" particleCellAssociation: ptr " << &particle << " dr " << dr << " useCaching " << useCaching);
+
+ // reset pointer
+ association = 0;
+ // check if link is already there
+ if( useCaching ){
+ if( particle.isAvailable< ParticleCellAssociation* >("cellAssociation") ){
+ ParticleCellAssociation* theAssociation = particle.auxdata< ParticleCellAssociation* >("cellAssociation");
+ if( theAssociation ){
+ // check whether the cached association is from the same container
+ if( container && theAssociation->container() != container ){
+ ATH_MSG_WARNING("Calling cached version with different container pointer");
+ return false;
+ }
+ // check if we need to resize the cone
+ if( dr > theAssociation->associationConeSize() ){
+ std::vector<const CaloCell*> cells;
+ ATH_MSG_DEBUG(" dr larger then cached dr: " << dr << " cached dr " << theAssociation->associationConeSize());
+ associateCells(*theAssociation->container(),theAssociation->caloExtension(),dr,cells);
+ theAssociation->updateData(std::move(cells),dr);
+ }
+ association = theAssociation;
+ ATH_MSG_DEBUG("Found existing calo extension");
+ return true;
+ }
+ }
+ }
+
// get the extrapolation into the calo
- const xAOD::ParticleCaloExtension* caloExtension = m_caloExtensionTool->caloExtension(trackParticle);
- if( !caloExtension ) return 0;
+ const Trk::CaloExtension* caloExtension = 0;
+ if( !m_caloExtensionTool->caloExtension(particle,caloExtension) ) {
+ ATH_MSG_DEBUG("Failed to get calo extension");
+ return false;
+ }
+ if( caloExtension->caloLayerIntersections().empty()){
+ ATH_MSG_DEBUG( "Received a caloExtension object without track extrapolation");
+ return false;
+ }
+
+ //retrieve the cell container if not provided, return false it retrieval failed
+ if( !container && !(container = getCellContainer()) ) {
+ ATH_MSG_DEBUG("Failed to get calo cell container");
+ return false;
+ }
+ std::vector<const CaloCell*> cells;
+ // update cone size in case it is smaller than the default
+ if( dr < m_coneSize ) dr = m_coneSize;
+ associateCells(*container,*caloExtension,dr,cells);
+
+ // get cell intersections
+ ParticleCellAssociation::CellIntersections cellIntersections;
+ getCellIntersections(*caloExtension,cells,cellIntersections);
+
+// for(auto it : cellIntersections){
+// double f_exp = (it.second)->pathLength();
+// double E_exp = (it.second)->expectedEnergyLoss();
+// ATH_MSG_DEBUG( " path " << f_exp << " expected Eloss " << E_exp );
+// }
+
+ ParticleCellAssociation* theAssocation = new ParticleCellAssociation( *caloExtension, std::move(cells), dr,
+ std::move(cellIntersections), container );
+
+ // now add the extension to the output collection so we are not causing any leaks
+ ParticleCellAssociationCollection* collection = 0;
+ if( !evtStore()->contains<ParticleCellAssociationCollection>(m_cellContainerName) ){
+ collection = new ParticleCellAssociationCollection();
+ if( evtStore()->record( collection, m_cellContainerName).isFailure() ) {
+ ATH_MSG_WARNING( "Failed to record output collection, will leak the ParticleCaloExtension");
+ delete collection;
+ collection = 0;
+ }
+ }else{
+ if(evtStore()->retrieve(collection,m_cellContainerName).isFailure()) {
+ ATH_MSG_WARNING( "Unable to retrieve " << m_cellContainerName << " will leak the ParticleCaloExtension" );
+ }
+ }
+ if( collection ) collection->push_back(theAssocation);
+ else{
+ ATH_MSG_WARNING( "No ParticleCaloCellAssociationCollection, failing extension to avoid memory leak");
+ delete theAssocation;
+ theAssocation = 0;
+ }
+
+ association = theAssocation;
+ if( useCaching ) particle.auxdecor< ParticleCellAssociation* >("cellAssociation") = theAssocation;
+
+
+ return true;
+
+ }
+
+ void ParticleCaloCellAssociationTool::associateCells( const CaloCellContainer& container,
+ const Trk::CaloExtension& caloExtension,
+ float dr,
+ std::vector<const CaloCell*>& cells ) const {
+
+ const Trk::TrackParameters* pars = caloExtension.caloEntryLayerIntersection();
+ if(!pars) {
+ ATH_MSG_WARNING( " NO TrackParameters caloExtension.caloEntryLayerIntersection() ");
+ return;
+ }
+
+ double eta = pars->position().eta();
+ double phi = pars->position().phi();
+
+ // Use Calorimeter list for CPU reasons
+ CaloCellList myList(&container);
+ myList.select(eta,phi,dr);
+ cells.reserve(myList.ncells());
+ cells.insert(cells.end(),myList.begin(),myList.end());
- // now collect the cells
-
- return caloExtension;
+ ATH_MSG_DEBUG("associated cells " << cells.size() << " using cone " << dr );
}
} // end of namespace Trk
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloCellAssociationTool.h b/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloCellAssociationTool.h
index 3fec827c3f93a926c431d6b5bc7a1bc928e089cc..a00fcacd3a6a42c90d089fd1809c6971b6c67cf4 100644
--- a/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloCellAssociationTool.h
+++ b/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloCellAssociationTool.h
@@ -15,12 +15,17 @@ authors : Niels van Eldik (CERN PH-ATC)
#include "AthenaBaseComps/AthAlgTool.h"
#include "GaudiKernel/ToolHandle.h"
-#include "xAODTracking/TrackParticle.h"
+#include "ParticleCaloExtension/ParticleCellAssociation.h"
+#include "TrackToCalo/CaloCellSelectorLayerdR.h"
+#include "PathLengthUtils.h"
namespace Trk {
-
class IParticleCaloExtensionTool;
+ class CaloExtension;
+}
+
+namespace Rec {
class ParticleCaloCellAssociationTool : virtual public IParticleCaloCellAssociationTool, public AthAlgTool {
public:
@@ -32,15 +37,44 @@ namespace Trk {
virtual StatusCode initialize();
virtual StatusCode finalize();
- /** Method to dress a TrackParticle with the calo cells crossed by its track
- @param trackParticle
- @return ParticleCaloExtension, caller does not obtain ownership of the pointer and should not delete it
+ /** Method to get the ParticleCellAssociation of a given TrackParticle
+ @param trackParticle input track particle
+ @param extension reference to a pointer to a ParticleCellAssociation, will be updated if call is successfull
+ NEVER delete the pointer, you will cause a crash!
+ @param dr cone size used for the association
+ If caching is enabled, the cells associated to the association contain at least all cells
+ in dr but could contain more. Users ALWAYS have to recalculate the associated cells in their cone.
+ @param container cell container to be used if provided
+ @param useCaching configure whether the tool caches the result on the track particle
+ The default behavior is 'true' to ensure optimal performance
+ If caching is enabled, the code will perform a consistency check on the container pointer
+ If the function is called twice on the same particle with different containers, the call will fail.
+ The same is true if the function is called once without container and once with on the same particle.
+ @return true if the call was successful
*/
- const xAOD::ParticleCaloExtension* caloAssociation( xAOD::TrackParticle& trackParticle ) const;
-
+ bool particleCellAssociation( const xAOD::IParticle& particle, const ParticleCellAssociation*& association, float dr,
+ const CaloCellContainer* container = 0, bool useCaching = true ) const final;
+
private:
- ToolHandle< IParticleCaloExtensionTool > m_caloExtensionTool;
+
+ void getCellIntersections( const Trk::CaloExtension& caloExtension,
+ const std::vector<const CaloCell*>& cells,
+ ParticleCellAssociation::CellIntersections& cellIntersections) const;
+
+ const CaloCellContainer* getCellContainer() const;
+ void associateCells( const CaloCellContainer& container, const Trk::CaloExtension& caloExtension, float dr,
+ std::vector<const CaloCell*>& cells ) const;
+
+ ToolHandle< Trk::IParticleCaloExtensionTool > m_caloExtensionTool;
+ std::string m_cellContainerName;
+ double m_coneSize;
+ mutable Trk::CaloCellSelectorLayerdR m_defaultSelector;
+
+ mutable PathLengthUtils pathLenUtil;
+
};
+
+
}
#endif
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloClusterAssociationTool.cxx b/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloClusterAssociationTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..ab136705f51f28206fb1de1620ad20ca12e38cb9
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloClusterAssociationTool.cxx
@@ -0,0 +1,159 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "ParticleCaloClusterAssociationTool.h"
+// forward declares
+#include "RecoToolInterfaces/IParticleCaloExtensionTool.h"
+
+#include "ParticleCaloExtension/ParticleClusterAssociationCollection.h"
+#include "ParticlesInConeTools/ICaloClustersInConeTool.h"
+
+#include "TrkCaloExtension/CaloExtension.h"
+#include "TrkCaloExtension/CaloExtensionHelpers.h"
+
+#include "FourMomUtils/P4Helpers.h"
+
+namespace Rec {
+
+ ParticleCaloClusterAssociationTool::ParticleCaloClusterAssociationTool(const std::string& t, const std::string& n, const IInterface* p )
+ : AthAlgTool(t,n,p),
+ m_caloExtensionTool("Trk::ParticleCaloExtensionTool/ParticleCaloExtensionTool"),
+ m_clustersInConeTool("xAOD::CaloClustersInConeTool/CaloClustersInConeTool")
+ {
+
+ declareInterface<IParticleCaloClusterAssociationTool>(this);
+ declareProperty("ParticleCaloExtensionTool", m_caloExtensionTool );
+ declareProperty("ClustersInConeTool", m_clustersInConeTool);
+
+ //coneSize for including calo cells around track
+ declareProperty("ConeSize", m_coneSize = 0.1);
+
+ }
+
+ ParticleCaloClusterAssociationTool::~ParticleCaloClusterAssociationTool() {}
+
+ StatusCode ParticleCaloClusterAssociationTool::initialize() {
+ /* Retrieve track extrapolator from ToolService */
+ ATH_CHECK( m_caloExtensionTool.retrieve() );
+
+ if (!m_clustersInConeTool.empty()) ATH_CHECK(m_clustersInConeTool.retrieve());
+
+ return StatusCode::SUCCESS;
+ }
+
+ StatusCode ParticleCaloClusterAssociationTool::finalize() {
+ return StatusCode::SUCCESS;
+ }
+
+ bool ParticleCaloClusterAssociationTool::particleClusterAssociation( const xAOD::IParticle& particle, const ParticleClusterAssociation*& association, float dr,
+ const xAOD::CaloClusterContainer* container, bool useCaching ) const {
+
+
+ ATH_MSG_DEBUG(" particleCellAssociation: ptr " << &particle << " dr " << dr << " useCaching " << useCaching);
+
+ // reset pointer
+ association = 0;
+ // check if link is already there
+ if( useCaching ){
+ if( particle.isAvailable< ParticleClusterAssociation* >("clusterAssociation") ){
+ ParticleClusterAssociation* theAssociation = particle.auxdata< ParticleClusterAssociation* >("clusterAssociation");
+ if( theAssociation ){
+ // check whether the cached association is from the same container
+ if( container && theAssociation->container() != container ){
+ ATH_MSG_WARNING("Calling cached version with different container pointer");
+ return false;
+ }
+ // check if we need to resize the cone
+ if( dr > theAssociation->associationConeSize() ){
+ ATH_MSG_DEBUG(" dr larger then cached dr: " << dr << " cached dr " << theAssociation->associationConeSize());
+ ParticleClusterAssociation::Data clusters;
+ associateClusters(container,theAssociation->caloExtension(),dr,clusters);
+ theAssociation->updateData(std::move(clusters),dr);
+ }
+ association = theAssociation;
+ ATH_MSG_DEBUG("Found existing calo extension");
+ return true;
+ }
+ }
+ }
+
+ // get the extrapolation into the calo
+ const Trk::CaloExtension* caloExtension = 0;
+ if( !m_caloExtensionTool->caloExtension(particle,caloExtension) ) {
+ ATH_MSG_DEBUG("Failed to get calo extension");
+ return false;
+ }
+ if( caloExtension->caloLayerIntersections().empty()){
+ ATH_MSG_DEBUG( "Received a caloExtension object without track extrapolation");
+ return false;
+ }
+
+
+ // update cone size in case it is smaller than the default
+ if( dr < m_coneSize ) dr = m_coneSize;
+ ParticleClusterAssociation::Data clusters;
+ associateClusters(container,*caloExtension,dr,clusters);
+
+ ParticleClusterAssociation* theAssocation = new ParticleClusterAssociation( *caloExtension, std::move(clusters), dr, container );
+
+ // now add the extension to the output collection so we are not causing any leaks
+ ParticleClusterAssociationCollection* collection = 0;
+ if( !evtStore()->contains<ParticleClusterAssociationCollection>(m_containerName) ){
+ collection = new ParticleClusterAssociationCollection();
+ if( evtStore()->record( collection, m_containerName).isFailure() ) {
+ ATH_MSG_WARNING( "Failed to record output collection, will leak the ParticleCaloExtension");
+ delete collection;
+ collection = 0;
+ }
+ }else{
+ if(evtStore()->retrieve(collection,m_containerName).isFailure()) {
+ ATH_MSG_WARNING( "Unable to retrieve " << m_containerName << " will leak the ParticleCaloExtension" );
+ }
+ }
+ if( collection ) collection->push_back(theAssocation);
+ else{
+ ATH_MSG_WARNING( "No ParticleCaloCellAssociationCollection, failing extension to avoid memory leak");
+ delete theAssocation;
+ theAssocation = 0;
+ }
+
+ association = theAssocation;
+ if( useCaching ) particle.auxdecor< ParticleClusterAssociation* >("clusterAssociation") = theAssocation;
+
+
+ return true;
+
+ }
+
+ void ParticleCaloClusterAssociationTool::associateClusters( const xAOD::CaloClusterContainer* container,
+ const Trk::CaloExtension& caloExtension,
+ float dr,
+ ParticleClusterAssociation::Data& clusters ) const {
+
+ const Trk::TrackParameters* pars = caloExtension.caloEntryLayerIntersection();
+ if(!pars) {
+ ATH_MSG_WARNING( " NO TrackParameters caloExtension.caloEntryLayerIntersection() ");
+ return;
+ }
+
+ float eta = pars->position().eta();
+ float phi = pars->position().phi();
+ if( container ){
+ float dr2Cut = dr*dr;
+ for( unsigned int i=0;i<container->size();++i ){
+
+ float dPhi = P4Helpers::deltaPhi( (*container)[i]->phi(), phi);
+ float dEta = (*container)[i]->eta()-eta;
+ float dr2 = dPhi*dPhi+ dEta*dEta;
+ if( dr2 < dr2Cut ) clusters.push_back( ElementLink<xAOD::CaloClusterContainer>(*container,i) );
+ }
+ }else{
+ if( !m_clustersInConeTool->particlesInCone(eta,phi,dr,clusters) ) {
+ ATH_MSG_WARNING("Failed to get clusters");
+ }
+ }
+ ATH_MSG_DEBUG("associated cells " << clusters.size() << " using cone " << dr );
+ }
+
+} // end of namespace Trk
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloClusterAssociationTool.h b/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloClusterAssociationTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..ed3b22a595acda5dca8d0c6f7bb964cd5110e081
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloClusterAssociationTool.h
@@ -0,0 +1,75 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/***************************************************************************
+ParticleCaloClusterAssociationTool.h - Description
+-------------------
+begin : Summer 2014
+authors : Niels van Eldik (CERN PH-ATC)
+***************************************************************************/
+#ifndef TRACKTOCALO_PARTICLECALOCLUSTERASSOCIATION_H
+#define TRACKTOCALO_PARTICLECALOCLUSTERASSOCIATION_H
+
+#include "RecoToolInterfaces/IParticleCaloClusterAssociationTool.h"
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+
+#include "ParticleCaloExtension/ParticleClusterAssociation.h"
+
+namespace Trk {
+ class IParticleCaloExtensionTool;
+ class CaloExtension;
+}
+namespace xAOD {
+ class ICaloClustersInConeTool;
+}
+namespace Rec {
+
+ class ParticleCaloClusterAssociationTool : virtual public IParticleCaloClusterAssociationTool, public AthAlgTool {
+ public:
+
+ ParticleCaloClusterAssociationTool(const std::string&,const std::string&,const IInterface*);
+
+ virtual ~ParticleCaloClusterAssociationTool();
+
+ virtual StatusCode initialize();
+ virtual StatusCode finalize();
+
+ /** Method to get the ParticleClusterAssociation of a given TrackParticle
+ @param particle input particle
+ @param extension reference to a pointer to a ParticleClusterAssociation, will be updated if call is successfull
+ NEVER delete the pointer, you will cause a crash!
+ @param dr cone size used for the association
+ If caching is enabled, the cells associated to the association contain at least all cells
+ in dr but could contain more. Users ALWAYS have to recalculate the associated cells in their cone.
+ @param container cluster container to be used if provided
+ @param useCaching configure whether the tool caches the result on the track particle
+ The default behavior is 'true' to ensure optimal performance
+ If caching is enabled, the code will perform a consistency check on the container pointer
+ If the function is called twice on the same particle with different containers, the call will fail.
+ The same is true if the function is called once without container and once with on the same particle.
+ @return true if the call was successful
+ */
+ bool particleClusterAssociation( const xAOD::IParticle& particle, const ParticleClusterAssociation*& association, float dr,
+ const xAOD::CaloClusterContainer* container = 0, bool useCaching = true ) const final;
+
+ private:
+
+
+ void associateClusters( const xAOD::CaloClusterContainer* container,
+ const Trk::CaloExtension& caloExtension,
+ float dr,
+ ParticleClusterAssociation::Data& clusters ) const;
+
+ ToolHandle< Trk::IParticleCaloExtensionTool > m_caloExtensionTool;
+ ToolHandle< xAOD::ICaloClustersInConeTool > m_clustersInConeTool;
+ std::string m_containerName;
+ double m_coneSize;
+
+ };
+
+
+}
+
+#endif
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloExtensionTool.cxx b/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloExtensionTool.cxx
index 2eac44ca3cd1204a1807b944c43c102fc8250926..5feb23a00d6d66a299b9d763cdd91a8162b15fd2 100644
--- a/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloExtensionTool.cxx
+++ b/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloExtensionTool.cxx
@@ -4,34 +4,34 @@
#include "ParticleCaloExtensionTool.h"
// forward declares
-#include "TrkExInterfaces/IExtrapolator.h"
-#include "TrkTrack/Track.h"
#include "TrkSurfaces/PerigeeSurface.h"
#include "TrkTrack/TrackStateOnSurface.h"
+#define private public
#include "TrkParameters/TrackParameters.h"
-
+#include "TrkTrack/Track.h"
+#include "TrkExInterfaces/IExtrapolator.h"
+#include "TrkCaloExtension/CaloExtension.h"
+#include "TrkCaloExtension/CaloExtensionCollection.h"
+#define public private
#include "xAODTracking/TrackingPrimitives.h"
#include "AtlasDetDescr/AtlasDetectorID.h"
-#include "MagFieldInterfaces/IMagFieldSvc.h"
-#include "EventPrimitives/EventPrimitivesHelpers.h"
-#include "TrkEventPrimitives/PropDirection.h"
-#include "TrkEventPrimitives/CurvilinearUVT.h"
-#include "TrkEventPrimitives/JacobianLocalToCurvilinear.h"
+#include "TrkParametersIdentificationHelpers/TrackParametersIdHelper.h"
+#include "xAODTruth/TruthVertex.h"
+#include "xAODMuon/Muon.h"
namespace Trk
{
ParticleCaloExtensionTool::ParticleCaloExtensionTool(const std::string& t, const std::string& n, const IInterface* p )
: AthAlgTool(t,n,p),
m_extrapolator("Trk::Extrapolator/AtlasExtrapolator"),
- m_magFieldSvc("AtlasFieldSvc", n),
- m_particleType(Trk::muon)
+ m_particleType(muon)
{
declareInterface<IParticleCaloExtensionTool>(this);
declareProperty("Extrapolator", m_extrapolator );
- declareProperty("MagFieldSvc", m_magFieldSvc);
declareProperty("ParticleType", m_particleTypeName = "muon" );
declareProperty("OutputContainerName", m_containerName = n );
+ declareProperty("StartFromPerigee", m_startFromPerigee = false);
}
ParticleCaloExtensionTool::~ParticleCaloExtensionTool() {}
@@ -40,11 +40,10 @@ namespace Trk
/* Retrieve track extrapolator from ToolService */
ATH_CHECK(detStore()->retrieve(m_detID, "AtlasID" ));
ATH_CHECK( m_extrapolator.retrieve() );
- ATH_CHECK( m_magFieldSvc.retrieve() );
- if( m_particleTypeName == "nonInteracting" ) m_particleType = Trk::nonInteracting;
- else if( m_particleTypeName == "muon" ) m_particleType = Trk::muon;
- else if( m_particleTypeName == "pion" ) m_particleType = Trk::pion;
+ if( m_particleTypeName == "nonInteracting" ) m_particleType = nonInteracting;
+ else if( m_particleTypeName == "muon" ) m_particleType = muon;
+ else if( m_particleTypeName == "pion" ) m_particleType = pion;
else ATH_MSG_WARNING("Unsupported particle type, using muon " << m_particleTypeName );
ATH_MSG_INFO(" Using particle type " << m_particleTypeName << " enum value " << m_particleType );
return StatusCode::SUCCESS;
@@ -54,26 +53,132 @@ namespace Trk
return StatusCode::SUCCESS;
}
- const xAOD::ParticleCaloExtension* ParticleCaloExtensionTool::caloExtension( xAOD::TrackParticle& trackParticle ) const {
+ const xAOD::TrackParticle* ParticleCaloExtensionTool::getTrackParticle(const xAOD::IParticle& particle ) const {
+
+ const xAOD::TrackParticle* trackParticle = dynamic_cast< const xAOD::TrackParticle*>(&particle);
+ if( trackParticle ) return trackParticle;
+
+ const xAOD::Muon* muon = dynamic_cast< const xAOD::Muon*>(&particle);
+ if( muon && muon->primaryTrackParticle() ) return muon->primaryTrackParticle();
+ return nullptr;
+ }
+
+ bool ParticleCaloExtensionTool::caloExtension( const xAOD::IParticle& particle, const Trk::CaloExtension*& extension,
+ bool useCaching ) const {
+ // reset input
+ extension = 0;
+
+ ATH_MSG_DEBUG(" caloExtension: ptr " << &particle << " useCaching " << useCaching);
+
// check if link is already there
- if( trackParticle.isAvailable< xAOD::ParticleCaloExtension* >("caloExtension") ){
- xAOD::ParticleCaloExtension* extension = trackParticle.auxdata< xAOD::ParticleCaloExtension* >("caloExtension");
+ if( useCaching && particle.isAvailable< CaloExtension* >("caloExtension") ){
+ extension = particle.auxdata< CaloExtension* >("caloExtension");
if( extension ){
ATH_MSG_DEBUG("Found existing calo extension");
- return trackParticle.auxdata< xAOD::ParticleCaloExtension* >("caloExtension");
+ return true;
}
}
+
+ // work out the type of particle and get the extension
+ CaloExtension* theExtension = 0;
+ const xAOD::TrackParticle* trackParticle = getTrackParticle(particle);
+ if( trackParticle ) theExtension = caloExtension(*trackParticle);
+ else{
+ const xAOD::NeutralParticle* neutralParticle = dynamic_cast< const xAOD::NeutralParticle*>(&particle);
+ if( neutralParticle ) theExtension = caloExtension(*neutralParticle);
+ else{
+ const xAOD::TruthParticle* truthParticle = dynamic_cast< const xAOD::TruthParticle*>(&particle);
+ if( truthParticle ) theExtension = caloExtension(*truthParticle);
+ else{
+ ATH_MSG_WARNING("Unsupported IParticle type");
+ return false;
+ }
+ }
+ }
+
+ // create and assign extension dressing
+ if( useCaching ) particle.auxdecor< CaloExtension* >("caloExtension") = theExtension;
+ extension = theExtension;
+
+ // return false is extension failed
+ return extension != nullptr;
+ }
+
+ Trk::CaloExtension* ParticleCaloExtensionTool::caloExtension( const xAOD::TruthParticle& particle ) const {
+
+ // get particle type
+ ParticleHypothesis particleType = muon; //nonInteracting;
+ if( abs(particle.pdgId()) == 11 ) particleType = muon; // we dont want the electron to loose energy when extrpolating in the calo
+ else if( abs(particle.pdgId()) == 13 ) particleType = muon;
+
+ // get start parameters
+ const xAOD::TruthVertex* pvtx = particle.prodVtx();
+ if ( pvtx == 0 ) return nullptr;
+
+ double charge = particle.charge();
+ Amg::Vector3D pos( pvtx->x() , pvtx->y() , pvtx->z() );
+ Amg::Vector3D mom( particle.px() , particle.py() , particle.pz() );
+
+ // hack, extrapolate neutral particles as infinit momentum for the time being
+ if( particle.isNeutral() ){
+ charge = 1.;
+ mom.normalize();
+ mom *= 1e10;
+ }
+ Trk::CurvilinearParameters startPars(pos,mom,charge);
+
+ // get extension
+ return caloExtension( startPars, alongMomentum, particleType );
+ }
+
+ Trk::CaloExtension* ParticleCaloExtensionTool::caloExtension( const xAOD::NeutralParticle& particle ) const {
+
+ // create start parameters
+ const Trk::NeutralPerigee& perigee = particle.perigeeParameters();
- // if not create it
- xAOD::ParticleCaloExtension*& caloExtension = trackParticle.auxdata< xAOD::ParticleCaloExtension* >("caloExtension");
+ double charge = 1.;
+ Amg::Vector3D pos( perigee.position() );
+ Amg::Vector3D mom( perigee.momentum() );
+ // hack, extrapolate neutral particles as infinit momentum for the time being
+ mom.normalize();
+ mom *= 1e10;
+ Trk::CurvilinearParameters startPars(pos,mom,charge);
- // check if there is a track
- if( !trackParticle.track() ){
- ATH_MSG_WARNING("No valid track, failing extension" );
- return 0;
+ // get extension
+ return caloExtension( startPars, alongMomentum, muon );
+ }
+
+ Trk::CaloExtension* ParticleCaloExtensionTool::caloExtension( const xAOD::TrackParticle& particle ) const {
+
+ //Determine if the track was fit electron hypothesis -- so extrapolate as if the particles is non interacting
+ ParticleHypothesis particleType = m_particleType;
+ if( particle.particleHypothesis() == xAOD::electron )
+ {
+ ATH_MSG_DEBUG("Fitting using electron hypothesis");
+ particleType = muon;//nonInteracting;
+
+ if(!m_startFromPerigee){
+ unsigned int index(0);
+ if (!particle.indexOfParameterAtPosition(index, xAOD::LastMeasurement)) {
+ ATH_MSG_WARNING("No TrackParticle or no have first measurement");
+ } else {
+ return caloExtension(particle.curvilinearParameters(index),alongMomentum,particleType);
+ }
+ }
+ }
+
+ if(m_startFromPerigee || !particle.track())
+ {
+ bool idExit = true;
+// Muon Entry is around z 6783 and r 4255
+ if(fabs(particle.perigeeParameters().position().z())>6700.) idExit = false;
+ if(particle.perigeeParameters().position().perp()>4200.) idExit = false;
+ PropDirection propDir = idExit ? alongMomentum : oppositeMomentum;
+ return caloExtension(particle.perigeeParameters(),propDir,particleType);
}
- const Trk::Track& track = *trackParticle.track();
+
+ const Track& track = *particle.track();
// look-up the parameters closest to the calorimeter in ID and muon system
ATH_MSG_DEBUG("trying to add calo layers" );
@@ -90,83 +195,101 @@ namespace Trk
// require at least one of them
if( !idExitParamers && !muonEntryParamers) {
- ATH_MSG_WARNING("Failed to find ID and Muon parameters");
- return 0;
+ idExitParamers = track.perigeeParameters();
}
// pick start parameters, start in ID if possible
const TrackParameters* startPars = idExitParamers ? idExitParamers : muonEntryParamers;
+ if( !startPars ){
+ ATH_MSG_WARNING("Failed to find start parameters");
+ return nullptr;
+ }
PropDirection propDir = idExitParamers ? alongMomentum : oppositeMomentum;
+
+ return caloExtension(*startPars,propDir,particleType);
+ }
+
- ATH_MSG_DEBUG("looking up calo states: r " << startPars->position().perp() << " z " << startPars->position().z()
- << " momentum " << startPars->momentum().mag() );
+ Trk::CaloExtension* ParticleCaloExtensionTool::caloExtension( const TrackParameters& startPars, PropDirection propDir,
+ ParticleHypothesis particleType ) const {
+
+ ATH_MSG_DEBUG("looking up calo states: r " << startPars.position().perp() << " z " << startPars.position().z()
+ << " momentum " << startPars.momentum().mag() );
// pointers to hold results and go
- std::vector<const TrackStateOnSurface*>* material = 0;
- const std::vector< std::pair< const TrackParameters*, int > >* caloParameters =
- m_extrapolator->extrapolate( *startPars, propDir, m_particleType, material, 3 );
- if( !caloParameters ) return 0;
+ std::vector<const TrackStateOnSurface*>* material = 0;//new std::vector<const TrackStateOnSurface*>();
+ const auto* caloParameters = m_extrapolator->extrapolate( startPars, propDir, particleType, material, 3 );
+ if( material ) {
+ ATH_MSG_DEBUG("Got material " << material->size() );
+ for( auto& m : *material ) {
+ if( msgLvl(MSG::DEBUG) ){
+ msg(MSG::DEBUG) << " layer ";
+ const TrackParameters* param = m->trackParameters();
+ if( param ) msg(MSG::DEBUG) << " param " << param << " pos: r " << param->position().perp() << " z " << param->position().z()
+ << " pt " << param->momentum().perp();
+ const Trk::MaterialEffectsBase* mat = m->materialEffectsOnTrack();
+ if( mat ) msg(MSG::DEBUG) << " mat: " << mat->thicknessInX0();
+ msg(MSG::DEBUG) << endreq;
+ }
+ delete m;
+ }
+ delete material;
+ }
+ if( !caloParameters ) return nullptr;
+ TrackParametersIdHelper parsIdHelper;
// create final object
- std::vector< std::vector<float> > parameters;
- std::vector< std::vector<float> > parametersCovariance;
- std::vector< int > identifiers;
- parameters.reserve(caloParameters->size());
- parametersCovariance.reserve(caloParameters->size());
+ const TrackParameters* caloEntry = 0;
+ const TrackParameters* muonEntry = 0;
+ std::vector<const CurvilinearParameters*> caloLayers;
+ caloLayers.reserve(caloParameters->size()-1);
ATH_MSG_DEBUG( " Found calo parameters: " << caloParameters->size() );
for( const auto& p : *caloParameters ){
- const Trk::TrackParameters* param = p.first;
+ const TrackParameters* param = p.first;
if( !param ) continue;
- ATH_MSG_DEBUG( " id " << p.second << " pos: r " << param->position().perp() << " z " << param->position().z()
- << " momentum " << param->momentum().perp() << " cov " << param->covariance() );
-
- std::vector<float> pars(6);
- pars[0] = param->position().x();
- pars[1] = param->position().y();
- pars[2] = param->position().z();
- pars[3] = param->momentum().x();
- pars[4] = param->momentum().y();
- pars[5] = param->momentum().z();
- parameters.push_back(pars);
- identifiers.push_back(p.second);
-
- std::vector<float> covVec;
- if( param->covariance() ){
- //now convert from to Curvilinear -- to be double checked for correctness
- Amg::Vector3D magnFieldVect; magnFieldVect.setZero();
- m_magFieldSvc->getField( ¶m->position(), &magnFieldVect );
-
- Trk::CurvilinearUVT curvilinearUVT(param->momentum());
- const Amg::Transform3D& localToGlobalTransform = param->associatedSurface().transform();
-
- Trk::JacobianLocalToCurvilinear jacobian(magnFieldVect,
- param->parameters()[Trk::qOverP],
- sin(param->parameters()[Trk::theta]),
- curvilinearUVT,
- localToGlobalTransform.rotation().col(0),
- localToGlobalTransform.rotation().col(1));
-
- AmgSymMatrix(5) covarianceMatrix = param->covariance()->similarity(jacobian);
- Amg::compress(covarianceMatrix,covVec);
- }
- parametersCovariance.push_back(covVec);
- delete param;
+ ATH_MSG_DEBUG( " param " << param << " id " << p.second /*<< " type " << p.first->type()*/ << " pos: r " << param->position().perp() << " z " << param->position().z()
+ << " pt " << param->momentum().perp() << " cov " << param->covariance() );
+
+ // assign parameters
+ if( p.second == 1 && propDir == Trk::alongMomentum) caloEntry = p.first;
+ else if( p.second == 3 && propDir == Trk::oppositeMomentum) caloEntry = p.first;
+ else if( p.second == 3 && propDir == Trk::alongMomentum) muonEntry = p.first;
+ else if( p.second == 4 && propDir == Trk::oppositeMomentum) muonEntry = p.first;
+ else{
+ bool isEntry = p.second > 0 ? true : false;
+ TrackParametersIdentifier id = parsIdHelper.encode( AtlasDetDescr::fFirstAtlasCaloTechnology,
+ static_cast<CaloSampling::CaloSample>( abs(p.second)%1000 ),
+ isEntry );
+ const CurvilinearParameters* cpars = dynamic_cast<const CurvilinearParameters*>(p.first);
+ if( !cpars ){ /*p.first->type() == Curvilinear*/
+ cpars = new CurvilinearParameters(p.first->position(),p.first->momentum(),p.first->charge(),nullptr,id);
+ delete p.first;
+ }else{
+ const_cast<CurvilinearParameters*>(cpars)->m_cIdentifier = id;
+ }
+ caloLayers.push_back( cpars );
+ }
+ }
+
+ if(!muonEntry && propDir == Trk::oppositeMomentum && fabs(startPars.position().perp()-4255.)<1.) {
+// muonEntry is right at the startPars position
+ muonEntry = &startPars;
+ }
+
+ if( muonEntry ) {
+ if( muonEntry->covariance() ) {
+ ATH_MSG_VERBOSE (" p at MuonEntry " << muonEntry->momentum().mag() << " cov 00 " << (*(muonEntry->covariance()))(0,0) << " cov 11 " << (*(muonEntry->covariance()))(1,1));
+ }
}
-
// clean-up memory
delete caloParameters;
- if( material ) {
- for( auto& m : *material ) delete m;
- delete material;
- }
- caloExtension = new xAOD::ParticleCaloExtension(trackParticle.charge(),std::move(parameters),
- std::move(parametersCovariance),std::move(identifiers));
+ CaloExtension* theExtension = new CaloExtension(caloEntry,muonEntry,std::move(caloLayers));
// now add the extension to the output collection so we are not causing any leaks
- ParticleCaloExtensionCollection* collection = 0;
- if( !evtStore()->contains<ParticleCaloExtensionCollection>(m_containerName) ){
- collection = new ParticleCaloExtensionCollection();
+ CaloExtensionCollection* collection = 0;
+ if( !evtStore()->contains<CaloExtensionCollection>(m_containerName) ){
+ collection = new CaloExtensionCollection();
if( evtStore()->record( collection, m_containerName).isFailure() ) {
ATH_MSG_WARNING( "Failed to record output collection, will leak the ParticleCaloExtension");
delete collection;
@@ -175,12 +298,17 @@ namespace Trk
}else{
if(evtStore()->retrieve(collection,m_containerName).isFailure()) {
ATH_MSG_WARNING( "Unable to retrieve " << m_containerName << " will leak the ParticleCaloExtension" );
+ collection = 0;
}
}
- if( collection ) collection->push_back(caloExtension);
-
- return caloExtension;
+ if( collection ) collection->push_back(theExtension);
+ else{
+ ATH_MSG_WARNING( "No CaloExtension Collection, failing extension to avoid memory leak");
+ delete theExtension;
+ theExtension = 0;
+ }
+ return theExtension;
}
-
-
} // end of namespace Trk
+
+
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloExtensionTool.h b/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloExtensionTool.h
index 7f45eaae5037b7292b9fd49886075755935db5e1..411ee8869eec456b6b0944047ef4995574b52984 100644
--- a/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloExtensionTool.h
+++ b/Reconstruction/RecoTools/TrackToCalo/src/ParticleCaloExtensionTool.h
@@ -11,26 +11,26 @@ authors : Niels van Eldik (CERN PH-ATC)
#ifndef TRKPARTICLECREATOR_PARTICLECALOEXTENSIONTOOL_H
#define TRKPARTICLECREATOR_PARTICLECALOEXTENSIONTOOL_H
-#include "RecoToolInterfaces/IParticleCaloExtensionTool.h"
+
#include "AthenaBaseComps/AthAlgTool.h"
#include "GaudiKernel/ToolHandle.h"
#include "GaudiKernel/ServiceHandle.h"
+#define private public
+#include "RecoToolInterfaces/IParticleCaloExtensionTool.h"
#include "xAODTracking/TrackParticle.h"
-#include "xAODTracking/ParticleCaloExtension.h"
+#define public private
#include "TrkEventPrimitives/ParticleHypothesis.h"
+#include "xAODTracking/NeutralParticle.h"
+#include "xAODTruth/TruthParticle.h"
class AtlasDetectorID;
-namespace MagField
-{
- class IMagFieldSvc;
-}
-
namespace Trk {
class IExtrapolator;
+ class CaloExtension;
class ParticleCaloExtensionTool : virtual public IParticleCaloExtensionTool, public AthAlgTool {
public:
@@ -42,20 +42,33 @@ namespace Trk {
virtual StatusCode initialize();
virtual StatusCode finalize();
- /** Method to dress a TrackParticle with the calo layers crossed by its track
- @param trackParticle
- @return ParticleCaloExtension, caller does not obtain ownership of the pointer and should not delete it
+ /** Method to dress a IParticle with the calo layers crossed by its track
+ @param IParticle reference to the particle
+ @param extension reference to a pointer to a CaloExtesion, will be updated if call is successfull
+ NEVER delete the pointer, you will cause a crash!
+ @param useCaching configure whether the tool caches the result on the track particle
+ The default behavior is 'true' to ensure optimal performance
+ @return true if the call was successful
*/
- const xAOD::ParticleCaloExtension* caloExtension( xAOD::TrackParticle& trackParticle ) const;
-
+ bool caloExtension( const xAOD::IParticle& particle, const Trk::CaloExtension*& extension, bool useCaching = true ) const final;
+ Trk::CaloExtension* caloExtension( const TrackParameters& startPars, PropDirection propDir, ParticleHypothesis particleType ) const;
private:
+
+ const xAOD::TrackParticle* getTrackParticle(const xAOD::IParticle& particle ) const;
+
+ Trk::CaloExtension* caloExtension( const xAOD::TruthParticle& particle ) const;
+ Trk::CaloExtension* caloExtension( const xAOD::NeutralParticle& particle ) const;
+ Trk::CaloExtension* caloExtension( const xAOD::TrackParticle& particle ) const;
+
+
+
+
const AtlasDetectorID* m_detID;
ToolHandle< IExtrapolator > m_extrapolator;
- ServiceHandle<MagField::IMagFieldSvc> m_magFieldSvc;
-
- Trk::ParticleHypothesis m_particleType;
- std::string m_particleTypeName;
- std::string m_containerName;
+ ParticleHypothesis m_particleType;
+ std::string m_particleTypeName;
+ std::string m_containerName;
+ bool m_startFromPerigee;
};
}
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/PathLengthUtils.cxx b/Reconstruction/RecoTools/TrackToCalo/src/PathLengthUtils.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..473be9d53b604ae8c58ec231a4bdbdf97db75c38
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/src/PathLengthUtils.cxx
@@ -0,0 +1,782 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "PathLengthUtils.h"
+#include "CaloGeoHelpers/CaloSampling.h"
+#include <iostream>
+
+//=================================
+PathLengthUtils::PathLengthUtils(){
+//=================================
+
+}
+
+//==================================
+PathLengthUtils::~PathLengthUtils(){
+//==================================
+
+}
+
+//=========================================================================================================================
+double PathLengthUtils::get3DPathLength(const CaloCell& cell, const Amg::Vector3D& entrance, const Amg::Vector3D& exit, double drFix, double dzFix){
+//=========================================================================================================================
+
+ int debugLevel = 0;
+
+ if(debugLevel>=1) std::cout << std::endl;
+
+ double dmin = 10.;
+ double dphimin = dmin/cell.caloDDE()->r();
+
+// Get cell parameters
+
+ double dphi = cell.caloDDE()->dphi();
+ if(dphi<dphimin) dphi = dphimin;
+ double CellPhimin = cell.caloDDE()->phi()-dphi*0.5;
+ double CellPhimax = cell.caloDDE()->phi()+dphi*0.5;
+
+// Order according to IP to follow track convention with entry and exit planes
+
+ int isign = 1;
+ if(cell.caloDDE()->z()<0) isign = -1;
+ double dr = cell.caloDDE()->dr();
+//
+// always use fixed values that correspond to the Calorimeter Tracking Geometry
+// these are different from the CaloCell values
+//
+ if(dr==0) dr = drFix;
+ dr = drFix;
+ if(dr<dmin) dr = dmin;
+ double dz = cell.caloDDE()->dz();
+ if(dz==0) dz = dzFix;
+ dz = dzFix;
+ if(dz<dmin) dz = dmin;
+
+ double CellZmin = cell.caloDDE()->z()-isign*dz*0.5;
+ double CellZmax = cell.caloDDE()->z()+isign*dz*0.5;
+ double CellRmin = cell.caloDDE()->r()-dr*0.5;
+ double CellRmax = cell.caloDDE()->r()+dr*0.5;
+
+// fix cells with zero R or Z width
+
+// if(cell.caloDDE()->dr() == 0){
+// // V = dr*dphi*r*dz
+// double CellVolume = cell.caloDDE()->volume();
+// double product = cell.caloDDE()->dphi() * cell.caloDDE()->r() * cell.caloDDE()->dz();
+// double drFormula = dmin;
+// if(product!=0){
+// drFormula = CellVolume / product;
+// }
+// if(drFormula<dmin) drFormula = dmin;
+// CellRmin = cell.caloDDE()->r() - drFormula*0.5;
+// CellRmax = cell.caloDDE()->r() + drFormula*0.5;
+// dr = drFormula;
+// if(debugLevel>=2) std::cout << "drFormula " << drFormula << " new CellRmin " << CellRmin << " new CellRmax " << CellRmax << std::endl;
+// }
+
+// if(cell.caloDDE()->dz() == 0){
+// // V = dr*dphi*r*dz
+// double CellVolume = cell.caloDDE()->volume();
+// double product = cell.caloDDE()->dphi() * cell.caloDDE()->r() * cell.caloDDE()->dr();
+// double dzFormula = dmin;
+// if(product!=0){
+// dzFormula = CellVolume / product;
+// }
+// if(dzFormula<dmin) dzFormula = dmin;
+// CellZmin = cell.caloDDE()->z() - isign*dzFormula*0.5;
+// CellZmax = cell.caloDDE()->z() + isign*dzFormula*0.5;
+// dz = dzFormula;
+// if(debugLevel>=2) std::cout << "dzFormula " << dzFormula << " new CellZmin " << CellZmin << " new CellZmax " << CellZmax << std::endl;
+// }
+//
+// For CPU reason do a fast check on crossing
+//
+// track direction
+//
+ Amg::Vector3D dir = exit - entrance;
+ if(dir.mag()==0) return 0.;
+ dir = dir/dir.mag();
+//
+// point of closest approach in x,y
+//
+ double r0 = (entrance.x()-cell.caloDDE()->x())*sin(dir.phi()) - (entrance.y()-cell.caloDDE()->y())*cos(dir.phi());
+ bool crossing = true;
+ if(fabs(r0)>sqrt(cell.caloDDE()->r()*cell.caloDDE()->r()*dphi*dphi + dr*dr)) crossing = false;
+//
+// point of closest approach in r,z
+//
+ double rz0 = (entrance.perp()-cell.caloDDE()->r())*cos(dir.theta()) - (entrance.z()-cell.caloDDE()->z())*sin(dir.theta());
+ if(fabs(rz0)>sqrt(dr*dr+dz*dz)) crossing = false;
+
+ if(debugLevel>=2&&crossing) {
+ std::cout << " get3DPathMatch distance xy " << r0 << " distance rz " << rz0 << std::endl;
+ std::cout << std::endl;
+ }
+
+ if(!crossing) return 0.;
+
+// construct 8 corners of the volume
+
+ Amg::Vector3D Corner0(CellRmin*cos(CellPhimin),CellRmin*sin(CellPhimin),CellZmin);
+ Amg::Vector3D Corner1(CellRmin*cos(CellPhimax),CellRmin*sin(CellPhimax),CellZmin);
+ Amg::Vector3D Corner2(CellRmin*cos(CellPhimin),CellRmin*sin(CellPhimin),CellZmax);
+ Amg::Vector3D Corner3(CellRmin*cos(CellPhimax),CellRmin*sin(CellPhimax),CellZmax);
+
+ Amg::Vector3D Corner4(CellRmax*cos(CellPhimin),CellRmax*sin(CellPhimin),CellZmin);
+ Amg::Vector3D Corner5(CellRmax*cos(CellPhimax),CellRmax*sin(CellPhimax),CellZmin);
+ Amg::Vector3D Corner6(CellRmax*cos(CellPhimin),CellRmax*sin(CellPhimin),CellZmax);
+ Amg::Vector3D Corner7(CellRmax*cos(CellPhimax),CellRmax*sin(CellPhimax),CellZmax);
+
+ if(debugLevel>=1) {
+ double volume_ratio = cell.caloDDE()->dr()*cell.caloDDE()->r()*cell.caloDDE()->dphi()*cell.caloDDE()->dz() /(cell.caloDDE()->volume() + 1.) ;
+ std::cout << " sampling " << cell.caloDDE()->getSampling() << " CellPhimin " << CellPhimin << " CellPhimax " << CellPhimax << " CellZmin " << CellZmin << " CellZmax " << CellZmax << " CellRmin " << CellRmin << " CellRmax " << CellRmax << " cell.caloDDE()->dz() " << cell.caloDDE()->dz() << " cell.caloDDE()->dr() " << cell.caloDDE()->dr() << " cell.caloDDE()->dx() " << cell.caloDDE()->dx() << " cell.caloDDE()->dy() " << cell.caloDDE()->dy() << " volume_ratio " << volume_ratio << std::endl;
+ }
+ if(debugLevel>=2) {
+ std::cout << " Corner0 x " << Corner0.x() << " y " << Corner0.y() << " z " << Corner0.z() << std::endl;
+ std::cout << " Corner1 x " << Corner1.x() << " y " << Corner1.y() << " z " << Corner1.z() << std::endl;
+ std::cout << " Corner2 x " << Corner2.x() << " y " << Corner2.y() << " z " << Corner2.z() << std::endl;
+ std::cout << " Corner3 x " << Corner3.x() << " y " << Corner3.y() << " z " << Corner3.z() << std::endl;
+ std::cout << " Corner4 x " << Corner4.x() << " y " << Corner4.y() << " z " << Corner4.z() << std::endl;
+ std::cout << " Corner5 x " << Corner5.x() << " y " << Corner5.y() << " z " << Corner5.z() << std::endl;
+ std::cout << " Corner6 x " << Corner6.x() << " y " << Corner6.y() << " z " << Corner6.z() << std::endl;
+ std::cout << " Corner7 x " << Corner7.x() << " y " << Corner7.y() << " z " << Corner7.z() << std::endl;
+ }
+
+// Entry plane vectors through Corner0
+
+ Amg::Vector3D dir01 = Corner1 - Corner0; // direction along Phi
+ Amg::Vector3D dir02 = Corner2 - Corner0; // direction along Z
+// Second (side) plane through Corner0 has dir20 and dir40
+ Amg::Vector3D dir04 = Corner4 - Corner0; // direction along R
+// Third plane through Corner0 has dir10 and dir40
+
+
+// Exit plane vectors through Corner6
+
+ Amg::Vector3D dir67 = Corner7 - Corner6; // direction along Phi
+ Amg::Vector3D dir64 = Corner4 - Corner6; // direction along Z
+// Second (side) plane through Corner dir62 and dir64
+ Amg::Vector3D dir62 = Corner2 - Corner6; // direction along R
+// Third plane through Corner dir67 and dir62
+
+// minus direction of track
+ Amg::Vector3D dirT = exit - entrance;
+ dirT = -dirT;
+
+// check positions
+
+ if(debugLevel>=2) {
+ if(entrance.perp()>Corner0.perp()) {
+ std::cout << " PROBLEM radius entrance larger than cell corner " << std::endl;
+ }
+ if(entrance.mag()>Corner0.mag()||fabs(exit.z())<fabs(Corner0.z())) {
+ std::cout << " PROBLEM 3D or z distance to IP entrance larger than cell corner " << std::endl;
+ }
+
+ if(exit.perp()<Corner6.perp()) {
+ std::cout << " PROBLEM radius exit small than cell corner " << std::endl;
+ }
+ if(entrance.mag()>Corner6.mag()||fabs(exit.z())<fabs(Corner6.z())) {
+ std::cout << " PROBLEM 3D or Z distance to IP exit smaller than cell corner " << std::endl;
+ }
+ }
+
+// dot products of track with plane vectors
+
+ if(debugLevel>=2) {
+ std::cout << " dir01 " << dir01.x() << " y " << dir01.y() << " z " << dir01.z() << " mag " << dir01.mag() << std::endl;
+ std::cout << " dir02 " << dir02.x() << " y " << dir02.y() << " z " << dir02.z() << " mag " << dir02.mag() << std::endl;
+ std::cout << " dir04 " << dir04.x() << " y " << dir04.y() << " z " << dir04.z() << " mag " << dir04.mag() << std::endl;
+ }
+
+ double dotp1 = fabs(dirT.dot(dir01)/dir01.mag()/dirT.mag());
+ double dotp2= fabs(dirT.dot(dir02)/dir02.mag()/dirT.mag());
+ double dotp4 = fabs(dirT.dot(dir04)/dir04.mag()/dirT.mag());
+//
+// order planes according to dotproduct and calculate Matrices
+// mEntry matrix
+// column(0) vector for first plane e.g. dir01
+// column(1) second vector for first plane e.g. dir02
+// column(2) = minus track direction
+// type = 124 here plane made by 1 and 2 = dir01 and dir02
+// later also plane made by 14 = dir01 and dir04 can be tried (second try)
+// or 24 = dir02 and dir04 (third try)
+//
+ Amg::MatrixX mEntry(3,3);
+ mEntry.col(2) = dirT;
+ Amg::MatrixX mExit(3,3);
+ mExit.col(2) = dirT;
+
+ int type = 0;
+ double dotmax = 0.;
+ if(dotp1<=dotp2&&dotp2<=dotp4) {
+ type = 124;
+ dotmax = dotp4;
+ mEntry.col(0) = dir01;
+ mEntry.col(1) = dir02;
+ mExit.col(0) = dir67;
+ mExit.col(1) = dir64;
+ } else if (dotp1<=dotp4&&dotp4<=dotp2) {
+ type = 142;
+ dotmax = dotp2;
+ mEntry.col(0) = dir01;
+ mEntry.col(1) = dir04;
+ mExit.col(0) = dir67;
+ mExit.col(1) = dir62;
+ } else if (dotp2<=dotp1&&dotp1<=dotp4) {
+ type = 214;
+ dotmax = dotp4;
+ mEntry.col(0) = dir02;
+ mEntry.col(1) = dir01;
+ mExit.col(0) = dir64;
+ mExit.col(1) = dir67;
+ } else if (dotp2<=dotp4&&dotp4<=dotp1) {
+ type = 241;
+ dotmax = dotp1;
+ mEntry.col(0) = dir02;
+ mEntry.col(1) = dir04;
+ mExit.col(0) = dir64;
+ mExit.col(1) = dir62;
+ } else if (dotp4<=dotp2&&dotp2<=dotp1) {
+ type = 421;
+ dotmax = dotp1;
+ mEntry.col(0) = dir04;
+ mEntry.col(1) = dir02;
+ mExit.col(0) = dir62;
+ mExit.col(1) = dir64;
+ } else if (dotp4<=dotp1&&dotp1<=dotp2) {
+ type = 412;
+ dotmax = dotp2;
+ mEntry.col(0) = dir04;
+ mEntry.col(1) = dir01;
+ mExit.col(0) = dir62;
+ mExit.col(1) = dir67;
+ }
+
+// position of track w.r.t. all planes
+
+ Amg::Vector3D posEn = entrance-(Corner1+Corner2+Corner4)/2.;
+ Amg::Vector3D posEx = entrance-(Corner2+Corner4+Corner7)/2.;
+//
+// put posEn en posEx in the middle of the plane
+//
+// this makes it possible to check that the Matrix solution lies within the bounds
+// |x|<0.5 and |y|<0.5
+// x corresponds to the first column(0) e.g. dir01
+// y corresponds to the second column(1) e.g. dir02
+//
+ int type0 = type-10*int(type/10);
+ if(type0==1) {
+// last plane 1 not used
+ posEn += Corner1/2.;
+ posEx += Corner7/2.;
+ }
+ if(type0==2) {
+// last plane 2 not used
+ posEn += Corner2/2.;
+ posEx += Corner4/2.;
+ }
+ if(type0==4) {
+// last plane 4 not used
+ posEn += Corner4/2.;
+ posEx += Corner2/2.;
+ }
+
+ if(debugLevel>=2) {
+ std::cout << " not normalized direction track " << dirT.x() << " y " << dirT.y() << " z " << dirT.z() << " eta " << dirT.eta() << " type " << type << " dotmax " << dotmax << std::endl;
+ std::cout << " dotp1 " << dotp1 << " dotp2 " << dotp2 << " dotp4 " << dotp4 << std::endl;
+ std::cout << " entry position " << entrance.x() << " y " << entrance.y() << " z " << entrance.z() << std::endl;
+ std::cout << " mEntry Matrix " << std::endl;
+ std::cout << mEntry << std::endl;
+ std::cout << " exit position " << exit.x() << " y " << exit.y() << " z " << exit.z() << std::endl;
+ std::cout << " mExit Matrix " << std::endl;
+ std::cout << mExit << std::endl;
+ }
+
+// cross with entry plane
+ int typeCheck = type0;
+
+ Amg::Vector3D pathEn;
+ if(debugLevel>=2) std::cout << " Try crossing Entry first try for type " << type << std::endl;
+ bool crossingEn = crossingMatrix(mEntry,posEn,pathEn);
+// if dotmax near to 1 Matrix is singular (track is parallel to the plane)
+ if(debugLevel>=2&&crossingEn) std::cout << " Succesfull crossing Entry first try for type " << type << " type 0 " << type0 << std::endl;
+ if(!crossingEn) {
+ int type1 = int((type-100*int(type/100))/10);
+ int type2 = int(type/100);
+ if(type1==1) posEn += Corner1/2.;
+ if(type1==2) posEn += Corner2/2.;
+ if(type1==4) posEn += Corner4/2.;
+ if(type0==1) {
+ mEntry.col(1) = dir01;
+ posEn -= Corner1/2.;
+ }
+ if(type0==2) {
+ mEntry.col(1) = dir02;
+ posEn -= Corner2/2.;
+ }
+ if(type0==4) {
+ mEntry.col(1) = dir04;
+ posEn -= Corner4/2.;
+ }
+ if(debugLevel>=2) std::cout << " Try crossing Entry second try for type " << type << std::endl;
+ crossingEn = crossingMatrix(mEntry,posEn,pathEn);
+ typeCheck = type1;
+ if(debugLevel>=2&&crossingEn) std::cout << " Succesfull crossing Entry second try for type " <<type << std::endl;
+ if(!crossingEn) {
+ if(type2==1) posEn += Corner1/2.;
+ if(type2==2) posEn += Corner2/2.;
+ if(type2==4) posEn += Corner4/2.;
+ if(type1==1) {
+ mEntry.col(0) = dir01;
+ posEn -= Corner1/2.;
+ }
+ if(type1==2) {
+ mEntry.col(0) = dir02;
+ posEn -= Corner2/2.;
+ }
+ if(type1==4) {
+ mEntry.col(0) = dir04;
+ posEn -= Corner4/2.;
+ }
+ if(debugLevel>=2) std::cout << " Try crossing Entry third try for type " << type << std::endl;
+ crossingEn = crossingMatrix(mEntry,posEn,pathEn);
+ typeCheck = type2;
+ if(debugLevel>=2&&crossingEn) std::cout << " Succesfull crossing Entry third try for type " << type << std::endl;
+ if(debugLevel>=2&&!crossingEn) std::cout << " No crossing Entry third try for type " << type << std::endl;
+ }
+ }
+
+ Amg::Vector3D posXEn = entrance;
+ if(crossingEn) {
+ posXEn = entrance - dirT*pathEn.z();
+ if(typeCheck==1) {
+ double dphiCheck = acos(cos(posXEn.phi())*cos((CellPhimax+CellPhimin)/2.)+sin(posXEn.phi())*sin((CellPhimax+CellPhimin)/2.));
+ if(dphiCheck>fabs(CellPhimax-CellPhimin)) {
+ crossingEn = false;
+ }
+ if(debugLevel>=1&&!crossingEn) std::cout << " check crossing in phi " << crossingEn << " dphiCheck " << dphiCheck << " CellPhimin " << CellPhimin << " CellPhimax " << CellPhimax << " position crossing phi " << posXEn.phi() << std::endl;
+ }
+ if(typeCheck==2) {
+ if(posXEn.perp()<CellRmin) crossingEn = false;
+ if(posXEn.perp()>CellRmax) crossingEn = false;
+ if(debugLevel>=1&&!crossingEn) std::cout << " check crossing in R " << crossingEn << " CellRmin " << CellRmin << " CellRmax " << CellRmax << " position crossing R " << posXEn.perp() << std::endl;
+ }
+ if(typeCheck==4) {
+ if(isign*posXEn.z()<isign*CellZmin) crossingEn = false;
+ if(isign*posXEn.z()>isign*CellZmax) crossingEn = false;
+ if(debugLevel>=1&&!crossingEn) std::cout << " check crossing in Z " << crossingEn << " CellZmin " << CellZmin << " CellZmax " << CellZmax << " position crossing Z " << posXEn.z() << std::endl;
+ }
+ }
+
+
+
+ if(debugLevel>=2) {
+ std::cout << " Entry pathE solution " << pathEn.x() << " y " << pathEn.y() << " z " << pathEn.z() << " crossing " << crossingEn << std::endl;
+ std::cout << " crossing entrance plane with bounds " << posXEn.x() << " y " << posXEn.y() << " z " << posXEn.z() << std::endl;
+ Amg::Vector3D check = entrance - dirT*pathEn.z();
+ std::cout << " crossing entrance plane check " << check.x() << " y " << check.y() << " z " << check.z() << std::endl;
+ }
+
+ if(!crossingEn) return 0.;
+
+
+// cross with exit plane
+
+ Amg::Vector3D pathEx;
+ if(debugLevel>=2) std::cout << " Try crossing Exit first try for type " << type << std::endl;
+ bool crossingEx = crossingMatrix(mExit,posEx,pathEx);
+ if(debugLevel>=2&&crossingEx) std::cout << " Succesfull crossing Exit first try for type " << type << " type0 " << std::endl;
+ typeCheck = type0;
+ if(!crossingEx) {
+ int type1 = int((type-100*int(type/100))/10);
+ int type2 = int(type/100);
+ if(type1==1) posEx += Corner7/2.;
+ if(type1==2) posEx += Corner4/2.;
+ if(type1==4) posEx += Corner2/2.;
+ if(type0==1) {
+ mExit.col(1) = dir67;
+ posEx -= Corner7/2.;
+ }
+ if(type0==2) {
+ mExit.col(1) = dir64;
+ posEx -= Corner4/2.;
+ }
+ if(type0==4) {
+ mExit.col(1) = dir62;
+ posEx -= Corner2/2.;
+ }
+ if(debugLevel>=2) std::cout << " Try crossing Exit second try for type " << type << std::endl;
+ crossingEx = crossingMatrix(mExit,posEx,pathEx);
+ typeCheck = type1;
+ if(debugLevel>=2&&crossingEx) std::cout << " Succesfull crossing Exit second try for type " << type << std::endl;
+ if(!crossingEx) {
+ if(type2==1) posEx += Corner7/2.;
+ if(type2==2) posEx += Corner4/2.;
+ if(type2==4) posEx += Corner2/2.;
+ if(type1==1) {
+ mExit.col(0) = dir67;
+ posEx -= Corner7/2.;
+ }
+ if(type1==2) {
+ mExit.col(0) = dir64;
+ posEx -= Corner4/2.;
+ }
+ if(type1==4) {
+ mExit.col(0) = dir62;
+ posEx -= Corner2/2.;
+ }
+ if(debugLevel>=2) std::cout << " Try crossing Exit third try for type " << type << std::endl;
+ crossingEx = crossingMatrix(mExit,posEx,pathEx);
+ typeCheck = type2;
+ if(debugLevel>=2&&crossingEx) std::cout << " Succesfull crossing Exit third try for type " << type << std::endl;
+ }
+ }
+
+ Amg::Vector3D posXEx = posXEn;
+ if(crossingEx) {
+ posXEx = entrance - dirT*pathEx.z();
+ if(typeCheck==1) {
+ double dphiCheck = acos(cos(posXEx.phi())*cos((CellPhimax+CellPhimin)/2.)+sin(posXEx.phi())*sin((CellPhimax+CellPhimin)/2.));
+ if(dphiCheck>fabs(CellPhimax-CellPhimin)) {
+ crossingEx = false;
+ }
+ if(debugLevel>=1&&!crossingEx) std::cout << " check crossing in phi " << crossingEx << " dphiCheck " << dphiCheck << " CellPhimin " << CellPhimin << " CellPhimax " << CellPhimax << " position crossing phi " << posXEx.phi() << std::endl;
+ }
+ if(typeCheck==2) {
+ if(posXEx.perp()<CellRmin) crossingEx = false;
+ if(posXEx.perp()>CellRmax) crossingEx = false;
+ if(debugLevel>=1&&!crossingEx) std::cout << " check crossing in R " << crossingEx << " CellRmin " << CellRmin << " CellRmax " << CellRmax << " position crossing R " << posXEx.perp() << std::endl;
+ }
+ if(typeCheck==4) {
+ if(isign*posXEx.z()<isign*CellZmin) crossingEx = false;
+ if(isign*posXEx.z()>isign*CellZmax) crossingEx = false;
+ if(debugLevel>=1&&!crossingEx) std::cout << " check crossing in Z " << crossingEx << " CellZmin " << CellZmin << " CellZmax " << CellZmax << " position crossing Z " << posXEx.z() << std::endl;
+ }
+ }
+
+ if(debugLevel>=2) {
+ std::cout << " Exit pathEx solution " << pathEx.x() << " y " << pathEx.y() << " z " << pathEx.z() << " crossing " << crossingEx << std::endl;
+ std::cout << " crossing exit plane with bounds " << posXEx.x() << " y " << posXEx.y() << " z " << posXEx.z() << std::endl;
+ Amg::Vector3D check = entrance - dirT*pathEx.z();
+ std::cout << " crossing exit plane check " << check.x() << " y " << check.y() << " z " << check.z() << std::endl;
+ }
+
+ double pathT = (entrance-exit).mag();
+ double path = 0.;
+ if(crossingEn&&crossingEx) {
+ path = (posXEx-posXEn).mag();
+ if(debugLevel>=1) {
+ std::cout << std::endl;
+ std::cout << " pathT " << pathT << " path " << path << std::endl;
+ std::cout << " entrance position " << entrance.x() << " y " << entrance.y() << " z " << entrance.z() << std::endl;
+ std::cout << " crossing entrance plane with bounds " << posXEn.x() << " y " << posXEn.y() << " z " << posXEn.z() << std::endl;
+ std::cout << " exit position " << exit.x() << " y " << exit.y() << " z " << exit.z() << std::endl;
+ std::cout << " crossing exit plane with bounds " << posXEx.x() << " y " << posXEx.y() << " z " << posXEx.z() << std::endl;
+ }
+ }
+ if(debugLevel>=2) {
+ std::cout << " pathT " << pathT << " path " << path << std::endl;
+ }
+
+ return path;
+
+}
+bool PathLengthUtils::crossingMatrix(Amg::MatrixX Matrix ,Amg::Vector3D entry, Amg::Vector3D& path) {
+
+// std::cout << " Matrix determinant " << Matrix.determinant() << std::endl;
+// std::cout << " Matrix " << std::endl;
+// std::cout << Matrix << std::endl;
+
+ if(Matrix.determinant()==0) {
+ // std::cout << " PathLengthUtils::crossingMatrix Matrix determinant ZERO " << std::endl;
+ return false;
+ }
+ Amg::MatrixX Minv = Matrix.inverse();
+
+// std::cout << " inverse Matrix " << std::endl;
+// std::cout << Minv << std::endl;
+
+ path = Minv*entry;
+
+ bool crossing = false;
+//
+// Inside middle of the plane ? range -0.5 to 0.5
+//
+ if(fabs(path.x())<0.5&&fabs(path.y())<0.5) crossing = true;
+
+ return crossing;
+}
+
+double PathLengthUtils::getPathLengthInTile(const CaloCell& cell, const Amg::Vector3D& entrance, const Amg::Vector3D& exit){
+//=========================================================================================================================
+ // OBTAIN LAYER INDICES FOR LINEAR INTERPOLATION
+ unsigned int SampleID = cell.caloDDE()->getSampling();
+
+ double CellZB[9] = {141.495, 424.49, 707.485, 999.605, 1300.855, 1615.8, 1949., 2300.46, 2651.52};
+ double CellDZB[9] = {282.99, 283., 282.99, 301.25, 301.25, 328.64, 337.76, 365.16, 336.96};
+ double CellZC[9] = {159.755, 483.83, 812.465, 1150.23, 1497.125, 1857.71, 2241.12, 2628.695,0};
+ double CellDZC[9] = {319.51, 328.64, 328.63, 346.9, 346.89, 374.28, 392.54, 382.61,0};
+
+ // OBTAIN TRACK AND CELL PARAMETERS
+ double pathl = 0.;
+
+ double Layer1X(exit.x()),Layer1Y(exit.y()),Layer1Z(exit.z());
+ double Layer2X(entrance.x()),Layer2Y(entrance.y()),Layer2Z(entrance.z());
+
+ double CellPhimin = cell.caloDDE()->phi()-cell.caloDDE()->dphi()*0.5;
+ double CellPhimax = cell.caloDDE()->phi()+cell.caloDDE()->dphi()*0.5;
+ double CellZmin = cell.caloDDE()->z()-cell.caloDDE()->dz()*0.5;
+ double CellZmax = cell.caloDDE()->z()+cell.caloDDE()->dz()*0.5;
+ double CellRmin = cell.caloDDE()->r()-cell.caloDDE()->dr()*0.5;
+ double CellRmax = cell.caloDDE()->r()+cell.caloDDE()->dr()*0.5;
+
+ // FIX FOR CELLS WITH ZERO WIDTH IN R
+ if(cell.caloDDE()->dr() == 0){
+ // V = dr*dphi*r*dz
+ double CellVolume = cell.caloDDE()->volume();
+ double product = cell.caloDDE()->dphi() * cell.caloDDE()->r() * cell.caloDDE()->dz();
+ double drFormula = 0;
+ if(product!=0){
+ drFormula = CellVolume / product;
+ } // IF
+ CellRmin = cell.caloDDE()->r() - drFormula*0.5;
+ CellRmax = cell.caloDDE()->r() + drFormula*0.5;
+ } // IF
+
+ //TileCellDim *cell_dim = m_tileDDM->get_cell_dim(cell->caloDDE()->identify());
+
+ double CellXimp[2], CellYimp[2], CellZimp[2];
+ double X(0), Y(0), Z(0), R(0), Phi(0);
+ double DeltaPhi;
+
+ // COMPUTE PATH
+ bool compute = true;
+ int lBC(0);
+ // LOOP IS USUALLY RUN ONCE, EXCEPT FOR LADDER SHAPED TILECAL CELLS
+ while(compute){
+ if(lBC==1) break;
+ int Np = 0;
+ if(sqrt((Layer1X-Layer2X)*(Layer1X-Layer2X)+(Layer1Y-Layer2Y)*(Layer1Y-Layer2Y)) < 3818.5){
+ if(SampleID == 13 && lBC == 0){
+ int cpos = fabs(cell.caloDDE()->eta())/0.1;
+
+ CellRmin = 2600.; //cell_dim->getRMin(0);
+ CellRmax = 2990.; //cell_dim->getRMax(2);
+ CellZmin = CellZB[cpos] - 0.5*CellDZB[cpos]; //cell_dim->getZMin(0);
+ CellZmax = CellZB[cpos] + 0.5*CellDZB[cpos]; //cell_dim->getZMax(0);
+ } // ELSE IF
+ else if(SampleID == 13 && lBC == 1){
+ int cpos = fabs(cell.caloDDE()->eta())/0.1;
+
+ CellRmin = 2990.; //cell_dim->getRMin(3);
+ CellRmax = 3440.; //cell_dim->getRMax(5);
+ CellZmin = CellZC[cpos] - 0.5*CellDZC[cpos]; //cell_dim->getZMin(3);
+ CellZmax = CellZC[cpos] + 0.5*CellDZC[cpos]; //cell_dim->getZMax(3);
+ } // ELSE IF
+
+ // CALCULATE GRADIENTS
+ double Sxy = (Layer2X-Layer1X)/(Layer2Y-Layer1Y);
+ double Sxz = (Layer2X-Layer1X)/(Layer2Z-Layer1Z);
+ double Syz = (Layer2Y-Layer1Y)/(Layer2Z-Layer1Z);
+
+ // CALCULATE POINTS OF INTERSECTION
+ // INTERSECTIONS R PLANES
+ double Radius(CellRmin);
+
+ double x0int(exit.x()), x1int(entrance.x()),
+ y0int(exit.y()), y1int(entrance.y()),
+ z0int(exit.z()), z1int(entrance.z());
+ double S((y1int-y0int)/(x1int-x0int));
+ double a(1+S*S), b(2*S*y0int-2*S*S*x0int), c(y0int*y0int-Radius*Radius+S*S*x0int*x0int-2*y0int*S*x0int);
+ double x1((-b+sqrt(b*b-4*a*c))/(2*a)), x2((-b-sqrt(b*b-4*a*c))/(2*a));
+ double y1(y0int+S*(x1-x0int)), y2(y0int+S*(x2-x0int));
+ double S1 = ((z1int-z0int)/(x1int-x0int));
+ double z1(z0int+S1*(x1-x0int)), z2(z0int+S1*(x2-x0int));
+
+ X = x1;
+ Y = y1;
+ Z = z1;
+
+ if( ((x1-x0int)*(x1-x0int)+(y1-y0int)*(y1-y0int)+(z1-z0int)*(z1-z0int)) >
+ ((x2-x0int)*(x2-x0int)+(y2-y0int)*(y2-y0int)+(z2-z0int)*(z2-z0int)) ){
+ X=x2;
+ Y=y2;
+ Z=z2;
+ } // IF
+
+ if(CellRmin != CellRmax){
+ Phi = acos(X/sqrt(X*X+Y*Y));
+ if(Y <= 0) Phi = -Phi;
+ R = CellRmin;
+
+ if(Z>=CellZmin && Z<=CellZmax && Phi>=CellPhimin && Phi<=CellPhimax){
+ CellXimp[Np]=X;
+ CellYimp[Np]=Y;
+ CellZimp[Np]=Z;
+ Np=Np+1;
+
+ } // IF
+
+ Radius = CellRmax;
+
+ c = y0int*y0int-Radius*Radius+S*S*x0int*x0int-2*y0int*S*x0int;
+ x1 = ((-b+sqrt(b*b-4*a*c))/(2*a));
+ x2 = ((-b-sqrt(b*b-4*a*c))/(2*a));
+ y1 = (y0int+S*(x1-x0int));
+ y2 = (y0int+S*(x2-x0int));
+ z1 = (z0int+S1*(x1-x0int));
+ z2 = (z0int+S1*(x2-x0int));
+ S1 = ((z1int-z0int)/(x1int-x0int));
+
+ X = x1;
+ Y = y1;
+ Z = z1;
+
+ if( ((x1-x0int)*(x1-x0int)+(y1-y0int)*(y1-y0int)+(z1-z0int)*(z1-z0int)) >
+ ((x2-x0int)*(x2-x0int)+(y2-y0int)*(y2-y0int)+(z2-z0int)*(z2-z0int)) ){
+ X = x2;
+ Y = y2;
+ Z = z2;
+ } // IF
+
+ Phi=std::acos(X/sqrt(X*X+Y*Y));
+ if (Y <= 0) Phi=-Phi;
+ R=CellRmax;
+
+ if(Z>=CellZmin && Z<=CellZmax && Phi>=CellPhimin && Phi<=CellPhimax){
+ CellXimp[Np]=X; CellYimp[Np]=Y; CellZimp[Np]=Z;
+ Np=Np+1;
+ } // IF
+ } // IF
+
+ // INTERSECTIONS Z PLANES
+ if(CellZmin != CellZmax){
+ if(Np < 2){
+ Z = CellZmin;
+ X = Layer1X+Sxz*(Z-Layer1Z);
+ Y = Layer1Y+Syz*(Z-Layer1Z);
+ R = sqrt(X*X+Y*Y);
+ Phi = std::acos(X/R);
+ if(Y <= 0) Phi=-Phi;
+ if(R>=CellRmin && R<=CellRmax && Phi>=CellPhimin && Phi<=CellPhimax){
+ CellXimp[Np]=X; CellYimp[Np]=Y; CellZimp[Np]=Z;
+ Np=Np+1;
+ } // IF
+ } // IF
+
+ if(Np < 2){
+ Z = CellZmax;
+ X = Layer1X+Sxz*(Z-Layer1Z);
+ Y = Layer1Y+Syz*(Z-Layer1Z);
+ R = sqrt(X*X+Y*Y);
+ Phi = std::acos(X/R);
+ if(Y <= 0) Phi=-Phi;
+ if(R>=CellRmin && R<=CellRmax && Phi>=CellPhimin && Phi<=CellPhimax){
+ CellXimp[Np]=X; CellYimp[Np]=Y; CellZimp[Np]=Z;
+ Np=Np+1;
+ } // IF
+ } // IF
+ } // IF
+
+ // INTERSECTIONS PHI PLANES
+ if(CellPhimin != CellPhimax){
+ if(Np < 2){
+ X = (Layer1X-Sxy*Layer1Y)/(1-Sxy*tan(CellPhimin));
+ Y = X*std::tan(CellPhimin);
+ Z = Layer1Z+(1/Sxz)*(X-Layer1X);
+ R = sqrt(X*X+Y*Y);
+ Phi = std::acos(X/R);
+ if(Y <= 0) Phi=-Phi;
+ DeltaPhi=fabs(Phi-CellPhimin);
+ if(DeltaPhi > 3.141593) DeltaPhi=fabs(Phi+CellPhimin);
+ if(R>=CellRmin && R<=CellRmax && Z>=CellZmin && Z<=CellZmax && DeltaPhi<0.0001){
+ CellXimp[Np]=X; CellYimp[Np]=Y; CellZimp[Np]=Z;
+ Np=Np+1;
+ } // IF
+ } // IF
+
+ if(Np < 2){
+ X = (Layer1X-Sxy*Layer1Y)/(1-Sxy*tan(CellPhimax));
+ Y = X*std::tan(CellPhimax);
+ Z = Layer1Z+(1/Sxz)*(X-Layer1X);
+ R = sqrt(X*X+Y*Y);
+ Phi=acos(X/R);
+ if(Y <= 0) Phi=-Phi;
+ DeltaPhi=fabs(Phi-CellPhimax);
+ if(DeltaPhi > 3.141593) DeltaPhi=fabs(Phi+CellPhimax);
+ if(R>=CellRmin && R<=CellRmax && Z>=CellZmin && Z<=CellZmax && DeltaPhi<0.0001){
+ CellXimp[Np]=X; CellYimp[Np]=Y; CellZimp[Np]=Z;
+ Np=Np+1;
+ } // IF
+ } // IF
+ } // IF
+
+ // CALCULATE PATH IF TWO INTERSECTIONS WERE FOUND
+ if(Np == 2){
+ pathl += sqrt( (CellXimp[0]-CellXimp[1])*(CellXimp[0]-CellXimp[1]) +
+ (CellYimp[0]-CellYimp[1])*(CellYimp[0]-CellYimp[1]) +
+ (CellZimp[0]-CellZimp[1])*(CellZimp[0]-CellZimp[1]) );
+ } // IF
+ } // IF
+ if(SampleID == 13 && lBC == 0) ++lBC;
+ else compute = false;
+ } // WHILE (FOR LBBC LAYER)
+
+ return pathl;
+} // PathLengthUtils::getPathLengthInTile
+
+
+
+CaloSampling::CaloSample PathLengthUtils::tileEntrance(CaloSampling::CaloSample sample){
+ if(sample==CaloSampling::TileBar0||sample==CaloSampling::TileBar1||sample==CaloSampling::TileBar2||sample==CaloSampling::TileGap2) return CaloSampling::TileBar0;
+ if(sample==CaloSampling::TileGap1) return CaloSampling::TileBar1;
+ if(sample==CaloSampling::TileGap3||sample==CaloSampling::TileExt0||sample==CaloSampling::TileExt1||sample==CaloSampling::TileExt2) return CaloSampling::TileExt0;
+ return CaloSampling::Unknown;
+ //return sample;
+} // PathLengthUtils::entrance
+
+CaloSampling::CaloSample PathLengthUtils::tileExit(CaloSampling::CaloSample sample){
+ if(sample==CaloSampling::TileBar0||sample==CaloSampling::TileBar1||sample==CaloSampling::TileBar2||sample==CaloSampling::TileGap1) return CaloSampling::TileBar2;
+ if(sample==CaloSampling::TileGap2) return CaloSampling::TileBar1;
+ if(sample==CaloSampling::TileGap3) return CaloSampling::TileExt1;
+ if(sample==CaloSampling::TileExt0||sample==CaloSampling::TileExt1||sample==CaloSampling::TileExt2) return CaloSampling::TileExt2;
+ return CaloSampling::Unknown;
+ //return sample;
+} // PathLengthUtils::exit
+
+/** Return the length(mm) of the path crossed inside the cell, given the parameters for the extrapolation at entrance and exit of the layer **/
+double PathLengthUtils::pathInsideCell(const CaloCell& cell, const CaloExtensionHelpers::EntryExitLayerMap& entryExitLayerMap) {
+
+ CaloSampling::CaloSample sample = cell.caloDDE()->getSampling();
+
+ //------------------------------------------------------------------------------------------
+ //special treatment for tile calo cells
+ CaloSampling::CaloSample tileEntranceID = tileEntrance(sample);
+ CaloSampling::CaloSample tileExitID = tileExit(sample);
+
+ auto tileEntrancePair = entryExitLayerMap.find(tileEntranceID);
+ auto tileExitPair = entryExitLayerMap.find(tileExitID);
+ if( tileEntrancePair != entryExitLayerMap.end() && tileExitPair != entryExitLayerMap.end() ){
+ return getPathLengthInTile( cell, tileEntrancePair->second.first, tileExitPair->second.second);
+ }
+ //------------------------------------------------------------------------------------------
+
+ auto entryExitPair = entryExitLayerMap.find(sample);
+ if( entryExitPair == entryExitLayerMap.end() ) return -1.;
+
+ const Amg::Vector3D& entry = entryExitPair->second.first;
+ const Amg::Vector3D& exit = entryExitPair->second.second;
+
+ if (!crossedPhi(cell, entry.phi(), exit.phi())) return -1.;
+ double pathCrossed = -1;
+ if (cell.caloDDE()->getSubCalo() != CaloCell_ID::TILE){
+ pathCrossed = getPathLengthInEta(cell, entry.eta(), exit.eta());
+ }
+ else{
+ pathCrossed = getPathLengthInZ(cell, entry.z(), exit.z());
+ }
+ if (pathCrossed <= 0) return -1.;
+ return pathCrossed * (exit-entry).mag();
+}
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/PathLengthUtils.h b/Reconstruction/RecoTools/TrackToCalo/src/PathLengthUtils.h
new file mode 100755
index 0000000000000000000000000000000000000000..17de47b05a3e1be4715ce8103fa6888791caa63b
--- /dev/null
+++ b/Reconstruction/RecoTools/TrackToCalo/src/PathLengthUtils.h
@@ -0,0 +1,101 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/***************************************************************************
+ PathLengthUtils.h - Description
+ -------------------
+ begin : Summer 2014
+ authors : Marius Cornelis van Woerden
+ email : mvanwoer@cern.ch
+ changes :
+
+ ***************************************************************************/
+
+#ifndef PATHLENGTUTILS_H
+#define PATHLENGTUTILS_H
+
+#include "xAODTracking/TrackParticle.h"
+#include "xAODTracking/TrackParticleContainer.h"
+
+#include "CaloEvent/CaloCellContainer.h"
+#include "CaloEvent/CaloCell.h"
+
+#include "CaloGeoHelpers/CaloSampling.h"
+#include "CaloGeoHelpers/CaloPhiRange.h"
+
+#include <cmath>
+
+#include "TrkCaloExtension/CaloExtensionHelpers.h"
+
+ class PathLengthUtils{
+ public:
+ PathLengthUtils();
+ ~PathLengthUtils();
+
+
+ double pathInsideCell(const CaloCell& cell, const CaloExtensionHelpers::EntryExitLayerMap& entryExitLayerMap);
+ double get3DPathLength(const CaloCell& cell, const Amg::Vector3D& entry, const Amg::Vector3D& exit, double drFix, double dzFix );
+
+ double getPathLengthInTile(const CaloCell& cell, const Amg::Vector3D& entry, const Amg::Vector3D& exit );
+ //double path(xAOD::TrackParticle& trackParticle, const CaloCell* cell);
+ //double path(xAOD::TrackParticle& trackParticle, const CaloCell_ID::CaloSample sample);
+
+ private:
+ CaloSampling::CaloSample tileEntrance(CaloSampling::CaloSample sample);
+ CaloSampling::CaloSample tileExit(CaloSampling::CaloSample sample);
+
+ double phiMean(double a, double b);
+ bool crossedPhi(const CaloCell& cell, double phi_entrance, double phi_exit);
+ double getPathLengthInEta(const CaloCell& cell, double eta_entrance, double eta_exit);
+ double getPathLengthInZ(double zMin, double zMax, double z_entrance, double z_exit);
+ double getPathLengthInZ(const CaloCell& cell, double z_entrance, double z_exit);
+ bool crossingMatrix(Amg::MatrixX Matrix ,Amg::Vector3D entry, Amg::Vector3D& path);
+ //static double CellZB[9];
+ //static double CellDZB[9];
+ //static double CellZC[9];
+ //static double CellDZC[9];
+ };
+
+inline double PathLengthUtils::phiMean(double a, double b) { return 0.5*(a+b) + (a*b < 0)*M_PI; }
+
+/** Return true if the cell crossed was crossed by the track in phi **/
+inline bool PathLengthUtils::crossedPhi(const CaloCell& cell, double phi_entrance, double phi_exit) {
+ double mean_phi = phiMean(phi_entrance, phi_exit);
+ double dphi = fabs( CaloPhiRange::diff(phi_entrance, phi_exit) );
+ double phi_min = mean_phi - dphi, phi_max = mean_phi + dphi;
+
+ return (CaloPhiRange::diff(cell.phi() + cell.caloDDE()->dphi()/2., phi_min) > 0 &&
+ CaloPhiRange::diff(phi_max, cell.phi() - cell.caloDDE()->dphi()/2.) > 0);
+}
+
+/** Return the % of path length crossed by the track inside a cell in eta **/
+inline double PathLengthUtils::getPathLengthInEta(const CaloCell& cell, double eta_entrance, double eta_exit) {
+ double etaMin = cell.eta() - 0.5*cell.caloDDE()->deta();
+ double etaMax = cell.eta() + 0.5*cell.caloDDE()->deta();
+ if ( fabs(eta_entrance - eta_exit) < 1e-6 ) // to avoid FPE
+ return eta_entrance > etaMin && eta_entrance < etaMax;
+
+ double etaMinTrack = std::min(eta_entrance, eta_exit);
+ double etaMaxTrack = std::max(eta_entrance, eta_exit);
+ return (std::min(etaMax, etaMaxTrack) - std::max(etaMin, etaMinTrack))/(etaMaxTrack - etaMinTrack);
+}
+
+/** Return the % of path length crossed by the track inside a cell in Z **/
+inline double PathLengthUtils::getPathLengthInZ(double zMin, double zMax, double z_entrance, double z_exit) {
+ if ( fabs(z_entrance - z_exit) < 1e-6 ) // to avoid FPE
+ return z_entrance > zMin && z_entrance < zMax;
+
+ double zMinTrack = std::min(z_entrance, z_exit);
+ double zMaxTrack = std::max(z_entrance, z_exit);
+ return (std::min(zMax, zMaxTrack) - std::max(zMin, zMinTrack))/(zMaxTrack - zMinTrack);
+}
+
+/** Return the % of path length crossed by the track inside a cell in Z **/
+inline double PathLengthUtils::getPathLengthInZ(const CaloCell& cell, double z_entrance, double z_exit) {
+ return getPathLengthInZ(cell.z() - 0.5*cell.caloDDE()->dz(), cell.z() + 0.5*cell.caloDDE()->dz(), z_entrance, z_exit);
+}
+
+
+#endif //> !PATHLENGTUTILS_H
+
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/TrackParticleCaloExtensionAlg.cxx b/Reconstruction/RecoTools/TrackToCalo/src/TrackParticleCaloExtensionAlg.cxx
index edf896bafc74d0e076bc286ee4a3514722400c3f..95a02e6102737cc6c2c81d87dea97619effa1078 100644
--- a/Reconstruction/RecoTools/TrackToCalo/src/TrackParticleCaloExtensionAlg.cxx
+++ b/Reconstruction/RecoTools/TrackToCalo/src/TrackParticleCaloExtensionAlg.cxx
@@ -61,28 +61,28 @@ StatusCode TrackParticleCaloExtensionAlg::execute() {
if( !m_trackSelector.empty() && !m_trackSelector->decision(*tp) ) continue;
++m_processedTracks;
- // call calo Extension
- if( !m_caloExtensionTool.empty() ){
- const xAOD::ParticleCaloExtension* extension = m_caloExtensionTool->caloExtension(const_cast<xAOD::TrackParticle&>(*tp));
- if( extension ) ++m_extensions;
+ // // call calo Extension
+ // if( !m_caloExtensionTool.empty() ){
+ // const Trk::CaloExtension* extension = m_caloExtensionTool->caloExtension(const_cast<xAOD::TrackParticle&>(*tp));
+ // if( extension ) ++m_extensions;
- if( msgLvl(MSG::DEBUG) && extension && extension->numberOfParameters() != 0 ){
- ATH_MSG_DEBUG("Got calo extension " << extension->numberOfParameters() << " cells " << extension->caloCells().size() );
- for( unsigned int i=0;i<extension->numberOfParameters();++i){
- xAOD::CurvilinearParameters_t pars = extension->trackParameters(i);
- unsigned int id = extension->parameterIdentifier(i);
- ATH_MSG_DEBUG( " id " << id << " pos: r " << sqrt(pars[0]*pars[0]+pars[1]*pars[1]) << " z " << pars[2]
- << " momentum " << sqrt(pars[3]*pars[3]+pars[4]*pars[4]+pars[5]*pars[5]) );
- }
- }
- }
- // call calo cell association
- if( !m_caloCellAssociationTool.empty() ){
- const xAOD::ParticleCaloExtension* extension = m_caloCellAssociationTool->caloAssociation(const_cast<xAOD::TrackParticle&>(*tp));
- if( msgLvl(MSG::DEBUG) && extension && extension->numberOfParameters() != 0 ){
- ATH_MSG_DEBUG("Got calo cell extension " << extension->numberOfParameters() << " cells " << extension->caloCells().size() );
- }
- }
+ // if( msgLvl(MSG::DEBUG) && extension && extension->numberOfParameters() != 0 ){
+ // ATH_MSG_DEBUG("Got calo extension " << extension->numberOfParameters() << " cells " << extension->caloCells().size() );
+ // for( unsigned int i=0;i<extension->numberOfParameters();++i){
+ // xAOD::CurvilinearParameters_t pars = extension->trackParameters(i);
+ // unsigned int id = extension->parameterIdentifier(i);
+ // ATH_MSG_DEBUG( " id " << id << " pos: r " << sqrt(pars[0]*pars[0]+pars[1]*pars[1]) << " z " << pars[2]
+ // << " momentum " << sqrt(pars[3]*pars[3]+pars[4]*pars[4]+pars[5]*pars[5]) );
+ // }
+ // }
+ // }
+ // // call calo cell association
+ // if( !m_caloCellAssociationTool.empty() ){
+ // const xAOD::ParticleCaloExtension* extension = m_caloCellAssociationTool->caloAssociation(const_cast<xAOD::TrackParticle&>(*tp));
+ // if( msgLvl(MSG::DEBUG) && extension && extension->numberOfParameters() != 0 ){
+ // ATH_MSG_DEBUG("Got calo cell extension " << extension->numberOfParameters() << " cells " << extension->caloCells().size() );
+ // }
+ // }
}
return StatusCode::SUCCESS;
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/TrackParticleCaloExtensionAlg.h b/Reconstruction/RecoTools/TrackToCalo/src/TrackParticleCaloExtensionAlg.h
index 270206409535466d5f483d08862daa2c2c8cf50b..07be8ae1ea356171d3da2e20773b611289cb19cb 100644
--- a/Reconstruction/RecoTools/TrackToCalo/src/TrackParticleCaloExtensionAlg.h
+++ b/Reconstruction/RecoTools/TrackToCalo/src/TrackParticleCaloExtensionAlg.h
@@ -20,6 +20,8 @@
namespace Trk {
class ITrackSelectorTool;
class IParticleCaloExtensionTool;
+}
+namespace Rec {
class IParticleCaloCellAssociationTool;
}
@@ -37,7 +39,7 @@ public:
private:
ToolHandle <Trk::IParticleCaloExtensionTool> m_caloExtensionTool; //!< Tool to make the step-wise extrapolation
- ToolHandle <Trk::IParticleCaloCellAssociationTool> m_caloCellAssociationTool; //!< Tool to make the step-wise extrapolation
+ ToolHandle <Rec::IParticleCaloCellAssociationTool> m_caloCellAssociationTool; //!< Tool to make the step-wise extrapolation
ToolHandle <Trk::ITrackSelectorTool> m_trackSelector; //!< Tool to select tracks
std::string m_trackParicleContainerName;
diff --git a/Reconstruction/RecoTools/TrackToCalo/src/components/TrackToCalo_entries.cxx b/Reconstruction/RecoTools/TrackToCalo/src/components/TrackToCalo_entries.cxx
index 521442589ce1c36f7d45c4962f28768fbcbdec13..602b2a965dfb7d88d790a5b6a15ec1d023bd0229 100755
--- a/Reconstruction/RecoTools/TrackToCalo/src/components/TrackToCalo_entries.cxx
+++ b/Reconstruction/RecoTools/TrackToCalo/src/components/TrackToCalo_entries.cxx
@@ -1,24 +1,25 @@
//====================================================================
#include "GaudiKernel/DeclareFactoryEntries.h"
-#include "TrackToCalo/ExtrapolateToCaloTool.h"
#include "../ParticleCaloExtensionTool.h"
#include "../ParticleCaloCellAssociationTool.h"
+#include "../ParticleCaloClusterAssociationTool.h"
#include "../TrackParticleCaloExtensionAlg.h"
-#include "TrackToCalo/ImpactInCaloAlg.h"
+#include "../MuonCaloEnergyTool.h"
using namespace Trk;
-DECLARE_TOOL_FACTORY( ExtrapolateToCaloTool )
+using namespace Rec;
DECLARE_TOOL_FACTORY( ParticleCaloExtensionTool )
DECLARE_TOOL_FACTORY( ParticleCaloCellAssociationTool )
-DECLARE_ALGORITHM_FACTORY( ImpactInCaloAlg )
+DECLARE_TOOL_FACTORY( ParticleCaloClusterAssociationTool )
+DECLARE_TOOL_FACTORY( MuonCaloEnergyTool )
DECLARE_ALGORITHM_FACTORY( TrackParticleCaloExtensionAlg )
DECLARE_FACTORY_ENTRIES ( TrackToCalo ) {
- DECLARE_TOOL( ExtrapolateToCaloTool );
DECLARE_TOOL( ParticleCaloExtensionTool );
DECLARE_TOOL( ParticleCaloCellAssociationTool );
- DECLARE_ALGORITHM( ImpactInCaloAlg );
+ DECLARE_TOOL( ParticleCaloClusterAssociationTool );
+ DECLARE_TOOL( MuonCaloEnergyTool );
DECLARE_ALGORITHM( TrackParticleCaloExtensionAlg );
}