]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGCF/Correlations/DPhi/AliAnalysisTaskPhiCorrelations.cxx
cuts for TOF requirements for MC corrections (Leonardo)
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / AliAnalysisTaskPhiCorrelations.cxx
index 4a9365b206395f27577ac95fa81751f8ccfb8bb0..92fe4281b3faf3cd2f3498ab14406c5d3cb88e88 100644 (file)
@@ -22,6 +22,7 @@
 #include <TMath.h>
 #include <TTree.h>
 #include <TH2F.h>
+#include <TH3F.h>
 #include <TRandom.h>
 
 #include "AliAnalysisTaskPhiCorrelations.h"
 #include "AliCentrality.h"
 #include "AliStack.h"
 #include "AliAODMCHeader.h"
+#include "AliGenCocktailEventHeader.h"
+#include "AliGenEventHeader.h"
 
 #include "AliEventPoolManager.h"
 
+#include "AliESDZDC.h"
+#include "AliESDtrackCuts.h"
+
+#include "AliHelperPID.h"
 
 ////////////////////////////////////////////////////////////////////////
 //
@@ -86,12 +93,16 @@ fCompareCentralities(kFALSE),
 fTwoTrackEfficiencyStudy(kFALSE),
 fTwoTrackEfficiencyCut(0),
 fUseVtxAxis(kFALSE),
+fCourseCentralityBinning(kFALSE),
 fSkipTrigger(kFALSE),
+fInjectedSignals(kFALSE),
 // pointers to UE classes
+fHelperPID(0x0),
 fAnalyseUE(0x0),
 fHistos(0x0),
 fHistosMixed(0),
-fkTrackingEfficiency(0x0),
+fEfficiencyCorrectionTriggers(0),
+fEfficiencyCorrectionAssociated(0),
 // handlers and events
 fAOD(0x0),
 fESD(0x0),
@@ -111,15 +122,26 @@ fTrackEtaCut(0.8),
 fOnlyOneEtaSide(0),
 fPtMin(0.5),
 fFilterBit(0xFF),
-fSelectBit(0),
+fTrackStatus(0),
+fSelectBit(AliVEvent::kMB|AliVEvent::kUserDefined),
 fUseChargeHadrons(kFALSE),
+fParticleSpeciesTrigger(-1),
+fParticleSpeciesAssociated(-1),
 fSelectCharge(0),
+fTriggerSelectCharge(0),
+fAssociatedSelectCharge(0),
 fTriggerRestrictEta(-1),
 fEtaOrdering(kFALSE),
 fCutConversions(kFALSE),
 fCutResonances(kFALSE),
 fFillOnlyStep0(kFALSE),
 fSkipStep6(kFALSE),
