[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / FourierDecomposition / AliMuonEffMC.cxx

// MUON track QA referring AliMuonEffMC.cxx
// Author : Saehanseul Oh

#include "AliMuonEffMC.h"

#include <TList.h>
#include <TH1D.h>
#include <TH2F.h>
#include <TH3F.h>
#include <THn.h>
#include <TChain.h>
#include <TFile.h>

#include "AliAnalysisManager.h"
#include "AliAnalysisTask.h"
#include "AliESDVertex.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliESDMuonTrack.h"
#include "AliAODTrack.h"
#include "AliESDVertex.h"
#include "AliAODVertex.h"
#include "AliCentrality.h"
#include "AliVParticle.h"
#include "AliMCParticle.h"
#include "AliAnalysisHelperJetTasks.h"
#include "AliMCEvent.h"
#include "AliAODMCParticle.h"

using std::cout;
using std::endl;

ClassImp(AliMuonEffMC)

//________________________________________________________________________
AliMuonEffMC::AliMuonEffMC() :
  AliAnalysisTaskSE(), fESD(0), fAOD(0), fMC(0), fCentrality(99), fZVertex(99), fOutputList(0x0),      
  fHEventStat(0), fHEvt(0x0), fIsMc(kTRUE), fMDProcess(kTRUE), fCentralityEstimator("V0M"),
  fNEtaBins(15), fNpTBins(100), fNCentBins(5), fNZvtxBins(10), fNPhiBins(24),
  fHMuonParGen(0x0), fHMuonDetGen(0x0), fHMuonDetRec(0x0), fHEtcDetRec(0x0)
{
  // Constructor
 
}

//________________________________________________________________________
AliMuonEffMC::AliMuonEffMC(const char *name) :
  AliAnalysisTaskSE(name), fESD(0), fAOD(0), fMC(0), fCentrality(99), fZVertex(99), fOutputList(0x0),      
  fHEventStat(0), fHEvt(0x0),  fIsMc(kTRUE), fMDProcess(kTRUE), fCentralityEstimator("V0M"),
  fNEtaBins(15), fNpTBins(100), fNCentBins(5), fNZvtxBins(10), fNPhiBins(24),
  fHMuonParGen(0x0), fHMuonDetGen(0x0), fHMuonDetRec(0x0), fHEtcDetRec(0x0)
{
  // Constructor
  DefineInput(0, TChain::Class());
  DefineOutput(1, TList::Class());

}

//________________________________________________________________________
AliMuonEffMC::~AliMuonEffMC()
{
  //Destructor
  if(fOutputList) delete fOutputList;
}

//________________________________________________________________________
void AliMuonEffMC::UserCreateOutputObjects()
{
  // Create histograms
  // Called once (per slave on PROOF!)

  fOutputList = new TList();
  fOutputList->SetOwner(1);
 
  for(Int_t i=0; i<4; i++)
  {
    fHMuMotherGenPt[i] = 0x0;
    fHMuMotherRecPt[i] = 0x0;
    fHMuMotherGenPhi[i] = 0x0;
    fHMuMotherRecPhi[i] = 0x0;
    fHMuMotherGenEta[i] = 0x0;
    fHMuMotherRecEta[i] = 0x0;
    fHMuDCA[i] = 0x0;
  }

  fHEventStat = new TH1D("fHEventStat","Event statistics for analysis",18,0,18);
  fHEventStat->GetXaxis()->SetBinLabel(1,"Event");
  fHEventStat->GetXaxis()->SetBinLabel(2,"SelectedEvent");
  fHEventStat->GetXaxis()->SetBinLabel(3,"File");
  fHEventStat->GetXaxis()->SetBinLabel(4,"fSPHighpt");  //!Global Trigger Single plus High p_T
  fHEventStat->GetXaxis()->SetBinLabel(5,"fSPAllpt");   //!Global Trigger Single plus All p_T
  fHEventStat->GetXaxis()->SetBinLabel(6,"fSMLowpt");   //!Global Trigger Single minus Low p_T
  fHEventStat->GetXaxis()->SetBinLabel(7,"fSMHighpt");  //!Global Trigger Single minus High p_T
  fHEventStat->GetXaxis()->SetBinLabel(8,"fSMAllpt");   //!Global Trigger Single minus All p_T
  fHEventStat->GetXaxis()->SetBinLabel(9,"fSULowpt");   //!Global Trigger Single undefined Low p_T
  fHEventStat->GetXaxis()->SetBinLabel(10,"fSUHighpt"); //!Global Trigger Single undefined High p_T
  fHEventStat->GetXaxis()->SetBinLabel(11,"fSUAllpt");  //!Global Trigger Single undefined All p_T
  fHEventStat->GetXaxis()->SetBinLabel(12,"fUSLowpt");  //!Global Trigger UnlikeSign pair Low p_T
  fHEventStat->GetXaxis()->SetBinLabel(13,"fUSHighpt"); //!Global Trigger UnlikeSign pair High p_T
  fHEventStat->GetXaxis()->SetBinLabel(14,"fUSAllpt");  //!Global Trigger UnlikeSign pair All p_T
  fHEventStat->GetXaxis()->SetBinLabel(15,"fLSLowpt");  //!Global Trigger LikeSign pair pair Low p_T
  fHEventStat->GetXaxis()->SetBinLabel(16,"fLSHighpt"); //!Global Trigger LikeSign pair pair High p_T
  fHEventStat->GetXaxis()->SetBinLabel(17,"fLSAllpt");  //!Global Trigger LikeSign pair pair All p_T
  fHEventStat->GetXaxis()->SetBinLabel(18,"fSPLowpt");   //!Global Trigger Single plus Low p_T
  fOutputList->Add(fHEventStat);

  fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", 30, -15, 15, 103, -2, 101);
  fOutputList->Add(fHEvt);


  Int_t iTrackBin[5];
  Double_t* trackBins[5];
  const char* trackAxisTitle[5];

  // eta
  Double_t etaBins[fNEtaBins+1];
  for (Int_t i=0; i<=fNEtaBins; i++) { etaBins[i] = (Double_t)(-4.0 + 1.5/fNEtaBins*i); }
  iTrackBin[0] = fNEtaBins;
  trackBins[0] = etaBins;
  trackAxisTitle[0] = "#eta";

  // p_T
  Double_t pTBins[fNpTBins+1];
  for (Int_t i=0; i<=fNpTBins; i++) { pTBins[i] = (Double_t)(10.0/fNpTBins*i); }
  iTrackBin[1] = fNpTBins;
  trackBins[1] = pTBins;
  trackAxisTitle[1] = "p_{T} (GeV/c)";

  // centrality
  Double_t CentBins[fNCentBins+1];
  for (Int_t i=0; i<=fNCentBins; i++) { CentBins[i] = (Double_t)(100.0/fNCentBins*i); }
  iTrackBin[2] = fNCentBins;
  trackBins[2] = CentBins;
  trackAxisTitle[2] = "Cent";

  // Z-vertex
  Double_t ZvtxBins[fNZvtxBins+1];
  for (Int_t i=0; i<=fNZvtxBins; i++) { ZvtxBins[i] = (Double_t)(-10.0 + 20.0/fNZvtxBins*i); }
  iTrackBin[3] = fNZvtxBins;
  trackBins[3] = ZvtxBins;
  trackAxisTitle[3] = "Zvtx";

  // phi
  Double_t phiBins[fNPhiBins+1];
  for (Int_t i=0; i<=fNPhiBins; i++) { phiBins[i] = (Double_t)(TMath::TwoPi()/fNPhiBins * i); }
  iTrackBin[4] = fNPhiBins;
  trackBins[4] = phiBins;
  trackAxisTitle[4] = "#phi";

  // THn for tracking efficiency
  fHMuonParGen = new THnF("fHMuonParGen", "", 5, iTrackBin, 0, 0);
  for (Int_t i=0; i<5; i++)
  {
    fHMuonParGen->SetBinEdges(i, trackBins[i]);
    fHMuonParGen->GetAxis(i)->SetTitle(trackAxisTitle[i]);
  }
  fHMuonParGen->Sumw2();
  fOutputList->Add(fHMuonParGen);

  fHMuonDetGen = (THnF*) fHMuonParGen->Clone("fHMuonDetGen");
  fHMuonDetGen->Sumw2();
  fOutputList->Add(fHMuonDetGen);

  fHMuonDetRec = (THnF*) fHMuonParGen->Clone("fHMuonDetRec");
  fHMuonDetRec->Sumw2();
  fOutputList->Add(fHMuonDetRec);

  fHEtcDetRec = (THnF*) fHMuonParGen->Clone("fHEtcDetRec");
  fHEtcDetRec->Sumw2();
  fOutputList->Add(fHEtcDetRec);

  if(fMDProcess)
  {
    const char* MotherSpecies[4] = {"Pion","Kaon","D", "Etc"};
    
    for(Int_t i=0; i<4; i++)
    {
      fHMuMotherGenPt[i] = new TH2F(Form("fHMuMotherGenPt_%s",MotherSpecies[i]),";p_{T,muon}^{gen} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
      fOutputList->Add(fHMuMotherGenPt[i]);
      
      fHMuMotherRecPt[i] = new TH2F(Form("fHMuMotherRecPt_%s",MotherSpecies[i]),";p_{T,muon}^{rec} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
      fOutputList->Add(fHMuMotherRecPt[i]);
      
      fHMuMotherGenPhi[i] = new TH2F(Form("fHMuMotherGenPhi_%s",MotherSpecies[i]),";#phi_{gen};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
      fOutputList->Add(fHMuMotherGenPhi[i]);
      
      fHMuMotherRecPhi[i] = new TH2F(Form("fHMuMotherRecPhi_%s",MotherSpecies[i]),";#phi_{rec};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
      fOutputList->Add(fHMuMotherRecPhi[i]);
      
      fHMuMotherGenEta[i] = new TH2F(Form("fHMuMotherGenEta_%s",MotherSpecies[i]),";#eta_{gen};mother #eta;",100, -5., -1., 100, -5., -1.);
      fOutputList->Add(fHMuMotherGenEta[i]);
      
      fHMuMotherRecEta[i] = new TH2F(Form("fHMuMotherRecEta_%s",MotherSpecies[i]),";#eta_{rec};mother #eta;",100, -5., -1., 100, -5., -1.);
      fOutputList->Add(fHMuMotherRecEta[i]);
      
      fHMuDCA[i] =  new TH1F(Form("fHMuDCA_%s",MotherSpecies[i]), ";DCA", 100, 0, 50);
      fOutputList->Add(fHMuDCA[i]);
    }
  }
  PostData(1, fOutputList);
}

//________________________________________________________________________
void AliMuonEffMC::UserExec(Option_t *) 
{
  // Main loop, Called for each event
  Int_t ntrks = 0; // number of tracks in an event

  if(((TString)InputEvent()->IsA()->GetName())=="AliAODEvent")
  {
    fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
    if (!fAOD) { AliError("AOD event not found. Nothing done!"); return; }
    ntrks = fAOD->GetNTracks();
  } 
  else
  {
    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
    if (!fESD) { AliError("ESD event not found. Nothing done!"); return; }
    ntrks = fESD->GetNumberOfMuonTracks();
  }
	
  fHEventStat->Fill(0.5);
  if(fIsMc) 
  {
    fMC = MCEvent();
    if (!fMC) { AliError("MC event not avaliable."); return; }
  }

  // Centrality, vertex, other event variables...
  if(fAOD)
  { 
    const AliAODVertex* vertex = fAOD->GetPrimaryVertex();
    fZVertex = vertex->GetZ();
    if(fAOD->GetCentrality())  fCentrality = fAOD->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
  }
  else if(fESD) 
  {
    const AliESDVertex* vertex = fESD->GetPrimaryVertex();
    fZVertex = vertex->GetZ();
    if(fESD->GetCentrality()) fCentrality = fESD->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
  }   

  if ((fESD && !VertexOk(fESD)) || (fAOD && !VertexOk(fAOD))) { AliInfo(Form("Event REJECTED. z = %.1f", fZVertex)); return; }
  if (fCentrality > 100. || fCentrality < -1.5) { AliInfo(Form("Event REJECTED. fCentrality = %.1f", fCentrality)); return; }
 
  if(fCentrality < 0) fCentrality = 0.5; //ad hoc centrality for pp
  // Fill Event histogram
  fHEvt->Fill(fZVertex, fCentrality);
  fHEventStat->Fill(1.5);
 
  ULong64_t trigword = 0;
  if(fAOD) trigword=fAOD->GetTriggerMask();
  else if(fESD) trigword=fESD->GetTriggerMask();
 
  if (trigword & 0x01) fHEventStat->Fill(17.5);
  if (trigword & 0x02) fHEventStat->Fill(3.5);
  if (trigword & 0x04) fHEventStat->Fill(4.5);
  if (trigword & 0x08) fHEventStat->Fill(5.5);      
  if (trigword & 0x010) fHEventStat->Fill(6.5);
  if (trigword & 0x020) fHEventStat->Fill(7.5);
  if (trigword & 0x040) fHEventStat->Fill(8.5);
  if (trigword & 0x080) fHEventStat->Fill(9.5);
  if (trigword & 0x100) fHEventStat->Fill(10.5);
  if (trigword & 0x200) fHEventStat->Fill(11.5);
  if (trigword & 0x400) fHEventStat->Fill(12.5);
  if (trigword & 0x800) fHEventStat->Fill(13.5);
  if (trigword & 0x1000) fHEventStat->Fill(14.5);
  if (trigword & 0x2000) fHEventStat->Fill(15.5);
  if (trigword & 0x4000) fHEventStat->Fill(16.5);
  
  // generated level loop
  for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
  {
    if(fAOD)
    {
      AliAODMCParticle *AodMcParticle  = (AliAODMCParticle*)fMC->GetTrack(ipart);
      if(TMath::Abs(AodMcParticle->PdgCode())==13) 
      {
	Double_t fillArrayParGen[5] = { AodMcParticle->Eta(), AodMcParticle->Pt(), fCentrality, fZVertex, AodMcParticle->Phi() };
	fHMuonParGen->Fill(fillArrayParGen);
      }
    } 
    else if(fESD)
    {
      AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(ipart);
      if(TMath::Abs(McParticle->PdgCode())==13) 
      {
	Double_t fillArrayParGen[5] = { McParticle->Eta(), McParticle->Pt(), fCentrality, fZVertex, McParticle->Phi() };
	fHMuonParGen->Fill(fillArrayParGen);
      }
    }
  }
  
  // reconstructed level loop
  for (Int_t iTrack = 0; iTrack<ntrks; iTrack++) 
  { 
    Int_t label = 0;
    Double_t trackpt = 0;
    Double_t tracketa = 0;
    Double_t trackphi = 0;
    Double_t dcavalue = 0;
    
    Double_t mcpt = 0;
    Double_t mceta = 0;
    Double_t mcphi = 0;
    
    Int_t motherlabel =0;
    Int_t motherpdg = 0;
    Double_t motherpt = 0;
    Double_t mothereta = 0;
    Double_t motherphi = 0;
    
    if(fAOD)
    {
      AliAODTrack* muonTrack = (AliAODTrack*)fAOD->GetTrack(iTrack);
      if(muonTrack)
      {
	if(!(IsGoodMUONtrack(*muonTrack)) || !(muonTrack->IsMuonTrack())) continue;
	trackpt = muonTrack->Pt();
	tracketa = muonTrack->Eta();
	trackphi = muonTrack->Phi();
	label =  TMath::Abs(muonTrack->GetLabel());
      }
    }
    else if(fESD)
    {
      AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iTrack);
      if(muonTrack)
      {
	if(!IsGoodMUONtrack(*muonTrack)) continue;
	trackpt = muonTrack->Pt();
	tracketa = muonTrack->Eta();
	trackphi = muonTrack->Phi();
	label =  TMath::Abs(muonTrack->GetLabel());
	dcavalue = muonTrack->GetDCA();
      }
    }
    
    Double_t fillArrayDetRec[5] = { tracketa, trackpt, fCentrality, fZVertex, trackphi };
    fHMuonDetRec->Fill(fillArrayDetRec); 
    
    if(fAOD)
    {
      AliAODMCParticle *aodMcParticle  = (AliAODMCParticle*)fMC->GetTrack(label);
      if(TMath::Abs(aodMcParticle->PdgCode()) != 13) 
      {
	fHEtcDetRec->Fill(fillArrayDetRec); 
	continue;
      }
      mcpt = aodMcParticle->Pt();
      mceta = aodMcParticle->Eta();
      mcphi = aodMcParticle->Phi();
      motherlabel = aodMcParticle->GetMother();
      if(motherlabel > 0) 
      { 
	AliAODMCParticle *aodMotherParticle  = (AliAODMCParticle*)fMC->GetTrack(motherlabel);
	motherpdg = TMath::Abs(aodMotherParticle->PdgCode());
	motherpt = aodMotherParticle->Pt();
	mothereta = aodMotherParticle->Eta();
	motherphi = aodMotherParticle->Phi();
      }      
    }
    else if(fESD)
    {
      AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(label);
      if(TMath::Abs(McParticle->PdgCode()) != 13) 
      {
	fHEtcDetRec->Fill(fillArrayDetRec); 
	continue;
      }
      mcpt = McParticle->Pt();
      mceta = McParticle->Eta();
      mcphi = McParticle->Phi();
      motherlabel = McParticle->GetMother();
      if(motherlabel > 0) 
      { 
	AliMCParticle *MotherParticle  = (AliMCParticle*)fMC->GetTrack(motherlabel);
	motherpdg = TMath::Abs(MotherParticle->PdgCode());
	motherpt = MotherParticle->Pt();
	mothereta = MotherParticle->Eta();
	motherphi = MotherParticle->Phi();
      }
    }
    Double_t fillArrayDetGen[5] = { mceta, mcpt, fCentrality, fZVertex, mcphi };
    fHMuonDetGen->Fill(fillArrayDetGen);
    
    // mother-daughter kinematic relation
    if(fMDProcess)
    {
      if(motherlabel > 0)
      {
	if(motherpdg==411 || motherpdg==413 || motherpdg==421 || motherpdg==423 || motherpdg==431 || motherpdg==433 || motherpdg==10413 || motherpdg==10411 || motherpdg==10423 || motherpdg==10421 || motherpdg==10433 || motherpdg==10431 || motherpdg==20413 || motherpdg==415 || motherpdg==20423 || motherpdg==425 || motherpdg==20433 || motherpdg==435)
	{  
	  fHMuMotherGenPt[2]->Fill(mcpt, motherpt);
	  fHMuMotherRecPt[2]->Fill(trackpt, motherpt);
	  fHMuMotherGenPhi[2]->Fill(mcphi, motherphi);
	  fHMuMotherRecPhi[2]->Fill(trackphi, motherphi);
	  fHMuMotherGenEta[2]->Fill(mceta, mothereta);
	  fHMuMotherRecEta[2]->Fill(tracketa, mothereta);
	  if(fESD) fHMuDCA[2]->Fill(dcavalue);
	}
	else if(motherpdg==211) 
	{
	  fHMuMotherGenPt[0]->Fill(mcpt, motherpt);
	  fHMuMotherRecPt[0]->Fill(trackpt, motherpt);
	  fHMuMotherGenPhi[0]->Fill(mcphi, motherphi);
	  fHMuMotherRecPhi[0]->Fill(trackphi, motherphi);
	  fHMuMotherGenEta[0]->Fill(mceta, mothereta);
	  fHMuMotherRecEta[0]->Fill(tracketa, mothereta);
	  if(fESD) fHMuDCA[0]->Fill(dcavalue);
	}
	else if(motherpdg==321)
	{
	  fHMuMotherGenPt[1]->Fill(mcpt, motherpt);
	  fHMuMotherRecPt[1]->Fill(trackpt, motherpt);
	  fHMuMotherGenPhi[1]->Fill(mcphi, motherphi);
	  fHMuMotherRecPhi[1]->Fill(trackphi, motherphi);
	  fHMuMotherGenEta[1]->Fill(mceta, mothereta);
	  fHMuMotherRecEta[1]->Fill(tracketa, mothereta);
	  if(fESD) fHMuDCA[1]->Fill(dcavalue);
	}	
	else 
	{
	  fHMuMotherGenPt[3]->Fill(mcpt, motherpt);
	  fHMuMotherRecPt[3]->Fill(trackpt, motherpt);
	  fHMuMotherGenPhi[3]->Fill(mcphi, motherphi);
	  fHMuMotherRecPhi[3]->Fill(trackphi, motherphi);
	  fHMuMotherGenEta[3]->Fill(mceta, mothereta);
	  fHMuMotherRecEta[3]->Fill(tracketa, mothereta);
	  if(fESD) fHMuDCA[3]->Fill(dcavalue);
	}
      }    
    }	 // (mother hadron) : (daughter muon) QA */
  } 
  
  PostData(1, fOutputList);
  return;
}

//________________________________________________________________________
void AliMuonEffMC::Terminate(Option_t *) 
{
  // Draw result to the screen
  // Called once at the end of the query

}

//________________________________________________________________________
Bool_t AliMuonEffMC::VertexOk(TObject* obj) const
{
  // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
  
  Int_t nContributors  = 0;
  Double_t zVertex     = 999;
  TString name("");
  
  if (obj->InheritsFrom("AliESDEvent")) {
    AliESDEvent* esdevt = (AliESDEvent*) obj;
    const AliESDVertex* vtx = esdevt->GetPrimaryVertex();
    if (!vtx)
      return 0;
    nContributors = vtx->GetNContributors();
    zVertex       = vtx->GetZ();
    name          = vtx->GetName();
  }
  else if (obj->InheritsFrom("AliAODEvent")) { 
    AliAODEvent* aodevt = (AliAODEvent*) obj;
    if (aodevt->GetNumberOfVertices() < 1)
      return 0;
    const AliAODVertex* vtx = aodevt->GetPrimaryVertex();
    nContributors = vtx->GetNContributors();
    zVertex       = vtx->GetZ();
    name          = vtx->GetName();
  }
  
  // Reject if TPC-only vertex
  if (name.CompareTo("TPCVertex")==0)
    return kFALSE;
  
  // Check # contributors and range...
  if( nContributors < 1 || TMath::Abs(zVertex) > 10 ) {
    return kFALSE;
  }
  
  return kTRUE;
}


//________________________________________________________________________
Bool_t AliMuonEffMC::IsGoodMUONtrack(AliESDMuonTrack &track)
{
  // Applying track cuts for MUON tracks
  if(!track.ContainTrackerData()) return kFALSE;

  Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
  Double_t eta = track.Eta();

  // Theta cut at absorber end
  if(thetaTrackAbsEnd <= 2. || thetaTrackAbsEnd >= 10.) return kFALSE;
  // Eta cut
  if(eta <= -4. || eta >= -2.5) return kFALSE;
  if(track.GetMatchTrigger() <= 0) return kFALSE;

  return kTRUE;
}

//________________________________________________________________________
Bool_t AliMuonEffMC::IsGoodMUONtrack(AliAODTrack &track)
{
  if (!track.IsMuonTrack()) return kFALSE;

  Double_t dThetaAbs = TMath::ATan(track.GetRAtAbsorberEnd()/505.)
                     * TMath::RadToDeg();
  if ((dThetaAbs<2.) || (dThetaAbs>10.)) return kFALSE;

  Double_t dEta = track.Eta();
  if ((dEta<-4.) || (dEta>2.5)) return kFALSE;

  if (track.GetMatchTrigger()<0.5) return kFALSE;

  return kTRUE;
}
Commit	Line	Data
7e990b20	1	// MUON track QA referring AliMuonEffMC.cxx
	2	// Author : Saehanseul Oh
	3
	4	#include "AliMuonEffMC.h"
	5
	6	#include <TList.h>
	7	#include <TH1D.h>
	8	#include <TH2F.h>
	9	#include <TH3F.h>
	10	#include <THn.h>
	11	#include <TChain.h>
	12	#include <TFile.h>
	13
	14	#include "AliAnalysisManager.h"
	15	#include "AliAnalysisTask.h"
	16	#include "AliESDVertex.h"
	17	#include "AliESDEvent.h"
	18	#include "AliAODEvent.h"
	19	#include "AliESDMuonTrack.h"
b673a083	20	#include "AliAODTrack.h"
7e990b20	21	#include "AliESDVertex.h"
	22	#include "AliAODVertex.h"
	23	#include "AliCentrality.h"
	24	#include "AliVParticle.h"
	25	#include "AliMCParticle.h"
	26	#include "AliAnalysisHelperJetTasks.h"
	27	#include "AliMCEvent.h"
fe76c976	28	#include "AliAODMCParticle.h"
7e990b20	29
	30	using std::cout;
	31	using std::endl;
	32
	33	ClassImp(AliMuonEffMC)
	34
	35	//________________________________________________________________________
	36	AliMuonEffMC::AliMuonEffMC() :
	37	AliAnalysisTaskSE(), fESD(0), fAOD(0), fMC(0), fCentrality(99), fZVertex(99), fOutputList(0x0),
	38	fHEventStat(0), fHEvt(0x0), fIsMc(kTRUE), fMDProcess(kTRUE), fCentralityEstimator("V0M"),
	39	fNEtaBins(15), fNpTBins(100), fNCentBins(5), fNZvtxBins(10), fNPhiBins(24),
	40	fHMuonParGen(0x0), fHMuonDetGen(0x0), fHMuonDetRec(0x0), fHEtcDetRec(0x0)
	41	{
	42	// Constructor
fe76c976	43
7e990b20	44	}
	45
	46	//________________________________________________________________________
	47	AliMuonEffMC::AliMuonEffMC(const char *name) :
	48	AliAnalysisTaskSE(name), fESD(0), fAOD(0), fMC(0), fCentrality(99), fZVertex(99), fOutputList(0x0),
	49	fHEventStat(0), fHEvt(0x0), fIsMc(kTRUE), fMDProcess(kTRUE), fCentralityEstimator("V0M"),
	50	fNEtaBins(15), fNpTBins(100), fNCentBins(5), fNZvtxBins(10), fNPhiBins(24),
	51	fHMuonParGen(0x0), fHMuonDetGen(0x0), fHMuonDetRec(0x0), fHEtcDetRec(0x0)
	52	{
	53	// Constructor
	54	DefineInput(0, TChain::Class());
	55	DefineOutput(1, TList::Class());
fe76c976	56
7e990b20	57	}
	58
	59	//________________________________________________________________________
	60	AliMuonEffMC::~AliMuonEffMC()
	61	{
	62	//Destructor
9370528f	63	if(fOutputList) delete fOutputList;
7e990b20	64	}
	65
	66	//________________________________________________________________________
	67	void AliMuonEffMC::UserCreateOutputObjects()
	68	{
	69	// Create histograms
	70	// Called once (per slave on PROOF!)
fe76c976	71
	72	fOutputList = new TList();
	73	fOutputList->SetOwner(1);
	74
53fc502b	75	for(Int_t i=0; i<4; i++)
	76	{
	77	fHMuMotherGenPt[i] = 0x0;
	78	fHMuMotherRecPt[i] = 0x0;
	79	fHMuMotherGenPhi[i] = 0x0;
	80	fHMuMotherRecPhi[i] = 0x0;
	81	fHMuMotherGenEta[i] = 0x0;
	82	fHMuMotherRecEta[i] = 0x0;
	83	fHMuDCA[i] = 0x0;
	84	}
7e990b20	85
7e990b20	86	fHEventStat = new TH1D("fHEventStat","Event statistics for analysis",18,0,18);
	87	fHEventStat->GetXaxis()->SetBinLabel(1,"Event");
	88	fHEventStat->GetXaxis()->SetBinLabel(2,"SelectedEvent");
	89	fHEventStat->GetXaxis()->SetBinLabel(3,"File");
	90	fHEventStat->GetXaxis()->SetBinLabel(4,"fSPHighpt"); //!Global Trigger Single plus High p_T
	91	fHEventStat->GetXaxis()->SetBinLabel(5,"fSPAllpt"); //!Global Trigger Single plus All p_T
	92	fHEventStat->GetXaxis()->SetBinLabel(6,"fSMLowpt"); //!Global Trigger Single minus Low p_T
	93	fHEventStat->GetXaxis()->SetBinLabel(7,"fSMHighpt"); //!Global Trigger Single minus High p_T
	94	fHEventStat->GetXaxis()->SetBinLabel(8,"fSMAllpt"); //!Global Trigger Single minus All p_T
	95	fHEventStat->GetXaxis()->SetBinLabel(9,"fSULowpt"); //!Global Trigger Single undefined Low p_T
	96	fHEventStat->GetXaxis()->SetBinLabel(10,"fSUHighpt"); //!Global Trigger Single undefined High p_T
	97	fHEventStat->GetXaxis()->SetBinLabel(11,"fSUAllpt"); //!Global Trigger Single undefined All p_T
	98	fHEventStat->GetXaxis()->SetBinLabel(12,"fUSLowpt"); //!Global Trigger UnlikeSign pair Low p_T
	99	fHEventStat->GetXaxis()->SetBinLabel(13,"fUSHighpt"); //!Global Trigger UnlikeSign pair High p_T
	100	fHEventStat->GetXaxis()->SetBinLabel(14,"fUSAllpt"); //!Global Trigger UnlikeSign pair All p_T
	101	fHEventStat->GetXaxis()->SetBinLabel(15,"fLSLowpt"); //!Global Trigger LikeSign pair pair Low p_T
	102	fHEventStat->GetXaxis()->SetBinLabel(16,"fLSHighpt"); //!Global Trigger LikeSign pair pair High p_T
	103	fHEventStat->GetXaxis()->SetBinLabel(17,"fLSAllpt"); //!Global Trigger LikeSign pair pair All p_T
	104	fHEventStat->GetXaxis()->SetBinLabel(18,"fSPLowpt"); //!Global Trigger Single plus Low p_T
	105	fOutputList->Add(fHEventStat);
	106
	107	fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", 30, -15, 15, 103, -2, 101);
	108	fOutputList->Add(fHEvt);
	109
	110
	111	Int_t iTrackBin[5];
	112	Double_t* trackBins[5];
	113	const char* trackAxisTitle[5];
	114
	115	// eta
	116	Double_t etaBins[fNEtaBins+1];
	117	for (Int_t i=0; i<=fNEtaBins; i++) { etaBins[i] = (Double_t)(-4.0 + 1.5/fNEtaBins*i); }
	118	iTrackBin[0] = fNEtaBins;
	119	trackBins[0] = etaBins;
	120	trackAxisTitle[0] = "#eta";
	121
	122	// p_T
	123	Double_t pTBins[fNpTBins+1];
	124	for (Int_t i=0; i<=fNpTBins; i++) { pTBins[i] = (Double_t)(10.0/fNpTBins*i); }
	125	iTrackBin[1] = fNpTBins;
	126	trackBins[1] = pTBins;
	127	trackAxisTitle[1] = "p_{T} (GeV/c)";
	128
	129	// centrality
	130	Double_t CentBins[fNCentBins+1];
	131	for (Int_t i=0; i<=fNCentBins; i++) { CentBins[i] = (Double_t)(100.0/fNCentBins*i); }
	132	iTrackBin[2] = fNCentBins;
	133	trackBins[2] = CentBins;
	134	trackAxisTitle[2] = "Cent";
	135
	136	// Z-vertex
	137	Double_t ZvtxBins[fNZvtxBins+1];
	138	for (Int_t i=0; i<=fNZvtxBins; i++) { ZvtxBins[i] = (Double_t)(-10.0 + 20.0/fNZvtxBins*i); }
	139	iTrackBin[3] = fNZvtxBins;
	140	trackBins[3] = ZvtxBins;
	141	trackAxisTitle[3] = "Zvtx";
	142
	143	// phi
	144	Double_t phiBins[fNPhiBins+1];
	145	for (Int_t i=0; i<=fNPhiBins; i++) { phiBins[i] = (Double_t)(TMath::TwoPi()/fNPhiBins * i); }
	146	iTrackBin[4] = fNPhiBins;
	147	trackBins[4] = phiBins;
	148	trackAxisTitle[4] = "#phi";
	149
150	// THn for tracking efficiency
151	fHMuonParGen = new THnF("fHMuonParGen", "", 5, iTrackBin, 0, 0);
152	for (Int_t i=0; i<5; i++)
153	{
154	fHMuonParGen->SetBinEdges(i, trackBins[i]);
155	fHMuonParGen->GetAxis(i)->SetTitle(trackAxisTitle[i]);
156	}
157	fHMuonParGen->Sumw2();
158	fOutputList->Add(fHMuonParGen);
159
160	fHMuonDetGen = (THnF*) fHMuonParGen->Clone("fHMuonDetGen");
161	fHMuonDetGen->Sumw2();
162	fOutputList->Add(fHMuonDetGen);
163
164	fHMuonDetRec = (THnF*) fHMuonParGen->Clone("fHMuonDetRec");
165	fHMuonDetRec->Sumw2();
166	fOutputList->Add(fHMuonDetRec);
167
168	fHEtcDetRec = (THnF*) fHMuonParGen->Clone("fHEtcDetRec");
169	fHEtcDetRec->Sumw2();
170	fOutputList->Add(fHEtcDetRec);
171
172	if(fMDProcess)
173	{
174	const char* MotherSpecies[4] = {"Pion","Kaon","D", "Etc"};
175
176	for(Int_t i=0; i<4; i++)
177	{
178	fHMuMotherGenPt[i] = new TH2F(Form("fHMuMotherGenPt_%s",MotherSpecies[i]),";p_{T,muon}^{gen} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
179	fOutputList->Add(fHMuMotherGenPt[i]);
180
181	fHMuMotherRecPt[i] = new TH2F(Form("fHMuMotherRecPt_%s",MotherSpecies[i]),";p_{T,muon}^{rec} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
182	fOutputList->Add(fHMuMotherRecPt[i]);
183
184	fHMuMotherGenPhi[i] = new TH2F(Form("fHMuMotherGenPhi_%s",MotherSpecies[i]),";#phi_{gen};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
185	fOutputList->Add(fHMuMotherGenPhi[i]);
186
187	fHMuMotherRecPhi[i] = new TH2F(Form("fHMuMotherRecPhi_%s",MotherSpecies[i]),";#phi_{rec};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
188	fOutputList->Add(fHMuMotherRecPhi[i]);
189
190	fHMuMotherGenEta[i] = new TH2F(Form("fHMuMotherGenEta_%s",MotherSpecies[i]),";#eta_{gen};mother #eta;",100, -5., -1., 100, -5., -1.);
191	fOutputList->Add(fHMuMotherGenEta[i]);
192
193	fHMuMotherRecEta[i] = new TH2F(Form("fHMuMotherRecEta_%s",MotherSpecies[i]),";#eta_{rec};mother #eta;",100, -5., -1., 100, -5., -1.);
194	fOutputList->Add(fHMuMotherRecEta[i]);
195
196	fHMuDCA[i] = new TH1F(Form("fHMuDCA_%s",MotherSpecies[i]), ";DCA", 100, 0, 50);
197	fOutputList->Add(fHMuDCA[i]);
198	}
199	}
200	PostData(1, fOutputList);
201	}
202
203	//________________________________________________________________________
204	void AliMuonEffMC::UserExec(Option_t *)
205	{
206	// Main loop, Called for each event
b673a083	207	Int_t ntrks = 0; // number of tracks in an event
	208
	209	if(((TString)InputEvent()->IsA()->GetName())=="AliAODEvent")
7e990b20	210	{
7e990b20	211	fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
b673a083	212	if (!fAOD) { AliError("AOD event not found. Nothing done!"); return; }
	213	ntrks = fAOD->GetNTracks();
	214	}
	215	else
7e990b20	216	{
b673a083	217	fESD = dynamic_cast<AliESDEvent*>(InputEvent());
	218	if (!fESD) { AliError("ESD event not found. Nothing done!"); return; }
	219	ntrks = fESD->GetNumberOfMuonTracks();
7e990b20	220	}
	221
	222	fHEventStat->Fill(0.5);
	223	if(fIsMc)
	224	{
	225	fMC = MCEvent();
b673a083	226	if (!fMC) { AliError("MC event not avaliable."); return; }
7e990b20	227	}
7e990b20	228
b673a083	229	// Centrality, vertex, other event variables...
	230	if(fAOD)
	231	{
	232	const AliAODVertex* vertex = fAOD->GetPrimaryVertex();
	233	fZVertex = vertex->GetZ();
fe76c976	234	if(fAOD->GetCentrality()) fCentrality = fAOD->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
7e990b20	235	}
b673a083	236	else if(fESD)
	237	{
	238	const AliESDVertex* vertex = fESD->GetPrimaryVertex();
	239	fZVertex = vertex->GetZ();
fe76c976	240	if(fESD->GetCentrality()) fCentrality = fESD->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
b673a083	241	}
7e990b20	242
b673a083	243	if ((fESD && !VertexOk(fESD)) \|\| (fAOD && !VertexOk(fAOD))) { AliInfo(Form("Event REJECTED. z = %.1f", fZVertex)); return; }
b673a083	244	if (fCentrality > 100. \|\| fCentrality < -1.5) { AliInfo(Form("Event REJECTED. fCentrality = %.1f", fCentrality)); return; }
7e990b20	245
fe76c976	246	if(fCentrality < 0) fCentrality = 0.5; //ad hoc centrality for pp
b673a083	247	// Fill Event histogram
b673a083	248	fHEvt->Fill(fZVertex, fCentrality);
7e990b20	249	fHEventStat->Fill(1.5);
7e990b20	250
b673a083	251	ULong64_t trigword = 0;
	252	if(fAOD) trigword=fAOD->GetTriggerMask();
	253	else if(fESD) trigword=fESD->GetTriggerMask();
7e990b20	254
	255	if (trigword & 0x01) fHEventStat->Fill(17.5);
	256	if (trigword & 0x02) fHEventStat->Fill(3.5);
	257	if (trigword & 0x04) fHEventStat->Fill(4.5);
	258	if (trigword & 0x08) fHEventStat->Fill(5.5);
	259	if (trigword & 0x010) fHEventStat->Fill(6.5);
	260	if (trigword & 0x020) fHEventStat->Fill(7.5);
	261	if (trigword & 0x040) fHEventStat->Fill(8.5);
	262	if (trigword & 0x080) fHEventStat->Fill(9.5);
	263	if (trigword & 0x100) fHEventStat->Fill(10.5);
	264	if (trigword & 0x200) fHEventStat->Fill(11.5);
	265	if (trigword & 0x400) fHEventStat->Fill(12.5);
	266	if (trigword & 0x800) fHEventStat->Fill(13.5);
	267	if (trigword & 0x1000) fHEventStat->Fill(14.5);
	268	if (trigword & 0x2000) fHEventStat->Fill(15.5);
	269	if (trigword & 0x4000) fHEventStat->Fill(16.5);
fe76c976	270
	271	// generated level loop
	272	for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
7e990b20	273	{
fe76c976	274	if(fAOD)
	275	{
	276	AliAODMCParticle AodMcParticle = (AliAODMCParticle)fMC->GetTrack(ipart);
	277	if(TMath::Abs(AodMcParticle->PdgCode())==13)
7e990b20	278	{
fe76c976	279	Double_t fillArrayParGen[5] = { AodMcParticle->Eta(), AodMcParticle->Pt(), fCentrality, fZVertex, AodMcParticle->Phi() };
	280	fHMuonParGen->Fill(fillArrayParGen);
	281	}
	282	}
	283	else if(fESD)
	284	{
	285	AliMCParticle McParticle = (AliMCParticle)fMC->GetTrack(ipart);
	286	if(TMath::Abs(McParticle->PdgCode())==13)
	287	{
	288	Double_t fillArrayParGen[5] = { McParticle->Eta(), McParticle->Pt(), fCentrality, fZVertex, McParticle->Phi() };
	289	fHMuonParGen->Fill(fillArrayParGen);
b673a083	290	}
fe76c976	291	}
	292	}
	293
	294	// reconstructed level loop
	295	for (Int_t iTrack = 0; iTrack<ntrks; iTrack++)
	296	{
	297	Int_t label = 0;
	298	Double_t trackpt = 0;
	299	Double_t tracketa = 0;
	300	Double_t trackphi = 0;
	301	Double_t dcavalue = 0;
	302
	303	Double_t mcpt = 0;
	304	Double_t mceta = 0;
	305	Double_t mcphi = 0;
	306
	307	Int_t motherlabel =0;
	308	Int_t motherpdg = 0;
	309	Double_t motherpt = 0;
	310	Double_t mothereta = 0;
	311	Double_t motherphi = 0;
	312
	313	if(fAOD)
	314	{
	315	AliAODTrack* muonTrack = (AliAODTrack*)fAOD->GetTrack(iTrack);
	316	if(muonTrack)
b673a083	317	{
fe76c976	318	if(!(IsGoodMUONtrack(*muonTrack)) \|\| !(muonTrack->IsMuonTrack())) continue;
	319	trackpt = muonTrack->Pt();
	320	tracketa = muonTrack->Eta();
	321	trackphi = muonTrack->Phi();
	322	label = TMath::Abs(muonTrack->GetLabel());
b673a083	323	}
fe76c976	324	}
	325	else if(fESD)
	326	{
	327	AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iTrack);
	328	if(muonTrack)
	329	{
	330	if(!IsGoodMUONtrack(*muonTrack)) continue;
	331	trackpt = muonTrack->Pt();
	332	tracketa = muonTrack->Eta();
	333	trackphi = muonTrack->Phi();
	334	label = TMath::Abs(muonTrack->GetLabel());
	335	dcavalue = muonTrack->GetDCA();
	336	}
	337	}
	338
	339	Double_t fillArrayDetRec[5] = { tracketa, trackpt, fCentrality, fZVertex, trackphi };
	340	fHMuonDetRec->Fill(fillArrayDetRec);
	341
	342	if(fAOD)
	343	{
	344	AliAODMCParticle aodMcParticle = (AliAODMCParticle)fMC->GetTrack(label);
	345	if(TMath::Abs(aodMcParticle->PdgCode()) != 13)
	346	{
	347	fHEtcDetRec->Fill(fillArrayDetRec);
	348	continue;
	349	}
	350	mcpt = aodMcParticle->Pt();
	351	mceta = aodMcParticle->Eta();
	352	mcphi = aodMcParticle->Phi();
	353	motherlabel = aodMcParticle->GetMother();
	354	if(motherlabel > 0)
	355	{
	356	AliAODMCParticle aodMotherParticle = (AliAODMCParticle)fMC->GetTrack(motherlabel);
	357	motherpdg = TMath::Abs(aodMotherParticle->PdgCode());
	358	motherpt = aodMotherParticle->Pt();
	359	mothereta = aodMotherParticle->Eta();
	360	motherphi = aodMotherParticle->Phi();
	361	}
	362	}
	363	else if(fESD)
	364	{
b673a083	365	AliMCParticle McParticle = (AliMCParticle)fMC->GetTrack(label);
	366	if(TMath::Abs(McParticle->PdgCode()) != 13)
	367	{
	368	fHEtcDetRec->Fill(fillArrayDetRec);
	369	continue;
	370	}
fe76c976	371	mcpt = McParticle->Pt();
	372	mceta = McParticle->Eta();
	373	mcphi = McParticle->Phi();
	374	motherlabel = McParticle->GetMother();
	375	if(motherlabel > 0)
	376	{
	377	AliMCParticle MotherParticle = (AliMCParticle)fMC->GetTrack(motherlabel);
	378	motherpdg = TMath::Abs(MotherParticle->PdgCode());
	379	motherpt = MotherParticle->Pt();
	380	mothereta = MotherParticle->Eta();
	381	motherphi = MotherParticle->Phi();
	382	}
7e990b20	383	}
fe76c976	384	Double_t fillArrayDetGen[5] = { mceta, mcpt, fCentrality, fZVertex, mcphi };
fe76c976	385	fHMuonDetGen->Fill(fillArrayDetGen);
7e990b20	386
fe76c976	387	// mother-daughter kinematic relation
fe76c976	388	if(fMDProcess)
7e990b20	389	{
fe76c976	390	if(motherlabel > 0)
7e990b20	391	{
fe76c976	392	if(motherpdg==411 \|\| motherpdg==413 \|\| motherpdg==421 \|\| motherpdg==423 \|\| motherpdg==431 \|\| motherpdg==433 \|\| motherpdg==10413 \|\| motherpdg==10411 \|\| motherpdg==10423 \|\| motherpdg==10421 \|\| motherpdg==10433 \|\| motherpdg==10431 \|\| motherpdg==20413 \|\| motherpdg==415 \|\| motherpdg==20423 \|\| motherpdg==425 \|\| motherpdg==20433 \|\| motherpdg==435)
	393	{
	394	fHMuMotherGenPt[2]->Fill(mcpt, motherpt);
	395	fHMuMotherRecPt[2]->Fill(trackpt, motherpt);
	396	fHMuMotherGenPhi[2]->Fill(mcphi, motherphi);
	397	fHMuMotherRecPhi[2]->Fill(trackphi, motherphi);
	398	fHMuMotherGenEta[2]->Fill(mceta, mothereta);
	399	fHMuMotherRecEta[2]->Fill(tracketa, mothereta);
	400	if(fESD) fHMuDCA[2]->Fill(dcavalue);
	401	}
	402	else if(motherpdg==211)
	403	{
	404	fHMuMotherGenPt[0]->Fill(mcpt, motherpt);
	405	fHMuMotherRecPt[0]->Fill(trackpt, motherpt);
	406	fHMuMotherGenPhi[0]->Fill(mcphi, motherphi);
	407	fHMuMotherRecPhi[0]->Fill(trackphi, motherphi);
	408	fHMuMotherGenEta[0]->Fill(mceta, mothereta);
	409	fHMuMotherRecEta[0]->Fill(tracketa, mothereta);
	410	if(fESD) fHMuDCA[0]->Fill(dcavalue);
	411	}
	412	else if(motherpdg==321)
	413	{
	414	fHMuMotherGenPt[1]->Fill(mcpt, motherpt);
	415	fHMuMotherRecPt[1]->Fill(trackpt, motherpt);
	416	fHMuMotherGenPhi[1]->Fill(mcphi, motherphi);
	417	fHMuMotherRecPhi[1]->Fill(trackphi, motherphi);
	418	fHMuMotherGenEta[1]->Fill(mceta, mothereta);
	419	fHMuMotherRecEta[1]->Fill(tracketa, mothereta);
	420	if(fESD) fHMuDCA[1]->Fill(dcavalue);
	421	}
	422	else
	423	{
	424	fHMuMotherGenPt[3]->Fill(mcpt, motherpt);
	425	fHMuMotherRecPt[3]->Fill(trackpt, motherpt);
	426	fHMuMotherGenPhi[3]->Fill(mcphi, motherphi);
	427	fHMuMotherRecPhi[3]->Fill(trackphi, motherphi);
	428	fHMuMotherGenEta[3]->Fill(mceta, mothereta);
	429	fHMuMotherRecEta[3]->Fill(tracketa, mothereta);
	430	if(fESD) fHMuDCA[3]->Fill(dcavalue);
	431	}
	432	}
	433	} // (mother hadron) : (daughter muon) QA */
	434	}
	435
7e990b20	436	PostData(1, fOutputList);
	437	return;
	438	}
	439
	440	//________________________________________________________________________
	441	void AliMuonEffMC::Terminate(Option_t *)
	442	{
	443	// Draw result to the screen
	444	// Called once at the end of the query
	445
7e990b20	446	}
	447
	448	//________________________________________________________________________
	449	Bool_t AliMuonEffMC::VertexOk(TObject* obj) const
	450	{
	451	// Modified from AliAnalyseLeadingTrackUE::VertexSelection()
	452
	453	Int_t nContributors = 0;
	454	Double_t zVertex = 999;
	455	TString name("");
	456
	457	if (obj->InheritsFrom("AliESDEvent")) {
	458	AliESDEvent* esdevt = (AliESDEvent*) obj;
	459	const AliESDVertex* vtx = esdevt->GetPrimaryVertex();
	460	if (!vtx)
	461	return 0;
	462	nContributors = vtx->GetNContributors();
	463	zVertex = vtx->GetZ();
	464	name = vtx->GetName();
	465	}
	466	else if (obj->InheritsFrom("AliAODEvent")) {
	467	AliAODEvent* aodevt = (AliAODEvent*) obj;
	468	if (aodevt->GetNumberOfVertices() < 1)
	469	return 0;
	470	const AliAODVertex* vtx = aodevt->GetPrimaryVertex();
	471	nContributors = vtx->GetNContributors();
	472	zVertex = vtx->GetZ();
	473	name = vtx->GetName();
	474	}
	475
	476	// Reject if TPC-only vertex
	477	if (name.CompareTo("TPCVertex")==0)
	478	return kFALSE;
	479
	480	// Check # contributors and range...
	481	if( nContributors < 1 \|\| TMath::Abs(zVertex) > 10 ) {
	482	return kFALSE;
	483	}
	484
	485	return kTRUE;
	486	}
	487
	488
	489	//________________________________________________________________________
	490	Bool_t AliMuonEffMC::IsGoodMUONtrack(AliESDMuonTrack &track)
	491	{
	492	// Applying track cuts for MUON tracks
	493	if(!track.ContainTrackerData()) return kFALSE;
7e990b20	494
	495	Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
	496	Double_t eta = track.Eta();
	497
	498	// Theta cut at absorber end
	499	if(thetaTrackAbsEnd <= 2. \|\| thetaTrackAbsEnd >= 10.) return kFALSE;
	500	// Eta cut
	501	if(eta <= -4. \|\| eta >= -2.5) return kFALSE;
	502	if(track.GetMatchTrigger() <= 0) return kFALSE;
	503
	504	return kTRUE;
	505	}
b673a083	506
	507	//________________________________________________________________________
	508	Bool_t AliMuonEffMC::IsGoodMUONtrack(AliAODTrack &track)
	509	{
	510	if (!track.IsMuonTrack()) return kFALSE;
	511
	512	Double_t dThetaAbs = TMath::ATan(track.GetRAtAbsorberEnd()/505.)
	513	* TMath::RadToDeg();
	514	if ((dThetaAbs<2.) \|\| (dThetaAbs>10.)) return kFALSE;
	515
	516	Double_t dEta = track.Eta();
	517	if ((dEta<-4.) \|\| (dEta>2.5)) return kFALSE;
	518
	519	if (track.GetMatchTrigger()<0.5) return kFALSE;
	520
	521	return kTRUE;
	522	}