ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MuonTools.cc
Revision: 1.6
Committed: Fri Nov 28 09:13:35 2008 UTC (16 years, 5 months ago) by loizides
Content type: text/plain
Branch: MAIN
Changes since 1.5: +2 -2 lines
Log Message:
More comments, rest to be done by Phil

File Contents

# User Rev Content
1 loizides 1.6 // $Id: MuonTools.cc,v 1.5 2008/11/27 16:28:58 loizides Exp $
2 loizides 1.5
3 pharris 1.2 #include "MitPhysics/Utils/interface/MuonTools.h"
4 loizides 1.5 #include <TFile.h>
5    
6 pharris 1.1 using namespace mithep;
7    
8 loizides 1.5 //--------------------------------------------------------------------------------------------------
9     MuonTools::MuonTools(const char *mutemp, const char *pitemp) :
10     fIsInit(kFALSE),
11     fmuon_em_etaEmi(0),
12     fmuon_had_etaEmi(0),
13     fmuon_had_etaTmi(0),
14     fmuon_em_etaB(0),
15     fmuon_had_etaB(0),
16     fmuon_ho_etaB(0),
17     fmuon_had_etaTpl(0),
18     fmuon_em_etaEpl(0),
19     fmuon_had_etaEpl(0),
20 pharris 1.1 fpion_em_etaEmi(0),
21     fpion_had_etaEmi(0),
22     fpion_had_etaTmi(0),
23     fpion_em_etaB(0),
24     fpion_had_etaB(0),
25     fpion_ho_etaB(0),
26     fpion_had_etaTpl(0),
27     fpion_em_etaEpl(0),
28 loizides 1.5 fpion_had_etaEpl(0)
29     {
30     // Constructor.
31    
32     if (mutemp && pitemp)
33     Init(mutemp, pitemp);
34     }
35    
36     //--------------------------------------------------------------------------------------------------
37     MuonTools::~MuonTools()
38     {
39     // Destructor.
40    
41     DeleteHistos();
42     }
43    
44     //--------------------------------------------------------------------------------------------------
45     void MuonTools::DeleteHistos()
46     {
47     // Delete histograms.
48    
49     if (fIsInit) {
50     delete fpion_em_etaEmi;
51     delete fpion_had_etaEmi;
52     delete fpion_had_etaTmi;
53     delete fpion_em_etaB;
54     delete fpion_had_etaB;
55     delete fpion_ho_etaB;
56     delete fpion_had_etaTpl;
57     delete fpion_em_etaEpl;
58     delete fpion_had_etaEpl;
59     delete fmuon_em_etaEmi;
60     delete fmuon_had_etaEmi;
61     delete fmuon_had_etaTmi;
62     delete fmuon_em_etaB;
63     delete fmuon_had_etaB;
64     delete fmuon_ho_etaB;
65     delete fmuon_had_etaTpl;
66     delete fmuon_em_etaEpl;
67     delete fmuon_had_etaEpl;
68     fpion_em_etaEmi = 0;
69     fpion_had_etaEmi = 0;
70     fpion_had_etaTmi = 0;
71     fpion_em_etaB = 0;
72     fpion_had_etaB = 0;
73     fpion_ho_etaB = 0;
74     fpion_had_etaTpl = 0;
75     fpion_em_etaEpl = 0;
76     fpion_had_etaEpl = 0;
77     fmuon_em_etaEmi = 0;
78     fmuon_had_etaEmi = 0;
79     fmuon_had_etaTmi = 0;
80     fmuon_em_etaB = 0;
81     fmuon_had_etaB = 0;
82     fmuon_ho_etaB = 0;
83     fmuon_had_etaTpl = 0;
84     fmuon_em_etaEpl = 0;
85     fmuon_had_etaEpl = 0;
86     fIsInit = kFALSE;
87     }
88     }
89    
90     //--------------------------------------------------------------------------------------------------
91     Double_t MuonTools::GetCaloCompatability(const Muon *iMuon,
92     Bool_t iEMSpecial, Bool_t iCorrectedHCAL) const
93     {
94     // todo
95    
96     Double_t lEta = iMuon->Eta();
97     Double_t aEta = TMath::Abs(lEta);
98     if (aEta > 2.5)
99     return 0.5;
100    
101     Double_t lP = iMuon->P();
102     if (lP >= 2000.)
103     lP = 1999.9;
104     if(lP < 0. )
105     return 0.5;
106    
107     Double_t lEM = -5.;
108     if (!iEMSpecial || iMuon->EmEnergy() != 0.)
109     lEM = iMuon->EmEnergy();
110    
111     Double_t lHad = iMuon->HadEnergy();
112     Double_t lHO = iMuon->HoEnergy();;
113    
114     TH2D *lTMuonHad = 0;
115     TH2D *lTPionHad = 0;
116     TH2D *lTMuonHo = 0;
117     TH2D *lTPionHo = 0;
118     TH2D *lTMuonEm = 0;
119     TH2D *lTPionEm = 0;
120    
121     if (aEta >= 1.27) {
122     if (iCorrectedHCAL)
123     lHad *= 1.8/2.2;
124     if (lEta > 0) {
125 pharris 1.1 lTPionHad = fpion_had_etaEpl;
126     lTMuonHad = fmuon_had_etaEpl;
127     } else {
128     lTPionHad = fpion_had_etaEmi;
129     lTMuonHad = fmuon_had_etaEmi;
130     }
131     }
132 loizides 1.5
133     if (aEta < 1.27 && aEta >= 1.1) {
134     if (iCorrectedHCAL)
135     lHad *= (1.8/(-2.2*aEta+5.5));
136     if (lEta > 0) {
137 pharris 1.1 lTPionHad = fpion_had_etaTpl;
138     lTMuonHad = fmuon_had_etaTpl;
139     } else {
140     lTPionHad = fpion_had_etaTmi;
141     lTMuonHad = fmuon_had_etaTmi;
142     }
143     }
144 loizides 1.5
145     if (aEta < 1.1) {
146     if(iCorrectedHCAL)
147     lHad *= TMath::Sin(2*TMath::ATan(TMath::Exp(lEta))); //todo ask!
148 pharris 1.1 lTPionHad = fpion_had_etaB;
149     lTMuonHad = fmuon_had_etaB;
150     }
151 loizides 1.5 if (lEta > 1.479) {
152     lTPionEm = fpion_em_etaEpl;
153     lTMuonEm = fmuon_em_etaEpl;
154 pharris 1.1 }
155 loizides 1.5 if (aEta <= 1.479) {
156 pharris 1.1 lTPionEm = fpion_em_etaB;
157     lTMuonEm = fmuon_em_etaB;
158     }
159 loizides 1.5 if (lEta < -1.479) {
160 pharris 1.1 lTPionEm = fpion_em_etaEmi;
161     lTMuonEm = fmuon_em_etaEmi;
162     }
163 loizides 1.5 if (aEta < 1.28) {
164 pharris 1.1 lTPionHo = fpion_ho_etaB;
165     lTMuonHo = fmuon_ho_etaB;
166     }
167    
168 loizides 1.5 Double_t lPBX = 1.;
169     Double_t lPSX = 1.;
170     Double_t lPBY = 1.;
171     Double_t lPSY = 1.;
172     Double_t lPBZ = 1.;
173     Double_t lPSZ = 1.;
174     if (!Overflow(lTPionEm, lP,lEM))
175     lPBX = lTPionEm ->GetBinContent(lTPionEm ->GetXaxis()->FindBin(lP),
176     lTPionEm ->GetYaxis()->FindBin(lEM));
177     if (!Overflow(lTPionHad,lP,lHad))
178     lPBY = lTPionHad->GetBinContent(lTPionHad->GetXaxis()->FindBin(lP),
179     lTPionHad->GetYaxis()->FindBin(lHad));
180     if (!Overflow(lTPionHo, lP,lHO))
181     lPBZ = lTPionHo ->GetBinContent(lTPionHo ->GetXaxis()->FindBin(lP),
182     lTPionHo ->GetYaxis()->FindBin(lHO));
183     if (!Overflow(lTMuonEm, lP,lEM ))
184     lPSX = lTMuonEm ->GetBinContent(lTMuonEm ->GetXaxis()->FindBin(lP),
185     lTMuonEm ->GetYaxis()->FindBin(lEM));
186     if (!Overflow(lTMuonHad,lP,lHad))
187     lPSY = lTMuonHad->GetBinContent(lTMuonHad->GetXaxis()->FindBin(lP),
188     lTMuonHad->GetYaxis()->FindBin(lHad));
189     if (!Overflow(lTMuonHo ,lP,lHO))
190     lPSZ = lTMuonHo ->GetBinContent(lTMuonHo ->GetXaxis()->FindBin(lP),
191     lTMuonHo ->GetYaxis()->FindBin(lHO));
192 pharris 1.1
193 loizides 1.5 if (lPSX == 0. || lPBX == 0. || (lEM <= 0. && !iEMSpecial)) {
194     lPSX = 1.;
195     lPBX = 1.;
196     }
197     if (lPSY == 0. || lPBY == 0. || lHad == 0.) {
198     lPSY = 1.;
199     lPBY = 1.;
200     }
201     if (lPSZ == 0. || lPBZ == 0. || lHO == 0.) {
202     lPSZ = 1.;
203     lPBZ = 1.;
204     }
205     if ((lPSX*lPSY*lPSZ+lPBX*lPBY*lPBZ) > 0.)
206     return lPSX*lPSY*lPSZ/(lPSX*lPSY*lPSZ+lPBX*lPBY*lPBZ);
207    
208 pharris 1.1 return 0.5;
209     }
210    
211 loizides 1.5 //--------------------------------------------------------------------------------------------------
212     Bool_t MuonTools::Init(const char *mutemp, const char *pitemp)
213     {
214     // Read histograms from given files.
215    
216     if (fIsInit) {
217     DeleteHistos();
218     }
219    
220     TDirectory::TContext context(0);
221    
222     TFile *muon_templates = TFile::Open(mutemp);
223     if (!muon_templates) {
224     Fatal("Init", "Could not open file %s", mutemp);
225     return kFALSE;
226     }
227     fmuon_em_etaEmi = LoadHisto("em_etaEmi", muon_templates);
228     fmuon_had_etaEmi = LoadHisto("had_etaEmi", muon_templates);
229     fmuon_had_etaTmi = LoadHisto("had_etaTmi", muon_templates);
230     fmuon_em_etaB = LoadHisto("em_etaB", muon_templates);
231     fmuon_had_etaB = LoadHisto("had_etaB", muon_templates);
232     fmuon_ho_etaB = LoadHisto("ho_etaB", muon_templates);
233     fmuon_had_etaTpl = LoadHisto("had_etaTpl", muon_templates);
234     fmuon_em_etaEpl = LoadHisto("em_etaEpl", muon_templates);
235     fmuon_had_etaEpl = LoadHisto("had_etaEpl", muon_templates);
236     muon_templates->Close();
237     delete muon_templates;
238    
239     TFile *pion_templates = TFile::Open(pitemp);
240     if (!pion_templates) {
241     Fatal("Init", "Could not open file %s", pitemp);
242     return kFALSE;
243     }
244    
245     fpion_em_etaEmi = LoadHisto("em_etaEmi", pion_templates);
246     fpion_had_etaEmi = LoadHisto("had_etaEmi", pion_templates);
247     fpion_had_etaTmi = LoadHisto("had_etaTmi", pion_templates);
248     fpion_em_etaB = LoadHisto("em_etaB", pion_templates);
249     fpion_had_etaB = LoadHisto("had_etaB", pion_templates);
250     fpion_ho_etaB = LoadHisto("ho_etaB", pion_templates);
251     fpion_had_etaTpl = LoadHisto("had_etaTpl", pion_templates);
252     fpion_em_etaEpl = LoadHisto("em_etaEpl", pion_templates);
253     fpion_had_etaEpl = LoadHisto("had_etaEpl", pion_templates);
254     pion_templates->Close();
255     delete pion_templates;
256    
257     fIsInit = kTRUE;
258     return kTRUE;
259     }
260    
261     //--------------------------------------------------------------------------------------------------
262     Bool_t MuonTools::IsGood(const mithep::Muon *iMuon, ESelType iSel) const
263     {
264     // Return true if given muon qualifies given selection criterium.
265    
266     Double_t tm2dcut = 0.;
267    
268     switch(iSel) {
269     case kAllArbitrated:
270     if (iMuon->StandaloneTrk() != 0 || iMuon->GlobalTrk()!= 0)
271     return kTRUE;
272     if (iMuon->NSegments() > 0)
273     return kTRUE;
274     break;
275     case kPromptTight:
276     return iMuon->PromptTight(Muon::kAny);
277     break;
278     case kTMOneStationLoose:
279     return iMuon->TMOneStation(99999,999999);
280     break;
281     case kTMOneStationTight:
282     return iMuon->TMOneStation();
283     break;
284     case kTMLastStationLoose:
285     return iMuon->TMLastStation(999999,999999);
286     break;
287     case kTMLastStationTight:
288     return iMuon->TMLastStation();
289     break;
290     case kTM2DCompatibilityLoose:
291     tm2dcut = 0.7;
292     break;
293     case kTM2DCompatibilityTight:
294     tm2dcut = 1.0;
295     break;
296     default:
297     return kFALSE;
298     break;
299     }
300    
301     Double_t lVal = 1.2*GetSegmentCompatability(iMuon);
302     if (lVal/1.2 == 0.5)
303     return kFALSE;
304    
305     lVal += 0.8*GetCaloCompatability(iMuon,kTRUE,kTRUE);
306     if (lVal > tm2dcut)
307     return kTRUE;
308     return kFALSE;
309     }
310    
311     //--------------------------------------------------------------------------------------------------
312     Double_t MuonTools::GetSegmentCompatability(const mithep::Muon *iMuon) const
313     {
314 loizides 1.6 // todo
315 loizides 1.5
316     Int_t lNStationsCrossed = 0;
317     Int_t lNStationsSegment = 0;
318    
319     Int_t lStSegmentmatch[8];
320     Int_t lStCrossed[8];
321     Double_t lStBoundary[8];
322    
323     Double_t lWeight = 0.;
324     Bool_t lAdjust = kTRUE;
325     for (Int_t i0 = 0; i0 < 8; ++i0) {
326     lStBoundary[i0] = 0.;
327     if(iMuon->GetTrackDist(i0) < 999999. ) {
328     lNStationsCrossed++;
329     lStCrossed[i0] = 1;
330     if (iMuon->GetTrackDist(i0) > -10. )
331     lStBoundary[i0] = iMuon->GetTrackDist(i0);
332     } else
333     lStCrossed[i0] = 0;
334    
335     if(iMuon->GetDX(i0) < 999999.) { //Use iMuon->GetSegmentX--> CHECK
336     lNStationsSegment++;
337     lStSegmentmatch[i0] = 1;
338     } else
339     lStSegmentmatch[i0] = 0;
340    
341     if(iMuon->GetDY(i0) < 999999.)
342     lAdjust = kFALSE;
343     }
344    
345     if (lNStationsCrossed == 0)
346     return 0.5;
347    
348     Double_t lStWeight[8];
349     Int_t lPCross = -1;
350     const Double_t lAtWeight = 0.5;
351     for (Int_t i0 = 0; i0< 8; ++i0) {
352     lStWeight[i0] = 0;
353     if (lStCrossed[i0] > 0) {
354     lPCross++;
355    
356     switch (lNStationsCrossed) {
357     case 1 :
358     lStWeight[i0] = 1.;
359     break;
360     case 2 :
361     if (lPCross == 0 )
362     lStWeight[i0] = 0.33;
363     else
364     lStWeight[i0] = 0.67;
365     break;
366     case 3 :
367     if (lPCross == 0)
368     lStWeight[i0] = 0.23;
369     else if (lPCross == 1)
370     lStWeight[i0] = 0.33;
371     else
372     lStWeight[i0] = 0.44;
373     break;
374     case 4 :
375     if (lPCross == 0)
376     lStWeight[i0] = 0.10;
377     else if (lPCross == 1)
378     lStWeight[i0] = 0.20;
379     else if (lPCross == 2)
380     lStWeight[i0] = 0.30;
381     else
382     lStWeight[i0] = 0.40;
383     break;
384     default :
385     lStWeight[i0] = 1./lNStationsCrossed;
386     }
387    
388     if (lStSegmentmatch[i0] <= 0 && lStBoundary[i0] != 0.)
389     lStWeight[i0] *= lAtWeight*0.5*(TMath::Erf(lStBoundary[i0]/6.)+1.);
390     else if (lStSegmentmatch[i0] <= 0 && lStBoundary[i0] == 0)
391     lStWeight[i0] = 0.;
392    
393     if (lStSegmentmatch[i0] > 0) {
394     Double_t lP2X = TMath::Power(iMuon->GetPullX(i0),2.);
395     Double_t lP2Y = TMath::Power(iMuon->GetPullY(i0),2.);
396     Double_t lD2X = TMath::Power(iMuon->GetDX(i0),2.);
397     Double_t lD2Y = TMath::Power(iMuon->GetDY(i0),2.);
398     if (iMuon->GetDY(i0) < 999999 && iMuon->GetDX(i0) < 999999)
399     lStWeight[i0] *= SigWeight(TMath::Sqrt(lD2X+lD2Y),TMath::Sqrt(lP2X+lP2Y));
400     else if (iMuon->GetDY(i0) >= 999999 && i0 < 4)
401     lStWeight[i0] *= SigWeight(iMuon->GetDX(i0),iMuon->GetPullX(i0));
402     else if(i0 < 4)
403     lStWeight[i0] *= SigWeight(iMuon->GetDY(i0),iMuon->GetPullY(i0));
404     }
405     }
406     lWeight += lStWeight[i0];
407     }
408    
409     return lWeight;
410     }
411    
412     //--------------------------------------------------------------------------------------------------
413     TH2D *MuonTools::LoadHisto(const char *name, TFile *file) const
414     {
415     // Load histogram with given name from given file and return it.
416    
417     TH2D *ret = dynamic_cast<TH2D*>(file->Get(name));
418     if (!ret) {
419     Fatal("LoadHisto", "Could not load histogram %s from file %s", name, file->GetName());
420     return 0;
421 pharris 1.1 }
422 loizides 1.5 ret->SetDirectory(0);
423     return ret;
424 pharris 1.1 }
425 loizides 1.5