+fRejectCentralityOutliers(kFALSE),
+fRemoveWeakDecays(kFALSE),
+fRemoveDuplicates(kFALSE),
+fSkipFastCluster(kFALSE),
+fWeightPerEvent(kFALSE),
+fCustomBinning(),
 fFillpT(kFALSE)
 {
   // Default constructor
@@ -192,8 +214,13 @@ void  AliAnalysisTaskPhiCorrelations::CreateOutputObjects()
   // Initialize class with main algorithms, event and track selection. 
   fAnalyseUE = new AliAnalyseLeadingTrackUE();
   fAnalyseUE->SetParticleSelectionCriteria(fFilterBit, fUseChargeHadrons, fTrackEtaCut, fPtMin);
+  fAnalyseUE->SetTrackStatus(fTrackStatus);
   fAnalyseUE->SetDebug(fDebug); 
   fAnalyseUE->DefineESDCuts(fFilterBit);
+  fAnalyseUE->SetEventSelection(fSelectBit);
+  fAnalyseUE->SetHelperPID(fHelperPID);
+  if ((fParticleSpeciesTrigger != -1 || fParticleSpeciesAssociated != -1) && !fHelperPID)
+    AliFatal("HelperPID object should be set in the steering macro");
 
   // Initialize output list of containers
   if (fListOfHistos != NULL){
@@ -206,15 +233,25 @@ void  AliAnalysisTaskPhiCorrelations::CreateOutputObjects()
        }
 
   // Initialize class to handle histograms 
-  const char* histType = "4";
-  if (fUseVtxAxis)
-    histType = "5";
-  fHistos = new AliUEHistograms("AliUEHistogramsSame", histType);
-  fHistosMixed = new AliUEHistograms("AliUEHistogramsMixed", histType);
+  TString histType = "4R";
+  if (fUseVtxAxis == 1)
+    histType = "5R";
+  else if (fUseVtxAxis == 2)
+    histType = "6R";
+  if (fCourseCentralityBinning)
+    histType += "C";
+  fHistos = new AliUEHistograms("AliUEHistogramsSame", histType, fCustomBinning);
+  fHistosMixed = new AliUEHistograms("AliUEHistogramsMixed", histType, fCustomBinning);
   
   fHistos->SetSelectCharge(fSelectCharge);
   fHistosMixed->SetSelectCharge(fSelectCharge);
   
+  fHistos->SetSelectTriggerCharge(fTriggerSelectCharge);
+  fHistosMixed->SetSelectTriggerCharge(fTriggerSelectCharge);
+  
+  fHistos->SetSelectAssociatedCharge(fAssociatedSelectCharge);
+  fHistosMixed->SetSelectAssociatedCharge(fAssociatedSelectCharge);
+
   fHistos->SetTriggerRestrictEta(fTriggerRestrictEta);
   fHistosMixed->SetTriggerRestrictEta(fTriggerRestrictEta);
   
@@ -227,15 +264,36 @@ void  AliAnalysisTaskPhiCorrelations::CreateOutputObjects()
   fHistos->SetPairCuts(fCutConversions, fCutResonances);
   fHistosMixed->SetPairCuts(fCutConversions, fCutResonances);
   
+  fHistos->SetTrackEtaCut(fTrackEtaCut);
+  fHistosMixed->SetTrackEtaCut(fTrackEtaCut);
+  
+  fHistos->SetWeightPerEvent(fWeightPerEvent);
+  fHistosMixed->SetWeightPerEvent(fWeightPerEvent);
+  
+  if (fEfficiencyCorrectionTriggers)
+   {
+    fHistos->SetEfficiencyCorrectionTriggers(fEfficiencyCorrectionTriggers);
+    fHistosMixed->SetEfficiencyCorrectionTriggers((THnF*) fEfficiencyCorrectionTriggers->Clone());
+   }
+  if (fEfficiencyCorrectionAssociated)
+  {
+    fHistos->SetEfficiencyCorrectionAssociated(fEfficiencyCorrectionAssociated);
+    fHistosMixed->SetEfficiencyCorrectionAssociated((THnF*) fEfficiencyCorrectionAssociated->Clone());
+  }
+  
   // add histograms to list
   fListOfHistos->Add(fHistos);
   fListOfHistos->Add(fHistosMixed);
+  // add HelperPID to list
+  if (fHelperPID)
+    fListOfHistos->Add(fHelperPID);
   
   fListOfHistos->Add(new TH2F("trackletsVsV0Cent", ";L1 clusters;v0 centrality", 100, -0.5, 9999.5, 101, 0, 101));
   fListOfHistos->Add(new TH2F("processIDs", ";#Delta#phi;process id", 100, -0.5 * TMath::Pi(), 1.5 * TMath::Pi(), kPNoProcess + 1, -0.5, kPNoProcess + 0.5));
   fListOfHistos->Add(new TH1F("eventStat", ";;events", 4, -0.5, 3.5));
   fListOfHistos->Add(new TH2F("mixedDist", ";centrality;tracks;events", 101, 0, 101, 200, 0, fMixingTracks * 1.5));
   fListOfHistos->Add(new TH1F("pids", ";pdg;tracks", 2001, -1000.5, 1000.5));
+  fListOfHistos->Add(new TH2F("referenceMultiplicity", ";centrality;tracks;events", 101, 0, 101, 200, 0, 200));
   
   PostData(0,fListOfHistos);
   
@@ -243,24 +301,30 @@ void  AliAnalysisTaskPhiCorrelations::CreateOutputObjects()
   AddSettingsTree();
 
   // event mixing
-  Int_t trackDepth = fMixingTracks; 
-  Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
+  Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemention of AliEventPoolManager
    
   Int_t nCentralityBins  = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(1);
   Double_t* centralityBins = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(1, 0)->GetXbins()->GetArray();
   
-  Int_t nZvtxBins  = 7+1+7;
-  // bins for second buffer are shifted by 100 cm
-  Double_t vertexBins[] = { -7, -5, -3, -1, 1, 3, 5, 7, 93, 95, 97, 99, 101, 103, 105, 107 };
+  const Int_t kNZvtxBins  = 10+(1+10)*4;
+  // 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, 
+                                      90,  92,  94,  96,  98, 100, 102, 104, 106, 108, 110, 
+                                     190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210, 
+                                     290, 292, 294, 296, 298, 300, 302, 304, 306, 308, 310, 
+                                     390, 392, 394, 396, 398, 400, 402, 404, 406, 408, 410 };
+
+  Int_t nZvtxBins  = kNZvtxBins;
   Double_t* zvtxbin = vertexBins;
 
-  if (fHistos->GetUEHist(2)->GetEventHist()->GetNVar() > 2)
+  if (fMode == 0 && fHistos->GetUEHist(2)->GetEventHist()->GetNVar() > 2)
   {
     nZvtxBins = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(2);
     zvtxbin = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(2, 0)->GetXbins()->GetArray();
   }
 
-  fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centralityBins, nZvtxBins, zvtxbin);
+  fPoolMgr = new AliEventPoolManager(poolsize, fMixingTracks, nCentralityBins, centralityBins, nZvtxBins, zvtxbin);
+  fPoolMgr->SetTargetValues(fMixingTracks, 0.1, 5);
 }
 
 //____________________________________________________________________
@@ -308,17 +372,28 @@ void  AliAnalysisTaskPhiCorrelations::AddSettingsTree()
   settingsTree->Branch("fOnlyOneEtaSide", &fOnlyOneEtaSide,"OnlyOneEtaSide/I");
   settingsTree->Branch("fPtMin", &fPtMin, "PtMin/D");
   settingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/I");
+  settingsTree->Branch("fTrackStatus", &fTrackStatus,"TrackStatus/I");
   settingsTree->Branch("fSelectBit", &fSelectBit,"EventSelectionBit/I");
   settingsTree->Branch("fUseChargeHadrons", &fUseChargeHadrons,"UseChHadrons/O");
+  settingsTree->Branch("fParticleSpeciesTrigger", &fParticleSpeciesTrigger,"ParticleSpeciesTrigger/I");
+  settingsTree->Branch("fParticleSpeciesAssociated", &fParticleSpeciesAssociated,"ParticleSpeciesAssociated/I");
   settingsTree->Branch("fSelectCharge", &fSelectCharge,"SelectCharge/I");
+  settingsTree->Branch("fTriggerSelectCharge", &fTriggerSelectCharge,"TriggerSelectCharge/I");
+  settingsTree->Branch("fAssociatedSelectCharge", &fAssociatedSelectCharge,"fAssociatedSelectCharge/I");
   settingsTree->Branch("fTriggerRestrictEta", &fTriggerRestrictEta,"TriggerRestrictEta/D");
   settingsTree->Branch("fEtaOrdering", &fEtaOrdering,"EtaOrdering/O");
   settingsTree->Branch("fCutConversions", &fCutConversions,"CutConversions/O");
   settingsTree->Branch("fCutResonances", &fCutResonances,"CutResonances/O");
   settingsTree->Branch("fFillpT", &fFillpT,"FillpT/O");
-  settingsTree->Branch("fkTrackingEfficiency", "TH1D", &fkTrackingEfficiency);
   settingsTree->Branch("fMixingTracks", &fMixingTracks,"MixingTracks/I");
   settingsTree->Branch("fSkipTrigger", &fSkipTrigger,"SkipTrigger/O");
+  settingsTree->Branch("fInjectedSignals", &fInjectedSignals,"SkipTrigger/O");
+  settingsTree->Branch("fRejectCentralityOutliers", &fRejectCentralityOutliers,"RejectCentralityOutliers/O");
+  settingsTree->Branch("fRemoveWeakDecays", &fRemoveWeakDecays,"RemoveWeakDecays/O");
+  settingsTree->Branch("fRemoveDuplicates", &fRemoveDuplicates,"RemoveDuplicates/O");
+  settingsTree->Branch("fSkipFastCluster", &fSkipFastCluster,"SkipFastCluster/O");
+  settingsTree->Branch("fWeightPerEvent", &fWeightPerEvent,"WeightPerEvent/O");
+  //fCustomBinning
   
   settingsTree->Fill();
   fListOfHistos->Add(settingsTree);
@@ -374,6 +449,9 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
   // count all events
   fHistos->FillEvent(centrality, -1);
   
+  if (centrality < 0)
+    return;
+
   // Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
   TObject* vertexSupplier = fMcEvent;
   if (fAOD) // AOD
@@ -390,10 +468,72 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
   Float_t weight = 1;
   if (fFillpT)
     weight = -1;
+  
+  // For productions with injected signals, figure out above which label to skip particles/tracks
+  Int_t skipParticlesAbove = 0;
+  if (fInjectedSignals)
+  {
+    AliGenEventHeader* eventHeader = 0;
+    Int_t headers = 0;
+
+    if (fMcEvent)
+    {
+      // ESD
+      AliHeader* header = (AliHeader*) fMcEvent->Header();
+      if (!header)
+       AliFatal("fInjectedSignals set but no MC header found");
+       
+      AliGenCocktailEventHeader* cocktailHeader = dynamic_cast<AliGenCocktailEventHeader*> (header->GenEventHeader());
+      if (!cocktailHeader)
+      {
+       header->Dump();
+       AliFatal("fInjectedSignals set but no MC cocktail header found");
+      }
+
+      headers = cocktailHeader->GetHeaders()->GetEntries();
+      eventHeader = dynamic_cast<AliGenEventHeader*> (cocktailHeader->GetHeaders()->First());
+    }
+    else
+    {
+      // AOD
+      AliAODMCHeader* header = (AliAODMCHeader*) fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
+      if (!header)
+       AliFatal("fInjectedSignals set but no MC header found");
+      
+      headers = header->GetNCocktailHeaders();
+      eventHeader = header->GetCocktailHeader(0);
+    }
     
+    if (!eventHeader)
+    {
+      // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing 
+      // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
+      AliError("First event header not found. Skipping this event.");
+      fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
+      return;
+    }
+    
+    skipParticlesAbove = eventHeader->NProduced();
+    AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
+  }
+  
   // Get MC primaries
-  TObjArray* tracksMC = fAnalyseUE->GetAcceptedParticles(mc, 0, kTRUE, -1, kTRUE);
+  // triggers
+  TObjArray* tmpList = fAnalyseUE->GetAcceptedParticles(mc, 0, kTRUE, fParticleSpeciesTrigger, kTRUE);
+  CleanUp(tmpList, mc, skipParticlesAbove);
+  TObjArray* tracksMC = CloneAndReduceTrackList(tmpList);
+  delete tmpList;
   
+  // associated
+  TObjArray* tracksCorrelateMC = tracksMC;
+  if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
+  {
+    // TODO for MC this uses to PDG of the mother of the particle
+    tracksCorrelateMC = fAnalyseUE->GetAcceptedParticles(mc, 0, kTRUE, fParticleSpeciesAssociated, kTRUE);
+    CleanUp(tracksCorrelateMC, mc, skipParticlesAbove);
+  }
+  
+  /*
   if (fAOD)
   {
     for (Int_t i=0; i<fArrayMC->GetEntriesFast(); i++)
@@ -404,20 +544,23 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
     for (Int_t i=0; i<fMcEvent->GetNumberOfTracks(); i++)
       ((TH1F*) fListOfHistos->FindObject("pids"))->Fill(fMcEvent->GetTrack(i)->PdgCode());
   }
+  */
+  
+  if (fFillOnlyStep0)
+    zVtx = 0;
   
   // (MC-true all particles)
   // STEP 0
-  fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepAll, tracksMC, 0, weight);
+  fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepAll, tracksMC, tracksCorrelateMC, weight);
   
   // mixed event
   if (fFillMixed)
   {
     AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx);
-    //pool->PrintInfo();
-    if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5) 
+    if (pool->IsReady())
       for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) 
        fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepAll, tracksMC, pool->GetEvent(jMix), 1.0 / pool->GetCurrentNEvents(), (jMix == 0));
-    pool->UpdatePool(CloneAndReduceTrackList(tracksMC));
+    pool->UpdatePool(CloneAndReduceTrackList(tracksCorrelateMC));
   }
   
 //   Printf("trigger: %d", ((AliInputEventHandler*)fInputHandler)->IsEventSelected());
@@ -428,7 +571,7 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
     // (MC-true all particles)
     // STEP 1
     if (!fReduceMemoryFootprint)
-      fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTriggered, tracksMC, 0, weight);
+      fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTriggered, tracksMC, tracksCorrelateMC, weight);
     else
       fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
       
@@ -448,8 +591,12 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
         TObjArray* primMCParticles = fAnalyseUE->GetAcceptedParticles(mc, 0x0, kTRUE, particleSpecies, kTRUE);
         TObjArray* primRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, particleSpecies, kTRUE);
         TObjArray* allRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, particleSpecies, kTRUE);
+       
+       CleanUp(primMCParticles, mc, skipParticlesAbove);
+       CleanUp(primRecoTracksMatched, mc, skipParticlesAbove);
+       CleanUp(allRecoTracksMatched, mc, skipParticlesAbove);
       
-        fHistos->FillTrackingEfficiency(primMCParticles, primRecoTracksMatched, allRecoTracksMatched, particleSpecies, centrality);
+        fHistos->FillTrackingEfficiency(primMCParticles, primRecoTracksMatched, allRecoTracksMatched, 0, particleSpecies, centrality, zVtx);
         
 //     Printf("%d --> %d %d %d", particleSpecies, primMCParticles->GetEntries(), primRecoTracksMatched->GetEntries(), allRecoTracksMatched->GetEntries());
 
@@ -457,58 +604,134 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
         delete primRecoTracksMatched;
         delete allRecoTracksMatched;
       }
+      
+      TObjArray* fakeParticles = fAnalyseUE->GetFakeParticles(inputEvent, mc, kFALSE, -1, kTRUE);
+      CleanUp((TObjArray*) fakeParticles->At(0), mc, skipParticlesAbove);
+      CleanUp((TObjArray*) fakeParticles->At(1), mc, skipParticlesAbove);
+
+      fHistos->FillTrackingEfficiency(0, 0, 0, (TObjArray*) fakeParticles->At(2), -1, centrality, zVtx);
+      fHistos->FillFakePt(fakeParticles, centrality);
+//       Printf(">>>>> %d %d %d fakes", ((TObjArray*) fakeParticles->At(0))->GetEntriesFast(), ((TObjArray*) fakeParticles->At(1))->GetEntriesFast(), ((TObjArray*) fakeParticles->At(2))->GetEntriesFast());
+      delete fakeParticles;
     
       // (MC-true all particles)
       // STEP 2
       if (!fReduceMemoryFootprint)
