Event Mixing for Balance Function Psi Task (PhiCorrelation style)
authormiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Jul 2012 11:46:01 +0000 (11:46 +0000)
committermiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Jul 2012 11:46:01 +0000 (11:46 +0000)
PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx
PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.h
PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx
PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h
PWGCF/EBYE/macros/AddTaskBalancePsiCentralityTrain.C

index 29c7b5e..30e055a 100755 (executable)
@@ -27,7 +27,9 @@
 #include "AliStack.h"\r
 #include "AliESDtrackCuts.h"\r
 #include "AliEventplane.h"\r
-#include "AliTHn.h"              \r
+#include "AliTHn.h"    \r
+\r
+#include "AliEventPoolManager.h"           \r
 \r
 #include "AliPID.h"                \r
 #include "AliPIDResponse.h"        \r
@@ -35,6 +37,7 @@
 \r
 #include "AliAnalysisTaskBFPsi.h"\r
 #include "AliBalancePsi.h"\r
+#include "AliAnalysisTaskTriggeredBF.h"\r
 \r
 \r
 // Analysis task for the BF vs Psi code\r
@@ -51,9 +54,14 @@ AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(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
   fListBF(0),\r
   fListBFS(0),\r
+  fListBFM(0),\r
   fHistListPIDQA(0),\r
   fHistEventStats(0),\r
   fHistCentStats(0),\r
@@ -138,6 +146,7 @@ AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name)
   DefineOutput(2, TList::Class());\r
   DefineOutput(3, TList::Class());\r
   DefineOutput(4, TList::Class());\r
+  DefineOutput(5, TList::Class());\r
 }\r
 \r
 //________________________________________________________________________\r
@@ -182,6 +191,12 @@ void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
       //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
     }\r
   }\r
+  if(fRunMixing) {\r
+    if(!fMixedBalance) {\r
+      fMixedBalance = new AliBalancePsi();\r
+      fMixedBalance->SetAnalysisLevel("ESD");\r
+    }\r
+  }\r
 \r
   //QA list\r
   fList = new TList();\r
@@ -199,6 +214,12 @@ void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
     fListBFS->SetOwner();\r
   }\r
 \r
+  if(fRunMixing) {\r
+    fListBFM = new TList();\r
+    fListBFM->SetName("listTriggeredBFMixed");\r
+    fListBFM->SetOwner();\r
+  }\r
+\r
   //PID QA list\r
   if(fUsePID) {\r
     fHistListPIDQA = new TList();\r
@@ -207,12 +228,12 @@ void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
   }\r
 \r
   //Event stats.\r
-  TString gCutName[4] = {"Total","Offline trigger",\r
-                         "Vertex","Analyzed"};\r
+  TString gCutName[5] = {"Total","Offline trigger",\r
+                         "Vertex","Analyzed","sel. Centrality"};\r
   fHistEventStats = new TH2F("fHistEventStats",\r
                              "Event statistics;;Centrality percentile;N_{events}",\r
-                             4,0.5,4.5,220,-5,105);\r
-  for(Int_t i = 1; i <= 4; i++)\r
+                             5,0.5,5.5,220,-5,105);\r
+  for(Int_t i = 1; i <= 5; i++)\r
     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
   fList->Add(fHistEventStats);\r
 \r
@@ -308,8 +329,33 @@ void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
     fListBFS->Add(fShuffledBalance->GetHistNpp());\r
     fListBFS->Add(fShuffledBalance->GetHistNnp());\r
   }  \r
+\r
+  if(fRunMixing) {\r
+    fListBFM->Add(fMixedBalance->GetHistNp());\r
+    fListBFM->Add(fMixedBalance->GetHistNn());\r
+    fListBFM->Add(fMixedBalance->GetHistNpn());\r
+    fListBFM->Add(fMixedBalance->GetHistNnn());\r
+    fListBFM->Add(fMixedBalance->GetHistNpp());\r
+    fListBFM->Add(fMixedBalance->GetHistNnp());\r
+  }  \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
+\r
+  fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);\r
+\r
   if(fESDtrackCuts) fList->Add(fESDtrackCuts);\r
 \r
   //====================PID========================//\r
@@ -365,279 +411,268 @@ void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {
   PostData(1, fList);\r
   PostData(2, fListBF);\r
   if(fRunShuffling) PostData(3, fListBFS);\r
-  if(fUsePID) PostData(4, fHistListPIDQA);       //PID\r
+  if(fRunMixing) PostData(4, fListBFM);\r
+  if(fUsePID) PostData(5, fHistListPIDQA);       //PID\r
 }\r
 \r
 //________________________________________________________________________\r
 void AliAnalysisTaskBFPsi::UserExec(Option_t *) {\r
   // Main loop\r
   // Called for each event\r
-  TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
-\r
-  AliESDtrack *track_TPC   = NULL;\r
 \r
+  TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
   Int_t gNumberOfAcceptedTracks = 0;\r
-  Float_t fCentrality           = 0.;\r
-  Double_t gReactionPlane       = 0.;\r
-\r
-  // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)\r
-  vector<Double_t> *chargeVectorShuffle[9];   // this will be shuffled\r
-  vector<Double_t> *chargeVector[9];          // original charge\r
-  for(Int_t i = 0; i < 9; i++){\r
-    chargeVectorShuffle[i] = new vector<Double_t>;\r
-    chargeVector[i]        = new vector<Double_t>;\r
-  }\r
+  Double_t fCentrality          = -1.;\r
+  Double_t gReactionPlane       = -1.; \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
+  // get the event (for generator level: MCEvent())\r
+  AliVEvent* eventMain = NULL;\r
+  if(gAnalysisLevel == "MC") {\r
+    eventMain = dynamic_cast<AliVEvent*>(MCEvent()); \r
+  }\r
+  else{\r
+    eventMain = dynamic_cast<AliVEvent*>(InputEvent()); \r
\r
+  }\r
+  if(!eventMain) {\r
+    AliError("eventMain not available");\r
+    return;\r
+  }\r
 \r
+  // PID Response task active?\r
   if(fUsePID) {\r
     fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();\r
     if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
   }\r
\r
-  //ESD analysis\r
-  if(gAnalysisLevel == "ESD") {\r
-    AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE\r
-    if (!gESD) {\r
-      Printf("ERROR: gESD 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
-    // store offline trigger bits\r
-    fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
+  // get the reaction plane\r
+  gReactionPlane = GetEventPlane(eventMain);\r
+  fHistEventPlane->Fill(gReactionPlane,fCentrality);\r
+  if(gReactionPlane < 0){\r
+    return;\r
+  }\r
+  \r
+  // get the accepted tracks in main event\r
+  TObjArray *tracksMain = GetAcceptedTracks(eventMain,fCentrality,gReactionPlane);\r
+  gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();\r
 \r
-    // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
-    fHistEventStats->Fill(1,fCentrality); //all events\r
-    Bool_t isSelected = kTRUE;\r
-    if(fUseOfflineTrigger)\r
-      isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
-    if(isSelected) {\r
-      fHistEventStats->Fill(2,fCentrality); //triggered events\r
+  //multiplicity cut (used in pp)\r
+  fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,fCentrality);\r
+  if(fUseMultiplicity) {\r
+    if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))\r
+      return;\r
+  }\r
 \r
-      if(fUseCentrality) {\r
-       //Centrality stuff\r
-       AliCentrality *centrality = gESD->GetCentrality();\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 \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
-       fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\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->CalculateBalance(gReactionPlane,tracksMain,tracksMixed); \r
+           }\r
+       }\r
        \r
-       // take only events inside centrality class\r
-       if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
-                                                fCentralityPercentileMax,\r
-                                                fCentralityEstimator.Data()))\r
-         return;\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->CalculateBalance(gReactionPlane,tracksMain,NULL);\r
+  \r
+  // calculate shuffled balance function\r
+  if(fRunShuffling && tracksShuffled != NULL) {\r
+    fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL);\r
+  }\r
+}      \r
+\r
+//________________________________________________________________________\r
+Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){\r
+  // Checks the Event cuts\r
+  // Fills Event statistics histograms\r
+  \r
+  Bool_t isSelectedMain = kTRUE;\r
+  Float_t fCentrality = -1.;\r
+  TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
+\r
+  fHistEventStats->Fill(1,fCentrality); //all events\r
+\r
+  // Event trigger bits\r
+  fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
+  if(fUseOfflineTrigger)\r
+    isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
+  \r
+  if(isSelectedMain) {\r
+    fHistEventStats->Fill(2,fCentrality); //triggered events\r
+    \r
+    //Centrality stuff \r
+    if(fUseCentrality) {\r
+      if(gAnalysisLevel == "AOD") { //centrality in AOD header\r
+       AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
+       if(header){\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
+       }//AOD header\r
+      }//AOD\r
+\r
+      else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
+       AliCentrality *centrality = event->GetCentrality();\r
+       fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
+\r
+       // QA for centrality estimators\r
+       fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));\r
+       fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("FMD"));\r
+       fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("TRK"));\r
+       fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("TKL"));\r
+       fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("CL0"));\r
+       fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("CL1"));\r
+       fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("V0MvsFMD"));\r
+       fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("TKLvsV0M"));\r
+       fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZEMvsZDC"));\r
+\r
        // centrality QA (V0M)\r
-       fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
+       fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
+      }//ESD\r
+      else if(gAnalysisLevel == "MC"){\r
+       Double_t gImpactParameter = 0.;\r
+       AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());\r
+       if(headerH){\r
+         gImpactParameter = headerH->ImpactParameter();\r
+         fCentrality      = gImpactParameter;\r
+       }//MC header\r
+      }//MC\r
+      else{\r
+       fCentrality = -1.;\r
       }\r
