Utils/src/MuonTools.cc

// $Id: MuonTools.cc,v 1.8 2009/03/23 14:48:08 loizides Exp $

#include "MitPhysics/Utils/interface/MuonTools.h"
#include <TFile.h>

using namespace mithep;

//--------------------------------------------------------------------------------------------------
MuonTools::MuonTools(const char *mutemp, const char *pitemp) :
  fIsInit(kFALSE),
  fmuon_em_etaEmi(0),
  fmuon_had_etaEmi(0),
  fmuon_had_etaTmi(0),
  fmuon_em_etaB(0),
  fmuon_had_etaB(0),
  fmuon_ho_etaB(0),
  fmuon_had_etaTpl(0),
  fmuon_em_etaEpl(0),
  fmuon_had_etaEpl(0),
  fpion_em_etaEmi(0),
  fpion_had_etaEmi(0),
  fpion_had_etaTmi(0),
  fpion_em_etaB(0),
  fpion_had_etaB(0),
  fpion_ho_etaB(0),
  fpion_had_etaTpl(0),
  fpion_em_etaEpl(0),
  fpion_had_etaEpl(0)
{
  // Constructor.

  if (mutemp && pitemp)
    Init(mutemp, pitemp);
}

//--------------------------------------------------------------------------------------------------
MuonTools::~MuonTools() 
{
  // Destructor.

  DeleteHistos();
}

//--------------------------------------------------------------------------------------------------
void MuonTools::DeleteHistos()
{
  // Delete histograms.

  if (fIsInit) {
    delete fpion_em_etaEmi; 
    delete fpion_had_etaEmi;
    delete fpion_had_etaTmi;
    delete fpion_em_etaB;
    delete fpion_had_etaB;
    delete fpion_ho_etaB;
    delete fpion_had_etaTpl;
    delete fpion_em_etaEpl;
    delete fpion_had_etaEpl;
    delete fmuon_em_etaEmi;
    delete fmuon_had_etaEmi;
    delete fmuon_had_etaTmi;
    delete fmuon_em_etaB;
    delete fmuon_had_etaB;
    delete fmuon_ho_etaB;
    delete fmuon_had_etaTpl;
    delete fmuon_em_etaEpl;
    delete fmuon_had_etaEpl;
    fpion_em_etaEmi  = 0;
    fpion_had_etaEmi = 0;
    fpion_had_etaTmi = 0;
    fpion_em_etaB    = 0;
    fpion_had_etaB   = 0;
    fpion_ho_etaB    = 0;
    fpion_had_etaTpl = 0;
    fpion_em_etaEpl  = 0;
    fpion_had_etaEpl = 0;
    fmuon_em_etaEmi  = 0;
    fmuon_had_etaEmi = 0;
    fmuon_had_etaTmi = 0;
    fmuon_em_etaB    = 0;
    fmuon_had_etaB   = 0;
    fmuon_ho_etaB    = 0;
    fmuon_had_etaTpl = 0;
    fmuon_em_etaEpl  = 0;
    fmuon_had_etaEpl = 0;
    fIsInit = kFALSE;
  }
}

