[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / MuonHadron / AliAnalysisTaskMuonHadronCorrelations.cxx

#include "AliLog.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TTree.h"
#include "TAxis.h"
#include "AliAODEvent.h"
#include "AliAODTrack.h"
#include "TMath.h"
#include "TString.h"

#include "AliAnalysisTaskMuonHadronCorrelations.h"
#include "AliAnalysisManager.h"
#include "AliInputEventHandler.h"
#include "AliEventPoolManager.h"

ClassImp(AliAnalysisTaskMuonHadronCorrelations)

//====================================================================================================================================================

AliAnalysisTaskMuonHadronCorrelations::AliAnalysisTaskMuonHadronCorrelations() : 
  AliAnalysisTaskSE(), 
  fAOD(0x0),
  fPoolMgr(0x0),
  fTrackCB(0x0),
  fTrackMA(0x0),
  fFilterBitCentralBarrel(0),
  fMaxEtaCentralBarrel(1.0),
  fMaxChi2Muon(9999999999.), 
  fMinRAbsMuon(0), 
  fMaxRAbsMuon(9999999999.),
  fTriggerMatchLevelMuon(0),
  fNbinsCent(1), 
  fNbinsPt(1),
  fCentAxis(0x0), 
  fPtAxis(0x0),
  fHistV0Multiplicity(0x0), 
  fHistITSMultiplicity(0x0),
  fHistCentrality(0x0),
  fCentMethod(0),
  fOutputList(0x0)
{

  // Default constructor

  for (Int_t iCent=0; iCent<fNMaxBinsCentrality; iCent++) {
    for (Int_t iPtBinCB=0; iPtBinCB<fNMaxBinsPt; iPtBinCB++) {
      for (Int_t iPtBinMA=0; iPtBinMA<fNMaxBinsPt; iPtBinMA++) {
	fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA]    = NULL;
	fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA] = NULL;
      }
    }
    fHistNTracksCB_vs_NTracksMA[iCent]  = NULL;
    fHistNTracksCB_vs_NTracksMAmixed[iCent]  = NULL;
    fHistSingleMuonsEtaVsPt[iCent]   = NULL;
    fHistSingleMuonsEtaVsRAbs[iCent] = NULL;
  }  
  
}


//====================================================================================================================================================

AliAnalysisTaskMuonHadronCorrelations::AliAnalysisTaskMuonHadronCorrelations(const char *name) : 
  AliAnalysisTaskSE(name), 
  fAOD(0x0),
  fPoolMgr(0x0),
  fTrackCB(0x0),
  fTrackMA(0x0),
  fFilterBitCentralBarrel(0),
  fMaxEtaCentralBarrel(1.0),
  fMaxChi2Muon(9999999999.), 
  fMinRAbsMuon(0), 
  fMaxRAbsMuon(9999999999.),
  fTriggerMatchLevelMuon(0),
  fNbinsCent(1), 
  fNbinsPt(1),
  fCentAxis(0x0), 
  fPtAxis(0x0),
  fHistV0Multiplicity(0x0), 
  fHistITSMultiplicity(0x0),
  fHistCentrality(0x0),
  fCentMethod(0),
  fOutputList(0x0)
{

  // Constructor

  for (Int_t iCent=0; iCent<fNMaxBinsCentrality; iCent++) {
    for (Int_t iPtBinCB=0; iPtBinCB<fNMaxBinsPt; iPtBinCB++) {
      for (Int_t iPtBinMA=0; iPtBinMA<fNMaxBinsPt; iPtBinMA++) {
	fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA]    = NULL;
	fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA] = NULL;
      }
    }
    fHistNTracksCB_vs_NTracksMA[iCent]  = NULL;
    fHistNTracksCB_vs_NTracksMAmixed[iCent]  = NULL;
    fHistSingleMuonsEtaVsPt[iCent]   = NULL;
    fHistSingleMuonsEtaVsRAbs[iCent] = NULL;
  }  
  
  // Define input and output slots here
  DefineOutput(1, TList::Class());
  
}

//====================================================================================================================================================

AliAnalysisTaskMuonHadronCorrelations::~AliAnalysisTaskMuonHadronCorrelations() {
  
  delete fCentAxis;
  delete fPtAxis;

  if (fOutputList  && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) 
    delete fOutputList;
}

//====================================================================================================================================================

void AliAnalysisTaskMuonHadronCorrelations::UserCreateOutputObjects() {
  
  fOutputList = new TList();
  fOutputList->SetOwner(kTRUE);

  for (Int_t iCent=0; iCent<fNbinsCent; iCent++) {
    for (Int_t iPtBinCB=0; iPtBinCB<fNbinsPt; iPtBinCB++) {
      for (Int_t iPtBinMA=0; iPtBinMA<fNbinsPt; iPtBinMA++) {
	
	fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA]    = new TH1D(Form("fHistDeltaPhi_Cent%02d_PtBin%02d_%02d",iCent,iPtBinCB,iPtBinMA), 
							       Form("%d-%d %%, %3.1f<p_{T}^{TPC}<%3.1f, %3.1f<p_{T}^{Muon}<%3.1f",
								    Int_t(fCentAxis->GetBinLowEdge(iCent+1)),
								    Int_t(fCentAxis->GetBinUpEdge(iCent+1)),
								    fPtAxis->GetBinLowEdge(iPtBinCB+1),
								    fPtAxis->GetBinUpEdge(iPtBinCB+1),
								    fPtAxis->GetBinLowEdge(iPtBinMA+1),
								    fPtAxis->GetBinUpEdge(iPtBinMA+1)),
							       100, -0.5*TMath::RadToDeg()*TMath::Pi(), 1.5*TMath::RadToDeg()*TMath::Pi());
	
	fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA]    = new TH1D(Form("fHistDeltaPhiMix_Cent%02d_PtBin%02d_%02d",iCent,iPtBinCB,iPtBinMA), 
								  Form("%d-%d %%, %3.1f<p_{T}^{TPC}<%3.1f, %3.1f<p_{T}^{Muon}<%3.1f MIXED",
								       Int_t(fCentAxis->GetBinLowEdge(iCent+1)),
								       Int_t(fCentAxis->GetBinUpEdge(iCent+1)),
								       fPtAxis->GetBinLowEdge(iPtBinCB+1),
								       fPtAxis->GetBinUpEdge(iPtBinCB+1),
								       fPtAxis->GetBinLowEdge(iPtBinMA+1),
								       fPtAxis->GetBinUpEdge(iPtBinMA+1)),
								  100, -0.5*TMath::RadToDeg()*TMath::Pi(), 1.5*TMath::RadToDeg()*TMath::Pi());
	
	fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA]    -> SetXTitle("#Delta#varphi  [degrees]");
	fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA] -> SetXTitle("#Delta#varphi  [degrees]");

	fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA]    -> Sumw2();
	fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA] -> Sumw2();
	
	fOutputList -> Add(fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA]);
	fOutputList -> Add(fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA]);

      }
    }

    fHistNTracksCB_vs_NTracksMA[iCent] = new TH2D(Form("fHistNTracksCB_vs_NTracksMA_Cent%02d",iCent),
						  Form("%d-%d %%",Int_t(fCentAxis->GetBinLowEdge(iCent+1)),Int_t(fCentAxis->GetBinUpEdge(iCent+1))),
						  100, 0, 500, 20, 0, 20);
    fHistNTracksCB_vs_NTracksMA[iCent] -> SetXTitle("N_{tracks} Central Barrel");
    fHistNTracksCB_vs_NTracksMA[iCent] -> SetYTitle("N_{tracks} Muon Arm");
    fHistNTracksCB_vs_NTracksMA[iCent] -> Sumw2();

    fHistNTracksCB_vs_NTracksMAmixed[iCent] = new TH2D(Form("fHistNTracksCB_vs_NTracksMAmixed_Cent%02d",iCent),
						  Form("%d-%d %% MIXED",Int_t(fCentAxis->GetBinLowEdge(iCent+1)),Int_t(fCentAxis->GetBinUpEdge(iCent+1))),
						  100, 0, 500, 20, 0, 20);
    fHistNTracksCB_vs_NTracksMAmixed[iCent] -> SetXTitle("N_{tracks} Central Barrel");
    fHistNTracksCB_vs_NTracksMAmixed[iCent] -> SetYTitle("N_{tracks} Muon Arm");
    fHistNTracksCB_vs_NTracksMAmixed[iCent] -> Sumw2();

    fOutputList -> Add(fHistNTracksCB_vs_NTracksMA[iCent]);
    fOutputList -> Add(fHistNTracksCB_vs_NTracksMAmixed[iCent]);

    fHistSingleMuonsEtaVsPt[iCent]   = new TH2D(Form("fHistSingleMuonsEtaVsPt_Cent%02d",iCent),
					       "#eta vs p_{T} for single muons",
					       100, -4.5, -2., 100, 0., 10.);
    fOutputList->Add(fHistSingleMuonsEtaVsPt[iCent]);
    fHistSingleMuonsEtaVsRAbs[iCent] = new TH2D(Form("fHistSingleMuonsEtaVsRAbs_Cent%02d",iCent),
					       "#eta vs R_{Abs} for single muons",
					       100, -4.5, -2., 100, 0, 100);
    fOutputList->Add(fHistSingleMuonsEtaVsRAbs[iCent]);
    
  }
  
  fHistV0Multiplicity = new TH1D("fHistV0Multiplicity", "V0 Multiplicity", 500, 0, 1000);
  fHistV0Multiplicity -> SetXTitle("Multiplicity");
  fHistV0Multiplicity -> Sumw2();

  fHistITSMultiplicity = new TH1D("fHistITSMultiplicity", "ITS Multiplicity", 500, 0, 500);
  fHistITSMultiplicity -> SetXTitle("N_{Clusters}");
  fHistITSMultiplicity -> Sumw2();

  fHistCentrality = new TH1D("fHistCentrality", Form("%s Centrality",fCentMethod.Data()), 300, -100, 200);
  fHistCentrality -> SetXTitle("Centrality  [%]");
  fHistCentrality -> Sumw2();

  fOutputList -> Add(fHistV0Multiplicity);
  fOutputList -> Add(fHistITSMultiplicity);
  fOutputList -> Add(fHistCentrality);

  const Int_t kNZvtxBins  = 10;
  // bins for further buffers are shifted by 100 cm
  Double_t vertexBins[kNZvtxBins+1] = { -10,   -8,  -6,  -4,  -2,   0,   2,   4,   6,   8,  10 };
  Int_t nZvtxBins  = kNZvtxBins;
  Double_t* zvtxbin = vertexBins;

  fPoolMgr = new AliEventPoolManager(1000, 20000, fNbinsCent, (Double_t*)fCentAxis->GetXbins()->GetArray(), nZvtxBins, zvtxbin);

  PostData(1, fOutputList); 

}

//====================================================================================================================================================

void AliAnalysisTaskMuonHadronCorrelations::UserExec(Option_t *) {

  fAOD = dynamic_cast<AliAODEvent *>(InputEvent());  
  if (!fAOD) return;  
  
  // Trigger selection
  if (!(IsTriggerFired())) return;
  
  fHistV0Multiplicity  -> Fill(GetV0Multiplicity());
  fHistITSMultiplicity -> Fill(GetITSMultiplicity());

  Int_t centBin = GetCentBin();
  if (centBin<0) return;
  Double_t percentile = fAOD->GetCentrality()->GetCentralityPercentile(fCentMethod.Data());
  fHistCentrality->Fill(percentile);

  const AliVVertex* vertex = fAOD->GetPrimaryVertex();
  Double_t zVtx = vertex->GetZ();
  if (TMath::Abs(zVtx) > 10.) return;
  
  TObjArray *tracksMuonArm = GetAcceptedTracksMuonArm(fAOD,centBin);
  if (tracksMuonArm->GetEntriesFast() == 0) {
    delete tracksMuonArm;
    return;
  }
  TObjArray *tracksCentralBarrel = GetAcceptedTracksCentralBarrel(fAOD);
  
  fHistNTracksCB_vs_NTracksMA[centBin] -> Fill(tracksCentralBarrel->GetEntries(), tracksMuonArm->GetEntries());

  AliDebug(1, Form("Single Event analysis : nTracksCB = %4d, nTracksMA = %4d", tracksCentralBarrel->GetEntries(), tracksMuonArm->GetEntries()));

  // Same event  
  for (Int_t iTrMA=0; iTrMA<tracksMuonArm->GetEntriesFast(); iTrMA++) {
    fTrackMA = (AliAODTrack*) tracksMuonArm->At(iTrMA);
    for (Int_t iTrCB=0; iTrCB<tracksCentralBarrel->GetEntriesFast(); iTrCB++) {
      fTrackCB = (AliAODTrack*) tracksCentralBarrel -> At(iTrCB);
      FillHistograms(centBin, kSingleEvent);
    }
  }

  // Mixed event
  {
    AliEventPool* pool = fPoolMgr->GetEventPool(percentile, zVtx);
    //pool->PrintInfo();
    if (pool->IsReady() || pool->NTracksInPool() > 2000 || pool->GetCurrentNEvents() >= 5) 
      for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) {
	TObjArray *mixedTracks = pool->GetEvent(jMix);
	fHistNTracksCB_vs_NTracksMAmixed[centBin]->Fill(mixedTracks->GetEntriesFast(), tracksMuonArm->GetEntriesFast());
	for (Int_t iTrMA=0; iTrMA<tracksMuonArm->GetEntriesFast(); iTrMA++) {
	  fTrackMA = (AliAODTrack*) tracksMuonArm->At(iTrMA);
	  for (Int_t iTrCB=0; iTrCB<mixedTracks->GetEntriesFast(); iTrCB++) {
	    fTrackCB = (AliAODTrack*) mixedTracks -> At(iTrCB);
	    FillHistograms(centBin, kMixedEvent);
	  }
	}
      }
    pool->UpdatePool(tracksCentralBarrel);
  }

  delete tracksMuonArm;

  PostData(1, fOutputList); 

}

//====================================================================================================================================================

void AliAnalysisTaskMuonHadronCorrelations::FillHistograms(Int_t centrality, Int_t option) {

  Int_t ptBinTrackCB = fPtAxis -> FindBin(fTrackCB->Pt());
  Int_t ptBinTrackMA = fPtAxis -> FindBin(fTrackMA->Pt());

  if (ptBinTrackCB<1 || ptBinTrackCB>fNbinsPt || ptBinTrackMA<1 || ptBinTrackMA>fNbinsPt) return;

  Double_t deltaPhi = fTrackCB->Phi() - fTrackMA->Phi();
  if (deltaPhi >  1.5*TMath::Pi()) deltaPhi -= TMath::TwoPi();
  if (deltaPhi < -0.5*TMath::Pi()) deltaPhi += TMath::TwoPi();

  if (option==kSingleEvent)     fHistDeltaPhi[centrality][ptBinTrackCB-1][ptBinTrackMA-1]    -> Fill(TMath::RadToDeg()*deltaPhi);
  else if (option==kMixedEvent) fHistDeltaPhiMix[centrality][ptBinTrackCB-1][ptBinTrackMA-1] -> Fill(TMath::RadToDeg()*deltaPhi);

}

//====================================================================================================================================================

Bool_t AliAnalysisTaskMuonHadronCorrelations::IsTriggerFired() {
  
  Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7); 

  return isSelected;
}

//====================================================================================================================================================

Float_t AliAnalysisTaskMuonHadronCorrelations::GetV0Multiplicity() {
  
  Float_t multiplicity=0;
  for (Int_t iChannel=0; iChannel<64; iChannel++) multiplicity += fAOD->GetVZEROData()->GetMultiplicity(iChannel);
  return multiplicity;

}

//====================================================================================================================================================

TObjArray* AliAnalysisTaskMuonHadronCorrelations::GetAcceptedTracksCentralBarrel(AliAODEvent *aodEvent) {

  // fills the array of central barrel tracks that pass the cuts

  TObjArray *tracks = new TObjArray;
  tracks->SetOwner(kTRUE);

  Int_t nTracks = aodEvent->GetNTracks();

  AliAODTrack *track = 0;

  for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
    track = aodEvent->GetTrack(iTrack);
    if (track->TestFilterBit(fFilterBitCentralBarrel) && TMath::Abs(track->Eta())<fMaxEtaCentralBarrel) {
      tracks->Add(new AliAODTrack(*track));
    }
  }

  return tracks;

}