-       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepVertex, tracksMC, 0, weight);
+       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepVertex, tracksMC, tracksCorrelateMC, weight);
       else
        fHistos->FillEvent(centrality, AliUEHist::kCFStepVertex);
       
       // Get MC primaries that match reconstructed track
-      TObjArray* tracksRecoMatchedPrim = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, -1, kTRUE);
-      
+      // triggers
+      TObjArray* tracksRecoMatchedPrim = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, fParticleSpeciesTrigger, kTRUE);
+      CleanUp(tracksRecoMatchedPrim, mc, skipParticlesAbove);
+      // associated
+      TObjArray* tracksCorrelateRecoMatchedPrim = tracksRecoMatchedPrim;
+      if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
+      {
+       tracksCorrelateRecoMatchedPrim = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, fParticleSpeciesAssociated, kTRUE);
+       CleanUp(tracksCorrelateRecoMatchedPrim, mc, skipParticlesAbove);
+      }
+
       // (RECO-matched (quantities from MC particle) primary particles)
       // STEP 4
-      fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTrackedOnlyPrim, tracksRecoMatchedPrim, 0, weight);
+      fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTrackedOnlyPrim, tracksRecoMatchedPrim, tracksCorrelateRecoMatchedPrim, weight);
+
+      // mixed event
+      if (fFillMixed)
+      {
+       AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx + 200);
+       if (pool->IsReady())
+         for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) 
+           fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTrackedOnlyPrim, tracksRecoMatchedPrim, pool->GetEvent(jMix), 1.0 / pool->GetCurrentNEvents(), (jMix == 0));
+       pool->UpdatePool(CloneAndReduceTrackList(tracksCorrelateRecoMatchedPrim));
+      }
       
       // Get MC primaries + secondaries that match reconstructed track
-      TObjArray* tracksRecoMatchedAll = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, -1, kTRUE);
+      // triggers
+      TObjArray* tracksRecoMatchedAll = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, fParticleSpeciesTrigger, kTRUE);
+      CleanUp(tracksRecoMatchedAll, mc, skipParticlesAbove);
+      // associated
+      TObjArray* tracksCorrelateRecoMatchedAll = tracksRecoMatchedAll;
+      if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
+      {
+       tracksCorrelateRecoMatchedAll = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, fParticleSpeciesAssociated, kTRUE);
+       CleanUp(tracksCorrelateRecoMatchedAll, mc, skipParticlesAbove);
+      }
       
       // (RECO-matched (quantities from MC particle) all particles)
       // STEP 5
-      fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTracked, tracksRecoMatchedAll, 0, weight);
+      fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTracked, tracksRecoMatchedAll, tracksCorrelateRecoMatchedAll, weight);
       
-      // Get RECO tracks
-      TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
+      // mixed event
+      if (fFillMixed)
+      {
+       AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx + 300);
+       if (pool->IsReady())
+         for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) 
+           fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTracked, tracksRecoMatchedAll, pool->GetEvent(jMix), 1.0 / pool->GetCurrentNEvents(), (jMix == 0));
+       pool->UpdatePool(CloneAndReduceTrackList(tracksCorrelateRecoMatchedAll));
+      }
       
+      // Get RECO tracks
+      // triggers
+      TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, fParticleSpeciesTrigger, kTRUE);
+      CleanUp(tracks, mc, skipParticlesAbove);
+      // associated
+      TObjArray* tracksCorrelate = tracks;
+      if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
+      {
+       tracksCorrelate = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, fParticleSpeciesAssociated, kTRUE);
+       CleanUp(tracksCorrelate, mc, skipParticlesAbove);
+      }
+     
       // (RECO all tracks)
       // STEP 6
       if (!fSkipStep6)
-       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, 0, weight);
+       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, tracksCorrelate, weight);
       
       // two track cut, STEP 8
       if (fTwoTrackEfficiencyCut > 0)
-       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracks, 0, weight, kTRUE, kTRUE, bSign, fTwoTrackEfficiencyCut);
+       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracks, tracksCorrelate, weight, kTRUE, kTRUE, bSign, fTwoTrackEfficiencyCut);
 
+      // apply correction efficiency, STEP 10
+      if (fEfficiencyCorrectionTriggers || fEfficiencyCorrectionAssociated)
+      {
+         // with or without two track efficiency depending on if fTwoTrackEfficiencyCut is set
+       Bool_t twoTrackCut = (fTwoTrackEfficiencyCut > 0);
+       
+       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepCorrected, tracks, tracksCorrelate, weight, kTRUE, twoTrackCut, bSign, fTwoTrackEfficiencyCut, kTRUE);
+      }
+      
       // mixed event
-      if (fFillMixed && !fSkipStep6)
+      if (fFillMixed)
       {
        AliEventPool* pool2 = fPoolMgr->GetEventPool(centrality, zVtx + 100);
-       //pool2->PrintInfo();
-       if (pool2->IsReady() || pool2->NTracksInPool() > fMixingTracks / 10 || pool2->GetCurrentNEvents() >= 5) 
+       if (pool2->IsReady())
+       {
          for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
          {
            // STEP 6
-           fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0));
+           if (!fSkipStep6)
+             fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0));
            
            // two track cut, STEP 8
            if (fTwoTrackEfficiencyCut > 0)
              fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0), kTRUE, bSign, fTwoTrackEfficiencyCut);
+           
+           // apply correction efficiency, STEP 10
+           if (fEfficiencyCorrectionTriggers || fEfficiencyCorrectionAssociated)
+           {
+             // with or without two track efficiency depending on if fTwoTrackEfficiencyCut is set
+             Bool_t twoTrackCut = (fTwoTrackEfficiencyCut > 0);
+             
+             fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepCorrected, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0), twoTrackCut, bSign, fTwoTrackEfficiencyCut, kTRUE);
+           }
          }
-       pool2->UpdatePool(CloneAndReduceTrackList(tracks));
+       }
+       pool2->UpdatePool(CloneAndReduceTrackList(tracksCorrelate));
       }
-
+      
       if (0 && !fReduceMemoryFootprint)
       {
         // make list of secondaries (matched with MC)
@@ -559,12 +782,22 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
         delete tracksRecoMatchedSecondaries;
       }
   
+      if (tracksCorrelateRecoMatchedPrim != tracksRecoMatchedPrim)
+       delete tracksCorrelateRecoMatchedPrim;
       delete tracksRecoMatchedPrim;
+
+      if (tracksCorrelateRecoMatchedAll != tracksRecoMatchedAll)
+       delete tracksCorrelateRecoMatchedAll;
       delete tracksRecoMatchedAll;
+      
+      if (tracksCorrelate != tracks)
+       delete tracksCorrelate;
       delete tracks;
     }
   }
   
+  if (tracksMC != tracksCorrelateMC)
+    delete tracksCorrelateMC;
   delete tracksMC;
 }
 
@@ -578,43 +811,90 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
 
   ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(0);
 
-  // skip not selected events here (the AOD is not updated for those)
   if (!fInputHandler)
     return;
     
-  if (!(fInputHandler->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kUserDefined)))
+  // skip not selected events here (the AOD is not updated for those)
+  if (!fSkipTrigger && !(fInputHandler->IsEventSelected() & fSelectBit))
+    return;
+  
+  // skip fast cluster events here if requested
+  if (fSkipFastCluster && (fInputHandler->IsEventSelected() & AliVEvent::kFastOnly))
     return;
 