//--------------------------------------------------------------------------------------------------
Double_t MuonTools::GetCaloCompatability(const Muon *iMuon,
                                         Bool_t iEMSpecial, Bool_t iCorrectedHCAL) const
{
  // Get calo compatibility value for given muon based on calorimeter templates.
  // If iEMSpecial is true, then a use different arrangement of ECAL for compatibility.

  Double_t lEta = iMuon->Eta();
  Double_t aEta = TMath::Abs(lEta);
  if (aEta > 2.5) 
    return 0.5; 

  Double_t lP = iMuon->P();
  if (lP >= 2000.) 
    lP = 1999.9;
  if(lP < 0. )           
    return 0.5; 

  Double_t lEM  = -5.;      
  if (!iEMSpecial || iMuon->EmEnergy() != 0.) 
    lEM  = iMuon->EmEnergy();

  Double_t lHad = iMuon->HadEnergy();
  Double_t lHO = iMuon->HoEnergy();;

  TH2D *lTMuonHad = 0;
  TH2D *lTPionHad = 0;
  TH2D *lTMuonHo  = 0;
  TH2D *lTPionHo  = 0;
  TH2D *lTMuonEm  = 0;
  TH2D *lTPionEm  = 0;
    
  if (aEta >= 1.27) {
  if (iCorrectedHCAL) 
      lHad *= 1.8/2.2;
    if (lEta > 0) {
      lTPionHad = fpion_had_etaEpl;
      lTMuonHad = fmuon_had_etaEpl;
    } else {
      lTPionHad = fpion_had_etaEmi;
      lTMuonHad = fmuon_had_etaEmi;
    }
  }

  if (aEta < 1.27 && aEta >= 1.1) {
    if (iCorrectedHCAL)    
      lHad *= (1.8/(-2.2*aEta+5.5));
    if (lEta > 0) {
      lTPionHad  = fpion_had_etaTpl;
      lTMuonHad  = fmuon_had_etaTpl;
    } else {
      lTPionHad  = fpion_had_etaTmi;
      lTMuonHad  = fmuon_had_etaTmi;
    }
  }

  if (aEta < 1.1) {
    if(iCorrectedHCAL)    
      lHad *= TMath::Sin(2*TMath::ATan(TMath::Exp(lEta)));
    lTPionHad  = fpion_had_etaB;
    lTMuonHad  = fmuon_had_etaB;
  }
  if (lEta > 1.479) {
    lTPionEm = fpion_em_etaEpl;
    lTMuonEm = fmuon_em_etaEpl;
  }
  if (aEta <= 1.479) {
    lTPionEm  = fpion_em_etaB;
    lTMuonEm  = fmuon_em_etaB;
  }
  if (lEta < -1.479) {
    lTPionEm  = fpion_em_etaEmi;
    lTMuonEm  = fmuon_em_etaEmi;
  }
  if (aEta < 1.28) {
    lTPionHo  = fpion_ho_etaB;
    lTMuonHo  = fmuon_ho_etaB;
  }
  
  Double_t lPBX = 1.;     
  Double_t lPSX = 1.; 
  Double_t lPBY = 1.;     
  Double_t lPSY = 1.; 
  Double_t lPBZ = 1.;     
  Double_t lPSZ = 1.; 
  if (!Overflow(lTPionEm, lP,lEM))  
    lPBX = lTPionEm ->GetBinContent(lTPionEm ->GetXaxis()->FindBin(lP),
                                     lTPionEm ->GetYaxis()->FindBin(lEM));
  if (!Overflow(lTPionHad,lP,lHad)) 
    lPBY = lTPionHad->GetBinContent(lTPionHad->GetXaxis()->FindBin(lP),
                                     lTPionHad->GetYaxis()->FindBin(lHad));
  if (!Overflow(lTPionHo, lP,lHO))  
    lPBZ = lTPionHo ->GetBinContent(lTPionHo ->GetXaxis()->FindBin(lP),
                                     lTPionHo ->GetYaxis()->FindBin(lHO));
  if (!Overflow(lTMuonEm, lP,lEM )) 
    lPSX = lTMuonEm ->GetBinContent(lTMuonEm ->GetXaxis()->FindBin(lP),
                                     lTMuonEm ->GetYaxis()->FindBin(lEM));
  if (!Overflow(lTMuonHad,lP,lHad)) 
    lPSY = lTMuonHad->GetBinContent(lTMuonHad->GetXaxis()->FindBin(lP),
                                    lTMuonHad->GetYaxis()->FindBin(lHad));
  if (!Overflow(lTMuonHo ,lP,lHO))  
    lPSZ =  lTMuonHo ->GetBinContent(lTMuonHo ->GetXaxis()->FindBin(lP),
                                     lTMuonHo ->GetYaxis()->FindBin(lHO));
  
  if (lPSX == 0. || lPBX == 0. || (lEM <= 0. && !iEMSpecial)) {
    lPSX = 1.; 
    lPBX = 1.;
  } 
  if (lPSY == 0. || lPBY == 0. || lHad == 0.) {
    lPSY = 1.; 
    lPBY = 1.;
  }
  if (lPSZ == 0. || lPBZ == 0. || lHO  == 0.) {
    lPSZ = 1.; 
    lPBZ = 1.;
  }
  if ((lPSX*lPSY*lPSZ+lPBX*lPBY*lPBZ) > 0.) 
    return lPSX*lPSY*lPSZ/(lPSX*lPSY*lPSZ+lPBX*lPBY*lPBZ);

  return 0.5;
}

