Reference Manual

Go to the documentation of this file.
00001 #ifndef RecoEcal_EgammaCoreTools_EcalClusterTools_h
00002 #define RecoEcal_EgammaCoreTools_EcalClusterTools_h
00003 
00016 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00017 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00018 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00019 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
00020 //#include "DataFormats/Math/interface/Point3D.h"
00021 #include "DataFormats/Math/interface/Vector3D.h"
00022 //includes for ShowerShape function to work
00023 #include <vector>
00024 #include <math.h>
00025 #include <TMath.h>
00026 #include <TMatrixT.h>
00027 #include <TMatrixD.h>
00028 #include <TVectorT.h>
00029 #include <TVectorD.h>
00030 
00031 
00032 class DetId;
00033 class CaloTopology;
00034 class CaloGeometry;
00035 
00036 struct Cluster2ndMoments {
00037 
00038   // major and minor cluster moments wrt principale axes:
00039   float sMaj;
00040   float sMin;
00041   // angle between sMaj and phi:
00042   float alpha;
00043 
00044 };
00045 
00046 class EcalClusterTools {
00047         public:
00048                 EcalClusterTools() {};
00049                 ~EcalClusterTools() {};
00050 
00051                 // various energies in the matrix nxn surrounding the maximum energy crystal of the input cluster
00052                 //we use an eta/phi coordinate system rather than phi/eta
00053                 //note e3x2 does not have a definate eta/phi geometry, it takes the maximum 3x2 block containing the 
00054                 //seed regardless of whether that 3 in eta or phi
00055                 static float e1x3( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00056 
00057 static float e1x3( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
00058 
00059                 static float e3x1( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00060 
00061 static float e3x1( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
00062 
00063                 static float e1x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00064                 //AA
00065                 //Take into account severities
00066                 static float e1x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
00067 
00068                 static float e5x1( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00069 static float e5x1( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
00070 
00071                 static float e2x2( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00072 static float e2x2( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
00073 
00074                 static float e3x2( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00075 static float e3x2( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
00076 
00077                 static float e3x3( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00078 static float e3x3( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
00079 
00080                 static float e4x4( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology);
00081 static float e4x4( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
00082 
00083                 static float e5x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00084 static float e5x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
00085 
00086                 // energy in the 2x5 strip right of the max crystal (does not contain max crystal)
00087                 // 2 crystals wide in eta, 5 wide in phi.
00088                 static float e2x5Right( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00089 static float e2x5Right( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
00090                 // energy in the 2x5 strip left of the max crystal (does not contain max crystal)
00091 
00092                 static float e2x5Left( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00093 static float e2x5Left( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
00094                 // energy in the 5x2 strip above the max crystal (does not contain max crystal)
00095                 // 5 crystals wide in eta, 2 wide in phi.
00096 
00097                 static float e2x5Top( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00098 static float e2x5Top( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
00099                 // energy in the 5x2 strip below the max crystal (does not contain max crystal)                
00100 
00101                 static float e2x5Bottom( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00102 static float e2x5Bottom( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
00103                 // energy in a 2x5 strip containing the seed (max) crystal.
00104                 // 2 crystals wide in eta, 5 wide in phi.
00105                 // it is the maximum of either (1x5left + 1x5center) or (1x5right + 1x5center)
00106                 static float e2x5Max( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00107 static float e2x5Max( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
00108 
00109                 // energies in the crystal left, right, top, bottom w.r.t. to the most energetic crystal
00110                 static float eLeft( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00111 static float eLeft( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
00112 
00113                 static float eRight( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00114 static float eRight( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
00115 
00116                 static float eTop( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00117 static float eTop( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
00118 
00119                 static float eBottom( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
00120 static float eBottom( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
00121                 // the energy of the most energetic crystal in the cluster
00122 
00123                 static float eMax( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits );
00124 static float eMax( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
00125 
00126                 // the energy of the second most energetic crystal in the cluster
00127                 static float e2nd( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits );
00128 static float e2nd( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
00129 
00130                 // get the DetId and the energy of the maximum energy crystal of the input cluster
00131                 static std::pair<DetId, float> getMaximum( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits);
00132 static std::pair<DetId, float> getMaximum( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
00133 
00134                 static std::vector<float> energyBasketFractionEta( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits );
00135 static std::vector<float> energyBasketFractionEta( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
00136 
00137                 static std::vector<float> energyBasketFractionPhi( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits);
00138 static std::vector<float> energyBasketFractionPhi( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
00139 
00140                 // return a vector v with v[0] = etaLat, v[1] = phiLat, v[2] = lat
00141                 static std::vector<float> lat( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, bool logW = true, float w0 = 4.7 );
00142 static std::vector<float> lat( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv, bool logW = true,float w0 = 4.7);
00143 
00144                 // return a vector v with v[0] = covEtaEta, v[1] = covEtaPhi, v[2] = covPhiPhi
00145 
00146                 static std::vector<float> covariances(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits, const CaloTopology *topology, const CaloGeometry* geometry, float w0 = 4.7);
00147 static std::vector<float> covariances(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits, const CaloTopology *topology, const CaloGeometry* geometry,const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv,float w0 = 4.7);
00148                 // return a vector v with v[0] = covIEtaIEta, v[1] = covIEtaIPhi, v[2] = covIPhiIPhi
00149                 //this function calculates differences in eta/phi in units of crystals not global eta/phi
00150                 //this is gives better performance in the crack regions of the calorimeter but gives otherwise identical results to covariances function
00151                 //   except that it doesnt need an eta based correction funtion in the endcap 
00152                 //it is multipled by an approprate crystal size to ensure it gives similar values to covariances(...)
00153                 //
00154                 //Warning: covIEtaIEta has been studied by egamma, but so far covIPhiIPhi hasnt been studied extensively so there could be a bug in 
00155                 //         the covIPhiIEta or covIPhiIPhi calculations. I dont think there is but as it hasnt been heavily used, there might be one
00156                 static std::vector<float> localCovariances(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits, const CaloTopology *topology, float w0 = 4.7);
00157 static std::vector<float> localCovariances(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits, const CaloTopology *topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv,float w0 = 4.7);
00158                 
00159                 static std::vector<float> scLocalCovariances(const reco::SuperCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology, float w0 = 4.7);
00160 static std::vector<float> scLocalCovariances(const reco::SuperCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv, float w0 = 4.7);
00161 
00162                 // cluster second moments with respect to principal axes:
00163                 static Cluster2ndMoments cluster2ndMoments( const reco::BasicCluster &basicCluster, const EcalRecHitCollection &recHits, double phiCorrectionFactor=0.8, double w0=4.7, bool useLogWeights=true);
00164 
00165                 static Cluster2ndMoments cluster2ndMoments( const reco::SuperCluster &superCluster, const EcalRecHitCollection &recHits, double phiCorrectionFactor=0.8, double w0=4.7, bool useLogWeights=true);
00166                 static Cluster2ndMoments cluster2ndMoments( const std::vector<const EcalRecHit*>& RH_ptrs, double  phiCorrectionFactor=0.8, double  w0=4.7, bool useLogWeights=true);
00167 
00168                 static double zernike20( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, double R0 = 6.6, bool logW = true, float w0 = 4.7 );
00169                 static double zernike42( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, double R0 = 6.6, bool logW = true, float w0 = 4.7 );
00170 
00171                 // get the detId's of a matrix centered in the maximum energy crystal = (0,0)
00172                 // the size is specified by ixMin, ixMax, iyMin, iyMax in unit of crystals
00173                 static std::vector<DetId> matrixDetId( const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax );
00174 static std::vector<DetId> matrixDetId( const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv  );
00175                 // get the energy deposited in a matrix centered in the maximum energy crystal = (0,0)
00176                 // the size is specified by ixMin, ixMax, iyMin, iyMax in unit of crystals
00177                 static float matrixEnergy( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax );
00178 
00179                 //AA
00180                 //Take into account the severities and flags
00181                 static float matrixEnergy( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
00182                 static float getFraction( const std::vector< std::pair<DetId, float> > &v_id, DetId id);
00183                 // get the DetId and the energy of the maximum energy crystal in a vector of DetId
00184                 static std::pair<DetId, float> getMaximum( const std::vector< std::pair<DetId, float> > &v_id, const EcalRecHitCollection *recHits);
00185 static std::pair<DetId, float> getMaximum( const std::vector< std::pair<DetId, float> > &v_id, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
00186                 // get the energy of a DetId, return 0 if the DetId is not in the collection
00187                 static float recHitEnergy(DetId id, const EcalRecHitCollection *recHits);
00188 
00189                 //AA
00190                 //Take into account severities and flags
00191                 static float recHitEnergy(DetId id, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
00192 
00193                 //Shower shape variables return vector <Roundness, Angle> of a photon
00194                 static std::vector<float> roundnessBarrelSuperClusters( const reco::SuperCluster &superCluster ,const EcalRecHitCollection &recHits, int weightedPositionMethod = 0, float energyThreshold = 0.0);
00195                 static std::vector<float> roundnessBarrelSuperClustersUserExtended( const reco::SuperCluster &superCluster ,const EcalRecHitCollection &recHits, int ieta_delta=0, int iphi_delta=0, float energyRHThresh=0.00000, int weightedPositionMethod=0);
00196                 static std::vector<float> roundnessSelectedBarrelRecHits(const std::vector<const EcalRecHit*>&rhVector, int weightedPositionMethod = 0);
00197         private:
00198                 struct EcalClusterEnergyDeposition
00199                 { 
00200                         double deposited_energy;
00201                         double r;
00202                         double phi;
00203                 };
00204 
00205                 static std::vector<EcalClusterEnergyDeposition> getEnergyDepTopology( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, bool logW, float w0 );
00206 
00207                 static math::XYZVector meanClusterPosition( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology, const CaloGeometry *geometry );
00208 static math::XYZVector meanClusterPosition( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology, const CaloGeometry *geometry, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
00209                 //return energy weighted mean distance from the seed crystal in number of crystals
00210                 //<iEta,iPhi>, iPhi is not defined for endcap and is returned as zero 
00211                 static std::pair<float,float>  mean5x5PositionInLocalCrysCoord(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology);
00212 static std::pair<float,float>  mean5x5PositionInLocalCrysCoord(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
00213 
00214                 static std::pair<float,float> mean5x5PositionInXY(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology);
00215 static std::pair<float,float> mean5x5PositionInXY(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
00216 
00217                 static double f00(double r) { return 1; }
00218                 static double f11(double r) { return r; }
00219                 static double f20(double r) { return 2.0*r*r-1.0; }
00220                 static double f22(double r) { return r*r; }
00221                 static double f31(double r) { return 3.0*r*r*r - 2.0*r; }
00222                 static double f33(double r) { return r*r*r; }
00223                 static double f40(double r) { return 6.0*r*r*r*r-6.0*r*r+1.0; }
00224                 static double f42(double r) { return 4.0*r*r*r*r-3.0*r*r; }
00225                 static double f44(double r) { return r*r*r*r; }
00226                 static double f51(double r) { return 10.0*pow(r,5)-12.0*pow(r,3)+3.0*r; }
00227                 static double f53(double r) { return 5.0*pow(r,5) - 4.0*pow(r,3); }
00228                 static double f55(double r) { return pow(r,5); }
00229 
00230                 static double absZernikeMoment( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, int n, int m, double R0, bool logW, float w0 );
00231                 static double fast_AbsZernikeMoment(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, int n, int m, double R0, bool logW, float w0 );
00232                 static double calc_AbsZernikeMoment(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, int n, int m, double R0, bool logW, float w0 );
00233 
00234                 static double factorial(int n) {
00235                         double res = 1.;
00236                         for (int i = 2; i <= n; ++i) res *= i;
00237                         return res;
00238                 }
00239 
00240                 //useful functions for the localCovariances function
00241                 static float getIEta(const DetId& id);
00242                 static float getIPhi(const DetId& id);
00243                 static float getNormedIX(const DetId& id);
00244                 static float getNormedIY(const DetId& id);
00245                 static float getDPhiEndcap(const DetId& crysId,float meanX,float meanY);
00246                 static float getNrCrysDiffInEta(const DetId& crysId,const DetId& orginId);
00247                 static float getNrCrysDiffInPhi(const DetId& crysId,const DetId& orginId);
00248 
00249                                 //useful functions for showerRoundnessBarrel function
00250                                 static int deltaIEta(int seed_ieta, int rh_ieta);
00251                                 static int deltaIPhi(int seed_iphi, int rh_iphi);
00252                                 static std::vector<int> getSeedPosition(const std::vector<const EcalRecHit*>&RH_ptrs);
00253                                 static float getSumEnergy(const std::vector<const EcalRecHit*>&RH_ptrs);
00254                                 static float computeWeight(float eRH, float energyTotal, int weightedPositionMethod);
00255 
00256 };
00257 
00258 #endif
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h
