ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/interface/MuonTools.h
(Generate patch)

Comparing UserCode/MitPhysics/Utils/interface/MuonTools.h (file contents):
Revision 1.2 by ceballos, Tue Nov 11 21:21:28 2008 UTC vs.
Revision 1.7 by loizides, Tue Apr 7 15:37:09 2009 UTC

# Line 1 | Line 1
1 + //--------------------------------------------------------------------------------------------------
2 + // $Id$
3 + //
4 + // MuonTools
5 + //
6 + // This class allows you to classify a given muon according to defined criteria. For this purpose
7 + // is loads histograms from two ROOT files, specified in the constructor. The main function then
8 + // is "IsGood(*muon, selection)" which returns true if the given muon fulfills the selection
9 + // criteria.
10 + //
11 + // Logically, the code has been put together by Phil who took most of the ideas from
12 + //  http://cmslxr.fnal.gov/lxr/source/RecoMuon/MuonIdentification/
13 + //  http://cmslxr.fnal.gov/lxr/source/DataFormats/MuonReco/
14 + //
15 + // Authors: P.Harris, C.Loizides
16 + //--------------------------------------------------------------------------------------------------
17 +
18   #ifndef MITPHYSICS_UTIL_MUONTOOLS_H
19   #define MITPHYSICS_UTIL_MUONTOOLS_H
3 #include <vector>
4 #include <iostream>
20  
21 < #include "TH2D.h"
7 < #include "TFile.h"
21 > #include "MitAna/DataTree/interface/Muon.h"
22   #include "MitAna/DataTree/interface/Collections.h"
23   #include "MitCommon/MathTools/interface/MathUtils.h"
24 + #include "TH2D.h"
25  
26   namespace mithep {
27    class MuonTools {
28 <  public:
29 <    MuonTools();
30 <    ~MuonTools();
31 <    double getCaloCompatability(mithep::Muon* iMuon,bool iEMSpecial, bool iCorrectedHCAL);
32 <    enum selection {
33 <      AllArbitrated,
34 <      PromptTight,
35 <      TMLastStationLoose,
36 <      TMLastStationTight,
37 <      TMOneStationLoose,
38 <      TMOneStationTight,
39 <      TM2DCompatibilityLoose,
40 <      TM2DCompatibilityTight
41 <    };
42 <    static inline double sigWeight(double iVal0,double iVal1) {
43 <      if(iVal1 < 1.) return 1;
44 <      if(iVal0 < 3. && iVal1 > 3.) {
45 <        double lVal = TMath::Max(iVal0,1.);
46 <        return 1/TMath::Power(lVal,0.25);
47 <      }
48 <      double lVal = TMath::Max(iVal1,1.);
49 <      return 1/TMath::Power(lVal,0.25);
50 <    }
51 <    static inline int lastHit(mithep::Muon *iMuon) {
52 <      int lId = -1;
53 <      for(int i0 = 0; i0 < 8; i0++) {
54 <        if(iMuon->GetDX(i0) < 99999 || iMuon->GetDY(i0) < 99999) lId = i0;
55 <      }
56 <      return lId;
57 <    }
58 <    static inline int maxChamberId(mithep::Muon* iMuon,double iMaxD,double iMaxP) {
59 <      int lId = -1;
60 <      for(int i0 = 0; i0 < 8; i0++) {
61 <        if(iMuon->GetTrackDist(i0)                            < iMaxD &&
62 <           iMuon->GetTrackDist(i0)/iMuon->GetTrackDistErr(i0) < iMaxP) {
63 <          if(iMuon->GetDX(i0) >= 999999) {lId = i0;} else {lId = -1;}
64 <        }
65 <      }
66 <      return lId;
67 <    }
68 <    static inline int lastStation(mithep::Muon* iMuon,double iMaxD,double iMaxP) {
69 <      int lId = -1;
70 <      for(int i0 = 0; i0 < 8; i0++) {
71 <        if((lId % 4) > (i0 % 4)) continue;
72 <        if(iMuon->GetTrackDist(i0)                            < iMaxD &&
73 <           iMuon->GetTrackDist(i0)/iMuon->GetTrackDistErr(i0) < iMaxP) lId = i0;
74 <      }
75 <      return lId;
76 <    }
77 <    static inline int lastStation(mithep::Muon* iMuon,int iMax=8) {
78 <      int lId = -1; if(iMax > 8) iMax = 8;
79 <      for(int i0 = 0; i0 < iMax; i0++) {if(iMuon->StationBit(i0) && ((lId % 4) < (i0 % 4)))lId = i0;}
65 <      return lId;
66 <    }
67 <    static inline bool overflow(TH2D *iHist,double lVal0,double lVal1) {
68 <      if(iHist == 0) return true;
69 <      if(iHist ->GetXaxis()->FindBin(lVal0) == 0               ||
70 <         iHist ->GetXaxis()->FindBin(lVal0) >  iHist->GetNbinsX() ||
71 <         iHist ->GetYaxis()->FindBin(lVal0) == 0               ||
72 <         iHist ->GetYaxis()->FindBin(lVal0) >  iHist->GetNbinsY()) return true;
73 <      return false;
74 <    }
75 <    static inline MCParticle* etaPhiMatch(MCParticleCol* iMCParts,Muon *iMuon) {
76 <      mithep::MCParticle *lMC = 0; double lDR = 1.; double lPt = -1.;
77 <      for(unsigned int i0 = 0; i0 < iMCParts->GetEntries(); i0++) {
78 <        mithep::MCParticle* pMC = iMCParts->At(i0);
79 <        if(pMC->Status() != 1) continue;
80 <        double pDR = mithep::MathUtils::DeltaR(pMC->Mom(),iMuon->Mom());
81 <        if(pDR > lDR && pDR > 0.01) continue;
82 <        if(fabs(pMC->Pt()) < lPt)   continue;
83 <        lDR = pDR;
84 <        lMC = pMC;
85 <        lPt = pMC->Pt();
86 <      }
87 <      return lMC;
88 <    }
89 <    static inline Muon* match(MuonCol* iMuons,Track *iTrack) {
90 <      mithep::Muon *lMuon = 0; double lDR = 1.; double lPt = -1.;
91 <      for(unsigned int i0 = 0; i0 < iMuons->GetEntries(); i0++) {
92 <        mithep::Muon* pMuon = iMuons->At(i0);
93 <        double pDR = mithep::MathUtils::DeltaR(pMuon->Mom(),iTrack->Mom4(0.109));
94 <        if(pDR > lDR && pDR > 0.01) continue;
95 <        if(fabs(pMuon->Pt()) < lPt)   continue;
96 <        lDR   = pDR;
97 <        lMuon = pMuon;
98 <        lPt   = pMuon->Pt();
99 <      }
100 <      return lMuon;
101 <    }
102 <    static inline int NSegments(mithep::Muon *iMuon) {
103 <      int lNSegs = 0;
104 <      for(int i0 = 0; i0 < 8; i0++) {
105 <        if(iMuon->StationBit(i0)) lNSegs++;
106 <      }
107 <      return lNSegs;
108 <    }
109 <    static inline bool promptTight(int iId,Muon *iMuon) {
110 <      const mithep::Track *lTrack = 0;
111 <      if(iId == 0                              ) lTrack = iMuon->GlobalTrk();
112 <      if(iId == 1 || (iId == -1 && lTrack == 0)) lTrack = iMuon->TrackerTrk();
113 <      if(iId == 2 || (iId == -1 && lTrack == 0)) lTrack = iMuon->StandaloneTrk();
114 <      if(lTrack == 0) return 0;
115 <      if(lTrack->NHits() < 11)                                        return false;
116 <      if(lTrack->Chi2()/lTrack->Ndof() > 10.)                         return false;
117 <      if(lTrack->D0() > 0.2)                                          return false;
118 <      if((lastHit(iMuon)% 4) == 0)                                    return false;
119 <      return true;
120 <    }
121 <    static inline bool TMLastStation(Muon *iMuon,double iDYMin = 3., double iPYMin = 3.,double iDXMin = 3., double iPXMin = 3.,int iN = 2) {
122 <      if(iMuon->NSegments()       < iN)               return false; //2Last 1One
123 <      int lLast = lastStation(iMuon,-3.,-3.);                       //Last Required Station
124 <      if(lLast < 0)                                   return false;
125 <      if(iMuon->GetDX(lLast)  > 9999.)                return false;
126 <      lLast = lastStation(iMuon);                                   //No Requirements Imply StationMask (Segment and Track Arbitrated)
127 <      if(lLast < 0)                                   return false;
128 <      if(!(fabs(iMuon->GetDX(lLast))      < iDXMin  ||
129 <           fabs(iMuon->GetPullX(lLast))   < iPXMin))  return false;
130 <      if(lLast == 3) lLast = lastStation(iMuon,3);
131 <      if(lLast < 0)                                   return false;
132 <      if(!(fabs(iMuon->GetDY(lLast))      < iDYMin ||               //Old Code it was 9999
133 <           fabs(iMuon->GetPullY(lLast))   < iPYMin))  return false;
134 <      return true;
135 <    }
136 <    static inline bool TMOneStation(Muon *iMuon,double iDYMin = 3., double iPYMin = 3.,double iDXMin = 3., double iPXMin = 3.,int iN = 1) {
137 <      if(iMuon->NSegments()       < iN)               return false; //2Last 1One
138 <      bool pGoodX  = false; bool pBadY = false;
139 <      for(int i0 = 0; i0 < 8; i0++) {
140 <        if((fabs(iMuon->GetDX(i0))      < iDXMin  ||
141 <            fabs(iMuon->GetPullX(i0))   < iPXMin)) pGoodX = true;
142 <        if(pGoodX                                 &&
143 <           (fabs(iMuon->GetDY(i0))      < iDYMin  ||
144 <            fabs(iMuon->GetPullY(i0))   < iPYMin))  return true;
145 <        if(fabs(iMuon->GetDY(i0)) < 999999) pBadY = true;
146 <        if(i0 == 3 && pGoodX && !pBadY)             return true;
147 <      }
148 <      return false;
149 <    }
150 <    static inline double getCaloCompatabilitySlowly(mithep::Muon* iMuon,bool iEMSpecial, bool iCorrectedHCAL) {
151 <      std::cout << "Warning Loading Compatability Slowly" << std::endl;
152 <      TFile* lPion_templates = new TFile("MitCommon/MuonTools/data/PionCaloTemplate.root","READ");
153 <      TFile* lMuon_templates = new TFile("MitCommon/MuonTools/data/MuonCaloTemplate.root","READ");
154 <      TH2D* lpion_em_etaEmi  = (TH2D*) lPion_templates->Get("em_etaEmi");
155 <      TH2D* lpion_had_etaEmi = (TH2D*) lPion_templates->Get("had_etaEmi");
156 <      TH2D* lpion_had_etaTmi = (TH2D*) lPion_templates->Get("had_etaTmi");
157 <      TH2D* lpion_em_etaB    = (TH2D*) lPion_templates->Get("em_etaB")   ;
158 <      TH2D* lpion_had_etaB   = (TH2D*) lPion_templates->Get("had_etaB")  ;
159 <      TH2D* lpion_ho_etaB    = (TH2D*) lPion_templates->Get("ho_etaB")   ;
160 <      TH2D* lpion_had_etaTpl = (TH2D*) lPion_templates->Get("had_etaTpl");
161 <      TH2D* lpion_em_etaEpl  = (TH2D*) lPion_templates->Get("em_etaEpl") ;
162 <      TH2D* lpion_had_etaEpl = (TH2D*) lPion_templates->Get("had_etaEpl");
163 <      TH2D* lmuon_em_etaEmi  = (TH2D*) lMuon_templates->Get("em_etaEmi") ;
164 <      TH2D* lmuon_had_etaEmi = (TH2D*) lMuon_templates->Get("had_etaEmi");
165 <      TH2D* lmuon_had_etaTmi = (TH2D*) lMuon_templates->Get("had_etaTmi");
166 <      TH2D* lmuon_em_etaB    = (TH2D*) lMuon_templates->Get("em_etaB")   ;
167 <      TH2D* lmuon_had_etaB   = (TH2D*) lMuon_templates->Get("had_etaB")  ;
168 <      TH2D* lmuon_ho_etaB    = (TH2D*) lMuon_templates->Get("ho_etaB")   ;
169 <      TH2D* lmuon_had_etaTpl = (TH2D*) lMuon_templates->Get("had_etaTpl");
170 <      TH2D* lmuon_em_etaEpl  = (TH2D*) lMuon_templates->Get("em_etaEpl");
171 <      TH2D* lmuon_had_etaEpl = (TH2D*) lMuon_templates->Get("had_etaEpl");
172 <
173 <      double lEta = -1.; double lP = -1;
174 <      double lEM  = -5.;      double lHad = 0;      double lHO = 0;
175 <      lEta = iMuon->Eta();
176 <      lP   = iMuon->P();
177 <      if(lP >= 2000.) lP = 1999.9;
178 <      if(!iEMSpecial || iMuon->EmEnergy() != 0.) lEM  = iMuon->EmEnergy();
179 <      lHad = iMuon->HadEnergy();
180 <      lHO  = iMuon->HoEnergy();
181 <      if(lP < 0. )           return 0.5;
182 <      if(fabs(lEta) >  2.5 ) return 0.5;
183 <      TH2D* lTMuonHad = NULL;
184 <      TH2D* lTPionHad = NULL;
185 <      TH2D* lTMuonHo  = NULL;
186 <      TH2D* lTPionHo  = NULL;
187 <      TH2D* lTMuonEm  = NULL;
188 <      TH2D* lTPionEm  = NULL;
189 <
190 <      if(fabs(lEta) >=  1.27) {
191 <        if(iCorrectedHCAL) lHad *= 1.8/2.2;
192 <        if(lEta > 0) {
193 <          lTPionHad = lpion_had_etaEpl;
194 <          lTMuonHad = lmuon_had_etaEpl;
195 <        } else {
196 <          lTPionHad = lpion_had_etaEmi;
197 <          lTMuonHad = lmuon_had_etaEmi;
198 <        }
199 <      }
200 <      if(fabs(lEta) <  1.27  && fabs(lEta) >=  1.1 ) {
201 <        if(iCorrectedHCAL)    lHad *= (1.8/(-2.2*fabs(lEta)+5.5));
202 <        if(lEta > 0) {
203 <          lTPionHad  = lpion_had_etaTpl;
204 <          lTMuonHad  = lmuon_had_etaTpl;
205 <        } else {
206 <          lTPionHad  = lpion_had_etaTmi;
207 <          lTMuonHad  = lmuon_had_etaTmi;
208 <        }
209 <      }
210 <      if(fabs(lEta) <  1.1) {
211 <        if(iCorrectedHCAL)    lHad *= sin(2*atan(exp(iMuon->Eta())));
212 <        lTPionHad  = lpion_had_etaB;
213 <        lTMuonHad  = lmuon_had_etaB;
214 <      }
215 <      if(lEta >  1.479  ) {
216 <        lTPionEm  = lpion_em_etaEpl;
217 <        lTMuonEm  = lmuon_em_etaEpl;
218 <      }
219 <      if(fabs(lEta) <=  1.479) {
220 <        lTPionEm  = lpion_em_etaB;
221 <        lTMuonEm  = lmuon_em_etaB;
222 <      }
223 <      if(lEta < -1.479 ) {
224 <        lTPionEm  = lpion_em_etaEmi;
225 <        lTMuonEm  = lmuon_em_etaEmi;
226 <      }
227 <      if(fabs(lEta) < 1.28) {
228 <        lTPionHo  = lpion_ho_etaB;
229 <        lTMuonHo  = lmuon_ho_etaB;
230 <      }
231 <
232 <      double lPBX = 1.;     double lPSX = 1.;
233 <      double lPBY = 1.;     double lPSY = 1.;
234 <      double lPBZ = 1.;     double lPSZ = 1.;
235 <      if(!overflow(lTPionEm, lP,lEM))  lPBX =  lTPionEm ->GetBinContent(lTPionEm ->GetXaxis()->FindBin(lP),lTPionEm ->GetYaxis()->FindBin(lEM) );
236 <      if(!overflow(lTPionHad,lP,lHad)) lPBY =  lTPionHad->GetBinContent(lTPionHad->GetXaxis()->FindBin(lP),lTPionHad->GetYaxis()->FindBin(lHad));
237 <      if(!overflow(lTPionHo, lP,lHO))  lPBZ =  lTPionHo ->GetBinContent(lTPionHo ->GetXaxis()->FindBin(lP),lTPionHo ->GetYaxis()->FindBin(lHO) );
238 <      if(!overflow(lTMuonEm, lP,lEM )) lPSX =  lTMuonEm ->GetBinContent(lTMuonEm ->GetXaxis()->FindBin(lP),lTMuonEm ->GetYaxis()->FindBin(lEM) );
239 <      if(!overflow(lTMuonHad,lP,lHad)) lPSY =  lTMuonHad->GetBinContent(lTMuonHad->GetXaxis()->FindBin(lP),lTMuonHad->GetYaxis()->FindBin(lHad));
240 <      if(!overflow(lTMuonHo ,lP,lHO))  lPSZ =  lTMuonHo ->GetBinContent(lTMuonHo ->GetXaxis()->FindBin(lP),lTMuonHo ->GetYaxis()->FindBin(lHO) );
241 <
242 <      if(lPSX == 0. || lPBX == 0. || (lEM <= 0. && !iEMSpecial)) {lPSX = 1.; lPBX = 1.;}
243 <      if(lPSY == 0. || lPBY == 0. || lHad == 0.) {lPSY = 1.; lPBY = 1.;}
244 <      if(lPSZ == 0. || lPBZ == 0. || lHO  == 0.) {lPSZ = 1.; lPBZ = 1.;}
245 <      lPion_templates->Close();
246 <      lMuon_templates->Close();
247 <      if((lPSX*lPSY*lPSZ+lPBX*lPBY*lPBZ) > 0.) return lPSX*lPSY*lPSZ/(lPSX*lPSY*lPSZ+lPBX*lPBY*lPBZ);
248 <      return 0.5;
249 <    }
250 <    static inline float getSegmentCompatability(const mithep::Muon* iMuon) {
251 <      int lNStationsCrossed = 0;
252 <      int lNStationsSegment = 0;
253 <      
254 <      std::vector<int>   lStTrack(8);
255 <      std::vector<int>   lStSegmentmatch(8);
256 <      std::vector<int>   lStCrossed(8);
257 <      std::vector<float> lStBoundary(8);
258 <      std::vector<float> lStWeight(8);
259 <      float lWeight    = 0.;
260 <      bool  lAdjust    = true;
261 <      for(int i0 = 0; i0 < 8; i0++) {
262 <        if(iMuon->GetTrackDist(i0) < 999999. ) {
263 <          lNStationsCrossed++;
264 <          lStCrossed[i0]  = 1;
265 <          lStBoundary[i0] = 0.;
266 <          if(iMuon->GetTrackDist(i0) > -10. ) lStBoundary[i0] = iMuon->GetTrackDist(i0);
267 <        }
268 <        if(iMuon->GetDX(i0) < 999999.) { //Use iMuon->GetSegmentX--> CHECK
269 <          lNStationsSegment++;
270 <          lStSegmentmatch[i0] = 1;
271 <        }
272 <        if(iMuon->GetDY(i0) < 999999.) lAdjust = false;
273 <      }
274 <      int lPCross = -1;
275 <      const float lAtWeight = 0.5;
276 <      for(int i0 = 0; i0< 8; i0++)     {
277 <        lStWeight[i0] = 0;
278 <        if(lStCrossed[i0] > 0 )        { lPCross++;
279 <        switch ( lNStationsCrossed ) {
280 <        case 1 :
281 <          lStWeight[i0] =  1.;
282 <          break;
283 <        case 2 :
284 <          if     ( lPCross  == 0 ) lStWeight[i0] =  0.33;
285 <          else   lStWeight[i0] =  0.67;
286 <          break;
287 <        case 3 :
288 <          if     ( lPCross == 0 ) lStWeight[i0] =  0.23;
289 <          else if( lPCross == 1 ) lStWeight[i0] =  0.33;
290 <          else                    lStWeight[i0] =  0.44;
291 <          break;
292 <        case 4 :
293 <          if     ( lPCross == 0 ) lStWeight[i0] =  0.10;
294 <          else if( lPCross == 1 ) lStWeight[i0] =  0.20;
295 <          else if( lPCross == 2 ) lStWeight[i0] =  0.30;
296 <          else                    lStWeight[i0] =  0.40;
297 <          break;
298 <          
299 <        default :
300 <          lStWeight[i0] = 1./lNStationsCrossed;
301 <        }
302 <        
303 <        if(lStSegmentmatch[i0] <= 0 && lStBoundary[i0] != 0. ) {
304 <          lStWeight[i0] *= lAtWeight*0.5*(TMath::Erf(lStBoundary[i0]/6.)+1.);
305 <        } else if(lStSegmentmatch[i0] <= 0 && lStBoundary[i0] == 0) {lStWeight[i0] = 0.;}
306 <        
307 <        if( lStSegmentmatch[i0] > 0) {
308 <          double lP2X = TMath::Power(iMuon->GetPullX(i0),2.);
309 <          double lP2Y = TMath::Power(iMuon->GetPullY(i0),2.);
310 <          double lD2X = TMath::Power(iMuon->GetDX(i0),2.);
311 <          double lD2Y = TMath::Power(iMuon->GetDY(i0),2.);
312 <          if( iMuon->GetDY(i0) < 999999 && iMuon->GetDX(i0) < 999999) {
313 <            lStWeight[i0] *= sigWeight(sqrt(lD2X+lD2Y),sqrt(lP2X+lP2Y));
314 <          } else if (iMuon->GetDY(i0) >= 999999 && i0 < 4) {
315 <            lStWeight[i0] *= sigWeight(iMuon->GetDX(i0),iMuon->GetPullX(i0));
316 <          } else if(i0 < 4) {
317 <            lStWeight[i0] *= sigWeight(iMuon->GetDY(i0),iMuon->GetPullY(i0));
318 <          }
319 <        }
320 <        }
321 <        lWeight += lStWeight[i0];
322 <      }
323 <      if(lNStationsCrossed == 0) lWeight = 0.5;
324 <      return lWeight;
325 <    }
326 <    static inline bool isGoodMuon(mithep::Muon *iMuon,selection iSelection) {
327 <      double lVal = 0;
328 <      switch(iSelection) {
329 <      case AllArbitrated:
330 <        if(iMuon->StandaloneTrk() != 0 || iMuon->GlobalTrk()!= 0)  return true;
331 <        if(iMuon->NSegments() > 0) return true;
332 <        return false;
333 <        break;
334 <      case PromptTight:
335 <        return promptTight(-1,iMuon);
336 <        break;
337 <      case TMOneStationLoose:
338 <        return TMOneStation(iMuon,99999,999999);
339 <      case TMOneStationTight:
340 <        return TMOneStation(iMuon);
341 <      case TMLastStationLoose:
342 <        return TMLastStation(iMuon,999999,999999);
343 <        break;
344 <      case TMLastStationTight:
345 <        return TMLastStation(iMuon);
346 <        break;
347 <      case TM2DCompatibilityLoose:
348 <        lVal             = 1.2*getSegmentCompatability(iMuon);
349 <        if(lVal/1.2 == 0.5) return false;
350 <        //if(fClass ) lVal += 0.8*getCaloCompatability(iMuon,true,true);
351 <        lVal += 0.8*getCaloCompatabilitySlowly(iMuon,true,true);
352 <        if(lVal > 0.7) return true;
353 <        return false;
354 <        break;
355 <      case TM2DCompatibilityTight:
356 <        lVal             = 1.2*getSegmentCompatability(iMuon);
357 <        if(lVal/1.2 == 0.5) return false;
358 <        //if(fClass ) lVal += 0.8*getCaloCompatability(iMuon,true,true);
359 <        lVal += 0.8*getCaloCompatabilitySlowly(iMuon,true,true);
360 <        if(lVal > 1.0) return true;
361 <        return false;
362 <        break;
363 <      default:
364 <        return false;
365 <      }
366 <      return false;
367 <    }
368 <    bool isGood(mithep::Muon *iMuon,selection iSelection);
369 <  private:
370 <    TFile *fMuon_templates;
371 <    TFile *fPion_templates;
372 <    TH2D* fpion_em_etaEmi;
373 <    TH2D* fpion_had_etaEmi;
374 <    TH2D* fpion_had_etaTmi;
375 <    TH2D* fpion_em_etaB;
376 <    TH2D* fpion_had_etaB;
377 <    TH2D* fpion_ho_etaB;
378 <    TH2D* fpion_had_etaTpl;
379 <    TH2D* fpion_em_etaEpl;
380 <    TH2D* fpion_had_etaEpl;
381 <    TH2D* fmuon_em_etaEmi;
382 <    TH2D* fmuon_had_etaEmi;
383 <    TH2D* fmuon_had_etaTmi;
384 <    TH2D* fmuon_em_etaB;
385 <    TH2D* fmuon_had_etaB;
386 <    TH2D* fmuon_ho_etaB;
387 <    TH2D* fmuon_had_etaTpl;
388 <    TH2D* fmuon_em_etaEpl;
389 <    TH2D* fmuon_had_etaEpl;
28 >    public:
29 >      MuonTools(const char *mutemp="$CMSSW_BASE/src/MitPhysics/data/MuonCaloTemplate.root",
30 >                const char *pitemp="$CMSSW_BASE/src/MitPhysics/data/PionCaloTemplate.root");
31 >      ~MuonTools();
32 >
33 >      enum ESelType {
34 >        kAllArbitrated,          //All arbitration (DT/CSC/RPC Hits) put on at least one
35 >                                 //  segments given a global Muon
36 >        kPromptTight,            //Standard global muon identification
37 >        kTMLastStationLoose,     //Loose matching requirements on lastmost muon station of reco
38 >        kTMLastStationTight,     //Tight matching requirements on lastmost muon station of reco
39 >        kTMOneStationLoose,      //Loose matching requirements on at least one muon station of reco
40 >        kTMOneStationTight,      //Tight matching requirements on at least one muon station of reco
41 >        kTM2DCompatibilityLoose, //Loose requirement on sum of compatabiliity variables
42 >                                 //  ===> 1.2 Segment compatability + 0.8 calo compatability > 0.8
43 >        kTM2DCompatibilityTight  //Tight requirement on sum of compatabiliity variables
44 >                                 //  ===> 1.2 Segment compatability + 0.8 calo compatability > 1.2
45 >      };
46 >
47 >      Bool_t      Init(const char *mutemp, const char *pitemp);
48 >      Bool_t      IsGood(const mithep::Muon *iMuon, ESelType iSel) const;
49 >      Double_t    GetCaloCompatability(const mithep::Muon *iMuon,
50 >                                       Bool_t iEMSpecial, Bool_t iCorrectedHCAL) const;
51 >      Double_t    GetSegmentCompatability(const mithep::Muon *iMuon)             const;
52 >
53 >    protected:
54 >      void        DeleteHistos();
55 >      Bool_t      Overflow(const TH2D *iHist, Double_t lVal0, Double_t lVal1)    const;
56 >      Double_t    SigWeight(Double_t iVal0, Double_t iVal1)                      const;
57 >  
58 >    private:
59 >      Bool_t      fIsInit;              //!=true if histograms are loaded
60 >      TH2D       *fmuon_em_etaEmi;      //!Neg Endcap EM       Calo Deposit Template for Muons
61 >      TH2D       *fmuon_had_etaEmi;     //!Neg Endcap Hadronic Calo Deposit Template for Muons
62 >      TH2D       *fmuon_had_etaTmi;     //!Neg Transition Hadronic Calo Deposit Template for Muons
63 >      TH2D       *fmuon_em_etaB;        //!Barrel EM       Calo Deposit Template for Muons
64 >      TH2D       *fmuon_had_etaB;       //!Barrel Hadronic Calo Deposit Template for Muons
65 >      TH2D       *fmuon_ho_etaB;        //!Barrel HO       Calo Deposit Template for Muons
66 >      TH2D       *fmuon_had_etaTpl;     //!Plus Transition Hadronic Calo Deposit Template for Muons
67 >      TH2D       *fmuon_em_etaEpl;      //!Plus Endcap EM       Calo Deposit Template for Muons
68 >      TH2D       *fmuon_had_etaEpl;     //!Plus Endcap Hadronic Calo Deposit Template for Muons
69 >      TH2D       *fpion_em_etaEmi;      //!Neg  Endcap EM       Calo Deposit Template for Pions
70 >      TH2D       *fpion_had_etaEmi;     //!Neg  Endcap Hadronic Calo Deposit Template for Pions
71 >      TH2D       *fpion_had_etaTmi;     //!Neg Transition Hadronic Calo Deposit Template for Pions
72 >      TH2D       *fpion_em_etaB;        //!Barrel EM       Calo Deposit Template for Pions
73 >      TH2D       *fpion_had_etaB;       //!Barrel Hadronic Calo Deposit Template for Pions
74 >      TH2D       *fpion_ho_etaB;        //!Barrel HO       Calo Deposit Template for Pions
75 >      TH2D       *fpion_had_etaTpl;     //!Plus Transition Hadronic Calo Deposit Template for Pions
76 >      TH2D       *fpion_em_etaEpl;      //!Plus Endcap EM       Calo Deposit Template for Pions
77 >      TH2D       *fpion_had_etaEpl;     //!Plus Endcap Hadronic Calo Deposit Template for Pions
78 >
79 >      TH2D       *LoadHisto(const char *fname, TFile *file)                      const;
80    };
81   }
82  
83 + //--------------------------------------------------------------------------------------------------
84 + inline Double_t mithep::MuonTools::SigWeight(Double_t iVal0, Double_t iVal1) const
85 + {
86 +  // Returns weighted uncertainty given segment matching uncertainty (iVal0) and
87 +  // segment matching pull (iVal1).
88 +
89 +  if (iVal1 < 1.)  //if pull defined and within range
90 +    return 1.;
91 +  if (iVal0 < 3. && iVal1 > 3.) {       //if pull not well defined and uncertainty defined
92 +    Double_t lVal = TMath::Max(iVal0,1.);
93 +    return 1./TMath::Power(lVal,0.25);
94 +  }
95 +
96 +  Double_t lVal = TMath::Max(iVal1,1.); //if pull > 1 and pull < 3 return 1/pull^4
97 +  return 1./TMath::Power(lVal,0.25);
98 + }
99 + //--------------------------------------------------------------------------------------------------
100 + inline Bool_t mithep::MuonTools::Overflow(const TH2D *iHist, Double_t lVal0, Double_t lVal1) const
101 + {
102 +  // Check if values are in overflow bins of given histogram.
103 +
104 +  if(iHist == 0)
105 +    return kTRUE;
106 +
107 +  if (iHist ->GetXaxis()->FindBin(lVal0) == 0                  ||
108 +      iHist ->GetXaxis()->FindBin(lVal0) >  iHist->GetNbinsX() ||
109 +      iHist ->GetYaxis()->FindBin(lVal0) == 0                  ||
110 +      iHist ->GetYaxis()->FindBin(lVal0) >  iHist->GetNbinsY()) {
111 +    return kTRUE;
112 +  }
113 +  return kFALSE;
114 + }
115   #endif

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines