Utils/src/MuonTools.cc

// $Id: MuonTools.cc,v 1.21 2011/10/02 10:05:42 ceballos Exp $

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

ClassImp(mithep::MuonTools)

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;
}
//--------------------------------------------------------------------------------------------------
Bool_t MuonTools::PassD0Cut(const Muon *mu, const VertexCol *vertices, Double_t fD0Cut, Int_t nVertex) 
{
  Bool_t d0cut = kFALSE;
  const Track *mt = mu->BestTrk();
  if (!mt) return kFALSE;

  Double_t d0_real = 1e30;
  if(nVertex >= 0) d0_real = TMath::Abs(mt->D0Corrected(*vertices->At(nVertex)));
  else            {
    Double_t distVtx = 999.0;
    Int_t closestVtx = 0;
    for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
      double dz = TMath::Abs(mt->DzCorrected(*vertices->At(nv)));
      if(dz < distVtx) {
        distVtx    = dz;
        closestVtx = nv;
      }
    }
    d0_real = TMath::Abs(mt->D0Corrected(*vertices->At(closestVtx)));
  }
  if(d0_real < fD0Cut) d0cut = kTRUE;
  
  return d0cut;
}

//--------------------------------------------------------------------------------------------------
Bool_t MuonTools::PassD0Cut(const Muon *mu, const BeamSpotCol *beamspots, Double_t fD0Cut) 
{
  Bool_t d0cut = kFALSE;
  const Track *mt = mu->BestTrk();
  if (!mt) return kFALSE;

  // d0 cut
  Double_t d0_real = 99999;
  for(UInt_t i0 = 0; i0 < beamspots->GetEntries(); i0++) {
    Double_t pD0 = mt->D0Corrected(*beamspots->At(i0));
    if(TMath::Abs(pD0) < TMath::Abs(d0_real)) d0_real = TMath::Abs(pD0);
  }
  if(d0_real < fD0Cut) d0cut = kTRUE;
  
  return d0cut;
}

//--------------------------------------------------------------------------------------------------
Bool_t MuonTools::PassDZCut(const Muon *mu, const VertexCol *vertices, Double_t fDZCut, Int_t nVertex) 
{
  Bool_t dzcut = kFALSE;
  const Track *mt = mu->BestTrk();
  if (!mt) return kFALSE;

  Double_t distVtx = 999.0;
  if(nVertex >= 0) distVtx = TMath::Abs(mt->DzCorrected(*vertices->At(nVertex)));
  else {
    for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
      double dz = TMath::Abs(mt->DzCorrected(*vertices->At(nv)));
      if(dz < distVtx) {
        distVtx = dz;
      }
    }
  }

  if(distVtx < fDZCut) dzcut = kTRUE;
  
  return dzcut;
}

