Adding the code for the bf analysis wrt Psi
authorpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 3 Apr 2012 22:23:09 +0000 (22:23 +0000)
committerpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 3 Apr 2012 22:23:09 +0000 (22:23 +0000)
PWGCF/CMakelibPWGCFebye.pkg
PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx [new file with mode: 0755]
PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.h [new file with mode: 0755]
PWGCF/EBYE/BalanceFunctions/AliBalance.cxx
PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx [new file with mode: 0644]
PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h [new file with mode: 0644]
PWGCF/EBYE/macros/AddTaskBalancePsiCentralityTrain.C [new file with mode: 0644]
PWGCF/EBYE/macros/configBalanceFunctionPsiAnalysis.C [new file with mode: 0644]
PWGCF/EBYE/macros/runBalanceFunctionPsi.C [new file with mode: 0755]
PWGCF/PWGCFebyeLinkDef.h

index ccbd2fd..f1604c1 100644 (file)
@@ -29,7 +29,9 @@
 set ( SRCS   
     EBYE/BalanceFunctions/AliAnalysisTaskToyModel.cxx 
     EBYE/BalanceFunctions/AliAnalysisTaskBF.cxx 
+    EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx 
     EBYE/BalanceFunctions/AliBalance.cxx 
+    EBYE/BalanceFunctions/AliBalancePsi.cxx 
     EBYE/LRC/AliLRCBase.cxx
     EBYE/LRC/AliAnalysisTaskLRC.cxx 
     EBYE/LRC/AliLRCAnalysis.cxx 
diff --git a/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx b/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx
new file mode 100755 (executable)
index 0000000..030e735
--- /dev/null
@@ -0,0 +1,1188 @@
+#include "TChain.h"\r
+#include "TList.h"\r
+#include "TCanvas.h"\r
+#include "TLorentzVector.h"\r
+#include "TGraphErrors.h"\r
+#include "TH1F.h"\r
+#include "TH2F.h"\r
+#include "TH2D.h"                  \r
+#include "TH3D.h"\r
+#include "TArrayF.h"\r
+#include "TF1.h"\r
+#include "TRandom.h"\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+#include "AliAnalysisManager.h"\r
+\r
+#include "AliESDVertex.h"\r
+#include "AliESDEvent.h"\r
+#include "AliESDInputHandler.h"\r
+#include "AliAODEvent.h"\r
+#include "AliAODTrack.h"\r
+#include "AliAODInputHandler.h"\r
+#include "AliGenEventHeader.h"\r
+#include "AliGenHijingEventHeader.h"\r
+#include "AliMCEventHandler.h"\r
+#include "AliMCEvent.h"\r
+#include "AliStack.h"\r
+#include "AliESDtrackCuts.h"\r
+#include "AliEventplane.h"\r
+\r
+#include "AliPID.h"                \r
+#include "AliPIDResponse.h"        \r
+#include "AliPIDCombined.h"        \r
+\r
+#include "AliAnalysisTaskBFPsi.h"\r
+#include "AliBalancePsi.h"\r
+\r
+\r
+// Analysis task for the BF vs Psi code\r
+// Authors: Panos.Christakoglou@nikhef.nl\r
+\r
+ClassImp(AliAnalysisTaskBFPsi)\r
+\r
+//________________________________________________________________________\r
+AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name) \r
+: AliAnalysisTaskSE(name), \r
+  fBalance(0),\r
+  fRunShuffling(kFALSE),\r
+  fShuffledBalance(0),\r
+  fList(0),\r
+  fListBF(0),\r
+  fListBFS(0),\r
+  fHistListPIDQA(0),\r
+  fHistEventStats(0),\r
+  fHistCentStats(0),\r
+  fHistTriggerStats(0),\r
+  fHistTrackStats(0),\r
+  fHistVx(0),\r
+  fHistVy(0),\r
+  fHistVz(0),\r
+  fHistEventPlane(0),\r
+  fHistClus(0),\r
+  fHistDCA(0),\r
+  fHistChi2(0),\r
+  fHistPt(0),\r
+  fHistEta(0),\r
+  fHistPhi(0),\r
+  fHistPhiBefore(0),\r
+  fHistPhiAfter(0),\r
+  fHistV0M(0),\r
+  fHistRefTracks(0),\r
+  fHistdEdxVsPTPCbeforePID(NULL),\r
+  fHistBetavsPTOFbeforePID(NULL), \r
+  fHistProbTPCvsPtbeforePID(NULL), \r
+  fHistProbTOFvsPtbeforePID(NULL), \r
+  fHistProbTPCTOFvsPtbeforePID(NULL),\r
+  fHistNSigmaTPCvsPtbeforePID(NULL), \r
+  fHistNSigmaTOFvsPtbeforePID(NULL), \r
+  fHistdEdxVsPTPCafterPID(NULL),\r
+  fHistBetavsPTOFafterPID(NULL), \r
+  fHistProbTPCvsPtafterPID(NULL), \r
+  fHistProbTOFvsPtafterPID(NULL), \r
+  fHistProbTPCTOFvsPtafterPID(NULL),\r
+  fHistNSigmaTPCvsPtafterPID(NULL), \r
+  fHistNSigmaTOFvsPtafterPID(NULL),  \r
+  fPIDResponse(0x0),\r
+  fPIDCombined(0x0),\r
+  fParticleOfInterest(kPion),\r
+  fPidDetectorConfig(kTPCTOF),\r
+  fUsePID(kFALSE),\r
+  fUsePIDnSigma(kTRUE),\r
+  fUsePIDPropabilities(kFALSE), \r
+  fPIDNSigma(3.),\r
+  fMinAcceptedPIDProbability(0.8),\r
+  fESDtrackCuts(0),\r
+  fCentralityEstimator("V0M"),\r
+  fUseCentrality(kFALSE),\r
+  fCentralityPercentileMin(0.), \r
+  fCentralityPercentileMax(5.),\r
+  fImpactParameterMin(0.),\r
+  fImpactParameterMax(20.),\r
+  fUseMultiplicity(kFALSE),\r
+  fNumberOfAcceptedTracksMin(0),\r
+  fNumberOfAcceptedTracksMax(10000),\r
+  fHistNumberOfAcceptedTracks(0),\r
+  fUseOfflineTrigger(kFALSE),\r
+  fVxMax(0.3),\r
+  fVyMax(0.3),\r
+  fVzMax(10.),\r
+  nAODtrackCutBit(128),\r
+  fPtMin(0.3),\r
+  fPtMax(1.5),\r
+  fEtaMin(-0.8),\r
+  fEtaMax(-0.8),\r
+  fDCAxyCut(-1),\r
+  fDCAzCut(-1),\r
+  fTPCchi2Cut(-1),\r
+  fNClustersTPCCut(-1),\r
+  fAcceptanceParameterization(0),\r
+  fDifferentialV2(0),\r
+  fUseFlowAfterBurner(kFALSE),\r
+  fExcludeResonancesInMC(kFALSE),\r
+  fUseMCPdgCode(kFALSE),\r
+  fPDGCodeToBeAnalyzed(-1) {\r
+  // Constructor\r
+  // Define input and output slots here\r
+  // Input slot #0 works with a TChain\r
+  DefineInput(0, TChain::Class());\r
+  // Output slot #0 writes into a TH1 container\r
+  DefineOutput(1, TList::Class());\r
+  DefineOutput(2, TList::Class());\r
+  DefineOutput(3, TList::Class());\r
+  DefineOutput(4, TList::Class());\r
+}\r
+\r
+//________________________________________________________________________\r
+AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {\r
+\r
+  // delete fBalance; \r
+  // delete fShuffledBalance; \r
+  // delete fList;\r
+  // delete fListBF; \r
+  // delete fListBFS;\r
+\r
+  // delete fHistEventStats; \r
+  // delete fHistTrackStats; \r
+  // delete fHistVx; \r
+  // delete fHistVy; \r
+  // delete fHistVz; \r
+\r
+  // delete fHistClus;\r
+  // delete fHistDCA;\r
+  // delete fHistChi2;\r
+  // delete fHistPt;\r
+  // delete fHistEta;\r
+  // delete fHistPhi;\r
+  // delete fHistV0M;\r
+}\r
+\r
+//________________________________________________________________________\r
+void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {\r
+  // Create histograms\r
+  // Called once\r
+  if(!fBalance) {\r
+    fBalance = new AliBalancePsi();\r
+    fBalance->SetAnalysisLevel("ESD");\r
+    //fBalance->SetNumberOfBins(-1,16);\r
+    fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
+  }\r
+  if(fRunShuffling) {\r
+    if(!fShuffledBalance) {\r
+      fShuffledBalance = new AliBalancePsi();\r
+      fShuffledBalance->SetAnalysisLevel("ESD");\r
+      //fShuffledBalance->SetNumberOfBins(-1,16);\r
+      fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
+    }\r
+  }\r
+\r
+  //QA list\r
+  fList = new TList();\r
+  fList->SetName("listQA");\r
+  fList->SetOwner();\r
+\r
+  //Balance Function list\r
+  fListBF = new TList();\r
+  fListBF->SetName("listBF");\r
+  fListBF->SetOwner();\r
+\r
+  if(fRunShuffling) {\r
+    fListBFS = new TList();\r
+    fListBFS->SetName("listBFShuffled");\r
+    fListBFS->SetOwner();\r
+  }\r
+\r
+  //PID QA list\r
+  if(fUsePID) {\r
+    fHistListPIDQA = new TList();\r
+    fHistListPIDQA->SetName("listQAPID");\r
+    fHistListPIDQA->SetOwner();\r
+  }\r
+\r
+  //Event stats.\r
+  TString gCutName[4] = {"Total","Offline trigger",\r
+                         "Vertex","Analyzed"};\r
+  fHistEventStats = new TH1F("fHistEventStats",\r
+                             "Event statistics;;N_{events}",\r
+                             4,0.5,4.5);\r
+  for(Int_t i = 1; i <= 4; i++)\r
+    fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
+  fList->Add(fHistEventStats);\r
+\r
+  TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r
+  fHistCentStats = new TH2F("fHistCentStats",\r
+                             "Centrality statistics;;Cent percentile",\r
+                           9,-0.5,8.5,220,-5,105);\r
+  for(Int_t i = 1; i <= 9; i++)\r
+    fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
+  fList->Add(fHistCentStats);\r
+\r
+  fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);\r
+  fList->Add(fHistTriggerStats);\r
+\r
+  fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);\r
+  fList->Add(fHistTrackStats);\r
+\r
+  fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);\r
+  fList->Add(fHistNumberOfAcceptedTracks);\r
+\r
+  // Vertex distributions\r
+  fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);\r
+  fList->Add(fHistVx);\r
+  fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);\r
+  fList->Add(fHistVy);\r
+  fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);\r
+  fList->Add(fHistVz);\r
+\r
+  //Event plane\r
+  fHistEventPlane = new TH1F("fHistEventPlane",";#Psi_{2} [rad];Counts",100,0,2.*TMath::Pi());\r
+  fList->Add(fHistEventPlane);\r
+\r
+  // QA histograms\r
+  fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);\r
+  fList->Add(fHistClus);\r
+  fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);\r
+  fList->Add(fHistChi2);\r
+  fHistDCA  = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); \r
+  fList->Add(fHistDCA);\r
+  fHistPt   = new TH1F("fHistPt","p_{T} distribution",200,0,10);\r
+  fList->Add(fHistPt);\r
+  fHistEta  = new TH1F("fHistEta","#eta distribution",200,-2,2);\r
+  fList->Add(fHistEta);\r
+  fHistPhi  = new TH1F("fHistPhi","#phi distribution",200,-20,380);\r
+  fList->Add(fHistPhi);\r
+  fHistPhiBefore  = new TH1F("fHistPhiBefore","#phi distribution",200,0.,2*TMath::Pi());\r
+  fList->Add(fHistPhiBefore);\r
+  fHistPhiAfter  = new TH1F("fHistPhiAfter","#phi distribution",200,0.,2*TMath::Pi());\r
+  fList->Add(fHistPhiAfter);\r
+  fHistV0M  = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);\r
+  fList->Add(fHistV0M);\r
+  TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};\r
+  fHistRefTracks  = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);\r
+  for(Int_t i = 1; i <= 6; i++)\r
+    fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
+  fList->Add(fHistRefTracks);\r
+\r
+  // Balance function histograms\r
+  // Initialize histograms if not done yet\r
+  if(!fBalance->GetHistNp(0)){\r
+    AliWarning("Histograms not yet initialized! --> Will be done now");\r
+    AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
+    fBalance->InitHistograms();\r
+  }\r
+\r
+  if(fRunShuffling) {\r
+    if(!fShuffledBalance->GetHistNp(0)) {\r
+      AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
+      AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
+      fShuffledBalance->InitHistograms();\r
+    }\r
+  }\r
+\r
+  for(Int_t a = 0; a < ANALYSIS_TYPES; a++){\r
+    fListBF->Add(fBalance->GetHistNp(a));\r
+    fListBF->Add(fBalance->GetHistNn(a));\r
+    fListBF->Add(fBalance->GetHistNpn(a));\r
+    fListBF->Add(fBalance->GetHistNnn(a));\r
+    fListBF->Add(fBalance->GetHistNpp(a));\r
+    fListBF->Add(fBalance->GetHistNnp(a));\r
+\r
+    if(fRunShuffling) {\r
+      fListBFS->Add(fShuffledBalance->GetHistNp(a));\r
+      fListBFS->Add(fShuffledBalance->GetHistNn(a));\r
+      fListBFS->Add(fShuffledBalance->GetHistNpn(a));\r
+      fListBFS->Add(fShuffledBalance->GetHistNnn(a));\r
+      fListBFS->Add(fShuffledBalance->GetHistNpp(a));\r
+      fListBFS->Add(fShuffledBalance->GetHistNnp(a));\r
+    }  \r
+  }\r
+\r
+  if(fESDtrackCuts) fList->Add(fESDtrackCuts);\r
+\r
+  //====================PID========================//\r
+  if(fUsePID) {\r
+    fPIDCombined = new AliPIDCombined();\r
+    fPIDCombined->SetDefaultTPCPriors();\r
+\r
+    fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r
+    fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition \r
+    \r
+    fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); \r
+    fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition\r
+    \r
+    fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); \r
+    fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition \r
+    \r
+    fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
+    fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition \r
+\r
+    fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
+    fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition \r
+    \r
+    fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
+    fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r
+    \r
+    fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
+    fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition \r
+    \r
+    fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r
+    fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition \r
+    \r
+    fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); \r
+    fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition \r
+    \r
+    fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
+    fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition \r
+  \r
+    fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000,  -10, 10, 1000, 0, 2); \r
+    fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition  \r
+    \r
+    fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); \r
+    fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition \r
+\r
+    fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
+    fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition  \r
+    \r
+    fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
+    fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition \r
+  }\r
+  //====================PID========================//\r
+\r
+  // Post output data.\r
+  PostData(1, fList);\r
+  PostData(2, fListBF);\r
+  if(fRunShuffling) PostData(3, fListBFS);\r
+  if(fUsePID) PostData(4, 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
+  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
+\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
+\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
+    // store offline trigger bits\r
+    fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
+\r
+    // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
+    fHistEventStats->Fill(1); //all events\r
+    Bool_t isSelected = kTRUE;\r
+    if(fUseOfflineTrigger)\r
+      isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
+    if(isSelected) {\r
+      fHistEventStats->Fill(2); //triggered events\r
+\r
+      if(fUseCentrality) {\r
+       //Centrality stuff\r
+       AliCentrality *centrality = gESD->GetCentrality();\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); //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); //analayzed events\r
+                 fHistVx->Fill(vertex->GetXv());\r
+                 fHistVy->Fill(vertex->GetYv());\r
+                 fHistVz->Fill(vertex->GetZv());\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
+                 fHistEventPlane->Fill(gVZEROEventPlane);\r
+                 gReactionPlane = gVZEROEventPlane;\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);\r
+                   fHistPt->Fill(v_pt);\r
+                   fHistEta->Fill(v_eta);\r
+                   fHistPhi->Fill(v_phi);\r
+\r
+                   // fill charge vector\r
+                   chargeVector[0]->push_back(v_charge);\r
+                   chargeVector[1]->push_back(v_y);\r
+                   chargeVector[2]->push_back(v_eta);\r
+                   chargeVector[3]->push_back(v_phi);\r
+                   chargeVector[4]->push_back(v_p[0]);\r
+                   chargeVector[5]->push_back(v_p[1]);\r
+                   chargeVector[6]->push_back(v_p[2]);\r
+                   chargeVector[7]->push_back(v_pt);\r
+                   chargeVector[8]->push_back(v_E);\r
+\r
+                   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
+               }//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
+  }//ESD analysis\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
+\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); //all events\r
+    Bool_t isSelected = kTRUE;\r
+    if(fUseOfflineTrigger)\r
+      isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
+    if(isSelected) {\r
+      fHistEventStats->Fill(2); //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
+      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); //events with a proper vertex\r
+           if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
+             if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
+               if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
+                 fHistEventStats->Fill(4); //analyzed events\r
+                 fHistVx->Fill(vertex->GetX());\r
+                 fHistVy->Fill(vertex->GetY());\r
+                 fHistVz->Fill(vertex->GetZ());\r
+       \r
+                 //========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
+                 fHistEventPlane->Fill(gVZEROEventPlane);\r
+                 gReactionPlane = gVZEROEventPlane;\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( v_eta < fEtaMin || v_eta > fEtaMax)  continue;\r
+                   \r
+                   // Extra DCA cuts (for systematic studies [!= -1])\r
+                   if( fDCAxyCut != -1 && fDCAxyCut != -1){\r
+                     if(TMath::Sqrt((DCAxy*DCAxy)/(fDCAxyCut*fDCAxyCut)+(DCAz*DCAz)/(fDCAzCut*fDCAzCut)) > 1 ){\r
+                       continue;  // 2D cut\r
+                     }\r
+                   }\r
+                   \r
+                   // Extra TPC cuts (for systematic studies [!= -1])\r
+                   if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
+                     continue;\r
+                   }\r
+                   if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
+                     continue;\r
+                   }\r
+                   \r
+                   // fill QA histograms\r
+                   fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
+                   fHistDCA->Fill(DCAz,DCAxy);\r
+                   fHistChi2->Fill(aodTrack->Chi2perNDF());\r
+                   fHistPt->Fill(v_pt);\r
+                   fHistEta->Fill(v_eta);\r
+                   fHistPhi->Fill(v_phi);\r
+\r
+                   // 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
+\r
+  //MC-ESD analysis\r
+  if(gAnalysisLevel == "MCESD") {\r
+    AliMCEvent*  mcEvent = MCEvent(); \r
+    if (!mcEvent) {\r
+      Printf("ERROR: mcEvent not available");\r
+      return;\r
+    }\r
+\r
+    AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE\r
+    if (!gESD) {\r
+      Printf("ERROR: gESD not available");\r
+      return;\r
+    }\r
+\r
+    // store offline trigger bits\r
+    fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
+\r
+    // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
+    fHistEventStats->Fill(1); //all events\r
+    Bool_t isSelected = kTRUE;\r
+    if(fUseOfflineTrigger)\r
+      isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
+    if(isSelected) {\r
+      fHistEventStats->Fill(2); //triggered events\r
+\r
+      if(fUseCentrality) {\r
+       //Centrality stuff\r
+       AliCentrality *centrality = gESD->GetCentrality();\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); //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); //analayzed events\r
+                 fHistVx->Fill(vertex->GetXv());\r
+                 fHistVy->Fill(vertex->GetYv());\r
+                 fHistVz->Fill(vertex->GetZv());\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
+                 fHistEventPlane->Fill(gVZEROEventPlane);\r
+                 gReactionPlane = gVZEROEventPlane;\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);\r
+                   fHistPt->Fill(v_pt);\r
+                   fHistEta->Fill(v_eta);\r
+                   fHistPhi->Fill(v_phi);\r
+\r
+                   // fill charge vector\r
+                   chargeVector[0]->push_back(v_charge);\r
+                   chargeVector[1]->push_back(v_y);\r
+                   chargeVector[2]->push_back(v_eta);\r
+                   chargeVector[3]->push_back(v_phi);\r
+                   chargeVector[4]->push_back(v_p[0]);\r
+                   chargeVector[5]->push_back(v_p[1]);\r
+                   chargeVector[6]->push_back(v_p[2]);\r
+                   chargeVector[7]->push_back(v_pt);\r
+                   chargeVector[8]->push_back(v_E);\r
+\r
+                   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
+\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); //total events\r
+    fHistEventStats->Fill(2); //offline trigger\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
+      fCentrality = gImpactParameter;\r
+      fHistEventPlane->Fill(gReactionPlane);\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
+    AliGenEventHeader *header = mcEvent->GenEventHeader();\r
+    if(!header) return;\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); //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); //analayzed events\r
+         fHistVx->Fill(gVertexArray.At(0));\r
+         fHistVy->Fill(gVertexArray.At(1));\r
+         fHistVz->Fill(gVertexArray.At(2));\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
+\r
+           v_eta    = track->Eta();\r
+           v_pt     = track->Pt();\r
+           \r
+           if( v_pt < fPtMin || v_pt > fPtMax)      \r
+             continue;\r
+           if( v_eta < fEtaMin || v_eta > fEtaMax)  \r
+             continue;\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 *>(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
+\r
+           v_charge = track->Charge();\r
+           v_y      = track->Y();\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);\r
+           fHistEta->Fill(v_eta);\r
+           fHistPhi->Fill(v_phi);\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);\r
+\r
+             Double_t v_DeltaphiAfter = v_phi - gReactionPlane;\r
+             if(v_DeltaphiAfter < 0) v_DeltaphiAfter += 2*TMath::Pi();\r
+             fHistPhiAfter->Fill(v_DeltaphiAfter);\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
+       }//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);\r
+  \r
+  // calculate balance function\r
+  if(fUseMultiplicity) \r
+    fBalance->CalculateBalance(gNumberOfAcceptedTracks,gReactionPlane,chargeVector);\r
+  else                 \r
+    fBalance->CalculateBalance(fCentrality,gReactionPlane,chargeVector);\r
+\r
+  if(fRunShuffling) {\r
+    // shuffle charges\r
+    random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );\r
+    if(fUseMultiplicity) \r
+      fShuffledBalance->CalculateBalance(gNumberOfAcceptedTracks,gReactionPlane,chargeVectorShuffle);\r
+    else                 \r
+      fShuffledBalance->CalculateBalance(fCentrality,gReactionPlane,chargeVectorShuffle);\r
+  }\r
+}      \r
+\r
+//________________________________________________________________________\r
+void  AliAnalysisTaskBFPsi::FinishTaskOutput(){\r
+  //Printf("END BF");\r
+\r
+  if (!fBalance) {\r
+    Printf("ERROR: fBalance not available");\r
+    return;\r
+  }  \r
+  if(fRunShuffling) {\r
+    if (!fShuffledBalance) {\r
+      Printf("ERROR: fShuffledBalance not available");\r
+      return;\r
+    }\r
+  }\r
+\r
+}\r
+\r
+//________________________________________________________________________\r
+void AliAnalysisTaskBFPsi::Terminate(Option_t *) {\r
+  // Draw result to the screen\r
+  // Called once at the end of the query\r
+\r
+  // not implemented ...\r
+\r
+}\r
diff --git a/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.h b/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.h
new file mode 100755 (executable)
index 0000000..5491da4
--- /dev/null
@@ -0,0 +1,241 @@
+#ifndef ALIANALYSISTASKBFPSI_CXX\r
+#define ALIANALYSISTASKBFPSI_CXX\r
+\r
+// Analysis task for the BF vs Psi code\r
+// Authors: Panos Cristakoglou@cern.ch\r
+\r
+class TList;\r
+class TH1F;\r
+class TH2F;\r
+class TF1;\r
+\r
+class AliBalancePsi;\r
+class AliESDtrackCuts;\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+#include "AliBalancePsi.h"\r
+\r
+#include "AliPID.h"  \r
+#include "AliPIDResponse.h"\r
+#include "AliPIDCombined.h"\r
\r
+\r
+class AliAnalysisTaskBFPsi : public AliAnalysisTaskSE {\r
+ public:\r
+  AliAnalysisTaskBFPsi(const char *name = "AliAnalysisTaskBFPsi");\r
+  virtual ~AliAnalysisTaskBFPsi(); \r
+  \r
+  \r
+  virtual void   UserCreateOutputObjects();\r
+  virtual void   UserExec(Option_t *option);\r
+  virtual void   FinishTaskOutput();\r
+  virtual void   Terminate(Option_t *);\r
+\r
+  void SetAnalysisObject(AliBalancePsi *const analysis) {\r
+    fBalance         = analysis;\r
+    }\r
+  void SetShufflingObject(AliBalancePsi *const analysisShuffled) {\r
+    fRunShuffling = kTRUE;\r
+    fShuffledBalance = analysisShuffled;\r
+  }\r
+  void SetAnalysisCutObject(AliESDtrackCuts *const trackCuts) {\r
+    fESDtrackCuts = trackCuts;}\r
+  void SetVertexDiamond(Double_t vx, Double_t vy, Double_t vz) {\r
+    fVxMax = vx;\r
+    fVyMax = vy;\r
+    fVzMax = vz;\r
+  }\r
+\r
+  //==============AOD analysis==============//\r
+  void SetAODtrackCutBit(Int_t bit){\r
+    nAODtrackCutBit = bit;\r
+  }\r
+\r
+  void SetKinematicsCutsAOD(Double_t ptmin, Double_t ptmax, Double_t etamin, Double_t etamax){\r
+    fPtMin  = ptmin;\r
+    fPtMax  = ptmax;\r
+    fEtaMin = etamin;\r
+    fEtaMax = etamax;\r
+\r
+  }\r
+\r
+  void SetExtraDCACutsAOD(Double_t DCAxy, Double_t DCAz){\r
+    fDCAxyCut  = DCAxy;\r
+    fDCAzCut = DCAz;\r
+  }\r
+\r
+   void SetExtraTPCCutsAOD(Double_t maxTPCchi2, Int_t minNClustersTPC){\r
+    fTPCchi2Cut      = maxTPCchi2;\r
+    fNClustersTPCCut = minNClustersTPC;\r
+  }\r
+\r
+  //==============MC analysis==============//\r
+  void SetKinematicsCutsMC(Double_t ptmin, Double_t ptmax,\r
+                           Double_t etamin, Double_t etamax){\r
+    fPtMin  = ptmin; fPtMax  = ptmax;\r
+    fEtaMin = etamin; fEtaMax = etamax;\r
+  }\r
+  void UseFlowAfterBurner(TF1 *gDifferentialV2) {\r
+    fDifferentialV2 = gDifferentialV2;\r
+    fUseFlowAfterBurner = kTRUE;\r
+  }\r
+  void ExcludeResonancesInMC() {fExcludeResonancesInMC = kTRUE;}\r
+\r
+  void SetPDGCode(Int_t gPdgCode) {\r
+    fUseMCPdgCode = kTRUE;\r
+    fPDGCodeToBeAnalyzed = gPdgCode;\r
+  }\r
+\r
+  //Centrality\r
+  void SetCentralityEstimator(const char* centralityEstimator) {fCentralityEstimator = centralityEstimator;}\r
+  const char* GetCentralityEstimator(void)                     {return fCentralityEstimator;}\r
+  void SetCentralityPercentileRange(Double_t min, Double_t max) { \r
+    fUseCentrality = kTRUE;\r
+    fCentralityPercentileMin=min;\r
+    fCentralityPercentileMax=max;\r
+  }\r
+  void SetImpactParameterRange(Double_t min, Double_t max) { \r
+    fUseCentrality = kTRUE;\r
+    fImpactParameterMin=min;\r
+    fImpactParameterMax=max;\r
+  }\r
+\r
+  //multiplicity\r
+  void SetMultiplicityRange(Int_t min, Int_t max) {\r
+    fUseMultiplicity = kTRUE;\r
+    fNumberOfAcceptedTracksMin = min;\r
+    fNumberOfAcceptedTracksMax = max;}\r
+  \r
+  void UseOfflineTrigger() {fUseOfflineTrigger = kTRUE;}\r
+  \r
+  //Acceptance filter\r
+  void SetAcceptanceParameterization(TF1 *parameterization) {\r
+    fAcceptanceParameterization = parameterization;}\r
+\r
+  //pid\r
+  enum kDetectorUsedForPID { kTPCpid, kTOFpid, kTPCTOF }; // default TPC & TOF pid (via GetTPCpid & GetTOFpid)  \r
+  enum kParticleOfInterest { kMuon, kElectron, kPion, kKaon, kProton };  \r
+\r
+  void SetUseBayesianPID(Double_t gMinProbabilityValue) {\r
+    fUsePID = kTRUE; fUsePIDnSigma = kFALSE; fUsePIDPropabilities = kTRUE;\r
+    fMinAcceptedPIDProbability = gMinProbabilityValue; }\r
+\r
+  void SetUseNSigmaPID(Double_t gMaxNSigma) {\r
+    fUsePID = kTRUE; fUsePIDPropabilities = kFALSE; fUsePIDnSigma = kTRUE;\r
+    fPIDNSigma = gMaxNSigma; }\r
+\r
+  void SetParticleOfInterest(kParticleOfInterest poi) {\r
+    fParticleOfInterest = poi;}\r
+  void SetDetectorUsedForPID(kDetectorUsedForPID detConfig) {\r
+    fPidDetectorConfig = detConfig;}\r
+\r
+ private:\r
+  AliBalancePsi *fBalance; //BF object\r
+  Bool_t fRunShuffling;//run shuffling or not\r
+  AliBalancePsi *fShuffledBalance; //BF object (shuffled)\r
+  TList *fList; //fList object\r
+  TList *fListBF; //fList object\r
+  TList *fListBFS; //fList object\r
+  TList *fHistListPIDQA;  //! list of histograms\r
+\r
+  TH1F *fHistEventStats; //event stats\r
+  TH2F *fHistCentStats; //centrality stats\r
+  TH1F *fHistTriggerStats; //trigger stats\r
+  TH1F *fHistTrackStats; //Track filter bit stats\r
+  TH1F *fHistVx; //x coordinate of the primary vertex\r
+  TH1F *fHistVy; //y coordinate of the primary vertex\r
+  TH1F *fHistVz; //z coordinate of the primary vertex\r
+\r
+  TH1F *fHistEventPlane; //event plane distribution\r
+\r
+  TH2F *fHistClus;//\r
+  TH2F *fHistDCA;//\r
+  TH1F *fHistChi2;//\r
+  TH1F *fHistPt;//\r
+  TH1F *fHistEta;//\r
+  TH1F *fHistPhi;//\r
+  TH1F *fHistPhiBefore;//\r
+  TH1F *fHistPhiAfter;//\r
+  TH2F *fHistV0M;//\r
+  TH2F *fHistRefTracks;//\r
+\r
+  //============PID============//\r
+  TH2D *fHistdEdxVsPTPCbeforePID;//\r
+  TH2D *fHistBetavsPTOFbeforePID;//\r
+  TH2D *fHistProbTPCvsPtbeforePID; //\r
+  TH2D *fHistProbTOFvsPtbeforePID;//\r
+  TH2D *fHistProbTPCTOFvsPtbeforePID;//\r
+  TH2D *fHistNSigmaTPCvsPtbeforePID;//\r
+  TH2D *fHistNSigmaTOFvsPtbeforePID;//\r
+  TH2D *fHistdEdxVsPTPCafterPID;//\r
+  TH2D *fHistBetavsPTOFafterPID;//\r
+  TH2D *fHistProbTPCvsPtafterPID;//\r
+  TH2D *fHistProbTOFvsPtafterPID;//\r
+  TH2D *fHistProbTPCTOFvsPtafterPID;//\r
+  TH2D *fHistNSigmaTPCvsPtafterPID;//\r
+  TH2D *fHistNSigmaTOFvsPtafterPID; //\r
+\r
+  AliPIDResponse *fPIDResponse;     //! PID response object\r
+  AliPIDCombined       *fPIDCombined;     //! combined PID object\r
+  \r
+  kParticleOfInterest  fParticleOfInterest;\r
+  kDetectorUsedForPID   fPidDetectorConfig;\r
+\r
+  Bool_t fUsePID; //\r
+  Bool_t fUsePIDnSigma;//\r
+  Bool_t fUsePIDPropabilities;//\r
+  Double_t fPIDNSigma;//\r
+  Double_t fMinAcceptedPIDProbability;//\r
+  //============PID============//\r
+\r
+  AliESDtrackCuts *fESDtrackCuts; //ESD track cuts\r
+\r
+  TString fCentralityEstimator;      //"V0M","TRK","TKL","ZDC","FMD"\r
+  Bool_t fUseCentrality;//use the centrality (PbPb) or not (pp)\r
+  Double_t fCentralityPercentileMin;//centrality percentile min\r
+  Double_t fCentralityPercentileMax;//centrality percentile max\r
+  Double_t fImpactParameterMin;//impact parameter min (used for MC)\r
+  Double_t fImpactParameterMax;//impact parameter max (used for MC)\r
+\r
+  Bool_t fUseMultiplicity;//use the multiplicity cuts\r
+  Int_t fNumberOfAcceptedTracksMin;//min. number of number of accepted tracks (used for the multiplicity dependence study - pp)\r
+  Int_t fNumberOfAcceptedTracksMax;//max. number of number of accepted tracks (used for the multiplicity dependence study - pp)\r
+  TH1F *fHistNumberOfAcceptedTracks;//hisot to store the number of accepted tracks\r
+\r
+  Bool_t fUseOfflineTrigger;//Usage of the offline trigger selection\r
+\r
+  Double_t fVxMax;//vxmax\r
+  Double_t fVyMax;//vymax\r
+  Double_t fVzMax;//vzmax\r
+\r
+  Int_t nAODtrackCutBit;//track cut bit from track selection (only used for AODs)\r
+\r
+  Double_t fPtMin;//only used for AODs\r
+  Double_t fPtMax;//only used for AODs\r
+  Double_t fEtaMin;//only used for AODs\r
+  Double_t fEtaMax;//only used for AODs\r
+\r
+  Double_t fDCAxyCut;//only used for AODs\r
+  Double_t fDCAzCut;//only used for AODs\r
+\r
+  Double_t fTPCchi2Cut;//only used for AODs\r
+  Int_t fNClustersTPCCut;//only used for AODs\r
+\r
+  TF1 *fAcceptanceParameterization;//acceptance filter used for MC\r
+\r
+  TF1 *fDifferentialV2;//pt-differential v2 (from real data)\r
+  Bool_t fUseFlowAfterBurner;//Usage of a flow after burner\r
+\r
+  Bool_t fExcludeResonancesInMC;//flag to exclude the resonances' decay products from the MC analysis\r
+  Bool_t fUseMCPdgCode; //Boolean to analyze a set of particles in MC\r
+  Int_t fPDGCodeToBeAnalyzed; //Analyze a set of particles in MC\r
+\r
+  \r
+\r
+  AliAnalysisTaskBFPsi(const AliAnalysisTaskBFPsi&); // not implemented\r
+  AliAnalysisTaskBFPsi& operator=(const AliAnalysisTaskBFPsi&); // not implemented\r
+  \r
+  ClassDef(AliAnalysisTaskBFPsi, 5); // example of analysis\r
+};\r
+\r
+#endif\r
index 30b6575..9b0ca0b 100644 (file)
@@ -254,7 +254,7 @@ void AliBalance::CalculateBalance(Float_t fCentrality,vector<Double_t> **chargeV
   //Printf("(AliBalance) Number of tracks: %d",gNtrack);
 
   for(i = 0; i < gNtrack;i++){
-      Short_t charge          = chargeVector[0]->at(i);
+    Double_t charge          = chargeVector[0]->at(i);
       Double_t rapidity       = chargeVector[1]->at(i);
       Double_t pseudorapidity = chargeVector[2]->at(i);
       Double_t phi            = chargeVector[3]->at(i);
@@ -306,14 +306,14 @@ void AliBalance::CalculateBalance(Float_t fCentrality,vector<Double_t> **chargeV
   Double_t qLong = 0., qOut = 0., qSide = 0., qInv = 0.;
   Double_t dphi = 0.;
 
-  Short_t charge1  = 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.;
 
-  Short_t charge2  = 0;
+  Double_t charge2  = 0;
   Double_t eta2 = 0., rap2 = 0.;
   Double_t px2 = 0., py2 = 0., pz2 = 0.;
   Double_t pt2 = 0.;
diff --git a/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx b/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx
new file mode 100644 (file)
index 0000000..3fbc348
--- /dev/null
@@ -0,0 +1,792 @@
+/**************************************************************************
+ * Author: Panos Christakoglou.                                           *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: AliBalancePsi.cxx 54125 2012-01-24 21:07:41Z miweber $ */
+
+//-----------------------------------------------------------------
+//           Balance Function class
+//   This is the class to deal with the Balance Function wrt Psi analysis
+//   Origin: Panos Christakoglou, Nikhef, Panos.Christakoglou@cern.ch
+//-----------------------------------------------------------------
+
+
+//ROOT
+#include <Riostream.h>
+#include <TMath.h>
+#include <TAxis.h>
+#include <TH2D.h>
+#include <TH3D.h>
+#include <TLorentzVector.h>
+#include <TObjArray.h>
+#include <TGraphErrors.h>
+#include <TString.h>
+
+#include "AliVParticle.h"
+#include "AliMCParticle.h"
+#include "AliESDtrack.h"
+#include "AliAODTrack.h"
+
+#include "AliBalancePsi.h"
+
+ClassImp(AliBalancePsi)
+
+//____________________________________________________________________//
+AliBalancePsi::AliBalancePsi() :
+  TObject(), 
+  bShuffle(kFALSE),
+  fAnalysisLevel("ESD"),
+  fAnalyzedEvents(0) ,
+  fCentralityId(0) ,
+  fCentStart(0.),
+  fCentStop(0.),
+  fPsiInterval(15.),
+  fPsiNumberOfBins(24) {
+  // Default constructor
+  for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
+    if(i == 6) {
+      fNumberOfBins[i] = 180;
+      fP1Start[i]      = -360.0;
+      fP1Stop[i]       = 360.0;
+      fP2Start[i]      = -360.0;
+      fP2Stop[i]       = 360.0;
+      fP2Step[i]       = 0.1;
+    }
+    else {
+      fNumberOfBins[i] = 20;
+      fP1Start[i]      = -1.0;
+      fP1Stop[i]       = 1.0;
+      fP2Start[i]      = 0.0;
+      fP2Stop[i]       = 2.0;
+    }
+    fP2Step[i] = TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins[i];
+    fCentStart = 0.;
+    fCentStop  = 0.;
+
+    for(Int_t k = 0; k < MAXIMUM_STEPS_IN_PSI; k++) {
+      fNn[i][k] = 0.0;
+      fNp[i][k] = 0.0;
+
+      for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
+       fNpp[i][k][j] = .0;
+       fNnn[i][k][j] = .0;
+       fNpn[i][k][j] = .0;
+       fNnp[i][k][j] = .0;
+       fB[i][k][j] = 0.0;
+       ferror[i][k][j] = 0.0;
+      }
+    }
+    fHistP[i]  = NULL;
+    fHistN[i]  = NULL;
+    fHistPP[i] = NULL;
+    fHistPN[i] = NULL;
+    fHistNP[i] = NULL;
+    fHistNN[i] = NULL;
+  }
+}
+
+
+//____________________________________________________________________//
+AliBalancePsi::AliBalancePsi(const AliBalancePsi& balance):
+  TObject(balance), bShuffle(balance.bShuffle), 
+  fAnalysisLevel(balance.fAnalysisLevel),
+  fAnalyzedEvents(balance.fAnalyzedEvents), 
+  fCentralityId(balance.fCentralityId),
+  fCentStart(balance.fCentStart),
+  fCentStop(balance.fCentStop),
+  fPsiInterval(balance.fPsiInterval),
+  fPsiNumberOfBins(balance.fPsiNumberOfBins) {
+  //copy constructor
+  for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
+    fP1Start[i]      = balance.fP1Start[i];
+    fP1Stop[i]       = balance.fP1Stop[i];
+    fNumberOfBins[i] = balance.fNumberOfBins[i];
+    fP2Start[i]      = balance.fP2Start[i];
+    fP2Stop[i]       = balance.fP2Stop[i];
+    fP2Step[i]       = balance.fP2Step[i];
+    fCentStart       = balance.fCentStart;
+    fCentStop        = balance.fCentStop; 
+
+    fHistP[i]        = balance.fHistP[i];
+    fHistN[i]        = balance.fHistN[i];
+    fHistPN[i]        = balance.fHistPN[i];
+    fHistNP[i]        = balance.fHistNP[i];
+    fHistPP[i]        = balance.fHistPP[i];
+    fHistNN[i]        = balance.fHistNN[i];
+
+    for(Int_t k = 0; k < MAXIMUM_STEPS_IN_PSI; k++) {
+      fNn[i][k] = balance.fNn[i][k];
+      fNp[i][k] = balance.fNp[i][k];
+      
+      for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
+        fNpp[i][k][j] = balance.fNpp[i][k][j];
+       fNnn[i][k][j] = balance.fNnn[i][k][j];
+       fNpn[i][k][j] = balance.fNpn[i][k][j];
+       fNnp[i][k][j] = balance.fNnp[i][k][j];
+       fB[i][k][j] = balance.fB[i][k][j];
+       ferror[i][k][j] = balance.ferror[i][k][j];
+      } 
+    }
+  }
+}
+
+//____________________________________________________________________//
+AliBalancePsi::~AliBalancePsi() {
+  // Destructor
+  for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
+    delete fHistP[i];
+    delete fHistN[i];
+    delete fHistPN[i];
+    delete fHistNP[i];
+    delete fHistPP[i];
+    delete fHistNN[i];
+  }
+}
+
+//____________________________________________________________________//
+void AliBalancePsi::SetInterval(Int_t iAnalysisType,
+                               Double_t p1Start, Double_t p1Stop,
+                               Int_t ibins, 
+                               Double_t p2Start, Double_t p2Stop,
+                               Double_t psiInterval) {
+  // Sets the analyzed interval. 
+  // Set the same Information for all analyses
+  fPsiInterval = psiInterval;
+  fPsiNumberOfBins = (Int_t)(360./fPsiInterval);
+
+  if(iAnalysisType == -1){             
+    for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
+      fP1Start[i] = p1Start;
+      fP1Stop[i] = p1Stop;
+      fNumberOfBins[i] = ibins;
+      fP2Start[i] = p2Start;
+      fP2Stop[i] = p2Stop;
+      fP2Step[i] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[i];
+    }
+  }
+  // Set the Information for one analysis
+  else if((iAnalysisType > -1) && (iAnalysisType < ANALYSIS_TYPES)) {
+    fP1Start[iAnalysisType] = p1Start;
+    fP1Stop[iAnalysisType] = p1Stop;
+    fNumberOfBins[iAnalysisType] = ibins;
+    fP2Start[iAnalysisType] = p2Start;
+    fP2Stop[iAnalysisType] = p2Stop;
+    fP2Step[iAnalysisType] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[iAnalysisType];
+  }
+  else {
+    AliError("Wrong ANALYSIS number!");
+  }
+}
+
+//____________________________________________________________________//
+void AliBalancePsi::InitHistograms() {
+  //Initialize the histograms
+  TString histName;
+  for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
+    histName = "fHistP"; histName += gBFPsiAnalysisType[iAnalysisType]; 
+    if(bShuffle) histName.Append("_shuffle");
+    if(fCentralityId) histName += fCentralityId.Data();
+    fHistP[iAnalysisType] = new TH3D(histName.Data(),"",(Int_t)(fCentStop-fCentStart),fCentStart,fCentStop,fPsiNumberOfBins,-fPsiInterval/2.,360.-fPsiInterval/2.,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
+
+    histName = "fHistN"; histName += gBFPsiAnalysisType[iAnalysisType]; 
+    if(bShuffle) histName.Append("_shuffle");
+    if(fCentralityId) histName += fCentralityId.Data();
+    fHistN[iAnalysisType] = new TH3D(histName.Data(),"",(Int_t)(fCentStop-fCentStart),fCentStart,fCentStop,fPsiNumberOfBins,-fPsiInterval/2.,360.-fPsiInterval/2.,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
+  
+    histName = "fHistPN"; histName += gBFPsiAnalysisType[iAnalysisType]; 
+    if(bShuffle) histName.Append("_shuffle");
+    if(fCentralityId) histName += fCentralityId.Data();
+    fHistPN[iAnalysisType] = new TH3D(histName.Data(),"",(Int_t)(fCentStop-fCentStart),fCentStart,fCentStop,fPsiNumberOfBins,-fPsiInterval/2.,360.-fPsiInterval/2.,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
+    
+    histName = "fHistNP"; histName += gBFPsiAnalysisType[iAnalysisType]; 
+    if(bShuffle) histName.Append("_shuffle");
+    if(fCentralityId) histName += fCentralityId.Data();
+    fHistNP[iAnalysisType] = new TH3D(histName.Data(),"",(Int_t)(fCentStop-fCentStart),fCentStart,fCentStop,fPsiNumberOfBins,-fPsiInterval/2.,360.-fPsiInterval/2.,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
+    
+    histName = "fHistPP"; histName += gBFPsiAnalysisType[iAnalysisType]; 
+    if(bShuffle) histName.Append("_shuffle");
+    if(fCentralityId) histName += fCentralityId.Data();
+    fHistPP[iAnalysisType] = new TH3D(histName.Data(),"",(Int_t)(fCentStop-fCentStart),fCentStart,fCentStop,fPsiNumberOfBins,-fPsiInterval/2.,360.-fPsiInterval/2.,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
+    
+    histName = "fHistNN"; histName += gBFPsiAnalysisType[iAnalysisType]; 
+    if(bShuffle) histName.Append("_shuffle");
+    if(fCentralityId) histName += fCentralityId.Data();
+    fHistNN[iAnalysisType] = new TH3D(histName.Data(),"",(Int_t)(fCentStop-fCentStart),fCentStart,fCentStop,fPsiNumberOfBins,-fPsiInterval/2.,360.-fPsiInterval/2.,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
+  }
+}
+
+//____________________________________________________________________//
+void AliBalancePsi::PrintAnalysisSettings() {
+  
+  Printf("======================================");
+  Printf("Analysis level: %s",fAnalysisLevel.Data());
+  Printf("======================================");
+  for(Int_t ibin = 0; ibin < ANALYSIS_TYPES; ibin++){
+    Printf("Interval info for variable %d",ibin);
+    Printf("Analyzed interval (min.): %lf",fP2Start[ibin]);
+    Printf("Analyzed interval (max.): %lf",fP2Stop[ibin]);
+    Printf("Number of bins: %d",fNumberOfBins[ibin]);
+    Printf("Step: %lf",fP2Step[ibin]);
+    Printf("          ");
+  }
+  Printf("======================================");
+}
+
+//____________________________________________________________________//
+void AliBalancePsi::CalculateBalance(Float_t fCentrality, 
+                                    Double_t gReactionPlane, 
+                                    vector<Double_t> **chargeVector) {
+  // Calculates the balance function
+  fAnalyzedEvents++;
+  Int_t i = 0 , j = 0;
+  Int_t iBin = 0;
+  
+  // Initialize histograms if not done yet
+  if(!fHistPN[0]){
+    AliWarning("Histograms not yet initialized! --> Will be done now");
+    AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
+    InitHistograms();
+  }
+
+  Int_t gNtrack = chargeVector[0]->size();
+  //Printf("(AliBalancePsi) Number of tracks: %d",gNtrack);
+  Double_t gPsiMinusPhi = 0., gPsiMinusPhiPrime = 0.;
+  Int_t gPsiBin = -10;
+  for(i = 0; i < gNtrack;i++){
+    Double_t charge          = chargeVector[0]->at(i);
+    Double_t rapidity       = chargeVector[1]->at(i);
+    Double_t pseudorapidity = chargeVector[2]->at(i);
+    Double_t phi            = chargeVector[3]->at(i);
+    gPsiMinusPhi   = gReactionPlane - phi;
+    gPsiMinusPhiPrime = gPsiMinusPhi;
+    if(gPsiMinusPhi < -fPsiInterval/2) gPsiMinusPhiPrime = 360. - gPsiMinusPhiPrime;
+    gPsiBin = (Int_t)((gPsiMinusPhiPrime + 360./fPsiNumberOfBins/2.)/(360./fPsiNumberOfBins));
+    //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
+    for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
+      if(iAnalysisType == kEta) {
+       if((pseudorapidity >= fP1Start[iAnalysisType]) && (pseudorapidity <= fP1Stop[iAnalysisType])) {
+         if(charge > 0) {
+           fNp[iAnalysisType][gPsiBin] += 1.;
+           fHistP[iAnalysisType]->Fill(fCentrality,gPsiMinusPhi,pseudorapidity);
+         }//charge > 0
+         if(charge < 0) {
+           fNn[iAnalysisType][gPsiBin] += 1.;
+           fHistN[iAnalysisType]->Fill(fCentrality,gPsiMinusPhi,pseudorapidity);
+         }//charge < 0
+       }//p1 interval check
+      }//analysis type: eta
+      else if(iAnalysisType == kPhi) {
+       if((phi >= fP1Start[iAnalysisType]) && (phi <= fP1Stop[iAnalysisType])) {
+         if(charge > 0) {
+           fNp[iAnalysisType][gPsiBin] += 1.;
+           fHistP[iAnalysisType]->Fill(fCentrality,gPsiMinusPhi,phi);
+         }//charge > 0
+         if(charge < 0) {
+           fNn[iAnalysisType][gPsiBin] += 1.;
+           fHistN[iAnalysisType]->Fill(fCentrality,gPsiMinusPhi,phi);
+         }//charge < 0
+       }//p1 interval check
+      }//analysis type: phi
+      else {
+       if((rapidity >= fP1Start[iAnalysisType]) && (rapidity <= fP1Stop[iAnalysisType])) {
+         if(charge > 0) {
+           fNp[iAnalysisType][gPsiBin] += 1.;
+           fHistP[iAnalysisType]->Fill(fCentrality,gPsiMinusPhi,rapidity);
+         }//charge > 0
+         if(charge < 0) {
+           fNn[iAnalysisType][gPsiBin] += 1.;
+           fHistN[iAnalysisType]->Fill(fCentrality,gPsiMinusPhi,rapidity);
+         }//charge < 0
+       }//p1 interval check
+      }//analysis type: y, qside, qout, qlong, qinv
+    }//analysis type loop
+  }
+  
+  //Printf("Np: %lf - Nn: %lf",fNp[0],fNn[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.;
+  //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
+  for(i = 1; i < gNtrack; i++) {
+
+      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   = gReactionPlane - phi1;
+      phi1 -= 360.;
+      gPsiMinusPhiPrime = gPsiMinusPhi;
+      if(gPsiMinusPhi < -fPsiInterval/2) gPsiMinusPhiPrime = 360. - gPsiMinusPhiPrime;
+      gPsiBin = (Int_t)((gPsiMinusPhiPrime + 360./fPsiNumberOfBins/2.)/(360./fPsiNumberOfBins));
+
+      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(i);
+       energy2 = chargeVector[8]->at(j);
+       phi2 -= 360.;
+       
+       // 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;
+       
+       //qside
+       qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
+       
+       //qinv
+       qInv = TMath::Sqrt(TMath::Abs(-q0Tot2 + qTot2 ));
+       
+       //phi
+       dphi = phi1 - phi2;
+       if(dphi < -180) dphi = 360 - dphi;  //dphi should be between -180 and 180!
+
+       //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
+       if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
+
+         // rapidity
+         if( dy > fP2Start[kRapidity] && dy < fP2Stop[kRapidity]){
+           iBin = Int_t((dy-fP2Start[kRapidity])/fP2Step[kRapidity]);
+           if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
+             
+             if((charge1 > 0)&&(charge2 > 0)) {
+               fNpp[kRapidity][gPsiBin][iBin] += 1.;
+               fHistPP[kRapidity]->Fill(fCentrality,gPsiMinusPhi,dy);
+             }
+             else if((charge1 < 0)&&(charge2 < 0)) {
+               fNnn[kRapidity][gPsiBin][iBin] += 1.;
+               fHistNN[kRapidity]->Fill(fCentrality,gPsiMinusPhi,dy);
+             }
+             else if((charge1 > 0)&&(charge2 < 0)) {
+               fNpn[kRapidity][gPsiBin][iBin] += 1.;
+               fHistPN[kRapidity]->Fill(fCentrality,gPsiMinusPhi,dy);
+             }
+             else if((charge1 < 0)&&(charge2 > 0)) {
+               fNpn[kRapidity][gPsiBin][iBin] += 1.;
+               fHistPN[kRapidity]->Fill(fCentrality,gPsiMinusPhi,dy);
+             }
+           }//BF binning check
+         }//p2 interval check
+       }//p1 interval check
+       
+       // pseudorapidity
+       if((eta1 >= fP1Start[kEta]) && (eta1 <= fP1Stop[kEta]) && (eta2 >= fP1Start[kEta]) && (eta2 <= fP1Stop[kEta])) {
+         if( deta > fP2Start[kEta] && deta < fP2Stop[kEta]){
+           iBin = Int_t((deta-fP2Start[kEta])/fP2Step[kEta]);  
+           if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
+             if((charge1 > 0)&&(charge2 > 0)) {
+               fNpp[kEta][gPsiBin][iBin] += 1.;
+               fHistPP[kEta]->Fill(fCentrality,gPsiMinusPhi,deta);
+             }
+             if((charge1 < 0)&&(charge2 < 0)) {
+               fNnn[kEta][gPsiBin][iBin] += 1.;
+               fHistNN[kEta]->Fill(fCentrality,gPsiMinusPhi,deta);
+             }
+             if((charge1 > 0)&&(charge2 < 0)) {
+               fNpn[kEta][gPsiBin][iBin] += 1.;
+               fHistPN[kEta]->Fill(fCentrality,gPsiMinusPhi,deta);
+             }
+             if((charge1 < 0)&&(charge2 > 0)) {
+               fNpn[kEta][gPsiBin][iBin] += 1.;
+               fHistPN[kEta]->Fill(fCentrality,gPsiMinusPhi,deta);
+             }
+           }//BF binning check
+         }//p2 interval check
+       }//p1 interval check
+       
+       // Qlong, out, side, inv
+       // Check the p1 intervall for rapidity here (like for single tracks above)
+       if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
+         if( qLong > fP2Start[kQlong] && qLong < fP2Stop[kQlong]){
+           iBin = Int_t((qLong-fP2Start[kQlong])/fP2Step[kQlong]);     
+           if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
+             if((charge1 > 0)&&(charge2 > 0)) {
+               fNpp[kQlong][gPsiBin][iBin] += 1.;
+               fHistPP[kQlong]->Fill(fCentrality,gPsiMinusPhi,qLong);
+             }
+             if((charge1 < 0)&&(charge2 < 0)) {
+               fNnn[kQlong][gPsiBin][iBin] += 1.;
+               fHistNN[kQlong]->Fill(fCentrality,gPsiMinusPhi,qLong);
+             }
+             if((charge1 > 0)&&(charge2 < 0)) {
+               fNpn[kQlong][gPsiBin][iBin] += 1.;
+               fHistPN[kQlong]->Fill(fCentrality,gPsiMinusPhi,qLong);
+             }
+             if((charge1 < 0)&&(charge2 > 0)) {
+               fNpn[kQlong][gPsiBin][iBin] += 1.;
+               fHistPN[kQlong]->Fill(fCentrality,gPsiMinusPhi,qLong);
+             }
+           }//BF binning check
+         }//p2 interval check
+       }//p1 interval check
+         
+       if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
+         if( qOut > fP2Start[kQout] && qOut < fP2Stop[kQout]){
+           iBin = Int_t((qOut-fP2Start[kQout])/fP2Step[kQout]);        
+           if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
+             if((charge1 > 0)&&(charge2 > 0)) {
+               fNpp[kQout][gPsiBin][iBin] += 1.;
+               fHistPP[kQout]->Fill(fCentrality,gPsiMinusPhi,qOut);
+                 }
+             if((charge1 < 0)&&(charge2 < 0)) {
+               fNnn[kQout][gPsiBin][iBin] += 1.;
+               fHistNN[kQout]->Fill(fCentrality,gPsiMinusPhi,qOut);
+             }
+             if((charge1 > 0)&&(charge2 < 0)) {
+               fNpn[kQout][gPsiBin][iBin] += 1.;
+               fHistPN[kQout]->Fill(fCentrality,gPsiMinusPhi,qOut);
+             }
+             if((charge1 < 0)&&(charge2 > 0)) {
+               fNpn[kQout][gPsiBin][iBin] += 1.;
+               fHistPN[kQout]->Fill(fCentrality,gPsiMinusPhi,qOut);
+             }
+           }//BF binning check
+         }//p2 interval check
+       }//p1 interval check    
+       
+       if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
+         if( qSide > fP2Start[kQside] && qSide < fP2Stop[kQside]){
+           iBin = Int_t((qSide-fP2Start[kQside])/fP2Step[kQside]);     
+           if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
+             if((charge1 > 0)&&(charge2 > 0)) {
+               fNpp[kQside][gPsiBin][iBin] += 1.;
+               fHistPP[kQside]->Fill(fCentrality,gPsiMinusPhi,qSide);
+             }
+             if((charge1 < 0)&&(charge2 < 0)) {
+               fNnn[kQside][gPsiBin][iBin] += 1.;
+               fHistNN[kQside]->Fill(fCentrality,gPsiMinusPhi,qSide);
+             }
+             if((charge1 > 0)&&(charge2 < 0)) {
+               fNpn[kQside][gPsiBin][iBin] += 1.;
+               fHistPN[kQside]->Fill(fCentrality,gPsiMinusPhi,qSide);
+             }
+             if((charge1 < 0)&&(charge2 > 0)) {
+               fNpn[kQside][gPsiBin][iBin] += 1.;
+               fHistPN[kQside]->Fill(fCentrality,gPsiMinusPhi,qSide);
+                         }
+           }//BF binning check
+         }//p2 interval check
+       }//p1 interval check
+       
+       if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
+         if( qInv > fP2Start[kQinv] && qInv < fP2Stop[kQinv]){
+           iBin = Int_t((qInv-fP2Start[kQinv])/fP2Step[kQinv]);        
+           if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
+             if((charge1 > 0)&&(charge2 > 0)) {
+               fNpp[kQinv][gPsiBin][iBin] += 1.;
+               fHistPP[kQinv]->Fill(fCentrality,gPsiMinusPhi,qInv);
+             }
+             if((charge1 < 0)&&(charge2 < 0)) {
+               fNnn[kQinv][gPsiBin][iBin] += 1.;
+               fHistNN[kQinv]->Fill(fCentrality,gPsiMinusPhi,qInv);
+             }
+             if((charge1 > 0)&&(charge2 < 0)) {
+               fNpn[kQinv][gPsiBin][iBin] += 1.;
+               fHistPN[kQinv]->Fill(fCentrality,gPsiMinusPhi,qInv);
+             }
+             if((charge1 < 0)&&(charge2 > 0)) {
+               fNpn[kQinv][gPsiBin][iBin] += 1.;
+               fHistPN[kQinv]->Fill(fCentrality,gPsiMinusPhi,qInv);
+             }
+           }//BF binning check
+         }//p2 interval check
+       }//p1 interval check
+       
+       // Phi
+       if((phi1 >= fP1Start[kPhi]) && (phi1 <= fP1Stop[kPhi]) && (phi2 >= fP1Start[kPhi]) && (phi2 <= fP1Stop[kPhi])) {
+         if( dphi > fP2Start[kPhi] && dphi < fP2Stop[kPhi]){
+           iBin = Int_t((dphi-fP2Start[kPhi])/fP2Step[kPhi]);  
+           if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
+             if((charge1 > 0)&&(charge2 > 0)) {
+               fNpp[kPhi][gPsiBin][iBin] += 1.;
+               fHistPP[kPhi]->Fill(fCentrality,gPsiMinusPhi,dphi);
+             }
+             if((charge1 < 0)&&(charge2 < 0)) {
+               fNnn[kPhi][gPsiBin][iBin] += 1.;
+               fHistNN[kPhi]->Fill(fCentrality,gPsiMinusPhi,dphi);
+             }
+             if((charge1 > 0)&&(charge2 < 0)) {
+               fNpn[kPhi][gPsiBin][iBin] += 1.;
+               fHistPN[kPhi]->Fill(fCentrality,gPsiMinusPhi,dphi);
+             }
+             if((charge1 < 0)&&(charge2 > 0)) {
+               fNpn[kPhi][gPsiBin][iBin] += 1.;
+               fHistPN[kPhi]->Fill(fCentrality,gPsiMinusPhi,dphi);
+             }
+           }//BF binning check
+         }//p2 interval check
+       }//p1 interval check
+    }//end of 2nd particle loop
+  }//end of 1st particle loop
+  //Printf("Number of analyzed events: %i",fAnalyzedEvents);
+  //Printf("DeltaEta NN[0] = %.0f, PP[0] = %.0f, NP[0] = %.0f, PN[0] = %.0f",fNnn[kEta][0],fNpp[kEta][0],fNnp[kEta][0],fNpn[kEta][0]);
+}  
+
+
+//____________________________________________________________________//
+Double_t AliBalancePsi::GetBalance(Int_t iAnalysisType, Int_t gPsiBin, 
+                                  Int_t p2) {
+  // Returns the value of the balance function in bin p2
+  fB[iAnalysisType][gPsiBin][p2] = 0.5*(((fNpn[iAnalysisType][gPsiBin][p2] - 2.*fNnn[iAnalysisType][gPsiBin][p2])/fNn[iAnalysisType][gPsiBin]) + ((fNpn[iAnalysisType][gPsiBin][p2] - 2.*fNpp[iAnalysisType][gPsiBin][p2])/fNp[iAnalysisType][gPsiBin]))/fP2Step[iAnalysisType];
+  
+  return fB[iAnalysisType][gPsiBin][p2];
+}
+    
+//____________________________________________________________________//
+Double_t AliBalancePsi::GetError(Int_t iAnalysisType, Int_t gPsiBin, 
+                                Int_t p2) {        
+  // Returns the error on the BF value for bin p2
+  // The errors for fNn and fNp are neglected here (0.1 % of total error)
+  ferror[iAnalysisType][gPsiBin][p2] = TMath::Sqrt( Double_t(fNpp[iAnalysisType][gPsiBin][p2])/(Double_t(fNp[iAnalysisType][gPsiBin])*Double_t(fNp[iAnalysisType][gPsiBin])) + 
+                                          Double_t(fNnn[iAnalysisType][gPsiBin][p2])/(Double_t(fNn[iAnalysisType][gPsiBin])*Double_t(fNn[iAnalysisType][gPsiBin])) + 
+                                          Double_t(fNpn[iAnalysisType][gPsiBin][p2])*TMath::Power((0.5/Double_t(fNp[iAnalysisType][gPsiBin]) + 0.5/Double_t(fNn[iAnalysisType][gPsiBin])),2))/fP2Step[iAnalysisType];
+  
+  return ferror[iAnalysisType][gPsiBin][p2];
+}
+//____________________________________________________________________//
+TGraphErrors *AliBalancePsi::DrawBalance(Int_t iAnalysisType, Int_t gPsiBin) {
+  // Draws the BF
+  Double_t x[MAXIMUM_NUMBER_OF_STEPS];
+  Double_t xer[MAXIMUM_NUMBER_OF_STEPS];
+  Double_t b[MAXIMUM_NUMBER_OF_STEPS];
+  Double_t ber[MAXIMUM_NUMBER_OF_STEPS];
+
+  if((fNp[iAnalysisType][gPsiBin] == 0)||(fNn[iAnalysisType][gPsiBin] == 0)) {
+    cerr<<"Couldn't find any particles in the analyzed interval!!!"<<endl;
+    return NULL;
+  }
+  
+  for(Int_t i = 0; i < fNumberOfBins[iAnalysisType]; i++) {
+    b[i] = GetBalance(iAnalysisType,gPsiBin,i);
+    ber[i] = GetError(iAnalysisType,gPsiBin,i);
+    x[i] = fP2Start[iAnalysisType] + fP2Step[iAnalysisType]*i + fP2Step[iAnalysisType]/2;
+    xer[i] = 0.0;
+  }
+  
+  TGraphErrors *gr = new TGraphErrors(fNumberOfBins[iAnalysisType],x,b,xer,ber);
+  gr->GetXaxis()->SetTitleColor(1);
+  if(iAnalysisType==0) {
+    gr->SetTitle("Balance function B(#Delta y)");
+    gr->GetXaxis()->SetTitle("#Delta y");
+    gr->GetYaxis()->SetTitle("B(#Delta y)");
+  }
+  if(iAnalysisType==1) {
+    gr->SetTitle("Balance function B(#Delta #eta)");
+    gr->GetXaxis()->SetTitle("#Delta #eta");
+    gr->GetYaxis()->SetTitle("B(#Delta #eta)");
+  }
+  if(iAnalysisType==2) {
+    gr->SetTitle("Balance function B(q_{long})");
+    gr->GetXaxis()->SetTitle("q_{long} (GeV/c)");
+    gr->GetYaxis()->SetTitle("B(q_{long}) ((GeV/c)^{-1})");
+  }
+  if(iAnalysisType==3) {
+    gr->SetTitle("Balance function B(q_{out})");
+    gr->GetXaxis()->SetTitle("q_{out} (GeV/c)");
+    gr->GetYaxis()->SetTitle("B(q_{out}) ((GeV/c)^{-1})");
+  }
+  if(iAnalysisType==4) {
+    gr->SetTitle("Balance function B(q_{side})");
+    gr->GetXaxis()->SetTitle("q_{side} (GeV/c)");
+    gr->GetYaxis()->SetTitle("B(q_{side}) ((GeV/c)^{-1})");
+  }
+  if(iAnalysisType==5) {
+    gr->SetTitle("Balance function B(q_{inv})");
+    gr->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
+    gr->GetYaxis()->SetTitle("B(q_{inv}) ((GeV/c)^{-1})");
+  }
+  if(iAnalysisType==6) {
+    gr->SetTitle("Balance function B(#Delta #phi)");
+    gr->GetXaxis()->SetTitle("#Delta #phi");
+    gr->GetYaxis()->SetTitle("B(#Delta #phi)");
+  }
+
+  return gr;
+}
+
+//____________________________________________________________________//
+void AliBalancePsi::PrintResults(Int_t iAnalysisType, TH1D *gHistBalance) {
+  //Prints the calculated width of the BF and its error
+  Double_t x[MAXIMUM_NUMBER_OF_STEPS];
+  Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
+  Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
+  Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
+  Double_t deltaBalP2 = 0.0, integral = 0.0;
+  Double_t deltaErrorNew = 0.0;
+  
+  cout<<"=================================================="<<endl;
+  for(Int_t i = 1; i <= fNumberOfBins[iAnalysisType]; i++) { 
+    x[i-1] = fP2Start[iAnalysisType] + fP2Step[iAnalysisType]*i + fP2Step[iAnalysisType]/2;
+    //cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
+  } 
+  //cout<<"=================================================="<<endl;
+  for(Int_t i = 2; i <= fNumberOfBins[iAnalysisType]; i++) {
+    gSumXi += gHistBalance->GetBinCenter(i);
+    gSumBi += gHistBalance->GetBinContent(i);
+    gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
+    gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
+    gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
+    gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
+    gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
+    
+    deltaBalP2 += fP2Step[iAnalysisType]*TMath::Power(gHistBalance->GetBinError(i),2);
+    integral += fP2Step[iAnalysisType]*gHistBalance->GetBinContent(i);
+  }
+  for(Int_t i = 1; i < fNumberOfBins[iAnalysisType]; i++)
+    deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
+  
+  Double_t integralError = TMath::Sqrt(deltaBalP2);
+  
+  Double_t delta = gSumBiXi / gSumBi;
+  Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
+  cout<<"Analysis type: "<<gBFPsiAnalysisType[iAnalysisType].Data()<<endl;
+  cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
+  cout<<"New error: "<<deltaErrorNew<<endl;
+  cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
+  cout<<"=================================================="<<endl;
+}
+//____________________________________________________________________//
+TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t psiMin, Double_t psiMax) {
+  //Returns the BF histogram, extracted from the 6 TH2D objects 
+  //(private members) of the AliBalancePsi class.
+  //
+  TString gAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"};
+  TString histName = "gHistBalanceFunctionHistogram";
+  histName += gAnalysisType[iAnalysisType];
+
+  SetInterval(iAnalysisType, fHistP[iAnalysisType]->GetZaxis()->GetXmin(),
+             fHistP[iAnalysisType]->GetZaxis()->GetXmin(),
+             fHistPP[iAnalysisType]->GetNbinsZ(),
+             fHistPP[iAnalysisType]->GetZaxis()->GetXmin(),
+             fHistPP[iAnalysisType]->GetZaxis()->GetXmax(),
+             psiMax-psiMin);
+
+  // determine the projection thresholds
+  Int_t binMinX, binMinY, binMinZ;
+  Int_t binMaxX, binMaxY, binMaxZ;
+
+  fHistPP[iAnalysisType]->GetBinXYZ(fHistPP[iAnalysisType]->FindBin(centrMin,psiMin),binMinX,binMinY,binMinZ);
+  fHistPP[iAnalysisType]->GetBinXYZ(fHistPP[iAnalysisType]->FindBin(centrMax,psiMax),binMaxX,binMaxY,binMaxZ);
+
+  TH1D *gHistBalanceFunctionHistogram = new TH1D(histName.Data(),"",fHistPP[iAnalysisType]->GetNbinsZ(),fHistPP[iAnalysisType]->GetZaxis()->GetXmin(),fHistPP[iAnalysisType]->GetZaxis()->GetXmax());
+  switch(iAnalysisType) {
+  case kRapidity:
+    gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta y");
+    gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta y)");
+    break;
+  case kEta:
+    gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #eta");
+    gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #eta)");
+    break;
+  case kQlong:
+    gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{long} (GeV/c)");
+    gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{long})");
+    break;
+  case kQout:
+    gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{out} (GeV/c)");
+    gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{out})");
+    break;
+  case kQside:
+    gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{side} (GeV/c)");
+    gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{side})");
+    break;
+  case kQinv:
+    gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
+    gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{inv})");
+    break;
+  case kPhi:
+    gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #phi (deg.)");
+    gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #phi)");
+    break;
+  default:
+    break;
+  }
+
+  TH1D *hTemp1 = dynamic_cast<TH1D *>(fHistPN[iAnalysisType]->ProjectionZ(Form("%s_Cent_%.0f_%.0f_Psi_%.0f_%.0f",fHistPN[iAnalysisType]->GetName(),centrMin,centrMax,psiMin,psiMax),binMinX,binMaxX,binMinY,binMaxY));
+  TH1D *hTemp2 = dynamic_cast<TH1D *>(fHistPN[iAnalysisType]->ProjectionZ(Form("%s_Cent_%.0f_%.0f_Psi_%.0f_%.0f_copy",fHistPN[iAnalysisType]->GetName(),centrMin,centrMax,psiMin,psiMax),binMinX,binMaxX,binMinY,binMaxY));
+  TH1D *hTemp3 = dynamic_cast<TH1D *>(fHistNN[iAnalysisType]->ProjectionZ(Form("%s_Cent_%.0f_%.0f_Psi_%.0f_%.0f",fHistNN[iAnalysisType]->GetName(),centrMin,centrMax,psiMin,psiMax),binMinX,binMaxX,binMinY,binMaxY));
+  TH1D *hTemp4 = dynamic_cast<TH1D *>(fHistPP[iAnalysisType]->ProjectionZ(Form("%s_Cent_%.0f_%.0f_Psi_%.0f_%.0f",fHistPP[iAnalysisType]->GetName(),centrMin,centrMax,psiMin,psiMax),binMinX,binMaxX,binMinY,binMaxY));
+  TH1D *hTemp5 = dynamic_cast<TH1D *>(fHistN[iAnalysisType]->ProjectionZ(Form("%s_Cent_%.0f_%.0f_Psi_%.0f_%.0f",fHistN[iAnalysisType]->GetName(),centrMin,centrMax,psiMin,psiMax),binMinX,binMaxX,binMinY,binMaxY));
+  TH1D *hTemp6 = dynamic_cast<TH1D *>(fHistP[iAnalysisType]->ProjectionZ(Form("%s_Cent_%.0f_%.0f_Psi_%.0f_%.0f",fHistP[iAnalysisType]->GetName(),centrMin,centrMax,psiMin,psiMax),binMinX,binMaxX,binMinY,binMaxY));
+
+  if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)) {
+    hTemp1->Sumw2();
+    hTemp2->Sumw2();
+    hTemp3->Sumw2();
+    hTemp4->Sumw2();
+    hTemp1->Add(hTemp3,-2.);
+    hTemp1->Scale(1./hTemp5->GetEntries());
+    hTemp2->Add(hTemp4,-2.);
+    hTemp2->Scale(1./hTemp6->GetEntries());
+    gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
+    gHistBalanceFunctionHistogram->Scale(0.5/fP2Step[iAnalysisType]);
+  }
+
+  PrintResults(iAnalysisType,gHistBalanceFunctionHistogram);
+
+  return gHistBalanceFunctionHistogram;
+}
diff --git a/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h b/PWGCF/EBYE/BalanceFunctions/AliBalancePsi.h
new file mode 100644 (file)
index 0000000..a5aa8a4
--- /dev/null
@@ -0,0 +1,166 @@
+#ifndef ALIBALANCEPSI_H
+#define ALIBALANCEPSI_H
+/*  See cxx source for full Copyright notice */
+
+
+/* $Id: AliBalancePsi.h 54125 2012-01-24 21:07:41Z miweber $ */
+
+//-------------------------------------------------------------------------
+//                          Class AliBalancePsi
+//   This is the class for the Balance Function writ Psi analysis
+//
+//    Origin: Panos Christakoglou, Nikhef, Panos.Christakoglou@cern.ch
+//-------------------------------------------------------------------------
+
+#include <TObject.h>
+#include "TString.h"
+
+#define ANALYSIS_TYPES 7
+#define MAXIMUM_NUMBER_OF_STEPS        1024
+#define MAXIMUM_STEPS_IN_PSI 360
+
+class TGraphErrors;
+class TH1D;
+class TH2D;
+class TH3D;
+
+const TString gBFPsiAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"};
+
+class AliBalancePsi : public TObject {
+ public:
+  enum EAnalysisType {
+    kRapidity = 0, 
+    kEta      = 1, 
+    kQlong    = 2, 
+    kQout     = 3, 
+    kQside    = 4, 
+    kQinv     = 5, 
+    kPhi      = 6 
+  };
+
+  AliBalancePsi();
+  AliBalancePsi(const AliBalancePsi& balance);
+  ~AliBalancePsi();
+
+  void SetPsiInterval(Double_t psiInterval) {
+    fPsiInterval = psiInterval; 
+    fPsiNumberOfBins = (Int_t)(360./fPsiInterval);}
+  void SetPsiNumberOfBins(Int_t psiNumberOfBins) {
+    fPsiNumberOfBins = psiNumberOfBins;
+    fPsiInterval = 360./fPsiNumberOfBins;}
+
+  void SetCentralityIdentifier(const char* centralityId) {
+    fCentralityId = centralityId;}
+  
+  void SetAnalysisLevel(const char* analysisLevel) {
+    fAnalysisLevel = analysisLevel;}
+  void SetShuffle(Bool_t shuffle) {bShuffle = shuffle;}
+  void SetInterval(Int_t iAnalysisType, Double_t p1Start, Double_t p1Stop,
+                  Int_t ibins, Double_t p2Start, Double_t p2Stop,
+                  Double_t psiInterval);
+  void SetCentralityInterval(Double_t cStart, Double_t cStop)  { fCentStart = cStart; fCentStop = cStop;};
+  void SetNp(Int_t analysisType, Int_t psiBin, Double_t NpSet)   { fNp[analysisType][psiBin] = NpSet; }
+  void SetNn(Int_t analysisType, Int_t psiBin, Double_t NnSet)   { fNn[analysisType][psiBin] = NnSet; }
+  void SetNpp(Int_t analysisType, Int_t psiBin, Int_t ibin, Double_t NppSet) { if(ibin > -1 && ibin < MAXIMUM_NUMBER_OF_STEPS) fNpp[analysisType][psiBin][ibin] = NppSet; }
+  void SetNpn(Int_t analysisType, Int_t psiBin, Int_t ibin, Double_t NpnSet) { if(ibin > -1 && ibin < MAXIMUM_NUMBER_OF_STEPS) fNpn[analysisType][psiBin][ibin] = NpnSet; }
+  void SetNnp(Int_t analysisType, Int_t psiBin, Int_t ibin, Double_t NnpSet) { if(ibin > -1 && ibin < MAXIMUM_NUMBER_OF_STEPS) fNnp[analysisType][psiBin][ibin] = NnpSet; }
+  void SetNnn(Int_t analysisType, Int_t psiBin, Int_t ibin, Double_t NnnSet) { if(ibin > -1 && ibin < MAXIMUM_NUMBER_OF_STEPS) fNnn[analysisType][psiBin][ibin] = NnnSet; }
+
+  void InitHistograms(void);
+
+  const char* GetAnalysisLevel() {return fAnalysisLevel.Data();}
+  Int_t GetNumberOfAnalyzedEvent() {return fAnalyzedEvents;}
+
+  Int_t GetNumberOfBins(Int_t ibin) {return fNumberOfBins[ibin];}
+  Double_t GetP1Start(Int_t ibin){return fP1Start[ibin];}
+  Double_t GetP1Stop(Int_t ibin){return fP1Stop[ibin];}   
+  Double_t GetP2Start(Int_t ibin){return fP2Start[ibin];}
+  Double_t GetP2Stop(Int_t ibin){return fP2Stop[ibin];}    
+  Double_t GetNp(Int_t analysisType, Int_t psiBin) const { return 1.0*fNp[analysisType][psiBin]; }
+  Double_t GetNn(Int_t analysisType, Int_t psiBin) const { return 1.0*fNn[analysisType][psiBin]; }
+  Double_t GetNnn(Int_t analysisType, Int_t psiBin, Int_t p2) const { 
+    return 1.0*fNnn[analysisType][psiBin][p2]; }
+  Double_t GetNpp(Int_t analysisType, Int_t psiBin, Int_t p2) const { 
+    return 1.0*fNpp[analysisType][psiBin][p2]; }
+  Double_t GetNpn(Int_t analysisType, Int_t psiBin, Int_t p2) const { 
+    return 1.0*fNpn[analysisType][psiBin][p2]; }  
+  Double_t GetNnp(Int_t analysisType, Int_t psiBin, Int_t p2) const { 
+    return 1.0*fNnp[analysisType][psiBin][p2]; }
+
+  void CalculateBalance(Float_t fCentrality, Double_t gReactionPlane, vector<Double_t> **chargeVector);
+  
+  Double_t GetBalance(Int_t iAnalysisType, Int_t psiBin, Int_t p2);
+  Double_t GetError(Int_t iAnalysisType, Int_t psiBin, Int_t p2);
+
+  TH3D *GetHistNp(Int_t iAnalysisType) { return fHistP[iAnalysisType];}
+  TH3D *GetHistNn(Int_t iAnalysisType) { return fHistN[iAnalysisType];}
+  TH3D *GetHistNpn(Int_t iAnalysisType) { return fHistPN[iAnalysisType];}
+  TH3D *GetHistNnp(Int_t iAnalysisType) { return fHistNP[iAnalysisType];}
+  TH3D *GetHistNpp(Int_t iAnalysisType) { return fHistPP[iAnalysisType];}
+  TH3D *GetHistNnn(Int_t iAnalysisType) { return fHistNN[iAnalysisType];}
+
+  void PrintAnalysisSettings();
+  TGraphErrors *DrawBalance(Int_t fAnalysisType, Int_t psiBin);
+
+  void SetHistNp(Int_t iAnalysisType, TH3D *gHist) { 
+    fHistP[iAnalysisType] = gHist;}
+  void SetHistNn(Int_t iAnalysisType, TH3D *gHist) { 
+    fHistN[iAnalysisType] = gHist;}
+  void SetHistNpn(Int_t iAnalysisType, TH3D *gHist) { 
+    fHistPN[iAnalysisType] = gHist;}
+  void SetHistNnp(Int_t iAnalysisType, TH3D *gHist) { 
+    fHistNP[iAnalysisType] = gHist;}
+  void SetHistNpp(Int_t iAnalysisType, TH3D *gHist) { 
+    fHistPP[iAnalysisType] = gHist;}
+  void SetHistNnn(Int_t iAnalysisType, TH3D *gHist) { 
+    fHistNN[iAnalysisType] = gHist;}
+
+  TH1D *GetBalanceFunctionHistogram(Int_t iAnalysisType,
+                                   Double_t centrMin, 
+                                   Double_t centrMax, 
+                                   Double_t psiMin, Double_t psiMax);
+  void PrintResults(Int_t iAnalysisType, TH1D *gHist);
+
+ private:
+  Bool_t bShuffle; //shuffled balance function object
+  TString fAnalysisLevel; //ESD, AOD or MC
+  Int_t fAnalyzedEvents; //number of events that have been analyzed
+
+  TString fCentralityId;//Centrality identifier to be used for the histo naming
+
+  Int_t fNumberOfBins[ANALYSIS_TYPES];//number of bins of the analyzed interval
+  Double_t fP1Start[ANALYSIS_TYPES];
+  Double_t fP1Stop[ANALYSIS_TYPES];
+  Double_t fP2Start[ANALYSIS_TYPES];
+  Double_t fP2Stop[ANALYSIS_TYPES];
+  Double_t fP2Step[ANALYSIS_TYPES]; 
+  Double_t fCentStart;
+  Double_t fCentStop;
+
+       
+  Double_t fNnn[ANALYSIS_TYPES][MAXIMUM_STEPS_IN_PSI][MAXIMUM_NUMBER_OF_STEPS]; //N(--)
+  Double_t fNpp[ANALYSIS_TYPES][MAXIMUM_STEPS_IN_PSI][MAXIMUM_NUMBER_OF_STEPS]; //N(++)
+  Double_t fNpn[ANALYSIS_TYPES][MAXIMUM_STEPS_IN_PSI][MAXIMUM_NUMBER_OF_STEPS]; //N(+-)
+  Double_t fNnp[ANALYSIS_TYPES][MAXIMUM_STEPS_IN_PSI][MAXIMUM_NUMBER_OF_STEPS]; //N(-+)
+  Double_t fNn[ANALYSIS_TYPES][MAXIMUM_STEPS_IN_PSI], fNp[ANALYSIS_TYPES][MAXIMUM_STEPS_IN_PSI]; //number of pos./neg. inside the analyzed interval
+  
+  Double_t fB[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS][MAXIMUM_STEPS_IN_PSI]; //BF matrix
+  Double_t ferror[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS][MAXIMUM_STEPS_IN_PSI]; //error of the BF
+  
+  TH3D *fHistP[ANALYSIS_TYPES]; //N+
+  TH3D *fHistN[ANALYSIS_TYPES]; //N-
+  TH3D *fHistPN[ANALYSIS_TYPES]; //N+-
+  TH3D *fHistNP[ANALYSIS_TYPES]; //N-+
+  TH3D *fHistPP[ANALYSIS_TYPES]; //N++
+  TH3D *fHistNN[ANALYSIS_TYPES]; //N--
+
+  Double_t fPsiInterval;// interval in Psi-phi1
+  Int_t fPsiNumberOfBins;// number of bins
+
+  AliBalancePsi & operator=(const AliBalancePsi & ) {return *this;}
+
+  ClassDef(AliBalancePsi, 1)
+};
+
+#endif
diff --git a/PWGCF/EBYE/macros/AddTaskBalancePsiCentralityTrain.C b/PWGCF/EBYE/macros/AddTaskBalancePsiCentralityTrain.C
new file mode 100644 (file)
index 0000000..8efa3cf
--- /dev/null
@@ -0,0 +1,181 @@
+// now in options\r
+//=============================================//\r
+//const char* centralityEstimator = "V0M";\r
+//const char* centralityEstimator = "CL1";\r
+//const char* centralityEstimator = "TRK";\r
+//=============================================//\r
+//Bool_t gRunShuffling = kFALSE;\r
+//Bool_t gRunShuffling = kTRUE;\r
+//=============================================//\r
+\r
+//PID config\r
+Bool_t kUseNSigmaPID = kFALSE;\r
+Double_t nSigmaMax = 3.0;\r
+Bool_t kUseBayesianPID = kTRUE;\r
+Double_t gMinAcceptedProbability = 0.7;\r
+\r
+//_________________________________________________________//\r
+AliAnalysisTaskBF *AddTaskBalancePsiCentralityTrain(Double_t centrMin=0.,\r
+                                                   Double_t centrMax=100.,\r
+                                                   Bool_t gRunShuffling=kFALSE,\r
+                                                   TString centralityEstimator="V0M",\r
+                                                   Double_t vertexZ=10.,\r
+                                                   Double_t DCAxy=-1,\r
+                                                   Double_t DCAz=-1,\r
+                                                   Double_t ptMin=0.3,\r
+                                                   Double_t ptMax=1.5,\r
+                                                   Double_t etaMin=-0.8,\r
+                                                   Double_t etaMax=0.8,\r
+                                                   Double_t maxTPCchi2 = -1, \r
+                                                   Int_t minNClustersTPC = -1,\r
+                                                   Bool_t kUsePID = kFALSE,\r
+                                                   Int_t AODfilterBit = 128,\r
+                                                   Bool_t bCentralTrigger = kFALSE,\r
+                                                   TString fileNameBase="AnalysisResults") {\r
+  // Creates a balance function analysis task and adds it to the analysis manager.\r
+  // Get the pointer to the existing analysis manager via the static access method.\r
+  TString centralityName("");\r
+  centralityName+=Form("%.0f",centrMin);\r
+  centralityName+="-";\r
+  centralityName+=Form("%.0f",centrMax);\r
+  centralityName+="_";\r
+  centralityName+=Form("%s",centralityEstimator.Data());\r
+  centralityName+="_";\r
+  centralityName+=Form("vZ%.1f",vertexZ);\r
+  centralityName+="_";\r
+  centralityName+=Form("DCAxy%.1f",DCAxy);\r
+  centralityName+="_";\r
+  centralityName+=Form("DCAz%.1f",DCAz);\r
+  centralityName+="_Pt";\r
+  centralityName+=Form("%.1f",ptMin);\r
+  centralityName+="-";\r
+  centralityName+=Form("%.1f",ptMax);\r
+  centralityName+="_Eta";\r
+  centralityName+=Form("%.1f",etaMin);\r
+  centralityName+="-";\r
+  centralityName+=Form("%.1f",etaMax);\r
+  centralityName+="_Chi";\r
+  centralityName+=Form("%.1f",maxTPCchi2);\r
+  centralityName+="_nClus";\r
+  centralityName+=Form("%d",minNClustersTPC);\r
+  centralityName+="_Bit";\r
+  centralityName+=Form("%d",AODfilterBit);\r
+  if(bCentralTrigger)   centralityName+="_withCentralTrigger";\r
+\r
+\r
+\r
+\r
+\r
+  TString outputFileName(fileNameBase);\r
+  outputFileName.Append(".root");\r
+\r
+  //===========================================================================\r
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
+  if (!mgr) {\r
+    ::Error("AddTaskBF", "No analysis manager to connect to.");\r
+    return NULL;\r
+  }\r
+\r
+  // Check the analysis type using the event handlers connected to the analysis manager.\r
+  //===========================================================================\r
+  if (!mgr->GetInputEventHandler()) {\r
+    ::Error("AddTaskBF", "This task requires an input event handler");\r
+    return NULL;\r
+  }\r
+  TString analysisType = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"\r
+  if(dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())) analysisType = "MC";\r
+\r
+  // for local changed BF configuration\r
+  //gROOT->LoadMacro("./configBalanceFunctionAnalysis.C");\r
+  gROOT->LoadMacro("$ALICE_ROOT/PWGCF/EBYE/macros/configBalanceFunctionAnalysis.C");\r
+  AliBalancePsi *bf  = 0;  // Balance Function object\r
+  AliBalancePsi *bfs = 0;  // shuffled 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
+  }\r
+  else if (analysisType=="AOD"){\r
+    bf  = GetBalanceFunctionObject("AOD",centralityEstimator,centrMin,centrMax);\r
+    if(gRunShuffling) bfs = GetBalanceFunctionObject("AOD",centralityEstimator,centrMin,centrMax,kTRUE);\r
+  }\r
+  else if (analysisType=="MC"){\r
+    bf  = GetBalanceFunctionObject("MC",centralityEstimator,centrMin,centrMax);\r
+    if(gRunShuffling) bfs = GetBalanceFunctionObject("MC",centralityEstimator,centrMin,centrMax,kTRUE);\r
+  }\r
+  else{\r
+    ::Error("AddTaskBF", "analysis type NOT known.");\r
+    return NULL;\r
+  }\r
+\r
+  // Create the task, add it to manager and configure it.\r
+  //===========================================================================\r
+  AliAnalysisTaskBFPsi *taskBF = new AliAnalysisTaskBFPsi("TaskBFPsi");\r
+  taskBF->SetAnalysisObject(bf);\r
+  if(gRunShuffling) taskBF->SetShufflingObject(bfs);\r
+\r
+  taskBF->SetCentralityPercentileRange(centrMin,centrMax);\r
+  if(analysisType == "ESD") {\r
+    AliESDtrackCuts *trackCuts = GetTrackCutsObject(ptMin,ptMax,etaMin,etaMax,maxTPCchi2,DCAxy,DCAz,minNClustersTPC);\r
+    taskBF->SetAnalysisCutObject(trackCuts);\r
+    if(kUsePID) {\r
+      if(kUseBayesianPID)\r
+       taskBF->SetUseBayesianPID(gMinAcceptedProbability);\r
+      else if(kUseNSigmaPID)\r
+       taskBF->SetUseNSigmaPID(nSigmaMax);\r
+      taskBF->SetParticleOfInterest(AliAnalysistaskBFPsi::kProton);\r
+      taskBF->SetDetectorUsedForPID(AliAnalysisTaskBFPsi::kTOFpid);\r
+    }\r
+  }\r
+  else if(analysisType == "AOD") {\r
+    // pt and eta cut (pt_min, pt_max, eta_min, eta_max)\r
+    taskBF->SetAODtrackCutBit(AODfilterBit);\r
+    taskBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax);\r
+\r
+    // set extra DCA cuts (-1 no extra cut)\r
+    taskBF->SetExtraDCACutsAOD(DCAxy,DCAz);\r
+\r
+    // set extra TPC chi2 / nr of clusters cut\r
+    taskBF->SetExtraTPCCutsAOD(maxTPCchi2, minNClustersTPC);\r
+    \r
+  }\r
+  else if(analysisType == "MC") {\r
+    taskBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax); \r
+  }\r
+\r
+  // offline trigger selection (AliVEvent.h)\r
+  // taskBF->UseOfflineTrigger(); // NOT used (selection is done with the AliAnalysisTaskSE::SelectCollisionCandidates()) \r
+  // with this only selected events are analyzed (first 2 bins in event QA histogram are the same))\r
+  // documentation in https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PWG1EvSelDocumentation\r
+  if(bCentralTrigger) taskBF->SelectCollisionCandidates(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral);\r
+  else                taskBF->SelectCollisionCandidates(AliVEvent::kMB);\r
+\r
+  // centrality estimator (default = V0M)\r
+  taskBF->SetCentralityEstimator(centralityEstimator);\r
+  \r
+  // vertex cut (x,y,z)\r
+  taskBF->SetVertexDiamond(.3,.3,vertexZ);\r
+  \r
+\r
+\r
+  //bf->PrintAnalysisSettings();\r
+  mgr->AddTask(taskBF);\r
+  \r
+  // Create ONLY the output containers for the data produced by the task.\r
+  // Get and connect other common input/output containers via the manager as below\r
+  //==============================================================================\r
+  TString outputFileName = AliAnalysisManager::GetCommonFileName();\r
+  outputFileName += ":PWGCFEbyE.outputBalanceFunctionAnalysis";\r
+  AliAnalysisDataContainer *coutQA = mgr->CreateContainer(Form("listQA_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  AliAnalysisDataContainer *coutBF = mgr->CreateContainer(Form("listBF_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  if(gRunShuffling) AliAnalysisDataContainer *coutBFS = mgr->CreateContainer(Form("listBFShuffled_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  if(kUsePID) AliAnalysisDataContainer *coutQAPID = mgr->CreateContainer(Form("listQAPID_%s",centralityName.Data()), 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
+\r
+  return taskBF;\r
+}\r
diff --git a/PWGCF/EBYE/macros/configBalanceFunctionPsiAnalysis.C b/PWGCF/EBYE/macros/configBalanceFunctionPsiAnalysis.C
new file mode 100644 (file)
index 0000000..ed2b8d3
--- /dev/null
@@ -0,0 +1,63 @@
+//__________________________________________________//\r
+AliBalancePsi *GetBalanceFunctionObject(const char* analysisLevel = "ESD", \r
+                                       const char* centralityName = 0x0,\r
+                                       Double_t centrMin = 0.,\r
+                                       Double_t centrMax = 100.,\r
+                                       Bool_t bShuffle = kFALSE) {\r
+  //Function to setup the AliBalance object and return it\r
+  AliBalancePsi *gBalance = new AliBalancePsi();\r
+  gBalance->SetAnalysisLevel(analysisLevel);\r
+  gBalance->SetShuffle(bShuffle);\r
+  if(centralityName) gBalance->SetCentralityIdentifier(centralityName);\r
+  gBalance->SetCentralityInterval(centrMin,centrMax);\r
+\r
+  //Set all analyses separately\r
+  //Rapidity\r
+  gBalance->SetInterval(AliBalance::kRapidity,-0.8,0.8,16,-0.8,0.8,15.);  \r
+  //Eta\r
+  gBalance->SetInterval(AliBalance::kEta,-0.8,0.8,16,-0.8,0.8,15);\r
+  //Qlong\r
+  gBalance->SetInterval(AliBalance::kQlong,-1,1,200,0.0,4.0,15);\r
+  //Qout\r
+  gBalance->SetInterval(AliBalance::kQout,-1,1,200,0.0,4.0,15);\r
+  //Qside\r
+  gBalance->SetInterval(AliBalance::kQside,-1,1,200,0.0,4.0,15);\r
+  //Qinv\r
+  gBalance->SetInterval(AliBalance::kQinv,-1,1,200,0.0,4.0,15);\r
+  //Phi\r
+  gBalance->SetInterval(AliBalance::kPhi,0.,360.,90,-180.,180.0,15);\r
+\r
+  //Init the histograms\r
+  gBalance->InitHistograms();\r
+  \r
+  return gBalance;\r
+}\r
+\r
+//__________________________________________________//\r
+AliESDtrackCuts *GetTrackCutsObject(Double_t ptMin, Double_t ptMax, Double_t etaMin, Double_t etaMax, Double_t maxTPCchi2, Double_t maxDCAz, Double_t maxDCAxy, Int_t minNClustersTPC) {\r
+\r
+  // only used for ESDs\r
+  // Function to setup the AliESDtrackCuts object and return it\r
+  AliESDtrackCuts *cuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();\r
+  cuts->SetRequireTPCStandAlone(kTRUE); // TPC only cuts!  \r
+\r
+  // extra TPC cuts (Syst studies)\r
+  if(minNClustersTPC != -1)  cuts->SetMinNClustersTPC(minNClustersTPC);\r
+  else cuts->SetMinNClustersTPC(70); // standard for filter bit 128\r
+  \r
+  if(maxTPCchi2 != -1) cuts->SetMaxChi2PerClusterTPC(maxTPCchi2);\r
+\r
+  // extra DCA cuts (Syst studies)  \r
+  if(maxDCAz!=-1 && maxDCAxy != -1){\r
+    cuts->SetMaxDCAToVertexZ(maxDCAz);\r
+    cuts->SetMaxDCAToVertexXY(maxDCAxy);\r
+  }\r
+\r
+  cuts->SetPtRange(ptMin,ptMax);\r
+  cuts->SetEtaRange(etaMin,etaMax);\r
+  cuts->DefineHistograms(1);\r
+  //cuts->SaveHistograms("trackCuts");\r
+\r
+  return cuts;\r
+}\r
+\r
diff --git a/PWGCF/EBYE/macros/runBalanceFunctionPsi.C b/PWGCF/EBYE/macros/runBalanceFunctionPsi.C
new file mode 100755 (executable)
index 0000000..4d927ff
--- /dev/null
@@ -0,0 +1,456 @@
+// run.C\r
+//\r
+// Template run macro for AliBasicTask.cxx/.h with example layout of\r
+// physics selections and options, in macro and task.\r
+//\r
+// Author: Arvinder Palaha\r
+//\r
+class AliAnalysisGrid;\r
+class AliAnalysisTaskBF;\r
+class AliBalance;\r
+\r
+//Centrality stuff\r
+Int_t binfirst = 0;  //where do we start numbering bins\r
+Int_t binlast = 8;  //where do we stop numbering bins\r
+const Int_t numberOfCentralityBins = 9;\r
+Float_t centralityArray[numberOfCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.}; // in centrality percentile\r
+\r
+//Systematic studies\r
+const Int_t numberOfSyst = 13;\r
+Float_t vZ[numberOfSyst]     = {10.,12.,6.,8.,10.,10.,10.,10.,10.,10.,10.,10.,10.};     // global Vertex Z cut\r
+Float_t DCAxy[numberOfSyst]  = {-1.,2.4,2.4,2.4,2.2,2.0,1.8,2.4,2.4,2.4,2.4,2.4,2.4};   // DCA xy cut (afterburner, -1 = w/o additional cut)\r
+Float_t DCAz[numberOfSyst]   = {-1.,3.2,3.2,3.2,3.0,2.8,2.6,3.2,3.2,3.2,3.2,3.2,3.2};   // DCA z cut (afterburner, -1 = w/o additional cut)\r
+Float_t ptMin[numberOfSyst]  = {0.3,0.3,0.3,0.3,0.3,0.3,0.3,1.5,5.0,0.3,0.3,0.3,0.3};   // pt cuts\r
+Float_t ptMax[numberOfSyst]  = {1.5,1.5,1.5,1.5,1.5,1.5,1.5,5.0,10.0,10.0,1.5,1.5,1.5}; // pt cuts\r
+Float_t etaMin[numberOfSyst] = {-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,-1.0,-0.6,-0.4}; // eta cuts\r
+Float_t etaMax[numberOfSyst] = {0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8,0.8,1.0,0.6,0.4};   // eta cuts\r
+\r
+Bool_t kUsePID = kFALSE;\r
+\r
+//______________________________________________________________________________\r
+void runBalanceFunctionPsi(\r
+         const char* runtype = "local", // local, proof or grid\r
+         const char *gridmode = "test", // Set the run mode (can be "full", "test", "offline", "submit" or "terminate"). Full & Test work for proof\r
+        const Int_t bunchN = 0,\r
+         const bool bAOD = 1, // 1 = AOD ANALYSIS, 0 = ESD ANALYSIS\r
+         const bool bMCtruth = 0, // 1 = MCEvent handler is on (MC truth), 0 = MCEvent handler is off (MC reconstructed/real data)\r
+         const bool bMCphyssel = 0, // 1 = looking at MC truth or reconstructed, 0 = looking at real data\r
+         const Long64_t nentries = 50000, // for local and proof mode, ignored in grid mode. Set to 1234567890 for all events.\r
+         const Long64_t firstentry = 0, // for local and proof mode, ignored in grid mode\r
+         TString proofdataset = "bunchPROOF", // path to dataset on proof cluster, for proof analysis\r
+         const char *proofcluster = "miweber@alice-caf.cern.ch", // which proof cluster to use in proof mode\r
+         const char *taskname = "BF_Syst_Test" // sets name of grid generated macros\r
+         )\r
+{\r
+    // check run type\r
+    if(runtype != "local" && runtype != "proof" && runtype != "grid"){\r
+        Printf("\n\tIncorrect run option, check first argument of run macro");\r
+        Printf("\tint runtype = local, proof or grid\n");\r
+        return;\r
+    }\r
+    Printf("%s analysis chosen",runtype);\r
+  \r
+    // load libraries\r
+    gSystem->Load("libCore.so");        \r
+    gSystem->Load("libGeom.so");\r
+    gSystem->Load("libVMC.so");\r
+    gSystem->Load("libPhysics.so");\r
+    gSystem->Load("libTree.so");\r
+    gSystem->Load("libSTEERBase.so");\r
+    gSystem->Load("libESD.so");\r
+    gSystem->Load("libAOD.so");\r
+    gSystem->Load("libANALYSIS.so");\r
+    gSystem->Load("libANALYSISalice.so");\r
+    gSystem->Load("libPWGCFebye.so");\r
+\r
+    // additional\r
+\r
+    // compile standalone stuff\r
+    //gROOT->LoadMacro("AliBalance.cxx++g");\r
+    //gROOT->LoadMacro("AliAnalysisTaskBF.cxx++g");\r
+\r
+    // add aliroot indlude path\r
+    //gROOT->ProcessLine(".include $PWD/.");\r
+    //gROOT->ProcessLine(Form(".include %s/include",gSystem->ExpandPathName("$ALICE_ROOT")));\r
+\r
+    gROOT->SetStyle("Plain");\r
+\r
+    // analysis manager\r
+    AliAnalysisManager* mgr = new AliAnalysisManager(Form("%s%i",taskname,bunchN));\r
+    \r
+    // create the alien handler and attach it to the manager\r
+    AliAnalysisGrid *plugin = CreateAlienHandler(bAOD,bunchN,Form("%s%i",taskname,bunchN), gridmode, proofcluster, Form("%s_%d.txt",proofdataset.Data(),bunchN)); \r
+    mgr->SetGridHandler(plugin);\r
+    \r
+\r
+    // input handler (ESD or AOD)\r
+    AliVEventHandler* inputH = NULL;\r
+    if(!bAOD){\r
+      inputH = new AliESDInputHandler();\r
+    }\r
+    else{\r
+      inputH = new AliAODInputHandler();\r
+    }\r
+    mgr->SetInputEventHandler(inputH);\r
+    \r
+    // mc event handler\r
+    if(bMCtruth) {\r
+        AliMCEventHandler* mchandler = new AliMCEventHandler();\r
+        // Not reading track references\r
+        mchandler->SetReadTR(kFALSE);\r
+        mgr->SetMCtruthEventHandler(mchandler);\r
+    }   \r
+\r
+    // AOD output handler\r
+    //AliAODHandler* aodoutHandler = new AliAODHandler();\r
+    //aodoutHandler->SetOutputFileName("aod.root");\r
+    //mgr->SetOutputEventHandler(aodoutHandler); \r
+    \r
+    // === Physics Selection Task ===\r
+    //\r
+    // In SelectCollisionCandidate(), default is kMB, so the task UserExec() \r
+    // function is only called for these events.\r
+    // Options are:\r
+    //    kMB             Minimum Bias trigger\r
+    //    kMBNoTRD        Minimum bias trigger where the TRD is not read out\r
+    //    kMUON           Muon trigger\r
+    //    kHighMult       High-Multiplicity Trigger\r
+    //    kUserDefined    For manually defined trigger selection\r
+    //\r
+    // Multiple options possible with the standard AND/OR operators && and ||\r
+    // These all have the usual offline SPD or V0 selections performed.\r
+    //\r
+    // With a pointer to the physics selection object using physSelTask->GetPhysicsSelection(),\r
+    // one can manually set the selected and background classes using:\r
+    //    AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL")\r
+    //    AddBGTriggerClass("+CINT1A-ABCE-NOPF-ALL");\r
+    //\r
+    // One can also specify multiple classes at once, or require a class to NOT\r
+    // trigger, for e.g.\r
+    //    AddBGTriggerClass("+CSMBA-ABCE-NOPF-ALL -CSMBB-ABCE-NOPF-ALL");\r
+    //\r
+    // NOTE that manually setting the physics selection overrides the standard\r
+    // selection, so it must be done in completeness.\r
+    //\r
+    // ALTERNATIVELY, one can make the physics selection inside the task\r
+    // UserExec().\r
+    // For this case, comment out the task->SelectCol.... line, \r
+    // and see AliBasicTask.cxx UserExec() function for details on this.\r
+\r
+    //gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");\r
+    //AliPhysicsSelectionTask *physSelTask = AddTaskPhysicsSelection(bMCphyssel);\r
+    //if(!physSelTask) { Printf("no physSelTask"); return; }\r
+    //AliPhysicsSelection *physSel = physSelTask->GetPhysicsSelection();\r
+    //physSel->AddCollisionTriggerClass("+CINT1B-ABCE-NOPF-ALL");// #3119 #769");\r
+                \r
+    // create task\r
+    //Add the centrality determination task and the physics selection \r
+    // (only on ESD level, in AODs centrality is already in header and events are selected)\r
+    if(!bAOD){\r
+      gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C");\r
+      AliCentralitySelectionTask *taskCentrality = AddTaskCentrality();\r
+\r
+      // Add physics selection task (NOT needed for AODs)\r
+      gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");\r
+      AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection(bMCphyssel);\r
+\r
+      //Add the PID response\r
+      if(kUsePID) {\r
+       gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPIDResponse.C");\r
+       AddTaskPIDResponse(bMCphyssel); \r
+      }\r
+    }\r
+\r
+    //Add the VZERO event plane task\r
+    gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskVZEROEPSelection.C"); \r
+    AliVZEROEPSelectionTask* epSelTask = AddTaskVZEROEPSelection();\r
+\r
+    //Add the BF task (all centralities)\r
+    gROOT->LoadMacro("AddTaskBalancePsiCentralityTrain.C"); \r
+    AliAnalysisTaskBFPsi *task = AddTaskBalancePsiCentralityTrain(0,100,0,"V0M",vZ[0],DCAxy[0],DCAz[0],ptMin[0],ptMax[0],etaMin[0],etaMax[0],-1,-1,kUsePID);\r
+    \r
+    // enable debug printouts\r
+    //mgr->SetDebugLevel(2);\r
+    //mgr->SetUseProgressBar(1,100);\r
+    if (!mgr->InitAnalysis()) return;\r
+    mgr->PrintStatus();\r
+  \r
+    // start analysis\r
+    Printf("Starting Analysis....");\r
+    mgr->StartAnalysis(runtype,nentries,firstentry);\r
+}\r
+\r
+//______________________________________________________________________________\r
+AliAnalysisGrid* CreateAlienHandler(Bool_t bAOD, Int_t bunchN, const char *taskname, const char *gridmode, const char *proofcluster, const char *proofdataset)\r
+{\r
+    AliAnalysisAlien *plugin = new AliAnalysisAlien();\r
+    // Set the run mode (can be "full", "test", "offline", "submit" or "terminate")\r
+    plugin->SetRunMode(gridmode);\r
+\r
+    // Set versions of used packages\r
+    plugin->SetAPIVersion("V1.1x");\r
+    plugin->SetROOTVersion("v5-28-00d");\r
+    plugin->SetAliROOTVersion("v5-02-05-AN");\r
+\r
+    // Declare input data to be processed.\r
+\r
+    // Method 1: Create automatically XML collections using alien 'find' command.\r
+    // Define production directory LFN\r
+    plugin->SetGridDataDir("/alice/data/2010/LHC10h/");\r
+    // On real reconstructed data:\r
+    // plugin->SetGridDataDir("/alice/data/2009/LHC09d");\r
+\r
+    // Set data search pattern\r
+    //plugin->SetDataPattern("*ESDs.root"); // THIS CHOOSES ALL PASSES\r
+    // Data pattern for reconstructed data\r
+    if(!bAOD){\r
+      plugin->SetDataPattern("*ESDs/pass2/*ESDs.root"); // CHECK LATEST PASS OF DATA SET IN ALIENSH\r
+    } \r
+    else{\r
+      plugin->SetDataPattern("*ESDs/pass2/AOD049/*/AliAOD.root");\r
+    }\r
+\r
+    plugin->SetRunPrefix("000");   // real data\r
+    // ...then add run numbers to be considered\r
+    //plugin->SetRunRange(114917,115322);\r
+\r
+    if(bunchN==0){\r
+      plugin->AddRunNumber(137366);\r
+    }\r
+    \r
+    //bunch1\r
+    else if(bunchN == 1){\r
+      plugin->AddRunNumber(139510);\r
+      plugin->AddRunNumber(139507);\r
+      plugin->AddRunNumber(139505);\r
+      plugin->AddRunNumber(139503); \r
+      plugin->AddRunNumber(139465); \r
+      plugin->AddRunNumber(139438);\r
+      plugin->AddRunNumber(139437);\r
+      plugin->AddRunNumber(139360); \r
+      plugin->AddRunNumber(139329);\r
+      plugin->AddRunNumber(139328); \r
+    }\r
+\r
+    //bunch2\r
+    else if(bunchN == 2){\r
+      plugin->AddRunNumber(139314); \r
+      plugin->AddRunNumber(139310);\r
+      plugin->AddRunNumber(139309); \r
+      plugin->AddRunNumber(139173); \r
+      plugin->AddRunNumber(139107); \r
+      plugin->AddRunNumber(139105); \r
+      plugin->AddRunNumber(139038); \r
+      plugin->AddRunNumber(139037); \r
+      plugin->AddRunNumber(139036); \r
+      plugin->AddRunNumber(139029); \r
+      plugin->AddRunNumber(139028); \r
+      plugin->AddRunNumber(138872); \r
+      plugin->AddRunNumber(138871); \r
+      plugin->AddRunNumber(138870); \r
+      plugin->AddRunNumber(138837); \r
+      plugin->AddRunNumber(138732); \r
+      plugin->AddRunNumber(138730);\r
+      plugin->AddRunNumber(138666);\r
+      plugin->AddRunNumber(138662); \r
+      plugin->AddRunNumber(138653); \r
+    }\r
+\r
+    else if(bunchN == 3){\r
+      plugin->AddRunNumber(138652);\r
+      plugin->AddRunNumber(138638);\r
+      plugin->AddRunNumber(138624); \r
+      plugin->AddRunNumber(138621); \r
+      plugin->AddRunNumber(138583); \r
+      plugin->AddRunNumber(138582); \r
+      plugin->AddRunNumber(138579); \r
+      plugin->AddRunNumber(138578);\r
+      plugin->AddRunNumber(138534);\r
+      plugin->AddRunNumber(138469); \r
+    }\r
+\r
+    else if(bunchN == 4){\r
+      \r
+      plugin->AddRunNumber(138442);\r
+      plugin->AddRunNumber(138439);\r
+      plugin->AddRunNumber(138438);\r
+      plugin->AddRunNumber(138396); \r
+      plugin->AddRunNumber(138364); \r
+      plugin->AddRunNumber(138275); \r
+      plugin->AddRunNumber(138225); \r
+      plugin->AddRunNumber(138201);\r
+      plugin->AddRunNumber(138197); \r
+      plugin->AddRunNumber(138192); \r
+    }\r
+\r
+    else if(bunchN == 5){\r
+\r
+      plugin->AddRunNumber(138190);\r
+      plugin->AddRunNumber(137848); \r
+      plugin->AddRunNumber(137844); \r
+      plugin->AddRunNumber(137752); \r
+      plugin->AddRunNumber(137751); \r
+      plugin->AddRunNumber(137724); \r
+      plugin->AddRunNumber(137722); \r
+      plugin->AddRunNumber(137718); \r
+      plugin->AddRunNumber(137704); \r
+      plugin->AddRunNumber(137693);\r
+    }\r
+\r
+    else if(bunchN == 6){\r
+\r
+      plugin->AddRunNumber(137692); \r
+      plugin->AddRunNumber(137691); \r
+      plugin->AddRunNumber(137686); \r
+      plugin->AddRunNumber(137685); \r
+      plugin->AddRunNumber(137639); \r
+      plugin->AddRunNumber(137638);\r
+      plugin->AddRunNumber(137608); \r
+      plugin->AddRunNumber(137595);\r
+      plugin->AddRunNumber(137549);\r
+      plugin->AddRunNumber(137546); \r
+\r
+    }\r
+\r
+    else if(bunchN == 7){\r
+\r
+      plugin->AddRunNumber(137544); \r
+      plugin->AddRunNumber(137541); \r
+      plugin->AddRunNumber(137539); \r
+      plugin->AddRunNumber(137531); \r
+      plugin->AddRunNumber(137530); \r
+      plugin->AddRunNumber(137443); \r
+      plugin->AddRunNumber(137441); \r
+      plugin->AddRunNumber(137440); \r
+      plugin->AddRunNumber(137439); \r
+      plugin->AddRunNumber(137434); \r
+\r
+    }\r
+\r
+    else if(bunchN == 8){\r
+\r
+      plugin->AddRunNumber(137432); \r
+      plugin->AddRunNumber(137431); \r
+      plugin->AddRunNumber(137430); \r
+      plugin->AddRunNumber(137366); \r
+      plugin->AddRunNumber(137243); \r
+      plugin->AddRunNumber(137236);\r
+      plugin->AddRunNumber(137235);\r
+      plugin->AddRunNumber(137232); \r
+      plugin->AddRunNumber(137231); \r
+      plugin->AddRunNumber(137162); \r
+      plugin->AddRunNumber(137161);\r
+    }\r
+\r
+    else{\r
+\r
+      stderr<<"BUNCH NOT THERE"<<endl;\r
+      return NULL;\r
+\r
+    }\r
+\r
+\r
+    //plugin->AddRunList("139510, 139507, 139505, 139503, 139465, 139438, 139437, 139360, 139329, 139328, 139314, 139310, 139309, 139173, 139107, 139105, 139038, 139037, 139036, 139029, 139028, 138872, 138871, 138870, 138837, 138732, 138730, 138666, 138662, 138653, 138652, 138638, 138624, 138621, 138583, 138582, 138579, 138578, 138534, 138469, 138442, 138439, 138438, 138396, 138364, 138275, 138225, 138201, 138197, 138192, 138190, 137848, 137844, 137752, 137751, 137724, 137722, 137718, 137704, 137693, 137692, 137691, 137686, 137685, 137639, 137638, 137608, 137595, 137549, 137546, 137544, 137541, 137539, 137531, 137530, 137443, 137441, 137440, 137439, 137434, 137432, 137431, 137430, 137366, 137243, 137236, 137235, 137232, 137231, 137162, 137161");\r
+\r
+\r
+\r
+\r
+\r
+    plugin->SetNrunsPerMaster(1);\r
+    plugin->SetOutputToRunNo();\r
+    // comment out the next line when using the "terminate" option, unless\r
+    // you want separate merged files for each run\r
+    plugin->SetMergeViaJDL();\r
+\r
+    // Method 2: Declare existing data files (raw collections, xml collections, root file)\r
+    // If no path mentioned data is supposed to be in the work directory (see SetGridWorkingDir())\r
+    // XML collections added via this method can be combined with the first method if\r
+    // the content is compatible (using or not tags)\r
+    //   plugin->AddDataFile("tag.xml");\r
+    //   plugin->AddDataFile("/alice/data/2008/LHC08c/000057657/raw/Run57657.Merged.RAW.tag.root");\r
+\r
+    // Define alien work directory where all files will be copied. Relative to alien $HOME.\r
+    plugin->SetGridWorkingDir(taskname);\r
+\r
+    // Declare alien output directory. Relative to working directory.\r
+    plugin->SetGridOutputDir("out"); // In this case will be $HOME/taskname/out\r
+\r
+   // Declare the analysis source files names separated by blancs. To be compiled runtime\r
+    // using ACLiC on the worker nodes.\r
+    plugin->SetAnalysisSource("AliBalance.cxx AliAnalysisTaskBF.cxx");\r
+\r
+    // Declare all libraries (other than the default ones for the framework. These will be\r
+    // loaded by the generated analysis macro. Add all extra files (task .cxx/.h) here.\r
+    //plugin->AddIncludePath("-I.");\r
+    //plugin->SetAdditionalLibs("libPWGCFebye.so");\r
+    plugin->SetAdditionalLibs("AliBalance.cxx AliBalance.h AliAnalysisTaskBF.cxx AliAnalysisTaskBF.h");\r
+\r
+     // Declare the output file names separated by blancs.\r
+    // (can be like: file.root or file.root@ALICE::Niham::File)\r
+    // To only save certain files, use SetDefaultOutputs(kFALSE), and then\r
+    // SetOutputFiles("list.root other.filename") to choose which files to save\r
+    plugin->SetDefaultOutputs();\r
+    //plugin->SetOutputFiles("list.root");\r
+\r
+    // Optionally set a name for the generated analysis macro (default MyAnalysis.C)\r
+    plugin->SetAnalysisMacro(Form("%s.C",taskname));\r
+\r
+    // Optionally set maximum number of input files/subjob (default 100, put 0 to ignore)\r
+    plugin->SetSplitMaxInputFileNumber(100);\r
+\r
+    // Optionally modify the executable name (default analysis.sh)\r
+    plugin->SetExecutable(Form("%s.sh",taskname));\r
+\r
+    // set number of test files to use in "test" mode\r
+    plugin->SetNtestFiles(1);\r
+\r
+    // Optionally resubmit threshold.\r
+    plugin->SetMasterResubmitThreshold(90);\r
+\r
+    // Optionally set time to live (default 30000 sec)\r
+    plugin->SetTTL(90000);\r
+\r
+    // Optionally set input format (default xml-single)\r
+    plugin->SetInputFormat("xml-single");\r
+\r
+    // Optionally modify the name of the generated JDL (default analysis.jdl)\r
+    plugin->SetJDLName(Form("%s.jdl",taskname));\r
+\r
+    // Optionally modify job price (default 1)\r
+    plugin->SetPrice(1);      \r
+\r
+    // Optionally modify split mode (default 'se')    \r
+    plugin->SetSplitMode("se");\r
+\r
+    //plugin->SetUseSubmitPolicy();\r
+    //plugin->SetKeepLogs();\r
+    \r
+    //----------------------------------------------------------\r
+    //---      PROOF MODE SPECIFIC SETTINGS         ------------\r
+    //---------------------------------------------------------- \r
+    // Proof cluster\r
+    plugin->SetProofCluster(proofcluster);\r
+    // Dataset to be used   \r
+    plugin->SetProofDataSet(proofdataset);\r
+    // May need to reset proof. Supported modes: 0-no reset, 1-soft, 2-hard\r
+    plugin->SetProofReset(0);\r
+    // May limit number of workers\r
+    plugin->SetNproofWorkers(0);\r
+    // May limit the number of workers per slave\r
+    plugin->SetNproofWorkersPerSlave(1);   \r
+    // May use a specific version of root installed in proof\r
+    plugin->SetRootVersionForProof("current");\r
+    // May set the aliroot mode. Check http://aaf.cern.ch/node/83 \r
+    plugin->SetAliRootMode("default"); // Loads AF libs by default\r
+    // May request ClearPackages (individual ClearPackage not supported)\r
+    plugin->SetClearPackages(kFALSE);\r
+    // Plugin test mode works only providing a file containing test file locations, used in "local" mode also\r
+    plugin->SetFileForTestMode("files.txt"); // file should contain path name to a local directory containg *ESDs.root etc\r
+    // Request connection to alien upon connection to grid\r
+    plugin->SetProofConnectGrid(kFALSE);\r
+\r
+    plugin->Print();\r
+\r
+    return plugin;\r
+}\r
+\r
index ceafcce..1c2ad3c 100644 (file)
@@ -5,7 +5,9 @@
 #pragma link off all functions;
 
 #pragma link C++ class AliBalance+;
+#pragma link C++ class AliBalancePsi+;
 #pragma link C++ class AliAnalysisTaskBF+;
+#pragma link C++ class AliAnalysisTaskBFPsi+;
 
 #pragma link C++ class AliAnalysisTaskToyModel+;