PhotonIdentification/src/PhotonMIPHaloTagger.cc

0001 /** \class PhotonMIPHaloTragger
0002  *  Determine and Set whether a haol or not
0003  *
0004  *  \author Sushil S. Chauhan, A. Nowack, M. Tripathi : (UCDavis)
0005  *  \author A. Askew (FSU)
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 PhotonMIPHaloTagger::PhotonMIPHaloTagger(const edm::ParameterSet& conf, edm::ConsumesCollector&& iC)
0016     : EBecalCollection_{iC.consumes<EcalRecHitCollection>(
0017           conf.getParameter<edm::InputTag>("barrelEcalRecHitCollection"))},
0018       yRangeFit_{conf.getParameter<double>("YRangeFit")},
0019       xRangeFit_{conf.getParameter<double>("XRangeFit")},
0020       residualWidthEnergy_{conf.getParameter<double>("ResidualWidth")},
0021       haloDiscThreshold_{conf.getParameter<double>("HaloDiscThreshold")} {}
0022
0023 reco::Photon::MIPVariables PhotonMIPHaloTagger::mipCalculate(const reco::Photon& pho,
0024                                                              const edm::Event& e,
0025                                                              const edm::EventSetup& es) const {
0026   // Get EcalRecHits
0027   edm::Handle<EcalRecHitCollection> barrelHitHandle;
0028   e.getByToken(EBecalCollection_, barrelHitHandle);
0029
0030   bool validEcalRecHits = barrelHitHandle.isValid();
0031
0032   reco::Photon::MIPVariables mipId;
0033   if (validEcalRecHits) {
0034     mipId =
0035         getMipTrailFit(pho, e, es, barrelHitHandle, yRangeFit_, xRangeFit_, residualWidthEnergy_, haloDiscThreshold_);
0036   } else {
0037     edm::LogError("MIPcalculate") << "Error! Can't get the barrel hits product, hence set some default values";
0038     mipId.mipChi2 = -999.;
0039     mipId.mipTotEnergy = -999.;
0040     mipId.mipSlope = -999.;
0041     mipId.mipIntercept = -999.;
0042     mipId.mipNhitCone = 0;
0043     mipId.mipIsHalo = false;
0044   }
0045
0046   // std::cout<<"Setting: isHalo = "<<ismipHalo_<<"    nhitcone=  "<<nhitCone_<<std::endl;
0047   // std::cout<<"Setting: Chi2  = "<<mipFitResults_[0]<<"  eT= "<<mipFitResults_[1]<<" slope = "<<mipFitResults_[2]<<"  Intercept ="<<mipFitResults_[3]<<std::endl;
0048   return mipId;
0049 }
0050
0051 // Get the MIP Variable NCone, Energy etc here
0052 reco::Photon::MIPVariables PhotonMIPHaloTagger::getMipTrailFit(const reco::Photon& photon,
0053                                                                const edm::Event& iEvent,
0054                                                                const edm::EventSetup& iSetup,
0055                                                                edm::Handle<EcalRecHitCollection> ecalhitsCollEB,
0056                                                                double inputRangeY,
0057                                                                double inputRangeX,
0058                                                                double inputResWidth,
0059                                                                double inputHaloDiscCut) const {
0060   constexpr bool printDebug = false;
0061
0062   if constexpr (printDebug)
0063     std::cout << " inside MipFitTrail " << std::endl;
0064   //getting them from cfi config
0065   const double yrange = inputRangeY;
0066   const double xrange = inputRangeX;
0067   const double res_width = inputResWidth;         //width of the residual distribution
0068   const double halo_disc_cut = inputHaloDiscCut;  //cut based on  Shower,Angle and MIP
0069
0070   //initilize them here
0071   const double m = yrange / xrange;  //slope of the lines which form the cone around the trail
0072
0073   //first get the seed cell index
0074   int seedIEta = -999;
0075   int seedIPhi = -999;
0076   double seedE = -999.;
0077
0078   //get seed propert.
0079   getSeedHighestE(photon, iEvent, iSetup, ecalhitsCollEB, seedIEta, seedIPhi, seedE);
0080
0081   if constexpr (printDebug)
0082     std::cout << "Seed E =" << seedE << "  Seed ieta = " << seedIEta << "   Seed iphi =" << seedIPhi << std::endl;
0083
0084   //create some vector
0085   std::vector<int> ieta_cell;
0086   std::vector<int> iphi_cell;
0087   std::vector<double> energy_cell;
0088
0089   for (auto const& hit : *ecalhitsCollEB) {
0090     const EBDetId dit = hit.detid();
0091
0092     const int iphicell = dit.iphi();
0093     int ietacell = dit.ieta();
0094
0095     if (ietacell < 0)
0096       ietacell++;
0097
0098     //Exclude all cells within +/- 5 ieta of seed cell
0099     if (std::abs(ietacell - seedIEta) >= 5 && hit.energy() > 0.) {
0100       int delt_ieta = ietacell - seedIEta;
0101       int delt_iphi = iphicell - seedIPhi;
0102
0103       //Phi wrapping inclusion
0104       if (delt_iphi > 180) {
0105         delt_iphi = delt_iphi - 360;
0106       }
0107       if (delt_iphi < -180) {
0108         delt_iphi = delt_iphi + 360;
0109       }
0110
0111       //Condition to be within the cones
0112       if (((delt_iphi >= (m * delt_ieta)) && (delt_iphi <= (-m * delt_ieta))) ||
0113           ((delt_iphi <= (m * delt_ieta)) && (delt_iphi >= (-m * delt_ieta)))) {
0114         ieta_cell.push_back(delt_ieta);
0115         iphi_cell.push_back(delt_iphi);
0116         energy_cell.push_back(hit.energy());
0117       }
0118
0119     }  //within certain range of seed cell
0120
0121   }  //loop over hits
0122
0123   //Iterations for improvements
0124
0125   int Npoints = ieta_cell.size();
0126   int throwaway_index = -1;
0127   double chi2 = 0.0;
0128   double eT = 0.;
0129   double a1 = 0.;
0130   double b1 = 0.;
0131   double hres = 99999.;
0132
0133   if constexpr (printDebug)
0134     std::cout << " starting npoints = " << Npoints << std::endl;
0135
0136   //start Iterations
0137   for (int it = 0; it < 200 && hres > (5.0 * res_width); it++) {  //Max iter. is 200
0138
0139     //Throw away previous highest residual, if not first loop
0140     if (throwaway_index != -1) {
0141       ieta_cell.erase(ieta_cell.begin() + throwaway_index);
0142       iphi_cell.erase(iphi_cell.begin() + throwaway_index);
0143       energy_cell.erase(energy_cell.begin() + throwaway_index);
0144       Npoints--;
0145     }
0146
0147     //defined some variable for iterative fitting the mip trail line
0148     double sx = 0.0;
0149     double sy = 0.0;
0150     double ss = 0.0;
0151     double sxx = 0.0;
0152     double sxy = 0.0;
0153     double m_chi2 = 0.0;
0154     double etot_cell = 0.0;
0155
0156     //Fit the line to trail
0157     for (int j = 0; j < Npoints; j++) {
0158       constexpr double wt = 1.0;
0159       ss += wt;
0160       sx += ieta_cell[j] * wt;
0161       sy += iphi_cell[j];
0162       sxx += ieta_cell[j] * ieta_cell[j] * wt;
0163       sxy += ieta_cell[j] * iphi_cell[j] * wt;
0164     }
0165
0166     const double delt = ss * sxx - (sx * sx);
0167     a1 = ((sxx * sy) - (sx * sxy)) / delt;  // INTERCEPT
0168     b1 = ((ss * sxy) - (sx * sy)) / delt;   // SLOPE
0169
0170     double highest_res = 0.;
0171     int highres_index = 0;
0172
0173     for (int j = 0; j < Npoints; j++) {
0174       const double res = 1.0 * iphi_cell[j] - a1 - b1 * ieta_cell[j];
0175       const double res_sq = res * res;
0176
0177       if (std::abs(res) > highest_res) {
0178         highest_res = std::abs(res);
0179         highres_index = j;
0180       }
0181
0182       m_chi2 += res_sq;
0183       etot_cell += energy_cell[j];
0184     }
0185
0186     throwaway_index = highres_index;
0187     hres = highest_res;
0188
0189     chi2 = m_chi2 / ((Npoints - 2));
0190     chi2 = chi2 / (res_width * res_width);
0191     eT = etot_cell;
0192
0193   }  //for loop for iterations
0194
0195   if constexpr (printDebug)
0196     std::cout << "hres = " << hres << std::endl;
0197
0198   //get roundness and angle for this photon candidate form EcalClusterTool
0199   const std::vector<float> showershapes_barrel =
0200       EcalClusterTools::roundnessBarrelSuperClusters(*(photon.superCluster()), (*ecalhitsCollEB.product()));
0201
0202   const double roundness = showershapes_barrel[0];
0203   const double angle = showershapes_barrel[1];
0204
0205   if constexpr (printDebug)
0206     std::cout << " eTot =" << eT << "     Rounness = " << roundness << "    angle  " << angle << std::endl;
0207
0208   //get the halo disc variable
0209   const double halo_disc = eT / (roundness * angle);
0210
0211   ///Now Fill the FitResults vector
0212   reco::Photon::MIPVariables results;
0213   results.mipChi2 = chi2;
0214   results.mipTotEnergy = eT;
0215   results.mipSlope = b1;
0216   results.mipIntercept = a1;
0217   results.mipNhitCone = Npoints;
0218   results.mipIsHalo = halo_disc > halo_disc_cut;
0219   //is halo?, yes if halo_disc > 70 by default
0220   // based on 2010 study
0221
0222   if constexpr (printDebug)
0223     std::cout << "Endof MIP Trail: halo_dic= " << halo_disc << "   nhitcone =" << results.mipNhitCone
0224               << "  isHalo= " << results.mipIsHalo << std::endl;
0225   if constexpr (printDebug)
0226     std::cout << "Endof MIP Trail: Chi2  = " << chi2 << "  eT= " << eT << " slope = " << b1 << "  Intercept =" << a1
0227               << std::endl;
0228
0229   return results;
0230 }
0231
0232 //get the seed crystal index
0233 void PhotonMIPHaloTagger::getSeedHighestE(const reco::Photon& photon,
0234                                           const edm::Event& iEvent,
0235                                           const edm::EventSetup& iSetup,
0236                                           edm::Handle<EcalRecHitCollection> Brechit,
0237                                           int& seedIeta,
0238                                           int& seedIphi,
0239                                           double& seedEnergy) const {
0240   constexpr bool printDebug = false;
0241
0242   if constexpr (printDebug)
0243     std::cout << "Inside GetSeed" << std::endl;
0244   //Get the Seed
0245   double SeedE = -999.;
0246
0247   //initilaze them here
0248   seedIeta = -999;
0249   seedIphi = -999;
0250   seedEnergy = -999.;
0251
0252   std::vector<std::pair<DetId, float> > const& PhotonHit_DetIds = photon.superCluster()->hitsAndFractions();
0253   for (auto pr : PhotonHit_DetIds) {
0254     if ((pr.first).det() == DetId::Ecal && (pr.first).subdetId() == EcalBarrel) {
0255       EcalRecHitCollection::const_iterator thishit = Brechit->find((pr.first));
0256
0257       if (thishit == Brechit->end()) {
0258         continue;
0259       }
0260
0261       const EBDetId detId{pr.first};
0262
0263       const double crysE = thishit->energy();
0264
0265       if (crysE > SeedE) {
0266         SeedE = crysE;
0267         seedIeta = detId.ieta();
0268         seedIphi = (detId.iphi());
0269         seedEnergy = SeedE;
0270
0271         if constexpr (printDebug)
0272           std::cout << "Current max Seed = " << SeedE << "   seedIphi = " << seedIphi << "  ieta= " << seedIeta
0273                     << std::endl;
0274       }
0275
0276     }  //check if in Barrel
0277
0278   }  //loop over EBrechits cells
0279
0280   if constexpr (printDebug)
0281     std::cout << "End of  GetSeed" << std::endl;
0282 }