File indexing completed on 2023-05-26 22:38:23
0001
0002
0003
0004
0005
0006
0007
0008 #include "RecoEgamma/PhotonIdentification/interface/PhotonMIPHaloTagger.h"
0009 #include "DataFormats/DetId/interface/DetId.h"
0010 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0011 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0014
0015 void PhotonMIPHaloTagger::setup(const edm::ParameterSet& conf, edm::ConsumesCollector&& iC) {
0016 EBecalCollection_ = iC.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("barrelEcalRecHitCollection"));
0017 EEecalCollection_ = iC.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("endcapEcalRecHitCollection"));
0018
0019 yRangeFit_ = conf.getParameter<double>("YRangeFit");
0020 xRangeFit_ = conf.getParameter<double>("XRangeFit");
0021 residualWidthEnergy_ = conf.getParameter<double>("ResidualWidth");
0022 haloDiscThreshold_ = conf.getParameter<double>("HaloDiscThreshold");
0023 }
0024
0025 void PhotonMIPHaloTagger::MIPcalculate(const reco::Photon* pho,
0026 const edm::Event& e,
0027 const edm::EventSetup& es,
0028 reco::Photon::MIPVariables& mipId) {
0029
0030 inputRangeY = yRangeFit_;
0031 inputRangeX = xRangeFit_;
0032 inputResWidth = residualWidthEnergy_;
0033 inputHaloDiscCut = haloDiscThreshold_;
0034
0035
0036 mipFitResults_.clear();
0037
0038 nhitCone_ = -99;
0039 ismipHalo_ = false;
0040
0041
0042 edm::Handle<EcalRecHitCollection> barrelHitHandle;
0043 e.getByToken(EBecalCollection_, barrelHitHandle);
0044
0045 bool validEcalRecHits = barrelHitHandle.isValid();
0046 if (validEcalRecHits) {
0047
0048 mipFitResults_ = GetMipTrailFit(
0049 pho, e, es, barrelHitHandle, inputRangeY, inputRangeX, inputResWidth, inputHaloDiscCut, nhitCone_, ismipHalo_);
0050
0051
0052 mipId.mipChi2 = mipFitResults_[0];
0053 mipId.mipTotEnergy = mipFitResults_[1];
0054 mipId.mipSlope = mipFitResults_[2];
0055 mipId.mipIntercept = mipFitResults_[3];
0056 mipId.mipNhitCone = nhitCone_;
0057 mipId.mipIsHalo = ismipHalo_;
0058 } else {
0059 edm::LogError("MIPcalculate") << "Error! Can't get the barrel hits product, hence set some default values";
0060 mipId.mipChi2 = -999.;
0061 mipId.mipTotEnergy = -999.;
0062 mipId.mipSlope = -999.;
0063 mipId.mipIntercept = -999.;
0064 mipId.mipNhitCone = 0;
0065 mipId.mipIsHalo = false;
0066 }
0067
0068
0069
0070 }
0071
0072
0073 std::vector<double> PhotonMIPHaloTagger::GetMipTrailFit(const reco::Photon* photon,
0074 const edm::Event& iEvent,
0075 const edm::EventSetup& iSetup,
0076 edm::Handle<EcalRecHitCollection> ecalhitsCollEB,
0077 double inputRangeY,
0078 double inputRangeX,
0079 double inputResWidth,
0080 double inputHaloDiscCut,
0081 int& NhitCone_,
0082 bool& ismipHalo_) {
0083 bool debug_ = false;
0084
0085 if (debug_)
0086 std::cout << " inside MipFitTrail " << std::endl;
0087
0088 double yrange = inputRangeY;
0089 double xrange = inputRangeX;
0090 double res_width = inputResWidth;
0091 double halo_disc_cut = inputHaloDiscCut;
0092
0093
0094 double m = 0.;
0095 m = yrange / xrange;
0096
0097
0098 int seedIEta = -999;
0099 int seedIPhi = -999;
0100 double seedE = -999.;
0101
0102
0103 GetSeedHighestE(photon, iEvent, iSetup, ecalhitsCollEB, seedIEta, seedIPhi, seedE);
0104
0105 if (debug_)
0106 std::cout << "Seed E =" << seedE << " Seed ieta = " << seedIEta << " Seed iphi =" << seedIPhi << std::endl;
0107
0108
0109 std::vector<double> FitResults_;
0110 FitResults_.clear();
0111
0112
0113 std::vector<int> ieta_cell;
0114 std::vector<int> iphi_cell;
0115 std::vector<double> energy_cell;
0116
0117 ieta_cell.clear();
0118 iphi_cell.clear();
0119 energy_cell.clear();
0120
0121 int ietacell = 0;
0122 int iphicell = 0;
0123 int kArray = 0;
0124
0125 int delt_ieta = 0;
0126 int delt_iphi = 0;
0127
0128 for (EcalRecHitCollection::const_iterator it = ecalhitsCollEB->begin(); it != ecalhitsCollEB->end(); ++it) {
0129 EBDetId dit = it->detid();
0130
0131 iphicell = dit.iphi();
0132 ietacell = dit.ieta();
0133
0134 if (ietacell < 0)
0135 ietacell++;
0136
0137
0138 if (std::abs(ietacell - seedIEta) >= 5 && it->energy() > 0.) {
0139 delt_ieta = ietacell - seedIEta;
0140 delt_iphi = iphicell - seedIPhi;
0141
0142
0143 if (delt_iphi > 180) {
0144 delt_iphi = delt_iphi - 360;
0145 }
0146 if (delt_iphi < -180) {
0147 delt_iphi = delt_iphi + 360;
0148 }
0149
0150
0151 if (((delt_iphi >= (m * delt_ieta)) && (delt_iphi <= (-m * delt_ieta))) ||
0152 ((delt_iphi <= (m * delt_ieta)) && (delt_iphi >= (-m * delt_ieta)))) {
0153 ieta_cell.push_back(delt_ieta);
0154 iphi_cell.push_back(delt_iphi);
0155 energy_cell.push_back(it->energy());
0156 kArray++;
0157 }
0158
0159 }
0160
0161 }
0162
0163
0164
0165 int Npoints = 0;
0166 int throwaway_index = -1;
0167 double chi2;
0168 double eT = 0.;
0169 double hres = 99999.;
0170
0171
0172 double Roundness_ = 999.;
0173 double Angle_ = 999.;
0174 double halo_disc_ = 0.;
0175
0176 Npoints = kArray;
0177
0178 if (debug_)
0179 std::cout << " starting npoing = " << Npoints << std::endl;
0180
0181
0182 double sx = 0.0;
0183 double sy = 0.0;
0184 double ss = 0.0;
0185 double sxx = 0.0;
0186 double sxy = 0.0;
0187 double a1 = 0.0;
0188 double b1 = 0.0;
0189 double m_chi2 = 0.0;
0190 double etot_cell = 0.0;
0191
0192
0193 for (int it = 0; it < 200 && hres > (5.0 * res_width); it++) {
0194
0195
0196 if (throwaway_index != -1) {
0197 ieta_cell.erase(ieta_cell.begin() + throwaway_index);
0198 iphi_cell.erase(iphi_cell.begin() + throwaway_index);
0199 energy_cell.erase(energy_cell.begin() + throwaway_index);
0200 Npoints--;
0201 }
0202
0203
0204 sx = 0.0;
0205 sy = 0.0;
0206 ss = 0.0;
0207 sxx = 0.0;
0208 sxy = 0.0;
0209 m_chi2 = 0.0;
0210 etot_cell = 0.0;
0211
0212
0213 for (int j = 0; j < Npoints; j++) {
0214 double wt = 1.0;
0215 ss += wt;
0216 sx += ieta_cell[j] * wt;
0217 sy += iphi_cell[j];
0218 sxx += ieta_cell[j] * ieta_cell[j] * wt;
0219 sxy += ieta_cell[j] * iphi_cell[j] * wt;
0220 }
0221
0222 double delt = ss * sxx - (sx * sx);
0223 a1 = ((sxx * sy) - (sx * sxy)) / delt;
0224 b1 = ((ss * sxy) - (sx * sy)) / delt;
0225
0226 double highest_res = 0.;
0227 int highres_index = 0;
0228
0229 for (int j = 0; j < Npoints; j++) {
0230 double res = 1.0 * iphi_cell[j] - a1 - b1 * ieta_cell[j];
0231 double res_sq = res * res;
0232
0233 if (std::abs(res) > highest_res) {
0234 highest_res = std::abs(res);
0235 highres_index = j;
0236 }
0237
0238 m_chi2 += res_sq;
0239 etot_cell += energy_cell[j];
0240 }
0241
0242 throwaway_index = highres_index;
0243 hres = highest_res;
0244
0245 chi2 = m_chi2 / ((Npoints - 2));
0246 chi2 = chi2 / (res_width * res_width);
0247 eT = etot_cell;
0248
0249 }
0250
0251 if (debug_)
0252 std::cout << "hres = " << hres << std::endl;
0253
0254
0255 std::vector<float> showershapes_barrel =
0256 EcalClusterTools::roundnessBarrelSuperClusters(*(photon->superCluster()), (*ecalhitsCollEB.product()));
0257
0258 Roundness_ = showershapes_barrel[0];
0259 Angle_ = showershapes_barrel[1];
0260
0261 if (debug_)
0262 std::cout << " eTot =" << eT << " Rounness = " << Roundness_ << " Angle_ " << Angle_ << std::endl;
0263
0264
0265 halo_disc_ = eT / (Roundness_ * Angle_);
0266
0267
0268 FitResults_.push_back(chi2);
0269 FitResults_.push_back(eT);
0270 FitResults_.push_back(b1);
0271 FitResults_.push_back(a1);
0272 NhitCone_ = Npoints;
0273 if (halo_disc_ > halo_disc_cut)
0274 ismipHalo_ = true;
0275
0276
0277 if (debug_)
0278 std::cout << "Endof MIP Trail: halo_dic= " << halo_disc_ << " nhitcone =" << NhitCone_
0279 << " isHalo= " << ismipHalo_ << std::endl;
0280 if (debug_)
0281 std::cout << "Endof MIP Trail: Chi2 = " << chi2 << " eT= " << eT << " slope = " << b1 << " Intercept =" << a1
0282 << std::endl;
0283
0284 return FitResults_;
0285 }
0286
0287
0288 void PhotonMIPHaloTagger::GetSeedHighestE(const reco::Photon* photon,
0289 const edm::Event& iEvent,
0290 const edm::EventSetup& iSetup,
0291 edm::Handle<EcalRecHitCollection> Brechit,
0292 int& seedIeta,
0293 int& seedIphi,
0294 double& seedEnergy) {
0295 bool debug_ = false;
0296
0297 if (debug_)
0298 std::cout << "Inside GetSeed" << std::endl;
0299
0300 double SeedE = -999.;
0301
0302
0303 seedIeta = -999;
0304 seedIphi = -999;
0305 seedEnergy = -999.;
0306
0307 std::vector<std::pair<DetId, float> > PhotonHit_DetIds = photon->superCluster()->hitsAndFractions();
0308 std::vector<std::pair<DetId, float> >::const_iterator detitr;
0309 for (detitr = PhotonHit_DetIds.begin(); detitr != PhotonHit_DetIds.end(); ++detitr) {
0310 if (((*detitr).first).det() == DetId::Ecal && ((*detitr).first).subdetId() == EcalBarrel) {
0311 EcalRecHitCollection::const_iterator j = Brechit->find(((*detitr).first));
0312 EcalRecHitCollection::const_iterator thishit;
0313
0314 if (j != Brechit->end())
0315 thishit = j;
0316 if (j == Brechit->end()) {
0317 continue;
0318 }
0319
0320 EBDetId detId = (EBDetId)((*detitr).first);
0321
0322 double crysE = thishit->energy();
0323
0324 if (crysE > SeedE) {
0325 SeedE = crysE;
0326 seedIeta = detId.ieta();
0327 seedIphi = (detId.iphi());
0328 seedEnergy = SeedE;
0329
0330 if (debug_)
0331 std::cout << "Current max Seed = " << SeedE << " seedIphi = " << seedIphi << " ieta= " << seedIeta
0332 << std::endl;
0333 }
0334
0335 }
0336
0337 }
0338
0339 if (debug_)
0340 std::cout << "End of GetSeed" << std::endl;
0341 }