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 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   //get the predefined variables
0030   inputRangeY = yRangeFit_;
0031   inputRangeX = xRangeFit_;
0032   inputResWidth = residualWidthEnergy_;
0033   inputHaloDiscCut = haloDiscThreshold_;
0034
0035   //First store in local variables
0036   mipFitResults_.clear();
0037
0038   nhitCone_ = -99;     //hit inside the cone
0039   ismipHalo_ = false;  // halo?
0040
0041   // Get EcalRecHits
0042   edm::Handle<EcalRecHitCollection> barrelHitHandle;
0043   e.getByToken(EBecalCollection_, barrelHitHandle);
0044
0045   bool validEcalRecHits = barrelHitHandle.isValid();
0046   if (validEcalRecHits) {
0047     // GetMIPTrailFit
0048     mipFitResults_ = GetMipTrailFit(
0049         pho, e, es, barrelHitHandle, inputRangeY, inputRangeX, inputResWidth, inputHaloDiscCut, nhitCone_, ismipHalo_);
0050
0051     //Now set the variable in "MIPVaraible"
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   // std::cout<<"Setting: isHalo = "<<ismipHalo_<<"    nhitcone=  "<<nhitCone_<<std::endl;
0069   // std::cout<<"Setting: Chi2  = "<<mipFitResults_[0]<<"  eT= "<<mipFitResults_[1]<<" slope = "<<mipFitResults_[2]<<"  Intercept ="<<mipFitResults_[3]<<std::endl;
0070 }
0071
0072 // Get the MIP Variable NCone, Energy etc here
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   //getting them from cfi config
0088   double yrange = inputRangeY;
0089   double xrange = inputRangeX;
0090   double res_width = inputResWidth;         //width of the residual distribution
0091   double halo_disc_cut = inputHaloDiscCut;  //cut based on  Shower,Angle and MIP
0092
0093   //initilize them here
0094   double m = 0.;
0095   m = yrange / xrange;  //slope of the lines which form the cone around the trail
0096
0097   //first get the seed cell index
0098   int seedIEta = -999;
0099   int seedIPhi = -999;
0100   double seedE = -999.;
0101
0102   //get seed propert.
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   //to store results
0109   std::vector<double> FitResults_;
0110   FitResults_.clear();
0111
0112   //create some vector and clear them
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     //Exclude all cells within +/- 5 ieta of seed cell
0138     if (std::abs(ietacell - seedIEta) >= 5 && it->energy() > 0.) {
0139       delt_ieta = ietacell - seedIEta;
0140       delt_iphi = iphicell - seedIPhi;
0141
0142       //Phi wrapping inclusion
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       //Condition to be within the cones
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     }  //within cerntain range of seed cell
0160
0161   }  //loop voer hits
0162
0163   //Iterations for imporovements
0164
0165   int Npoints = 0;
0166   int throwaway_index = -1;
0167   double chi2;
0168   double eT = 0.;
0169   double hres = 99999.;
0170
0171   //some tmp local variale
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   //defined some variable for iterative fitting the mip trail line
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   //start Iterations
0193   for (int it = 0; it < 200 && hres > (5.0 * res_width); it++) {  //Max iter. is 200
0194
0195     //Throw away previous highest residual, if not first loop
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     //Lets Initialize them first for each iteration
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     //Fit the line to trail
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;  // INTERCEPT
0224     b1 = ((ss * sxy) - (sx * sy)) / delt;   // SLOPE
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   }  //for loop for iterations
0250
0251   if (debug_)
0252     std::cout << "hres = " << hres << std::endl;
0253
0254   //get roundness and angle for this photon candidate form EcalClusterTool
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   //get the halo disc variable
0265   halo_disc_ = eT / (Roundness_ * Angle_);
0266
0267   ///Now Filll the FitResults vector
0268   FitResults_.push_back(chi2);  //chi2
0269   FitResults_.push_back(eT);    //total energy
0270   FitResults_.push_back(b1);    //slope
0271   FitResults_.push_back(a1);    //intercept
0272   NhitCone_ = Npoints;          //nhit in cone
0273   if (halo_disc_ > halo_disc_cut)
0274     ismipHalo_ = true;  //is halo?, yes if halo_disc > 70 by default
0275                         // based on 2010 study
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 //get the seed crystal index
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   //Get the Seed
0300   double SeedE = -999.;
0301
0302   //initilaze them here
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     }  //check if in Barrel
0336
0337   }  //loop over EBrechits cells
0338
0339   if (debug_)
0340     std::cout << "End of  GetSeed" << std::endl;
0341 }