+  // Support for ESD and AOD based analysis
+  AliVEvent* inputEvent = fAOD;
+  if (!inputEvent)
+    inputEvent = fESD;
+
   Double_t centrality = 0;
   
   AliCentrality *centralityObj = 0;
   if (fCentralityMethod.Length() > 0)
   {
-    if (fAOD)
-      centralityObj = fAOD->GetHeader()->GetCentralityP();
-    else if (fESD)
-      centralityObj = fESD->GetCentrality();
-    
-    if (centralityObj)
-      centrality = centralityObj->GetCentralityPercentile(fCentralityMethod);
-      //centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
-    else
-      centrality = -1;
+    if (fCentralityMethod == "ZNA_MANUAL")
+    {
+      Bool_t zna = kFALSE;
+      for(Int_t j = 0; j < 4; ++j) {
+       if (fESD->GetZDCData()->GetZDCTDCData(12,j) != 0) {
+         zna = kTRUE;
+       }
+      }
 
-    if (fAOD)
+//       Printf("%d %f", zna, fZNAtower[0]);
+      if (zna)
+      {
+       // code from Chiara O (23.10.12)
+       const Double_t *fZNAtower = fESD->GetZDCData()->GetZN2TowerEnergy();
+       Float_t znacut[4] = {681., 563., 413., 191.};
+       
+       if(fZNAtower[0]>znacut[0]) centrality = 1;
+       else if(fZNAtower[0]>znacut[1]) centrality = 21;
+       else if(fZNAtower[0]>znacut[2]) centrality = 41;
+       else if(fZNAtower[0]>znacut[3]) centrality = 61;
+       else centrality = 81;
+      }
+      else
+       centrality = -1;
+    }
+    else if (fCentralityMethod == "TRACKS_MANUAL")
     {
-      // remove outliers
-      if (centrality == 0)
+      // for pp
+      TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
+      centrality = tracks->GetEntriesFast();
+      if (centrality > 40)
+       centrality = 41;
+//       Printf("%d %f", tracks->GetEntriesFast(), centrality);
+      delete tracks;
+    }
+    else
+    {
+      if (fAOD)
+       centralityObj = fAOD->GetHeader()->GetCentralityP();
+      else if (fESD)
+       centralityObj = fESD->GetCentrality();
+      
+      if (centralityObj)
+       centrality = centralityObj->GetCentralityPercentile(fCentralityMethod);
+       //centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
+      else
+       centrality = -1;
+
+      if (fAOD)
       {
-       if (fAOD->GetVZEROData())
+       // remove outliers
+       if (centrality == 0)
        {
-         Float_t multV0 = 0;
-         for (Int_t i=0; i<64; i++)
-           multV0 += fAOD->GetVZEROData()->GetMultiplicity(i);
-         if (multV0 < 19500)
+         if (fAOD->GetVZEROData())
          {
-           centrality = -1;
-           AliInfo("Rejecting event due to too small V0 multiplicity");
+           Float_t multV0 = 0;
+           for (Int_t i=0; i<64; i++)
+             multV0 += fAOD->GetVZEROData()->GetMultiplicity(i);
+           if (multV0 < 19500)
+           {
+             centrality = -1;
+             AliInfo("Rejecting event due to too small V0 multiplicity");
+           }
          }
        }
       }
@@ -623,11 +903,6 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
     AliInfo(Form("Centrality is %f", centrality));
   }
   
-  // Support for ESD and AOD based analysis
-  AliVEvent* inputEvent = fAOD;
-  if (!inputEvent)
-    inputEvent = fESD;
-
   Float_t bSign = (inputEvent->GetMagneticField() > 0) ? 1 : -1;
 
   fHistos->SetRunNumber(inputEvent->GetRunNumber());
@@ -636,7 +911,7 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
   fHistos->FillEvent(centrality, AliUEHist::kCFStepAll);
   
   // Trigger selection ************************************************
-  if (!fAnalyseUE->TriggerSelection(fInputHandler)) return;
+  if (!fSkipTrigger && !fAnalyseUE->TriggerSelection(fInputHandler)) return;
   
   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
   fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
@@ -651,13 +926,50 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
   if (centrality < 0 && !fCompareCentralities)
     return;
 
-  TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
-  // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
-  TObjArray* tracksClone = CloneAndReduceTrackList(tracks);
-  delete tracks;
-  
+  TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, fParticleSpeciesTrigger, kTRUE);
   //Printf("Accepted %d tracks", tracks->GetEntries());
   
+  // check for outlier in centrality vs number of tracks (rough constants extracted from correlation histgram)
+  Bool_t reject = kFALSE;
+  if (fRejectCentralityOutliers)
+  {
+    if (centrality > 40 && centrality <= 50 && tracks->GetEntriesFast() > 1160)
+      reject = kTRUE;
+    if (centrality > 50 && centrality <= 60 && tracks->GetEntriesFast() > 650)
+      reject = kTRUE;
+    if (centrality > 60 && centrality <= 70 && tracks->GetEntriesFast() > 370)
+      reject = kTRUE;
+    if (centrality > 70 && centrality <= 80 && tracks->GetEntriesFast() > 220)
+      reject = kTRUE;
+    if (centrality > 80 && centrality <= 90 && tracks->GetEntriesFast() > 130)
+      reject = kTRUE;
+    if (centrality > 90 && tracks->GetEntriesFast() > 75)
+      reject = kTRUE;
+  }
+  
+  if (reject)
+  {
+    AliInfo(Form("Rejecting event due to centrality vs tracks correlation: %f %d", centrality, tracks->GetEntriesFast()));
+    fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
+    delete tracks;
+    return;
+  }
+  
+  // correlate particles with...
+  TObjArray* tracksCorrelate = 0;
+  if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
+    tracksCorrelate = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, fParticleSpeciesAssociated, kTRUE);
+  
+  // reference multiplicity
+  Int_t referenceMultiplicity = -1;
+  if (fESD)
+    referenceMultiplicity = AliESDtrackCuts::GetReferenceMultiplicity(fESD);
+  else if (fAOD)
+    referenceMultiplicity = tracks->GetEntriesFast(); // TODO to be replaced by the estimator once available in the AOD
+//    referenceMultiplicity = fAOD->GetHeader()->GetRefMultiplicityComb05();
+
+  ((TH2F*) fListOfHistos->FindObject("referenceMultiplicity"))->Fill(centrality, referenceMultiplicity);
+  
   const AliVVertex* vertex = inputEvent->GetPrimaryVertex();
   Double_t zVtx = vertex->GetZ();
   
@@ -674,12 +986,12 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
   if (centrality >= 0)
   {
     if (!fSkipStep6)
-      fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracksClone, 0, weight);
+      fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, tracksCorrelate, weight, kTRUE, kFALSE, 0, 0.02, kTRUE);
     
     ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(1);
     
     if (fTwoTrackEfficiencyCut > 0)
-      fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksClone, 0, weight, kTRUE, kTRUE, bSign, fTwoTrackEfficiencyCut);
+      fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracks, tracksCorrelate, weight, kTRUE, kTRUE, bSign, fTwoTrackEfficiencyCut, kTRUE);
   }
 
   // fill second time with SPD centrality
@@ -687,9 +999,13 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
   {
     centrality = centralityObj->GetCentralityPercentile("CL1");
     if (centrality >= 0 && !fSkipStep6)
-      fHistos->FillCorrelations(centrality, 2, AliUEHist::kCFStepReconstructed, tracksClone, 0, weight);
+      fHistos->FillCorrelations(centrality, 2, AliUEHist::kCFStepReconstructed, tracks, tracksCorrelate, weight, kFALSE, kFALSE, 0, 0.02, kTRUE);
   }
   
+  // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
+  TObjArray* tracksClone = CloneAndReduceTrackList(tracks);
+  delete tracks;
+  
   if (fFillMixed)
   {
     // event mixing
@@ -715,9 +1031,9 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
     if (!pool)
       AliFatal(Form("No pool found for centrality = %f, zVtx = %f", centrality, zVtx));
     
-    //pool->SetDebug(1);
+//     pool->SetDebug(1);
      
-    if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5
+    if (pool->IsReady()) 
     {
       
       Int_t nMix = pool->GetCurrentNEvents();
@@ -734,19 +1050,29 @@ void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
        TObjArray* bgTracks = pool->GetEvent(jMix);
        
        if (!fSkipStep6)
-         fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0));
+         fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0), kFALSE, 0, 0.02, kTRUE);
 
        if (fTwoTrackEfficiencyCut > 0)
-         fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0), kTRUE, bSign, fTwoTrackEfficiencyCut);
+         fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0), kTRUE, bSign, fTwoTrackEfficiencyCut, kTRUE);
       }
     }
     
     // ownership is with the pool now
-    pool->UpdatePool(tracksClone);
+    if (tracksCorrelate)
+    {
+      pool->UpdatePool(CloneAndReduceTrackList(tracksCorrelate));
+      delete tracksClone;
+    }
+    else
+      pool->UpdatePool(tracksClone);
     //pool->PrintInfo();
   }
   else
+  {
     delete tracksClone;
+    if (tracksCorrelate)
+      delete tracksCorrelate;
+  }
 }
 
 TObjArray* AliAnalysisTaskPhiCorrelations::CloneAndReduceTrackList(TObjArray* tracks)
@@ -774,3 +1100,47 @@ void  AliAnalysisTaskPhiCorrelations::Initialize()
   // MC handler
   fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 }
+
+//____________________________________________________________________
+void AliAnalysisTaskPhiCorrelations::RemoveDuplicates(TObjArray* tracks)
+{
+  // remove particles with the same label
+  
+  Int_t before = tracks->GetEntriesFast();
+
+  for (Int_t i=0; i<before; ++i) 
+  {
+    AliVParticle* part = (AliVParticle*) tracks->At(i);
+    
+    for (Int_t j=i+1; j<before; ++j) 
+    {
+      AliVParticle* part2 = (AliVParticle*) tracks->At(j);
+      
+      if (part->GetLabel() == part2->GetLabel())
+      {
+       Printf("Removing %d with label %d (duplicated in %d)", i, part->GetLabel(), j); part->Dump(); part2->Dump();
+       TObject* object = tracks->RemoveAt(i);
+       if (tracks->IsOwner())
+         delete object;
+       break;
+      }
+    }
+  }
+  tracks->Compress();
+  
+  if (before > tracks->GetEntriesFast())
+    AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast())); 
+}
+
+void AliAnalysisTaskPhiCorrelations::CleanUp(TObjArray* tracks, TObject* mcObj, Int_t maxLabel)
+{
+  // calls RemoveInjectedSignals, RemoveWeakDecays and RemoveDuplicates
+  
+  if (fInjectedSignals)
+    fAnalyseUE->RemoveInjectedSignals(tracks, mcObj, maxLabel);
+  if (fRemoveWeakDecays)
+    fAnalyseUE->RemoveWeakDecays(tracks, mcObj);
+  if (fRemoveDuplicates)
+    RemoveDuplicates(tracks);
+}