//--------------------------------------------------------------------------------------------------
Bool_t MuonTools::Init(const char *mutemp, const char *pitemp)
{
  // Read histograms from given files.

  if (fIsInit) {
    DeleteHistos();
  }

  TDirectory::TContext context(0);

  TFile *muon_templates = TFile::Open(mutemp);
  if (!muon_templates) {
    Fatal("Init", "Could not open file %s", mutemp);
    return kFALSE;
  }
  fmuon_em_etaEmi  = LoadHisto("em_etaEmi",  muon_templates);
  fmuon_had_etaEmi = LoadHisto("had_etaEmi", muon_templates);
  fmuon_had_etaTmi = LoadHisto("had_etaTmi", muon_templates);
  fmuon_em_etaB    = LoadHisto("em_etaB",    muon_templates);
  fmuon_had_etaB   = LoadHisto("had_etaB",   muon_templates);
  fmuon_ho_etaB    = LoadHisto("ho_etaB",    muon_templates);
  fmuon_had_etaTpl = LoadHisto("had_etaTpl", muon_templates);
  fmuon_em_etaEpl  = LoadHisto("em_etaEpl",  muon_templates);
  fmuon_had_etaEpl = LoadHisto("had_etaEpl", muon_templates);
  muon_templates->Close();
  delete muon_templates;

  TFile *pion_templates = TFile::Open(pitemp);
  if (!pion_templates) {
    Fatal("Init", "Could not open file %s", pitemp);
    return kFALSE;
  }

  fpion_em_etaEmi  = LoadHisto("em_etaEmi",  pion_templates);
  fpion_had_etaEmi = LoadHisto("had_etaEmi", pion_templates);
  fpion_had_etaTmi = LoadHisto("had_etaTmi", pion_templates);
  fpion_em_etaB    = LoadHisto("em_etaB",    pion_templates);
  fpion_had_etaB   = LoadHisto("had_etaB",   pion_templates);
  fpion_ho_etaB    = LoadHisto("ho_etaB",    pion_templates);
  fpion_had_etaTpl = LoadHisto("had_etaTpl", pion_templates);
  fpion_em_etaEpl  = LoadHisto("em_etaEpl",  pion_templates);
  fpion_had_etaEpl = LoadHisto("had_etaEpl", pion_templates);
  pion_templates->Close();
  delete pion_templates;

  fIsInit = kTRUE;
  return kTRUE;
}

//--------------------------------------------------------------------------------------------------
Bool_t MuonTools::IsGood(const mithep::Muon *iMuon, ESelType iSel) const
{
  // Return true if given muon qualifies given selection criterium.

  Double_t tm2dcut = 0.;

  switch(iSel) {
    case kAllArbitrated:
      if (iMuon->StandaloneTrk() != 0 || iMuon->GlobalTrk()!= 0)  
        return kTRUE;
      if (iMuon->NSegments() > 0) 
        return kTRUE;
      break;
    case kPromptTight: 
      return iMuon->PromptTight(Muon::kAny);
      break;
    case kTMOneStationLoose:
      return iMuon->TMOneStation(999999,999999);
      break;
    case kTMOneStationTight:
      return iMuon->TMOneStation();
      break;
    case kTMLastStationLoose: 
      return iMuon->TMLastStation(999999,999999);
      break;
    case kTMLastStationTight: 
      return iMuon->TMLastStation();
      break;
    case kTM2DCompatibilityLoose: 
      tm2dcut = 0.7;
      break;
    case kTM2DCompatibilityTight: 
      tm2dcut = 1.0;
      break;
    default:
      return kFALSE;
      break;
  }

  Double_t lVal = GetSegmentCompatability(iMuon); 
  if (lVal == 0.5) // exclude this border case
    return kFALSE;

  lVal *= 1.2;
  lVal += 0.8*GetCaloCompatability(iMuon,kTRUE,kTRUE);
  if (lVal > tm2dcut) 
    return kTRUE;

  return kFALSE;
}

