EgammaTools/src/EcalRegressionData.cc

0001 #include "RecoEgamma/EgammaTools/interface/EcalRegressionData.h"
0002
0003 #include "RecoEcal/EgammaCoreTools/interface/EcalTools.h"
0004 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0005 #include "DataFormats/Math/interface/deltaR.h"
0006 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0007 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0008 #include "RecoEgamma/EgammaTools/interface/EcalClusterLocal.h"
0009
0010 #include "DataFormats/VertexReco/interface/Vertex.h"
0011 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0012
0013 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0014 #include "FWCore/Utilities/interface/isFinite.h"
0015
0016 float EcalRegressionData::seedLeftRightAsym() const {
0017   float eLeftRightSum = eLeft() + eRight();
0018   float eLeftRightDiff = eLeft() - eRight();
0019   return eLeftRightSum != 0 ? eLeftRightDiff / eLeftRightSum : 0.;
0020 }
0021
0022 float EcalRegressionData::seedTopBottomAsym() const {
0023   float eTopBottomSum = eTop() + eBottom();
0024   float eTopBottomDiff = eTop() - eBottom();
0025   return eTopBottomSum != 0 ? eTopBottomDiff / eTopBottomSum : 0.;
0026 }
0027
0028 float EcalRegressionData::subClusRawEnergy(size_t clusNr) const {
0029   if (clusNr < subClusRawEnergy_.size())
0030     return subClusRawEnergy_[clusNr];
0031   else
0032     return 0.;
0033 }
0034
0035 float EcalRegressionData::subClusDEta(size_t clusNr) const {
0036   if (clusNr < subClusDEta_.size())
0037     return subClusDEta_[clusNr];
0038   else
0039     return 0.;
0040 }
0041
0042 float EcalRegressionData::subClusDPhi(size_t clusNr) const {
0043   if (clusNr < subClusDPhi_.size())
0044     return subClusDPhi_[clusNr];
0045   else
0046     return 0.;
0047 }
0048
0049 void EcalRegressionData::fill(const reco::SuperCluster& superClus,
0050                               const EcalRecHitCollection* ebRecHits,
0051                               const EcalRecHitCollection* eeRecHits,
0052                               const CaloGeometry* geom,
0053                               const CaloTopology* topology,
0054                               int nrVertices) {
0055   clear();
0056
0057   const DetId& seedid = superClus.seed()->hitsAndFractions().at(0).first;
0058   isEB_ = (seedid.subdetId() == EcalBarrel);
0059
0060   // skip HGCal
0061   if (EcalTools::isHGCalDet(seedid.det()))
0062     return;
0063
0064   const EcalRecHitCollection* recHits = isEB_ ? ebRecHits : eeRecHits;
0065
0066   scRawEnergy_ = superClus.rawEnergy();
0067   scCalibEnergy_ = superClus.correctedEnergy();
0068   scPreShowerEnergy_ = superClus.preshowerEnergy();
0069   scEta_ = superClus.eta();
0070   scPhi_ = superClus.phi();
0071   scEtaWidth_ = superClus.etaWidth();
0072   scPhiWidth_ = superClus.phiWidth();
0073   scNrAdditionalClusters_ = superClus.clustersSize() - 1;
0074
0075   seedClusEnergy_ = superClus.seed()->energy();
0076   eMax_ = EcalClusterTools::eMax(*superClus.seed(), recHits);
0077   e2nd_ = EcalClusterTools::e2nd(*superClus.seed(), recHits);
0078   e3x3_ = EcalClusterTools::e3x3(*superClus.seed(), recHits, topology);
0079   eTop_ = EcalClusterTools::eTop(*superClus.seed(), recHits, topology);
0080   eBottom_ = EcalClusterTools::eBottom(*superClus.seed(), recHits, topology);
0081   eLeft_ = EcalClusterTools::eLeft(*superClus.seed(), recHits, topology);
0082   eRight_ = EcalClusterTools::eRight(*superClus.seed(), recHits, topology);
0083   const auto& localCovs = EcalClusterTools::localCovariances(*superClus.seed(), recHits, topology);
0084   sigmaIEtaIEta_ = edm::isNotFinite(localCovs[0]) ? 0. : std::sqrt(localCovs[0]);
0085   sigmaIPhiIPhi_ = edm::isNotFinite(localCovs[2]) ? 0. : std::sqrt(localCovs[2]);
0086
0087   if (sigmaIEtaIEta_ * sigmaIPhiIPhi_ > 0)
0088     sigmaIEtaIPhi_ = localCovs[1] / (sigmaIEtaIEta_ * sigmaIPhiIPhi_);
0089   else if (localCovs[1] > 0)
0090     sigmaIEtaIPhi_ = 1.;
0091   else
0092     sigmaIEtaIPhi_ = -1.;
0093
0094   auto localCoordsFunc = isEB() ? egammaTools::localEcalClusterCoordsEB : egammaTools::localEcalClusterCoordsEE;
0095
0096   float dummy;
0097   localCoordsFunc(
0098       *superClus.seed(), *geom, seedCrysEtaOrX_, seedCrysPhiOrY_, seedCrysIEtaOrIX_, seedCrysIPhiOrIY_, dummy, dummy);
0099
0100   for (auto clus = superClus.clustersBegin() + 1; clus != superClus.clustersEnd(); ++clus) {
0101     const float dEta = (*clus)->eta() - superClus.seed()->eta();
0102     const float dPhi = reco::deltaPhi((*clus)->phi(), superClus.seed()->phi());
0103     const float dR2 = dEta * dEta + dPhi * dPhi;
0104     if (dR2 > maxSubClusDR2_ || maxSubClusDR2_ == 998001.) {
0105       maxSubClusDR2_ = dR2;
0106       maxSubClusDRDEta_ = dEta;
0107       maxSubClusDRDPhi_ = dPhi;
0108       maxSubClusDRRawEnergy_ = (*clus)->energy();
0109     }
0110     subClusRawEnergy_.push_back((*clus)->energy());
0111     subClusDEta_.push_back(dEta);
0112     subClusDPhi_.push_back(dPhi);
0113   }
0114
0115   nrVtx_ = nrVertices;
0116 }
0117
0118 void EcalRegressionData::clear() {
0119   isEB_ = false;
0120   scRawEnergy_ = 0.;
0121   scCalibEnergy_ = 0.;
0122   scPreShowerEnergy_ = 0.;
0123   scEta_ = 0.;
0124   scPhi_ = 0.;
0125   scEtaWidth_ = 0.;
0126   scPhiWidth_ = 0.;
0127   scNrAdditionalClusters_ = 0;
0128
0129   seedClusEnergy_ = 0.;
0130   eMax_ = 0.;
0131   e2nd_ = 0.;
0132   e3x3_ = 0.;
0133   eTop_ = 0.;
0134   eBottom_ = 0.;
0135   eLeft_ = 0.;
0136   eRight_ = 0.;
0137   sigmaIEtaIEta_ = 0.;
0138   sigmaIEtaIPhi_ = 0.;
0139   sigmaIPhiIPhi_ = 0.;
0140
0141   seedCrysPhiOrY_ = 0.;
0142   seedCrysEtaOrX_ = 0.;
0143   seedCrysIEtaOrIX_ = 0;
0144   seedCrysIPhiOrIY_ = 0;
0145
0146   maxSubClusDR2_ = 998001.;
0147   maxSubClusDRDPhi_ = 999.;
0148   maxSubClusDRDEta_ = 999;
0149   maxSubClusDRRawEnergy_ = 0.;
0150
0151   subClusRawEnergy_.clear();
0152   subClusDPhi_.clear();
0153   subClusDEta_.clear();
0154
0155   nrVtx_ = 0;
0156 }
0157
0158 void EcalRegressionData::fillVec(std::vector<float>& inputVec) const {
0159   if (isEB())
0160     fillVecEB_(inputVec);
0161   else
0162     fillVecEE_(inputVec);
0163 }
0164
0165 void EcalRegressionData::fillVecEB_(std::vector<float>& inputVec) const {
0166   inputVec.clear();
0167   inputVec.resize(33);
0168   inputVec[0] = nrVtx();                                          //nVtx
0169   inputVec[1] = scEta();                                          //scEta
0170   inputVec[2] = scPhi();                                          //scPhi
0171   inputVec[3] = scEtaWidth();                                     //scEtaWidth
0172   inputVec[4] = scPhiWidth();                                     //scPhiWidth
0173   inputVec[5] = scSeedR9();                                       //scSeedR9
0174   inputVec[6] = seedClusEnergyOverSCRawEnergy();                  //scSeedRawEnergy/scRawEnergy
0175   inputVec[7] = eMaxOverSCRawEnergy();                            //scSeedEmax/scRawEnergy
0176   inputVec[8] = e2ndOverSCRawEnergy();                            //scSeedE2nd/scRawEnergy
0177   inputVec[9] = seedLeftRightAsym();                              //scSeedLeftRightAsym
0178   inputVec[10] = seedTopBottomAsym();                             //scSeedTopBottomAsym
0179   inputVec[11] = sigmaIEtaIEta();                                 //scSeedSigmaIetaIeta
0180   inputVec[12] = sigmaIEtaIPhi();                                 //scSeedSigmaIetaIphi
0181   inputVec[13] = sigmaIPhiIPhi();                                 //scSeedSigmaIphiIphi
0182   inputVec[14] = scNrAdditionalClusters();                        //N_ECALClusters
0183   inputVec[15] = maxSubClusDR();                                  //clusterMaxDR
0184   inputVec[16] = maxSubClusDRDPhi();                              //clusterMaxDRDPhi
0185   inputVec[17] = maxSubClusDRDEta();                              //clusterMaxDRDEta
0186   inputVec[18] = maxSubClusDRRawEnergyOverSCRawEnergy();          //clusMaxDRRawEnergy/scRawEnergy
0187   inputVec[19] = subClusRawEnergyOverSCRawEnergy(SubClusNr::C1);  //clusterRawEnergy[0]/scRawEnergy
0188   inputVec[20] = subClusRawEnergyOverSCRawEnergy(
0189       SubClusNr::C2);  //clusterRawEnergy[1]/scRawEnergy  float scPreShowerEnergy()const{return scPreShowerEnergy_;}
0190
0191   inputVec[21] = subClusRawEnergyOverSCRawEnergy(SubClusNr::C3);  //clusterRawEnergy[2]/scRawEnergy
0192   inputVec[22] = subClusDPhi(SubClusNr::C1);                      //clusterDPhiToSeed[0]
0193   inputVec[23] = subClusDPhi(SubClusNr::C2);                      //clusterDPhiToSeed[1]
0194   inputVec[24] = subClusDPhi(SubClusNr::C3);                      //clusterDPhiToSeed[2]
0195   inputVec[25] = subClusDEta(SubClusNr::C1);                      //clusterDEtaToSeed[0]
0196   inputVec[26] = subClusDEta(SubClusNr::C2);                      //clusterDEtaToSeed[1]
0197   inputVec[27] = subClusDEta(SubClusNr::C3);                      //clusterDEtaToSeed[2]
0198   inputVec[28] = seedCrysEtaOrX();                                //scSeedCryEta
0199   inputVec[29] = seedCrysPhiOrY();                                //scSeedCryPhi
0200   inputVec[30] = seedCrysIEtaOrIX();                              //scSeedCryIeta
0201   inputVec[31] = seedCrysIPhiOrIY();                              //scSeedCryIphi
0202   inputVec[32] = scCalibEnergy();                                 //scCalibratedEnergy
0203 }
0204
0205 void EcalRegressionData::fillVecEE_(std::vector<float>& inputVec) const {
0206   inputVec.clear();
0207   inputVec.resize(33);    //this emulates the old behaviour of RegressionHelper, even if past 29 we dont use elements
0208   inputVec[0] = nrVtx();  //nVtx
0209   inputVec[1] = scEta();  //scEta
0210   inputVec[2] = scPhi();  //scPhi
0211   inputVec[3] = scEtaWidth();                                     //scEtaWidth
0212   inputVec[4] = scPhiWidth();                                     //scPhiWidth
0213   inputVec[5] = scSeedR9();                                       //scSeedR9
0214   inputVec[6] = seedClusEnergyOverSCRawEnergy();                  //scSeedRawEnergy/scRawEnergy
0215   inputVec[7] = eMaxOverSCRawEnergy();                            //scSeedEmax/scRawEnergy
0216   inputVec[8] = e2ndOverSCRawEnergy();                            //scSeedE2nd/scRawEnergy
0217   inputVec[9] = seedLeftRightAsym();                              //scSeedLeftRightAsym
0218   inputVec[10] = seedTopBottomAsym();                             //scSeedTopBottomAsym
0219   inputVec[11] = sigmaIEtaIEta();                                 //scSeedSigmaIetaIeta
0220   inputVec[12] = sigmaIEtaIPhi();                                 //scSeedSigmaIetaIphi
0221   inputVec[13] = sigmaIPhiIPhi();                                 //scSeedSigmaIphiIphi
0222   inputVec[14] = scNrAdditionalClusters();                        //N_ECALClusters
0223   inputVec[15] = maxSubClusDR();                                  //clusterMaxDR
0224   inputVec[16] = maxSubClusDRDPhi();                              //clusterMaxDRDPhi
0225   inputVec[17] = maxSubClusDRDEta();                              //clusterMaxDRDEta
0226   inputVec[18] = maxSubClusDRRawEnergyOverSCRawEnergy();          //clusMaxDRRawEnergy/scRawEnergy
0227   inputVec[19] = subClusRawEnergyOverSCRawEnergy(SubClusNr::C1);  //clusterRawEnergy[0]/scRawEnergy
0228   inputVec[20] = subClusRawEnergyOverSCRawEnergy(SubClusNr::C2);  //clusterRawEnergy[1]/scRawEnergy
0229   inputVec[21] = subClusRawEnergyOverSCRawEnergy(SubClusNr::C3);  //clusterRawEnergy[2]/scRawEnergy
0230   inputVec[22] = subClusDPhi(SubClusNr::C1);                      //clusterDPhiToSeed[0]
0231   inputVec[23] = subClusDPhi(SubClusNr::C2);                      //clusterDPhiToSeed[1]
0232   inputVec[24] = subClusDPhi(SubClusNr::C3);                      //clusterDPhiToSeed[2]
0233   inputVec[25] = subClusDEta(SubClusNr::C1);                      //clusterDEtaToSeed[0]
0234   inputVec[26] = subClusDEta(SubClusNr::C2);                      //clusterDEtaToSeed[1]
0235   inputVec[27] = subClusDEta(SubClusNr::C3);                      //clusterDEtaToSeed[2]
0236   inputVec[28] = scPreShowerEnergyOverSCRawEnergy();
0237   inputVec[29] = scCalibEnergy();  //scCalibratedEnergy
0238 }