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

Comparing UserCode/MitPhysics/Utils/src/MuonTools.cc (file contents):
Revision 1.3 by pharris, Mon Nov 3 16:55:48 2008 UTC vs.
Revision 1.14 by ceballos, Fri Jan 21 11:25:29 2011 UTC

# Line 1 | Line 1
1 + // $Id$
2 +
3   #include "MitPhysics/Utils/interface/MuonTools.h"
4 + #include <TFile.h>
5 +
6 + ClassImp(mithep::MuonTools)
7 +
8   using namespace mithep;
9  
10 < MuonTools::MuonTools() :
11 <  fpion_em_etaEmi(0),
12 <  fpion_had_etaEmi(0),
7 <  fpion_had_etaTmi(0),
8 <  fpion_em_etaB(0),
9 <  fpion_had_etaB(0),
10 <  fpion_ho_etaB(0),
11 <  fpion_had_etaTpl(0),
12 <  fpion_em_etaEpl(0),
13 <  fpion_had_etaEpl(0),
10 > //--------------------------------------------------------------------------------------------------
11 > MuonTools::MuonTools(const char *mutemp, const char *pitemp) :
12 >  fIsInit(kFALSE),
13    fmuon_em_etaEmi(0),
14    fmuon_had_etaEmi(0),
15    fmuon_had_etaTmi(0),
# Line 19 | Line 18 | MuonTools::MuonTools() :
18    fmuon_ho_etaB(0),
19    fmuon_had_etaTpl(0),
20    fmuon_em_etaEpl(0),
21 <  fmuon_had_etaEpl(0) {}
22 < MuonTools::~MuonTools() {
23 <  delete fpion_em_etaEmi;
24 <  delete fpion_had_etaEmi;
25 <  delete fpion_had_etaTmi;
26 <  delete fpion_em_etaB;
27 <  delete fpion_had_etaB;
28 <  delete fpion_ho_etaB;
29 <  delete fpion_had_etaTpl;
30 <  delete fpion_em_etaEpl;
31 <  delete fpion_had_etaEpl;
32 <  delete fmuon_em_etaEmi;
33 <  delete fmuon_had_etaEmi;
34 <  delete fmuon_had_etaTmi;
35 <  delete fmuon_em_etaB;
36 <  delete fmuon_had_etaB;
37 <  delete fmuon_ho_etaB;
38 <  delete fmuon_had_etaTpl;
39 <  delete fmuon_em_etaEpl;
40 <  delete fmuon_had_etaEpl;
41 < }
42 <
43 < double MuonTools::getCaloCompatability(mithep::Muon* iMuon,bool iEMSpecial, bool iCorrectedHCAL) {
44 <  if(fpion_em_etaEmi == 0) {
45 <    TFile* fPion_templates = new TFile("$CMSSW_BASE/src/MitPhysics/Init/PionCaloTemplate.root","READ");
46 <    TFile* fMuon_templates = new TFile("$CMSSW_BASE/src/MitPhysics/Init/MuonCaloTemplate.root","READ");
47 <    fpion_em_etaEmi  = (TH2D*) fPion_templates->Get("em_etaEmi");
48 <    fpion_had_etaEmi = (TH2D*) fPion_templates->Get("had_etaEmi");
49 <    fpion_had_etaTmi = (TH2D*) fPion_templates->Get("had_etaTmi");
50 <    fpion_em_etaB    = (TH2D*) fPion_templates->Get("em_etaB")   ;
51 <    fpion_had_etaB   = (TH2D*) fPion_templates->Get("had_etaB")  ;
52 <    fpion_ho_etaB    = (TH2D*) fPion_templates->Get("ho_etaB")   ;
53 <    fpion_had_etaTpl = (TH2D*) fPion_templates->Get("had_etaTpl");
54 <    fpion_em_etaEpl  = (TH2D*) fPion_templates->Get("em_etaEpl") ;
55 <    fpion_had_etaEpl = (TH2D*) fPion_templates->Get("had_etaEpl");
56 <    fmuon_em_etaEmi  = (TH2D*) fMuon_templates->Get("em_etaEmi") ;
57 <    fmuon_had_etaEmi = (TH2D*) fMuon_templates->Get("had_etaEmi");
58 <    fmuon_had_etaTmi = (TH2D*) fMuon_templates->Get("had_etaTmi");
59 <    fmuon_em_etaB    = (TH2D*) fMuon_templates->Get("em_etaB")   ;
60 <    fmuon_had_etaB   = (TH2D*) fMuon_templates->Get("had_etaB")  ;
61 <    fmuon_ho_etaB    = (TH2D*) fMuon_templates->Get("ho_etaB")   ;
62 <    fmuon_had_etaTpl = (TH2D*) fMuon_templates->Get("had_etaTpl");
63 <    fmuon_em_etaEpl  = (TH2D*) fMuon_templates->Get("em_etaEpl");
64 <    fmuon_had_etaEpl = (TH2D*) fMuon_templates->Get("had_etaEpl");
65 <
66 <    fpion_em_etaEmi ->SetDirectory(0);
67 <    fpion_had_etaEmi->SetDirectory(0);
68 <    fpion_had_etaTmi->SetDirectory(0);
69 <    fpion_em_etaB   ->SetDirectory(0);
70 <    fpion_had_etaB  ->SetDirectory(0);
71 <    fpion_ho_etaB   ->SetDirectory(0);
72 <    fpion_had_etaTpl->SetDirectory(0);
73 <    fpion_em_etaEpl ->SetDirectory(0);  
74 <    fpion_had_etaEpl->SetDirectory(0);  
75 <    fmuon_em_etaEmi ->SetDirectory(0);  
76 <    fmuon_had_etaEmi->SetDirectory(0);  
77 <    fmuon_had_etaTmi->SetDirectory(0);  
78 <    fmuon_em_etaB   ->SetDirectory(0);  
79 <    fmuon_had_etaB  ->SetDirectory(0);  
80 <    fmuon_ho_etaB   ->SetDirectory(0);  
81 <    fmuon_had_etaTpl->SetDirectory(0);  
82 <    fmuon_em_etaEpl ->SetDirectory(0);  
83 <    fmuon_had_etaEpl->SetDirectory(0);  
84 <    fPion_templates->Close();
85 <    fMuon_templates->Close();
86 <  }
87 <  double lEta = -1.; double lP = -1;
88 <  double lEM  = -5.;      double lHad = 0;      double lHO = 0;
89 <  lEta = iMuon->Eta();
90 <  lP   = iMuon->P();
91 <  if(lP >= 2000.) lP = 1999.9;
92 <  if(!iEMSpecial || iMuon->EmEnergy() != 0.) lEM  = iMuon->EmEnergy();
93 <  lHad = iMuon->HadEnergy();
94 <  lHO  = iMuon->HoEnergy();
95 <  if(lP < 0. )           return 0.5;
96 <  if(fabs(lEta) >  2.5 ) return 0.5;
97 <  TH2D* lTMuonHad = NULL;
98 <  TH2D* lTPionHad = NULL;
99 <  TH2D* lTMuonHo  = NULL;
100 <  TH2D* lTPionHo  = NULL;
101 <  TH2D* lTMuonEm  = NULL;
102 <  TH2D* lTPionEm  = NULL;
103 <  
104 <  if(fabs(lEta) >=  1.27) {
105 <    if(iCorrectedHCAL) lHad *= 1.8/2.2;
106 <    if(lEta > 0) {
21 >  fmuon_had_etaEpl(0),
22 >  fpion_em_etaEmi(0),
23 >  fpion_had_etaEmi(0),
24 >  fpion_had_etaTmi(0),
25 >  fpion_em_etaB(0),
26 >  fpion_had_etaB(0),
27 >  fpion_ho_etaB(0),
28 >  fpion_had_etaTpl(0),
29 >  fpion_em_etaEpl(0),
30 >  fpion_had_etaEpl(0)
31 > {
32 >  // Constructor.
33 >
34 >  if (mutemp && pitemp)
35 >    Init(mutemp, pitemp);
36 > }
37 >
38 > //--------------------------------------------------------------------------------------------------
39 > MuonTools::~MuonTools()
40 > {
41 >  // Destructor.
42 >
43 >  DeleteHistos();
44 > }
45 >
46 > //--------------------------------------------------------------------------------------------------
47 > void MuonTools::DeleteHistos()
48 > {
49 >  // Delete histograms.
50 >
51 >  if (fIsInit) {
52 >    delete fpion_em_etaEmi;
53 >    delete fpion_had_etaEmi;
54 >    delete fpion_had_etaTmi;
55 >    delete fpion_em_etaB;
56 >    delete fpion_had_etaB;
57 >    delete fpion_ho_etaB;
58 >    delete fpion_had_etaTpl;
59 >    delete fpion_em_etaEpl;
60 >    delete fpion_had_etaEpl;
61 >    delete fmuon_em_etaEmi;
62 >    delete fmuon_had_etaEmi;
63 >    delete fmuon_had_etaTmi;
64 >    delete fmuon_em_etaB;
65 >    delete fmuon_had_etaB;
66 >    delete fmuon_ho_etaB;
67 >    delete fmuon_had_etaTpl;
68 >    delete fmuon_em_etaEpl;
69 >    delete fmuon_had_etaEpl;
70 >    fpion_em_etaEmi  = 0;
71 >    fpion_had_etaEmi = 0;
72 >    fpion_had_etaTmi = 0;
73 >    fpion_em_etaB    = 0;
74 >    fpion_had_etaB   = 0;
75 >    fpion_ho_etaB    = 0;
76 >    fpion_had_etaTpl = 0;
77 >    fpion_em_etaEpl  = 0;
78 >    fpion_had_etaEpl = 0;
79 >    fmuon_em_etaEmi  = 0;
80 >    fmuon_had_etaEmi = 0;
81 >    fmuon_had_etaTmi = 0;
82 >    fmuon_em_etaB    = 0;
83 >    fmuon_had_etaB   = 0;
84 >    fmuon_ho_etaB    = 0;
85 >    fmuon_had_etaTpl = 0;
86 >    fmuon_em_etaEpl  = 0;
87 >    fmuon_had_etaEpl = 0;
88 >    fIsInit = kFALSE;
89 >  }
90 > }
91 >
92 > //--------------------------------------------------------------------------------------------------
93 > Double_t MuonTools::GetCaloCompatability(const Muon *iMuon,
94 >                                         Bool_t iEMSpecial, Bool_t iCorrectedHCAL) const
95 > {
96 >  // Get calo compatibility value for given muon based on calorimeter templates.
97 >  // If iEMSpecial is true, then a use different arrangement of ECAL for compatibility.
98 >
99 >  Double_t lEta = iMuon->Eta();
100 >  Double_t aEta = TMath::Abs(lEta);
101 >  if (aEta > 2.5)
102 >    return 0.5;
103 >
104 >  Double_t lP = iMuon->P();
105 >  if (lP >= 2000.)
106 >    lP = 1999.9;
107 >  if(lP < 0. )          
108 >    return 0.5;
109 >
110 >  Double_t lEM  = -5.;      
111 >  if (!iEMSpecial || iMuon->EmEnergy() != 0.)
112 >    lEM  = iMuon->EmEnergy();
113 >
114 >  Double_t lHad = iMuon->HadEnergy();
115 >  Double_t lHO = iMuon->HoEnergy();;
116 >
117 >  TH2D *lTMuonHad = 0;
118 >  TH2D *lTPionHad = 0;
119 >  TH2D *lTMuonHo  = 0;
120 >  TH2D *lTPionHo  = 0;
121 >  TH2D *lTMuonEm  = 0;
122 >  TH2D *lTPionEm  = 0;
123 >    
124 >  if (aEta >= 1.27) {
125 >  if (iCorrectedHCAL)
126 >      lHad *= 1.8/2.2;
127 >    if (lEta > 0) {
128        lTPionHad = fpion_had_etaEpl;
129        lTMuonHad = fmuon_had_etaEpl;
130      } else {
# Line 112 | Line 132 | double MuonTools::getCaloCompatability(m
132        lTMuonHad = fmuon_had_etaEmi;
133      }
134    }
135 <  if(fabs(lEta) <  1.27  && fabs(lEta) >=  1.1 ) {
136 <    if(iCorrectedHCAL)    lHad *= (1.8/(-2.2*fabs(lEta)+5.5));
137 <    if(lEta > 0) {
135 >
136 >  if (aEta < 1.27 && aEta >= 1.1) {
137 >    if (iCorrectedHCAL)    
138 >      lHad *= (1.8/(-2.2*aEta+5.5));
139 >    if (lEta > 0) {
140        lTPionHad  = fpion_had_etaTpl;
141        lTMuonHad  = fmuon_had_etaTpl;
142      } else {
# Line 122 | Line 144 | double MuonTools::getCaloCompatability(m
144        lTMuonHad  = fmuon_had_etaTmi;
145      }
146    }
147 <  if(fabs(lEta) <  1.1) {
148 <    if(iCorrectedHCAL)    lHad *= sin(2*atan(exp(iMuon->Eta())));
147 >
148 >  if (aEta < 1.1) {
149 >    if(iCorrectedHCAL)    
150 >      lHad *= TMath::Sin(2*TMath::ATan(TMath::Exp(lEta)));
151      lTPionHad  = fpion_had_etaB;
152      lTMuonHad  = fmuon_had_etaB;
153    }
154 <  if(lEta >  1.479  ) {
155 <    lTPionEm  = fpion_em_etaEpl;
156 <    lTMuonEm  = fmuon_em_etaEpl;
154 >  if (lEta > 1.479) {
155 >    lTPionEm = fpion_em_etaEpl;
156 >    lTMuonEm = fmuon_em_etaEpl;
157    }
158 <  if(fabs(lEta) <=  1.479) {
158 >  if (aEta <= 1.479) {
159      lTPionEm  = fpion_em_etaB;
160      lTMuonEm  = fmuon_em_etaB;
161    }
162 <  if(lEta < -1.479 ) {
162 >  if (lEta < -1.479) {
163      lTPionEm  = fpion_em_etaEmi;
164      lTMuonEm  = fmuon_em_etaEmi;
165    }
166 <  if(fabs(lEta) < 1.28) {
166 >  if (aEta < 1.28) {
167      lTPionHo  = fpion_ho_etaB;
168      lTMuonHo  = fmuon_ho_etaB;
169    }
170    
171 <  double lPBX = 1.;     double lPSX = 1.;
172 <  double lPBY = 1.;     double lPSY = 1.;
173 <  double lPBZ = 1.;     double lPSZ = 1.;
174 <  if(!overflow(lTPionEm, lP,lEM))  lPBX =  lTPionEm ->GetBinContent(lTPionEm ->GetXaxis()->FindBin(lP),lTPionEm ->GetYaxis()->FindBin(lEM) );
175 <  if(!overflow(lTPionHad,lP,lHad)) lPBY =  lTPionHad->GetBinContent(lTPionHad->GetXaxis()->FindBin(lP),lTPionHad->GetYaxis()->FindBin(lHad));
176 <  if(!overflow(lTPionHo, lP,lHO))  lPBZ =  lTPionHo ->GetBinContent(lTPionHo ->GetXaxis()->FindBin(lP),lTPionHo ->GetYaxis()->FindBin(lHO) );
177 <  if(!overflow(lTMuonEm, lP,lEM )) lPSX =  lTMuonEm ->GetBinContent(lTMuonEm ->GetXaxis()->FindBin(lP),lTMuonEm ->GetYaxis()->FindBin(lEM) );
178 <  if(!overflow(lTMuonHad,lP,lHad)) lPSY =  lTMuonHad->GetBinContent(lTMuonHad->GetXaxis()->FindBin(lP),lTMuonHad->GetYaxis()->FindBin(lHad));
179 <  if(!overflow(lTMuonHo ,lP,lHO))  lPSZ =  lTMuonHo ->GetBinContent(lTMuonHo ->GetXaxis()->FindBin(lP),lTMuonHo ->GetYaxis()->FindBin(lHO) );
180 <  
181 <  if(lPSX == 0. || lPBX == 0. || (lEM <= 0. && !iEMSpecial)) {lPSX = 1.; lPBX = 1.;}
182 <  if(lPSY == 0. || lPBY == 0. || lHad == 0.) {lPSY = 1.; lPBY = 1.;}
183 <  if(lPSZ == 0. || lPBZ == 0. || lHO  == 0.) {lPSZ = 1.; lPBZ = 1.;}
184 <  if((lPSX*lPSY*lPSZ+lPBX*lPBY*lPBZ) > 0.) return lPSX*lPSY*lPSZ/(lPSX*lPSY*lPSZ+lPBX*lPBY*lPBZ);
171 >  Double_t lPBX = 1.;    
172 >  Double_t lPSX = 1.;
173 >  Double_t lPBY = 1.;    
174 >  Double_t lPSY = 1.;
175 >  Double_t lPBZ = 1.;    
176 >  Double_t lPSZ = 1.;
177 >  if (!Overflow(lTPionEm, lP,lEM))  
178 >    lPBX = lTPionEm ->GetBinContent(lTPionEm ->GetXaxis()->FindBin(lP),
179 >                                     lTPionEm ->GetYaxis()->FindBin(lEM));
180 >  if (!Overflow(lTPionHad,lP,lHad))
181 >    lPBY = lTPionHad->GetBinContent(lTPionHad->GetXaxis()->FindBin(lP),
182 >                                     lTPionHad->GetYaxis()->FindBin(lHad));
183 >  if (!Overflow(lTPionHo, lP,lHO))  
184 >    lPBZ = lTPionHo ->GetBinContent(lTPionHo ->GetXaxis()->FindBin(lP),
185 >                                     lTPionHo ->GetYaxis()->FindBin(lHO));
186 >  if (!Overflow(lTMuonEm, lP,lEM ))
187 >    lPSX = lTMuonEm ->GetBinContent(lTMuonEm ->GetXaxis()->FindBin(lP),
188 >                                     lTMuonEm ->GetYaxis()->FindBin(lEM));
189 >  if (!Overflow(lTMuonHad,lP,lHad))
190 >    lPSY = lTMuonHad->GetBinContent(lTMuonHad->GetXaxis()->FindBin(lP),
191 >                                    lTMuonHad->GetYaxis()->FindBin(lHad));
192 >  if (!Overflow(lTMuonHo ,lP,lHO))  
193 >    lPSZ =  lTMuonHo ->GetBinContent(lTMuonHo ->GetXaxis()->FindBin(lP),
194 >                                     lTMuonHo ->GetYaxis()->FindBin(lHO));
195 >  
196 >  if (lPSX == 0. || lPBX == 0. || (lEM <= 0. && !iEMSpecial)) {
197 >    lPSX = 1.;
198 >    lPBX = 1.;
199 >  }
200 >  if (lPSY == 0. || lPBY == 0. || lHad == 0.) {
201 >    lPSY = 1.;
202 >    lPBY = 1.;
203 >  }
204 >  if (lPSZ == 0. || lPBZ == 0. || lHO  == 0.) {
205 >    lPSZ = 1.;
206 >    lPBZ = 1.;
207 >  }
208 >  if ((lPSX*lPSY*lPSZ+lPBX*lPBY*lPBZ) > 0.)
209 >    return lPSX*lPSY*lPSZ/(lPSX*lPSY*lPSZ+lPBX*lPBY*lPBZ);
210 >
211    return 0.5;
212   }
213  
214 < bool MuonTools::isGood(mithep::Muon *iMuon,selection iSelection) {
215 <  double lVal = 0;
216 <  switch(iSelection) {
217 <  case AllArbitrated:
218 <    if(iMuon->StandaloneTrk() != 0 || iMuon->GlobalTrk()!= 0)  return true;
219 <    if(iMuon->NSegments() > 0) return true;
220 <    return false;
221 <    break;
222 <  case PromptTight:
223 <    return promptTight(-1,iMuon);
224 <    break;
225 <   case TMOneStationLoose:
226 <     return TMOneStation(iMuon,99999,999999);
227 <     break;
228 <  case TMOneStationTight:
229 <    return TMOneStation(iMuon);
230 <    break;
231 <  case TMLastStationLoose:
232 <    return TMLastStation(iMuon,999999,999999);
233 <    break;
234 <  case TMLastStationTight:
235 <    return TMLastStation(iMuon);
236 <    break;
237 <  case TM2DCompatibilityLoose:
238 <    lVal             = 1.2*getSegmentCompatability(iMuon);
239 <    if(lVal/1.2 == 0.5) return false;
240 <    lVal += 0.8*getCaloCompatability(iMuon,true,true);
241 <    if(lVal > 0.7) return true;
242 <    return false;
243 <    break;
244 <  case TM2DCompatibilityTight:
245 <    lVal             = 1.2*getSegmentCompatability(iMuon);
196 <    if(lVal/1.2 == 0.5) return false;
197 <    lVal += 0.8*getCaloCompatability(iMuon,true,true);
198 <    if(lVal > 1.0) return true;
199 <    return false;
200 <    break;
201 <  default:
202 <    return false;
214 > //--------------------------------------------------------------------------------------------------
215 > Bool_t MuonTools::Init(const char *mutemp, const char *pitemp)
216 > {
217 >  // Read histograms from given files.
218 >
219 >  if (fIsInit) {
220 >    DeleteHistos();
221 >  }
222 >
223 >  TDirectory::TContext context(0);
224 >
225 >  TFile *muon_templates = TFile::Open(mutemp);
226 >  if (!muon_templates) {
227 >    Fatal("Init", "Could not open file %s", mutemp);
228 >    return kFALSE;
229 >  }
230 >  fmuon_em_etaEmi  = LoadHisto("em_etaEmi",  muon_templates);
231 >  fmuon_had_etaEmi = LoadHisto("had_etaEmi", muon_templates);
232 >  fmuon_had_etaTmi = LoadHisto("had_etaTmi", muon_templates);
233 >  fmuon_em_etaB    = LoadHisto("em_etaB",    muon_templates);
234 >  fmuon_had_etaB   = LoadHisto("had_etaB",   muon_templates);
235 >  fmuon_ho_etaB    = LoadHisto("ho_etaB",    muon_templates);
236 >  fmuon_had_etaTpl = LoadHisto("had_etaTpl", muon_templates);
237 >  fmuon_em_etaEpl  = LoadHisto("em_etaEpl",  muon_templates);
238 >  fmuon_had_etaEpl = LoadHisto("had_etaEpl", muon_templates);
239 >  muon_templates->Close();
240 >  delete muon_templates;
241 >
242 >  TFile *pion_templates = TFile::Open(pitemp);
243 >  if (!pion_templates) {
244 >    Fatal("Init", "Could not open file %s", pitemp);
245 >    return kFALSE;
246    }
247 <  return false;
247 >
248 >  fpion_em_etaEmi  = LoadHisto("em_etaEmi",  pion_templates);
249 >  fpion_had_etaEmi = LoadHisto("had_etaEmi", pion_templates);
250 >  fpion_had_etaTmi = LoadHisto("had_etaTmi", pion_templates);
251 >  fpion_em_etaB    = LoadHisto("em_etaB",    pion_templates);
252 >  fpion_had_etaB   = LoadHisto("had_etaB",   pion_templates);
253 >  fpion_ho_etaB    = LoadHisto("ho_etaB",    pion_templates);
254 >  fpion_had_etaTpl = LoadHisto("had_etaTpl", pion_templates);
255 >  fpion_em_etaEpl  = LoadHisto("em_etaEpl",  pion_templates);
256 >  fpion_had_etaEpl = LoadHisto("had_etaEpl", pion_templates);
257 >  pion_templates->Close();
258 >  delete pion_templates;
259 >
260 >  fIsInit = kTRUE;
261 >  return kTRUE;
262 > }
263 >
264 > //--------------------------------------------------------------------------------------------------
265 > Bool_t MuonTools::IsGood(const mithep::Muon *iMuon, ESelType iSel) const
266 > {
267 >  // Return true if given muon qualifies given selection criterium.
268 >
269 >  Double_t tm2dcut = 0.;
270 >
271 >  switch(iSel) {
272 >    case kAllArbitrated:
273 >      if (iMuon->StandaloneTrk() != 0 || iMuon->GlobalTrk()!= 0)  
274 >        return kTRUE;
275 >      if (iMuon->NSegments() > 0)
276 >        return kTRUE;
277 >      break;
278 >    case kPromptTight:
279 >      return iMuon->PromptTight(Muon::kAny);
280 >      break;
281 >    case kTMOneStationLoose:
282 >      return iMuon->TMOneStation(999999,999999);
283 >      break;
284 >    case kTMOneStationTight:
285 >      return iMuon->TMOneStation();
286 >      break;
287 >    case kTMLastStationLoose:
288 >      return iMuon->TMLastStation(999999,999999);
289 >      break;
290 >    case kTMLastStationTight:
291 >      return iMuon->TMLastStation();
292 >      break;
293 >    case kTM2DCompatibilityLoose:
294 >      tm2dcut = 0.7;
295 >      break;
296 >    case kTM2DCompatibilityTight:
297 >      tm2dcut = 1.0;
298 >      break;
299 >    default:
300 >      return kFALSE;
301 >      break;
302 >  }
303 >
304 >  Double_t lVal = GetSegmentCompatability(iMuon);
305 >  if (lVal == 0.5) // exclude this border case
306 >    return kFALSE;
307 >
308 >  lVal *= 1.2;
309 >  lVal += 0.8*GetCaloCompatability(iMuon,kTRUE,kTRUE);
310 >  if (lVal > tm2dcut)
311 >    return kTRUE;
312 >
313 >  return kFALSE;
314 > }
315 >
316 > //--------------------------------------------------------------------------------------------------
317 > Double_t MuonTools::GetSegmentCompatability(const mithep::Muon *iMuon) const
318 > {
319 >  // Get segment compatability for given muon based on likelihood of well defined
320 >  // track through chambers.
321 >
322 >  Int_t lNStationsCrossed = 0;
323 >  Int_t lNStationsSegment = 0;
324 >  
325 >  Int_t lStSegmentmatch[8];
326 >  Int_t lStCrossed[8];
327 >  Double_t lStBoundary[8];
328 >
329 >  Double_t lWeight  = 0.;
330 >  Bool_t lAdjust    = kTRUE;
331 >  for (Int_t i0 = 0; i0 < 8; ++i0) {
332 >    lStBoundary[i0] = 0.;
333 >    if(iMuon->GetTrackDist(i0) < 999999. ) {
334 >      lNStationsCrossed++;
335 >      lStCrossed[i0]  = 1;
336 >      if (iMuon->GetTrackDist(i0) > -10. )
337 >        lStBoundary[i0] = iMuon->GetTrackDist(i0);
338 >    } else
339 >      lStCrossed[i0]  = 0;
340 >
341 >    if(iMuon->GetDX(i0) < 999999.) {
342 >      lNStationsSegment++;
343 >      lStSegmentmatch[i0] = 1;
344 >    } else
345 >      lStSegmentmatch[i0] = 0;
346 >
347 >    if(iMuon->GetDY(i0) < 999999.)
348 >      lAdjust = kFALSE;
349 >  }
350 >
351 >  if (lNStationsCrossed == 0)
352 >    return 0.5;
353 >
354 >  Double_t lStWeight[8];
355 >  Int_t lPCross = -1;
356 >  const Double_t lAtWeight = 0.5;
357 >  for (Int_t i0 = 0; i0< 8; ++i0) {
358 >    lStWeight[i0] = 0;
359 >    if (lStCrossed[i0] > 0) {
360 >      lPCross++;
361 >
362 >      switch (lNStationsCrossed) {
363 >        case 1 :
364 >          lStWeight[i0] =  1.;
365 >          break;
366 >        case 2 :
367 >          if (lPCross == 0 )
368 >            lStWeight[i0] = 0.33;
369 >          else
370 >            lStWeight[i0] = 0.67;
371 >          break;
372 >        case 3 :
373 >          if (lPCross == 0)
374 >            lStWeight[i0] = 0.23;
375 >          else if (lPCross == 1)
376 >            lStWeight[i0] = 0.33;
377 >          else
378 >            lStWeight[i0] = 0.44;
379 >          break;
380 >        case 4 :
381 >          if (lPCross == 0)
382 >            lStWeight[i0] = 0.10;
383 >          else if (lPCross == 1)
384 >            lStWeight[i0] = 0.20;
385 >          else if (lPCross == 2)
386 >            lStWeight[i0] = 0.30;
387 >          else                    
388 >            lStWeight[i0] = 0.40;
389 >          break;
390 >        default :
391 >          lStWeight[i0] = 1./lNStationsCrossed;
392 >      }
393 >    
394 >      if (lStSegmentmatch[i0] <= 0 && lStBoundary[i0] != 0.)
395 >        lStWeight[i0] *= lAtWeight*0.5*(TMath::Erf(lStBoundary[i0]/6.)+1.);
396 >      else if (lStSegmentmatch[i0] <= 0 && lStBoundary[i0] == 0)
397 >        lStWeight[i0] = 0.;
398 >      
399 >      if (lStSegmentmatch[i0] > 0) {
400 >        Double_t lP2X = TMath::Power(iMuon->GetPullX(i0),2.);
401 >        Double_t lP2Y = TMath::Power(iMuon->GetPullY(i0),2.);
402 >        Double_t lD2X = TMath::Power(iMuon->GetDX(i0),2.);
403 >        Double_t lD2Y = TMath::Power(iMuon->GetDY(i0),2.);
404 >        if (iMuon->GetDY(i0) < 999999 && iMuon->GetDX(i0) < 999999)
405 >          lStWeight[i0] *= SigWeight(TMath::Sqrt(lD2X+lD2Y),TMath::Sqrt(lP2X+lP2Y));
406 >        else if (iMuon->GetDY(i0) >= 999999 && i0 < 4)
407 >          lStWeight[i0] *= SigWeight(iMuon->GetDX(i0),iMuon->GetPullX(i0));
408 >        else if(i0 < 4)
409 >          lStWeight[i0] *= SigWeight(iMuon->GetDY(i0),iMuon->GetPullY(i0));
410 >      }
411 >    }
412 >    lWeight += lStWeight[i0];
413 >  }
414 >
415 >  return lWeight;
416 > }
417 >
418 > //--------------------------------------------------------------------------------------------------
419 > TH2D *MuonTools::LoadHisto(const char *name, TFile *file) const
420 > {
421 >  // Load histogram with given name from given file and return it.
422 >
423 >  TH2D *ret = dynamic_cast<TH2D*>(file->Get(name));
424 >  if (!ret) {
425 >    Fatal("LoadHisto", "Could not load histogram %s from file %s", name, file->GetName());
426 >    return 0;
427 >  }
428 >  ret->SetDirectory(0);
429 >  return ret;
430 > }
431 > //--------------------------------------------------------------------------------------------------
432 > Bool_t MuonTools::PassD0Cut(const Muon *mu, const VertexCol *vertices, Double_t fD0Cut)
433 > {
434 >  Bool_t d0cut = kFALSE;
435 >  const Track *mt = mu->BestTrk();
436 >  if (!mt) return kFALSE;
437 >
438 >  Double_t d0_real = 1e30;
439 >  for(UInt_t i0 = 0; i0 < vertices->GetEntries(); i0++) {
440 >    if(vertices->At(i0)->NTracks() > 0){
441 >      Double_t pD0 = mt->D0Corrected(*vertices->At(i0));
442 >      d0_real = TMath::Abs(pD0);
443 >      break;
444 >    }
445 >  }
446 >  if(d0_real < fD0Cut) d0cut = kTRUE;
447 >  
448 >  return d0cut;
449 > }
450 >
451 > //--------------------------------------------------------------------------------------------------
452 > Bool_t MuonTools::PassD0Cut(const Muon *mu, const BeamSpotCol *beamspots, Double_t fD0Cut)
453 > {
454 >  Bool_t d0cut = kFALSE;
455 >  const Track *mt = mu->BestTrk();
456 >  if (!mt) return kFALSE;
457 >
458 >  // d0 cut
459 >  Double_t d0_real = 99999;
460 >  for(UInt_t i0 = 0; i0 < beamspots->GetEntries(); i0++) {
461 >    Double_t pD0 = mt->D0Corrected(*beamspots->At(i0));
462 >    if(TMath::Abs(pD0) < TMath::Abs(d0_real)) d0_real = TMath::Abs(pD0);
463 >  }
464 >  if(d0_real < fD0Cut) d0cut = kTRUE;
465 >  
466 >  return d0cut;
467 > }
468 >
469 > //--------------------------------------------------------------------------------------------------
470 > Bool_t MuonTools::PassSoftMuonCut(const Muon *mu, const VertexCol *vertices)
471 > {
472 >  if(mu->Pt() <= 3.0) return kFALSE;
473 >  
474 >  if(!mu->IsTrackerMuon()) return kFALSE;
475 >  
476 >  if(!mu->Quality().Quality(MuonQuality::TMLastStationAngTight)) return kFALSE;
477 >  
478 >  if(mu->BestTrk()->NHits() <= 10) return kFALSE;
479 >
480 >  if(!PassD0Cut(mu, vertices, 0.2)) return kFALSE;
481 >  
482 >  Double_t totalIso = 1.0 * mu->IsoR03SumPt() +
483 >                      1.0 * mu->IsoR03EmEt() +
484 >                      1.0 * mu->IsoR03HadEt();
485 >  if (totalIso < (mu->Pt()*0.10) && mu->Pt() > 20.0) return kFALSE;
486 >
487 >  return kTRUE;
488   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines