]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGCF/Correlations/DPhi/FourierDecomposition/AliDhcTask.cxx
update DHC task (Constantin Loizides <cloizides@lbl.gov>)
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / FourierDecomposition / AliDhcTask.cxx
index 59fc84650687fec56aa3f8eb19af6abcd860a092..77aa89f142231f551150703a7d2d2ee77b4684e1 100644 (file)
@@ -7,6 +7,7 @@
 #include "TFormula.h"
 #include "TH1.h"
 #include "TH2.h"
+#include "TH3.h"
 #include "TAxis.h"
 #include "TProfile2D.h"
 #include "TROOT.h"
@@ -20,6 +21,7 @@
 #include "AliESDEvent.h"
 #include "AliESDInputHandler.h"
 #include "AliESDtrackCuts.h"
+#include "AliESDMuonTrack.h"
 #include "AliPool.h"
 #include "AliVParticle.h"
 
@@ -28,17 +30,39 @@ using std::endl;
 
 ClassImp(AliDhcTask)
 
+
+//________________________________________________________________________
+AliDhcTask::AliDhcTask()
+: AliAnalysisTaskSE(), fVerbosity(0), fEtaMax(1), fZVtxMax(10), fPtMin(0.25), fPtMax(15),
+  fTrackDepth(1000), fPoolSize(200), fTracksName(), fDoWeights(kFALSE), fFillMuons(kFALSE),
+  fEtaTLo(-1.0), fEtaTHi(1.0), fEtaALo(-1.0), fEtaAHi(1.0),
+  fESD(0x0), fAOD(0x0), fOutputList(0x0), fHEvt(0x0), fHTrk(0x0),
+  fHPtAss(0x0), fHPtTrg(0x0), fHPtTrgEvt(0x0),
+  fHPtTrgNorm1S(0x0), fHPtTrgNorm1M(0x0), fHPtTrgNorm2S(0x0), fHPtTrgNorm2M(0x0),
+  fHCent(0x0), fHZvtx(0x0), fNbins(0), fHSs(0x0), fHMs(0x0), fHPts(0x0), fIndex(0x0),
+  fCentrality(99), fZVertex(99), fEsdTPCOnly(0), fPoolMgr(0),
+  fCentMethod("V0M"), fNBdeta(20), fNBdphi(36),
+  fBPtT(0x0), fBPtA(0x0), fBCent(0x0), fBZvtx(0x0),
+  fMixBCent(0x0), fMixBZvtx(0x0),
+  fHEffT(0x0), fHEffA(0x0)
+{
+  
+}
+
 //________________________________________________________________________
 AliDhcTask::AliDhcTask(const char *name) 
 : AliAnalysisTaskSE(name), fVerbosity(0), fEtaMax(1), fZVtxMax(10), fPtMin(0.25), fPtMax(15),
-  fTrackDepth(1000), fPoolSize(200), fTracksName(), fDoWeights(0),
-  fESD(0), fAOD(0), fOutputList(0), fHistPt(0), fHEvt(0), fHTrk(0), fHPtAss(0), fHPtTrg(0),
-  fHCent(0), fHZvtx(0), fNbins(0), fHSs(0), fHMs(0),
-  fIndex(0), fMeanPtTrg(0), fMeanPtAss(0), fMean2PtTrg(0), fMean2PtAss(0),
-  fCentrality(99), fZVertex(99), fEsdTrackCutsTPCOnly(0), fPoolMgr(0),
+  fTrackDepth(1000), fPoolSize(200), fTracksName(), fDoWeights(kFALSE), fFillMuons(kFALSE),
+  fEtaTLo(-1.0), fEtaTHi(1.0), fEtaALo(-1.0), fEtaAHi(1.0),
+  fESD(0x0), fAOD(0x0), fOutputList(0x0), fHEvt(0x0), fHTrk(0x0),
+  fHPtAss(0x0), fHPtTrg(0x0), fHPtTrgEvt(0x0),
+  fHPtTrgNorm1S(0x0), fHPtTrgNorm1M(0x0), fHPtTrgNorm2S(0x0), fHPtTrgNorm2M(0x0),
+  fHCent(0x0), fHZvtx(0x0), fNbins(0), fHSs(0x0), fHMs(0x0), fHPts(0x0), fIndex(0x0),
+  fCentrality(99), fZVertex(99), fEsdTPCOnly(0), fPoolMgr(0),
   fCentMethod("V0M"), fNBdeta(20), fNBdphi(36),
   fBPtT(0x0), fBPtA(0x0), fBCent(0x0), fBZvtx(0x0),
-  fMixBCent(0x0), fMixBZvtx(0x0)
+  fMixBCent(0x0), fMixBZvtx(0x0),
+  fHEffT(0x0), fHEffA(0x0)
 {
   // Constructor
 
@@ -71,14 +95,22 @@ void AliDhcTask::UserCreateOutputObjects()
 {
   // Create histograms
   // Called once (per slave on PROOF!)
+  if (fVerbosity > 1) {
+    AliInfo("Initialize Dhc Task");
+    AliInfo(Form(" efficiency correct triggers? %d", fHEffT!=0));
+    AliInfo(Form(" efficiency correct associates? %d", fHEffA!=0));
+    AliInfo(Form(" fill muons? %d", fFillMuons));
+    AliInfo(Form(" trigger eta range %f .. %f", fEtaTLo, fEtaTHi));
+    AliInfo(Form(" associate eta range %f .. %f", fEtaALo, fEtaAHi));
+  }
 
   fOutputList = new TList();
   fOutputList->SetOwner(1);
 
-  fEsdTrackCutsTPCOnly = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
-  //fEsdTrackCutsTPCOnly->SetMinNClustersTPC(70);
-  fEsdTrackCutsTPCOnly->SetMinNCrossedRowsTPC(70);
-  fEsdTrackCutsTPCOnly->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
+  fEsdTPCOnly = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
+  //fEsdTPCOnly->SetMinNClustersTPC(70);
+  fEsdTPCOnly->SetMinNCrossedRowsTPC(70);
+  fEsdTPCOnly->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
 
   BookHistos();
   InitEventMixer(); 
@@ -91,6 +123,10 @@ void AliDhcTask::BookHistos()
   // Book histograms.
 
   if (fVerbosity > 1) {
+    AliInfo(Form("Number of pt(t) bins: %d", fBPtT->GetNbins()));
+    for (Int_t i=1; i<=fBPtT->GetNbins(); i++) {
+      AliInfo(Form("pt(t) bin %d, %f to %f", i, fBPtT->GetBinLowEdge(i), fBPtT->GetBinUpEdge(i)));
+    }
     AliInfo(Form("Number of pt(a) bins: %d", fBPtA->GetNbins()));
     for (Int_t i=1; i<=fBPtA->GetNbins(); i++) {
       AliInfo(Form("pt(a) bin %d, %f to %f", i, fBPtA->GetBinLowEdge(i), fBPtA->GetBinUpEdge(i)));
@@ -122,77 +158,76 @@ void AliDhcTask::BookHistos()
     zvtx[i] = fBZvtx->GetBinUpEdge(i);
   }
 
-  // Event histo
-  fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", 30, -15, 15, 101, 0, 101);
+  fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", nZvtx, zvtx, nCent, cent);
   fOutputList->Add(fHEvt);
-  // Track histo
   fHTrk = new TH2F("fHTrk", "Track-level variables", 
                   100, 0, TMath::TwoPi(), 100, -fEtaMax, +fEtaMax);
   fOutputList->Add(fHTrk);
   
-  // Left over from the tutorial :)
-  fHistPt = new TH1F("fHistPt", "P_{T} distribution", 200, 0., fPtMax);
-  fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
-  fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
-  fHistPt->SetMarkerStyle(kFullCircle);
-  fOutputList->Add(fHistPt);
-
   fHPtAss = new TH1F("fHPtAss","PtAssoc;P_{T} (GeV/c) [GeV/c]",nPtAssc,pta);
   fOutputList->Add(fHPtAss);
   fHPtTrg = new TH1F("fHPtTrg","PtTrig;P_{T} (GeV/c) [GeV/c]",nPtTrig,ptt);
   fOutputList->Add(fHPtTrg);
+  fHPtTrgEvt = new TH1F("fHPtTrgEvt","PtTrig;P_{T} (GeV/c) [GeV/c]",nPtTrig,ptt);
+  fHPtTrgNorm1S = new TH3F("fHPtTrgNorm1S","PtTrig;P_{T} (GeV/c) [GeV/c];centrality;z_{vtx}",nPtTrig,ptt,nCent,cent,nZvtx,zvtx);
+  fOutputList->Add(fHPtTrgNorm1S);
+  fHPtTrgNorm1M = (TH3*) fHPtTrgNorm1S->Clone("fHPtTrgNorm1M");
+  fOutputList->Add(fHPtTrgNorm1M);
+  fHPtTrgNorm2S = (TH3*) fHPtTrgNorm1S->Clone("fHPtTrgNorm2S");
+  fOutputList->Add(fHPtTrgNorm2S);
+  fHPtTrgNorm2M = (TH3*) fHPtTrgNorm1S->Clone("fHPtTrgNorm2M");
+  fOutputList->Add(fHPtTrgNorm2M);
+
   fHCent = new TH1F("fHCent","Cent;bins",nCent,cent);
   fOutputList->Add(fHCent);
   fHZvtx = new TH1F("fHZvtx","Zvertex;bins",nZvtx,zvtx);
   fOutputList->Add(fHZvtx);
 
   fNbins = nPtTrig*nPtAssc*nCent*nZvtx;
-  fHSs = new TH2*[fNbins];
-  fHMs = new TH2*[fNbins];
+  fHSs   = new TH2*[fNbins];
+  fHMs   = new TH2*[fNbins];
+  fHPts  = new TH1*[fNbins];
 
   fIndex = new TFormula("GlobIndex","(t-1)*[0]*[1]*[2]+(z-1)*[0]*[1]+(x-1)*[0]+(y-1)+0*[4]");
   fIndex->SetParameters(nPtTrig,nPtAssc,nZvtx,nCent);
   fIndex->SetParNames("NTrigBins","NAssocBins", "NZvertexBins", "NCentBins");
-  //fOutputList->Add(fIndex);
-
-  fMeanPtTrg = new TProfile2D*[nPtTrig*nPtAssc];
-  fMeanPtAss = new TProfile2D*[nPtTrig*nPtAssc];
-  fMean2PtTrg = new TProfile2D*[nPtTrig*nPtAssc];
-  fMean2PtAss = new TProfile2D*[nPtTrig*nPtAssc];
-  for (Int_t c=1; c<=nCent; ++c) {
-    TString title(Form("cen=%d (%.1f)",c,fHCent->GetBinCenter(c)));
-    fMeanPtTrg[c-1]  = new TProfile2D(Form("hMeanPtTrgCen%d",c),title,nPtTrig,ptt,nPtAssc,pta);
-    fMeanPtAss[c-1]  = new TProfile2D(Form("hMeanPtAssCen%d",c),title,nPtTrig,ptt,nPtAssc,pta);
-    fMean2PtTrg[c-1] = new TProfile2D(Form("hMean2PtTrgCen%d",c),title,nPtTrig,ptt,nPtAssc,pta);
-    fMean2PtAss[c-1] = new TProfile2D(Form("hMean2PtAssCen%d",c),title,nPtTrig,ptt,nPtAssc,pta);
-    fOutputList->Add(fMeanPtTrg[c-1]);
-    fOutputList->Add(fMeanPtAss[c-1]);
-    fOutputList->Add(fMean2PtTrg[c-1]);
-    fOutputList->Add(fMean2PtAss[c-1]);
-  }
+  fOutputList->Add(fIndex);
 
   Int_t count = 0;
   for (Int_t c=1; c<=nCent; ++c) {
     for (Int_t z=1; z<=nZvtx; ++z) {
       for (Int_t t=1; t<=nPtTrig; ++t) {
        for (Int_t a=1; a<=nPtAssc; ++a) {
-         fHSs[count] = 0;
-         fHMs[count] = 0;
+         fHSs[count]  = 0;
+         fHMs[count]  = 0;
+          fHPts[count] = 0;
          if (a>t) {
            ++count;
            continue;
          }
-         TString title(Form("cen=%d (%.1f), zVtx=%d (%.1f), trig=%d (%.1f), assc=%d (%.1f)",
-                            c,fHCent->GetBinCenter(c), z,fHZvtx->GetBinCenter(z),
-                            t,fHPtTrg->GetBinCenter(t),a, fHPtAss->GetBinCenter(a)));
-         fHSs[count] = new TH2F(Form("hS%d",count), Form("Signal %s",title.Data()),
+          if (t==1 && a==1) {
+            TString title(Form("cen=%d (%.1f to %.1f), zVtx=%d (%.1f to %.1f)",
+                               c, fBCent->GetBinLowEdge(c), fBCent->GetBinUpEdge(c), 
+                               z, fBZvtx->GetBinLowEdge(z), fBZvtx->GetBinUpEdge(z)));
+            fHPts[count] = new TH1F(Form("hPt%d",count), Form("Ptdist %s;p_{T} (GeV/c)",title.Data()), 200,0,20);
+            fOutputList->Add(fHPts[count]);
+          }
+         TString title(Form("cen=%d (%.1f to %.1f), zVtx=%d (%.1f to %.1f), trig=%d (%.1f to %.1f), assoc=%d (%.1f to %.1f)",
+                            c, fBCent->GetBinLowEdge(c), fBCent->GetBinUpEdge(c), 
+                             z, fBZvtx->GetBinLowEdge(z), fBZvtx->GetBinUpEdge(z),
+                            t, fBPtT->GetBinLowEdge(t),  fBPtT->GetBinUpEdge(t),
+                             a, fBPtA->GetBinLowEdge(a),  fBPtA->GetBinUpEdge(a)));
+         fHSs[count] = new TH2F(Form("hS%d",count), Form("Signal %s;#Delta #varphi;#Delta #eta",title.Data()),
                                 fNBdphi,-0.5*TMath::Pi(),1.5*TMath::Pi(),fNBdeta,-2*fEtaMax,2*fEtaMax);
-         fHMs[count] = new TH2F(Form("hM%d",count), Form("Signal %s",title.Data()),
+         fHMs[count] = new TH2F(Form("hM%d",count), Form("Mixed %s;#Delta #varphi;#Delta #eta",title.Data()),
                                 fNBdphi,-0.5*TMath::Pi(),1.5*TMath::Pi(),fNBdeta,-2*fEtaMax,2*fEtaMax);
          fOutputList->Add(fHSs[count]);
          fOutputList->Add(fHMs[count]);
+          Int_t tcount = (Int_t)fIndex->Eval(t,a,z,c);
          if (fVerbosity>5)
-           cout << count << " " << fIndex->Eval(t,a,z,c) << ": " << title << endl;
+           cout << count << " " << tcount << ": " << title << endl;
+          if (count != tcount)
+            AliFatal(Form("Index mismatch: %d %d", count, tcount));
          ++count;
        }
       }
@@ -200,6 +235,50 @@ void AliDhcTask::BookHistos()
   }
 }
 
+//________________________________________________________________________
+void AliDhcTask::SetAnaMode(Int_t iAna=kHH)
+{
+  if (iAna==kHH) {
+    SetFillMuons(kFALSE);
+    fEtaTLo = -1.0;
+    fEtaTHi = 1.0;
+    fEtaALo = -1.0;
+    fEtaAHi = 1.0;
+  } else if (iAna==kMuH) {
+    SetFillMuons(kTRUE);
+    fEtaTLo = -5.0;
+    fEtaTHi = -2.0;
+    fEtaALo = -1.0;
+    fEtaAHi = 1.0;
+  } else if (iAna==kHMu) {
+    SetFillMuons(kTRUE);
+    fEtaTLo = -1.0;
+    fEtaTHi = 1.0;
+    fEtaALo = -5.0;
+    fEtaAHi = -2.0;
+  } else if (iAna==kMuMu) {
+    SetFillMuons(kTRUE);
+    fEtaTLo = -5.0;
+    fEtaTHi = -2.0;
+    fEtaALo = -5.0;
+    fEtaAHi = -2.0;
+  } else if (iAna==kPSide) {
+    SetFillMuons(kFALSE);
+    fEtaTLo = -1.0;
+    fEtaTHi = -0.465;
+    fEtaALo = -1.0;
+    fEtaAHi = -0.465;
+  } else if (iAna==kASide) {
+    SetFillMuons(kFALSE);
+    fEtaTLo = -0.465;
+    fEtaTHi = 1.0;
+    fEtaALo = -0.465;
+    fEtaAHi = 1.0;
+  } else {
+    AliInfo(Form("Unrecognized analysis option: %d", iAna));
+  }
+}
+
 //________________________________________________________________________
 void AliDhcTask::InitEventMixer()
 {
@@ -235,7 +314,6 @@ void AliDhcTask::InitEventMixer()
 void AliDhcTask::UserExec(Option_t *) 
 {
   // Main loop, called for each event.
-
   static int iEvent = -1; ++iEvent;
 
   if (fVerbosity>2) {
@@ -268,7 +346,9 @@ void AliDhcTask::UserExec(Option_t *)
   if (mcgen) {
     fZVertex = 0;
     TList *list = InputEvent()->GetList();
-    TClonesArray *tcaTracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracksName));
+    TClonesArray *tcaTracks = 0;
+    if (list) 
+      tcaTracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracksName));
     if (tcaTracks) 
       fCentrality = tcaTracks->GetEntries();
   } else {
@@ -302,9 +382,9 @@ void AliDhcTask::UserExec(Option_t *)
 
   // Fill Event histogram
   fHEvt->Fill(fZVertex, fCentrality);
-
   if (fCentrality > fBCent->GetXmax() || fCentrality < fBCent->GetXmin()) {
-    AliInfo(Form("Event REJECTED (centrality out of range). fCentrality = %.1f", fCentrality));
+    if (fVerbosity > 1)
+      AliInfo(Form("Event REJECTED (centrality out of range). fCentrality = %.1f", fCentrality));
     return;
   }
 
@@ -383,7 +463,7 @@ void AliDhcTask::GetESDTracks(MiniEvent* miniEvt)
         AliError(Form("Couldn't get ESD track %d\n", i));
         continue;
       }
-      Bool_t trkOK = fEsdTrackCutsTPCOnly->AcceptTrack(esdtrack);
+      Bool_t trkOK = fEsdTPCOnly->AcceptTrack(esdtrack);
       if (!trkOK)
         continue;
       Double_t pt = esdtrack->Pt();
@@ -427,30 +507,73 @@ void AliDhcTask::GetESDTracks(MiniEvent* miniEvt)
       arr.Add(newtrack);
       nSelTrax++;
     }
-    return;
-  }
 
-  TList *list = InputEvent()->GetList();
-  TClonesArray *tcaTracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracksName));
-  const Int_t ntracks = tcaTracks->GetEntries();
-  if (miniEvt)
-    miniEvt->reserve(ntracks);
-  else {
-    AliError("Ptr to miniEvt zero");
-    return;
+    for(Int_t itrack = 0; itrack < nSelTrax; itrack++) {
+      AliVTrack *esdtrack = static_cast<AliESDtrack*>(arr.At(itrack));
+      if(!esdtrack) {
+        AliError(Form("ERROR: Could not retrieve esdtrack %d",itrack));
+        continue;
+      }
+      Double_t pt = esdtrack->Pt();
+      Double_t eta  = esdtrack->Eta();
+      Double_t phi  = esdtrack->Phi();
+      Int_t    sign = esdtrack->Charge() > 0 ? 1 : -1;
+      miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
+    }
+  } else {
+    TList *list = InputEvent()->GetList();
+    TClonesArray *tcaTracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracksName));
+
+    if(!tcaTracks){
+      AliError("Ptr to tcaTracks zero");
+      return;
+    }
+
+    const Int_t ntracks = tcaTracks->GetEntries();
+    if (miniEvt)
+      miniEvt->reserve(ntracks);
+    else {
+      AliError("Ptr to miniEvt zero");
+      return;
+    }
+
+    for (Int_t itrack = 0; itrack < ntracks; itrack++) {
+      AliVTrack *vtrack = static_cast<AliVTrack*>(tcaTracks->At(itrack));
+      if (!vtrack) {
+        AliError(Form("ERROR: Could not retrieve track %d",itrack));
+        continue;
+      }
+      Double_t pt   = vtrack->Pt();
+      Double_t eta  = vtrack->Eta();
+      Double_t phi  = vtrack->Phi();
+      Int_t    sign = vtrack->Charge() > 0 ? 1 : -1;
+      miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
+    }
   }
 
-  for(Int_t itrack = 0; itrack < ntracks; itrack++) {
-    AliVTrack *esdtrack = static_cast<AliESDtrack*>(tcaTracks->At(itrack));
-    if(!esdtrack) {
-      AliError(Form("ERROR: Could not retrieve esdtrack %d",itrack));
-      continue;
+  if (fFillMuons) {
+    // count good muons
+    Int_t nGoodMuons = 0;
+    for (Int_t iMu = 0; iMu<fESD->GetNumberOfMuonTracks(); iMu++) {
+      AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iMu);
+      if (muonTrack) {
+          if (IsGoodMUONtrack(*muonTrack)) nGoodMuons++;
+      }
+    }
+    miniEvt->reserve(miniEvt->size()+nGoodMuons);
+    // fill them into the mini event
+    for (Int_t iMu = 0; iMu<fESD->GetNumberOfMuonTracks(); iMu++) {
+      AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iMu);
+      if (muonTrack) {
+        if (!IsGoodMUONtrack(*muonTrack)) 
+          continue;
+        Double_t ptMu   = muonTrack->Pt();
+        Double_t etaMu  = muonTrack->Eta();
+        Double_t phiMu  = muonTrack->Phi();
+        Double_t signMu = muonTrack->Charge() > 0 ? 1 : -1;
+        miniEvt->push_back(AliMiniTrack(ptMu, etaMu, phiMu, signMu));
+      }
     }
-    Double_t pt = esdtrack->Pt();
-    Double_t eta  = esdtrack->Eta();
-    Double_t phi  = esdtrack->Phi();
-    Int_t    sign = esdtrack->Charge() > 0 ? 1 : -1;
-    miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
   }
 }
 
@@ -459,65 +582,122 @@ void AliDhcTask::GetAODTracks(MiniEvent* miniEvt)
 {
   // Loop twice: 1. Count sel. tracks. 2. Fill vector.
 
-  Int_t nTrax = fAOD->GetNumberOfTracks();
-  Int_t nSelTrax = 0;
+  if (fTracksName.IsNull()) {
+    Int_t nTrax = fAOD->GetNumberOfTracks();
+    Int_t nSelTrax = 0;
 
-  if (fVerbosity > 2)
-    AliInfo(Form("%d tracks in event",nTrax));
+    if (fVerbosity > 2)
+      AliInfo(Form("%d tracks in event",nTrax));
 
-  // Loop 1.
-  for (Int_t i = 0; i < nTrax; ++i) {
-    AliAODTrack* aodtrack = fAOD->GetTrack(i);
-    if (!aodtrack) {
-      AliError(Form("Couldn't get AOD track %d\n", i));
-      continue;
+    // Loop 1.
+    for (Int_t i = 0; i < nTrax; ++i) {
+      AliAODTrack* aodtrack = fAOD->GetTrack(i);
+      if (!aodtrack) {
+        AliError(Form("Couldn't get AOD track %d\n", i));
+        continue;
+      }
+      // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
+      UInt_t tpcOnly = 1 << 7;
+      Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly);
+      if (!trkOK)
+        continue;
+      Double_t pt = aodtrack->Pt();
+      Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
+      if (!ptOK)
+        continue;
+      Double_t eta = aodtrack->Eta();
+      if (TMath::Abs(eta) > fEtaMax)
+        continue;
+      nSelTrax++;
     }
-    // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
-    UInt_t tpcOnly = 1 << 7;
-    Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly);
-    if (!trkOK)
-      continue;
-    Double_t pt = aodtrack->Pt();
-    Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
-    if (!ptOK)
-      continue;
-    Double_t eta = aodtrack->Eta();
-    if (TMath::Abs(eta) > fEtaMax)
-      continue;
-    nSelTrax++;
-  }
 
-  if (miniEvt)
-    miniEvt->reserve(nSelTrax);
-  else {
-    AliError("!miniEvt");
-    return;
-  }
-  
-  // Loop 2.  
-  for (Int_t i = 0; i < nTrax; ++i) {
-    AliAODTrack* aodtrack = fAOD->GetTrack(i);
-    if (!aodtrack) {
-      AliError(Form("Couldn't get AOD track %d\n", i));
-      continue;
+    if (miniEvt)
+      miniEvt->reserve(nSelTrax);
+    else {
+      AliError("!miniEvt");
+      return;
     }
+  
+    // Loop 2.  
+    for (Int_t i = 0; i < nTrax; ++i) {
+      AliAODTrack* aodtrack = fAOD->GetTrack(i);
+      if (!aodtrack) {
+        AliError(Form("Couldn't get AOD track %d\n", i));
+        continue;
+      }
     
-    // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
-    UInt_t tpcOnly = 1 << 7;
-    Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly);
-    if (!trkOK)
-      continue;
-    Double_t pt = aodtrack->Pt();
-    Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
-    if (!ptOK)
-      continue;
-    Double_t eta  = aodtrack->Eta();
-    if (TMath::Abs(eta) > fEtaMax)
-      continue;
+      // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
+      UInt_t tpcOnly = 1 << 7;
+      Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly);
+      if (!trkOK)
+        continue;
+      Double_t pt = aodtrack->Pt();
+      Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
+      if (!ptOK)
+        continue;
+      Double_t eta  = aodtrack->Eta();
+      if (TMath::Abs(eta) > fEtaMax)
+        continue;
+      
+      Double_t phi  = aodtrack->Phi();
+      Int_t    sign = aodtrack->Charge() > 0 ? 1 : -1;
+      miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
+    }
+  } else {
+    TList *list = InputEvent()->GetList();
+    TClonesArray *tcaTracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracksName));
+
+    if (!tcaTracks){
+      AliError("Ptr to tcaTracks zero");
+      return;
+    }
+
+    const Int_t ntracks = tcaTracks->GetEntries();
+    if (miniEvt)
+      miniEvt->reserve(ntracks);
+    else {
+      AliError("Ptr to miniEvt zero");
+      return;
+    }
 
-    Double_t phi  = aodtrack->Phi();
-    Int_t    sign = aodtrack->Charge() > 0 ? 1 : -1;
-    miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
+    for (Int_t itrack = 0; itrack < ntracks; itrack++) {
+      AliVTrack *vtrack = static_cast<AliVTrack*>(tcaTracks->At(itrack));
+      if (!vtrack) {
+        AliError(Form("ERROR: Could not retrieve vtrack %d",itrack));
+        continue;
+      }
+      Double_t pt   = vtrack->Pt();
+      Double_t eta  = vtrack->Eta();
+      Double_t phi  = vtrack->Phi();
+      Int_t    sign = vtrack->Charge() > 0 ? 1 : -1;
+      miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
+    }
+  }
+
+  if (fFillMuons) {
+    // count good muons
+    Int_t nGoodMuons = 0;
+    for (Int_t iMu = 0; iMu<fAOD->GetNumberOfTracks(); iMu++) {
+      AliAODTrack* muonTrack = fAOD->GetTrack(iMu);
+      if(muonTrack) {
+        if (IsGoodMUONtrack(*muonTrack)) 
+          nGoodMuons++;
+      }
+    }
+    miniEvt->reserve(miniEvt->size()+nGoodMuons);
+    // fill them into the mini event
+    for (Int_t iMu = 0; iMu<fAOD->GetNumberOfTracks(); iMu++) {
+      AliAODTrack* muonTrack = fAOD->GetTrack(iMu);
+      if (muonTrack) {
+        if (!IsGoodMUONtrack(*muonTrack)) 
+          continue;
+        Double_t ptMu   = muonTrack->Pt();
+        Double_t etaMu  = muonTrack->Eta();
+        Double_t phiMu  = muonTrack->Phi();
+        Double_t signMu = muonTrack->Charge() > 0 ? 1 : -1;
+        miniEvt->push_back(AliMiniTrack(ptMu, etaMu, phiMu, signMu));
+      }
+    }
   }
 }
 
@@ -562,8 +742,8 @@ Int_t AliDhcTask::Correlate(const MiniEvent &evt1, const MiniEvent &evt2, Int_t
   TH2  **hist = fHMs;
   if (pairing == kSameEvt) {
     hist = fHSs;
-    fHCent->AddBinContent(cbin);
-    fHZvtx->AddBinContent(zbin);
+    fHCent->Fill(fCentrality);
+    fHZvtx->Fill(fZVertex);
   }
 
   Int_t nZvtx = fHZvtx->GetNbinsX();
@@ -572,69 +752,128 @@ Int_t AliDhcTask::Correlate(const MiniEvent &evt1, const MiniEvent &evt2, Int_t
 
   Int_t globIndex = (cbin-1)*nZvtx*nPtTrig*nPtAssc+(zbin-1)*nPtTrig*nPtAssc;
 
-  Double_t weight = 1;
-  if (fDoWeights) { // Count trigger particles in this event
-    for (Int_t i=0; i<iMax; ++i) {
-      const AliMiniTrack &a(evt1.at(i));
-      Float_t pta  = a.Pt();
-      Int_t abin = fHPtTrg->FindBin(pta);
-      if (fHPtTrg->IsBinOverflow(abin) ||
-          fHPtTrg->IsBinUnderflow(abin))
-        continue;
-      ++weight;
+  Int_t ptindex = (Int_t)fIndex->Eval(1,1,zbin,cbin);
+
+  fHPtTrgEvt->Reset();
+  for (Int_t i=0; i<iMax; ++i) {
+    const AliMiniTrack &a(evt1.at(i));
+    Float_t pta  = a.Pt();
+    fHPtTrgEvt->Fill(pta);
+    if (pairing == kSameEvt) {
+      fHPts[ptindex]->Fill(pta);
     }
   }
-  
-  for (Int_t i=0; i<iMax; ++i) {
 
+  Bool_t bCountTrg; // only count the trigger if an associated particle is found
+
+  for (Int_t i=0; i<iMax; ++i) {
     // Trigger particles
     const AliMiniTrack &a(evt1.at(i));
 
     Float_t pta  = a.Pt();
+    Float_t etaa = a.Eta();
+    Float_t phia = a.Phi();
     Int_t abin = fHPtTrg->FindBin(pta);
-    if (fHPtTrg->IsBinOverflow(abin) ||
-       fHPtTrg->IsBinUnderflow(abin))
+    if (fHPtTrg->IsBinOverflow(abin) || fHPtTrg->IsBinUnderflow(abin))
+      continue;
+
+    if (etaa<fEtaTLo || etaa>fEtaTHi)
       continue;
 
+    // efficiency weighting
+    Double_t effWtT = 1.0;
+    if (fHEffT) {
+      // trigger particle
+      const Int_t nEffDimT = fHEffT->GetNdimensions();
+      Int_t effBinT[nEffDimT];
+      effBinT[0] = fHEffT->GetAxis(0)->FindBin(etaa);
+      effBinT[1] = fHEffT->GetAxis(1)->FindBin(pta);
+      effBinT[2] = fHEffT->GetAxis(2)->FindBin(fCentrality);
+      effBinT[3] = fHEffT->GetAxis(3)->FindBin(fZVertex);
+      if (nEffDimT>4) {
+        effBinT[4] = fHEffT->GetAxis(4)->FindBin(phia);
+      }
+      effWtT = fHEffT->GetBinContent(effBinT);
+    }
+    
     if (pairing == kSameEvt) {
-      fHistPt->Fill(pta);
-      fHTrk->Fill(a.Phi(),a.Eta());
-      fHPtTrg->AddBinContent(abin);
+      fHTrk->Fill(phia,etaa);
+      fHPtTrgNorm1S->Fill(pta,fCentrality,fZVertex,effWtT);
+    } else {
+      fHPtTrgNorm1M->Fill(pta,fCentrality,fZVertex,effWtT);
     }
+    bCountTrg = kFALSE;
 
     for (Int_t j=0; j<jMax; ++j) {
       // Associated particles
       if (pairing == kSameEvt && i==j)
-       continue;
+        continue;
 
       const AliMiniTrack &b(evt2.at(j));
       
       Float_t ptb  = b.Pt();
-      if (pta < ptb) 
-       continue;
+      Float_t etab = b.Eta();
+      Float_t phib = b.Phi();
+      if (pta < ptb)
+           continue;
+
+      Int_t bbin = fHPtAss->FindBin(ptb);
+      if (fHPtAss->IsBinOverflow(bbin) || fHPtAss->IsBinUnderflow(bbin))
+        continue;
 
-      Int_t bbin = fHPtTrg->FindBin(ptb);
-      if (fHPtAss->IsBinOverflow(bbin) ||
-         fHPtAss->IsBinUnderflow(bbin))
-       continue;
+      if (etab<fEtaALo || etab>fEtaAHi)
+        continue;
 
-      Float_t dphi = DeltaPhi(a.Phi(), b.Phi());
-      Float_t deta = a.Eta() - b.Eta();
+      Float_t dphi = DeltaPhi(phia, phib);
+      Float_t deta = etaa - etab;
 
       Int_t index = globIndex+(abin-1)*nPtAssc+(bbin-1);
-      hist[index]->Fill(dphi,deta,1./weight);
+      Double_t weight = 1.0;
+      // number of trigger weight event-by-event (CMS method)
+      if (fDoWeights) {
+        Double_t nTrgs = fHPtTrgEvt->GetBinContent(abin);
+        if (nTrgs==0.0) {
+          AliError(Form("Trying to work with number of triggers weight = %f",nTrgs));
+          return 0;
+        }
+        weight *= 1./nTrgs;
+      }
+
+      // efficiency weighting
+      if (fHEffT) {
+        // trigger particle
+        weight *= effWtT;
+      }
+      if (fHEffA) {
+        // associated particle
+        const Int_t nEffDimA = fHEffA->GetNdimensions();
+        Int_t effBinA[nEffDimA];
+        effBinA[0] = fHEffA->GetAxis(0)->FindBin(etab);
+        effBinA[1] = fHEffA->GetAxis(1)->FindBin(ptb);
+        effBinA[2] = fHEffA->GetAxis(2)->FindBin(fCentrality);
+        effBinA[3] = fHEffA->GetAxis(3)->FindBin(fZVertex);
+        if (nEffDimA>4) {
+          effBinA[4] = fHEffA->GetAxis(4)->FindBin(phib);
+        }
+        weight *= fHEffA->GetBinContent(effBinA);
+      }
+      hist[index]->Fill(dphi,deta,weight);
+      bCountTrg = kTRUE;
 
       if (pairing == kSameEvt) {
-       fHPtAss->AddBinContent(bbin);
-       fMeanPtTrg[cbin-1]->Fill(pta,ptb,pta);
-       fMeanPtAss[cbin-1]->Fill(pta,ptb,ptb);
-       fMean2PtTrg[cbin-1]->Fill(pta,ptb,pta*pta);
-       fMean2PtAss[cbin-1]->Fill(pta,ptb,ptb*ptb);
+        fHPtAss->Fill(ptb);
+      }
+    }
+    if (bCountTrg) {
+      if (pairing == kSameEvt) {
+        fHPtTrgNorm2S->Fill(pta,fCentrality,fZVertex,effWtT);
+      } else {
+        fHPtTrgNorm2M->Fill(pta,fCentrality,fZVertex,effWtT);
       }
     }
   }
 
-  return 0;
+  return 1;
 }
 
 //________________________________________________________________________
@@ -642,31 +881,6 @@ void AliDhcTask::Terminate(Option_t *)
 {
   // Draw result to the screen
   // Called once at the end of the query
-
-  delete fPoolMgr;
-
-  fHCent->SetEntries(fHCent->Integral());
-  fHZvtx->SetEntries(fHZvtx->Integral());
-  fHPtTrg->SetEntries(fHPtTrg->Integral());
-  fHPtAss->SetEntries(fHPtAss->Integral());
-
-  if (gROOT->IsBatch())
-    return;
-
-  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
-  if (!fOutputList) {
-    AliError("Output list not available");
-    return;
-  }
-  
-  fHistPt = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistPt"));
-  if (!fHistPt) {
-    AliError("ERROR: fHistPt not available\n");
-    return;
-  }
-  TCanvas *c1 = new TCanvas("AliDhcTask","Pt",10,10,510,510);
-  c1->cd(1)->SetLogy();
-  fHistPt->DrawCopy("E");
 }
 
 //________________________________________________________________________
@@ -708,3 +922,49 @@ Bool_t AliDhcTask::VertexOk(TObject* obj) const
   
   return kTRUE;
 }
+
+//________________________________________________________________________
+Bool_t AliDhcTask::IsGoodMUONtrack(AliESDMuonTrack &track)
+{
+  // Applying track cuts for MUON tracks
+
+  if (!track.ContainTrackerData()) 
+    return kFALSE;
+  if (!track.ContainTriggerData()) 
+    return kFALSE;
+  Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
+  if ((thetaTrackAbsEnd < 2.) || (thetaTrackAbsEnd > 10.)) 
+    return kFALSE;
+  Double_t eta = track.Eta();
+  if ((eta < -4.) || (eta > -2.5))
+    return kFALSE;
+  if (track.GetMatchTrigger() < 0.5) 
+    return kFALSE;
+  return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliDhcTask::IsGoodMUONtrack(AliAODTrack &track)
+{
+  // Applying track cuts for MUON tracks
+
+  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;
+}
+
+//________________________________________________________________________
+AliDhcTask::~AliDhcTask()
+{
+  //Destructor
+  if (fPoolMgr) 
+    delete fPoolMgr;
+}