/ElectroWeakAnalysis/ZEE/interface/SCEnergyCorrections.h

  1#include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
  2
  3reco::CaloClusterPtrVector CaloClusterVectorCopier(const reco::SuperCluster& sc);
  4
  5reco::SuperCluster fEtaScCorr(const reco::SuperCluster& sc)
  6{
  7	reco::CaloClusterPtrVector bcs = CaloClusterVectorCopier(sc);		
  8	double ieta = fabs(sc.eta())*(5/0.087);
  9	double p0 = 40.2198;
 10	double p1 = -3.03103e-6;
 11	double newE;
 12//                      std::cout << "Corrected E = Raw E * (1+ p1*(ieta - p0)*(ieta - p0))"<< std::endl;
 13	if ( ieta < p0 ) newE = sc.rawEnergy();
 14	else newE = sc.rawEnergy()/(1 + p1*(ieta - p0)*(ieta - p0));
 15
 16	reco::SuperCluster corrSc(newE, sc.position(), sc.seed(), bcs, sc.preshowerEnergy(), 0., 0.);
 17	return corrSc;
 18}
 19
 20reco::SuperCluster fBremScCorr(const reco::SuperCluster& sc, const edm::ParameterSet& ps)
 21{
 22	std::vector<double> fBrem = ps.getParameter<std::vector<double> >("fBremVec");
 23	double bremFrLowThr = ps.getParameter<double>("brLinearLowThr");
 24	double bremFrHighThr = ps.getParameter<double>("brLinearHighThr");
 25	
 26	reco::CaloClusterPtrVector bcs = CaloClusterVectorCopier(sc);		
 27	double bremFrac = sc.phiWidth()/sc.etaWidth();	
 28	double newE = sc.energy();
 29	if(fabs(sc.eta()) < 1.479)
 30	{
 31		reco::SuperCluster fEtaSC = fEtaScCorr(sc);
 32		reco::CaloClusterPtrVector bcs = CaloClusterVectorCopier(sc);		
 33		newE = fEtaSC.energy();
 34	}
 35
 36	if(bremFrac < bremFrLowThr)  bremFrac = bremFrLowThr;
 37	if(bremFrac < bremFrHighThr)  bremFrac = bremFrHighThr;
 38	
 39	double p0 = fBrem[0]; 
 40	double p1 = fBrem[1]; 
 41	double p2 = fBrem[2]; 
 42	double p3 = fBrem[3]; 
 43	double p4 = fBrem[4]; 
 44	//
 45	double threshold = p4; 
 46	  
 47	double y = p0*threshold*threshold + p1*threshold + p2;
 48	double yprime = 2*p0*threshold + p1;
 49	double a = p3;
 50	double b = yprime - 2*a*threshold;
 51	double c = y - a*threshold*threshold - b*threshold;
 52	 
 53	double fCorr = 1;
 54	if( bremFrac < threshold ) 
 55	    fCorr = p0*bremFrac*bremFrac + p1*bremFrac + p2;
 56	else 
 57	    fCorr = a*bremFrac*bremFrac + b*bremFrac + c;
 58
 59	newE /= fCorr;
 60	reco::SuperCluster corrSc(newE, sc.position(), sc.seed(), bcs, sc.preshowerEnergy(), 0., 0.);
 61	return corrSc;
 62}
 63
 64reco::SuperCluster fEtEtaCorr(const reco::SuperCluster& sc, const edm::ParameterSet& ps)
 65{
 66  // et -- Et of the SuperCluster (with respect to (0,0,0))
 67  // eta -- eta of the SuperCluster
 68  	std::vector<double> fEtEtaParams = ps.getParameter<std::vector<double> >("fEtEtaParamsVec");
 69
 70	reco::SuperCluster fBremSC = fBremScCorr(sc, ps);
 71	reco::CaloClusterPtrVector bcs = CaloClusterVectorCopier(sc);		
 72
 73	double eta = sc.eta();
 74	double et = fBremSC.energy()/cosh(eta);
 75	double fCorr = 0.;
 76
 77	double p0 = fEtEtaParams[0] + fEtEtaParams[1]/(et + fEtEtaParams[ 2]) + fEtEtaParams[ 3]/(et*et);
 78	double p1 = fEtEtaParams[4] + fEtEtaParams[5]/(et + fEtEtaParams[ 6]) + fEtEtaParams[ 7]/(et*et);
 79	double p2 = fEtEtaParams[8] + fEtEtaParams[9]/(et + fEtEtaParams[10]) + fEtEtaParams[11]/(et*et);
 80
 81	fCorr = 
 82	p0 + 
 83	p1 * atan(fEtEtaParams[12]*(fEtEtaParams[13]-fabs(eta))) + fEtEtaParams[14] * fabs(eta) + 
 84	p1 * fEtEtaParams[15] * fabs(eta) +
 85	p2 * fEtEtaParams[16] * eta * eta; 
 86
 87	if ( fCorr < 0.5 ) fCorr = 0.5;
 88
 89	double newE = et/(fCorr*cosh(eta));
 90	reco::SuperCluster corrSc(newE, sc.position(), sc.seed(), bcs, sc.preshowerEnergy(), 0., 0.);
 91	return corrSc;
 92}
 93
 94reco::SuperCluster fEAddScCorr(const reco::SuperCluster& sc, double Ecorr)
 95{
 96	reco::CaloClusterPtrVector bcs = CaloClusterVectorCopier(sc);		
 97	
 98	double newE = sc.rawEnergy()+Ecorr; 
 99	reco::SuperCluster corrSc(newE, sc.position(), sc.seed(), bcs, sc.preshowerEnergy(), 0., 0.);
100	return corrSc;
101}
102	
103
104reco::CaloClusterPtrVector CaloClusterVectorCopier(const reco::SuperCluster& sc)
105{
106  	reco::CaloClusterPtrVector clusters_v;
107
108  	for(reco::CaloCluster_iterator cluster = sc.clustersBegin(); cluster != sc.clustersEnd(); cluster ++)
109 	{
110		clusters_v.push_back(*cluster);
111	}
112
113	return clusters_v;
114}