//====================================================================================================================================================

TObjArray* AliAnalysisTaskMuonHadronCorrelations::GetAcceptedTracksMuonArm(AliAODEvent *aodEvent, Int_t centBin) {

  // fills the array of muon tracks that pass the cuts

  TObjArray *tracks = new TObjArray;
  tracks->SetOwner(kFALSE);

  Int_t nTracks = aodEvent->GetNTracks();

  AliAODTrack *track = 0;
  
  for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
    track = aodEvent->GetTrack(iTrack);
    if (track->IsMuonTrack() && track->GetMatchTrigger()>=fTriggerMatchLevelMuon) {
      if (track->Chi2perNDF() < 5.) {
	// histos
	fHistSingleMuonsEtaVsPt[centBin]->Fill(track->Eta(),track->Pt());
	fHistSingleMuonsEtaVsRAbs[centBin]->Fill(track->Eta(),track->GetRAtAbsorberEnd());
	if (track->Eta() > -4. && track->Eta() < -2.5) {
	  Double_t rabs = track->GetRAtAbsorberEnd();
	  if (rabs > 17.6 && rabs < 89.5) {
	    tracks->Add(new AliAODTrack(*track));
	  }
	}
      }
    }
  }

  return tracks;

}

//====================================================================================================================================================

void AliAnalysisTaskMuonHadronCorrelations::SetCentBinning(Int_t nBins, Double_t *limits) {

  if (nBins>fNMaxBinsCentrality) {
    AliInfo(Form("WARNING : only %d centrality bins (out of the %d proposed) will be considered",fNMaxBinsCentrality,nBins));
    nBins = fNMaxBinsCentrality;
  }
  if (nBins<=0) {
    AliInfo("WARNING : at least one centrality bin must be considered");
    nBins = 1;
  }
  
  fNbinsCent = nBins;
  fCentAxis  = new TAxis(fNbinsCent, limits);

}

//====================================================================================================================================================

void AliAnalysisTaskMuonHadronCorrelations::SetPtBinning(Int_t nBins, Double_t *limits) {

  if (nBins>fNMaxBinsPt) {
    AliInfo(Form("WARNING : only %d pt bins (out of the %d proposed) will be considered",fNMaxBinsPt,nBins));
    nBins = fNMaxBinsPt;
  }
  if (nBins<=0) {
    AliInfo("WARNING : at least one pt bin must be considered");
    nBins = 1;
  }
  
  fNbinsPt = nBins;
  fPtAxis  = new TAxis(fNbinsPt, limits);

}

//====================================================================================================================================================

Int_t AliAnalysisTaskMuonHadronCorrelations::GetCentBin() {

  Double_t percentile = fAOD->GetCentrality()->GetCentralityPercentile(fCentMethod.Data());

  Int_t bin = fCentAxis->FindBin(percentile) - 1;
  if (bin >= fNbinsCent) bin = -1;
  return bin;
  
}

//====================================================================================================================================================

Double_t AliAnalysisTaskMuonHadronCorrelations::GetITSMultiplicity() {

  Double_t multiplicity = fAOD->GetHeader()->GetNumberOfITSClusters(1);

  return multiplicity;

}

//====================================================================================================================================================

void AliAnalysisTaskMuonHadronCorrelations::Terminate(Option_t *) {

  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
  if (!fOutputList) {
    printf("ERROR: Output list not available\n");
    return;
  }  

}

//====================================================================================================================================================
Commit	Line	Data
07c2615a	1	#include "AliLog.h"
	2	#include "TH1D.h"
	3	#include "TH2D.h"
	4	#include "TTree.h"
	5	#include "TAxis.h"
	6	#include "AliAODEvent.h"
	7	#include "AliAODTrack.h"
	8	#include "TMath.h"
	9	#include "TString.h"
	10
	11	#include "AliAnalysisTaskMuonHadronCorrelations.h"
6de71f1d	12	#include "AliAnalysisManager.h"
	13	#include "AliInputEventHandler.h"
	14	#include "AliEventPoolManager.h"
07c2615a	15
	16	ClassImp(AliAnalysisTaskMuonHadronCorrelations)
	17
	18	//====================================================================================================================================================
	19
	20	AliAnalysisTaskMuonHadronCorrelations::AliAnalysisTaskMuonHadronCorrelations() :
	21	AliAnalysisTaskSE(),
	22	fAOD(0x0),
6de71f1d	23	fPoolMgr(0x0),
07c2615a	24	fTrackCB(0x0),
	25	fTrackMA(0x0),
	26	fFilterBitCentralBarrel(0),
	27	fMaxEtaCentralBarrel(1.0),
29d66ae8	28	fMaxChi2Muon(9999999999.),
07c2615a	29	fMinRAbsMuon(0),
29d66ae8	30	fMaxRAbsMuon(9999999999.),
07c2615a	31	fTriggerMatchLevelMuon(0),
07c2615a	32	fNbinsCent(1),
	33	fNbinsPt(1),
	34	fCentAxis(0x0),
07c2615a	35	fPtAxis(0x0),
	36	fHistV0Multiplicity(0x0),
	37	fHistITSMultiplicity(0x0),
	38	fHistCentrality(0x0),
	39	fCentMethod(0),
07c2615a	40	fOutputList(0x0)
	41	{
	42
	43	// Default constructor
	44
	45	for (Int_t iCent=0; iCent<fNMaxBinsCentrality; iCent++) {
	46	for (Int_t iPtBinCB=0; iPtBinCB<fNMaxBinsPt; iPtBinCB++) {
	47	for (Int_t iPtBinMA=0; iPtBinMA<fNMaxBinsPt; iPtBinMA++) {
	48	fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA] = NULL;
	49	fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA] = NULL;
	50	}
	51	}
	52	fHistNTracksCB_vs_NTracksMA[iCent] = NULL;
57a9cc15	53	fHistNTracksCB_vs_NTracksMAmixed[iCent] = NULL;
	54	fHistSingleMuonsEtaVsPt[iCent] = NULL;
	55	fHistSingleMuonsEtaVsRAbs[iCent] = NULL;
07c2615a	56	}
	57
	58	}
	59
	60
	61	//====================================================================================================================================================
	62
	63	AliAnalysisTaskMuonHadronCorrelations::AliAnalysisTaskMuonHadronCorrelations(const char *name) :
	64	AliAnalysisTaskSE(name),
	65	fAOD(0x0),
6de71f1d	66	fPoolMgr(0x0),
07c2615a	67	fTrackCB(0x0),
	68	fTrackMA(0x0),
	69	fFilterBitCentralBarrel(0),
	70	fMaxEtaCentralBarrel(1.0),
29d66ae8	71	fMaxChi2Muon(9999999999.),
07c2615a	72	fMinRAbsMuon(0),
29d66ae8	73	fMaxRAbsMuon(9999999999.),
07c2615a	74	fTriggerMatchLevelMuon(0),
07c2615a	75	fNbinsCent(1),
	76	fNbinsPt(1),
	77	fCentAxis(0x0),
07c2615a	78	fPtAxis(0x0),
	79	fHistV0Multiplicity(0x0),
	80	fHistITSMultiplicity(0x0),
	81	fHistCentrality(0x0),
	82	fCentMethod(0),
07c2615a	83	fOutputList(0x0)
	84	{
	85
	86	// Constructor
	87
	88	for (Int_t iCent=0; iCent<fNMaxBinsCentrality; iCent++) {
	89	for (Int_t iPtBinCB=0; iPtBinCB<fNMaxBinsPt; iPtBinCB++) {
	90	for (Int_t iPtBinMA=0; iPtBinMA<fNMaxBinsPt; iPtBinMA++) {
	91	fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA] = NULL;
	92	fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA] = NULL;
	93	}
	94	}
	95	fHistNTracksCB_vs_NTracksMA[iCent] = NULL;
57a9cc15	96	fHistNTracksCB_vs_NTracksMAmixed[iCent] = NULL;
	97	fHistSingleMuonsEtaVsPt[iCent] = NULL;
	98	fHistSingleMuonsEtaVsRAbs[iCent] = NULL;
07c2615a	99	}
	100
	101	// Define input and output slots here
	102	DefineOutput(1, TList::Class());
	103
	104	}
	105
	106	//====================================================================================================================================================
	107
	108	AliAnalysisTaskMuonHadronCorrelations::~AliAnalysisTaskMuonHadronCorrelations() {
	109
07c2615a	110	delete fCentAxis;
	111	delete fPtAxis;
	112
6de71f1d	113	if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())
6de71f1d	114	delete fOutputList;
07c2615a	115	}
	116
	117	//====================================================================================================================================================
	118
	119	void AliAnalysisTaskMuonHadronCorrelations::UserCreateOutputObjects() {
	120
	121	fOutputList = new TList();
	122	fOutputList->SetOwner(kTRUE);
	123
	124	for (Int_t iCent=0; iCent<fNbinsCent; iCent++) {
	125	for (Int_t iPtBinCB=0; iPtBinCB<fNbinsPt; iPtBinCB++) {
	126	for (Int_t iPtBinMA=0; iPtBinMA<fNbinsPt; iPtBinMA++) {
	127
	128	fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA] = new TH1D(Form("fHistDeltaPhi_Cent%02d_PtBin%02d_%02d",iCent,iPtBinCB,iPtBinMA),
	129	Form("%d-%d %%, %3.1f<p_{T}^{TPC}<%3.1f, %3.1f<p_{T}^{Muon}<%3.1f",
	130	Int_t(fCentAxis->GetBinLowEdge(iCent+1)),
	131	Int_t(fCentAxis->GetBinUpEdge(iCent+1)),
	132	fPtAxis->GetBinLowEdge(iPtBinCB+1),
	133	fPtAxis->GetBinUpEdge(iPtBinCB+1),
	134	fPtAxis->GetBinLowEdge(iPtBinMA+1),
	135	fPtAxis->GetBinUpEdge(iPtBinMA+1)),
	136	100, -0.5TMath::RadToDeg()TMath::Pi(), 1.5TMath::RadToDeg()TMath::Pi());
	137
	138	fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA] = new TH1D(Form("fHistDeltaPhiMix_Cent%02d_PtBin%02d_%02d",iCent,iPtBinCB,iPtBinMA),
	139	Form("%d-%d %%, %3.1f<p_{T}^{TPC}<%3.1f, %3.1f<p_{T}^{Muon}<%3.1f MIXED",
	140	Int_t(fCentAxis->GetBinLowEdge(iCent+1)),
	141	Int_t(fCentAxis->GetBinUpEdge(iCent+1)),
	142	fPtAxis->GetBinLowEdge(iPtBinCB+1),
	143	fPtAxis->GetBinUpEdge(iPtBinCB+1),
	144	fPtAxis->GetBinLowEdge(iPtBinMA+1),
	145	fPtAxis->GetBinUpEdge(iPtBinMA+1)),
	146	100, -0.5TMath::RadToDeg()TMath::Pi(), 1.5TMath::RadToDeg()TMath::Pi());
	147
	148	fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA] -> SetXTitle("#Delta#varphi [degrees]");
	149	fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA] -> SetXTitle("#Delta#varphi [degrees]");
	150
	151	fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA] -> Sumw2();
	152	fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA] -> Sumw2();
	153
	154	fOutputList -> Add(fHistDeltaPhi[iCent][iPtBinCB][iPtBinMA]);
	155	fOutputList -> Add(fHistDeltaPhiMix[iCent][iPtBinCB][iPtBinMA]);
	156
	157	}
	158	}
	159
	160	fHistNTracksCB_vs_NTracksMA[iCent] = new TH2D(Form("fHistNTracksCB_vs_NTracksMA_Cent%02d",iCent),
	161	Form("%d-%d %%",Int_t(fCentAxis->GetBinLowEdge(iCent+1)),Int_t(fCentAxis->GetBinUpEdge(iCent+1))),
57a9cc15	162	100, 0, 500, 20, 0, 20);
07c2615a	163	fHistNTracksCB_vs_NTracksMA[iCent] -> SetXTitle("N_{tracks} Central Barrel");
07c2615a	164	fHistNTracksCB_vs_NTracksMA[iCent] -> SetYTitle("N_{tracks} Muon Arm");
07c2615a	165	fHistNTracksCB_vs_NTracksMA[iCent] -> Sumw2();
07c2615a	166
57a9cc15	167	fHistNTracksCB_vs_NTracksMAmixed[iCent] = new TH2D(Form("fHistNTracksCB_vs_NTracksMAmixed_Cent%02d",iCent),
	168	Form("%d-%d %% MIXED",Int_t(fCentAxis->GetBinLowEdge(iCent+1)),Int_t(fCentAxis->GetBinUpEdge(iCent+1))),
	169	100, 0, 500, 20, 0, 20);
	170	fHistNTracksCB_vs_NTracksMAmixed[iCent] -> SetXTitle("N_{tracks} Central Barrel");
	171	fHistNTracksCB_vs_NTracksMAmixed[iCent] -> SetYTitle("N_{tracks} Muon Arm");
	172	fHistNTracksCB_vs_NTracksMAmixed[iCent] -> Sumw2();
	173
