Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-07-07 22:33:56

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   double energy_total = 0.;
0126   int delt_ieta = 0;
0127   int delt_iphi = 0;
0128 
0129   for (EcalRecHitCollection::const_iterator it = ecalhitsCollEB->begin(); it != ecalhitsCollEB->end(); ++it) {
0130     EBDetId dit = it->detid();
0131 
0132     iphicell = dit.iphi();
0133     ietacell = dit.ieta();
0134 
0135     if (ietacell < 0)
0136       ietacell++;
0137 
0138     //Exclude all cells within +/- 5 ieta of seed cell
0139     if (std::abs(ietacell - seedIEta) >= 5 && it->energy() > 0.) {
0140       delt_ieta = ietacell - seedIEta;
0141       delt_iphi = iphicell - seedIPhi;
0142 
0143       //Phi wrapping inclusion
0144       if (delt_iphi > 180) {
0145         delt_iphi = delt_iphi - 360;
0146       }
0147       if (delt_iphi < -180) {
0148         delt_iphi = delt_iphi + 360;
0149       }
0150 
0151       //Condition to be within the cones
0152       if (((delt_iphi >= (m * delt_ieta)) && (delt_iphi <= (-m * delt_ieta))) ||
0153           ((delt_iphi <= (m * delt_ieta)) && (delt_iphi >= (-m * delt_ieta)))) {
0154         ieta_cell.push_back(delt_ieta);
0155         iphi_cell.push_back(delt_iphi);
0156         energy_cell.push_back(it->energy());
0157         energy_total += it->energy();
0158         kArray++;
0159       }
0160 
0161     }  //within cerntain range of seed cell
0162 
0163   }  //loop voer hits
0164 
0165   //Iterations for imporovements
0166 
0167   int Npoints = 0;
0168   int throwaway_index = -1;
0169   double chi2;
0170   double eT = 0.;
0171   double hres = 99999.;
0172 
0173   //some tmp local variale
0174   double Roundness_ = 999.;
0175   double Angle_ = 999.;
0176   double halo_disc_ = 0.;
0177 
0178   Npoints = kArray;
0179 
0180   if (debug_)
0181     std::cout << " starting npoing = " << Npoints << std::endl;
0182 
0183   //defined some variable for iterative fitting the mip trail line
0184   double sx = 0.0;
0185   double sy = 0.0;
0186   double ss = 0.0;
0187   double sxx = 0.0;
0188   double sxy = 0.0;
0189   double a1 = 0.0;
0190   double b1 = 0.0;
0191   double m_chi2 = 0.0;
0192   double etot_cell = 0.0;
0193 
0194   //start Iterations
0195   for (int it = 0; it < 200 && hres > (5.0 * res_width); it++) {  //Max iter. is 200
0196 
0197     //Throw away previous highest residual, if not first loop
0198     if (throwaway_index != -1) {
0199       ieta_cell.erase(ieta_cell.begin() + throwaway_index);
0200       iphi_cell.erase(iphi_cell.begin() + throwaway_index);
0201       energy_cell.erase(energy_cell.begin() + throwaway_index);
0202       Npoints--;
0203     }
0204 
0205     //Lets Initialize them first for each iteration
0206     sx = 0.0;
0207     sy = 0.0;
0208     ss = 0.0;
0209     sxx = 0.0;
0210     sxy = 0.0;
0211     m_chi2 = 0.0;
0212     etot_cell = 0.0;
0213 
0214     //Fit the line to trail
0215     for (int j = 0; j < Npoints; j++) {
0216       double wt = 1.0;
0217       ss += wt;
0218       sx += ieta_cell[j] * wt;
0219       sy += iphi_cell[j];
0220       sxx += ieta_cell[j] * ieta_cell[j] * wt;
0221       sxy += ieta_cell[j] * iphi_cell[j] * wt;
0222     }
0223 
0224     double delt = ss * sxx - (sx * sx);
0225     a1 = ((sxx * sy) - (sx * sxy)) / delt;  // INTERCEPT
0226     b1 = ((ss * sxy) - (sx * sy)) / delt;   // SLOPE
0227 
0228     double highest_res = 0.;
0229     int highres_index = 0;
0230 
0231     for (int j = 0; j < Npoints; j++) {
0232       double res = 1.0 * iphi_cell[j] - a1 - b1 * ieta_cell[j];
0233       double res_sq = res * res;
0234 
0235       if (std::abs(res) > highest_res) {
0236         highest_res = std::abs(res);
0237         highres_index = j;
0238       }
0239 
0240       m_chi2 += res_sq;
0241       etot_cell += energy_cell[j];
0242     }
0243 
0244     throwaway_index = highres_index;
0245     hres = highest_res;
0246 
0247     chi2 = m_chi2 / ((Npoints - 2));
0248     chi2 = chi2 / (res_width * res_width);
0249     eT = etot_cell;
0250 
0251   }  //for loop for iterations
0252 
0253   if (debug_)
0254     std::cout << "hres = " << hres << std::endl;
0255 
0256   //get roundness and angle for this photon candidate form EcalClusterTool
0257   std::vector<float> showershapes_barrel =
0258       EcalClusterTools::roundnessBarrelSuperClusters(*(photon->superCluster()), (*ecalhitsCollEB.product()));
0259 
0260   Roundness_ = showershapes_barrel[0];
0261   Angle_ = showershapes_barrel[1];
0262 
0263   if (debug_)
0264     std::cout << " eTot =" << eT << "     Rounness = " << Roundness_ << "    Angle_  " << Angle_ << std::endl;
0265 
0266   //get the halo disc variable
0267   halo_disc_ = eT / (Roundness_ * Angle_);
0268 
0269   ///Now Filll the FitResults vector
0270   FitResults_.push_back(chi2);  //chi2
0271   FitResults_.push_back(eT);    //total energy
0272   FitResults_.push_back(b1);    //slope
0273   FitResults_.push_back(a1);    //intercept
0274   NhitCone_ = Npoints;          //nhit in cone
0275   if (halo_disc_ > halo_disc_cut)
0276     ismipHalo_ = true;  //is halo?, yes if halo_disc > 70 by default
0277                         // based on 2010 study
0278 
0279   if (debug_)
0280     std::cout << "Endof MIP Trail: halo_dic= " << halo_disc_ << "   nhitcone =" << NhitCone_
0281               << "  isHalo= " << ismipHalo_ << std::endl;
0282   if (debug_)
0283     std::cout << "Endof MIP Trail: Chi2  = " << chi2 << "  eT= " << eT << " slope = " << b1 << "  Intercept =" << a1
0284               << std::endl;
0285 
0286   return FitResults_;
0287 }
0288 
0289 //get the seed crystal index
0290 void PhotonMIPHaloTagger::GetSeedHighestE(const reco::Photon* photon,
0291                                           const edm::Event& iEvent,
0292                                           const edm::EventSetup& iSetup,
0293                                           edm::Handle<EcalRecHitCollection> Brechit,
0294                                           int& seedIeta,
0295                                           int& seedIphi,
0296                                           double& seedEnergy) {
0297   bool debug_ = false;
0298 
0299   if (debug_)
0300     std::cout << "Inside GetSeed" << std::endl;
0301   //Get the Seed
0302   double SeedE = -999.;
0303 
0304   //initilaze them here
0305   seedIeta = -999;
0306   seedIphi = -999;
0307   seedEnergy = -999.;
0308 
0309   std::vector<std::pair<DetId, float> > PhotonHit_DetIds = photon->superCluster()->hitsAndFractions();
0310   int ncrys = 0;
0311 
0312   std::vector<std::pair<DetId, float> >::const_iterator detitr;
0313   for (detitr = PhotonHit_DetIds.begin(); detitr != PhotonHit_DetIds.end(); ++detitr) {
0314     if (((*detitr).first).det() == DetId::Ecal && ((*detitr).first).subdetId() == EcalBarrel) {
0315       EcalRecHitCollection::const_iterator j = Brechit->find(((*detitr).first));
0316       EcalRecHitCollection::const_iterator thishit;
0317 
0318       if (j != Brechit->end())
0319         thishit = j;
0320       if (j == Brechit->end()) {
0321         continue;
0322       }
0323 
0324       EBDetId detId = (EBDetId)((*detitr).first);
0325 
0326       double crysE = thishit->energy();
0327 
0328       if (crysE > SeedE) {
0329         SeedE = crysE;
0330         seedIeta = detId.ieta();
0331         seedIphi = (detId.iphi());
0332         seedEnergy = SeedE;
0333 
0334         if (debug_)
0335           std::cout << "Current max Seed = " << SeedE << "   seedIphi = " << seedIphi << "  ieta= " << seedIeta
0336                     << std::endl;
0337       }
0338 
0339       ncrys++;
0340 
0341     }  //check if in Barrel
0342 
0343   }  //loop over EBrechits cells
0344 
0345   if (debug_)
0346     std::cout << "End of  GetSeed" << std::endl;
0347 }