//--------------------------------------------------------------------------------------------------
Double_t MuonTools::GetSegmentCompatability(const mithep::Muon *iMuon) const
{
  // Get segment compatability for given muon based on likelihood of well defined 
  // track through chambers.

  Int_t lNStationsCrossed = 0;
  Int_t lNStationsSegment = 0;
  
  Int_t lStSegmentmatch[8];
  Int_t lStCrossed[8];
  Double_t lStBoundary[8];

  Double_t lWeight  = 0.;
  Bool_t lAdjust    = kTRUE;
  for (Int_t i0 = 0; i0 < 8; ++i0) {
    lStBoundary[i0] = 0.;
    if(iMuon->GetTrackDist(i0) < 999999. ) { 
      lNStationsCrossed++;
      lStCrossed[i0]  = 1;
      if (iMuon->GetTrackDist(i0) > -10. ) 
        lStBoundary[i0] = iMuon->GetTrackDist(i0); 
    } else
      lStCrossed[i0]  = 0;

    if(iMuon->GetDX(i0) < 999999.) {
      lNStationsSegment++;
      lStSegmentmatch[i0] = 1;
    } else
      lStSegmentmatch[i0] = 0;

    if(iMuon->GetDY(i0) < 999999.) 
      lAdjust = kFALSE;
  }

  if (lNStationsCrossed == 0)
    return 0.5;

  Double_t lStWeight[8];
  Int_t lPCross = -1;
  const Double_t lAtWeight = 0.5; 
  for (Int_t i0 = 0; i0< 8; ++i0) { 
    lStWeight[i0] = 0;
    if (lStCrossed[i0] > 0) { 
      lPCross++;

      switch (lNStationsCrossed) { 
        case 1 : 
          lStWeight[i0] =  1.;
          break;
        case 2 :
          if (lPCross == 0 ) 
            lStWeight[i0] = 0.33;
          else 
            lStWeight[i0] = 0.67;
          break;
        case 3 : 
          if (lPCross == 0) 
            lStWeight[i0] = 0.23;
          else if (lPCross == 1) 
            lStWeight[i0] = 0.33;
          else 
            lStWeight[i0] = 0.44;
          break;
        case 4 : 
          if (lPCross == 0) 
            lStWeight[i0] = 0.10;
          else if (lPCross == 1) 
            lStWeight[i0] = 0.20;
          else if (lPCross == 2) 
            lStWeight[i0] = 0.30;
          else                    
            lStWeight[i0] = 0.40;
          break;
        default : 
          lStWeight[i0] = 1./lNStationsCrossed;
      }
    
      if (lStSegmentmatch[i0] <= 0 && lStBoundary[i0] != 0.)
        lStWeight[i0] *= lAtWeight*0.5*(TMath::Erf(lStBoundary[i0]/6.)+1.); 
      else if (lStSegmentmatch[i0] <= 0 && lStBoundary[i0] == 0) 
        lStWeight[i0] = 0.;
      
      if (lStSegmentmatch[i0] > 0) { 
        Double_t lP2X = TMath::Power(iMuon->GetPullX(i0),2.);
        Double_t lP2Y = TMath::Power(iMuon->GetPullY(i0),2.);
        Double_t lD2X = TMath::Power(iMuon->GetDX(i0),2.);
        Double_t lD2Y = TMath::Power(iMuon->GetDY(i0),2.);
        if (iMuon->GetDY(i0) < 999999 && iMuon->GetDX(i0) < 999999) 
          lStWeight[i0] *= SigWeight(TMath::Sqrt(lD2X+lD2Y),TMath::Sqrt(lP2X+lP2Y));
        else if (iMuon->GetDY(i0) >= 999999 && i0 < 4)
          lStWeight[i0] *= SigWeight(iMuon->GetDX(i0),iMuon->GetPullX(i0));
        else if(i0 < 4)
          lStWeight[i0] *= SigWeight(iMuon->GetDY(i0),iMuon->GetPullY(i0));
      }
    } 
    lWeight += lStWeight[i0];
  }

  return lWeight;
}

