00001
00002
00003
00004
00005 #include "JetMETCorrections/Type1MET/interface/TauMETAlgo.h"
00006 #include "DataFormats/Math/interface/deltaR.h"
00007
00008 using namespace std;
00009 using namespace reco;
00010
00011 typedef math::XYZTLorentzVector LorentzVector;
00012 typedef math::XYZPoint Point;
00013
00014 TauMETAlgo::TauMETAlgo() {}
00015 TauMETAlgo::~TauMETAlgo() {}
00016
00017 void TauMETAlgo::run(edm::Event& iEvent,const edm::EventSetup& iSetup,
00018 edm::Handle<PFTauCollection> tauHandle,edm::Handle<CaloJetCollection> calojetHandle,
00019 double jetPTthreshold,double jetEMfracLimit,
00020 const JetCorrector& correctedjets,const std::vector<CaloMET>& uncorMET,double jetMatchDeltaR,
00021 double tauMinEt,double tauEtaMax,bool useSeedTrack,double seedTrackPt,bool useTrackIsolation,double trackIsolationMinPt,
00022 bool useECALIsolation,double gammaIsolationMinPt,bool useProngStructure,std::vector<CaloMET>* corrMET) {
00023
00024
00025
00026 double DeltaPx = 0.0;
00027 double DeltaPy = 0.0;
00028 double DeltaSumET = 0.0;
00029 for(PFTauCollection::size_type iPFTau=0;iPFTau<tauHandle->size();iPFTau++) {
00030 PFTauRef thePFTau(tauHandle,iPFTau);
00031 bool matchFlag = false;
00032 bool goodTau = false;
00033 if((fabs(thePFTau->eta()) <= tauEtaMax) && (thePFTau->et() >= tauMinEt)) {
00034 goodTau = true;
00035 if(useSeedTrack) {
00036 if(!(thePFTau->leadPFChargedHadrCand())) {goodTau = false;}
00037 else {
00038 if(thePFTau->leadPFChargedHadrCand()->et() < seedTrackPt) {goodTau = false;}
00039 }
00040 }
00041 if(useTrackIsolation) {
00042 PFCandidateRefVector PFTauProdIsolCHCands = (*thePFTau).isolationPFChargedHadrCands();
00043 for(PFCandidateRefVector::const_iterator iChargedHadrCand=PFTauProdIsolCHCands.begin();
00044 iChargedHadrCand!=PFTauProdIsolCHCands.end();++iChargedHadrCand) {
00045 if((**iChargedHadrCand).pt() > trackIsolationMinPt) {goodTau = false;}
00046 }
00047 }
00048 if(useECALIsolation) {
00049 PFCandidateRefVector PFTauProdIsolGammaCands = (*thePFTau).isolationPFGammaCands();
00050 for(PFCandidateRefVector::const_iterator iGammaCand=PFTauProdIsolGammaCands.begin();
00051 iGammaCand!=PFTauProdIsolGammaCands.end();++iGammaCand) {
00052 if((**iGammaCand).et() > gammaIsolationMinPt) {goodTau = false;}
00053 }
00054 }
00055 if(useProngStructure) {
00056 if((thePFTau->signalPFChargedHadrCands().size() != 1) && (thePFTau->signalPFChargedHadrCands().size() != 3)) {goodTau = false;}
00057 }
00058 }
00059 if(goodTau) {
00060 for(CaloJetCollection::const_iterator calojetIter = calojetHandle->begin();calojetIter != calojetHandle->end();++calojetIter) {
00061 if(deltaR(calojetIter->p4(),thePFTau->p4())<jetMatchDeltaR) {
00062 if( calojetIter->pt() > jetPTthreshold && calojetIter->emEnergyFraction() < jetEMfracLimit ) {
00063 double correct = correctedjets.correction (*calojetIter);
00064 DeltaPx += ((correct * calojetIter->px()) - thePFTau->px());
00065 DeltaPy += ((correct * calojetIter->py()) - thePFTau->py());
00066 DeltaSumET += ((correct * calojetIter->et()) - thePFTau->et());
00067 } else {
00068 DeltaPx += (calojetIter->px() - thePFTau->px());
00069 DeltaPy += (calojetIter->py() - thePFTau->py());
00070 DeltaSumET += (calojetIter->et() - thePFTau->et());
00071 }
00072 if(matchFlag) {std::cerr << "### TauMETAlgo - ERROR: Multiple jet matches!!!! " << std::endl;}
00073 matchFlag = true;
00074 }
00075 }
00076 }
00077 }
00078 CorrMETData delta;
00079 delta.mex = DeltaPx;
00080 delta.mey = DeltaPy;
00081 delta.sumet = - DeltaSumET;
00082 const CaloMET* u = &(uncorMET.front());
00083 double corrMetPx = u->px() + delta.mex;
00084 double corrMetPy = u->py() + delta.mey;
00085 MET::LorentzVector correctedMET4vector( corrMetPx, corrMetPy, 0., sqrt(corrMetPx*corrMetPx + corrMetPy*corrMetPy));
00086 std::vector<CorrMETData> corrections = u->mEtCorr();
00087 corrections.push_back(delta);
00088 CaloMET result;
00089 result = CaloMET(u->getSpecific(), (u->sumEt() + delta.sumet), corrections, correctedMET4vector, u->vertex());
00090 corrMET->push_back(result);
00091 }
00092
00093 void TauMETAlgo::run(edm::Event& iEvent,const edm::EventSetup& iSetup,
00094 edm::Handle<PFTauCollection> tauHandle,edm::Handle<CaloJetCollection> calojetHandle,
00095 double jetPTthreshold,double jetEMfracLimit,
00096 const JetCorrector& correctedjets,const std::vector<MET>& uncorMET,double jetMatchDeltaR,
00097 double tauMinEt,double tauEtaMax,bool useSeedTrack,double seedTrackPt,bool useTrackIsolation,double trackIsolationMinPt,
00098 bool useECALIsolation,double gammaIsolationMinPt,bool useProngStructure,std::vector<MET>* corrMET) {
00099
00100 std::cerr << "TauMETAlgo::run -> Test.. " << std::endl;
00101
00102 double DeltaPx = 0.0;
00103 double DeltaPy = 0.0;
00104 double DeltaSumET = 0.0;
00105 for(PFTauCollection::size_type iPFTau=0;iPFTau<tauHandle->size();iPFTau++) {
00106 PFTauRef thePFTau(tauHandle,iPFTau);
00107 bool matchFlag = false;
00108 bool goodTau = false;
00109 if((fabs(thePFTau->eta()) <= tauEtaMax) && (thePFTau->et() >= tauMinEt)) {
00110 goodTau = true;
00111 if(useSeedTrack) {
00112 if(!(thePFTau->leadPFChargedHadrCand())) {goodTau = false;}
00113 else {
00114 if(thePFTau->leadPFChargedHadrCand()->et() < seedTrackPt) {goodTau = false;}
00115 }
00116 }
00117 if(useTrackIsolation) {
00118 PFCandidateRefVector PFTauProdIsolCHCands = (*thePFTau).isolationPFChargedHadrCands();
00119 for(PFCandidateRefVector::const_iterator iChargedHadrCand=PFTauProdIsolCHCands.begin();
00120 iChargedHadrCand!=PFTauProdIsolCHCands.end();++iChargedHadrCand) {
00121 if((**iChargedHadrCand).pt() > trackIsolationMinPt) {goodTau = false;}
00122 }
00123 }
00124 if(useECALIsolation) {
00125 PFCandidateRefVector PFTauProdIsolGammaCands = (*thePFTau).isolationPFGammaCands();
00126 for(PFCandidateRefVector::const_iterator iGammaCand=PFTauProdIsolGammaCands.begin();
00127 iGammaCand!=PFTauProdIsolGammaCands.end();++iGammaCand) {
00128 if((**iGammaCand).et() > gammaIsolationMinPt) {goodTau = false;}
00129 }
00130 }
00131 if(useProngStructure) {
00132 if((thePFTau->signalPFChargedHadrCands().size() != 1) && (thePFTau->signalPFChargedHadrCands().size() != 3)) {goodTau = false;}
00133 }
00134 }
00135 if(goodTau) {
00136 for(CaloJetCollection::const_iterator calojetIter = calojetHandle->begin();calojetIter != calojetHandle->end();++calojetIter) {
00137 if(deltaR(calojetIter->p4(),thePFTau->p4())<jetMatchDeltaR) {
00138 if( calojetIter->pt() > jetPTthreshold && calojetIter->emEnergyFraction() < jetEMfracLimit ) {
00139 double correct = correctedjets.correction (*calojetIter);
00140 DeltaPx += ((correct * calojetIter->px()) - thePFTau->px());
00141 DeltaPy += ((correct * calojetIter->py()) - thePFTau->py());
00142 DeltaSumET += ((correct * calojetIter->et()) - thePFTau->et());
00143 } else {
00144 DeltaPx += (calojetIter->px() - thePFTau->px());
00145 DeltaPy += (calojetIter->py() - thePFTau->py());
00146 DeltaSumET += (calojetIter->et() - thePFTau->et());
00147 }
00148 if(matchFlag) {std::cerr << "### TauMETAlgo - ERROR: Multiple jet matches!!!! " << std::endl;}
00149 matchFlag = true;
00150 }
00151 }
00152 }
00153 }
00154 CorrMETData delta;
00155 delta.mex = DeltaPx;
00156 delta.mey = DeltaPy;
00157 delta.sumet = - DeltaSumET;
00158 const MET* u = &(uncorMET.front());
00159 double corrMetPx = u->px() + delta.mex;
00160 double corrMetPy = u->py() + delta.mey;
00161 MET::LorentzVector correctedMET4vector( corrMetPx, corrMetPy, 0., sqrt(corrMetPx*corrMetPx + corrMetPy*corrMetPy));
00162 std::vector<CorrMETData> corrections = u->mEtCorr();
00163 corrections.push_back(delta);
00164 MET result;
00165 result = MET((u->sumEt() + delta.sumet), corrections, correctedMET4vector, u->vertex());
00166 corrMET->push_back(result);
00167 }
00168
1.7.3