07c2615a	174	fOutputList -> Add(fHistNTracksCB_vs_NTracksMA[iCent]);
57a9cc15	175	fOutputList -> Add(fHistNTracksCB_vs_NTracksMAmixed[iCent]);
	176
	177	fHistSingleMuonsEtaVsPt[iCent] = new TH2D(Form("fHistSingleMuonsEtaVsPt_Cent%02d",iCent),
	178	"#eta vs p_{T} for single muons",
	179	100, -4.5, -2., 100, 0., 10.);
	180	fOutputList->Add(fHistSingleMuonsEtaVsPt[iCent]);
	181	fHistSingleMuonsEtaVsRAbs[iCent] = new TH2D(Form("fHistSingleMuonsEtaVsRAbs_Cent%02d",iCent),
	182	"#eta vs R_{Abs} for single muons",
	183	100, -4.5, -2., 100, 0, 100);
	184	fOutputList->Add(fHistSingleMuonsEtaVsRAbs[iCent]);
07c2615a	185
	186	}
	187
	188	fHistV0Multiplicity = new TH1D("fHistV0Multiplicity", "V0 Multiplicity", 500, 0, 1000);
	189	fHistV0Multiplicity -> SetXTitle("Multiplicity");
	190	fHistV0Multiplicity -> Sumw2();
	191
	192	fHistITSMultiplicity = new TH1D("fHistITSMultiplicity", "ITS Multiplicity", 500, 0, 500);
	193	fHistITSMultiplicity -> SetXTitle("N_{Clusters}");
	194	fHistITSMultiplicity -> Sumw2();
	195
	196	fHistCentrality = new TH1D("fHistCentrality", Form("%s Centrality",fCentMethod.Data()), 300, -100, 200);
	197	fHistCentrality -> SetXTitle("Centrality [%]");
	198	fHistCentrality -> Sumw2();
	199
	200	fOutputList -> Add(fHistV0Multiplicity);
	201	fOutputList -> Add(fHistITSMultiplicity);
	202	fOutputList -> Add(fHistCentrality);
	203
6de71f1d	204	const Int_t kNZvtxBins = 10;
	205	// bins for further buffers are shifted by 100 cm
	206	Double_t vertexBins[kNZvtxBins+1] = { -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10 };
	207	Int_t nZvtxBins = kNZvtxBins;
	208	Double_t* zvtxbin = vertexBins;
	209
	210	fPoolMgr = new AliEventPoolManager(1000, 20000, fNbinsCent, (Double_t*)fCentAxis->GetXbins()->GetArray(), nZvtxBins, zvtxbin);
	211
07c2615a	212	PostData(1, fOutputList);
	213
	214	}
	215
	216	//====================================================================================================================================================
	217
	218	void AliAnalysisTaskMuonHadronCorrelations::UserExec(Option_t *) {
	219
07c2615a	220	fAOD = dynamic_cast<AliAODEvent *>(InputEvent());
	221	if (!fAOD) return;
	222
	223	// Trigger selection
	224	if (!(IsTriggerFired())) return;
	225
	226	fHistV0Multiplicity -> Fill(GetV0Multiplicity());
	227	fHistITSMultiplicity -> Fill(GetITSMultiplicity());
6de71f1d	228
07c2615a	229	Int_t centBin = GetCentBin();
07c2615a	230	if (centBin<0) return;
6de71f1d	231	Double_t percentile = fAOD->GetCentrality()->GetCentralityPercentile(fCentMethod.Data());
	232	fHistCentrality->Fill(percentile);
	233
	234	const AliVVertex* vertex = fAOD->GetPrimaryVertex();
	235	Double_t zVtx = vertex->GetZ();
	236	if (TMath::Abs(zVtx) > 10.) return;
07c2615a	237
57a9cc15	238	TObjArray *tracksMuonArm = GetAcceptedTracksMuonArm(fAOD,centBin);
6de71f1d	239	if (tracksMuonArm->GetEntriesFast() == 0) {
	240	delete tracksMuonArm;
	241	return;
	242	}
	243	TObjArray *tracksCentralBarrel = GetAcceptedTracksCentralBarrel(fAOD);
07c2615a	244
6de71f1d	245	fHistNTracksCB_vs_NTracksMA[centBin] -> Fill(tracksCentralBarrel->GetEntries(), tracksMuonArm->GetEntries());
07c2615a	246
6de71f1d	247	AliDebug(1, Form("Single Event analysis : nTracksCB = %4d, nTracksMA = %4d", tracksCentralBarrel->GetEntries(), tracksMuonArm->GetEntries()));
	248
	249	// Same event
	250	for (Int_t iTrMA=0; iTrMA<tracksMuonArm->GetEntriesFast(); iTrMA++) {
	251	fTrackMA = (AliAODTrack*) tracksMuonArm->At(iTrMA);
	252	for (Int_t iTrCB=0; iTrCB<tracksCentralBarrel->GetEntriesFast(); iTrCB++) {
	253	fTrackCB = (AliAODTrack*) tracksCentralBarrel -> At(iTrCB);
07c2615a	254	FillHistograms(centBin, kSingleEvent);
	255	}
	256	}
	257
6de71f1d	258	// Mixed event
	259	{
	260	AliEventPool* pool = fPoolMgr->GetEventPool(percentile, zVtx);
	261	//pool->PrintInfo();
	262	if (pool->IsReady() \|\| pool->NTracksInPool() > 2000 \|\| pool->GetCurrentNEvents() >= 5)
	263	for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) {
57a9cc15	264	TObjArray *mixedTracks = pool->GetEvent(jMix);
57a9cc15	265	fHistNTracksCB_vs_NTracksMAmixed[centBin]->Fill(mixedTracks->GetEntriesFast(), tracksMuonArm->GetEntriesFast());
6de71f1d	266	for (Int_t iTrMA=0; iTrMA<tracksMuonArm->GetEntriesFast(); iTrMA++) {
6de71f1d	267	fTrackMA = (AliAODTrack*) tracksMuonArm->At(iTrMA);
6de71f1d	268	for (Int_t iTrCB=0; iTrCB<mixedTracks->GetEntriesFast(); iTrCB++) {
	269	fTrackCB = (AliAODTrack*) mixedTracks -> At(iTrCB);
	270	FillHistograms(centBin, kMixedEvent);
	271	}
	272	}
	273	}
	274	pool->UpdatePool(tracksCentralBarrel);
	275	}
	276
	277	delete tracksMuonArm;
07c2615a	278
	279	PostData(1, fOutputList);
	280
	281	}
	282
	283	//====================================================================================================================================================
	284
	285	void AliAnalysisTaskMuonHadronCorrelations::FillHistograms(Int_t centrality, Int_t option) {
	286
	287	Int_t ptBinTrackCB = fPtAxis -> FindBin(fTrackCB->Pt());
	288	Int_t ptBinTrackMA = fPtAxis -> FindBin(fTrackMA->Pt());
	289
	290	if (ptBinTrackCB<1 \|\| ptBinTrackCB>fNbinsPt \|\| ptBinTrackMA<1 \|\| ptBinTrackMA>fNbinsPt) return;
	291
	292	Double_t deltaPhi = fTrackCB->Phi() - fTrackMA->Phi();
	293	if (deltaPhi > 1.5*TMath::Pi()) deltaPhi -= TMath::TwoPi();
	294	if (deltaPhi < -0.5*TMath::Pi()) deltaPhi += TMath::TwoPi();
	295
	296	if (option==kSingleEvent) fHistDeltaPhi[centrality][ptBinTrackCB-1][ptBinTrackMA-1] -> Fill(TMath::RadToDeg()*deltaPhi);
	297	else if (option==kMixedEvent) fHistDeltaPhiMix[centrality][ptBinTrackCB-1][ptBinTrackMA-1] -> Fill(TMath::RadToDeg()*deltaPhi);
	298
	299	}
	300
	301	//====================================================================================================================================================
	302
	303	Bool_t AliAnalysisTaskMuonHadronCorrelations::IsTriggerFired() {
	304
6de71f1d	305	Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7);
07c2615a	306
6de71f1d	307	return isSelected;
07c2615a	308	}
	309
	310	//====================================================================================================================================================
	311
	312	Float_t AliAnalysisTaskMuonHadronCorrelations::GetV0Multiplicity() {
	313
	314	Float_t multiplicity=0;
	315	for (Int_t iChannel=0; iChannel<64; iChannel++) multiplicity += fAOD->GetVZEROData()->GetMultiplicity(iChannel);
	316	return multiplicity;
	317
	318	}
	319
	320	//====================================================================================================================================================
	321
