end-of-line normalization
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskTriggeredBF.cxx
index 261891c..603d085 100755 (executable)
-#include <vector>\r
-#include "TChain.h"\r
-#include "TList.h"\r
-#include "TCanvas.h"\r
-#include "TLorentzVector.h"\r
-#include "TGraphErrors.h"\r
-#include "TH1F.h"\r
-#include "TH2F.h"\r
-#include "TArrayF.h"\r
-#include "TF1.h"\r
-#include "TRandom.h"\r
-\r
-#include "AliLog.h"\r
-\r
-#include "AliAnalysisTaskSE.h"\r
-#include "AliAnalysisManager.h"\r
-\r
-#include "AliESDVertex.h"\r
-#include "AliESDEvent.h"\r
-#include "AliESDInputHandler.h"\r
-#include "AliAODEvent.h"\r
-#include "AliAODTrack.h"\r
-#include "AliAODInputHandler.h"\r
-#include "AliGenEventHeader.h"\r
-#include "AliGenHijingEventHeader.h"\r
-#include "AliMCEventHandler.h"\r
-#include "AliMCEvent.h"\r
-#include "AliMixInputEventHandler.h"\r
-#include "AliStack.h"\r
-\r
-#include "TH2D.h"    \r
-#include "AliTHn.h"             \r
-\r
-#include "AliEventPoolManager.h" \r
-\r
-#include "AliAnalysisTaskTriggeredBF.h"\r
-#include "AliBalanceTriggered.h"\r
-\r
-\r
-// Analysis task for the TriggeredBF code\r
-// Authors: Panos.Christakoglou@nikhef.nl, m.weber@cern.ch\r
-\r
-// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+\r
-//\r
-// For the V0 part:\r
-// --> AliAnalysisTaskExtractV0AOD (by david.chinellato@gmail.com)\r
-//\r
-// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+\r
-\r
-using std::cout;\r
-using std::endl;\r
-using std::vector;\r
-\r
-ClassImp(AliAnalysisTaskTriggeredBF)\r
-\r
-//________________________________________________________________________\r
-AliAnalysisTaskTriggeredBF::AliAnalysisTaskTriggeredBF(const char *name) \r
-: AliAnalysisTaskSE(name), \r
-  fBalance(0),\r
-  fRunShuffling(kFALSE),\r
-  fShuffledBalance(0),\r
-  fRunMixing(kFALSE),\r
-  fMixingTracks(50000),\r
-  fMixedBalance(0),\r
-  fPoolMgr(0),\r
-  fRunV0(kFALSE),\r
-  fPIDResponse(0),\r
-  fPIDCombined(0),\r
-  fList(0),\r
-  fListTriggeredBF(0),\r
-  fListTriggeredBFS(0),\r
-  fListTriggeredBFM(0),\r
-  fHistListPIDQA(0),\r
-  fHistListV0(0),\r
-  fHistEventStats(0),\r
-  fHistCentStats(0),\r
-  fHistTriggerStats(0),\r
-  fHistTrackStats(0),\r
-  fHistVx(0),\r
-  fHistVy(0),\r
-  fHistVz(0),\r
-  fHistClus(0),\r
-  fHistDCA(0),\r
-  fHistChi2(0),\r
-  fHistPt(0),\r
-  fHistEta(0),\r
-  fHistPhi(0),\r
-  fHistPhiBefore(0),\r
-  fHistPhiAfter(0),\r
-  fHistV0M(0),\r
-  fHistRefTracks(0),\r
-  fHistV0MultiplicityBeforeTrigSel(0),\r
-  fHistV0MultiplicityForTrigEvt(0),\r
-  fHistV0MultiplicityForSelEvt(0),\r
-  fHistV0MultiplicityForSelEvtNoTPCOnly(0),\r
-  fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),\r
-  fHistMultiplicityBeforeTrigSel(0),\r
-  fHistMultiplicityForTrigEvt(0),\r
-  fHistMultiplicity(0),\r
-  fHistMultiplicityNoTPCOnly(0),\r
-  fHistMultiplicityNoTPCOnlyNoPileup(0),\r
-  fHistV0InvMassK0(0),\r
-  fHistV0InvMassLambda(0),\r
-  fHistV0InvMassAntiLambda(0),\r
-  fHistV0Armenteros(0),\r
-  fHistV0SelInvMassK0(0),\r
-  fHistV0SelInvMassLambda(0),\r
-  fHistV0SelInvMassAntiLambda(0),\r
-  fHistV0SelArmenteros(0),\r
-  fCentralityEstimator("V0M"),\r
-  fUseCentrality(kFALSE),\r
-  fCentralityPercentileMin(0.), \r
-  fCentralityPercentileMax(5.),\r
-  fImpactParameterMin(0.),\r
-  fImpactParameterMax(20.),\r
-  fUseMultiplicity(kFALSE),\r
-  fNumberOfAcceptedTracksMin(0),\r
-  fNumberOfAcceptedTracksMax(10000),\r
-  fHistNumberOfAcceptedTracks(0),\r
-  fUseOfflineTrigger(kFALSE),\r
-  fVxMax(0.3),\r
-  fVyMax(0.3),\r
-  fVzMax(10.),\r
-  nAODtrackCutBit(128),\r
-  fPtMin(0.3),\r
-  fPtMax(1.5),\r
-  fEtaMin(-0.8),\r
-  fEtaMax(-0.8),\r
-  fDCAxyCut(-1),\r
-  fDCAzCut(-1),\r
-  fTPCchi2Cut(-1),\r
-  fNClustersTPCCut(-1)\r
-{\r
-  // Constructor\r
-  // Define input and output slots here\r
-  // Input slot #0 works with a TChain\r
-  DefineInput(0, TChain::Class());\r
-  // Output slot #0 writes into a TH1 container\r
-  DefineOutput(1, TList::Class());\r
-  DefineOutput(2, TList::Class());\r
-  DefineOutput(3, TList::Class());\r
-  DefineOutput(4, TList::Class());\r
-  DefineOutput(5, TList::Class());\r
-}\r
-\r
-//________________________________________________________________________\r
-AliAnalysisTaskTriggeredBF::~AliAnalysisTaskTriggeredBF() {\r
-\r
-  // Destructor\r
-\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskTriggeredBF::UserCreateOutputObjects() {\r
-  // Create histograms\r
-  // Called once\r
-\r
-  // global switch disabling the reference \r
-  // (to avoid "Replacing existing TH1" if several wagons are created in train)\r
-  Bool_t oldStatus = TH1::AddDirectoryStatus();\r
-  TH1::AddDirectory(kFALSE);\r
-\r
-  if(!fBalance) {\r
-    fBalance = new AliBalanceTriggered();\r
-    fBalance->SetAnalysisLevel("AOD");\r
-  }\r
-  if(fRunShuffling) {\r
-    if(!fShuffledBalance) {\r
-      fShuffledBalance = new AliBalanceTriggered();\r
-      fShuffledBalance->SetAnalysisLevel("AOD");\r
-    }\r
-  }\r
-  if(fRunMixing) {\r
-    if(!fMixedBalance) {\r
-      fMixedBalance = new AliBalanceTriggered();\r
-      fMixedBalance->SetAnalysisLevel("AOD");\r
-    }\r
-  }\r
-\r
-  //QA list\r
-  fList = new TList();\r
-  fList->SetName("listQA");\r
-  fList->SetOwner();\r
-\r
-  //Balance Function list\r
-  fListTriggeredBF = new TList();\r
-  fListTriggeredBF->SetName("listTriggeredBF");\r
-  fListTriggeredBF->SetOwner();\r
-\r
-  if(fRunShuffling) {\r
-    fListTriggeredBFS = new TList();\r
-    fListTriggeredBFS->SetName("listTriggeredBFShuffled");\r
-    fListTriggeredBFS->SetOwner();\r
-  }\r
-  if(fRunMixing) {\r
-    fListTriggeredBFM = new TList();\r
-    fListTriggeredBFM->SetName("listTriggeredBFMixed");\r
-    fListTriggeredBFM->SetOwner();\r
-  }\r
-  \r
-  \r
-  //Event stats.\r
-  TString gCutName[4] = {"Total","Offline trigger",\r
-                         "Vertex","Analyzed"};\r
-  fHistEventStats = new TH1F("fHistEventStats",\r
-                             "Event statistics;;N_{events}",\r
-                             4,0.5,4.5);\r
-  for(Int_t i = 1; i <= 4; i++)\r
-    fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
-  fList->Add(fHistEventStats);\r
-  \r
-  TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r
-  fHistCentStats = new TH2F("fHistCentStats",\r
-                           "Centrality statistics;;Cent percentile",\r
-                           9,-0.5,8.5,220,-5,105);\r
-  for(Int_t i = 1; i <= 9; i++)\r
-    fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
-  fList->Add(fHistCentStats);\r
-  \r
-  fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);\r
-  fList->Add(fHistTriggerStats);\r
-  \r
-  fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);\r
-  fList->Add(fHistTrackStats);\r
-\r
-  fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);\r
-  fList->Add(fHistNumberOfAcceptedTracks);\r
-\r
-  // Vertex distributions\r
-  fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);\r
-  fList->Add(fHistVx);\r
-  fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);\r
-  fList->Add(fHistVy);\r
-  fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);\r
-  fList->Add(fHistVz);\r
-\r
-  // QA histograms\r
-  fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);\r
-  fList->Add(fHistClus);\r
-  fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);\r
-  fList->Add(fHistChi2);\r
-  fHistDCA  = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); \r
-  fList->Add(fHistDCA);\r
-  fHistPt   = new TH1F("fHistPt","p_{T} distribution",200,0,10);\r
-  fList->Add(fHistPt);\r
-  fHistEta  = new TH1F("fHistEta","#eta distribution",200,-2,2);\r
-  fList->Add(fHistEta);\r
-  fHistPhi  = new TH1F("fHistPhi","#phi distribution",200,-20,380);\r
-  fList->Add(fHistPhi);\r
-  fHistPhiBefore  = new TH1F("fHistPhiBefore","#phi distribution",200,0.,2*TMath::Pi());\r
-  fList->Add(fHistPhiBefore);\r
-  fHistPhiAfter  = new TH1F("fHistPhiAfter","#phi distribution",200,0.,2*TMath::Pi());\r
-  fList->Add(fHistPhiAfter);\r
-  fHistV0M  = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);\r
-  fList->Add(fHistV0M);\r
-  TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};\r
-  fHistRefTracks  = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);\r
-  for(Int_t i = 1; i <= 6; i++)\r
-    fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
-  fList->Add(fHistRefTracks);\r
-\r
-  //------------------------------------------------\r
-  // V0 Multiplicity Histograms\r
-  //------------------------------------------------\r
-  if(fRunV0){\r
-    fHistListV0 = new TList();\r
-    fHistListV0->SetOwner();  // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner\r
-    \r
-    if(! fHistV0MultiplicityBeforeTrigSel) {\r
-      fHistV0MultiplicityBeforeTrigSel = new TH1F("fHistV0MultiplicityBeforeTrigSel", \r
-                                                 "V0s per event (before Trig. Sel.);Nbr of V0s/Evt;Events", \r
-                                                 25, 0, 25);\r
-      fHistListV0->Add(fHistV0MultiplicityBeforeTrigSel);\r
-    }\r
-    \r
-    if(! fHistV0MultiplicityForTrigEvt) {\r
-      fHistV0MultiplicityForTrigEvt = new TH1F("fHistV0MultiplicityForTrigEvt", \r
-                                              "V0s per event (for triggered evt);Nbr of V0s/Evt;Events", \r
-                                              25, 0, 25);\r
-      fHistListV0->Add(fHistV0MultiplicityForTrigEvt);\r
-    }\r
-    \r
-    if(! fHistV0MultiplicityForSelEvt) {\r
-      fHistV0MultiplicityForSelEvt = new TH1F("fHistV0MultiplicityForSelEvt", \r
-                                             "V0s per event;Nbr of V0s/Evt;Events", \r
-                                             25, 0, 25);\r
-      fHistListV0->Add(fHistV0MultiplicityForSelEvt);\r
-    }\r
-    \r
-    if(! fHistV0MultiplicityForSelEvtNoTPCOnly) {\r
-    fHistV0MultiplicityForSelEvtNoTPCOnly = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnly", \r
-                                                    "V0s per event;Nbr of V0s/Evt;Events", \r
-                                                    25, 0, 25);\r
-    fHistListV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnly);\r
-    }\r
-    \r
-    if(! fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup) {\r
-      fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup", \r
-                                                            "V0s per event;Nbr of V0s/Evt;Events", \r
-                                                              25, 0, 25);\r
-      fHistListV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup);\r
-    }\r
-    \r
-    //------------------------------------------------\r
-    // Track Multiplicity Histograms\r
-    //------------------------------------------------\r
-    \r
-    if(! fHistMultiplicityBeforeTrigSel) {\r
-      fHistMultiplicityBeforeTrigSel = new TH1F("fHistMultiplicityBeforeTrigSel", \r
-                                             "Tracks per event;Nbr of Tracks;Events", \r
-                                               200, 0, 200);           \r
-      fHistListV0->Add(fHistMultiplicityBeforeTrigSel);\r
-    }\r
-    if(! fHistMultiplicityForTrigEvt) {\r
-      fHistMultiplicityForTrigEvt = new TH1F("fHistMultiplicityForTrigEvt", \r
-                                            "Tracks per event;Nbr of Tracks;Events", \r
-                                            200, 0, 200);              \r
-      fHistListV0->Add(fHistMultiplicityForTrigEvt);\r
-    }\r
-    if(! fHistMultiplicity) {\r
-      fHistMultiplicity = new TH1F("fHistMultiplicity", \r
-                                  "Tracks per event;Nbr of Tracks;Events", \r
-                                  200, 0, 200);                \r
-      fHistListV0->Add(fHistMultiplicity);\r
-    }\r
-    if(! fHistMultiplicityNoTPCOnly) {\r
-      fHistMultiplicityNoTPCOnly = new TH1F("fHistMultiplicityNoTPCOnly", \r
-                                           "Tracks per event;Nbr of Tracks;Events", \r
-                                           200, 0, 200);               \r
-    fHistListV0->Add(fHistMultiplicityNoTPCOnly);\r
-    }\r
-    if(! fHistMultiplicityNoTPCOnlyNoPileup) {\r
-      fHistMultiplicityNoTPCOnlyNoPileup = new TH1F("fHistMultiplicityNoTPCOnlyNoPileup", \r
-                                                   "Tracks per event;Nbr of Tracks;Events", \r
-                                                   200, 0, 200);               \r
-      fHistListV0->Add(fHistMultiplicityNoTPCOnlyNoPileup);\r
-    }\r
-\r
-    //------------------------------------------------\r
-    // V0 selection Histograms (before)\r
-    //------------------------------------------------\r
-    if(!fHistV0InvMassK0) {\r
-      fHistV0InvMassK0 = new TH1F("fHistV0InvMassK0",\r
-                                 "Invariant Mass for K0;Mass (GeV/c^{2});Events",\r
-                                 200,0,2);\r
-      fHistListV0->Add(fHistV0InvMassK0);\r
-    }\r
-    if(!fHistV0InvMassLambda) {\r
-      fHistV0InvMassLambda = new TH1F("fHistV0InvMassLambda",\r
-                                 "Invariant Mass for Lambda;Mass (GeV/c^{2});Events",\r
-                                 200,0,2);\r
-      fHistListV0->Add(fHistV0InvMassLambda);\r
-    }\r
-    if(!fHistV0InvMassAntiLambda) {\r
-      fHistV0InvMassAntiLambda = new TH1F("fHistV0InvMassAntiLambda",\r
-                                 "Invariant Mass for AntiLambda;Mass (GeV/c^{2});Events",\r
-                                 200,0,2);\r
-      fHistListV0->Add(fHistV0InvMassAntiLambda);\r
-    }\r
-    if(!fHistV0Armenteros) {\r
-      fHistV0Armenteros = new TH2F("fHistV0Armenteros",\r
-                                 "Armenteros plot;#alpha;q_{t}",\r
-                                  200,-1,1,200,0,0.5);\r
-      fHistListV0->Add(fHistV0Armenteros);\r
-    }\r
-    \r
-    //------------------------------------------------\r
-    // V0 selection Histograms (after)\r
-    //------------------------------------------------\r
-    if(!fHistV0SelInvMassK0) {\r
-      fHistV0SelInvMassK0 = new TH1F("fHistV0SelInvMassK0",\r
-                                 "Invariant Mass for K0;Mass (GeV/c^{2});Events",\r
-                                 200,0,2);\r
-      fHistListV0->Add(fHistV0SelInvMassK0);\r
-    }\r
-    if(!fHistV0SelInvMassLambda) {\r
-      fHistV0SelInvMassLambda = new TH1F("fHistV0SelInvMassLambda",\r
-                                 "Invariant Mass for Lambda;Mass (GeV/c^{2});Events",\r
-                                 200,0,2);\r
-      fHistListV0->Add(fHistV0SelInvMassLambda);\r
-    }\r
-    if(!fHistV0SelInvMassAntiLambda) {\r
-      fHistV0SelInvMassAntiLambda = new TH1F("fHistV0SelInvMassAntiLambda",\r
-                                 "Invariant Mass for AntiLambda;Mass (GeV/c^{2});Events",\r
-                                 200,0,2);\r
-      fHistListV0->Add(fHistV0SelInvMassAntiLambda);\r
-    }\r
-    if(!fHistV0SelArmenteros) {\r
-      fHistV0SelArmenteros = new TH2F("fHistV0SelArmenteros",\r
-                                 "Armenteros plot;#alpha;q_{t}",\r
-                                  200,-1,1,200,0,0.5);\r
-      fHistListV0->Add(fHistV0SelArmenteros);\r
-    }\r
-  }//V0\r
-    \r
-  // Balance function histograms\r
-  // Initialize histograms if not done yet\r
-  if(!fBalance->GetHistNp()){\r
-    AliWarning("Histograms not yet initialized! --> Will be done now");\r
-    AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
-    fBalance->InitHistograms();\r
-  }\r
-\r
-  if(fRunShuffling) {\r
-    if(!fShuffledBalance->GetHistNp()) {\r
-      AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
-      AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
-      fShuffledBalance->InitHistograms();\r
-    }\r
-  }\r
-\r
-  if(fRunMixing) {\r
-    if(!fMixedBalance->GetHistNp()) {\r
-      AliWarning("Histograms (mixing) not yet initialized! --> Will be done now");\r
-      AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
-      fMixedBalance->InitHistograms();\r
-    }\r
-  }\r
-\r
-  fListTriggeredBF->Add(fBalance->GetHistNp());\r
-  fListTriggeredBF->Add(fBalance->GetHistNn());\r
-  fListTriggeredBF->Add(fBalance->GetHistNpn());\r
-  fListTriggeredBF->Add(fBalance->GetHistNnn());\r
-  fListTriggeredBF->Add(fBalance->GetHistNpp());\r
-  fListTriggeredBF->Add(fBalance->GetHistNnp());\r
-  \r
-  if(fRunShuffling) {\r
-    fListTriggeredBFS->Add(fShuffledBalance->GetHistNp());\r
-    fListTriggeredBFS->Add(fShuffledBalance->GetHistNn());\r
-    fListTriggeredBFS->Add(fShuffledBalance->GetHistNpn());\r
-    fListTriggeredBFS->Add(fShuffledBalance->GetHistNnn());\r
-    fListTriggeredBFS->Add(fShuffledBalance->GetHistNpp());\r
-    fListTriggeredBFS->Add(fShuffledBalance->GetHistNnp());\r
-  }  \r
-\r
-  if(fRunMixing) {\r
-    fListTriggeredBFM->Add(fMixedBalance->GetHistNp());\r
-    fListTriggeredBFM->Add(fMixedBalance->GetHistNn());\r
-    fListTriggeredBFM->Add(fMixedBalance->GetHistNpn());\r
-    fListTriggeredBFM->Add(fMixedBalance->GetHistNnn());\r
-    fListTriggeredBFM->Add(fMixedBalance->GetHistNpp());\r
-    fListTriggeredBFM->Add(fMixedBalance->GetHistNnp());\r
-  }  \r
-\r
-  // PID Response task active?\r
-  if(fRunV0) {\r
-    fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();\r
-    if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
-  }\r
-\r
-  // Event Mixing\r
-  Int_t trackDepth = fMixingTracks; \r
-  Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager\r
-   \r
-  Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,90.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
-  Double_t* centbins        = centralityBins;\r
-  Int_t nCentralityBins     = sizeof(centralityBins) / sizeof(Double_t) - 1;\r
-  \r
-  // bins for second buffer are shifted by 100 cm\r
-  Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
-  Double_t* vtxbins     = vertexBins;\r
-  Int_t nVertexBins     = sizeof(vertexBins) / sizeof(Double_t) - 1;\r
-\r
-  fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);\r
-\r
-\r
-  // Post output data.\r
-  PostData(1, fList);\r
-  PostData(2, fListTriggeredBF);\r
-  if(fRunShuffling) PostData(3, fListTriggeredBFS);\r
-  if(fRunMixing) PostData(4, fListTriggeredBFM);\r
-  if(fRunV0) PostData(5,fHistListV0);\r
-\r
-  TH1::AddDirectory(oldStatus);\r
-\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskTriggeredBF::UserExec(Option_t *) {\r
-  // Main loop\r
-  // Called for each event\r
-\r
-  TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
-  Float_t fCentrality = -1.;  \r
-\r
-  // -------------------------------------------------------------                  \r
-  // AOD analysis (vertex and track cuts also here!!!!)\r
-  if(gAnalysisLevel == "AOD") {\r
-    AliVEvent* eventMain = dynamic_cast<AliVEvent*>(InputEvent()); \r
-    if(!eventMain) {\r
-      AliError("eventMain not available");\r
-      return;\r
-    }\r
-\r
-    // check event cuts and fill event histograms\r
-    if((fCentrality = IsEventAccepted(eventMain)) < 0){\r
-      return;\r
-    }\r
-    \r
-    // get the accepted tracks in main event\r
-    TObjArray *tracksMain = NULL;\r
-    if(fRunV0) tracksMain = GetAcceptedV0s(eventMain);\r
-    else       tracksMain = GetAcceptedTracks(eventMain);\r
-\r
-    // store charges of all accepted tracks, shuffle and reassign (two extra loops!)\r
-    TObjArray* tracksShuffled = NULL;\r
-    if(fRunShuffling){\r
-      tracksShuffled = GetShuffledTracks(tracksMain);\r
-    }\r
-    \r
-    // Event mixing --> UPDATE POOL IS MISSING!!!\r
-    if (fRunMixing)\r
-      {\r
-        // 1. First get an event pool corresponding in mult (cent) and\r
-        //    zvertex to the current event. Once initialized, the pool\r
-        //    should contain nMix (reduced) events. This routine does not\r
-        //    pre-scan the chain. The first several events of every chain\r
-        //    will be skipped until the needed pools are filled to the\r
-        //    specified depth. If the pool categories are not too rare, this\r
-        //    should not be a problem. If they are rare, you could lose`\r
-       //    statistics.\r
-       \r
-       // 2. Collect the whole pool's content of tracks into one TObjArray\r
-       //    (bgTracks), which is effectively a single background super-event.\r
-       \r
-       // 3. The reduced and bgTracks arrays must both be passed into\r
-       //    FillCorrelations(). Also nMix should be passed in, so a weight\r
-       //    of 1./nMix can be applied.\r
-       \r
-       AliEventPool* pool = fPoolMgr->GetEventPool(fCentrality, eventMain->GetPrimaryVertex()->GetZ());\r
-       \r
-       if (!pool){\r
-         AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCentrality, eventMain->GetPrimaryVertex()->GetZ()));\r
-       }\r
-       else{\r
-\r
-       //pool->SetDebug(1);\r
-       \r
-         if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){ \r
-           \r
-           \r
-           Int_t nMix = pool->GetCurrentNEvents();\r
-           //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;\r
-           \r
-           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);\r
-           //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());\r
-           //if (pool->IsReady())\r
-           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);\r
-           \r
-           // Fill mixed-event histos here  \r
-           for (Int_t jMix=0; jMix<nMix; jMix++) \r
-             {\r
-               TObjArray* tracksMixed = pool->GetEvent(jMix);\r
-               fMixedBalance->FillBalance(fCentrality,tracksMain,tracksMixed); \r
-             }\r
-         }\r
-         \r
-         // Update the Event pool\r
-         pool->UpdatePool(tracksMain);\r
-         //pool->PrintInfo();\r
-         \r
-       }//pool NULL check  \r
-      }//run mixing\r
-    \r
-    // calculate balance function\r
-    fBalance->FillBalance(fCentrality,tracksMain,NULL);\r
-    \r
-    // calculate shuffled balance function\r
-    if(fRunShuffling && tracksShuffled != NULL) {\r
-      fShuffledBalance->FillBalance(fCentrality,tracksShuffled,NULL);\r
-    }\r
-    \r
-  }//AOD analysis\r
-  else{\r
-    AliError("Triggered Balance Function analysis only for AODs!");\r
-  }\r
-}     \r
-\r
-//________________________________________________________________________\r
-Float_t AliAnalysisTaskTriggeredBF::IsEventAccepted(AliVEvent *event){\r
-  // Checks the Event cuts\r
-  // Fills Event statistics histograms\r
-  \r
-  // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
-  fHistEventStats->Fill(1); //all events\r
-\r
-  Bool_t isSelectedMain = kTRUE;\r
-  Float_t fCentrality = -1.;\r
-  Int_t nV0s          = event->GetNumberOfV0s();\r
-  TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
-  \r
-  if(fUseOfflineTrigger)\r
-    isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
-  \r
-  //V0 QA histograms (before trigger selection)\r
-  if(fRunV0){\r
-    fHistMultiplicityBeforeTrigSel->Fill ( -1 );\r
-    fHistV0MultiplicityBeforeTrigSel->Fill ( nV0s );\r
-  }\r
-  \r
-  if(isSelectedMain) {\r
-    fHistEventStats->Fill(2); //triggered events\r
-    \r
-    //Centrality stuff \r
-    if(fUseCentrality) {\r
-      if(gAnalysisLevel == "AOD") { //centrality in AOD header\r
-       AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
-       fCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
-\r
-       // QA for centrality estimators\r
-       fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));\r
-       fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("FMD"));\r
-       fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("TRK"));\r
-       fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("TKL"));\r
-       fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("CL0"));\r
-       fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("CL1"));\r
-       fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
-       fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
-       fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
-       \r
-       // centrality QA (V0M)\r
-       fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
-       \r
-       // centrality QA (reference tracks)\r
-       fHistRefTracks->Fill(0.,header->GetRefMultiplicity());\r
-       fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());\r
-       fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());\r
-       fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());\r
-       fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));\r
-       fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));\r
-       fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));\r
-       fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));\r
-       fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));\r
-\r
-       //V0 QA histograms (after trigger selection)\r
-       if(fRunV0){\r
-         fHistMultiplicityForTrigEvt->Fill ( fCentrality );\r
-         fHistV0MultiplicityForTrigEvt->Fill ( nV0s );\r
-       }\r
-      }\r
-    }\r
-    \r
-    \r
-    const AliVVertex *vertex = event->GetPrimaryVertex();\r
-    \r
-    if(vertex) {\r
-      Double32_t fCov[6];\r
-      vertex->GetCovarianceMatrix(fCov);\r
-      if(vertex->GetNContributors() > 0) {\r
-       if(fCov[5] != 0) {\r
-         fHistEventStats->Fill(3); //events with a proper vertex\r
-         if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
-           if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
-             if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
-               fHistEventStats->Fill(4); //analyzed events\r
-               fHistVx->Fill(vertex->GetX());\r
-               fHistVy->Fill(vertex->GetY());\r
-               fHistVz->Fill(vertex->GetZ());\r
-\r
-\r
-\r
-               //V0 QA histograms (vertex Z check)\r
-               if(fRunV0){\r
-                 fHistV0MultiplicityForSelEvt ->Fill( nV0s );\r
-                 fHistMultiplicity->Fill(fCentrality);\r
-\r
-                 //V0 QA histograms (Only look at events with well-established PV)\r
-                 const AliAODVertex *lPrimarySPDVtx = ((AliAODEvent*)event)->GetPrimaryVertexSPD();\r
-                 if(lPrimarySPDVtx){\r
-                   fHistMultiplicityNoTPCOnly->Fill ( fCentrality );\r
-                   fHistV0MultiplicityForSelEvtNoTPCOnly->Fill ( nV0s );\r
-                   \r
-                   //V0 QA histograms (Pileup Rejection)\r
-                   // FIXME : quality selection regarding pile-up rejection \r
-                   fHistMultiplicityNoTPCOnlyNoPileup->Fill(fCentrality);\r
-                   fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup ->Fill( nV0s );\r
-\r
-                 }\r
-                 else{\r
-                   return -1;\r
-                 }\r
-               }\r
-               \r
-               \r
-               // take only events inside centrality class\r
-               if((fCentrality > fCentralityPercentileMin) && (fCentrality < fCentralityPercentileMax)){\r
-                 return fCentrality;           \r
-               }//centrality class\r
-             }//Vz cut\r
-           }//Vy cut\r
-         }//Vx cut\r
-       }//proper vertex resolution\r
-      }//proper number of contributors\r
-    }//vertex object valid\r
-  }//triggered event \r
-  \r
-  // in all other cases return -1 (event not accepted)\r
-  return -1;\r
-}\r
-\r
-//________________________________________________________________________\r
-TObjArray* AliAnalysisTaskTriggeredBF::GetAcceptedTracks(AliVEvent *event){\r
-  // Returns TObjArray with tracks after all track cuts (only for AOD!)\r
-  // Fills QA histograms\r
-\r
-  //output TObjArray holding all good tracks\r
-  TObjArray* tracksAccepted = new TObjArray;\r
-  tracksAccepted->SetOwner(kTRUE);\r
-\r
-  Short_t vCharge = 0;\r
-  Double_t vEta    = 0.;\r
-  Double_t vPhi    = 0.;\r
-  Double_t vPt     = 0.;\r
-  \r
-  // Loop over tracks in event\r
-  for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
-    AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
-    if (!aodTrack) {\r
-      AliError(Form("Could not receive track %d", iTracks));\r
-      continue;\r
-    }\r
-    \r
-    // AOD track cuts\r
-    \r
-    // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
-    // take only TPC only tracks \r
-    fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
-    if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
-    \r
-    vCharge = aodTrack->Charge();\r
-    vEta    = aodTrack->Eta();\r
-    vPhi    = aodTrack->Phi() * TMath::RadToDeg();\r
-    vPt     = aodTrack->Pt();\r
-    \r
-    Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
-    Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
-    \r
-    \r
-    // Kinematics cuts from ESD track cuts\r
-    if( vPt < fPtMin || vPt > fPtMax)      continue;\r
-    if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
-    \r
-    // Extra DCA cuts (for systematic studies [!= -1])\r
-    if( fDCAxyCut != -1 && fDCAzCut != -1){\r
-      if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
-       continue;  // 2D cut\r
-      }\r
-    }\r
-    \r
-    // Extra TPC cuts (for systematic studies [!= -1])\r
-    if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
-      continue;\r
-    }\r
-    if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
-      continue;\r
-    }\r
-    \r
-    // fill QA histograms\r
-    fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
-    fHistDCA->Fill(dcaZ,dcaXY);\r
-    fHistChi2->Fill(aodTrack->Chi2perNDF());\r
-    fHistPt->Fill(vPt);\r
-    fHistEta->Fill(vEta);\r
-    fHistPhi->Fill(vPhi);\r
-    \r
-    // add the track to the TObjArray\r
-    tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge,1.));\r
-  }\r
-\r
-  return tracksAccepted;\r
-}\r
-\r
-//________________________________________________________________________\r
-TObjArray* AliAnalysisTaskTriggeredBF::GetAcceptedV0s(AliVEvent *event){\r
-  // Returns TObjArray with tracks after all track cuts (only for AOD!)\r
-  // Fills QA histograms\r
-\r
-  //output TObjArray holding all good tracks\r
-  TObjArray* tracksAccepted = new TObjArray;\r
-  tracksAccepted->SetOwner(kTRUE);\r
-\r
-  Short_t vCharge = 0;\r
-  Double_t vEta    = 0.;\r
-  Double_t vPhi    = 0.;\r
-  Double_t vPt     = 0.;\r
-  \r
-  //------------------------------------------------\r
-  // MAIN LAMBDA LOOP STARTS HERE (basically a copy of AliAnalysisTaskExtractV0AOD)\r
-  //------------------------------------------------\r
-\r
-  // parameters (for the time being hard coded here) --> from David for EbyE Lambdas\r
-  Bool_t fkUseOnTheFly = kFALSE;\r
-  Double_t fRapidityBoundary  = 0.5; \r
-  Double_t fCutDaughterEta    = 0.8;\r
-  Double_t fCutV0Radius       = 0.9;\r
-  Double_t fCutDCANegToPV     = 0.1;\r
-  Double_t fCutDCAPosToPV     = 0.1;\r
-  Double_t fCutDCAV0Daughters = 1.0;\r
-  Double_t fCutV0CosPA        = 0.9995;\r
-  Double_t fMassLambda        = 1.115683;\r
-  Double_t fCutMassLambda     = 0.007;\r
-  Double_t fCutProperLifetime = 3*7.9;\r
-  Double_t fCutLeastNumberOfCrossedRows = 70;\r
-  Double_t fCutLeastNumberOfCrossedRowsOverFindable = 0.8;\r
-  Double_t fCutTPCPIDNSigmasProton  = 3.0;\r
-  Double_t fCutTPCPIDNSigmasPion    = 5.0;\r
-\r
-\r
-  //Variable definition\r
-  Int_t    lOnFlyStatus = 0;// nv0sOn = 0, nv0sOff = 0;\r
-  Double_t lChi2V0 = 0;\r
-  Double_t lDcaV0Daughters = 0, lDcaV0ToPrimVertex = 0;\r
-  Double_t lDcaPosToPrimVertex = 0, lDcaNegToPrimVertex = 0;\r
-  Double_t lV0CosineOfPointingAngle = 0;\r
-  Double_t lV0Radius = 0, lPt = 0;\r
-  Double_t lEta = 0, lPhi = 0;\r
-  Double_t lRap = 0, lRapK0Short = 0, lRapLambda = 0;\r
-  Double_t lInvMassK0s = 0, lInvMassLambda = 0, lInvMassAntiLambda = 0;\r
-  Double_t lAlphaV0 = 0, lPtArmV0 = 0;\r
-  \r
-  Double_t fMinV0Pt = 0; \r
-  Double_t fMaxV0Pt = 100; \r
-  \r
-\r
-  \r
-  // some event observables\r
-  Int_t nv0s = event->GetNumberOfV0s();\r
-  Double_t tPrimaryVtxPosition[3];\r
-  const AliVVertex *primaryVtx = event->GetPrimaryVertex();\r
-  tPrimaryVtxPosition[0] = primaryVtx->GetX();\r
-  tPrimaryVtxPosition[1] = primaryVtx->GetY();\r
-  tPrimaryVtxPosition[2] = primaryVtx->GetZ();\r
-\r
-\r
-  //loop over V0s  \r
-  for (Int_t iV0 = 0; iV0 < nv0s; iV0++) \r
-    {// This is the begining of the V0 loop\r
-      AliAODv0 *v0 = ((AliAODEvent*)event)->GetV0(iV0);\r
-      if (!v0) continue;\r
-\r
-      //Obsolete at AOD level... \r
-      //---> Fix On-the-Fly candidates, count how many swapped\r
-      //if( v0->GetParamN()->Charge() > 0 && v0->GetParamP()->Charge() < 0 ){\r
-      //  fHistSwappedV0Counter -> Fill( 1 );\r
-      //}else{\r
-      //  fHistSwappedV0Counter -> Fill( 0 ); \r
-      //}\r
-      //if ( fkUseOnTheFly ) CheckChargeV0(v0); \r
-      \r
-      Double_t tDecayVertexV0[3]; v0->GetXYZ(tDecayVertexV0); \r
-      Double_t tV0mom[3];\r
-      v0->GetPxPyPz( tV0mom ); \r
-      Double_t lV0TotalMomentum = TMath::Sqrt(\r
-                                             tV0mom[0]*tV0mom[0]+tV0mom[1]*tV0mom[1]+tV0mom[2]*tV0mom[2] );\r
-      \r
-      lV0Radius = TMath::Sqrt(tDecayVertexV0[0]*tDecayVertexV0[0]+tDecayVertexV0[1]*tDecayVertexV0[1]);\r
-      lPt = v0->Pt();\r
-      lEta = v0->Eta();\r
-      lPhi = v0->Phi()*TMath::RadToDeg();\r
-      lRapK0Short = v0->RapK0Short();\r
-      lRapLambda  = v0->RapLambda();\r
-      lRap        = lRapLambda;//v0->Y(); //FIXME!!!\r
-      if ((lPt<fMinV0Pt)||(fMaxV0Pt<lPt)) continue;\r
-      \r
-      //UInt_t lKeyPos = (UInt_t)TMath::Abs(v0->GetPosID());\r
-      //UInt_t lKeyNeg = (UInt_t)TMath::Abs(v0->GetPosID());\r
-\r
-      Double_t lMomPos[3]; //v0->GetPPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);\r
-      Double_t lMomNeg[3]; //v0->GetNPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);\r
-      lMomPos[0] = v0->MomPosX();\r
-      lMomPos[1] = v0->MomPosY();\r
-      lMomPos[2] = v0->MomPosZ();\r
-      lMomNeg[0] = v0->MomNegX();\r
-      lMomNeg[1] = v0->MomNegY();\r
-      lMomNeg[2] = v0->MomNegZ();\r
-      \r
-      AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter\r
-      AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter\r
-      if (!pTrack || !nTrack) {\r
-       AliError("ERROR: Could not retreive one of the daughter track");\r
-       continue;\r
-      }\r
-\r
-      //Daughter Eta for Eta selection, afterwards\r
-      Double_t lNegEta = nTrack->Eta();\r
-      Double_t lPosEta = pTrack->Eta();\r
-      \r
-      // Filter like-sign V0 (next: add counter and distribution)\r
-      if ( pTrack->Charge() == nTrack->Charge()){\r
-       continue;\r
-      } \r
-      \r
-      //Quick test this far! \r
-      \r
-\r
-      //________________________________________________________________________\r
-      // Track quality cuts \r
-      Float_t lPosTrackCrossedRows = pTrack->GetTPCClusterInfo(2,1);\r
-      Float_t lNegTrackCrossedRows = nTrack->GetTPCClusterInfo(2,1);\r
-      Float_t lLeastNbrCrossedRows =  (lPosTrackCrossedRows>lNegTrackCrossedRows) ? lNegTrackCrossedRows : lPosTrackCrossedRows;\r
-\r
-      // TPC refit condition (done during reconstruction for Offline but not for On-the-fly)\r
-      if( !(pTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;\r
-      if( !(nTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;\r
-      \r
-      if ( ( ( pTrack->GetTPCClusterInfo(2,1) ) < 70 ) || ( ( nTrack->GetTPCClusterInfo(2,1) ) < 70 ) ) continue;\r
-      \r
-      //Findable clusters > 0 condition\r
-      if( pTrack->GetTPCNclsF()<=0 || nTrack->GetTPCNclsF()<=0 ) continue;\r
-      \r
-      //Compute ratio Crossed Rows / Findable clusters\r
-      //Note: above test avoids division by zero! \r
-      Float_t lPosTrackCrossedRowsOverFindable = lPosTrackCrossedRows / ((double)(pTrack->GetTPCNclsF())); \r
-      Float_t lNegTrackCrossedRowsOverFindable = lNegTrackCrossedRows / ((double)(nTrack->GetTPCNclsF())); \r
-      Float_t lLeastNbrCrossedRowsOverFindable = (lPosTrackCrossedRowsOverFindable>lNegTrackCrossedRowsOverFindable) ? lNegTrackCrossedRowsOverFindable : lPosTrackCrossedRowsOverFindable;\r
-\r
-      //Lowest Cut Level for Ratio Crossed Rows / Findable = 0.8, set here\r
-      if ( lLeastNbrCrossedRowsOverFindable < 0.8) continue;\r
-      \r
-      //End track Quality Cuts\r
-      //________________________________________________________________________\r
-      \r
-      \r
-      lDcaPosToPrimVertex = v0->DcaPosToPrimVertex();\r
-      lDcaNegToPrimVertex = v0->DcaNegToPrimVertex();\r
-          \r
-      lOnFlyStatus = v0->GetOnFlyStatus();\r
-      lChi2V0 = v0->Chi2V0();\r
-      lDcaV0Daughters = v0->DcaV0Daughters();\r
-      lDcaV0ToPrimVertex = v0->DcaV0ToPrimVertex();\r
-      lV0CosineOfPointingAngle = v0->CosPointingAngle(tPrimaryVtxPosition);\r
-      \r
-      // Distance over total momentum\r
-      Double_t lDistOverTotMom = TMath::Sqrt(\r
-                                   TMath::Power( tDecayVertexV0[0] - tPrimaryVtxPosition[0] , 2) +\r
-                                   TMath::Power( tDecayVertexV0[1] - tPrimaryVtxPosition[1] , 2) +\r
-                                   TMath::Power( tDecayVertexV0[2] - tPrimaryVtxPosition[2] , 2)\r
-                                   );\r
-      lDistOverTotMom /= (lV0TotalMomentum+1e-10); //avoid division by zero, to be sure\r
-      \r
-      \r
-      // Getting invariant mass infos directly from ESD\r
-      lInvMassK0s        = v0->MassK0Short();\r
-      lInvMassLambda     = v0->MassLambda();\r
-      lInvMassAntiLambda = v0->MassAntiLambda();\r
-      lAlphaV0 = v0->AlphaV0();\r
-      lPtArmV0 = v0->PtArmV0();\r
-\r
-      //Official means of acquiring N-sigmas \r
-      Double_t lNSigmasPosProton = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kProton );\r
-      Double_t lNSigmasPosPion   = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kPion );\r
-      Double_t lNSigmasNegProton = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kProton );\r
-      Double_t lNSigmasNegPion   = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kPion );\r
-\r
-      //V0 QA histograms (before V0 selection)\r
-      fHistV0InvMassK0->Fill(lInvMassK0s);\r
-      fHistV0InvMassLambda->Fill(lInvMassLambda);\r
-      fHistV0InvMassAntiLambda->Fill(lInvMassAntiLambda);\r
-      fHistV0Armenteros->Fill(lAlphaV0,lPtArmV0);\r
-      \r
-      \r
-      //First Selection: Reject OnFly\r
-      if( (lOnFlyStatus == 0 && fkUseOnTheFly == kFALSE) || (lOnFlyStatus != 0 && fkUseOnTheFly == kTRUE ) ){\r
-       \r
-\r
-       //Second Selection: rough 20-sigma band, parametric.    \r
-       //K0Short: Enough to parametrize peak broadening with linear function.    \r
-       Double_t lUpperLimitK0Short = (5.63707e-01) + (1.14979e-02)*lPt; \r
-       Double_t lLowerLimitK0Short = (4.30006e-01) - (1.10029e-02)*lPt;\r
-       \r
-       //Lambda: Linear (for higher pt) plus exponential (for low-pt broadening)\r
-       //[0]+[1]*x+[2]*TMath::Exp(-[3]*x)\r
-       Double_t lUpperLimitLambda = (1.13688e+00) + (5.27838e-03)*lPt + (8.42220e-02)*TMath::Exp(-(3.80595e+00)*lPt); \r
-       Double_t lLowerLimitLambda = (1.09501e+00) - (5.23272e-03)*lPt - (7.52690e-02)*TMath::Exp(-(3.46339e+00)*lPt);\r
-       \r
-       //Do Selection      \r
-       if( (lInvMassLambda     < lUpperLimitLambda  && lInvMassLambda     > lLowerLimitLambda     ) || \r
-           (lInvMassAntiLambda < lUpperLimitLambda  && lInvMassAntiLambda > lLowerLimitLambda     ) || \r
-           (lInvMassK0s        < lUpperLimitK0Short && lInvMassK0s        > lLowerLimitK0Short    ) ){\r
-\r
-\r
-      //         //Pre-selection in case this is AA...\r
-      //         //if( fkIsNuclear == kFALSE ) fTree->Fill();\r
-      //         //if( fkIsNuclear == kTRUE){ \r
-      //         //If this is a nuclear collision___________________\r
-      //         // ... pre-filter with TPC, daughter eta selection\r
-\r
-         \r
-         if( (lInvMassLambda     < lUpperLimitLambda  && lInvMassLambda     > lLowerLimitLambda \r
-              && TMath::Abs(lNSigmasPosProton) < 6.0 && TMath::Abs(lNSigmasNegPion) < 6.0 ) || \r
-             (lInvMassAntiLambda < lUpperLimitLambda  && lInvMassAntiLambda > lLowerLimitLambda \r
-              && TMath::Abs(lNSigmasNegProton) < 6.0 && TMath::Abs(lNSigmasPosPion) < 6.0 ) ||  \r
-             (lInvMassK0s        < lUpperLimitK0Short && lInvMassK0s        > lLowerLimitK0Short \r
-              && TMath::Abs(lNSigmasNegPion)   < 6.0 && TMath::Abs(lNSigmasPosPion) < 6.0 ) ){\r
-           \r
-           //insane test\r
-           if ( TMath::Abs(lNegEta)<0.8 && TMath::Abs(lPosEta)<0.8 ){\r
-\r
-             // start the fine selection (usually done in post processing, but we don't have time to waste) --> Lambdas!\r
-             if(\r
-                TMath::Abs(lRap)<fRapidityBoundary &&\r
-                TMath::Abs(lNegEta)       <= fCutDaughterEta               &&                   \r
-                TMath::Abs(lPosEta)       <= fCutDaughterEta               &&\r
-                lV0Radius                 >= fCutV0Radius                  &&\r
-                lDcaNegToPrimVertex       >= fCutDCANegToPV                &&\r
-                lDcaPosToPrimVertex       >= fCutDCAPosToPV                &&\r
-                lDcaV0Daughters           <= fCutDCAV0Daughters            &&\r
-                lV0CosineOfPointingAngle  >= fCutV0CosPA                   && \r
-                fMassLambda*lDistOverTotMom    <= fCutProperLifetime       &&\r
-                lLeastNbrCrossedRows             >= fCutLeastNumberOfCrossedRows             &&\r
-                lLeastNbrCrossedRowsOverFindable >= fCutLeastNumberOfCrossedRowsOverFindable &&\r
-                lPtArmV0 * 5 < TMath::Abs(lAlphaV0)                        && \r
-                ((TMath::Abs(lNSigmasNegPion)   <= fCutTPCPIDNSigmasPion     &&\r
-                 TMath::Abs(lNSigmasPosProton) <= fCutTPCPIDNSigmasProton) ||\r
-                 (TMath::Abs(lNSigmasPosPion)   <= fCutTPCPIDNSigmasPion     &&\r
-                  TMath::Abs(lNSigmasNegProton) <= fCutTPCPIDNSigmasProton))            \r
-                )\r
-               {\r
-\r
-                 //V0 QA histograms (after V0 selection)\r
-                 fHistV0SelInvMassK0->Fill(lInvMassK0s);\r
-                 fHistV0SelInvMassLambda->Fill(lInvMassLambda);\r
-                 fHistV0SelInvMassAntiLambda->Fill(lInvMassAntiLambda);\r
-\r
-                 // this means a V0 candidate is found\r
-                 if(TMath::Abs(lInvMassLambda-fMassLambda) < fCutMassLambda ||\r
-                    TMath::Abs(lInvMassAntiLambda-fMassLambda) < fCutMassLambda){\r
-\r
-                   fHistV0SelArmenteros->Fill(lAlphaV0,lPtArmV0);                \r
-\r
-                   vEta    = lEta;\r
-                   vPhi    = lPhi;\r
-                   vPt     = lPt;\r
-                   if(lAlphaV0 > 0) vCharge = 1;\r
-                   if(lAlphaV0 < 0) vCharge = -1;\r
-\r
-                   // fill QA histograms\r
-                   fHistPt->Fill(vPt);\r
-                   fHistEta->Fill(vEta);\r
-                   fHistPhi->Fill(vPhi);\r
-                   \r
-                   // add the track to the TObjArray\r
-                   tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge,1.));\r
-                 }\r
-               }\r
-             }\r
-         }\r
-         //}//end nuclear_____________________________________\r
-       }\r
-      }\r
-    }//V0 loop\r
-  \r
-  return tracksAccepted;\r
-}\r
-\r
-//________________________________________________________________________\r
-TObjArray* AliAnalysisTaskTriggeredBF::GetShuffledTracks(TObjArray *tracks){\r
-  // Clones TObjArray and returns it with tracks after shuffling the charges\r
-\r
-  TObjArray* tracksShuffled = new TObjArray;\r
-  tracksShuffled->SetOwner(kTRUE);\r
-\r
-  vector<Short_t> *chargeVector = new vector<Short_t>;   //original charge of accepted tracks \r
-\r
-  for (Int_t i=0; i<tracks->GetEntriesFast(); i++)\r
-  {\r
-    AliVParticle* track = (AliVParticle*) tracks->At(i);\r
-    chargeVector->push_back(track->Charge());\r
-  }  \r
\r
-  random_shuffle(chargeVector->begin(), chargeVector->end());\r
-  \r
-  for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){\r
-    AliVParticle* track = (AliVParticle*) tracks->At(i);\r
-    tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i),1.));\r
-  }\r
-\r
-  delete chargeVector;\r
-   \r
-  return tracksShuffled;\r
-}\r
-\r
-//________________________________________________________________________\r
-void  AliAnalysisTaskTriggeredBF::FinishTaskOutput(){\r
-  //checks if Balance Function objects are there (needed to write the histograms)\r
-  if (!fBalance) {\r
-    AliError("fBalance not available");\r
-    return;\r
-  }  \r
-  if(fRunShuffling) {\r
-    if (!fShuffledBalance) {\r
-      AliError("fShuffledBalance not available");\r
-      return;\r
-    }\r
-  }\r
-\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskTriggeredBF::Terminate(Option_t *) {\r
-  // Called once at the end of the query\r
-\r
-  // not implemented ...\r
-\r
-}\r
-\r
-void AliAnalysisTaskTriggeredBF::UserExecMix(Option_t *)\r
-{\r
-\r
-  // not yet done for event mixing!\r
-  return;\r
-\r
-}\r
-\r
+#include <vector>
+#include "TChain.h"
+#include "TList.h"
+#include "TCanvas.h"
+#include "TLorentzVector.h"
+#include "TGraphErrors.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TArrayF.h"
+#include "TF1.h"
+#include "TRandom.h"
+
+#include "AliLog.h"
+
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisManager.h"
+
+#include "AliESDVertex.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODTrack.h"
+#include "AliAODInputHandler.h"
+#include "AliGenEventHeader.h"
+#include "AliGenHijingEventHeader.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliMixInputEventHandler.h"
+#include "AliStack.h"
+
+#include "TH2D.h"    
+#include "AliTHn.h"             
+
+#include "AliEventPoolManager.h" 
+
+#include "AliAnalysisTaskTriggeredBF.h"
+#include "AliBalanceTriggered.h"
+
+
+// Analysis task for the TriggeredBF code
+// Authors: Panos.Christakoglou@nikhef.nl, m.weber@cern.ch
+
+// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
+//
+// For the V0 part:
+// --> AliAnalysisTaskExtractV0AOD (by david.chinellato@gmail.com)
+//
+// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
+
+using std::cout;
+using std::endl;
+using std::vector;
+
+ClassImp(AliAnalysisTaskTriggeredBF)
+
+//________________________________________________________________________
+AliAnalysisTaskTriggeredBF::AliAnalysisTaskTriggeredBF(const char *name) 
+: AliAnalysisTaskSE(name), 
+  fBalance(0),
+  fRunShuffling(kFALSE),
+  fShuffledBalance(0),
+  fRunMixing(kFALSE),
+  fMixingTracks(50000),
+  fMixedBalance(0),
+  fPoolMgr(0),
+  fRunV0(kFALSE),
+  fPIDResponse(0),
+  fPIDCombined(0),
+  fList(0),
+  fListTriggeredBF(0),
+  fListTriggeredBFS(0),
+  fListTriggeredBFM(0),
+  fHistListPIDQA(0),
+  fHistListV0(0),
+  fHistEventStats(0),
+  fHistCentStats(0),
+  fHistTriggerStats(0),
+  fHistTrackStats(0),
+  fHistVx(0),
+  fHistVy(0),
+  fHistVz(0),
+  fHistClus(0),
+  fHistDCA(0),
+  fHistChi2(0),
+  fHistPt(0),
+  fHistEta(0),
+  fHistPhi(0),
+  fHistPhiBefore(0),
+  fHistPhiAfter(0),
+  fHistV0M(0),
+  fHistRefTracks(0),
+  fHistV0MultiplicityBeforeTrigSel(0),
+  fHistV0MultiplicityForTrigEvt(0),
+  fHistV0MultiplicityForSelEvt(0),
+  fHistV0MultiplicityForSelEvtNoTPCOnly(0),
+  fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
+  fHistMultiplicityBeforeTrigSel(0),
+  fHistMultiplicityForTrigEvt(0),
+  fHistMultiplicity(0),
+  fHistMultiplicityNoTPCOnly(0),
+  fHistMultiplicityNoTPCOnlyNoPileup(0),
+  fHistV0InvMassK0(0),
+  fHistV0InvMassLambda(0),
+  fHistV0InvMassAntiLambda(0),
+  fHistV0Armenteros(0),
+  fHistV0SelInvMassK0(0),
+  fHistV0SelInvMassLambda(0),
+  fHistV0SelInvMassAntiLambda(0),
+  fHistV0SelArmenteros(0),
+  fCentralityEstimator("V0M"),
+  fUseCentrality(kFALSE),
+  fCentralityPercentileMin(0.), 
+  fCentralityPercentileMax(5.),
+  fImpactParameterMin(0.),
+  fImpactParameterMax(20.),
+  fUseMultiplicity(kFALSE),
+  fNumberOfAcceptedTracksMin(0),
+  fNumberOfAcceptedTracksMax(10000),
+  fHistNumberOfAcceptedTracks(0),
+  fUseOfflineTrigger(kFALSE),
+  fVxMax(0.3),
+  fVyMax(0.3),
+  fVzMax(10.),
+  nAODtrackCutBit(128),
+  fPtMin(0.3),
+  fPtMax(1.5),
+  fEtaMin(-0.8),
+  fEtaMax(-0.8),
+  fDCAxyCut(-1),
+  fDCAzCut(-1),
+  fTPCchi2Cut(-1),
+  fNClustersTPCCut(-1)
+{
+  // Constructor
+  // Define input and output slots here
+  // Input slot #0 works with a TChain
+  DefineInput(0, TChain::Class());
+  // Output slot #0 writes into a TH1 container
+  DefineOutput(1, TList::Class());
+  DefineOutput(2, TList::Class());
+  DefineOutput(3, TList::Class());
+  DefineOutput(4, TList::Class());
+  DefineOutput(5, TList::Class());
+}
+
+//________________________________________________________________________
+AliAnalysisTaskTriggeredBF::~AliAnalysisTaskTriggeredBF() {
+
+  // Destructor
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskTriggeredBF::UserCreateOutputObjects() {
+  // Create histograms
+  // Called once
+
+  // global switch disabling the reference 
+  // (to avoid "Replacing existing TH1" if several wagons are created in train)
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+
+  if(!fBalance) {
+    fBalance = new AliBalanceTriggered();
+    fBalance->SetAnalysisLevel("AOD");
+  }
+  if(fRunShuffling) {
+    if(!fShuffledBalance) {
+      fShuffledBalance = new AliBalanceTriggered();
+      fShuffledBalance->SetAnalysisLevel("AOD");
+    }
+  }
+  if(fRunMixing) {
+    if(!fMixedBalance) {
+      fMixedBalance = new AliBalanceTriggered();
+      fMixedBalance->SetAnalysisLevel("AOD");
+    }
+  }
+
+  //QA list
+  fList = new TList();
+  fList->SetName("listQA");
+  fList->SetOwner();
+
+  //Balance Function list
+  fListTriggeredBF = new TList();
+  fListTriggeredBF->SetName("listTriggeredBF");
+  fListTriggeredBF->SetOwner();
+
+  if(fRunShuffling) {
+    fListTriggeredBFS = new TList();
+    fListTriggeredBFS->SetName("listTriggeredBFShuffled");
+    fListTriggeredBFS->SetOwner();
+  }
+  if(fRunMixing) {
+    fListTriggeredBFM = new TList();
+    fListTriggeredBFM->SetName("listTriggeredBFMixed");
+    fListTriggeredBFM->SetOwner();
+  }
+  
+  
+  //Event stats.
+  TString gCutName[4] = {"Total","Offline trigger",
+                         "Vertex","Analyzed"};
+  fHistEventStats = new TH1F("fHistEventStats",
+                             "Event statistics;;N_{events}",
+                             4,0.5,4.5);
+  for(Int_t i = 1; i <= 4; i++)
+    fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
+  fList->Add(fHistEventStats);
+  
+  TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};
+  fHistCentStats = new TH2F("fHistCentStats",
+                           "Centrality statistics;;Cent percentile",
+                           9,-0.5,8.5,220,-5,105);
+  for(Int_t i = 1; i <= 9; i++)
+    fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
+  fList->Add(fHistCentStats);
+  
+  fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);
+  fList->Add(fHistTriggerStats);
+  
+  fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);
+  fList->Add(fHistTrackStats);
+
+  fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);
+  fList->Add(fHistNumberOfAcceptedTracks);
+
+  // Vertex distributions
+  fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);
+  fList->Add(fHistVx);
+  fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);
+  fList->Add(fHistVy);
+  fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);
+  fList->Add(fHistVz);
+
+  // QA histograms
+  fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);
+  fList->Add(fHistClus);
+  fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);
+  fList->Add(fHistChi2);
+  fHistDCA  = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); 
+  fList->Add(fHistDCA);
+  fHistPt   = new TH1F("fHistPt","p_{T} distribution",200,0,10);
+  fList->Add(fHistPt);
+  fHistEta  = new TH1F("fHistEta","#eta distribution",200,-2,2);
+  fList->Add(fHistEta);
+  fHistPhi  = new TH1F("fHistPhi","#phi distribution",200,-20,380);
+  fList->Add(fHistPhi);
+  fHistPhiBefore  = new TH1F("fHistPhiBefore","#phi distribution",200,0.,2*TMath::Pi());
+  fList->Add(fHistPhiBefore);
+  fHistPhiAfter  = new TH1F("fHistPhiAfter","#phi distribution",200,0.,2*TMath::Pi());
+  fList->Add(fHistPhiAfter);
+  fHistV0M  = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);
+  fList->Add(fHistV0M);
+  TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};
+  fHistRefTracks  = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);
+  for(Int_t i = 1; i <= 6; i++)
+    fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());
+  fList->Add(fHistRefTracks);
+
+  //------------------------------------------------
+  // V0 Multiplicity Histograms
+  //------------------------------------------------
+  if(fRunV0){
+    fHistListV0 = new TList();
+    fHistListV0->SetOwner();  // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner
+    
+    if(! fHistV0MultiplicityBeforeTrigSel) {
+      fHistV0MultiplicityBeforeTrigSel = new TH1F("fHistV0MultiplicityBeforeTrigSel", 
+                                                 "V0s per event (before Trig. Sel.);Nbr of V0s/Evt;Events", 
+                                                 25, 0, 25);
+      fHistListV0->Add(fHistV0MultiplicityBeforeTrigSel);
+    }
+    
+    if(! fHistV0MultiplicityForTrigEvt) {
+      fHistV0MultiplicityForTrigEvt = new TH1F("fHistV0MultiplicityForTrigEvt", 
+                                              "V0s per event (for triggered evt);Nbr of V0s/Evt;Events", 
+                                              25, 0, 25);
+      fHistListV0->Add(fHistV0MultiplicityForTrigEvt);
+    }
+    
+    if(! fHistV0MultiplicityForSelEvt) {
+      fHistV0MultiplicityForSelEvt = new TH1F("fHistV0MultiplicityForSelEvt", 
+                                             "V0s per event;Nbr of V0s/Evt;Events", 
+                                             25, 0, 25);
+      fHistListV0->Add(fHistV0MultiplicityForSelEvt);
+    }
+    
+    if(! fHistV0MultiplicityForSelEvtNoTPCOnly) {
+    fHistV0MultiplicityForSelEvtNoTPCOnly = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnly", 
+                                                    "V0s per event;Nbr of V0s/Evt;Events", 
+                                                    25, 0, 25);
+    fHistListV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnly);
+    }
+    
+    if(! fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup) {
+      fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup", 
+                                                            "V0s per event;Nbr of V0s/Evt;Events", 
+                                                              25, 0, 25);
+      fHistListV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup);
+    }
+    
+    //------------------------------------------------
+    // Track Multiplicity Histograms
+    //------------------------------------------------
+    
+    if(! fHistMultiplicityBeforeTrigSel) {
+      fHistMultiplicityBeforeTrigSel = new TH1F("fHistMultiplicityBeforeTrigSel", 
+                                             "Tracks per event;Nbr of Tracks;Events", 
+                                               200, 0, 200);           
+      fHistListV0->Add(fHistMultiplicityBeforeTrigSel);
+    }
+    if(! fHistMultiplicityForTrigEvt) {
+      fHistMultiplicityForTrigEvt = new TH1F("fHistMultiplicityForTrigEvt", 
+                                            "Tracks per event;Nbr of Tracks;Events", 
+                                            200, 0, 200);              
+      fHistListV0->Add(fHistMultiplicityForTrigEvt);
+    }
+    if(! fHistMultiplicity) {
+      fHistMultiplicity = new TH1F("fHistMultiplicity", 
+                                  "Tracks per event;Nbr of Tracks;Events", 
+                                  200, 0, 200);                
+      fHistListV0->Add(fHistMultiplicity);
+    }
+    if(! fHistMultiplicityNoTPCOnly) {
+      fHistMultiplicityNoTPCOnly = new TH1F("fHistMultiplicityNoTPCOnly", 
+                                           "Tracks per event;Nbr of Tracks;Events", 
+                                           200, 0, 200);               
+    fHistListV0->Add(fHistMultiplicityNoTPCOnly);
+    }
+    if(! fHistMultiplicityNoTPCOnlyNoPileup) {
+      fHistMultiplicityNoTPCOnlyNoPileup = new TH1F("fHistMultiplicityNoTPCOnlyNoPileup", 
+                                                   "Tracks per event;Nbr of Tracks;Events", 
+                                                   200, 0, 200);               
+      fHistListV0->Add(fHistMultiplicityNoTPCOnlyNoPileup);
+    }
+
+    //------------------------------------------------
+    // V0 selection Histograms (before)
+    //------------------------------------------------
+    if(!fHistV0InvMassK0) {
+      fHistV0InvMassK0 = new TH1F("fHistV0InvMassK0",
+                                 "Invariant Mass for K0;Mass (GeV/c^{2});Events",
+                                 200,0,2);
+      fHistListV0->Add(fHistV0InvMassK0);
+    }
+    if(!fHistV0InvMassLambda) {
+      fHistV0InvMassLambda = new TH1F("fHistV0InvMassLambda",
+                                 "Invariant Mass for Lambda;Mass (GeV/c^{2});Events",
+                                 200,0,2);
+      fHistListV0->Add(fHistV0InvMassLambda);
+    }
+    if(!fHistV0InvMassAntiLambda) {
+      fHistV0InvMassAntiLambda = new TH1F("fHistV0InvMassAntiLambda",
+                                 "Invariant Mass for AntiLambda;Mass (GeV/c^{2});Events",
+                                 200,0,2);
+      fHistListV0->Add(fHistV0InvMassAntiLambda);
+    }
+    if(!fHistV0Armenteros) {
+      fHistV0Armenteros = new TH2F("fHistV0Armenteros",
+                                 "Armenteros plot;#alpha;q_{t}",
+                                  200,-1,1,200,0,0.5);
+      fHistListV0->Add(fHistV0Armenteros);
+    }
+    
+    //------------------------------------------------
+    // V0 selection Histograms (after)
+    //------------------------------------------------
+    if(!fHistV0SelInvMassK0) {
+      fHistV0SelInvMassK0 = new TH1F("fHistV0SelInvMassK0",
+                                 "Invariant Mass for K0;Mass (GeV/c^{2});Events",
+                                 200,0,2);
+      fHistListV0->Add(fHistV0SelInvMassK0);
+    }
+    if(!fHistV0SelInvMassLambda) {
+      fHistV0SelInvMassLambda = new TH1F("fHistV0SelInvMassLambda",
+                                 "Invariant Mass for Lambda;Mass (GeV/c^{2});Events",
+                                 200,0,2);
+      fHistListV0->Add(fHistV0SelInvMassLambda);
+    }
+    if(!fHistV0SelInvMassAntiLambda) {
+      fHistV0SelInvMassAntiLambda = new TH1F("fHistV0SelInvMassAntiLambda",
+                                 "Invariant Mass for AntiLambda;Mass (GeV/c^{2});Events",
+                                 200,0,2);
+      fHistListV0->Add(fHistV0SelInvMassAntiLambda);
+    }
+    if(!fHistV0SelArmenteros) {
+      fHistV0SelArmenteros = new TH2F("fHistV0SelArmenteros",
+                                 "Armenteros plot;#alpha;q_{t}",
+                                  200,-1,1,200,0,0.5);
+      fHistListV0->Add(fHistV0SelArmenteros);
+    }
+  }//V0
+    
+  // Balance function histograms
+  // Initialize histograms if not done yet
+  if(!fBalance->GetHistNp()){
+    AliWarning("Histograms not yet initialized! --> Will be done now");
+    AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
+    fBalance->InitHistograms();
+  }
+
+  if(fRunShuffling) {
+    if(!fShuffledBalance->GetHistNp()) {
+      AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");
+      AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
+      fShuffledBalance->InitHistograms();
+    }
+  }
+
+  if(fRunMixing) {
+    if(!fMixedBalance->GetHistNp()) {
+      AliWarning("Histograms (mixing) not yet initialized! --> Will be done now");
+      AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
+      fMixedBalance->InitHistograms();
+    }
+  }
+
+  fListTriggeredBF->Add(fBalance->GetHistNp());
+  fListTriggeredBF->Add(fBalance->GetHistNn());
+  fListTriggeredBF->Add(fBalance->GetHistNpn());
+  fListTriggeredBF->Add(fBalance->GetHistNnn());
+  fListTriggeredBF->Add(fBalance->GetHistNpp());
+  fListTriggeredBF->Add(fBalance->GetHistNnp());
+  
+  if(fRunShuffling) {
+    fListTriggeredBFS->Add(fShuffledBalance->GetHistNp());
+    fListTriggeredBFS->Add(fShuffledBalance->GetHistNn());
+    fListTriggeredBFS->Add(fShuffledBalance->GetHistNpn());
+    fListTriggeredBFS->Add(fShuffledBalance->GetHistNnn());
+    fListTriggeredBFS->Add(fShuffledBalance->GetHistNpp());
+    fListTriggeredBFS->Add(fShuffledBalance->GetHistNnp());
+  }  
+
+  if(fRunMixing) {
+    fListTriggeredBFM->Add(fMixedBalance->GetHistNp());
+    fListTriggeredBFM->Add(fMixedBalance->GetHistNn());
+    fListTriggeredBFM->Add(fMixedBalance->GetHistNpn());
+    fListTriggeredBFM->Add(fMixedBalance->GetHistNnn());
+    fListTriggeredBFM->Add(fMixedBalance->GetHistNpp());
+    fListTriggeredBFM->Add(fMixedBalance->GetHistNnp());
+  }  
+
+  // PID Response task active?
+  if(fRunV0) {
+    fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
+    if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");
+  }
+
+  // Event Mixing
+  Int_t trackDepth = fMixingTracks; 
+  Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
+   
+  Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,90.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
+  Double_t* centbins        = centralityBins;
+  Int_t nCentralityBins     = sizeof(centralityBins) / sizeof(Double_t) - 1;
+  
+  // bins for second buffer are shifted by 100 cm
+  Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
+  Double_t* vtxbins     = vertexBins;
+  Int_t nVertexBins     = sizeof(vertexBins) / sizeof(Double_t) - 1;
+
+  fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
+
+
+  // Post output data.
+  PostData(1, fList);
+  PostData(2, fListTriggeredBF);
+  if(fRunShuffling) PostData(3, fListTriggeredBFS);
+  if(fRunMixing) PostData(4, fListTriggeredBFM);
+  if(fRunV0) PostData(5,fHistListV0);
+
+  TH1::AddDirectory(oldStatus);
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskTriggeredBF::UserExec(Option_t *) {
+  // Main loop
+  // Called for each event
+
+  TString gAnalysisLevel = fBalance->GetAnalysisLevel();
+  Float_t fCentrality = -1.;  
+
+  // -------------------------------------------------------------                  
+  // AOD analysis (vertex and track cuts also here!!!!)
+  if(gAnalysisLevel == "AOD") {
+    AliVEvent* eventMain = dynamic_cast<AliVEvent*>(InputEvent()); 
+    if(!eventMain) {
+      AliError("eventMain not available");
+      return;
+    }
+
+    // check event cuts and fill event histograms
+    if((fCentrality = IsEventAccepted(eventMain)) < 0){
+      return;
+    }
+    
+    // get the accepted tracks in main event
+    TObjArray *tracksMain = NULL;
+    if(fRunV0) tracksMain = GetAcceptedV0s(eventMain);
+    else       tracksMain = GetAcceptedTracks(eventMain);
+
+    // store charges of all accepted tracks, shuffle and reassign (two extra loops!)
+    TObjArray* tracksShuffled = NULL;
+    if(fRunShuffling){
+      tracksShuffled = GetShuffledTracks(tracksMain);
+    }
+    
+    // Event mixing --> UPDATE POOL IS MISSING!!!
+    if (fRunMixing)
+      {
+        // 1. First get an event pool corresponding in mult (cent) and
+        //    zvertex to the current event. Once initialized, the pool
+        //    should contain nMix (reduced) events. This routine does not
+        //    pre-scan the chain. The first several events of every chain
+        //    will be skipped until the needed pools are filled to the
+        //    specified depth. If the pool categories are not too rare, this
+        //    should not be a problem. If they are rare, you could lose`
+       //    statistics.
+       
+       // 2. Collect the whole pool's content of tracks into one TObjArray
+       //    (bgTracks), which is effectively a single background super-event.
+       
+       // 3. The reduced and bgTracks arrays must both be passed into
+       //    FillCorrelations(). Also nMix should be passed in, so a weight
+       //    of 1./nMix can be applied.
+       
+       AliEventPool* pool = fPoolMgr->GetEventPool(fCentrality, eventMain->GetPrimaryVertex()->GetZ());
+       
+       if (!pool){
+         AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCentrality, eventMain->GetPrimaryVertex()->GetZ()));
+       }
+       else{
+
+       //pool->SetDebug(1);
+       
+         if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){ 
+           
+           
+           Int_t nMix = pool->GetCurrentNEvents();
+           //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
+           
+           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
+           //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
+           //if (pool->IsReady())
+           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
+           
+           // Fill mixed-event histos here  
+           for (Int_t jMix=0; jMix<nMix; jMix++) 
+             {
+               TObjArray* tracksMixed = pool->GetEvent(jMix);
+               fMixedBalance->FillBalance(fCentrality,tracksMain,tracksMixed); 
+             }
+         }
+         
+         // Update the Event pool
+         pool->UpdatePool(tracksMain);
+         //pool->PrintInfo();
+         
+       }//pool NULL check  
+      }//run mixing
+    
+    // calculate balance function
+    fBalance->FillBalance(fCentrality,tracksMain,NULL);
+    
+    // calculate shuffled balance function
+    if(fRunShuffling && tracksShuffled != NULL) {
+      fShuffledBalance->FillBalance(fCentrality,tracksShuffled,NULL);
+    }
+    
+  }//AOD analysis
+  else{
+    AliError("Triggered Balance Function analysis only for AODs!");
+  }
+}     
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskTriggeredBF::IsEventAccepted(AliVEvent *event){
+  // Checks the Event cuts
+  // Fills Event statistics histograms
+  
+  // event selection done in AliAnalysisTaskSE::Exec() --> this is not used
+  fHistEventStats->Fill(1); //all events
+
+  Bool_t isSelectedMain = kTRUE;
+  Float_t fCentrality = -1.;
+  Int_t nV0s          = event->GetNumberOfV0s();
+  TString gAnalysisLevel = fBalance->GetAnalysisLevel();
+  
+  if(fUseOfflineTrigger)
+    isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+  
+  //V0 QA histograms (before trigger selection)
+  if(fRunV0){
+    fHistMultiplicityBeforeTrigSel->Fill ( -1 );
+    fHistV0MultiplicityBeforeTrigSel->Fill ( nV0s );
+  }
+  
+  if(isSelectedMain) {
+    fHistEventStats->Fill(2); //triggered events
+    
+    //Centrality stuff 
+    if(fUseCentrality) {
+      if(gAnalysisLevel == "AOD") { //centrality in AOD header
+       AliAODHeader *header = (AliAODHeader*) event->GetHeader();
+       fCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
+
+       // QA for centrality estimators
+       fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));
+       fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("FMD"));
+       fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("TRK"));
+       fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("TKL"));
+       fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("CL0"));
+       fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("CL1"));
+       fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
+       fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
+       fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
+       
+       // centrality QA (V0M)
+       fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());
+       
+       // centrality QA (reference tracks)
+       fHistRefTracks->Fill(0.,header->GetRefMultiplicity());
+       fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());
+       fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());
+       fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());
+       fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));
+       fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));
+       fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));
+       fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));
+       fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));
+
+       //V0 QA histograms (after trigger selection)
+       if(fRunV0){
+         fHistMultiplicityForTrigEvt->Fill ( fCentrality );
+         fHistV0MultiplicityForTrigEvt->Fill ( nV0s );
+       }
+      }
+    }
+    
+    
+    const AliVVertex *vertex = event->GetPrimaryVertex();
+    
+    if(vertex) {
+      Double32_t fCov[6];
+      vertex->GetCovarianceMatrix(fCov);
+      if(vertex->GetNContributors() > 0) {
+       if(fCov[5] != 0) {
+         fHistEventStats->Fill(3); //events with a proper vertex
+         if(TMath::Abs(vertex->GetX()) < fVxMax) {
+           if(TMath::Abs(vertex->GetY()) < fVyMax) {
+             if(TMath::Abs(vertex->GetZ()) < fVzMax) {
+               fHistEventStats->Fill(4); //analyzed events
+               fHistVx->Fill(vertex->GetX());
+               fHistVy->Fill(vertex->GetY());
+               fHistVz->Fill(vertex->GetZ());
+
+
+
+               //V0 QA histograms (vertex Z check)
+               if(fRunV0){
+                 fHistV0MultiplicityForSelEvt ->Fill( nV0s );
+                 fHistMultiplicity->Fill(fCentrality);
+
+                 //V0 QA histograms (Only look at events with well-established PV)
+                 const AliAODVertex *lPrimarySPDVtx = ((AliAODEvent*)event)->GetPrimaryVertexSPD();
+                 if(lPrimarySPDVtx){
+                   fHistMultiplicityNoTPCOnly->Fill ( fCentrality );
+                   fHistV0MultiplicityForSelEvtNoTPCOnly->Fill ( nV0s );
+                   
+                   //V0 QA histograms (Pileup Rejection)
+                   // FIXME : quality selection regarding pile-up rejection 
+                   fHistMultiplicityNoTPCOnlyNoPileup->Fill(fCentrality);
+                   fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup ->Fill( nV0s );
+
+                 }
+                 else{
+                   return -1;
+                 }
+               }
+               
+               
+               // take only events inside centrality class
+               if((fCentrality > fCentralityPercentileMin) && (fCentrality < fCentralityPercentileMax)){
+                 return fCentrality;           
+               }//centrality class
+             }//Vz cut
+           }//Vy cut
+         }//Vx cut
+       }//proper vertex resolution
+      }//proper number of contributors
+    }//vertex object valid
+  }//triggered event 
+  
+  // in all other cases return -1 (event not accepted)
+  return -1;
+}
+
+//________________________________________________________________________
+TObjArray* AliAnalysisTaskTriggeredBF::GetAcceptedTracks(AliVEvent *event){
+  // Returns TObjArray with tracks after all track cuts (only for AOD!)
+  // Fills QA histograms
+
+  //output TObjArray holding all good tracks
+  TObjArray* tracksAccepted = new TObjArray;
+  tracksAccepted->SetOwner(kTRUE);
+
+  Short_t vCharge = 0;
+  Double_t vEta    = 0.;
+  Double_t vPhi    = 0.;
+  Double_t vPt     = 0.;
+  
+  // Loop over tracks in event
+  for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {
+    AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));
+    if (!aodTrack) {
+      AliError(Form("Could not receive track %d", iTracks));
+      continue;
+    }
+    
+    // AOD track cuts
+    
+    // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
+    // take only TPC only tracks 
+    fHistTrackStats->Fill(aodTrack->GetFilterMap());
+    if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;
+    
+    vCharge = aodTrack->Charge();
+    vEta    = aodTrack->Eta();
+    vPhi    = aodTrack->Phi() * TMath::RadToDeg();
+    vPt     = aodTrack->Pt();
+    
+    Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)
+    Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)
+    
+    
+    // Kinematics cuts from ESD track cuts
+    if( vPt < fPtMin || vPt > fPtMax)      continue;
+    if( vEta < fEtaMin || vEta > fEtaMax)  continue;
+    
+    // Extra DCA cuts (for systematic studies [!= -1])
+    if( fDCAxyCut != -1 && fDCAzCut != -1){
+      if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
+       continue;  // 2D cut
+      }
+    }
+    
+    // Extra TPC cuts (for systematic studies [!= -1])
+    if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
+      continue;
+    }
+    if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
+      continue;
+    }
+    
+    // fill QA histograms
+    fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
+    fHistDCA->Fill(dcaZ,dcaXY);
+    fHistChi2->Fill(aodTrack->Chi2perNDF());
+    fHistPt->Fill(vPt);
+    fHistEta->Fill(vEta);
+    fHistPhi->Fill(vPhi);
+    
+    // add the track to the TObjArray
+    tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge,1.));
+  }
+
+  return tracksAccepted;
+}
+
+//________________________________________________________________________
+TObjArray* AliAnalysisTaskTriggeredBF::GetAcceptedV0s(AliVEvent *event){
+  // Returns TObjArray with tracks after all track cuts (only for AOD!)
+  // Fills QA histograms
+
+  //output TObjArray holding all good tracks
+  TObjArray* tracksAccepted = new TObjArray;
+  tracksAccepted->SetOwner(kTRUE);
+
+  Short_t vCharge = 0;
+  Double_t vEta    = 0.;
+  Double_t vPhi    = 0.;
+  Double_t vPt     = 0.;
+  
+  //------------------------------------------------
+  // MAIN LAMBDA LOOP STARTS HERE (basically a copy of AliAnalysisTaskExtractV0AOD)
+  //------------------------------------------------
+
+  // parameters (for the time being hard coded here) --> from David for EbyE Lambdas
+  Bool_t fkUseOnTheFly = kFALSE;
+  Double_t fRapidityBoundary  = 0.5; 
+  Double_t fCutDaughterEta    = 0.8;
+  Double_t fCutV0Radius       = 0.9;
+  Double_t fCutDCANegToPV     = 0.1;
+  Double_t fCutDCAPosToPV     = 0.1;
+  Double_t fCutDCAV0Daughters = 1.0;
+  Double_t fCutV0CosPA        = 0.9995;
+  Double_t fMassLambda        = 1.115683;
+  Double_t fCutMassLambda     = 0.007;
+  Double_t fCutProperLifetime = 3*7.9;
+  Double_t fCutLeastNumberOfCrossedRows = 70;
+  Double_t fCutLeastNumberOfCrossedRowsOverFindable = 0.8;
+  Double_t fCutTPCPIDNSigmasProton  = 3.0;
+  Double_t fCutTPCPIDNSigmasPion    = 5.0;
+
+
+  //Variable definition
+  Int_t    lOnFlyStatus = 0;// nv0sOn = 0, nv0sOff = 0;
+  Double_t lChi2V0 = 0;
+  Double_t lDcaV0Daughters = 0, lDcaV0ToPrimVertex = 0;
+  Double_t lDcaPosToPrimVertex = 0, lDcaNegToPrimVertex = 0;
+  Double_t lV0CosineOfPointingAngle = 0;
+  Double_t lV0Radius = 0, lPt = 0;
+  Double_t lEta = 0, lPhi = 0;
+  Double_t lRap = 0, lRapK0Short = 0, lRapLambda = 0;
+  Double_t lInvMassK0s = 0, lInvMassLambda = 0, lInvMassAntiLambda = 0;
+  Double_t lAlphaV0 = 0, lPtArmV0 = 0;
+  
+  Double_t fMinV0Pt = 0; 
+  Double_t fMaxV0Pt = 100; 
+  
+
+  
+  // some event observables
+  Int_t nv0s = event->GetNumberOfV0s();
+  Double_t tPrimaryVtxPosition[3];
+  const AliVVertex *primaryVtx = event->GetPrimaryVertex();
+  tPrimaryVtxPosition[0] = primaryVtx->GetX();
+  tPrimaryVtxPosition[1] = primaryVtx->GetY();
+  tPrimaryVtxPosition[2] = primaryVtx->GetZ();
+
+
+  //loop over V0s  
+  for (Int_t iV0 = 0; iV0 < nv0s; iV0++) 
+    {// This is the begining of the V0 loop
+      AliAODv0 *v0 = ((AliAODEvent*)event)->GetV0(iV0);
+      if (!v0) continue;
+
+      //Obsolete at AOD level... 
+      //---> Fix On-the-Fly candidates, count how many swapped
+      //if( v0->GetParamN()->Charge() > 0 && v0->GetParamP()->Charge() < 0 ){
+      //  fHistSwappedV0Counter -> Fill( 1 );
+      //}else{
+      //  fHistSwappedV0Counter -> Fill( 0 ); 
+      //}
+      //if ( fkUseOnTheFly ) CheckChargeV0(v0); 
+      
+      Double_t tDecayVertexV0[3]; v0->GetXYZ(tDecayVertexV0); 
+      Double_t tV0mom[3];
+      v0->GetPxPyPz( tV0mom ); 
+      Double_t lV0TotalMomentum = TMath::Sqrt(
+                                             tV0mom[0]*tV0mom[0]+tV0mom[1]*tV0mom[1]+tV0mom[2]*tV0mom[2] );
+      
+      lV0Radius = TMath::Sqrt(tDecayVertexV0[0]*tDecayVertexV0[0]+tDecayVertexV0[1]*tDecayVertexV0[1]);
+      lPt = v0->Pt();
+      lEta = v0->Eta();
+      lPhi = v0->Phi()*TMath::RadToDeg();
+      lRapK0Short = v0->RapK0Short();
+      lRapLambda  = v0->RapLambda();
+      lRap        = lRapLambda;//v0->Y(); //FIXME!!!
+      if ((lPt<fMinV0Pt)||(fMaxV0Pt<lPt)) continue;
+      
+      //UInt_t lKeyPos = (UInt_t)TMath::Abs(v0->GetPosID());
+      //UInt_t lKeyNeg = (UInt_t)TMath::Abs(v0->GetPosID());
+
+      Double_t lMomPos[3]; //v0->GetPPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);
+      Double_t lMomNeg[3]; //v0->GetNPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
+      lMomPos[0] = v0->MomPosX();
+      lMomPos[1] = v0->MomPosY();
+      lMomPos[2] = v0->MomPosZ();
+      lMomNeg[0] = v0->MomNegX();
+      lMomNeg[1] = v0->MomNegY();
+      lMomNeg[2] = v0->MomNegZ();
+      
+      AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
+      AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
+      if (!pTrack || !nTrack) {
+       AliError("ERROR: Could not retreive one of the daughter track");
+       continue;
+      }
+
+      //Daughter Eta for Eta selection, afterwards
+      Double_t lNegEta = nTrack->Eta();
+      Double_t lPosEta = pTrack->Eta();
+      
+      // Filter like-sign V0 (next: add counter and distribution)
+      if ( pTrack->Charge() == nTrack->Charge()){
+       continue;
+      } 
+      
+      //Quick test this far! 
+      
+
+      //________________________________________________________________________
+      // Track quality cuts 
+      Float_t lPosTrackCrossedRows = pTrack->GetTPCClusterInfo(2,1);
+      Float_t lNegTrackCrossedRows = nTrack->GetTPCClusterInfo(2,1);
+      Float_t lLeastNbrCrossedRows =  (lPosTrackCrossedRows>lNegTrackCrossedRows) ? lNegTrackCrossedRows : lPosTrackCrossedRows;
+
+      // TPC refit condition (done during reconstruction for Offline but not for On-the-fly)
+      if( !(pTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;
+      if( !(nTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;
+      
+      if ( ( ( pTrack->GetTPCClusterInfo(2,1) ) < 70 ) || ( ( nTrack->GetTPCClusterInfo(2,1) ) < 70 ) ) continue;
+      
+      //Findable clusters > 0 condition
+      if( pTrack->GetTPCNclsF()<=0 || nTrack->GetTPCNclsF()<=0 ) continue;
+      
+      //Compute ratio Crossed Rows / Findable clusters
+      //Note: above test avoids division by zero! 
+      Float_t lPosTrackCrossedRowsOverFindable = lPosTrackCrossedRows / ((double)(pTrack->GetTPCNclsF())); 
+      Float_t lNegTrackCrossedRowsOverFindable = lNegTrackCrossedRows / ((double)(nTrack->GetTPCNclsF())); 
+      Float_t lLeastNbrCrossedRowsOverFindable = (lPosTrackCrossedRowsOverFindable>lNegTrackCrossedRowsOverFindable) ? lNegTrackCrossedRowsOverFindable : lPosTrackCrossedRowsOverFindable;
+
+      //Lowest Cut Level for Ratio Crossed Rows / Findable = 0.8, set here
+      if ( lLeastNbrCrossedRowsOverFindable < 0.8) continue;
+      
+      //End track Quality Cuts
+      //________________________________________________________________________
+      
+      
+      lDcaPosToPrimVertex = v0->DcaPosToPrimVertex();
+      lDcaNegToPrimVertex = v0->DcaNegToPrimVertex();
+          
+      lOnFlyStatus = v0->GetOnFlyStatus();
+      lChi2V0 = v0->Chi2V0();
+      lDcaV0Daughters = v0->DcaV0Daughters();
+      lDcaV0ToPrimVertex = v0->DcaV0ToPrimVertex();
+      lV0CosineOfPointingAngle = v0->CosPointingAngle(tPrimaryVtxPosition);
+      
+      // Distance over total momentum
+      Double_t lDistOverTotMom = TMath::Sqrt(
+                                   TMath::Power( tDecayVertexV0[0] - tPrimaryVtxPosition[0] , 2) +
+                                   TMath::Power( tDecayVertexV0[1] - tPrimaryVtxPosition[1] , 2) +
+                                   TMath::Power( tDecayVertexV0[2] - tPrimaryVtxPosition[2] , 2)
+                                   );
+      lDistOverTotMom /= (lV0TotalMomentum+1e-10); //avoid division by zero, to be sure
+      
+      
+      // Getting invariant mass infos directly from ESD
+      lInvMassK0s        = v0->MassK0Short();
+      lInvMassLambda     = v0->MassLambda();
+      lInvMassAntiLambda = v0->MassAntiLambda();
+      lAlphaV0 = v0->AlphaV0();
+      lPtArmV0 = v0->PtArmV0();
+
+      //Official means of acquiring N-sigmas 
+      Double_t lNSigmasPosProton = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kProton );
+      Double_t lNSigmasPosPion   = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kPion );
+      Double_t lNSigmasNegProton = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kProton );
+      Double_t lNSigmasNegPion   = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kPion );
+
+      //V0 QA histograms (before V0 selection)
+      fHistV0InvMassK0->Fill(lInvMassK0s);
+      fHistV0InvMassLambda->Fill(lInvMassLambda);
+      fHistV0InvMassAntiLambda->Fill(lInvMassAntiLambda);
+      fHistV0Armenteros->Fill(lAlphaV0,lPtArmV0);
+      
+      
+      //First Selection: Reject OnFly
+      if( (lOnFlyStatus == 0 && fkUseOnTheFly == kFALSE) || (lOnFlyStatus != 0 && fkUseOnTheFly == kTRUE ) ){
+       
+
+       //Second Selection: rough 20-sigma band, parametric.    
+       //K0Short: Enough to parametrize peak broadening with linear function.    
+       Double_t lUpperLimitK0Short = (5.63707e-01) + (1.14979e-02)*lPt; 
+       Double_t lLowerLimitK0Short = (4.30006e-01) - (1.10029e-02)*lPt;
+       
+       //Lambda: Linear (for higher pt) plus exponential (for low-pt broadening)
+       //[0]+[1]*x+[2]*TMath::Exp(-[3]*x)
+       Double_t lUpperLimitLambda = (1.13688e+00) + (5.27838e-03)*lPt + (8.42220e-02)*TMath::Exp(-(3.80595e+00)*lPt); 
+       Double_t lLowerLimitLambda = (1.09501e+00) - (5.23272e-03)*lPt - (7.52690e-02)*TMath::Exp(-(3.46339e+00)*lPt);
+       
+       //Do Selection      
+       if( (lInvMassLambda     < lUpperLimitLambda  && lInvMassLambda     > lLowerLimitLambda     ) || 
+           (lInvMassAntiLambda < lUpperLimitLambda  && lInvMassAntiLambda > lLowerLimitLambda     ) || 
+           (lInvMassK0s        < lUpperLimitK0Short && lInvMassK0s        > lLowerLimitK0Short    ) ){
+
+
+      //         //Pre-selection in case this is AA...
+      //         //if( fkIsNuclear == kFALSE ) fTree->Fill();
+      //         //if( fkIsNuclear == kTRUE){ 
+      //         //If this is a nuclear collision___________________
+      //         // ... pre-filter with TPC, daughter eta selection
+
+         
+         if( (lInvMassLambda     < lUpperLimitLambda  && lInvMassLambda     > lLowerLimitLambda 
+              && TMath::Abs(lNSigmasPosProton) < 6.0 && TMath::Abs(lNSigmasNegPion) < 6.0 ) || 
+             (lInvMassAntiLambda < lUpperLimitLambda  && lInvMassAntiLambda > lLowerLimitLambda 
+              && TMath::Abs(lNSigmasNegProton) < 6.0 && TMath::Abs(lNSigmasPosPion) < 6.0 ) ||  
+             (lInvMassK0s        < lUpperLimitK0Short && lInvMassK0s        > lLowerLimitK0Short 
+              && TMath::Abs(lNSigmasNegPion)   < 6.0 && TMath::Abs(lNSigmasPosPion) < 6.0 ) ){
+           
+           //insane test
+           if ( TMath::Abs(lNegEta)<0.8 && TMath::Abs(lPosEta)<0.8 ){
+
+             // start the fine selection (usually done in post processing, but we don't have time to waste) --> Lambdas!
+             if(
+                TMath::Abs(lRap)<fRapidityBoundary &&
+                TMath::Abs(lNegEta)       <= fCutDaughterEta               &&                   
+                TMath::Abs(lPosEta)       <= fCutDaughterEta               &&
+                lV0Radius                 >= fCutV0Radius                  &&
+                lDcaNegToPrimVertex       >= fCutDCANegToPV                &&
+                lDcaPosToPrimVertex       >= fCutDCAPosToPV                &&
+                lDcaV0Daughters           <= fCutDCAV0Daughters            &&
+                lV0CosineOfPointingAngle  >= fCutV0CosPA                   && 
+                fMassLambda*lDistOverTotMom    <= fCutProperLifetime       &&
+                lLeastNbrCrossedRows             >= fCutLeastNumberOfCrossedRows             &&
+                lLeastNbrCrossedRowsOverFindable >= fCutLeastNumberOfCrossedRowsOverFindable &&
+                lPtArmV0 * 5 < TMath::Abs(lAlphaV0)                        && 
+                ((TMath::Abs(lNSigmasNegPion)   <= fCutTPCPIDNSigmasPion     &&
+                 TMath::Abs(lNSigmasPosProton) <= fCutTPCPIDNSigmasProton) ||
+                 (TMath::Abs(lNSigmasPosPion)   <= fCutTPCPIDNSigmasPion     &&
+                  TMath::Abs(lNSigmasNegProton) <= fCutTPCPIDNSigmasProton))            
+                )
+               {
+
+                 //V0 QA histograms (after V0 selection)
+                 fHistV0SelInvMassK0->Fill(lInvMassK0s);
+                 fHistV0SelInvMassLambda->Fill(lInvMassLambda);
+                 fHistV0SelInvMassAntiLambda->Fill(lInvMassAntiLambda);
+
+                 // this means a V0 candidate is found
+                 if(TMath::Abs(lInvMassLambda-fMassLambda) < fCutMassLambda ||
+                    TMath::Abs(lInvMassAntiLambda-fMassLambda) < fCutMassLambda){
+
+                   fHistV0SelArmenteros->Fill(lAlphaV0,lPtArmV0);                
+
+                   vEta    = lEta;
+                   vPhi    = lPhi;
+                   vPt     = lPt;
+                   if(lAlphaV0 > 0) vCharge = 1;
+                   if(lAlphaV0 < 0) vCharge = -1;
+
+                   // fill QA histograms
+                   fHistPt->Fill(vPt);
+                   fHistEta->Fill(vEta);
+                   fHistPhi->Fill(vPhi);
+                   
+                   // add the track to the TObjArray
+                   tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge,1.));
+                 }
+               }
+             }
+         }
+         //}//end nuclear_____________________________________
+       }
+      }
+    }//V0 loop
+  
+  return tracksAccepted;
+}
+
+//________________________________________________________________________
+TObjArray* AliAnalysisTaskTriggeredBF::GetShuffledTracks(TObjArray *tracks){
+  // Clones TObjArray and returns it with tracks after shuffling the charges
+
+  TObjArray* tracksShuffled = new TObjArray;
+  tracksShuffled->SetOwner(kTRUE);
+
+  vector<Short_t> *chargeVector = new vector<Short_t>;   //original charge of accepted tracks 
+
+  for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
+  {
+    AliVParticle* track = (AliVParticle*) tracks->At(i);
+    chargeVector->push_back(track->Charge());
+  }  
+  random_shuffle(chargeVector->begin(), chargeVector->end());
+  
+  for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){
+    AliVParticle* track = (AliVParticle*) tracks->At(i);
+    tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i),1.));
+  }
+
+  delete chargeVector;
+   
+  return tracksShuffled;
+}
+
+//________________________________________________________________________
+void  AliAnalysisTaskTriggeredBF::FinishTaskOutput(){
+  //checks if Balance Function objects are there (needed to write the histograms)
+  if (!fBalance) {
+    AliError("fBalance not available");
+    return;
+  }  
+  if(fRunShuffling) {
+    if (!fShuffledBalance) {
+      AliError("fShuffledBalance not available");
+      return;
+    }
+  }
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskTriggeredBF::Terminate(Option_t *) {
+  // Called once at the end of the query
+
+  // not implemented ...
+
+}
+
+void AliAnalysisTaskTriggeredBF::UserExecMix(Option_t *)
+{
+
+  // not yet done for event mixing!
+  return;
+
+}
+