//--------------------------------------------------------------------------------------------------
TH2D *MuonTools::LoadHisto(const char *name, TFile *file) const
{
  // Load histogram with given name from given file and return it.

  TH2D *ret = dynamic_cast<TH2D*>(file->Get(name));
  if (!ret) {
    Fatal("LoadHisto", "Could not load histogram %s from file %s", name, file->GetName());
    return 0;
  }
  ret->SetDirectory(0);
  return ret;
}

Revision:	1.9
Committed:	Tue Apr 7 15:37:10 2009 UTC (16 years ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_009c, Mit_009b, Mit_009a, Mit_009
Changes since 1.8:	+5 -3 lines
Log Message:	Cleanup
#	User	Rev	Content
1	loizides	1.9	// $Id: MuonTools.cc,v 1.8 2009/03/23 14:48:08 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	loizides	1.9	// Get calo compatibility value for given muon based on calorimeter templates.
95			// If iEMSpecial is true, then a use different arrangement of ECAL for compatibility.
96	loizides	1.5
97			Double_t lEta = iMuon->Eta();
98			Double_t aEta = TMath::Abs(lEta);
99			if (aEta > 2.5)
100			return 0.5;
101
102			Double_t lP = iMuon->P();
103			if (lP >= 2000.)
104			lP = 1999.9;
105			if(lP < 0. )
106			return 0.5;
107
108			Double_t lEM = -5.;
109			if (!iEMSpecial \|\| iMuon->EmEnergy() != 0.)
110			lEM = iMuon->EmEnergy();
111
112			Double_t lHad = iMuon->HadEnergy();
113			Double_t lHO = iMuon->HoEnergy();;
114
115			TH2D *lTMuonHad = 0;
116			TH2D *lTPionHad = 0;
117			TH2D *lTMuonHo = 0;
118			TH2D *lTPionHo = 0;
119			TH2D *lTMuonEm = 0;
120			TH2D *lTPionEm = 0;
121
122			if (aEta >= 1.27) {
123			if (iCorrectedHCAL)
124			lHad *= 1.8/2.2;
125			if (lEta > 0) {
126	pharris	1.1	lTPionHad = fpion_had_etaEpl;
127			lTMuonHad = fmuon_had_etaEpl;
128			} else {
129			lTPionHad = fpion_had_etaEmi;
130			lTMuonHad = fmuon_had_etaEmi;
131			}
132			}
133	loizides	1.5
134			if (aEta < 1.27 && aEta >= 1.1) {
135			if (iCorrectedHCAL)
136			lHad = (1.8/(-2.2aEta+5.5));
137			if (lEta > 0) {
138	pharris	1.1	lTPionHad = fpion_had_etaTpl;
139			lTMuonHad = fmuon_had_etaTpl;
140			} else {
141			lTPionHad = fpion_had_etaTmi;
142			lTMuonHad = fmuon_had_etaTmi;
143			}
144			}
145	loizides	1.5
146			if (aEta < 1.1) {
147			if(iCorrectedHCAL)
148	loizides	1.7	lHad = TMath::Sin(2TMath::ATan(TMath::Exp(lEta)));
149	pharris	1.1	lTPionHad = fpion_had_etaB;
150			lTMuonHad = fmuon_had_etaB;
151			}
152	loizides	1.5	if (lEta > 1.479) {
153			lTPionEm = fpion_em_etaEpl;
154			lTMuonEm = fmuon_em_etaEpl;
155	pharris	1.1	}
156	loizides	1.5	if (aEta <= 1.479) {
157	pharris	1.1	lTPionEm = fpion_em_etaB;
158			lTMuonEm = fmuon_em_etaB;
159			}
160	loizides	1.5	if (lEta < -1.479) {
161	pharris	1.1	lTPionEm = fpion_em_etaEmi;
162			lTMuonEm = fmuon_em_etaEmi;
163			}
164	loizides	1.5	if (aEta < 1.28) {
165	pharris	1.1	lTPionHo = fpion_ho_etaB;
166			lTMuonHo = fmuon_ho_etaB;
167			}
168
169	loizides	1.5	Double_t lPBX = 1.;
170			Double_t lPSX = 1.;
171			Double_t lPBY = 1.;
172			Double_t lPSY = 1.;
173			Double_t lPBZ = 1.;
174			Double_t lPSZ = 1.;
175			if (!Overflow(lTPionEm, lP,lEM))
176			lPBX = lTPionEm ->GetBinContent(lTPionEm ->GetXaxis()->FindBin(lP),
177			lTPionEm ->GetYaxis()->FindBin(lEM));
178			if (!Overflow(lTPionHad,lP,lHad))
179			lPBY = lTPionHad->GetBinContent(lTPionHad->GetXaxis()->FindBin(lP),
180			lTPionHad->GetYaxis()->FindBin(lHad));
181			if (!Overflow(lTPionHo, lP,lHO))
182			lPBZ = lTPionHo ->GetBinContent(lTPionHo ->GetXaxis()->FindBin(lP),
183			lTPionHo ->GetYaxis()->FindBin(lHO));
184			if (!Overflow(lTMuonEm, lP,lEM ))
185			lPSX = lTMuonEm ->GetBinContent(lTMuonEm ->GetXaxis()->FindBin(lP),
186			lTMuonEm ->GetYaxis()->FindBin(lEM));
187			if (!Overflow(lTMuonHad,lP,lHad))
188			lPSY = lTMuonHad->GetBinContent(lTMuonHad->GetXaxis()->FindBin(lP),
189			lTMuonHad->GetYaxis()->FindBin(lHad));
190			if (!Overflow(lTMuonHo ,lP,lHO))
191			lPSZ = lTMuonHo ->GetBinContent(lTMuonHo ->GetXaxis()->FindBin(lP),
192			lTMuonHo ->GetYaxis()->FindBin(lHO));
193	pharris	1.1
194	loizides	1.5	if (lPSX == 0. \|\| lPBX == 0. \|\| (lEM <= 0. && !iEMSpecial)) {
195			lPSX = 1.;
196			lPBX = 1.;
197			}
198			if (lPSY == 0. \|\| lPBY == 0. \|\| lHad == 0.) {
199			lPSY = 1.;
200			lPBY = 1.;
201			}
202			if (lPSZ == 0. \|\| lPBZ == 0. \|\| lHO == 0.) {
203			lPSZ = 1.;
204			lPBZ = 1.;
205			}
206			if ((lPSXlPSYlPSZ+lPBXlPBYlPBZ) > 0.)
207			return lPSXlPSYlPSZ/(lPSXlPSYlPSZ+lPBXlPBYlPBZ);
208
209	pharris	1.1	return 0.5;
210			}
211
212	loizides	1.5	//--------------------------------------------------------------------------------------------------
213			Bool_t MuonTools::Init(const char mutemp, const char pitemp)
214			{
215			// Read histograms from given files.
216
217			if (fIsInit) {
218			DeleteHistos();
219			}
220
221			TDirectory::TContext context(0);
222
223			TFile *muon_templates = TFile::Open(mutemp);
224			if (!muon_templates) {
225			Fatal("Init", "Could not open file %s", mutemp);
226			return kFALSE;
227			}
228			fmuon_em_etaEmi = LoadHisto("em_etaEmi", muon_templates);
229			fmuon_had_etaEmi = LoadHisto("had_etaEmi", muon_templates);
230			fmuon_had_etaTmi = LoadHisto("had_etaTmi", muon_templates);
231			fmuon_em_etaB = LoadHisto("em_etaB", muon_templates);
232			fmuon_had_etaB = LoadHisto("had_etaB", muon_templates);
233			fmuon_ho_etaB = LoadHisto("ho_etaB", muon_templates);
234			fmuon_had_etaTpl = LoadHisto("had_etaTpl", muon_templates);
235			fmuon_em_etaEpl = LoadHisto("em_etaEpl", muon_templates);
236			fmuon_had_etaEpl = LoadHisto("had_etaEpl", muon_templates);
237			muon_templates->Close();
238			delete muon_templates;
239
240			TFile *pion_templates = TFile::Open(pitemp);
241			if (!pion_templates) {
242			Fatal("Init", "Could not open file %s", pitemp);
243			return kFALSE;
244			}
245
246			fpion_em_etaEmi = LoadHisto("em_etaEmi", pion_templates);
247			fpion_had_etaEmi = LoadHisto("had_etaEmi", pion_templates);
248			fpion_had_etaTmi = LoadHisto("had_etaTmi", pion_templates);
249			fpion_em_etaB = LoadHisto("em_etaB", pion_templates);
250			fpion_had_etaB = LoadHisto("had_etaB", pion_templates);
251			fpion_ho_etaB = LoadHisto("ho_etaB", pion_templates);
252			fpion_had_etaTpl = LoadHisto("had_etaTpl", pion_templates);
253			fpion_em_etaEpl = LoadHisto("em_etaEpl", pion_templates);
254			fpion_had_etaEpl = LoadHisto("had_etaEpl", pion_templates);
255			pion_templates->Close();
256			delete pion_templates;
257
258			fIsInit = kTRUE;
259			return kTRUE;
260			}
261
262			//--------------------------------------------------------------------------------------------------
263			Bool_t MuonTools::IsGood(const mithep::Muon *iMuon, ESelType iSel) const
264			{
265			// Return true if given muon qualifies given selection criterium.
266
267			Double_t tm2dcut = 0.;
268
269			switch(iSel) {
270			case kAllArbitrated:
271			if (iMuon->StandaloneTrk() != 0 \|\| iMuon->GlobalTrk()!= 0)
272			return kTRUE;
273			if (iMuon->NSegments() > 0)
274			return kTRUE;
275			break;
276			case kPromptTight:
277			return iMuon->PromptTight(Muon::kAny);
278			break;
279			case kTMOneStationLoose:
280	loizides	1.8	return iMuon->TMOneStation(999999,999999);
281	loizides	1.5	break;
282			case kTMOneStationTight:
283			return iMuon->TMOneStation();
284			break;
285			case kTMLastStationLoose:
286			return iMuon->TMLastStation(999999,999999);
287			break;
288			case kTMLastStationTight:
289			return iMuon->TMLastStation();
290			break;
291			case kTM2DCompatibilityLoose:
292			tm2dcut = 0.7;
293			break;
294			case kTM2DCompatibilityTight:
295			tm2dcut = 1.0;
296			break;
297			default:
298			return kFALSE;
299			break;
300			}
301
302	loizides	1.8	Double_t lVal = GetSegmentCompatability(iMuon);
303			if (lVal == 0.5) // exclude this border case
304	loizides	1.5	return kFALSE;
305
306	loizides	1.8	lVal *= 1.2;
307	loizides	1.5	lVal += 0.8*GetCaloCompatability(iMuon,kTRUE,kTRUE);
308			if (lVal > tm2dcut)
309			return kTRUE;
310	loizides	1.8
311	loizides	1.5	return kFALSE;
312			}
313
314			//--------------------------------------------------------------------------------------------------
315			Double_t MuonTools::GetSegmentCompatability(const mithep::Muon *iMuon) const
316			{
317	loizides	1.9	// Get segment compatability for given muon based on likelihood of well defined
318			// track through chambers.
319	loizides	1.5
320			Int_t lNStationsCrossed = 0;
321			Int_t lNStationsSegment = 0;
322
323			Int_t lStSegmentmatch[8];
324			Int_t lStCrossed[8];
325			Double_t lStBoundary[8];
326
327			Double_t lWeight = 0.;
328			Bool_t lAdjust = kTRUE;
329			for (Int_t i0 = 0; i0 < 8; ++i0) {
330			lStBoundary[i0] = 0.;
331			if(iMuon->GetTrackDist(i0) < 999999. ) {
332			lNStationsCrossed++;
333			lStCrossed[i0] = 1;
334			if (iMuon->GetTrackDist(i0) > -10. )
335			lStBoundary[i0] = iMuon->GetTrackDist(i0);
336			} else
337			lStCrossed[i0] = 0;
338
339	loizides	1.7	if(iMuon->GetDX(i0) < 999999.) {
340	loizides	1.5	lNStationsSegment++;
341			lStSegmentmatch[i0] = 1;
342			} else
343			lStSegmentmatch[i0] = 0;
344
345			if(iMuon->GetDY(i0) < 999999.)
346			lAdjust = kFALSE;
347			}
348
349			if (lNStationsCrossed == 0)
350			return 0.5;
351
352			Double_t lStWeight[8];
353			Int_t lPCross = -1;
354			const Double_t lAtWeight = 0.5;
355			for (Int_t i0 = 0; i0< 8; ++i0) {
356			lStWeight[i0] = 0;
357			if (lStCrossed[i0] > 0) {
358			lPCross++;
359
360			switch (lNStationsCrossed) {
361			case 1 :
362			lStWeight[i0] = 1.;
363			break;
364			case 2 :
365			if (lPCross == 0 )
366			lStWeight[i0] = 0.33;
367			else
368			lStWeight[i0] = 0.67;
369			break;
370			case 3 :
371			if (lPCross == 0)
372			lStWeight[i0] = 0.23;
373			else if (lPCross == 1)
374			lStWeight[i0] = 0.33;
375			else
376			lStWeight[i0] = 0.44;
377			break;
378			case 4 :
379			if (lPCross == 0)
380			lStWeight[i0] = 0.10;
381			else if (lPCross == 1)
382			lStWeight[i0] = 0.20;
383			else if (lPCross == 2)
384			lStWeight[i0] = 0.30;
385			else
386			lStWeight[i0] = 0.40;
387			break;
388			default :
389			lStWeight[i0] = 1./lNStationsCrossed;
390			}
391
392			if (lStSegmentmatch[i0] <= 0 && lStBoundary[i0] != 0.)
393			lStWeight[i0] = lAtWeight0.5*(TMath::Erf(lStBoundary[i0]/6.)+1.);
394			else if (lStSegmentmatch[i0] <= 0 && lStBoundary[i0] == 0)
395			lStWeight[i0] = 0.;
396
397			if (lStSegmentmatch[i0] > 0) {
398			Double_t lP2X = TMath::Power(iMuon->GetPullX(i0),2.);
399			Double_t lP2Y = TMath::Power(iMuon->GetPullY(i0),2.);
400			Double_t lD2X = TMath::Power(iMuon->GetDX(i0),2.);
401			Double_t lD2Y = TMath::Power(iMuon->GetDY(i0),2.);
402			if (iMuon->GetDY(i0) < 999999 && iMuon->GetDX(i0) < 999999)
403			lStWeight[i0] *= SigWeight(TMath::Sqrt(lD2X+lD2Y),TMath::Sqrt(lP2X+lP2Y));
404			else if (iMuon->GetDY(i0) >= 999999 && i0 < 4)
405			lStWeight[i0] *= SigWeight(iMuon->GetDX(i0),iMuon->GetPullX(i0));
406			else if(i0 < 4)
407			lStWeight[i0] *= SigWeight(iMuon->GetDY(i0),iMuon->GetPullY(i0));
408			}
409			}
410			lWeight += lStWeight[i0];
411			}
412
413			return lWeight;
414			}
415
416			//--------------------------------------------------------------------------------------------------
417			TH2D MuonTools::LoadHisto(const char name, TFile *file) const
418			{
419			// Load histogram with given name from given file and return it.
420
421			TH2D ret = dynamic_cast<TH2D>(file->Get(name));
422			if (!ret) {
423			Fatal("LoadHisto", "Could not load histogram %s from file %s", name, file->GetName());
424			return 0;
425	pharris	1.1	}
426	loizides	1.5	ret->SetDirectory(0);
427			return ret;
428	pharris	1.1	}
429	loizides	1.5