6de71f1d	322	TObjArray* AliAnalysisTaskMuonHadronCorrelations::GetAcceptedTracksCentralBarrel(AliAODEvent *aodEvent) {
07c2615a	323
	324	// fills the array of central barrel tracks that pass the cuts
	325
6de71f1d	326	TObjArray *tracks = new TObjArray;
6de71f1d	327	tracks->SetOwner(kTRUE);
07c2615a	328
	329	Int_t nTracks = aodEvent->GetNTracks();
	330
	331	AliAODTrack *track = 0;
	332
	333	for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
	334	track = aodEvent->GetTrack(iTrack);
	335	if (track->TestFilterBit(fFilterBitCentralBarrel) && TMath::Abs(track->Eta())<fMaxEtaCentralBarrel) {
6de71f1d	336	tracks->Add(new AliAODTrack(*track));
07c2615a	337	}
	338	}
	339
	340	return tracks;
	341
	342	}
	343
	344	//====================================================================================================================================================
	345
57a9cc15	346	TObjArray* AliAnalysisTaskMuonHadronCorrelations::GetAcceptedTracksMuonArm(AliAODEvent *aodEvent, Int_t centBin) {
07c2615a	347
	348	// fills the array of muon tracks that pass the cuts
	349
6de71f1d	350	TObjArray *tracks = new TObjArray;
6de71f1d	351	tracks->SetOwner(kFALSE);
07c2615a	352
	353	Int_t nTracks = aodEvent->GetNTracks();
	354
	355	AliAODTrack *track = 0;
	356
	357	for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
	358	track = aodEvent->GetTrack(iTrack);
	359	if (track->IsMuonTrack() && track->GetMatchTrigger()>=fTriggerMatchLevelMuon) {
57a9cc15	360	if (track->Chi2perNDF() < 5.) {
	361	// histos
	362	fHistSingleMuonsEtaVsPt[centBin]->Fill(track->Eta(),track->Pt());
	363	fHistSingleMuonsEtaVsRAbs[centBin]->Fill(track->Eta(),track->GetRAtAbsorberEnd());
	364	if (track->Eta() > -4. && track->Eta() < -2.5) {
	365	Double_t rabs = track->GetRAtAbsorberEnd();
	366	if (rabs > 17.6 && rabs < 89.5) {
	367	tracks->Add(new AliAODTrack(*track));
	368	}
	369	}
	370	}
07c2615a	371	}
	372	}
	373
	374	return tracks;
	375
	376	}
	377
	378	//====================================================================================================================================================
	379
	380	void AliAnalysisTaskMuonHadronCorrelations::SetCentBinning(Int_t nBins, Double_t *limits) {
	381
	382	if (nBins>fNMaxBinsCentrality) {
	383	AliInfo(Form("WARNING : only %d centrality bins (out of the %d proposed) will be considered",fNMaxBinsCentrality,nBins));
	384	nBins = fNMaxBinsCentrality;
	385	}
	386	if (nBins<=0) {
	387	AliInfo("WARNING : at least one centrality bin must be considered");
	388	nBins = 1;
	389	}
	390
	391	fNbinsCent = nBins;
	392	fCentAxis = new TAxis(fNbinsCent, limits);
	393
	394	}
	395
	396	//====================================================================================================================================================
	397
07c2615a	398	void AliAnalysisTaskMuonHadronCorrelations::SetPtBinning(Int_t nBins, Double_t *limits) {
	399
	400	if (nBins>fNMaxBinsPt) {
	401	AliInfo(Form("WARNING : only %d pt bins (out of the %d proposed) will be considered",fNMaxBinsPt,nBins));
	402	nBins = fNMaxBinsPt;
	403	}
	404	if (nBins<=0) {
	405	AliInfo("WARNING : at least one pt bin must be considered");
	406	nBins = 1;
	407	}
	408
	409	fNbinsPt = nBins;
	410	fPtAxis = new TAxis(fNbinsPt, limits);
	411
	412	}
	413
	414	//====================================================================================================================================================
	415
	416	Int_t AliAnalysisTaskMuonHadronCorrelations::GetCentBin() {
07c2615a	417
6de71f1d	418	Double_t percentile = fAOD->GetCentrality()->GetCentralityPercentile(fCentMethod.Data());
07c2615a	419
6de71f1d	420	Int_t bin = fCentAxis->FindBin(percentile) - 1;
	421	if (bin >= fNbinsCent) bin = -1;
	422	return bin;
07c2615a	423
	424	}
	425
	426	//====================================================================================================================================================
	427
	428	Double_t AliAnalysisTaskMuonHadronCorrelations::GetITSMultiplicity() {
	429
6de71f1d	430	Double_t multiplicity = fAOD->GetHeader()->GetNumberOfITSClusters(1);
07c2615a	431
	432	return multiplicity;
	433
	434	}
	435
	436	//====================================================================================================================================================
	437
	438	void AliAnalysisTaskMuonHadronCorrelations::Terminate(Option_t *) {
	439
	440	fOutputList = dynamic_cast<TList*> (GetOutputData(1));
	441	if (!fOutputList) {
	442	printf("ERROR: Output list not available\n");
	443	return;
	444	}
	445
	446	}
	447
	448	//====================================================================================================================================================
	449