ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/RecoilTools.cc
Revision: 1.1
Committed: Wed Mar 21 18:56:26 2012 UTC (13 years, 1 month ago) by pharris
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_025e, Mit_025d
Log Message:
Adding MET regression

File Contents

# User Rev Content
1 pharris 1.1 #include "MitPhysics/Utils/interface/RecoilTools.h"
2     #include "MitPhysics/Utils/interface/JetTools.h"
3     #include <TFile.h>
4    
5     ClassImp(mithep::RecoilTools)
6    
7     using namespace mithep;
8    
9     RecoilTools::RecoilTools(TString iJetMVAFile) {
10     fJetIDMVA = new JetIDMVA();
11     fJetIDMVA->Initialize( "JetIDMVA",iJetMVAFile,JetIDMVA::kBaseline);
12     }
13     //--------------------------------------------------------------------------------------------------
14     RecoilTools::~RecoilTools() {
15     delete fJetIDMVA;
16     }
17     //--------------------------------------------------------------------------------------------------
18     bool RecoilTools::filter(const PFJet *iJet,Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2) {
19     double pDEta1 = iJet->Eta() - iEta1;
20     double pDPhi1 = fabs(iJet->Phi() - iPhi1); if(pDPhi1 > 2.*TMath::Pi()-pDPhi1) pDPhi1 = 2.*TMath::Pi()-pDPhi1;
21     double pDR1 = sqrt(pDEta1*pDEta1 + pDPhi1*pDPhi1);
22     if(pDR1 < 0.5) return false;
23     double pDEta2 = iJet->Eta() - iEta2;
24     double pDPhi2 = fabs(iJet->Phi() - iPhi2); if(pDPhi2 > 2.*TMath::Pi()-pDPhi2) pDPhi2 = 2.*TMath::Pi()-pDPhi2;
25     double pDR2 = sqrt(pDEta2*pDEta2 + pDPhi2*pDPhi2);
26     if(pDR2 < 0.5) return false;
27     return true;
28     }
29    
30     //--------------------------------------------------------------------------------------------------
31     Met RecoilTools::pfRecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
32     const PFMet *iMet) {
33     Met lPFMet(iMet->Px(),iMet->Py());
34     lPFMet.SetSumEt(iMet->SumEt());
35     lPFMet.SetMex (lPFMet.Mex()+iVisPt*cos(iVisPhi));
36     lPFMet.SetMey (lPFMet.Mey()+iVisPt*sin(iVisPhi));
37     lPFMet.SetSumEt(lPFMet.SumEt()-iVisSumEt);
38     return lPFMet;
39     }
40     //--------------------------------------------------------------------------------------------------
41     Met RecoilTools::trackMet(const PFCandidateCol *iCands,const Vertex *iVertex,Double_t iDZCut) {
42     double trkMetx = 0;
43     double trkMety = 0;
44     double trkSumEt = 0;
45     for(UInt_t i=0; i<iCands->GetEntries(); ++i) {
46     const PFCandidate *pfcand = iCands->At(i);
47     if( (pfcand->HasTrackerTrk() && (fabs(pfcand->TrackerTrk()->DzCorrected(*iVertex))< iDZCut)) ||
48     (pfcand->HasGsfTrk() && (fabs(pfcand->GsfTrk()->DzCorrected(*iVertex)) < iDZCut)) ) {
49     trkMetx -= pfcand->Px();
50     trkMety -= pfcand->Py();
51     trkSumEt += pfcand->Pt();
52     }
53     }
54     Met lMet(trkMetx,trkMety);
55     lMet.SetSumEt(trkSumEt);
56     return lMet;
57     }
58    
59     //--------------------------------------------------------------------------------------------------
60     //Compute the recoil => here this requires the vector sum of the visible components
61     //VisPt => Vector sum pT of the visible non-recoiling components
62     //VisPhi => Vector sum Phi of the visible non-recoiling components
63     //visSumEt => Vector sum Et =f the visible non-recoiling components
64     Met RecoilTools::trackRecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
65     const PFCandidateCol *iCands,const Vertex *iVertex,double iDZCut) {
66     Met lTrkMet = trackMet(iCands,iVertex,iDZCut);
67     lTrkMet.SetMex (lTrkMet.Mex()+iVisPt*cos(iVisPhi));
68     lTrkMet.SetMey (lTrkMet.Mey()+iVisPt*sin(iVisPhi));
69     lTrkMet.SetSumEt(lTrkMet.SumEt()-iVisSumEt);
70     return lTrkMet;
71     }
72     //--------------------------------------------------------------------------------------------------
73     void RecoilTools::addNeut(const PFJet *iJet,FourVectorM &iVec,Double_t &iSumEt,
74     FactorizedJetCorrector *iJetCorrector,const PileupEnergyDensityCol *iPUEnergyDensity,
75     int iSign) {
76     FourVectorM lVec(0,0,0,0);
77     double lPt = fJetIDMVA->correctedPt(iJet,iJetCorrector,iPUEnergyDensity);
78     lPt *= (iJet->NeutralEmEnergy()/iJet->E() + iJet->NeutralHadronEnergy()/iJet->E());
79     lVec.SetPt(lPt); lVec.SetEta(iJet->Eta()); lVec.SetPhi(iJet->Phi()); lVec.SetM(iJet->Mass());
80     if(iSign > 0) iVec -= lVec;
81     if(iSign < 0) iVec += lVec;
82     iSumEt += lPt;
83     }
84     //--------------------------------------------------------------------------------------------------
85     Met RecoilTools::NoPUMet( const PFJetCol *iJets,FactorizedJetCorrector *iJetCorrector,
86     const PileupEnergyDensityCol *iPileupEnergyDensity,
87     const PFCandidateCol *iCands,const Vertex *iVertex,
88     Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,Double_t iDZCut) {
89    
90     FourVectorM lVec (0,0,0,0); double lSumEt = 0;
91     for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) {
92     const PFCandidate *pPF = iCands->At(i0);
93     const Track* pTrack = pPF->TrackerTrk();
94     if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
95     if(pTrack == 0 ) continue;
96     if( !((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) ||
97     (pPF->HasGsfTrk() && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex)) <iDZCut)))) continue;
98     lVec -= pPF->Mom();
99     lSumEt += pPF->Pt();
100     }
101     for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
102     const PFJet *pJet = iJets->At(i0);
103     if(!fJetIDMVA->pass(pJet,iVertex,iJetCorrector,iPileupEnergyDensity)) continue;
104     if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2)) continue; //Quick cleaning==> if not done already
105     addNeut(pJet,lVec,lSumEt,iJetCorrector,iPileupEnergyDensity);
106     }
107     Met lMet(lVec.Px(),lVec.Py());
108     lMet.SetSumEt(lSumEt);
109     return lMet;
110     }
111     //--------------------------------------------------------------------------------------------------
112     Met RecoilTools::NoPURecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
113     const PFJetCol *iJets,FactorizedJetCorrector *iJetCorrector,
114     const PileupEnergyDensityCol *iPileupEnergyDensity,
115     const PFCandidateCol *iCands,const Vertex *iVertex,
116     Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
117     Double_t iDZCut) {
118    
119     Met lNoPUMet = NoPUMet(iJets,iJetCorrector,iPileupEnergyDensity,iCands,iVertex,iPhi1,iEta1,iPhi2,iEta2,iDZCut);
120     lNoPUMet.SetMex (lNoPUMet.Mex()+iVisPt*cos(iVisPhi));
121     lNoPUMet.SetMey (lNoPUMet.Mey()+iVisPt*sin(iVisPhi));
122     lNoPUMet.SetSumEt(lNoPUMet.SumEt()-iVisSumEt);
123     return lNoPUMet;
124     }
125     //--------------------------------------------------------------------------------------------------
126     Met RecoilTools::PUCMet( const PFJetCol *iJets,FactorizedJetCorrector *iJetCorrector,
127     const PileupEnergyDensityCol *iPileupEnergyDensity,
128     const PFCandidateCol *iCands,const Vertex *iVertex,
129     Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
130     Double_t iDZCut) {
131    
132     FourVectorM lVec (0,0,0,0); double lSumEt = 0;
133     for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) {
134     const PFCandidate *pPF = iCands->At(i0);
135     const Track* pTrack = pPF->TrackerTrk();
136     if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
137     if(pTrack == 0 &&
138     (pPF->PFType() == PFCandidate::eGamma ||
139     pPF->PFType() == PFCandidate::eEGammaHF ||
140     pPF->PFType() == PFCandidate::eNeutralHadron ||
141     pPF->PFType() == PFCandidate::eHadronHF ))
142     {lVec -= pPF->Mom(); lSumEt += pPF->Pt();}
143     if(pTrack == 0 ) continue;
144     if( !((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) ||
145     (pPF->HasGsfTrk() && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex)) <iDZCut)))) continue;
146     lVec -= pPF->Mom();
147     lSumEt += pPF->Pt();
148     }
149     for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
150     const PFJet *pJet = iJets->At(i0);
151     if(fJetIDMVA->correctedPt(pJet,iJetCorrector,iPileupEnergyDensity) < fJetIDMVA->fJetPtMin
152     && pJet->TrackCountingHighEffBJetTagsDisc() == -100) continue; //This line is a bug in the Met training//
153     if(!JetTools::passPFLooseId(pJet)) continue;
154     if(fJetIDMVA->pass(pJet,iVertex,iJetCorrector,iPileupEnergyDensity)) continue;
155     if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2)) continue; //Quick cleaning==> if not done already
156     addNeut(pJet,lVec,lSumEt,iJetCorrector,iPileupEnergyDensity,-1);
157     }
158     Met lMet(lVec.Px(),lVec.Py());
159     lMet.SetSumEt(lSumEt);
160     return lMet;
161     }
162     //----> This MET is a bug need to fix it
163     //--------------------------------------------------------------------------------------------------
164     Met RecoilTools::PUCRecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
165     const PFJetCol *iJets,FactorizedJetCorrector *iJetCorrector,
166     const PileupEnergyDensityCol *iPileupEnergyDensity,
167     const PFCandidateCol *iCands,const Vertex *iVertex,
168     Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
169     Double_t iDZCut) {
170     Met lPUCMet = PUCMet(iJets,iJetCorrector,iPileupEnergyDensity,iCands,iVertex,iPhi1,iEta1,iPhi2,iEta2,iDZCut);
171     lPUCMet.SetMex (lPUCMet.Mex()+iVisPt*cos(iVisPhi));
172     lPUCMet.SetMey (lPUCMet.Mey()+iVisPt*sin(iVisPhi));
173     lPUCMet.SetSumEt(lPUCMet.SumEt()-iVisSumEt);
174     return lPUCMet;
175     }
176     //--------------------------------------------------------------------------------------------------
177     Met RecoilTools::PUMet( const PFJetCol *iJets,FactorizedJetCorrector *iJetCorrector,
178     const PileupEnergyDensityCol *iPileupEnergyDensity,
179     const PFCandidateCol *iCands,const Vertex *iVertex,
180     Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
181     Double_t iDZCut) {
182    
183     FourVectorM lVec (0,0,0,0); double lSumEt = 0;
184     for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) {
185     const PFCandidate *pPF = iCands->At(i0);
186     const Track* pTrack = pPF->TrackerTrk();
187     if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
188     if(pTrack == 0 ) continue;
189     if( ((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) ||
190     (pPF->HasGsfTrk() && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex)) <iDZCut)))) continue;
191     lVec -= pPF->Mom();
192     lSumEt += pPF->Pt();
193     }
194     for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
195     const PFJet *pJet = iJets->At(i0);
196     if(fJetIDMVA->correctedPt(pJet,iJetCorrector,iPileupEnergyDensity) < fJetIDMVA->fJetPtMin
197     && pJet->TrackCountingHighEffBJetTagsDisc() == -100) continue; //This line is a bug in the Met training//
198     if(!JetTools::passPFLooseId(pJet)) continue;
199     if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2)) continue; //Quick cleaning
200     if(fJetIDMVA->pass(pJet,iVertex,iJetCorrector,iPileupEnergyDensity)) continue;
201     addNeut(pJet,lVec,lSumEt,iJetCorrector,iPileupEnergyDensity);
202     }
203     Met lMet(lVec.Px(),lVec.Py());
204     lMet.SetSumEt(lSumEt);
205     return lMet;
206     }
207     //--------------------------------------------------------------------------------------------------
208     Met RecoilTools::PURecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
209     const PFJetCol *iJets,FactorizedJetCorrector *iJetCorrector,
210     const PileupEnergyDensityCol *iPileupEnergyDensity,
211     const PFCandidateCol *iCands,const Vertex *iVertex,
212     Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
213     Double_t iDZCut) {
214     Met lPUMet = PUMet(iJets,iJetCorrector,iPileupEnergyDensity,iCands,iVertex,iPhi1,iEta1,iPhi2,iEta2,iDZCut);
215     lPUMet.SetMex (lPUMet.Mex()+iVisPt*cos(iVisPhi));
216     lPUMet.SetMey (lPUMet.Mey()+iVisPt*sin(iVisPhi));
217     lPUMet.SetSumEt(lPUMet.SumEt()-iVisSumEt);
218     return lPUMet;
219     }