-       \r
-      const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
+    }\r
+    \r
+    // Event Vertex MC\r
+    if(gAnalysisLevel == "MC"){\r
+      AliGenEventHeader *header = dynamic_cast<AliMCEvent*>(event)->GenEventHeader();\r
+      if(header){  \r
+       TArrayF gVertexArray;\r
+       header->PrimaryVertex(gVertexArray);\r
+       //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r
+       //gVertexArray.At(0),\r
+       //gVertexArray.At(1),\r
+       //gVertexArray.At(2));\r
+       fHistEventStats->Fill(3,fCentrality); //events with a proper vertex\r
+       if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
+         if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
+           if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
+             fHistEventStats->Fill(4,fCentrality); //analayzed events\r
+             fHistVx->Fill(gVertexArray.At(0));\r
+             fHistVy->Fill(gVertexArray.At(1));\r
+             fHistVz->Fill(gVertexArray.At(2),fCentrality);\r
+             \r
+             // take only events inside centrality class\r
+             if((fImpactParameterMin < fCentrality) && (fImpactParameterMax > fCentrality)){\r
+               return fCentrality;         \r
+             }//centrality class\r
+           }//Vz cut\r
+         }//Vy cut\r
+       }//Vx cut\r
+      }//header    \r
+    }//MC\r
+    \r
+    // Event Vertex AOD, ESD, ESDMC\r
+    else{\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(vertex->GetZRes() != 0) {\r
+         if(fCov[5] != 0) {\r
            fHistEventStats->Fill(3,fCentrality); //events with a proper vertex\r
-           if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
-             if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
-               if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
-                 fHistEventStats->Fill(4,fCentrality); //analayzed events\r
-                 fHistVx->Fill(vertex->GetXv());\r
-                 fHistVy->Fill(vertex->GetYv());\r
-                 fHistVz->Fill(vertex->GetZv(),fCentrality);\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,fCentrality); //analyzed events\r
+                 fHistVx->Fill(vertex->GetX());\r
+                 fHistVy->Fill(vertex->GetY());\r
+                 fHistVz->Fill(vertex->GetZ(),fCentrality);\r
                  \r
-                 //========Get the VZERO event plane========//\r
-                 Double_t gVZEROEventPlane = -10.0;\r
-                 Double_t qxTot = 0.0, qyTot = 0.0;\r
-                 AliEventplane *ep = gESD->GetEventplane();\r
-                 if(ep) \r
-                   gVZEROEventPlane = ep->CalculateVZEROEventPlane(gESD,10,2,qxTot,qyTot);\r
-                 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
-                 gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
-                 fHistEventPlane->Fill(gReactionPlane,fCentrality);\r
-                 //========Get the VZERO event plane========//\r
-\r
-                 //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());\r
-                 for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {\r
-                   AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));\r
-                   if (!track) {\r
-                     Printf("ERROR: Could not receive track %d", iTracks);\r
-                     continue;\r
-                   }   \r
-                   \r
-                   // take only TPC only tracks\r
-                   track_TPC   = new AliESDtrack();\r
-                   if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
-                   \r
-                   //ESD track cuts\r
-                   if(fESDtrackCuts) \r
-                     if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
-                   \r
-                   // fill QA histograms\r
-                   Float_t b[2];\r
-                   Float_t bCov[3];\r
-                   track_TPC->GetImpactParameters(b,bCov);\r
-                   if (bCov[0]<=0 || bCov[2]<=0) {\r
-                     AliDebug(1, "Estimated b resolution lower or equal zero!");\r
-                     bCov[0]=0; bCov[2]=0;\r
-                   }\r
-                   \r
-                   Int_t nClustersTPC = -1;\r
-                   nClustersTPC = track_TPC->GetTPCNclsIter1();   // TPC standalone\r
-                   //nClustersTPC = track->GetTPCclusters(0);   // global track\r
-                   Float_t chi2PerClusterTPC = -1;\r
-                   if (nClustersTPC!=0) {\r
-                     chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
-                     //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
-                   }\r
-\r
-                   //===========================PID===============================//               \r
-                   if(fUsePID) {\r
-                     Double_t prob[AliPID::kSPECIES]={0.};\r
-                     Double_t probTPC[AliPID::kSPECIES]={0.};\r
-                     Double_t probTOF[AliPID::kSPECIES]={0.};\r
-                     Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
-\r
-                     Double_t nSigma = 0.;\r
-                      UInt_t detUsedTPC = 0;\r
-                     UInt_t detUsedTOF = 0;\r
-                      UInt_t detUsedTPCTOF = 0;\r
-\r
-                     //Decide what detector configuration we want to use\r
-                     switch(fPidDetectorConfig) {\r
-                     case kTPCpid:\r
-                       fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
-                       nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r
-                       detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
-                       for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
-                         prob[iSpecies] = probTPC[iSpecies];\r
-                       break;\r
-                     case kTOFpid:\r
-                       fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
-                       nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r
-                       detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
-                       for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
-                         prob[iSpecies] = probTOF[iSpecies];\r
-                       break;\r
-                     case kTPCTOF:\r
-                       fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
-                       detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
-                       for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
-                         prob[iSpecies] = probTPCTOF[iSpecies];\r
-                       break;\r
-                     default:\r
-                       break;\r
-                     }//end switch: define detector mask\r
-                     \r
-                     //Filling the PID QA\r
-                     Double_t tofTime = -999., length = 999., tof = -999.;\r
-                     Double_t c = TMath::C()*1.E-9;// m/ns\r
-                     Double_t beta = -999.;\r
-                     Double_t  nSigmaTOFForParticleOfInterest = -999.;\r
-                     if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r
-                          (track->IsOn(AliESDtrack::kTIME))  ) { \r
-                       tofTime = track->GetTOFsignal();//in ps\r
-                       length = track->GetIntegratedLength();\r
-                       tof = tofTime*1E-3; // ns       \r
-                       \r
-                       if (tof <= 0) {\r
-                         //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
-                         continue;\r
-                       }\r
-                       if (length <= 0){\r
-                         //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
-                         continue;\r
-                       }\r
-                       \r
-                       length = length*0.01; // in meters\r
-                       tof = tof*c;\r
-                       beta = length/tof;\r
-                       \r
-                       nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r
-                       fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r
-                       fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
-                       fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
-                     }//TOF signal \r
-                     \r
-                     \r
-                     Double_t  nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r
-                     fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
-                     fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
-                     fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
-                     fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
-                     //end of QA-before pid\r
-                     \r
-                     if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
-                       //Make the decision based on the n-sigma\r
-                       if(fUsePIDnSigma) {\r
-                         if(nSigma > fPIDNSigma) continue;}\r
-                       \r
-                       //Make the decision based on the bayesian\r
-                       else if(fUsePIDPropabilities) {\r
-                         if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
-                         if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;      \r
-                       }\r
-                       \r
-                       //Fill QA after the PID\r
-                       fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r
-                       fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
-                       fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
-                       \r
-                       fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
-                       fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
-                       fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
-                       fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
-                     }\r
-                     \r
-                     PostData(4, fHistListPIDQA);\r
-                   }\r
-                    //===========================PID===============================//\r
-                   v_charge = track_TPC->Charge();\r
-                   v_y      = track_TPC->Y();\r
-                   v_eta    = track_TPC->Eta();\r
-                   v_phi    = track_TPC->Phi() * TMath::RadToDeg();\r
-                   v_E      = track_TPC->E();\r
-                   v_pt     = track_TPC->Pt();\r
-                   track_TPC->PxPyPz(v_p);\r
-                   fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
-                   fHistDCA->Fill(b[1],b[0]);\r
-                   fHistChi2->Fill(chi2PerClusterTPC,fCentrality);\r
-                   fHistPt->Fill(v_pt,fCentrality);\r
-                   fHistEta->Fill(v_eta,fCentrality);\r
-                   fHistPhi->Fill(v_phi,fCentrality);\r
-                   fHistRapidity->Fill(v_y,fCentrality);\r
-                   if(v_charge > 0) fHistPhiPos->Fill(v_phi,fCentrality);\r
-                   else if(v_charge < 0) fHistPhiNeg->Fill(v_phi,fCentrality);\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
-                   if(fRunShuffling) {\r
-                     chargeVectorShuffle[0]->push_back(v_charge);\r
-                     chargeVectorShuffle[1]->push_back(v_y);\r
-                     chargeVectorShuffle[2]->push_back(v_eta);\r
-                     chargeVectorShuffle[3]->push_back(v_phi);\r
-                     chargeVectorShuffle[4]->push_back(v_p[0]);\r
-                     chargeVectorShuffle[5]->push_back(v_p[1]);\r
-                     chargeVectorShuffle[6]->push_back(v_p[2]);\r
-                     chargeVectorShuffle[7]->push_back(v_pt);\r
-                     chargeVectorShuffle[8]->push_back(v_E);\r
-                   }\r
-                   \r
-                   delete track_TPC;\r
-                   \r
-                 } //track loop\r
+                 // take only events inside centrality class\r
+                 if((fCentrality > fCentralityPercentileMin) && (fCentrality < fCentralityPercentileMax)){\r
+                   fHistEventStats->Fill(5,fCentrality); //events with correct centrality\r
+                   return fCentrality;         \r
+                 }//centrality class\r
                }//Vz cut\r
              }//Vy cut\r
            }//Vx cut\r
@@ -645,569 +680,460 @@ void AliAnalysisTaskBFPsi::UserExec(Option_t *) {
        }//proper number of contributors\r
       }//vertex object valid\r
     }//triggered event \r
-  }//ESD analysis\r
+  }//AOD,ESD,ESDMC\r
   \r
-  //AOD analysis (vertex and track cuts also here!!!!)\r
-  else if(gAnalysisLevel == "AOD") {\r
-    AliAODEvent* gAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE\r
-    if(!gAOD) {\r
-      Printf("ERROR: gAOD not available");\r
-      return;\r
-    }\r
-\r
-    AliAODHeader *aodHeader = gAOD->GetHeader();\r
+  // in all other cases return -1 (event not accepted)\r
+  return -1;\r
+}\r
 \r
-    // store offline trigger bits\r
-    fHistTriggerStats->Fill(aodHeader->GetOfflineTrigger());\r
-    \r
-    // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
-    fHistEventStats->Fill(1,fCentrality); //all events\r
-    Bool_t isSelected = kTRUE;\r
-    if(fUseOfflineTrigger)\r
-      isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
-    if(isSelected) {\r
-      fHistEventStats->Fill(2,fCentrality); //triggered events\r
-                 \r
-      //Centrality stuff (centrality in AOD header)\r
-      if(fUseCentrality) {\r
-       fCentrality = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
-       \r
-       // QA for centrality estimators\r
-       fHistCentStats->Fill(0.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0M"));\r
-       fHistCentStats->Fill(1.,aodHeader->GetCentralityP()->GetCentralityPercentile("FMD"));\r
-       fHistCentStats->Fill(2.,aodHeader->GetCentralityP()->GetCentralityPercentile("TRK"));\r
-       fHistCentStats->Fill(3.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKL"));\r
-       fHistCentStats->Fill(4.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL0"));\r
-       fHistCentStats->Fill(5.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL1"));\r
-       fHistCentStats->Fill(6.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
-       fHistCentStats->Fill(7.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
-       fHistCentStats->Fill(8.,aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
-       \r
-       // take only events inside centrality class\r
-       if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax)) \r
-         return;\r
-       \r
-       // centrality QA (V0M)\r
-       fHistV0M->Fill(gAOD->GetVZEROData()->GetMTotV0A(), gAOD->GetVZEROData()->GetMTotV0C());\r
-       \r
-       // centrality QA (reference tracks)\r
-       fHistRefTracks->Fill(0.,aodHeader->GetRefMultiplicity());\r
-       fHistRefTracks->Fill(1.,aodHeader->GetRefMultiplicityPos());\r
-       fHistRefTracks->Fill(2.,aodHeader->GetRefMultiplicityNeg());\r
-       fHistRefTracks->Fill(3.,aodHeader->GetTPConlyRefMultiplicity());\r
-       fHistRefTracks->Fill(4.,aodHeader->GetNumberOfITSClusters(0));\r
-       fHistRefTracks->Fill(5.,aodHeader->GetNumberOfITSClusters(1));\r
-       fHistRefTracks->Fill(6.,aodHeader->GetNumberOfITSClusters(2));\r
-       fHistRefTracks->Fill(7.,aodHeader->GetNumberOfITSClusters(3));\r
-       fHistRefTracks->Fill(8.,aodHeader->GetNumberOfITSClusters(4));\r
-      }\r
+//________________________________________________________________________\r
+Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){\r
+  // Get the event plane\r
 \r
-      const AliAODVertex *vertex = gAOD->GetPrimaryVertex();\r
-      \r
-      if(vertex) {\r
-       Double32_t fCov[6];\r
-       vertex->GetCovarianceMatrix(fCov);\r
-       \r
-       if(vertex->GetNContributors() > 0) {\r
-         if(fCov[5] != 0) {\r
-           fHistEventStats->Fill(3,fCentrality); //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,fCentrality); //analyzed events\r
-                 fHistVx->Fill(vertex->GetX());\r
-                 fHistVy->Fill(vertex->GetY());\r
-                 fHistVz->Fill(vertex->GetZ(),fCentrality);\r
-       \r
-                 //========Get the VZERO event plane========//\r
-                 Double_t gVZEROEventPlane = -10.0;\r
-                 Double_t qxTot = 0.0, qyTot = 0.0;\r
-                 AliEventplane *ep = gAOD->GetEventplane();\r
-                 if(ep) \r
-                   gVZEROEventPlane = ep->CalculateVZEROEventPlane(gAOD,10,2,qxTot,qyTot);\r
-                 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
-                 gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
-                 fHistEventPlane->Fill(gReactionPlane,fCentrality);\r
-                 //========Get the VZERO event plane========//\r
-\r
-                 //Printf("There are %d tracks in this event", gAOD->GetNumberOfTracks());\r
-                 for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {\r
-                   AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));\r
-                   if (!aodTrack) {\r
-                     Printf("ERROR: 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_y      = aodTrack->Y();\r
-                   v_eta    = aodTrack->Eta();\r
-                   v_phi    = aodTrack->Phi() * TMath::RadToDeg();\r
-                   v_E      = aodTrack->E();\r
-                   v_pt     = aodTrack->Pt();\r
-                   aodTrack->PxPyPz(v_p);\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 (!fUsePID) {\r
-                     if( v_eta < fEtaMin || v_eta > fEtaMax)  continue;\r
-                   }\r
-                   else if (fUsePID){\r
-                     if( v_y < fEtaMin || v_y > fEtaMax)  continue;\r
-                   }\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(),fCentrality);\r
-                   fHistPt->Fill(v_pt,fCentrality);\r
-                   fHistEta->Fill(v_eta,fCentrality);\r
-                   fHistPhi->Fill(v_phi,fCentrality);\r
-                   fHistRapidity->Fill(v_y,fCentrality);\r
-                   if(v_charge > 0) fHistPhiPos->Fill(v_phi,fCentrality);\r
-                   else if(v_charge < 0) fHistPhiNeg->Fill(v_phi,fCentrality);\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
-                   if(fRunShuffling) {\r
-                     chargeVectorShuffle[0]->push_back(v_charge);\r
-                     chargeVectorShuffle[1]->push_back(v_y);\r
-                     chargeVectorShuffle[2]->push_back(v_eta);\r
-                     chargeVectorShuffle[3]->push_back(v_phi);\r
-                     chargeVectorShuffle[4]->push_back(v_p[0]);\r
-                     chargeVectorShuffle[5]->push_back(v_p[1]);\r
-                     chargeVectorShuffle[6]->push_back(v_p[2]);\r
-                     chargeVectorShuffle[7]->push_back(v_pt);\r
-                     chargeVectorShuffle[8]->push_back(v_E);\r
-                   }\r
-                                   \r
-                   gNumberOfAcceptedTracks += 1;\r
-                   \r
-                 } //track loop\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
+  TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 \r
-  //MC-ESD analysis\r
-  if(gAnalysisLevel == "MCESD") {\r
-    AliMCEvent*  mcEvent = MCEvent(); \r
-    if (!mcEvent) {\r
-      Printf("ERROR: mcEvent not available");\r
-      return;\r
+  Float_t gVZEROEventPlane    = -10.;\r
+  Float_t gReactionPlane      = -10.;\r
+  Double_t qxTot = 0.0, qyTot = 0.0;\r
+\r
+  //MC: from reaction plane\r
+  if(gAnalysisLevel == "MC"){\r
+   \r
+    AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(dynamic_cast<AliMCEvent*>(event)->GenEventHeader());\r
+    if (headerH) {\r
+      gReactionPlane = headerH->ReactionPlaneAngle();\r
+      gReactionPlane *= TMath::RadToDeg();\r
     }\r
-\r
-    AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE\r
-    if (!gESD) {\r
-      Printf("ERROR: gESD not available");\r
-      return;\r
+  }//MC\r
+  \r
+  // AOD,ESD,ESDMC: from VZERO Event Plane\r
+  else{  \r
+   \r
+    AliEventplane *ep = event->GetEventplane();\r
+    if(ep){ \r
+      gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);\r
+      if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
+      gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
     }\r
+    }//AOD,ESD,ESDMC\r
 \r
-    // store offline trigger bits\r
-    fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
+  return gReactionPlane;\r
+}\r
 \r
-    // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
-    fHistEventStats->Fill(1,fCentrality); //all events\r
-    Bool_t isSelected = kTRUE;\r
-    if(fUseOfflineTrigger)\r
-      isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
-    if(isSelected) {\r
-      fHistEventStats->Fill(2,fCentrality); //triggered events\r
+//________________________________________________________________________\r
+TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t fCentrality, Double_t gReactionPlane){\r
+  // Returns TObjArray with tracks after all track cuts (only for AOD!)\r
+  // Fills QA histograms\r
 \r
-      if(fUseCentrality) {\r
-       //Centrality stuff\r
-       AliCentrality *centrality = gESD->GetCentrality();\r
+  TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 \r
-       fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
-       \r
-       // take only events inside centrality class\r
-       if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
-                                                fCentralityPercentileMax,\r
-                                                fCentralityEstimator.Data()))\r
-         return;\r
-       \r
-       // centrality QA (V0M)\r
-       fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
-      }\r
-       \r
-      const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
-      if(vertex) {\r
-       if(vertex->GetNContributors() > 0) {\r
-         if(vertex->GetZRes() != 0) {\r
-           fHistEventStats->Fill(3,fCentrality); //events with a proper vertex\r
-           if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
-             if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
-               if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
-                 fHistEventStats->Fill(4,fCentrality); //analayzed events\r
-                 fHistVx->Fill(vertex->GetXv());\r
-                 fHistVy->Fill(vertex->GetYv());\r
-                 fHistVz->Fill(vertex->GetZv(),fCentrality);\r
-                 \r
-                 //========Get the VZERO event plane========//\r
-                 Double_t gVZEROEventPlane = -10.0;\r
-                 Double_t qxTot = 0.0, qyTot = 0.0;\r
-                 AliEventplane *ep = gESD->GetEventplane();\r
-                 if(ep) \r
-                   gVZEROEventPlane = ep->CalculateVZEROEventPlane(gESD,10,2,qxTot,qyTot);\r
-                 if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
-                 gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
-                 fHistEventPlane->Fill(gReactionPlane,fCentrality);\r
-                 //========Get the VZERO event plane========//\r
-\r
-                 //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());\r
-                 for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {\r
-                   AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));\r
-                   if (!track) {\r
-                     Printf("ERROR: Could not receive track %d", iTracks);\r
-                     continue;\r
-                   }   \r
-                   \r
-                   Int_t label = TMath::Abs(track->GetLabel());\r
-                   if(label > mcEvent->GetNumberOfTracks()) continue;\r
-                   if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
-                   \r
-                   AliMCParticle* mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
-                   if(!mcTrack) continue;\r
-\r
-                   // take only TPC only tracks\r
-                   track_TPC   = new AliESDtrack();\r
-                   if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
-                   \r
-                   //ESD track cuts\r
-                   if(fESDtrackCuts) \r
-                     if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
-                   \r
-                   // fill QA histograms\r
-                   Float_t b[2];\r
-                   Float_t bCov[3];\r
-                   track_TPC->GetImpactParameters(b,bCov);\r
-                   if (bCov[0]<=0 || bCov[2]<=0) {\r
-                     AliDebug(1, "Estimated b resolution lower or equal zero!");\r
-                     bCov[0]=0; bCov[2]=0;\r
-                   }\r
-                   \r
-                   Int_t nClustersTPC = -1;\r
-                   nClustersTPC = track_TPC->GetTPCNclsIter1();   // TPC standalone\r
-                   //nClustersTPC = track->GetTPCclusters(0);   // global track\r
-                   Float_t chi2PerClusterTPC = -1;\r
-                   if (nClustersTPC!=0) {\r
-                     chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
-                     //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
-                   }\r
-                   \r
-                   v_charge = track_TPC->Charge();\r
-                   v_y      = track_TPC->Y();\r
-                   v_eta    = track_TPC->Eta();\r
-                   v_phi    = track_TPC->Phi() * TMath::RadToDeg();\r
-                   v_E      = track_TPC->E();\r
-                   v_pt     = track_TPC->Pt();\r
-                   track_TPC->PxPyPz(v_p);\r
-\r
-                   fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
-                   fHistDCA->Fill(b[1],b[0]);\r
-                   fHistChi2->Fill(chi2PerClusterTPC,fCentrality);\r
-                   fHistPt->Fill(v_pt,fCentrality);\r
-                   fHistEta->Fill(v_eta,fCentrality);\r
-                   fHistPhi->Fill(v_phi,fCentrality);\r
-                   fHistRapidity->Fill(v_y,fCentrality);\r
-                   if(v_charge > 0) fHistPhiPos->Fill(v_phi,fCentrality);\r
-                   else if(v_charge < 0) fHistPhiNeg->Fill(v_phi,fCentrality);\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
-                   if(fRunShuffling) {\r
-                     chargeVectorShuffle[0]->push_back(v_charge);\r
-                     chargeVectorShuffle[1]->push_back(v_y);\r
-                     chargeVectorShuffle[2]->push_back(v_eta);\r
-                     chargeVectorShuffle[3]->push_back(v_phi);\r
-                     chargeVectorShuffle[4]->push_back(v_p[0]);\r
-                     chargeVectorShuffle[5]->push_back(v_p[1]);\r
-                     chargeVectorShuffle[6]->push_back(v_p[2]);\r
-                     chargeVectorShuffle[7]->push_back(v_pt);\r
-                     chargeVectorShuffle[8]->push_back(v_E);\r
-                   }\r
-                   \r
-                   delete track_TPC;\r
-                   gNumberOfAcceptedTracks += 1;\r
-                   \r
-                 } //track loop\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
-  }//MC-ESD analysis\r
+  //output TObjArray holding all good tracks\r
+  TObjArray* tracksAccepted = new TObjArray;\r
+  tracksAccepted->SetOwner(kTRUE);\r
 \r
-  //MC analysis\r
-  else if(gAnalysisLevel == "MC") {\r
-    AliMCEvent*  mcEvent = MCEvent(); \r
-    if (!mcEvent) {\r
-      Printf("ERROR: mcEvent not available");\r
-      return;\r
-    }\r
-    fHistEventStats->Fill(1,fCentrality); //total events\r
-    fHistEventStats->Fill(2,fCentrality); //offline trigger\r
+  Double_t v_charge;\r
+  Double_t v_eta;\r
+  Double_t v_y;\r
+  Double_t v_phi;\r
+  Double_t v_pt;\r
 \r
-    Double_t gImpactParameter = 0.;\r
-    if(fUseCentrality) {\r
-      //Get the MC header\r
-      AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());\r
-      if (headerH) {\r
-       //Printf("=====================================================");\r
-       //Printf("Reaction plane angle: %lf",headerH->ReactionPlaneAngle());\r
-       //Printf("Impact parameter: %lf",headerH->ImpactParameter());\r
-       //Printf("=====================================================");\r
-       gReactionPlane = headerH->ReactionPlaneAngle();\r
-       gImpactParameter = headerH->ImpactParameter();\r
-       fCentrality = gImpactParameter;\r
+\r
+  if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD\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_y      = aodTrack->Y();\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
-      fCentrality = gImpactParameter;\r
-      fHistEventPlane->Fill(gReactionPlane*TMath::RadToDeg(),fCentrality);\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(),fCentrality);\r
+      fHistPt->Fill(v_pt,fCentrality);\r
+      fHistEta->Fill(v_eta,fCentrality);\r
+      fHistRapidity->Fill(v_y,fCentrality);\r
+      if(v_charge > 0) fHistPhiPos->Fill(v_phi,fCentrality);\r
+      else if(v_charge < 0) fHistPhiNeg->Fill(v_phi,fCentrality);\r
+      fHistPhi->Fill(v_phi,fCentrality);\r
+      \r
+      // add the track to the TObjArray\r
+      tracksAccepted->Add(new AliBFBasicParticle(v_eta, v_phi, v_pt, v_charge));\r
+    }//track loop\r
+  }// AOD analysis\r
 \r
-      // take only events inside centrality class (DIDN'T CHANGE THIS UP TO NOW)\r
-      if((fImpactParameterMin > gImpactParameter) || (fImpactParameterMax < gImpactParameter))\r
-       return;\r
-    }\r
+\r
+  else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD\r
+\r
+    AliESDtrack *track_TPC   = NULL;\r
+    AliMCParticle *mcTrack   = NULL;\r
+\r
+    AliMCEvent*  mcEvent     = NULL;\r
     \r
-    AliGenEventHeader *header = mcEvent->GenEventHeader();\r
-    if(!header) return;\r
+    // for MC ESDs use also MC information\r
+    if(gAnalysisLevel == "MCESD"){\r
+      mcEvent = MCEvent(); \r
+      if (!mcEvent) {\r
+       AliError("mcEvent not available");\r
+       return tracksAccepted;\r
+      }\r
+    }\r
     \r
-    TArrayF gVertexArray;\r
-    header->PrimaryVertex(gVertexArray);\r
-    //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r
-    //gVertexArray.At(0),\r
-    //gVertexArray.At(1),\r
-    //gVertexArray.At(2));\r
-    fHistEventStats->Fill(3,fCentrality); //events with a proper vertex\r
-    if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
-      if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
-       if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
-         fHistEventStats->Fill(4,fCentrality); //analayzed events\r
-         fHistVx->Fill(gVertexArray.At(0));\r
-         fHistVy->Fill(gVertexArray.At(1));\r
-         fHistVz->Fill(gVertexArray.At(2),fCentrality);\r
-         \r
-         Printf("There are %d tracks in this event", mcEvent->GetNumberOfPrimaries());\r
-         for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfPrimaries(); iTracks++) {\r
-           AliMCParticle* track = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iTracks));\r
-           if (!track) {\r
-             Printf("ERROR: Could not receive particle %d", iTracks);\r
-             continue;\r
-           }\r
-           \r
-           //exclude non stable particles\r
-           if(!mcEvent->IsPhysicalPrimary(iTracks)) continue;\r
+    // Loop over tracks in event\r
+    for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
+      AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));\r
+      if (!track) {\r
+       AliError(Form("Could not receive track %d", iTracks));\r
+       continue;\r
+      }        \r
+\r
+      // for MC ESDs use also MC information --> MC track not used further???\r
+      if(gAnalysisLevel == "MCESD"){\r
+       Int_t label = TMath::Abs(track->GetLabel());\r
+       if(label > mcEvent->GetNumberOfTracks()) continue;\r
+       if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
+       \r
+       mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
+       if(!mcTrack) continue;\r
+      }\r
 \r
-           v_eta    = track->Eta();\r
-           v_pt     = track->Pt();\r
-           v_y      = track->Y();\r
+      // take only TPC only tracks\r
+      track_TPC   = new AliESDtrack();\r
+      if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
+      \r
+      //ESD track cuts\r
+      if(fESDtrackCuts) \r
+       if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
+      \r
+      // fill QA histograms\r
+      Float_t b[2];\r
+      Float_t bCov[3];\r
+      track_TPC->GetImpactParameters(b,bCov);\r
+      if (bCov[0]<=0 || bCov[2]<=0) {\r
+       AliDebug(1, "Estimated b resolution lower or equal zero!");\r
+       bCov[0]=0; bCov[2]=0;\r
+      }\r
+      \r
+      Int_t nClustersTPC = -1;\r
+      nClustersTPC = track_TPC->GetTPCNclsIter1();   // TPC standalone\r
+      //nClustersTPC = track->GetTPCclusters(0);   // global track\r
+      Float_t chi2PerClusterTPC = -1;\r
+      if (nClustersTPC!=0) {\r
+       chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
+       //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
+      }\r
+      \r
+      //===========================PID===============================//                    \r
+      if(fUsePID) {\r
+       Double_t prob[AliPID::kSPECIES]={0.};\r
+       Double_t probTPC[AliPID::kSPECIES]={0.};\r
+       Double_t probTOF[AliPID::kSPECIES]={0.};\r
+       Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
+       \r
+       Double_t nSigma = 0.;\r
+       UInt_t detUsedTPC = 0;\r
+       UInt_t detUsedTOF = 0;\r
+       UInt_t detUsedTPCTOF = 0;\r
+       \r
+       //Decide what detector configuration we want to use\r
+       switch(fPidDetectorConfig) {\r
+       case kTPCpid:\r
+         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
+         nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r
+         detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
+         for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
+           prob[iSpecies] = probTPC[iSpecies];\r
+         break;\r
+       case kTOFpid:\r
+         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
+         nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r
+         detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
+         for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
+           prob[iSpecies] = probTOF[iSpecies];\r
+         break;\r
+       case kTPCTOF:\r
+         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
+         detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
+         for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
+           prob[iSpecies] = probTPCTOF[iSpecies];\r
+         break;\r
+       default:\r
+         break;\r
+       }//end switch: define detector mask\r
+       \r
+       //Filling the PID QA\r
+       Double_t tofTime = -999., length = 999., tof = -999.;\r
+       Double_t c = TMath::C()*1.E-9;// m/ns\r
+       Double_t beta = -999.;\r
+       Double_t  nSigmaTOFForParticleOfInterest = -999.;\r
+       if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r
+            (track->IsOn(AliESDtrack::kTIME))  ) { \r
+         tofTime = track->GetTOFsignal();//in ps\r
+         length = track->GetIntegratedLength();\r
+         tof = tofTime*1E-3; // ns     \r
+         \r
+         if (tof <= 0) {\r
+           //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
+           continue;\r
+         }\r
+         if (length <= 0){\r
+           //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
+           continue;\r
+         }\r
+         \r
+         length = length*0.01; // in meters\r
+         tof = tof*c;\r
+         beta = length/tof;\r
+         \r
+         nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r
+         fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r
+         fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
+         fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
+       }//TOF signal \r
+       \r
+       \r
+       Double_t  nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r
+       fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
+       fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
+       fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
+       fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
+       //end of QA-before pid\r
+       \r
+       if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
+         //Make the decision based on the n-sigma\r
+         if(fUsePIDnSigma) {\r
+           if(nSigma > fPIDNSigma) continue;}\r
+         \r
+         //Make the decision based on the bayesian\r
+         else if(fUsePIDPropabilities) {\r
+           if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
+           if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;      \r
+         }\r
+         \r
+         //Fill QA after the PID\r
+         fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r
+         fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
+         fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
+         \r
+         fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
+         fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
+         fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
+         fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
+       }\r
+      }\r
+      //===========================PID===============================//\r
+      v_charge = track_TPC->Charge();\r
+      v_y      = track_TPC->Y();\r
+      v_eta    = track_TPC->Eta();\r
+      v_phi    = track_TPC->Phi() * TMath::RadToDeg();\r
+      v_pt     = track_TPC->Pt();\r
+      fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
+      fHistDCA->Fill(b[1],b[0]);\r
+      fHistChi2->Fill(chi2PerClusterTPC,fCentrality);\r
+      fHistPt->Fill(v_pt,fCentrality);\r
+      fHistEta->Fill(v_eta,fCentrality);\r
+      fHistPhi->Fill(v_phi,fCentrality);\r
+      fHistRapidity->Fill(v_y,fCentrality);\r
+      if(v_charge > 0) fHistPhiPos->Fill(v_phi,fCentrality);\r
+      else if(v_charge < 0) fHistPhiNeg->Fill(v_phi,fCentrality);\r
+      \r
+      // add the track to the TObjArray\r
+      tracksAccepted->Add(new AliBFBasicParticle(v_eta, v_phi, v_pt, v_charge));      \r
 \r
-           if( v_pt < fPtMin || v_pt > fPtMax)      \r
-             continue;\r
-           if (!fUsePID) {\r
-             if( v_eta < fEtaMin || v_eta > fEtaMax)  continue;\r
-           }\r
-           else if (fUsePID){\r
-             if( v_y < fEtaMin || v_y > fEtaMax)  continue;\r
-           }\r
+      delete track_TPC;\r
+    }//track loop\r
+  }// ESD analysis\r
 \r
-           //analyze one set of particles\r
-           if(fUseMCPdgCode) {\r
-             TParticle *particle = track->Particle();\r
-             if(!particle) continue;\r
-             \r
-             Int_t gPdgCode = particle->GetPdgCode();\r
-             if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r
-               continue;\r
-           }\r
-           \r
-           //Use the acceptance parameterization\r
-           if(fAcceptanceParameterization) {\r
-             Double_t gRandomNumber = gRandom->Rndm();\r
-             if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r
-               continue;\r
-           }\r
+  else if(gAnalysisLevel = "MC"){\r
+    \r
+    // Loop over tracks in event\r
+    for (Int_t iTracks = 0; iTracks < dynamic_cast<AliMCEvent*>(event)->GetNumberOfPrimaries(); iTracks++) {\r
+      AliMCParticle* track = dynamic_cast<AliMCParticle *>(event->GetTrack(iTracks));\r
+      if (!track) {\r
+       AliError(Form("Could not receive particle %d", iTracks));\r
+       continue;\r
+      }\r
            \r
-           //Exclude resonances\r
-           if(fExcludeResonancesInMC) {\r
-             TParticle *particle = track->Particle();\r
-             if(!particle) continue;\r
-             \r
-             Bool_t kExcludeParticle = kFALSE;\r
-             Int_t gMotherIndex = particle->GetFirstMother();\r
-             if(gMotherIndex != -1) {\r
-               AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r
-               if(motherTrack) {\r
-                 TParticle *motherParticle = motherTrack->Particle();\r
-                 if(motherParticle) {\r
-                   Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r
-                   //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
-                   if(pdgCodeOfMother == 113) {\r
-                     kExcludeParticle = kTRUE;\r
-                   }\r
-                 }\r
-               }\r
-             }\r
-             \r
-             //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
-             if(kExcludeParticle) continue;\r
-           }\r
+      //exclude non stable particles\r
+      if(!(dynamic_cast<AliMCEvent*>(event)->IsPhysicalPrimary(iTracks))) continue;\r
 \r
-           v_charge = track->Charge();\r
-           v_phi    = track->Phi();\r
-           v_E      = track->E();\r
-           track->PxPyPz(v_p);\r
-           //Printf("phi (before): %lf",v_phi);\r
-\r
-           fHistPt->Fill(v_pt,fCentrality);\r
-           fHistEta->Fill(v_eta,fCentrality);\r
-           fHistPhi->Fill(v_phi*TMath::RadToDeg(),fCentrality);\r
-           fHistRapidity->Fill(v_y,fCentrality);\r
-           if(v_charge > 0) fHistPhiPos->Fill(v_phi*TMath::RadToDeg(),fCentrality);\r
-           else if(v_charge < 0) fHistPhiNeg->Fill(v_phi*TMath::RadToDeg(),fCentrality);\r
-\r
-           //Flow after burner\r
-           if(fUseFlowAfterBurner) {\r
-             Double_t precisionPhi = 0.001;\r
-             Int_t maxNumberOfIterations = 100;\r
-\r
-             Double_t phi0 = v_phi;\r
-             Double_t gV2 = fDifferentialV2->Eval(v_pt);\r
-\r
-             for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r
-               Double_t phiprev = v_phi;\r
-               Double_t fl = v_phi - phi0 + gV2*TMath::Sin(2.*(v_phi - gReactionPlane));\r
-               Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(v_phi - gReactionPlane)); \r
-               v_phi -= fl/fp;\r
-               if (TMath::AreEqualAbs(phiprev,v_phi,precisionPhi)) break;\r
+      v_eta    = track->Eta();\r
+      v_pt     = track->Pt();\r
+      v_y      = track->Y();\r
+      \r
+      if( v_pt < fPtMin || v_pt > fPtMax)      \r
+       continue;\r
+      if (!fUsePID) {\r
+       if( v_eta < fEtaMin || v_eta > fEtaMax)  continue;\r
+      }\r
+      else if (fUsePID){\r
+       if( v_y < fEtaMin || v_y > fEtaMax)  continue;\r
+      }\r
+      \r
+      //analyze one set of particles\r
+      if(fUseMCPdgCode) {\r
+       TParticle *particle = track->Particle();\r
+       if(!particle) continue;\r
+       \r
+       Int_t gPdgCode = particle->GetPdgCode();\r
+       if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r
+         continue;\r
+      }\r
+      \r
+      //Use the acceptance parameterization\r
+      if(fAcceptanceParameterization) {\r
+       Double_t gRandomNumber = gRandom->Rndm();\r
+       if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r
+         continue;\r
+      }\r
+      \r
+      //Exclude resonances\r
+      if(fExcludeResonancesInMC) {\r
+       TParticle *particle = track->Particle();\r
+       if(!particle) continue;\r
+       \r
+       Bool_t kExcludeParticle = kFALSE;\r
+       Int_t gMotherIndex = particle->GetFirstMother();\r
+       if(gMotherIndex != -1) {\r
+         AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));\r
+         if(motherTrack) {\r
+           TParticle *motherParticle = motherTrack->Particle();\r
+           if(motherParticle) {\r
+             Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r
+             //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
+             if(pdgCodeOfMother == 113) {\r
+               kExcludeParticle = kTRUE;\r
              }\r
-             //Printf("phi (after): %lf\n",v_phi);\r
-                     Double_t v_DeltaphiBefore = phi0 - gReactionPlane;\r
-             if(v_DeltaphiBefore < 0) v_DeltaphiBefore += 2*TMath::Pi();\r
-             fHistPhiBefore->Fill(v_DeltaphiBefore,fCentrality);\r
-\r
-             Double_t v_DeltaphiAfter = v_phi - gReactionPlane;\r
-             if(v_DeltaphiAfter < 0) v_DeltaphiAfter += 2*TMath::Pi();\r
-             fHistPhiAfter->Fill(v_DeltaphiAfter,fCentrality);\r
            }\r
-           \r
-           v_phi *= TMath::RadToDeg();\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
-           if(fRunShuffling) {\r
-             chargeVectorShuffle[0]->push_back(v_charge);\r
-             chargeVectorShuffle[1]->push_back(v_y);\r
-             chargeVectorShuffle[2]->push_back(v_eta);\r
-             chargeVectorShuffle[3]->push_back(v_phi);\r
-             chargeVectorShuffle[4]->push_back(v_p[0]);\r
-             chargeVectorShuffle[5]->push_back(v_p[1]);\r
-             chargeVectorShuffle[6]->push_back(v_p[2]);\r
-             chargeVectorShuffle[7]->push_back(v_pt);\r
-             chargeVectorShuffle[8]->push_back(v_E);\r
-           }\r
-           gNumberOfAcceptedTracks += 1;\r
-                   \r
-         } //track loop\r
-         gReactionPlane *= TMath::RadToDeg();\r
-       }//Vz cut\r
-      }//Vy cut\r
-    }//Vx cut\r
-  }//MC analysis\r
-  \r
-  //multiplicity cut (used in pp)\r
-  if(fUseMultiplicity) {\r
-    if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))\r
-      return;\r
-  }\r
-  fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,fCentrality);\r
+         }\r
+       }\r
+       \r
+       //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
+       if(kExcludeParticle) continue;\r
+      }\r
+      \r
+      v_charge = track->Charge();\r
+      v_phi    = track->Phi();\r
+      //Printf("phi (before): %lf",v_phi);\r
+      \r
+      fHistPt->Fill(v_pt,fCentrality);\r
+      fHistEta->Fill(v_eta,fCentrality);\r
+      fHistPhi->Fill(v_phi*TMath::RadToDeg(),fCentrality);\r
+      fHistRapidity->Fill(v_y,fCentrality);\r
+      if(v_charge > 0) fHistPhiPos->Fill(v_phi*TMath::RadToDeg(),fCentrality);\r
+      else if(v_charge < 0) fHistPhiNeg->Fill(v_phi*TMath::RadToDeg(),fCentrality);\r
+      \r
+      //Flow after burner\r
+      if(fUseFlowAfterBurner) {\r
+       Double_t precisionPhi = 0.001;\r
+       Int_t maxNumberOfIterations = 100;\r
+       \r
+       Double_t phi0 = v_phi;\r
+       Double_t gV2 = fDifferentialV2->Eval(v_pt);\r
+       \r
+       for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r
+         Double_t phiprev = v_phi;\r
+         Double_t fl = v_phi - phi0 + gV2*TMath::Sin(2.*(v_phi - gReactionPlane));\r
+         Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(v_phi - gReactionPlane)); \r
+         v_phi -= fl/fp;\r
+         if (TMath::AreEqualAbs(phiprev,v_phi,precisionPhi)) break;\r
+       }\r
+       //Printf("phi (after): %lf\n",v_phi);\r
+       Double_t v_DeltaphiBefore = phi0 - gReactionPlane;\r
+       if(v_DeltaphiBefore < 0) v_DeltaphiBefore += 2*TMath::Pi();\r
+       fHistPhiBefore->Fill(v_DeltaphiBefore,fCentrality);\r
+       \r
+       Double_t v_DeltaphiAfter = v_phi - gReactionPlane;\r
+       if(v_DeltaphiAfter < 0) v_DeltaphiAfter += 2*TMath::Pi();\r
+       fHistPhiAfter->Fill(v_DeltaphiAfter,fCentrality);\r
+      }\r
+      \r
+      v_phi *= TMath::RadToDeg();\r
+          \r
+      tracksAccepted->Add(new AliBFBasicParticle(v_eta, v_phi, v_pt, v_charge));\r
+      \r
+    } //track loop\r
+  }//MC\r
   \r
-  // calculate balance function\r
-  if(fUseMultiplicity) \r
-    //fBalance->CalculateBalance(gNumberOfAcceptedTracks,gReactionPlane,chargeVector);\r
-    fBalance->CalculateBalance(gReactionPlane,chargeVector);\r
-  else                 \r
-    fBalance->CalculateBalance(gReactionPlane,chargeVector);\r
+  return tracksAccepted;  \r
+}\r
+//________________________________________________________________________\r
+TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks){\r
+  // Clones TObjArray and returns it with tracks after shuffling the charges\r
 \r
-  if(fRunShuffling) {\r
-    // shuffle charges\r
-    random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );\r
-    if(fUseMultiplicity) \r
-      fShuffledBalance->CalculateBalance(gReactionPlane,chargeVectorShuffle);\r
-    else                 \r
-      fShuffledBalance->CalculateBalance(gReactionPlane,chargeVectorShuffle);\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
+\r
+  delete chargeVector;\r
+   \r
+  return tracksShuffled;\r
+}\r
+\r
 \r
 //________________________________________________________________________\r
 void  AliAnalysisTaskBFPsi::FinishTaskOutput(){\r
   //Printf("END BF");\r
 \r
   if (!fBalance) {\r
-    Printf("ERROR: fBalance not available");\r
+    AliError("fBalance not available");\r
     return;\r
   }  \r
   if(fRunShuffling) {\r
     if (!fShuffledBalance) {\r
-      Printf("ERROR: fShuffledBalance not available");\r
+      AliError("fShuffledBalance not available");\r
       return;\r
     }\r
   }\r
index 1be5ee9..b1f1e46 100755 (executable)
@@ -11,6 +11,8 @@ class TF1;
 \r
 class AliBalancePsi;\r
 class AliESDtrackCuts;\r
+class AliEventPoolManager;\r
+\r
 \r
 #include "AliAnalysisTaskSE.h"\r
 #include "AliBalancePsi.h"\r
@@ -38,6 +40,11 @@ class AliAnalysisTaskBFPsi : public AliAnalysisTaskSE {
     fRunShuffling = kTRUE;\r
     fShuffledBalance = analysisShuffled;\r
   }\r
+  void SetMixingObject(AliBalancePsi *const analysisMixed) {\r
+    fRunMixing = kTRUE;\r
+    fMixedBalance = analysisMixed;\r
+  }\r
+  void SetMixingTracks(Int_t tracks) { fMixingTracks = tracks; }\r
   void SetAnalysisCutObject(AliESDtrackCuts *const trackCuts) {\r
     fESDtrackCuts = trackCuts;}\r
   void SetVertexDiamond(Double_t vx, Double_t vy, Double_t vz) {\r
@@ -130,12 +137,23 @@ class AliAnalysisTaskBFPsi : public AliAnalysisTaskSE {
     fPidDetectorConfig = detConfig;}\r
 \r
  private:\r
+  Double_t    IsEventAccepted(AliVEvent* event);\r
+  Double_t    GetEventPlane(AliVEvent* event);\r
+  TObjArray* GetAcceptedTracks(AliVEvent* event, Double_t fCentrality, Double_t gReactionPlane);\r
+  TObjArray* GetShuffledTracks(TObjArray* tracks);\r
+\r
   AliBalancePsi *fBalance; //BF object\r
   Bool_t fRunShuffling;//run shuffling or not\r
   AliBalancePsi *fShuffledBalance; //BF object (shuffled)\r
+  Bool_t fRunMixing;//run mixing or not\r
+  Int_t  fMixingTracks;\r
+  AliBalancePsi *fMixedBalance; //TriggeredBF object (mixed)\r
+  AliEventPoolManager*     fPoolMgr;         //! event pool manager\r
+\r
   TList *fList; //fList object\r
   TList *fListBF; //fList object\r
   TList *fListBFS; //fList object\r
+  TList *fListBFM; //fList object\r
   TList *fHistListPIDQA;  //! list of histograms\r
 \r
   TH2F *fHistEventStats; //event stats\r
@@ -241,4 +259,6 @@ class AliAnalysisTaskBFPsi : public AliAnalysisTaskSE {
   ClassDef(AliAnalysisTaskBFPsi, 5); // example of analysis\r
 };\r
 \r
+\r
+\r
 #endif\r
index 4acee47..cf54283 100644 (file)
@@ -220,15 +220,15 @@ void AliBalancePsi::InitHistograms() {
     fHistNN->SetBinLimits(j, dBinsPair[j]);
     fHistNN->SetVarTitle(j, axisTitlePair[j]);
   }
-  Printf("Finished setting up the AliTHn");
+  AliInfo("Finished setting up the AliTHn");
 }
 
 //____________________________________________________________________//
-void AliBalancePsi::CalculateBalance(Double_t gReactionPlane, 
-                                    vector<Double_t> **chargeVector) {
+void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
+                                    TObjArray *particles, 
+                                    TObjArray *particlesMixed ) {
   // Calculates the balance function
   fAnalyzedEvents++;
-  Int_t i = 0 , j = 0;
   
   // Initialize histograms if not done yet
   if(!fHistPN){
@@ -237,153 +237,106 @@ void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
     InitHistograms();
   }
 
-  Double_t trackVarsSingle[nTrackVariablesSingle];
-  Double_t trackVarsPair[nTrackVariablesPair];
-
-  Int_t gNtrack = chargeVector[0]->size();
-  //Printf("(AliBalancePsi) Number of tracks: %d",gNtrack);
-  Double_t gPsiMinusPhi = 0.;
-  Double_t dy = 0., deta = 0.;
-  Double_t qLong = 0., qOut = 0., qSide = 0., qInv = 0.;
-  Double_t dphi = 0.;
-
-  Double_t charge1  = 0;
-  Double_t eta1 = 0., rap1 = 0.;
-  Double_t px1 = 0., py1 = 0., pz1 = 0.;
-  Double_t pt1 = 0.;
-  Double_t energy1 = 0.;
-  Double_t phi1    = 0.;
-
-  Double_t charge2  = 0;
-  Double_t eta2 = 0., rap2 = 0.;
-  Double_t px2 = 0., py2 = 0., pz2 = 0.;
-  Double_t pt2 = 0.;
-  Double_t energy2 = 0.;
-  Double_t phi2    = 0.;
-  Double_t gPsiMinusPhiBin = -10.;
+  Double_t trackVariablesSingle[nTrackVariablesSingle];
+  Double_t trackVariablesPair[nTrackVariablesPair];
+
+  if (!particles){
+    AliWarning("particles TObjArray is NULL pointer --> return");
+    return;
+  }
   
-  for(i = 0; i < gNtrack; i++) {
-    gPsiMinusPhiBin = -10.;
-    charge1 = chargeVector[0]->at(i);
-    rap1    = chargeVector[1]->at(i);
-    eta1    = chargeVector[2]->at(i);
-    phi1    = chargeVector[3]->at(i);
-    px1     = chargeVector[4]->at(i);
-    py1     = chargeVector[5]->at(i);
-    pz1     = chargeVector[6]->at(i);
-    pt1     = chargeVector[7]->at(i);
-    energy1 = chargeVector[8]->at(i);
-    gPsiMinusPhi   = TMath::Abs(phi1 - gReactionPlane);
-    if((gPsiMinusPhi <= 7.5)||
-       ((172.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5)))
-      gPsiMinusPhiBin = 0.0;
-    else if(((37.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5))||
-           ((127.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5))||
-           ((217.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5))||
-           ((307.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5)))
-      gPsiMinusPhiBin = 1.0;
-    else if(((82.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5))||
-           ((262.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5)))
-      gPsiMinusPhiBin = 2.0;
-    else continue;
-    //cout<<"phi: "<<phi1<<" - Psi: "<<gReactionPlane<<" - phi - Psi: "<<gPsiMinusPhi<<" - Bin: "<<gPsiMinusPhiBin<<endl; 
-    //if(gPsiMinusPhi > 180.) gPsiMinusPhi = 360. - gPsiMinusPhi;
-    //if(gPsiMinusPhi < -fPsiInterval/2) gPsiMinusPhi = 360. + gPsiMinusPhi;
-    
-    //trackVarsSingle[0]    =  fCentrality;
-    trackVarsSingle[0]    =  gPsiMinusPhiBin;
-    //trackVarsSingle[2]    =  eta1;
-    //trackVarsSingle[3]    =  rap1;
-    //trackVarsSingle[3]    =  phi1;  
-    trackVarsSingle[1]    =  pt1;  
-
-    //Printf("Track(a) %d - phi-Psi: %lf",i+1,trackVarsSingle[1]);
-    //fill single particle histograms
-    if(charge1 > 0)  
-      fHistP->Fill(trackVarsSingle,0,1.); 
-    else if(charge1 < 0)            
-      fHistN->Fill(trackVarsSingle,0,1.); 
-
-    for(j = 0; j < i; j++) {
-      charge2 = chargeVector[0]->at(j);
-      rap2    = chargeVector[1]->at(j);
-      eta2    = chargeVector[2]->at(j);
-      phi2    = chargeVector[3]->at(j);
-      px2     = chargeVector[4]->at(j);
-      py2     = chargeVector[5]->at(j);
-      pz2     = chargeVector[6]->at(j);
-      pt2     = chargeVector[7]->at(j);
-      energy2 = chargeVector[8]->at(j);
-      //Printf("Track(b) %d - pt: %lf",j+1,pt2);
-      
-      // filling the arrays
-      // RAPIDITY 
-      dy = rap1 - rap2;
-      
-      // Eta
-      deta = eta1 - eta2;
-      
-      //qlong
-      Double_t eTot = energy1 + energy2;
-      Double_t pxTot = px1 + px2;
-      Double_t pyTot = py1 + py2;
-      Double_t pzTot = pz1 + pz2;
-      Double_t q0Tot = energy1 - energy2;
-      Double_t qxTot = px1 - px2;
-      Double_t qyTot = py1 - py2;
-      Double_t qzTot = pz1 - pz2;
-      
-      Double_t eTot2 = eTot*eTot;
-      Double_t pTot2 = pxTot*pxTot + pyTot*pyTot + pzTot*pzTot;
-      Double_t pzTot2 = pzTot*pzTot;
-      
-      Double_t q0Tot2 = q0Tot*q0Tot;
-      Double_t qTot2  = qxTot*qxTot + qyTot*qyTot + qzTot*qzTot;
-      
-      Double_t snn    = eTot2 - pTot2;
-      Double_t ptTot2 = pTot2 - pzTot2 ;
-      Double_t ptTot  = TMath::Sqrt( ptTot2 );
-      
-      qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + ptTot2);
-      
-      //qout
-      qOut = TMath::Sqrt(snn/(snn + ptTot2)) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
+  // define end of particle loops
+  Int_t iMax = particles->GetEntriesFast();
+  Int_t jMax = iMax;
+  if (particlesMixed)
+    jMax = particlesMixed->GetEntriesFast();
+
+  // Eta() is extremely time consuming, therefore cache it for the inner loop here:
+  TObjArray* particlesSecond = (particlesMixed) ? particlesMixed : particles;
+
+  TArrayF secondEta(jMax);
+  TArrayF secondPhi(jMax);
+  TArrayF secondPt(jMax);
+
+  for (Int_t i=0; i<jMax; i++){
+    secondEta[i] = ((AliVParticle*) particlesSecond->At(i))->Eta();
+    secondPhi[i] = ((AliVParticle*) particlesSecond->At(i))->Phi();
+    secondPt[i]  = ((AliVParticle*) particlesSecond->At(i))->Pt();
+  }
+  
+  // 1st particle loop
+  for (Int_t i=0; i<iMax; i++)
+    {
       
-      //qside
-      qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
+      AliVParticle* firstParticle = (AliVParticle*) particles->At(i);
       
-      //qinv
-      qInv = TMath::Sqrt(TMath::Abs(-q0Tot2 + qTot2 ));
+      // some optimization
+      Float_t firstEta = firstParticle->Eta();
+      Float_t firstPhi = firstParticle->Phi();
+      Float_t firstPt  = firstParticle->Pt();
+
+      // Event plane (determine psi bin)
+      Double_t gPsiMinusPhi    =   0.;
+      Double_t gPsiMinusPhiBin = -10.;
+      gPsiMinusPhi   = TMath::Abs(firstPhi - gReactionPlane);
+      if((gPsiMinusPhi <= 7.5)||
+        ((172.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5)))
+       gPsiMinusPhiBin = 0.0;
+      else if(((37.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5))||
+             ((127.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5))||
+             ((217.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5))||
+             ((307.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5)))
+       gPsiMinusPhiBin = 1.0;
+      else if(((82.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5))||
+             ((262.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5)))
+       gPsiMinusPhiBin = 2.0;
+      else continue;
+
       
-      //phi
-      dphi = phi2 - phi1;
-      if(dphi < -180.) dphi += 360.;  //dphi should be between -180 and 180!
-      else if(dphi > 180.) dphi -= 360.; //dphi should be between -180 and 180!
+      Short_t  charge = (Short_t) firstParticle->Charge();
+
+      trackVariablesSingle[0]    =  gPsiMinusPhiBin;
+      trackVariablesSingle[1]    =  firstPt;  
       
-      //trackVarsPair[0]    =  fCentrality;             
-      trackVarsPair[0]    =  gPsiMinusPhiBin;
-      //trackVarsPair[2]    =  eta1;
-      //trackVarsPair[3]    =  rap1;
-      //trackVarsPair[4]    =  phi1;
-      trackVarsPair[1]    =  deta;
-      //trackVarsPair[8]    =  dy;
-      trackVarsPair[2]    =  dphi;
-      trackVarsPair[3]    =  pt1;
-      trackVarsPair[4]    =  pt2;
-      //trackVarsPair[10]   =  qSide;
-      //trackVarsPair[11]   =  qOut;
-      //trackVarsPair[12]   =  qLong;
-      //trackVarsPair[13]   =  qInv;
+      //fill single particle histograms
+      if(charge > 0)      fHistP->Fill(trackVariablesSingle,0,1.); 
+      else if(charge < 0) fHistN->Fill(trackVariablesSingle,0,1.);  
       
-      if( charge1 > 0 && charge2 < 0)  fHistPN->Fill(trackVarsPair,0,1.); 
-      else if( charge1 < 0 && charge2 > 0)  fHistNP->Fill(trackVarsPair,0,1.); 
-      else if( charge1 > 0 && charge2 > 0)  fHistPP->Fill(trackVarsPair,0,1.); 
-      else if( charge1 < 0 && charge2 < 0)  fHistNN->Fill(trackVarsPair,0,1.); 
-      else AliWarning("Wrong charge combination!");
+
       
-    }//end of 2nd particle loop
-  }//end of 1st particle loop
+      // 2nd particle loop (only for j < i for non double counting in the same pT region)
+      // --> SAME pT region for trigger and assoc: NO double counting with this
+      // --> DIFF pT region for trigger and assoc: Missing assoc. particles with j > i to a trigger i 
+      //                          --> can be handled afterwards by using assoc. as trigger as well ?!     
+      for(Int_t j = 0; j < i; j++) {   // or go to full here (everything prepared)?
+       
+       if (particlesMixed && jMax < i)  // if the mixed track number is smaller than the main event one (could be done better if one loops over all tracks)
+         break;
+
+       AliVParticle* secondParticle = (AliVParticle*) particlesSecond->At(j);
+
+       Short_t charge2 = (Short_t) secondParticle->Charge();
+       
+       trackVariablesPair[0]    =  gPsiMinusPhiBin;
+       trackVariablesPair[1]    =  firstEta - secondEta[j];  // delta eta
+       trackVariablesPair[2]    =  firstPhi - secondPhi[j];  // delta phi
+       if (trackVariablesPair[2] > 180.)   // delta phi between -180 and 180 
+         trackVariablesPair[2] -= 360.;
+       if (trackVariablesPair[2] <  - 180.) 
+         trackVariablesPair[2] += 360.;
+       
+       trackVariablesPair[3]    =  firstPt;      // pt trigger
+       trackVariablesPair[4]    =  secondPt[j];  // pt
+       //      trackVariablesPair[5]    =  fCentrality;  // centrality
+       
+       if( charge > 0 && charge2 < 0)  fHistPN->Fill(trackVariablesPair,0,1.); 
+       else if( charge < 0 && charge2 > 0)  fHistNP->Fill(trackVariablesPair,0,1.); 
+       else if( charge > 0 && charge2 > 0)  fHistPP->Fill(trackVariablesPair,0,1.); 
+       else if( charge < 0 && charge2 < 0)  fHistNN->Fill(trackVariablesPair,0,1.); 
+       else AliWarning(Form("Wrong charge combination: charge1 = %d and charge2 = %d",charge,charge2));
+       
+      }//end of 2nd particle loop
+    }//end of 1st particle loop
 }  
 
 //____________________________________________________________________//
@@ -468,7 +421,7 @@ TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin,
   fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax); 
   fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax); 
 
-  Printf("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0));
+  AliInfo(Form("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0)));
 
   // Project into the wanted space (1st: analysis step, 2nd: axis)
   TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
@@ -529,9 +482,9 @@ TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
     return gHist;
   }
 
-  Printf("Entries (test): %lf",(Double_t)(gHistTest->GetEntries()));
-  Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries()));
-  Printf("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries()));
+  AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
+  AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
+  AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
   
   TCanvas *c2 = new TCanvas("c2","");
   c2->cd();
index b45d53b..5e76601 100644 (file)
@@ -62,7 +62,8 @@ class AliBalancePsi : public TObject {
   Int_t GetNumberOfAnalyzedEvent() {return fAnalyzedEvents;}
 
   void CalculateBalance(Double_t gReactionPlane, 
-                       vector<Double_t> **chargeVector);
+                       TObjArray* particles,
+                       TObjArray* particlesMixed);
   
   TH2D   *GetCorrelationFunctionPN(Double_t psiMin, Double_t psiMax);
   TH2D   *GetCorrelationFunctionNP(Double_t psiMin, Double_t psiMax);
index 2b5ce5b..3f8cc01 100644 (file)
@@ -18,6 +18,7 @@ Double_t gMinAcceptedProbability = 0.7;
 AliAnalysisTaskBFPsi *AddTaskBalancePsiCentralityTrain(Double_t centrMin=0.,\r
                                                       Double_t centrMax=100.,\r
                                                       Bool_t gRunShuffling=kFALSE,\r
+                                                      Bool_t gRunMixing=kFALSE,\r
                                                       TString centralityEstimator="V0M",\r
                                                       Double_t vertexZ=10.,\r
                                                       Double_t DCAxy=-1,\r
@@ -58,18 +59,22 @@ AliAnalysisTaskBFPsi *AddTaskBalancePsiCentralityTrain(Double_t centrMin=0.,
   gROOT->LoadMacro("$ALICE_ROOT/PWGCF/EBYE/macros/configBalanceFunctionPsiAnalysis.C");\r
   AliBalancePsi *bf  = 0;  // Balance Function object\r
   AliBalancePsi *bfs = 0;  // shuffled Balance function object\r
+  AliBalancePsi *bfm = 0;  // mixing Balance function object\r
 \r
   if (analysisType=="ESD"){\r
     bf  = GetBalanceFunctionObject("ESD",centralityEstimator,centrMin,centrMax);\r
     if(gRunShuffling) bfs = GetBalanceFunctionObject("ESD",centralityEstimator,centrMin,centrMax,kTRUE);\r
+    if(gRunMixing)    bfm = GetBalanceFunctionObject("ESD",centralityEstimator,centrMin,centrMax,kFALSE);\r
   }\r
   else if (analysisType=="AOD"){\r
     bf  = GetBalanceFunctionObject("AOD",centralityEstimator,centrMin,centrMax);\r
     if(gRunShuffling) bfs = GetBalanceFunctionObject("AOD",centralityEstimator,centrMin,centrMax,kTRUE);\r
+    if(gRunMixing)    bfm = GetBalanceFunctionObject("AOD",centralityEstimator,centrMin,centrMax,kFALSE);\r
   }\r
   else if (analysisType=="MC"){\r
     bf  = GetBalanceFunctionObject("MC",centralityEstimator,centrMin,centrMax);\r
     if(gRunShuffling) bfs = GetBalanceFunctionObject("MC",centralityEstimator,centrMin,centrMax,kTRUE);\r
+    if(gRunMixing)    bfm = GetBalanceFunctionObject("MC",centralityEstimator,centrMin,centrMax,kFALSE);\r
   }\r
   else{\r
     ::Error("AddTaskBF", "analysis type NOT known.");\r
@@ -81,6 +86,10 @@ AliAnalysisTaskBFPsi *AddTaskBalancePsiCentralityTrain(Double_t centrMin=0.,
   AliAnalysisTaskBFPsi *taskBF = new AliAnalysisTaskBFPsi("TaskBFPsi");\r
   taskBF->SetAnalysisObject(bf);\r
   if(gRunShuffling) taskBF->SetShufflingObject(bfs);\r
+  if(gRunMixing){\r
+    taskBF->SetMixingObject(bfm);\r
+    taskBF->SetMixingTracks(50000);\r
+  }\r
 \r
   taskBF->SetCentralityPercentileRange(centrMin,centrMax);\r
   if(analysisType == "ESD") {\r
@@ -137,13 +146,15 @@ AliAnalysisTaskBFPsi *AddTaskBalancePsiCentralityTrain(Double_t centrMin=0.,
   AliAnalysisDataContainer *coutQA = mgr->CreateContainer("listQAPsi", TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
   AliAnalysisDataContainer *coutBF = mgr->CreateContainer("listBFPsi", TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
   if(gRunShuffling) AliAnalysisDataContainer *coutBFS = mgr->CreateContainer("listBFPsiShuffled", TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  if(gRunMixing) AliAnalysisDataContainer *coutBFM = mgr->CreateContainer("listBFPsiMixed", TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
   if(kUsePID) AliAnalysisDataContainer *coutQAPID = mgr->CreateContainer("listQAPIDPsi", TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
 \r
   mgr->ConnectInput(taskBF, 0, mgr->GetCommonInputContainer());\r
   mgr->ConnectOutput(taskBF, 1, coutQA);\r
   mgr->ConnectOutput(taskBF, 2, coutBF);\r
   if(gRunShuffling) mgr->ConnectOutput(taskBF, 3, coutBFS);\r
-  if(kUsePID && analysisType == "ESD") mgr->ConnectOutput(taskBF, 4, coutQAPID);\r
+  if(gRunMixing) mgr->ConnectOutput(taskBF, 4, coutBFM);\r
+  if(kUsePID && analysisType == "ESD") mgr->ConnectOutput(taskBF, 5, coutQAPID);\r
 \r
   return taskBF;\r
 }\r