Event Mixing for Triggered BF analysis
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliAnalysisTaskTriggeredBF.cxx
index e781e7c354c07b0ce1522b0ac32b7c7616bbbab8..c1c6153cffe3441e6e4fe72783f6dd0ff5f3a41d 100755 (executable)
@@ -28,7 +28,9 @@
 #include "AliStack.h"\r
 \r
 #include "TH2D.h"    \r
-#include "AliTHn.h"              \r
+#include "AliTHn.h"             \r
+\r
+#include "AliEventPoolManager.h" \r
 \r
 #include "AliAnalysisTaskTriggeredBF.h"\r
 #include "AliBalanceTriggered.h"\r
@@ -45,9 +47,14 @@ AliAnalysisTaskTriggeredBF::AliAnalysisTaskTriggeredBF(const char *name)
   fBalance(0),\r
   fRunShuffling(kFALSE),\r
   fShuffledBalance(0),\r
+  fRunMixing(kFALSE),\r
+  fMixingTracks(50000),\r
+  fMixedBalance(0),\r
+  fPoolMgr(0),\r
   fList(0),\r
   fListTriggeredBF(0),\r
   fListTriggeredBFS(0),\r
+  fListTriggeredBFM(0),\r
   fHistListPIDQA(0),\r
   fHistEventStats(0),\r
   fHistCentStats(0),\r
@@ -98,6 +105,7 @@ AliAnalysisTaskTriggeredBF::AliAnalysisTaskTriggeredBF(const char *name)
   DefineOutput(1, TList::Class());\r
   DefineOutput(2, TList::Class());\r
   DefineOutput(3, TList::Class());\r
+  DefineOutput(4, TList::Class());\r
 }\r
 \r
 //________________________________________________________________________\r
@@ -121,6 +129,12 @@ void AliAnalysisTaskTriggeredBF::UserCreateOutputObjects() {
       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
@@ -137,6 +151,11 @@ void AliAnalysisTaskTriggeredBF::UserCreateOutputObjects() {
     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
@@ -217,6 +236,14 @@ void AliAnalysisTaskTriggeredBF::UserCreateOutputObjects() {
     }\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
@@ -233,13 +260,38 @@ void AliAnalysisTaskTriggeredBF::UserCreateOutputObjects() {
     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
+\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
+  cout<<nCentralityBins<<" "<<nVertexBins<<endl;\r
+\r
+  fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);\r
 \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
 }\r
 \r
 //________________________________________________________________________\r
@@ -248,183 +300,276 @@ void AliAnalysisTaskTriggeredBF::UserExec(Option_t *) {
   // Called for each event\r
 \r
   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
-\r
-  Float_t fCentrality           = 0.;\r
-  \r
-  // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)\r
-  vector<Double_t> *chargeVector[9];          // original charge\r
-  for(Int_t i = 0; i < 9; i++){\r
-    chargeVector[i]        = new vector<Double_t>;\r
-  }\r
-  \r
-  Double_t v_charge;\r
-  Double_t v_y;\r
-  Double_t v_eta;\r
-  Double_t v_phi;\r
-  Double_t v_p[3];\r
-  Double_t v_pt;\r
-  Double_t v_E;\r
+  Float_t fCentrality = 0.;  \r
 \r
   // -------------------------------------------------------------                  \r
   // AOD analysis (vertex and track cuts also here!!!!)\r
   if(gAnalysisLevel == "AOD") {\r
-    AliAODEvent* aodEventMain = dynamic_cast<AliAODEvent*>(InputEvent()); \r
-    if(!aodEventMain) {\r
-      AliError("aodEventMain not available");\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 = 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
-    AliAODHeader *aodHeaderMain = aodEventMain->GetHeader();\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
+       //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
+      }//run mixing\r
     \r
-    // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
-    fHistEventStats->Fill(1); //all events\r
+    // calculate balance function\r
+    fBalance->FillBalance(fCentrality,tracksMain,NULL);\r
     \r
-    Bool_t isSelectedMain = kTRUE;\r
+    // calculate shuffled balance function\r
+    if(fRunShuffling && tracksShuffled != NULL) {\r
+      fShuffledBalance->FillBalance(fCentrality,tracksShuffled,NULL);\r
+    }\r
     \r
-    if(fUseOfflineTrigger)\r
-      isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\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
+  TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
+  \r
+  if(fUseOfflineTrigger)\r
+    isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
+  \r
+  if(isSelectedMain) {\r
+    fHistEventStats->Fill(2); //triggered events\r
     \r
-    if(isSelectedMain) {\r
-      fHistEventStats->Fill(2); //triggered events\r
-      \r
-      //Centrality stuff (centrality in AOD header)\r
-      if(fUseCentrality) {\r
-       fCentrality = aodHeaderMain->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\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.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0M"));\r
-       fHistCentStats->Fill(1.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("FMD"));\r
-       fHistCentStats->Fill(2.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TRK"));\r
-       fHistCentStats->Fill(3.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKL"));\r
-       fHistCentStats->Fill(4.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL0"));\r
-       fHistCentStats->Fill(5.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL1"));\r
-       fHistCentStats->Fill(6.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
-       fHistCentStats->Fill(7.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
-       fHistCentStats->Fill(8.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
-       \r
-       // take only events inside centrality class\r
-       if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax)) \r
-         return;\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(aodEventMain->GetVZEROData()->GetMTotV0A(), aodEventMain->GetVZEROData()->GetMTotV0C());\r
+       fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
        \r
        // centrality QA (reference tracks)\r
-       fHistRefTracks->Fill(0.,aodHeaderMain->GetRefMultiplicity());\r
-       fHistRefTracks->Fill(1.,aodHeaderMain->GetRefMultiplicityPos());\r
-       fHistRefTracks->Fill(2.,aodHeaderMain->GetRefMultiplicityNeg());\r
-       fHistRefTracks->Fill(3.,aodHeaderMain->GetTPConlyRefMultiplicity());\r
-       fHistRefTracks->Fill(4.,aodHeaderMain->GetNumberOfITSClusters(0));\r
-       fHistRefTracks->Fill(5.,aodHeaderMain->GetNumberOfITSClusters(1));\r
-       fHistRefTracks->Fill(6.,aodHeaderMain->GetNumberOfITSClusters(2));\r
-       fHistRefTracks->Fill(7.,aodHeaderMain->GetNumberOfITSClusters(3));\r
-       fHistRefTracks->Fill(8.,aodHeaderMain->GetNumberOfITSClusters(4));\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
-      \r
-      const AliAODVertex *vertexMain = aodEventMain->GetPrimaryVertex();\r
-      \r
-      if(vertexMain) {\r
-       Double32_t fCovMain[6];\r
-       vertexMain->GetCovarianceMatrix(fCovMain);\r
-       \r
-       if(vertexMain->GetNContributors() > 0) {\r
-         if(fCovMain[5] != 0) {\r
-           fHistEventStats->Fill(3); //events with a proper vertex\r
-           if(TMath::Abs(vertexMain->GetX()) < fVxMax) {\r
-             if(TMath::Abs(vertexMain->GetY()) < fVyMax) {\r
-               if(TMath::Abs(vertexMain->GetZ()) < fVzMax) {\r
-                 fHistEventStats->Fill(4); //analyzed events\r
-                 fHistVx->Fill(vertexMain->GetX());\r
-                 fHistVy->Fill(vertexMain->GetY());\r
-                 fHistVz->Fill(vertexMain->GetZ());\r
-                 \r
-                 // Loop over tracks in main event\r
-                 for (Int_t iTracksMain = 0; iTracksMain < aodEventMain->GetNumberOfTracks(); iTracksMain++) {\r
-                   AliAODTrack* aodTrackMain = dynamic_cast<AliAODTrack *>(aodEventMain->GetTrack(iTracksMain));\r
-                   if (!aodTrackMain) {\r
-                     AliError(Form("Could not receive track %d", iTracksMain));\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(aodTrackMain->GetFilterMap());\r
-                   if(!aodTrackMain->TestFilterBit(nAODtrackCutBit)) continue;\r
-                   \r
-                   v_charge = aodTrackMain->Charge();\r
-                   v_y      = aodTrackMain->Y();\r
-                   v_eta    = aodTrackMain->Eta();\r
-                   v_phi    = aodTrackMain->Phi() * TMath::RadToDeg();\r
-                   v_E      = aodTrackMain->E();\r
-                   v_pt     = aodTrackMain->Pt();\r
-                   aodTrackMain->PxPyPz(v_p);\r
-                   \r
-                   Float_t DCAxy = aodTrackMain->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
-                   Float_t DCAz  = aodTrackMain->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( v_pt < fPtMin || v_pt > fPtMax)      continue;\r
-                   if( v_eta < fEtaMin || v_eta > 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 && aodTrackMain->Chi2perNDF() > fTPCchi2Cut){\r
-                     continue;\r
-                   }\r
-                   if( fNClustersTPCCut != -1 && aodTrackMain->GetTPCNcls() < fNClustersTPCCut){\r
-                     continue;\r
-                   }\r
-                   \r
-                   // fill QA histograms\r
-                   fHistClus->Fill(aodTrackMain->GetITSNcls(),aodTrackMain->GetTPCNcls());\r
-                   fHistDCA->Fill(DCAz,DCAxy);\r
-                   fHistChi2->Fill(aodTrackMain->Chi2perNDF());\r
-                   fHistPt->Fill(v_pt);\r
-                   fHistEta->Fill(v_eta);\r
-                   fHistPhi->Fill(v_phi);\r
-                   \r
-                   // fill charge vector\r
-                   chargeVector[0]->push_back(v_charge);\r
-                   chargeVector[1]->push_back(v_y);\r
-                   chargeVector[2]->push_back(v_eta);\r
-                   chargeVector[3]->push_back(v_phi);\r
-                   chargeVector[4]->push_back(v_p[0]);\r
-                   chargeVector[5]->push_back(v_p[1]);\r
-                   chargeVector[6]->push_back(v_p[2]);\r
-                   chargeVector[7]->push_back(v_pt);\r
-                   chargeVector[8]->push_back(v_E);\r
-                   \r
-                 } //track loop\r
-                 \r
-                   // calculate balance function\r
-                 fBalance->FillBalance(fCentrality,chargeVector);\r
-                 // clean charge vector afterwards\r
-                 for(Int_t i = 0; i < 9; i++){                \r
-                   chargeVector[i]->clear();\r
-                 }\r
-\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
-  }//AOD analysis\r
-  else{\r
-    AliError("Triggered Balance Function analysis only for AODs!");\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
+               // 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
+  Double_t v_charge;\r
+  Double_t v_eta;\r
+  Double_t v_phi;\r
+  Double_t v_pt;\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
+    v_charge = aodTrack->Charge();\r
+    v_eta    = aodTrack->Eta();\r
+    v_phi    = aodTrack->Phi() * TMath::RadToDeg();\r
+    v_pt     = 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( v_pt < fPtMin || v_pt > fPtMax)      continue;\r
+    if( v_eta < fEtaMin || v_eta > 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(v_pt);\r
+    fHistEta->Fill(v_eta);\r
+    fHistPhi->Fill(v_phi);\r
+    \r
+    // add the track to the TObjArray\r
+    tracksAccepted->Add(new AliBFBasicParticle(v_eta, v_phi, v_pt, v_charge));\r
   }\r
-}     \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)));\r
+  }\r
+   \r
+  return tracksShuffled;\r
+}\r
 \r
 //________________________________________________________________________\r
 void  AliAnalysisTaskTriggeredBF::FinishTaskOutput(){\r