Utils/src/MuonTools.cc

// $Id: MuonTools.cc,v 1.7 2008/12/04 13:53:34 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.

  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.

  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.8
Committed:	Mon Mar 23 14:48:08 2009 UTC (16 years, 1 month ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_008
Changes since 1.7:	+8 -6 lines
Log Message:	Cosmetics
#	User	Rev	Content
1	loizides	1.8	// $Id: MuonTools.cc,v 1.7 2008/12/04 13:53:34 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.8	// Get calo compatibility value for given muon.
95	loizides	1.5
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.2aEta+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	loizides	1.7	lHad = TMath::Sin(2TMath::ATan(TMath::Exp(lEta)));
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 ((lPSXlPSYlPSZ+lPBXlPBYlPBZ) > 0.)
206			return lPSXlPSYlPSZ/(lPSXlPSYlPSZ+lPBXlPBYlPBZ);
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	loizides	1.8	return iMuon->TMOneStation(999999,999999);
280	loizides	1.5	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	loizides	1.8	Double_t lVal = GetSegmentCompatability(iMuon);
302			if (lVal == 0.5) // exclude this border case
303	loizides	1.5	return kFALSE;
304
305	loizides	1.8	lVal *= 1.2;
306	loizides	1.5	lVal += 0.8*GetCaloCompatability(iMuon,kTRUE,kTRUE);
307			if (lVal > tm2dcut)
308			return kTRUE;
309	loizides	1.8
310	loizides	1.5	return kFALSE;
311			}
312
313			//--------------------------------------------------------------------------------------------------
314			Double_t MuonTools::GetSegmentCompatability(const mithep::Muon *iMuon) const
315			{
316	loizides	1.8	// Get segment compatability for given muon.
317	loizides	1.5
318			Int_t lNStationsCrossed = 0;
319			Int_t lNStationsSegment = 0;
320
321			Int_t lStSegmentmatch[8];
322			Int_t lStCrossed[8];
323			Double_t lStBoundary[8];
324
325			Double_t lWeight = 0.;
326			Bool_t lAdjust = kTRUE;
327			for (Int_t i0 = 0; i0 < 8; ++i0) {
328			lStBoundary[i0] = 0.;
329			if(iMuon->GetTrackDist(i0) < 999999. ) {
330			lNStationsCrossed++;
331			lStCrossed[i0] = 1;
332			if (iMuon->GetTrackDist(i0) > -10. )
333			lStBoundary[i0] = iMuon->GetTrackDist(i0);
334			} else
335			lStCrossed[i0] = 0;
336
337	loizides	1.7	if(iMuon->GetDX(i0) < 999999.) {
338	loizides	1.5	lNStationsSegment++;
339			lStSegmentmatch[i0] = 1;
340			} else
341			lStSegmentmatch[i0] = 0;
342
343			if(iMuon->GetDY(i0) < 999999.)
344			lAdjust = kFALSE;
345			}
346
347			if (lNStationsCrossed == 0)
348			return 0.5;
349
350			Double_t lStWeight[8];
351			Int_t lPCross = -1;
352			const Double_t lAtWeight = 0.5;
353			for (Int_t i0 = 0; i0< 8; ++i0) {
354			lStWeight[i0] = 0;
355			if (lStCrossed[i0] > 0) {
356			lPCross++;
357
358			switch (lNStationsCrossed) {
359			case 1 :
360			lStWeight[i0] = 1.;
361			break;
362			case 2 :
363			if (lPCross == 0 )
364			lStWeight[i0] = 0.33;
365			else
366			lStWeight[i0] = 0.67;
367			break;
368			case 3 :
369			if (lPCross == 0)
370			lStWeight[i0] = 0.23;
371			else if (lPCross == 1)
372			lStWeight[i0] = 0.33;
373			else
374			lStWeight[i0] = 0.44;
375			break;
376			case 4 :
377			if (lPCross == 0)
378			lStWeight[i0] = 0.10;
379			else if (lPCross == 1)
380			lStWeight[i0] = 0.20;
381			else if (lPCross == 2)
382			lStWeight[i0] = 0.30;
383			else
384			lStWeight[i0] = 0.40;
385			break;
386			default :
387			lStWeight[i0] = 1./lNStationsCrossed;
388			}
389
390			if (lStSegmentmatch[i0] <= 0 && lStBoundary[i0] != 0.)
391			lStWeight[i0] = lAtWeight0.5*(TMath::Erf(lStBoundary[i0]/6.)+1.);
392			else if (lStSegmentmatch[i0] <= 0 && lStBoundary[i0] == 0)
393			lStWeight[i0] = 0.;
394
395			if (lStSegmentmatch[i0] > 0) {
396			Double_t lP2X = TMath::Power(iMuon->GetPullX(i0),2.);
397			Double_t lP2Y = TMath::Power(iMuon->GetPullY(i0),2.);
398			Double_t lD2X = TMath::Power(iMuon->GetDX(i0),2.);
399			Double_t lD2Y = TMath::Power(iMuon->GetDY(i0),2.);
400			if (iMuon->GetDY(i0) < 999999 && iMuon->GetDX(i0) < 999999)
401			lStWeight[i0] *= SigWeight(TMath::Sqrt(lD2X+lD2Y),TMath::Sqrt(lP2X+lP2Y));
402			else if (iMuon->GetDY(i0) >= 999999 && i0 < 4)
403			lStWeight[i0] *= SigWeight(iMuon->GetDX(i0),iMuon->GetPullX(i0));
404			else if(i0 < 4)
405			lStWeight[i0] *= SigWeight(iMuon->GetDY(i0),iMuon->GetPullY(i0));
406			}
407			}
408			lWeight += lStWeight[i0];
409			}
410
411			return lWeight;
412			}
413
414			//--------------------------------------------------------------------------------------------------
415			TH2D MuonTools::LoadHisto(const char name, TFile *file) const
416			{
417			// Load histogram with given name from given file and return it.
418
419			TH2D ret = dynamic_cast<TH2D>(file->Get(name));
420			if (!ret) {
421			Fatal("LoadHisto", "Could not load histogram %s from file %s", name, file->GetName());
422			return 0;
423	pharris	1.1	}
424	loizides	1.5	ret->SetDirectory(0);
425			return ret;
426	pharris	1.1	}
427	loizides	1.5