CMS 3D CMS Logo

CMSSW_4_4_3_patch1/src/RecoEgamma/EgammaHFProducers/src/HFRecoEcalCandidateAlgo.cc

Go to the documentation of this file.
00001 
00009 //#includes
00010 #include "DataFormats/Math/interface/LorentzVector.h"
00011 #include "DataFormats/Math/interface/deltaR.h"
00012 #include "DataFormats/Math/interface/deltaPhi.h"
00013 #include "CLHEP/Vector/LorentzVector.h"
00014 #include <iostream>
00015 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
00016 #include "RecoEgamma/EgammaHFProducers/interface/HFRecoEcalCandidateAlgo.h"
00017 using namespace std;
00018 using namespace reco;
00019 
00020 HFRecoEcalCandidateAlgo::HFRecoEcalCandidateAlgo(bool correct,double e9e25Cut,double intercept2DCut,
00021                                                  const std::vector<double>& e1e9Cut,
00022                                                  const std::vector<double>& eCOREe9Cut,
00023                                                  const std::vector<double>& eSeLCut) :
00024 
00025   m_correct(correct), 
00026   m_e9e25Cut(e9e25Cut),
00027   m_intercept2DCut(intercept2DCut),
00028   m_e1e9Cuthi(e1e9Cut[1]),
00029   m_eCOREe9Cuthi(eCOREe9Cut[1]),
00030   m_eSeLCuthi(eSeLCut[1]),
00031   m_e1e9Cutlo(e1e9Cut[0]),
00032   m_eCOREe9Cutlo(eCOREe9Cut[0]),
00033   m_eSeLCutlo(eSeLCut[0]){
00034 
00035 }
00036 
00037 RecoEcalCandidate HFRecoEcalCandidateAlgo::correctEPosition(const SuperCluster& original , const HFEMClusterShape& shape) {
00038   double energyCorrect=0.7397;//.7515;
00039   double etaCorrect=.00938422+0.00682824*sin(6.28318531*shape.CellEta());//.0144225-.00484597*sin(6.17851*shape.CellEta());//0.01139;
00040   double phiAmpCorrect=0.00644091;//-0.006483;
00041   double phiFreqCorrect=6.28318531;//6.45377;
00042 
00043   double corEnergy= original.energy()/energyCorrect;
00044   double corEta=original.eta();
00045   corEta+=(original.eta()>0)?(etaCorrect):(-etaCorrect);
00046   double corPhi=original.phi()+phiAmpCorrect*sin(phiFreqCorrect*shape.CellPhi());
00047   double corPx=corEnergy*cos(corPhi)/cosh(corEta);
00048   double corPy=corEnergy*sin(corPhi)/cosh(corEta);
00049   double corPz=corEnergy*tanh(corEta);
00050     RecoEcalCandidate corCand(0,
00051                               math::XYZTLorentzVector(corPx,corPy,corPz,corEnergy),
00052                               math::XYZPoint(0,0,0));
00053 
00054  return corCand;
00055 }
00056 
00057 void HFRecoEcalCandidateAlgo::produce(const edm::Handle<SuperClusterCollection>& SuperClusters,
00058                                       const HFEMClusterShapeAssociationCollection& AssocShapes,
00059                                       RecoEcalCandidateCollection& RecoECand) {
00060   
00061   
00062   
00063   //get super's and cluster shapes and associations 
00064   for (unsigned int i=0; i < SuperClusters->size(); ++i) {
00065     const SuperCluster& supClus=(*SuperClusters)[i];    
00066     reco::SuperClusterRef theClusRef=edm::Ref<SuperClusterCollection>(SuperClusters,i);
00067   const HFEMClusterShapeRef clusShapeRef=AssocShapes.find(theClusRef)->val;
00068     const HFEMClusterShape& clusShape=*clusShapeRef;
00069 
00070     // basic candidate
00071     double px=supClus.energy()*cos(supClus.phi())/cosh(supClus.eta());
00072     double py=supClus.energy()*sin(supClus.phi())/cosh(supClus.eta());
00073     double pz=supClus.energy()*tanh(supClus.eta());
00074     RecoEcalCandidate theCand(0,
00075                               math::XYZTLorentzVector(px,py,pz,supClus.energy()),
00076                               math::XYZPoint(0,0,0));
00077 
00078     // correct it?
00079     if (m_correct)
00080       theCand=correctEPosition(supClus,clusShape);
00081 
00082     double e9e25=clusShape.eLong3x3()/clusShape.eLong5x5();
00083     double e1e9=clusShape.eLong1x1()/clusShape.eLong3x3();
00084     // EMID cuts...  
00085     //if((clusShape.e9e25()> m_e9e25Cut)&&((clusShape.eCOREe9()-(clusShape.eSeL()*1.125)) > m_intercept2DCut)){
00086     //  if((e9e25> m_e9e25Cut)&&((clusShape.eCOREe9()-(clusShape.eSeL()*1.125)) > m_intercept2DCut)){
00087     double var2d=(clusShape.eCOREe9()-(clusShape.eSeL()*1.125));
00088  
00089     bool isAcceptable=true;
00090     isAcceptable=isAcceptable && (e9e25> m_e9e25Cut);
00091     isAcceptable=isAcceptable && (var2d > m_intercept2DCut);
00092     isAcceptable=isAcceptable && ((e1e9< m_e1e9Cuthi)&&(e1e9> m_e1e9Cutlo));
00093     isAcceptable=isAcceptable && ((clusShape.eCOREe9()< m_eCOREe9Cuthi)&&(clusShape.eCOREe9()>  m_eCOREe9Cutlo));
00094     isAcceptable=isAcceptable && ((clusShape.eSeL()<m_eSeLCuthi)&&(clusShape.eSeL()>  m_eSeLCutlo));
00095    
00096     
00097     if(isAcceptable){
00098 
00099       theCand.setSuperCluster(theClusRef);
00100       RecoECand.push_back(theCand);
00101     }
00102   }
00103 }
00104  
00105 
00106 //&&((clusShape.e1e9()< m_e1e9Cuthi)&&(clusShape.e1e9()> m_e1e9Cutlo))&&((clusShape.eCOREe9()< m_eCOREe9Cuthi)&&(clusShape.eCOREe9()>  m_eCOREe9Cutlo))&&((clusShape.eSeL()<m_eSeLCuthi)&&(clusShape.eSeL()>  m_eSeLCutlo))