//--------------------------------------------------------------------------------------------------
Bool_t MuonTools::PassSoftMuonCut(const Muon *mu, const VertexCol *vertices, const Double_t fDZCut,
                                  const Bool_t applyIso) 
{
  if(mu->Pt() <= 3.0) return kFALSE;

  if(!mu->IsTrackerMuon()) return kFALSE;
  
  if(!mu->Quality().Quality(MuonQuality::TMLastStationAngTight)) return kFALSE;

  if(mu->BestTrk()->NHits() <= 10) return kFALSE;

  if(!PassD0Cut(mu, vertices, 0.2, 0)) return kFALSE;

  if(!PassDZCut(mu, vertices, fDZCut, 0)) return kFALSE;

  if(applyIso == kTRUE){
    Double_t totalIso = 1.0 * mu->IsoR03SumPt() + 
                        1.0 * mu->IsoR03EmEt() + 
                        1.0 * mu->IsoR03HadEt();
    if (totalIso < (mu->Pt()*0.10) && mu->Pt() > 20.0) return kFALSE;
  }

  return kTRUE;
}
Revision:	1.22
Committed:	Wed Oct 5 16:11:10 2011 UTC (13 years, 7 months ago) by ceballos
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_025c, Mit_025b, Mit_025a, Mit_025
Changes since 1.21:	+9 -6 lines
Log Message:	new
#	User	Rev	Content
1	ceballos	1.22	// $Id: MuonTools.cc,v 1.21 2011/10/02 10:05:42 ceballos Exp $
2	loizides	1.5
3	pharris	1.2	#include "MitPhysics/Utils/interface/MuonTools.h"
4	loizides	1.5	#include <TFile.h>
5
6	loizides	1.10	ClassImp(mithep::MuonTools)
7
8	pharris	1.1	using namespace mithep;
9
10	loizides	1.5	//--------------------------------------------------------------------------------------------------
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),
16			fmuon_em_etaB(0),
17			fmuon_had_etaB(0),
18			fmuon_ho_etaB(0),
19			fmuon_had_etaTpl(0),
20			fmuon_em_etaEpl(0),
21			fmuon_had_etaEpl(0),
22	pharris	1.1	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	loizides	1.5	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	loizides	1.9	// 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	loizides	1.5
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	pharris	1.1	lTPionHad = fpion_had_etaEpl;
129			lTMuonHad = fmuon_had_etaEpl;
130			} else {
131			lTPionHad = fpion_had_etaEmi;
132			lTMuonHad = fmuon_had_etaEmi;
133			}
134			}
135	loizides	1.5
136			if (aEta < 1.27 && aEta >= 1.1) {
137			if (iCorrectedHCAL)
138			lHad = (1.8/(-2.2aEta+5.5));
139			if (lEta > 0) {
140	pharris	1.1	lTPionHad = fpion_had_etaTpl;
141			lTMuonHad = fmuon_had_etaTpl;
142			} else {
143			lTPionHad = fpion_had_etaTmi;
144			lTMuonHad = fmuon_had_etaTmi;
145			}
146			}
147	loizides	1.5
148			if (aEta < 1.1) {
149			if(iCorrectedHCAL)
150	loizides	1.7	lHad = TMath::Sin(2TMath::ATan(TMath::Exp(lEta)));
151	pharris	1.1	lTPionHad = fpion_had_etaB;
152			lTMuonHad = fmuon_had_etaB;
153			}
154	loizides	1.5	if (lEta > 1.479) {
155			lTPionEm = fpion_em_etaEpl;
156			lTMuonEm = fmuon_em_etaEpl;
157	pharris	1.1	}
158	loizides	1.5	if (aEta <= 1.479) {
159	pharris	1.1	lTPionEm = fpion_em_etaB;
160			lTMuonEm = fmuon_em_etaB;
161			}
162	loizides	1.5	if (lEta < -1.479) {
163	pharris	1.1	lTPionEm = fpion_em_etaEmi;
164			lTMuonEm = fmuon_em_etaEmi;
165			}
166	loizides	1.5	if (aEta < 1.28) {
167	pharris	1.1	lTPionHo = fpion_ho_etaB;
168			lTMuonHo = fmuon_ho_etaB;
169			}
170
171	loizides	1.5	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	pharris	1.1
196	loizides	1.5	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 ((lPSXlPSYlPSZ+lPBXlPBYlPBZ) > 0.)
209			return lPSXlPSYlPSZ/(lPSXlPSYlPSZ+lPBXlPBYlPBZ);
210
211	pharris	1.1	return 0.5;
212			}
213
214	loizides	1.5	//--------------------------------------------------------------------------------------------------
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
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	loizides	1.8	return iMuon->TMOneStation(999999,999999);
283	loizides	1.5	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	loizides	1.8	Double_t lVal = GetSegmentCompatability(iMuon);
305			if (lVal == 0.5) // exclude this border case
306	loizides	1.5	return kFALSE;
307
308	loizides	1.8	lVal *= 1.2;
309	loizides	1.5	lVal += 0.8*GetCaloCompatability(iMuon,kTRUE,kTRUE);
310			if (lVal > tm2dcut)
311			return kTRUE;
312	loizides	1.8
313	loizides	1.5	return kFALSE;
314			}
315
316			//--------------------------------------------------------------------------------------------------
317			Double_t MuonTools::GetSegmentCompatability(const mithep::Muon *iMuon) const
318			{
319	loizides	1.9	// Get segment compatability for given muon based on likelihood of well defined
320			// track through chambers.
321	loizides	1.5
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	loizides	1.7	if(iMuon->GetDX(i0) < 999999.) {
342	loizides	1.5	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] = lAtWeight0.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	pharris	1.1	}
428	loizides	1.5	ret->SetDirectory(0);
429			return ret;
430	pharris	1.1	}
431	ceballos	1.11	//--------------------------------------------------------------------------------------------------
432	ceballos	1.15	Bool_t MuonTools::PassD0Cut(const Muon mu, const VertexCol vertices, Double_t fD0Cut, Int_t nVertex)
433	ceballos	1.11	{
434			Bool_t d0cut = kFALSE;
435			const Track *mt = mu->BestTrk();
436			if (!mt) return kFALSE;
437
438			Double_t d0_real = 1e30;
439	ceballos	1.15	if(nVertex >= 0) d0_real = TMath::Abs(mt->D0Corrected(*vertices->At(nVertex)));
440			else {
441			Double_t distVtx = 999.0;
442			Int_t closestVtx = 0;
443			for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
444			double dz = TMath::Abs(mt->DzCorrected(*vertices->At(nv)));
445			if(dz < distVtx) {
446			distVtx = dz;
447			closestVtx = nv;
448			}
449	ceballos	1.11	}
450	ceballos	1.15	d0_real = TMath::Abs(mt->D0Corrected(*vertices->At(closestVtx)));
451	ceballos	1.11	}
452			if(d0_real < fD0Cut) d0cut = kTRUE;
453
454			return d0cut;
455			}
456
457			//--------------------------------------------------------------------------------------------------
458	ceballos	1.14	Bool_t MuonTools::PassD0Cut(const Muon mu, const BeamSpotCol beamspots, Double_t fD0Cut)
459	ceballos	1.11	{
460			Bool_t d0cut = kFALSE;
461			const Track *mt = mu->BestTrk();
462			if (!mt) return kFALSE;
463
464			// d0 cut
465			Double_t d0_real = 99999;
466			for(UInt_t i0 = 0; i0 < beamspots->GetEntries(); i0++) {
467			Double_t pD0 = mt->D0Corrected(*beamspots->At(i0));
468			if(TMath::Abs(pD0) < TMath::Abs(d0_real)) d0_real = TMath::Abs(pD0);
469			}
470			if(d0_real < fD0Cut) d0cut = kTRUE;
471
472			return d0cut;
473			}
474	loizides	1.5
475	ceballos	1.12	//--------------------------------------------------------------------------------------------------
476	ceballos	1.16	Bool_t MuonTools::PassDZCut(const Muon mu, const VertexCol vertices, Double_t fDZCut, Int_t nVertex)
477	ceballos	1.15	{
478			Bool_t dzcut = kFALSE;
479			const Track *mt = mu->BestTrk();
480			if (!mt) return kFALSE;
481
482			Double_t distVtx = 999.0;
483	ceballos	1.16	if(nVertex >= 0) distVtx = TMath::Abs(mt->DzCorrected(*vertices->At(nVertex)));
484			else {
485			for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
486			double dz = TMath::Abs(mt->DzCorrected(*vertices->At(nv)));
487			if(dz < distVtx) {
488			distVtx = dz;
489			}
490	ceballos	1.15	}
491			}
492
493			if(distVtx < fDZCut) dzcut = kTRUE;
494
495			return dzcut;
496			}
497
498			//--------------------------------------------------------------------------------------------------
499	ceballos	1.22	Bool_t MuonTools::PassSoftMuonCut(const Muon mu, const VertexCol vertices, const Double_t fDZCut,
500			const Bool_t applyIso)
501	ceballos	1.12	{
502			if(mu->Pt() <= 3.0) return kFALSE;
503	sixie	1.18
504	ceballos	1.12	if(!mu->IsTrackerMuon()) return kFALSE;
505
506			if(!mu->Quality().Quality(MuonQuality::TMLastStationAngTight)) return kFALSE;
507	ceballos	1.21
508	ceballos	1.12	if(mu->BestTrk()->NHits() <= 10) return kFALSE;
509
510	ceballos	1.17	if(!PassD0Cut(mu, vertices, 0.2, 0)) return kFALSE;
511
512	ceballos	1.20	if(!PassDZCut(mu, vertices, fDZCut, 0)) return kFALSE;
513	ceballos	1.21
514	ceballos	1.22	if(applyIso == kTRUE){
515			Double_t totalIso = 1.0 * mu->IsoR03SumPt() +
516			1.0 * mu->IsoR03EmEt() +
517			1.0 * mu->IsoR03HadEt();
518			if (totalIso < (mu->Pt()*0.10) && mu->Pt() > 20.0) return kFALSE;
519			}
520	ceballos	1.13
521	ceballos	1.12	return kTRUE;
522			}