ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/RecoilTools.cc
(Generate patch)

Comparing UserCode/MitPhysics/Utils/src/RecoilTools.cc (file contents):
Revision 1.13 by pharris, Tue May 1 18:17:54 2012 UTC vs.
Revision 1.17 by pharris, Fri May 25 14:10:22 2012 UTC

# Line 26 | Line 26 | bool RecoilTools::filter(const PFJet *iJ
26    if(pDR2 < 0.5) return false;
27    return true;
28   }
29 <
29 > //--------------------------------------------------------------------------------------------------
30 > bool RecoilTools::filter(const PFCandidate *iCand,Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2) {
31 >  double pDEta1 = iCand->Eta() - iEta1;
32 >  double pDPhi1 = fabs(iCand->Phi() - iPhi1); if(pDPhi1 > 2.*TMath::Pi()-pDPhi1) pDPhi1 = 2.*TMath::Pi()-pDPhi1;
33 >  double pDR1   = sqrt(pDEta1*pDEta1 + pDPhi1*pDPhi1);
34 >  if(pDR1 < 0.3) return false;
35 >  double pDEta2 = iCand->Eta() - iEta2;
36 >  double pDPhi2 = fabs(iCand->Phi() - iPhi2); if(pDPhi2 > 2.*TMath::Pi()-pDPhi2) pDPhi2 = 2.*TMath::Pi()-pDPhi2;
37 >  double pDR2   = sqrt(pDEta2*pDEta2 + pDPhi2*pDPhi2);
38 >  if(pDR2 < 0.3) return false;
39 >  return true;
40 > }
41   //--------------------------------------------------------------------------------------------------
42   Met RecoilTools::pfRecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
43                            const PFCandidateCol *iCands) {
44    double lSumEt = 0;
45    FourVectorM lVec(0,0,0,0);
46 <  for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) { lVec -= iCands->At(i0)->Mom(); lSumEt += iCands->At(i0)->Et();}
46 >  for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) {
47 >    lVec -= iCands->At(i0)->Mom();
48 >    lSumEt += iCands->At(i0)->Et();
49 >  }
50    Met lPFMet(lVec.Px(),lVec.Py());
51    lPFMet.SetSumEt(lSumEt);
52    lPFMet.SetMex  (lPFMet.Mex()+iVisPt*cos(iVisPhi));  
# Line 41 | Line 55 | Met RecoilTools::pfRecoil(Double_t iVisP
55    return lPFMet;
56   }
57   //--------------------------------------------------------------------------------------------------
58 + Met RecoilTools::pfRecoil(double iPhi1,double iEta1,double iPhi2,double iEta2,
59 +                          const PFCandidateCol *iCands) {
60 +  double lSumEt = 0;
61 +  FourVectorM lVec(0,0,0,0);
62 +  for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) {
63 +    const PFCandidate *pfcand = iCands->At(i0);
64 +    if(!filter(pfcand,iPhi1,iEta1,iPhi2,iEta2)) continue;
65 +    lVec   -= pfcand->Mom();
66 +    lSumEt += pfcand->Et();
67 +  }
68 +  Met lPFMet(lVec.Px(),lVec.Py());
69 +  lPFMet.SetSumEt(lSumEt);
70 +  return lPFMet;
71 + }
72 + //--------------------------------------------------------------------------------------------------
73   Met RecoilTools::trackMet(const PFCandidateCol *iCands,const Vertex *iVertex,Double_t iDZCut) {
74    double trkMetx  = 0;
75    double trkMety  = 0;
# Line 58 | Line 87 | Met RecoilTools::trackMet(const PFCandid
87    lMet.SetSumEt(trkSumEt);
88    return lMet;
89   }
90 <
90 > //--------------------------------------------------------------------------------------------------
91 > Met RecoilTools::trackRecoil(double iPhi1,double iEta1,double iPhi2,double iEta2,const PFCandidateCol *iCands,const Vertex *iVertex,Double_t iDZCut) {
92 >  double trkMetx  = 0;
93 >  double trkMety  = 0;
94 >  double trkSumEt = 0;
95 >  for(UInt_t i0=0; i0<iCands->GetEntries(); i0++) {
96 >    const PFCandidate *pfcand = iCands->At(i0);
97 >    if(!filter(pfcand,iPhi1,iEta1,iPhi2,iEta2)) continue;
98 >    if( (pfcand->HasTrackerTrk() && (fabs(pfcand->TrackerTrk()->DzCorrected(*iVertex))< iDZCut)) ||
99 >        (pfcand->HasGsfTrk()     && (fabs(pfcand->GsfTrk()->DzCorrected(*iVertex))    < iDZCut)) ) {
100 >      trkMetx  -= pfcand->Px();
101 >      trkMety  -= pfcand->Py();
102 >      trkSumEt += pfcand->Pt();
103 >    }
104 >  }
105 >  Met lMet(trkMetx,trkMety);
106 >  lMet.SetSumEt(trkSumEt);
107 >  return lMet;
108 > }
109   //--------------------------------------------------------------------------------------------------
110   //Compute the recoil => here this requires the vector sum of the visible components
111   //VisPt    => Vector sum pT  of the visible non-recoiling components
# Line 78 | Line 125 | void RecoilTools::addNeut(const PFJet *i
125                            int iSign) {
126    FourVectorM lVec(0,0,0,0);
127    double lPt = fJetIDMVA->correctedPt(iJet,iJetCorrector,iPUEnergyDensity);
128 <  if(iJet->RawMom().Pt() < 10) lPt = TMath::Max(iJet->RawMom().Pt()-iJet->JetArea()*iPUEnergyDensity->At(0)->Rho(),0.);
128 >  //if(iJet->RawMom().Pt() < 10) lPt = TMath::Max(iJet->RawMom().Pt()-iJet->JetArea()*iPUEnergyDensity->At(0)->Rho(),0.);
129    lPt *= (iJet->NeutralEmEnergy()/iJet->RawMom().E() + iJet->NeutralHadronEnergy()/iJet->RawMom().E());
130    lVec.SetPt(lPt); lVec.SetEta(iJet->Eta()); lVec.SetPhi(iJet->Phi()); lVec.SetM(iJet->Mass());
131    if(iSign > 0) iVec -= lVec;
132    if(iSign < 0) iVec += lVec;
133 <  iSumEt += lPt;
133 >  //iSumEt += lPt;
134    //=== Above was a bug in the training
135 <  //if(iSign > 0) iSumEt += lPt;
136 <  //if(iSign < 0) iSumEt -= lPt;
135 >  if(iSign > 0) iSumEt += lPt;
136 >  if(iSign < 0) iSumEt -= lPt;
137   }
138  
139   //--------------------------------------------------------------------------------------------------
# Line 94 | Line 141 | void RecoilTools::addNeut(const PFJet *i
141   void RecoilTools::addNeut(const PFJet *iJet,FourVectorM &iVec,Double_t &iSumEt,double iRho,int iSign) {
142    FourVectorM lVec(0,0,0,0);
143    double lPt = iJet->Pt();
144 <  if(iJet->RawMom().Pt() < 10) lPt = TMath::Max(iJet->RawMom().Pt()-iJet->JetArea()*iRho,0.);//to be fixed
144 >  //if(iJet->RawMom().Pt() < 10) lPt = TMath::Max(iJet->RawMom().Pt()-iJet->JetArea()*iRho,0.);//to be fixed
145    //if(iJet->RawMom().Pt() < 10) lPt = iJet->RawMom().Pt()*iJet->L1OffsetCorrectionScale();
146    lPt *= (iJet->NeutralEmEnergy()/iJet->RawMom().E() + iJet->NeutralHadronEnergy()/iJet->RawMom().E());
147    lVec.SetPt(lPt); lVec.SetEta(iJet->Eta()); lVec.SetPhi(iJet->Phi()); lVec.SetM(iJet->Mass());
148    if(iSign > 0) iVec   -= lVec;
149    if(iSign < 0) iVec   += lVec;
150 <  iSumEt += lPt;
150 >  //iSumEt += lPt;
151    //=== Above was a bug in the training
152 <  //if(iSign > 0) iSumEt += lPt;
153 <  //if(iSign < 0) iSumEt -= lPt;
152 >  if(iSign > 0) iSumEt += lPt;
153 >  if(iSign < 0) iSumEt -= lPt;
154   }
155  
156   //--------------------------------------------------------------------------------------------------
# Line 118 | Line 165 | Met RecoilTools::NoPUMet( const PFJetCol
165      const Track* pTrack = pPF->TrackerTrk();
166      if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
167      if(pTrack        ==  0                           ) continue;
168 <     lSumEt   += pPF->Pt();  //=> another bug
122 <     if(         !((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) ||
168 >    if(  !((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) ||
169             (pPF->HasGsfTrk()     && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex))    <iDZCut)))) continue;
170 +    lSumEt   += pPF->Pt();  
171      lVec     -= pPF->Mom();
172    }
173    int lNPass = 0;
174    for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
175      const PFJet *pJet = iJets->At(i0);
176 <    if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2))                                       continue; //Quick cleaning==> if not done already
130 <    double lPt = fJetIDMVA->correctedPt(pJet,iJetCorrector,iPileupEnergyDensity);
131 <    //if(fJetIDMVA->pass(pJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity))
132 <    //  std::cout << " =====> Jet Passes Id : " << lPt << " -- " << pJet->Eta() << " --- " << fJetIDMVA->pass(pJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity) << std::endl;
133 <    //if(!fJetIDMVA->pass(pJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity))
134 <    //  std::cout << " =====> Jet Fails  Id :  " << lPt << " -- " << pJet->Eta() << " --- " << fJetIDMVA->pass(pJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity) << std::endl;
176 >    if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2))                                       continue;
177      if(!fJetIDMVA->pass(pJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity)) continue;
178      addNeut(pJet,lVec,lSumEt,iJetCorrector,iPileupEnergyDensity);
179      lNPass++;
# Line 140 | Line 182 | Met RecoilTools::NoPUMet( const PFJetCol
182    lMet.SetSumEt( lSumEt);
183    return lMet;
184   }
185 < //--------------------------------------------------------------------------------------------------
185 > //-------------------------------------------------------------------------------------------------
186   //Corrected Jets
187   Met RecoilTools::NoPUMet( const PFJetCol       *iJets,
188                            const PFCandidateCol *iCands,const Vertex *iVertex,const VertexCol *iVertices,Double_t iRho,
# Line 152 | Line 194 | Met RecoilTools::NoPUMet( const PFJetCol
194      const Track* pTrack = pPF->TrackerTrk();
195      if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
196      if(pTrack        ==  0                           ) continue;
155 –    lSumEt   += pPF->Pt();  //=> another bug
197      if(  !((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) ||
198             (pPF->HasGsfTrk()     && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex))    <iDZCut)))) continue;
199      lVec     -= pPF->Mom();
200 +    lSumEt   += pPF->Pt();  
201    }
202    for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
203      const PFJet *pJet = iJets->At(i0);
204 <    if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2))                             continue; //Quick cleaning==> if not done already
163 <    //if(fJetIDMVA->pass(pJet,iVertex,iVertices))
164 <    //  std::cout << " =====> Jet Passes Id : " << pJet->Pt() << " -- " << pJet->Eta() << " --- " << fJetIDMVA->pass(pJet,iVertex,iVertices) << std::endl;
165 <    //if(!fJetIDMVA->pass(pJet,iVertex,iVertices))
166 <    //  std::cout << " =====> Jet Fails  Id :  " << pJet->Pt() << " -- " << pJet->Eta() << " --- " << fJetIDMVA->pass(pJet,iVertex,iVertices) << std::endl;
204 >    if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2))                             continue;
205      if(!fJetIDMVA->pass(pJet,iVertex,iVertices))                          continue;
206      addNeut(pJet,lVec,lSumEt,iRho);
207    }
# Line 172 | Line 210 | Met RecoilTools::NoPUMet( const PFJetCol
210    return lMet;
211   }
212   //--------------------------------------------------------------------------------------------------
213 + Met RecoilTools::NoPURecoil(Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
214 +                            const PFJetCol            *iJets,
215 +                            const PFCandidateCol   *iCands   ,
216 +                            const Vertex *iVertex,const VertexCol *iVertices,
217 +                            FactorizedJetCorrector *iJetCorrector,
218 +                            const PileupEnergyDensityCol *iPileupEnergyDensity,Double_t iDZCut) {
219 +
220 +  FourVectorM lVec        (0,0,0,0); double lSumEt          = 0;
221 +  for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) {
222 +    const PFCandidate *pPF = iCands->At(i0);
223 +    if(!filter(pPF,iPhi1,iEta1,iPhi2,iEta2)) continue;
224 +    const Track* pTrack = pPF->TrackerTrk();
225 +    if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
226 +    if(pTrack        ==  0                           ) continue;
227 +    if(  !((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) ||
228 +           (pPF->HasGsfTrk()     && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex))    <iDZCut)))) continue;
229 +    lSumEt   += pPF->Pt();  
230 +    lVec     -= pPF->Mom();
231 +  }
232 +  int lNPass = 0;
233 +  for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
234 +    const PFJet *pJet = iJets->At(i0);
235 +    if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2))                                       continue;
236 +    if(iJetCorrector != 0 && iPileupEnergyDensity != 0) if(!fJetIDMVA->pass(pJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity)) continue;
237 +    if(iJetCorrector == 0 || iPileupEnergyDensity == 0) if(!fJetIDMVA->pass(pJet,iVertex,iVertices)) continue;
238 +    if(iJetCorrector != 0 && iPileupEnergyDensity != 0) addNeut(pJet,lVec,lSumEt,iJetCorrector,iPileupEnergyDensity);
239 +    if(iJetCorrector == 0 || iPileupEnergyDensity == 0) addNeut(pJet,lVec,lSumEt,1.);
240 +    lNPass++;
241 +  }
242 +  Met lMet(lVec.Px(),lVec.Py());
243 +  lMet.SetSumEt( lSumEt);
244 +  return lMet;
245 + }
246 +
247 + //--------------------------------------------------------------------------------------------------
248   Met RecoilTools::NoPURecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
249                              const PFJetCol            *iJets,FactorizedJetCorrector *iJetCorrector,
250                              const PileupEnergyDensityCol *iPileupEnergyDensity,
# Line 271 | Line 344 | Met RecoilTools::PUCMet( const PFJetCol
344    lMet.SetSumEt(lSumEt);
345    return lMet;
346   }
347 + //--------------------------------------------------------------------------------------------------
348 + Met RecoilTools::PUCRecoil( Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
349 +                            const PFJetCol            *iJets,
350 +                            const PFCandidateCol   *iCands   ,
351 +                            const Vertex *iVertex,const VertexCol *iVertices,
352 +                            FactorizedJetCorrector *iJetCorrector,
353 +                           const PileupEnergyDensityCol *iPileupEnergyDensity,Double_t iDZCut) {
354 +
355 +  
356 +  FourVectorM lVec        (0,0,0,0); double lSumEt          = 0;
357 +  for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) {
358 +    const PFCandidate *pPF = iCands->At(i0);
359 +    if(!filter(pPF,iPhi1,iEta1,iPhi2,iEta2)) continue;    
360 +    const Track* pTrack = pPF->TrackerTrk();
361 +    if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
362 +    if(pTrack == 0                                   &&
363 +       (pPF->PFType() == PFCandidate::eGamma         ||
364 +        pPF->PFType() == PFCandidate::eEGammaHF      ||
365 +        pPF->PFType() == PFCandidate::eNeutralHadron ||
366 +        pPF->PFType() == PFCandidate::eHadronHF      ))
367 +      {lVec -= pPF->Mom(); lSumEt += pPF->Pt();}
368 +    if(pTrack        ==  0                           ) continue;
369 +    if(  !((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) ||
370 +           (pPF->HasGsfTrk()     && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex))    <iDZCut)))) continue;
371 +    lVec     -= pPF->Mom();
372 +    lSumEt   += pPF->Pt();
373 +  }
374 +  for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
375 +    const PFJet *pJet = iJets->At(i0);
376 +    if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2))                                      continue; //Quick cleaning==> if not done already
377 +    if(!JetTools::passPFLooseId(pJet))                                             continue;
378 +    if(iJetCorrector != 0 && iPileupEnergyDensity != 0) if(fJetIDMVA->pass(pJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity)) continue;
379 +    if(iJetCorrector == 0 || iPileupEnergyDensity == 0) if(fJetIDMVA->pass(pJet,iVertex,iVertices)) continue;
380 +    if(iJetCorrector != 0 && iPileupEnergyDensity != 0) addNeut(pJet,lVec,lSumEt,iJetCorrector,iPileupEnergyDensity,-1);
381 +    if(iJetCorrector == 0 || iPileupEnergyDensity == 0) addNeut(pJet,lVec,lSumEt,1,-1);
382 +  }
383 +  Met lMet(lVec.Px(),lVec.Py());
384 +  lMet.SetSumEt(lSumEt);
385 +  return lMet;
386 + }
387   //----> This MET is a bug need to fix it
388   //--------------------------------------------------------------------------------------------------
389   Met RecoilTools::PUCRecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
# Line 357 | Line 470 | Met RecoilTools::PUMet( const PFJetCol
470    }
471    Met lMet(lVec.Px(),lVec.Py());
472    lMet.SetSumEt(lSumEt);
473 +  return lMet;
474 + }
475 + //--------------------------------------------------------------------------------------------------
476 + Met RecoilTools::PUMet( Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
477 +                        const PFJetCol            *iJets,
478 +                        const PFCandidateCol   *iCands   ,
479 +                        const Vertex *iVertex,const VertexCol *iVertices,
480 +                        FactorizedJetCorrector *iJetCorrector,
481 +                        const PileupEnergyDensityCol *iPileupEnergyDensity,Double_t iDZCut) {
482 +  
483 +  FourVectorM lVec        (0,0,0,0); double lSumEt          = 0;
484 +  for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) {
485 +    const PFCandidate *pPF = iCands->At(i0);
486 +    if(!filter(pPF,iPhi1,iEta1,iPhi2,iEta2)) continue;    
487 +    const Track* pTrack = pPF->TrackerTrk();
488 +    if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
489 +    if(pTrack        ==  0                           ) continue;
490 +    if(   ((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) ||
491 +           (pPF->HasGsfTrk()     && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex))    <iDZCut)))) continue;
492 +    lVec     -= pPF->Mom();
493 +    lSumEt   += pPF->Pt();
494 +  }
495 +  for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
496 +    const PFJet *pJet = iJets->At(i0);
497 +    if(!JetTools::passPFLooseId(pJet))                                             continue;
498 +    if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2))                                      continue; //Quick cleaning
499 +    if(fJetIDMVA->pass(pJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity)) continue;
500 +    addNeut(pJet,lVec,lSumEt,iJetCorrector,iPileupEnergyDensity);
501 +  }
502 +  Met lMet(lVec.Px(),lVec.Py());
503 +  lMet.SetSumEt(lSumEt);
504    return lMet;
505   }
506   //--------------------------------------------------------------------------------------------------

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines