89 bool validEcalRecHits=
true;
94 if (!barrelHitHandle.
isValid()) {
95 edm::LogError(
"MIPcalculate") <<
"Error! Can't get the barrel hits product ";
96 validEcalRecHits=
false;
100 if(validEcalRecHits){
118 else{
std::cout<<
"MIPCalculate::Ecal Collection is not there hence set some default values"<<std::endl;
142 double inputResWidth,
143 double inputHaloDiscCut,
149 if(debug_)
std::cout<<
" inside MipFitTrail "<<std::endl;
167 double seedE = -999.;
179 if(debug_)
std::cout<<
"Seed E ="<<seedE<<
" Seed ieta = "<<seedIEta<<
" Seed iphi ="<< seedIPhi<<std::endl;
182 std::vector<double> FitResults_;
187 std::vector<int> ieta_cell;
188 std::vector<int> iphi_cell;
189 std::vector<double> energy_cell;
200 double energy_total = 0.;
210 iphicell = dit.
iphi();
211 ietacell = dit.
ieta();
213 if(ietacell < 0)ietacell++;
216 if(
TMath::Abs(ietacell - seedIEta) >= 5 && it->energy() > 0.)
219 delt_ieta = ietacell - seedIEta;
220 delt_iphi = iphicell - seedIPhi;
223 if(delt_iphi > 180){ delt_iphi = delt_iphi - 360; }
224 if(delt_iphi < -180){ delt_iphi = delt_iphi + 360; }
227 if( ( (delt_iphi >= (m*delt_ieta)) && (delt_iphi <= (-m*delt_ieta)) ) ||
228 ( (delt_iphi <= (m*delt_ieta)) && (delt_iphi >= (-m*delt_ieta)) )
231 ieta_cell.push_back(delt_ieta);
232 iphi_cell.push_back(delt_iphi);
233 energy_cell.push_back(it->energy());
234 energy_total += it->energy();
247 int throwaway_index = -1;
250 double hres = 99999.;
253 double Roundness_ =999.;
255 double halo_disc_ =0.;
259 if(debug_)
std::cout<<
" starting npoing = "<<Npoints<<std::endl;
274 double etot_cell=0.0;
278 for(
int it=0; it < 200 && hres >(5.0*res_width); it++)
283 if (throwaway_index!=-1)
285 ieta_cell.erase(ieta_cell.begin()+throwaway_index);
286 iphi_cell.erase(iphi_cell.begin()+throwaway_index);
287 energy_cell.erase(energy_cell.begin()+throwaway_index);
309 for(
int j=0;
j<Npoints;
j++)
313 sx += ieta_cell[
j]*wt;
315 sxx += ieta_cell[
j]*ieta_cell[
j]*wt;
316 sxy += ieta_cell[
j]*iphi_cell[
j]*wt;
319 delt = ss*sxx - (sx*sx);
320 a1 = ((sxx*sy)-(sx*sxy))/delt;
321 b1 = ((ss*sxy)-(sx*sy))/delt;
324 double highest_res = 0.;
325 int highres_index = 0;
328 for(
int j=0;
j<Npoints;
j++)
330 res = 1.0*iphi_cell[
j] - a1 - b1*ieta_cell[
j];
340 etot_cell += energy_cell[
j];
345 throwaway_index = highres_index;
348 chi2 = m_chi2 /((Npoints-2));
349 chi2 = chi2/(res_width*res_width);
354 if(debug_)
std::cout<<
"hres = "<<hres<<std::endl;
357 std::vector<float> showershapes_barrel = EcalClusterTools::roundnessBarrelSuperClusters(*(photon->
superCluster()), (*ecalhitsCollEB.
product()));
359 Roundness_ = showershapes_barrel[0];
360 Angle_ = showershapes_barrel[1];
362 if(debug_)
std::cout<<
" eTot ="<<eT<<
" Rounness = "<<Roundness_<<
" Angle_ "<<Angle_ <<std::endl;
366 halo_disc_ = eT/(Roundness_* Angle_);
370 FitResults_.push_back(chi2);
371 FitResults_.push_back(eT);
372 FitResults_.push_back(b1);
373 FitResults_.push_back(a1);
375 if(halo_disc_ > halo_disc_cut)ismipHalo_ =
true;
379 if(debug_)
std::cout<<
"Endof MIP Trail: halo_dic= "<<halo_disc_<<
" nhitcone ="<<NhitCone_<<
" isHalo= "<<ismipHalo_<<std::endl;
380 if(debug_)
std::cout<<
"Endof MIP Trail: Chi2 = "<<chi2<<
" eT= "<<eT<<
" slope = "<<b1<<
" Intercept ="<<a1<<std::endl;
402 if(debug_)
std::cout<<
"Inside GetSeed"<<std::endl;
404 double SeedE = -999.;
414 std::vector< std::pair<DetId, float> > PhotonHit_DetIds = photon->
superCluster()->hitsAndFractions();
417 std::vector< std::pair<DetId, float> >::const_iterator detitr;
418 for(detitr = PhotonHit_DetIds.begin(); detitr != PhotonHit_DetIds.end(); ++detitr)
426 if ( j!= Brechit->end()) thishit = j;
427 if ( j== Brechit->end()){
continue;}
432 double crysE = thishit->energy();
436 seedIeta = detId.
ieta();
437 seedIphi = (detId.
iphi());
440 if(debug_)
std::cout<<
"Current max Seed = "<<SeedE<<
" seedIphi = "<<seedIphi<<
" ieta= "<<seedIeta<<std::endl;
453 if(debug_)
std::cout<<
"End of GetSeed"<<std::endl;
T getParameter(std::string const &) const
reco::SuperClusterRef superCluster() const
Ref to SuperCluster.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
void MIPcalculate(const reco::Photon *, const edm::Event &, const edm::EventSetup &es, reco::Photon::MIPVariables &mipId)
std::vector< EcalRecHit >::const_iterator const_iterator
double haloDiscThreshold_
void GetSeedHighestE(const reco::Photon *photon, const edm::Event &iEvent, const edm::EventSetup &iSetup, edm::Handle< EcalRecHitCollection > Brechit, int &seedIEta, int &seedIPhi, double &seedE)
int iphi() const
get the crystal iphi
edm::EDGetToken EBecalCollection_
std::vector< double > mipFitResults_
int ieta() const
get the crystal ieta
void setup(const edm::ParameterSet &conf, edm::ConsumesCollector &&iC)
std::vector< double > GetMipTrailFit(const reco::Photon *photon, const edm::Event &iEvent, const edm::EventSetup &iSetup, edm::Handle< EcalRecHitCollection > ecalhitsCollEB, double inputRangeY, double inputRangeX, double inputResWidth, double inputHaloDiscCut, int &NhitCone_, bool &ismipHalo_)
double residualWidthEnergy_
T const * product() const
edm::EDGetToken EEecalCollection_