Event Mixing for Balance Function Analysis
authormiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 23 Apr 2012 10:18:22 +0000 (10:18 +0000)
committermiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 23 Apr 2012 10:18:22 +0000 (10:18 +0000)
PWGCF/CMakelibPWGCFebye.pkg
PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEventMixingBF.cxx [new file with mode: 0755]
PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEventMixingBF.h [new file with mode: 0755]
PWGCF/EBYE/BalanceFunctions/AliBalanceEventMixing.cxx [new file with mode: 0644]
PWGCF/EBYE/BalanceFunctions/AliBalanceEventMixing.h [new file with mode: 0644]
PWGCF/EBYE/macros/AddMixingHandler.C [new file with mode: 0644]
PWGCF/EBYE/macros/AddTaskBalanceEventMixing.C [new file with mode: 0644]
PWGCF/EBYE/macros/configBalanceFunctionAnalysisEventMixing.C [new file with mode: 0644]
PWGCF/PWGCFebyeLinkDef.h

index f1604c1..0837d00 100644 (file)
@@ -30,8 +30,10 @@ set ( SRCS
     EBYE/BalanceFunctions/AliAnalysisTaskToyModel.cxx 
     EBYE/BalanceFunctions/AliAnalysisTaskBF.cxx 
     EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx 
+    EBYE/BalanceFunctions/AliAnalysisTaskEventMixingBF.cxx 
     EBYE/BalanceFunctions/AliBalance.cxx 
     EBYE/BalanceFunctions/AliBalancePsi.cxx 
+    EBYE/BalanceFunctions/AliBalanceEventMixing.cxx 
     EBYE/LRC/AliLRCBase.cxx
     EBYE/LRC/AliAnalysisTaskLRC.cxx 
     EBYE/LRC/AliLRCAnalysis.cxx 
diff --git a/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEventMixingBF.cxx b/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEventMixingBF.cxx
new file mode 100755 (executable)
index 0000000..9b9aa62
--- /dev/null
@@ -0,0 +1,683 @@
+#include "TChain.h"\r
+#include "TList.h"\r
+#include "TCanvas.h"\r
+#include "TLorentzVector.h"\r
+#include "TGraphErrors.h"\r
+#include "TH1F.h"\r
+#include "TH2F.h"\r
+#include "TArrayF.h"\r
+#include "TF1.h"\r
+#include "TRandom.h"\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+#include "AliAnalysisManager.h"\r
+\r
+#include "AliESDVertex.h"\r
+#include "AliESDEvent.h"\r
+#include "AliESDInputHandler.h"\r
+#include "AliAODEvent.h"\r
+#include "AliAODTrack.h"\r
+#include "AliAODInputHandler.h"\r
+#include "AliGenEventHeader.h"\r
+#include "AliGenHijingEventHeader.h"\r
+#include "AliMCEventHandler.h"\r
+#include "AliMCEvent.h"\r
+#include "AliMixInputEventHandler.h"\r
+#include "AliStack.h"\r
+#include "AliESDtrackCuts.h"\r
+\r
+#include "TH2D.h"                  \r
+#include "AliPID.h"                \r
+#include "AliPIDResponse.h"        \r
+#include "AliPIDCombined.h"        \r
+\r
+#include "AliAnalysisTaskEventMixingBF.h"\r
+#include "AliBalanceEventMixing.h"\r
+\r
+\r
+// Analysis task for the EventMixingBF code\r
+// Authors: Panos.Christakoglou@nikhef.nl, m.weber@cern.ch\r
+\r
+ClassImp(AliAnalysisTaskEventMixingBF)\r
+\r
+//________________________________________________________________________\r
+AliAnalysisTaskEventMixingBF::AliAnalysisTaskEventMixingBF(const char *name) \r
+: AliAnalysisTaskSE(name), \r
+  fBalance(0),\r
+  fRunShuffling(kFALSE),\r
+  fShuffledBalance(0),\r
+  fList(0),\r
+  fListEventMixingBF(0),\r
+  fListEventMixingBFS(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
+  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
+AliAnalysisTaskEventMixingBF::~AliAnalysisTaskEventMixingBF() {\r
+\r
+  // delete fBalance; \r
+  // delete fShuffledBalance; \r
+  // delete fList;\r
+  // delete fListEventMixingBF; \r
+  // delete fListEventMixingBFS;\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 AliAnalysisTaskEventMixingBF::UserCreateOutputObjects() {\r
+  // Create histograms\r
+  // Called once\r
+  if(!fBalance) {\r
+    fBalance = new AliBalanceEventMixing();\r
+    fBalance->SetAnalysisLevel("ESD");\r
+    //fBalance->SetNumberOfBins(-1,16);\r
+    fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);\r
+  }\r
+  if(fRunShuffling) {\r
+    if(!fShuffledBalance) {\r
+      fShuffledBalance = new AliBalanceEventMixing();\r
+      fShuffledBalance->SetAnalysisLevel("ESD");\r
+      //fShuffledBalance->SetNumberOfBins(-1,16);\r
+      fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);\r
+    }\r
+  }\r
+\r
+  //QA list\r
+  fList = new TList();\r
+  fList->SetName("listQA");\r
+  fList->SetOwner();\r
+\r
+  //Balance Function list\r
+  fListEventMixingBF = new TList();\r
+  fListEventMixingBF->SetName("listEventMixingBF");\r
+  fListEventMixingBF->SetOwner();\r
+\r
+  if(fRunShuffling) {\r
+    fListEventMixingBFS = new TList();\r
+    fListEventMixingBFS->SetName("listEventMixingBFShuffled");\r
+    fListEventMixingBFS->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
+  // 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
+    fListEventMixingBF->Add(fBalance->GetHistNp(a));\r
+    fListEventMixingBF->Add(fBalance->GetHistNn(a));\r
+    fListEventMixingBF->Add(fBalance->GetHistNpn(a));\r
+    fListEventMixingBF->Add(fBalance->GetHistNnn(a));\r
+    fListEventMixingBF->Add(fBalance->GetHistNpp(a));\r
+    fListEventMixingBF->Add(fBalance->GetHistNnp(a));\r
+\r
+    if(fRunShuffling) {\r
+      fListEventMixingBFS->Add(fShuffledBalance->GetHistNp(a));\r
+      fListEventMixingBFS->Add(fShuffledBalance->GetHistNn(a));\r
+      fListEventMixingBFS->Add(fShuffledBalance->GetHistNpn(a));\r
+      fListEventMixingBFS->Add(fShuffledBalance->GetHistNnn(a));\r
+      fListEventMixingBFS->Add(fShuffledBalance->GetHistNpp(a));\r
+      fListEventMixingBFS->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, fListEventMixingBF);\r
+  if(fRunShuffling) PostData(3, fListEventMixingBFS);\r
+  if(fUsePID) PostData(4, fHistListPIDQA);       //PID\r
+}\r
+\r
+//________________________________________________________________________\r
+void AliAnalysisTaskEventMixingBF::UserExec(Option_t *) {\r
+  // Main loop\r
+  // Called for each event\r
+  // NOTHING TO DO for event mixing!\r
+}      \r
+\r
+//________________________________________________________________________\r
+void  AliAnalysisTaskEventMixingBF::FinishTaskOutput(){\r
+  //Printf("END EventMixingBF");\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 AliAnalysisTaskEventMixingBF::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
+\r
+void AliAnalysisTaskEventMixingBF::UserExecMix(Option_t *)\r
+{\r
+\r
+  TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
+\r
+  TObjArray *array               = new TObjArray();\r
+\r
+  AliMixInputEventHandler *mixIEH = SetupEventsForMixing();\r
+\r
+  Int_t gNumberOfAcceptedTracks = 0;\r
+  Float_t fCentrality           = 0.;\r
+  \r
+  // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)\r
+  vector<Double_t> *chargeVector[9];          // original charge\r
+  for(Int_t i = 0; i < 9; i++){\r
+    chargeVector[i]        = new vector<Double_t>;\r
+  }\r
+  \r
+  Double_t v_charge;\r
+  Double_t v_y;\r
+  Double_t v_eta;\r
+  Double_t v_phi;\r
+  Double_t v_p[3];\r
+  Double_t v_pt;\r
+  Double_t v_E;\r
+\r
+  Int_t iMainTrackUsed = -1;\r
+\r
+  // -------------------------------------------------------------                  \r
+  // At the moment MIXING only for AODs\r
+  if(mixIEH){\r
+\r
+    //AOD analysis (vertex and track cuts also here!!!!)\r
+    if(gAnalysisLevel == "AOD") {\r
+      AliAODEvent* aodEventMain = dynamic_cast<AliAODEvent*>(fMainEvent); \r
+      if(!aodEventMain) {\r
+       Printf("ERROR: aodEventMain not available");\r
+       return;\r
+      }\r
+      AliAODEvent *aodEventMix  = dynamic_cast<AliAODEvent *>(fMixEvent); \r
+     if(!aodEventMix) {\r
+       Printf("ERROR: aodEventMix not available");\r
+       return;\r
+      }\r
+      \r
+      AliAODHeader *aodHeaderMain = aodEventMain->GetHeader();\r
+      AliAODHeader *aodHeaderMix  = aodEventMix->GetHeader();    \r
+  \r
+\r
+      // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
+      fHistEventStats->Fill(1); //all events\r
+\r
+      // Bool_t isSelectedMain = kTRUE;\r
+      // Bool_t isSelectedMix = kTRUE;\r
+\r
+      // if(fUseOfflineTrigger)\r
+      //       isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
+      // isSelectedMix = ((AliInputEventHandler*)((AliMultiInputEventHandler *)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetFirstMultiInputHandler())->IsEventSelected();\r
+      \r
+      // if(isSelectedMain && isSelectedMix) {\r
+      //       fHistEventStats->Fill(2); //triggered events\r
+       \r
+      //       //Centrality stuff (centrality in AOD header)\r
+      //       if(fUseCentrality) {\r
+      //         fCentrality = aodHeaderMain->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
+         \r
+      //         // QA for centrality estimators\r
+      //         fHistCentStats->Fill(0.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0M"));\r
+      //         fHistCentStats->Fill(1.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("FMD"));\r
+      //         fHistCentStats->Fill(2.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TRK"));\r
+      //         fHistCentStats->Fill(3.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKL"));\r
+      //         fHistCentStats->Fill(4.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL0"));\r
+      //         fHistCentStats->Fill(5.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL1"));\r
+      //         fHistCentStats->Fill(6.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
+      //         fHistCentStats->Fill(7.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
+      //         fHistCentStats->Fill(8.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
+         \r
+      //         // take only events inside centrality class\r
+      //         if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax)) \r
+      //           return;\r
+         \r
+      //         // centrality QA (V0M)\r
+      //         fHistV0M->Fill(aodEventMain->GetVZEROData()->GetMTotV0A(), aodEventMain->GetVZEROData()->GetMTotV0C());\r
+         \r
+      //         // centrality QA (reference tracks)\r
+      //         fHistRefTracks->Fill(0.,aodHeaderMain->GetRefMultiplicity());\r
+      //         fHistRefTracks->Fill(1.,aodHeaderMain->GetRefMultiplicityPos());\r
+      //         fHistRefTracks->Fill(2.,aodHeaderMain->GetRefMultiplicityNeg());\r
+      //         fHistRefTracks->Fill(3.,aodHeaderMain->GetTPConlyRefMultiplicity());\r
+      //         fHistRefTracks->Fill(4.,aodHeaderMain->GetNumberOfITSClusters(0));\r
+      //         fHistRefTracks->Fill(5.,aodHeaderMain->GetNumberOfITSClusters(1));\r
+      //         fHistRefTracks->Fill(6.,aodHeaderMain->GetNumberOfITSClusters(2));\r
+      //         fHistRefTracks->Fill(7.,aodHeaderMain->GetNumberOfITSClusters(3));\r
+      //         fHistRefTracks->Fill(8.,aodHeaderMain->GetNumberOfITSClusters(4));\r
+      //       }\r
+       \r
+      //       const AliAODVertex *vertexMain = aodEventMain->GetPrimaryVertex();\r
+      //       const AliAODVertex *vertexMix  = aodEventMix->GetPrimaryVertex();\r
+      //       \r
+      //       if(vertexMain && vertexMix) {\r
+      //         Double32_t fCovMain[6];\r
+      //         Double32_t fCovMix[6];\r
+      //         vertexMain->GetCovarianceMatrix(fCovMain);\r
+      //         vertexMix->GetCovarianceMatrix(fCovMix);\r
+         \r
+      //         if(vertexMain->GetNContributors() > 0 && vertexMix->GetNContributors() > 0) {\r
+      //           if(fCovMain[5] != 0 && fCovMix[5] != 0) {\r
+      //             fHistEventStats->Fill(3); //events with a proper vertex\r
+      //             if(TMath::Abs(vertexMain->GetX()) < fVxMax && TMath::Abs(vertexMix->GetX()) < fVxMax ) {\r
+      //               if(TMath::Abs(vertexMain->GetY()) < fVyMax && TMath::Abs(vertexMix->GetY()) < fVyMax) {\r
+      //                 if(TMath::Abs(vertexMain->GetZ()) < fVzMax && TMath::Abs(vertexMix->GetZ()) < fVzMax) {\r
+      //                   fHistEventStats->Fill(4); //analyzed events\r
+      //                   fHistVx->Fill(vertexMain->GetX());\r
+      //                   fHistVy->Fill(vertexMain->GetY());\r
+      //                   fHistVz->Fill(vertexMain->GetZ());\r
+\r
+                   // Loop over tracks in main event\r
+                   for (Int_t iTracksMain = 0; iTracksMain < aodEventMain->GetNumberOfTracks(); iTracksMain++) {\r
+                     AliAODTrack* aodTrackMain = dynamic_cast<AliAODTrack *>(aodEventMain->GetTrack(iTracksMain));\r
+                     if (!aodTrackMain) {\r
+                       Printf("ERROR: Could not receive track %d", iTracksMain);\r
+                       continue;\r
+                     }\r
+                     \r
+                     // AOD track cuts\r
+                     \r
+                     // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
+                     // take only TPC only tracks \r
+                     fHistTrackStats->Fill(aodTrackMain->GetFilterMap());\r
+                     if(!aodTrackMain->TestFilterBit(nAODtrackCutBit)) continue;\r
+                     \r
+                     v_charge = aodTrackMain->Charge();\r
+                     v_y      = aodTrackMain->Y();\r
+                     v_eta    = aodTrackMain->Eta();\r
+                     v_phi    = aodTrackMain->Phi() * TMath::RadToDeg();\r
+                     v_E      = aodTrackMain->E();\r
+                     v_pt     = aodTrackMain->Pt();\r
+                     aodTrackMain->PxPyPz(v_p);\r
+                     \r
+                     Float_t DCAxy = aodTrackMain->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
+                     Float_t DCAz  = aodTrackMain->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
+                     \r
+                     \r
+                     // Kinematics cuts from ESD track cuts\r
+                     if( v_pt < fPtMin || v_pt > fPtMax)      continue;\r
+                     if( v_eta < fEtaMin || v_eta > fEtaMax)  continue;\r
+                     \r
+                     // Extra DCA cuts (for systematic studies [!= -1])\r
+                     if( fDCAxyCut != -1 && 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 && aodTrackMain->Chi2perNDF() > fTPCchi2Cut){\r
+                       continue;\r
+                     }\r
+                     if( fNClustersTPCCut != -1 && aodTrackMain->GetTPCNcls() < fNClustersTPCCut){\r
+                       continue;\r
+                     }\r
+                     \r
+                     // fill QA histograms\r
+                     fHistClus->Fill(aodTrackMain->GetITSNcls(),aodTrackMain->GetTPCNcls());\r
+                     fHistDCA->Fill(DCAz,DCAxy);\r
+                     fHistChi2->Fill(aodTrackMain->Chi2perNDF());\r
+                     fHistPt->Fill(v_pt);\r
+                     fHistEta->Fill(v_eta);\r
+                     fHistPhi->Fill(v_phi);\r
+                     \r
+                     // fill charge vector\r
+                     chargeVector[0]->push_back(v_charge);\r
+                     chargeVector[1]->push_back(v_y);\r
+                     chargeVector[2]->push_back(v_eta);\r
+                     chargeVector[3]->push_back(v_phi);\r
+                     chargeVector[4]->push_back(v_p[0]);\r
+                     chargeVector[5]->push_back(v_p[1]);\r
+                     chargeVector[6]->push_back(v_p[2]);\r
+                     chargeVector[7]->push_back(v_pt);\r
+                     chargeVector[8]->push_back(v_E);\r
+\r
+                     // -------------------------------------------------------------               \r
+                     // for each track in main event loop over all tracks in mix event\r
+                     for (Int_t iTracksMix = 0; iTracksMix < aodEventMix->GetNumberOfTracks(); iTracksMix++) {\r
+\r
+                       AliAODTrack* aodTrackMix = dynamic_cast<AliAODTrack *>(aodEventMix->GetTrack(iTracksMix));\r
+                       if (!aodTrackMix) {\r
+                         Printf("ERROR: Could not receive track %d", iTracksMix);\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(aodTrackMix->GetFilterMap());\r
+                       if(!aodTrackMix->TestFilterBit(nAODtrackCutBit)) continue;\r
+                       \r
+                       v_charge = aodTrackMix->Charge();\r
+                       v_y      = aodTrackMix->Y();\r
+                       v_eta    = aodTrackMix->Eta();\r
+                       v_phi    = aodTrackMix->Phi() * TMath::RadToDeg();\r
+                       v_E      = aodTrackMix->E();\r
+                       v_pt     = aodTrackMix->Pt();\r
+                       aodTrackMix->PxPyPz(v_p);\r
+                     \r
+                       Float_t DCAxy = aodTrackMix->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
+                       Float_t DCAz  = aodTrackMix->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 && aodTrackMix->Chi2perNDF() > fTPCchi2Cut){\r
+                         continue;\r
+                       }\r
+                       if( fNClustersTPCCut != -1 && aodTrackMix->GetTPCNcls() < fNClustersTPCCut){\r
+                         continue;\r
+                       }\r
+                       \r
+                       // fill QA histograms\r
+                       fHistClus->Fill(aodTrackMix->GetITSNcls(),aodTrackMix->GetTPCNcls());\r
+                       fHistDCA->Fill(DCAz,DCAxy);\r
+                       fHistChi2->Fill(aodTrackMix->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
+                       \r
+                       \r
+                     } //mix track loop\r
+\r
+                     // calculate balance function for each track in main event\r
+                     iMainTrackUsed++; // is needed to do no double counting in Balance Function calculation   \r
+                     if(iMainTrackUsed >= chargeVector[0]->size()) break; //do not allow more tracks than in mixed event!\r
+                     fBalance->CalculateBalance(fCentrality,chargeVector,iMainTrackUsed);\r
+                     // clean charge vector afterwards\r
+                     for(Int_t i = 0; i < 9; i++){                    \r
+                       chargeVector[i]->clear();\r
+                     }\r
+                     \r
+\r
+                   } //main 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
+}\r
+\r
+AliMixInputEventHandler *AliAnalysisTaskEventMixingBF::SetupEventsForMixing() {\r
+\r
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
+  AliMultiInputEventHandler *inEvHMain = dynamic_cast<AliMultiInputEventHandler *>(mgr->GetInputEventHandler());\r
+  if (inEvHMain) {\r
+\r
+      AliMixInputEventHandler *mixEH = dynamic_cast<AliMixInputEventHandler *>(inEvHMain->GetFirstMultiInputHandler());\r
+      if (!mixEH) return kFALSE;\r
+      if (mixEH->CurrentBinIndex() < 0) {\r
+         AliDebug(AliLog::kDebug + 1, "Current event mixEH->CurrentEntry() == -1");\r
+         return kFALSE ;\r
+      }\r
+      AliDebug(AliLog::kDebug, Form("Mixing %lld %d [%lld,%lld] %d", mixEH->CurrentEntry(), mixEH->CurrentBinIndex(), mixEH->CurrentEntryMain(), mixEH->CurrentEntryMix(), mixEH->NumberMixed()));\r
+\r
+      AliInputEventHandler      *ihMainCurrent     = inEvHMain->GetFirstInputEventHandler();\r
+      fMainEvent = ihMainCurrent->GetEvent();\r
+\r
+      AliMultiInputEventHandler *inEvHMixedCurrent = mixEH->GetFirstMultiInputHandler(); // for buffer = 1\r
+      AliInputEventHandler      *ihMixedCurrent    = inEvHMixedCurrent->GetFirstInputEventHandler();\r
+      fMixEvent                                    = ihMixedCurrent->GetEvent();\r
+      \r
+      return mixEH;\r
+  }\r
+  return NULL;\r
+} \r
diff --git a/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEventMixingBF.h b/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEventMixingBF.h
new file mode 100755 (executable)
index 0000000..b0c06f5
--- /dev/null
@@ -0,0 +1,245 @@
+#ifndef ALIANALYSISTASKEVENTMIXINGBF_CXX\r
+#define ALIANALYSISTASKEVENTMIXINGBF_CXX\r
+\r
+// Analysis task for the EventMixingBF code\r
+// Authors: Panos Cristakoglou@cern.ch, m.weber@cern.ch\r
+\r
+class TList;\r
+class TH1F;\r
+class TH2F;\r
+class TF1;\r
+\r
+class AliMixInputEventHandler;\r
+class AliBalanceEventMixing;\r
+class AliESDtrackCuts;\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+#include "AliBalanceEventMixing.h"\r
+\r
+#include "AliPID.h"  \r
+#include "AliPIDResponse.h"\r
+#include "AliPIDCombined.h"\r
\r
+\r
+class AliAnalysisTaskEventMixingBF : public AliAnalysisTaskSE {\r
+ public:\r
+  AliAnalysisTaskEventMixingBF(const char *name = "AliAnalysisTaskEventMixingBF");\r
+  virtual ~AliAnalysisTaskEventMixingBF(); \r
+  \r
+  \r
+  virtual void   UserCreateOutputObjects();\r
+  virtual void   UserExec(Option_t *option);\r
+  virtual void   UserExecMix(Option_t*);\r
+  virtual void   FinishTaskOutput();\r
+  virtual void   Terminate(Option_t *);\r
+\r
+  void SetAnalysisObject(AliBalanceEventMixing *const analysis) {\r
+    fBalance         = analysis;\r
+    }\r
+  void SetShufflingObject(AliBalanceEventMixing *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
+  AliBalanceEventMixing *fBalance; //EventMixingBF object\r
+  Bool_t fRunShuffling;//run shuffling or not\r
+  AliBalanceEventMixing *fShuffledBalance; //EventMixingBF object (shuffled)\r
+  TList *fList; //fList object\r
+  TList *fListEventMixingBF; //fList object\r
+  TList *fListEventMixingBFS; //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
+  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
+  // Event Mixing\r
+  AliVEvent       *fMainEvent;\r
+  AliVEvent       *fMixEvent;\r
+\r
+  AliMixInputEventHandler *SetupEventsForMixing();  \r
+\r
+  AliAnalysisTaskEventMixingBF(const AliAnalysisTaskEventMixingBF&); // not implemented\r
+  AliAnalysisTaskEventMixingBF& operator=(const AliAnalysisTaskEventMixingBF&); // not implemented\r
+  \r
+  ClassDef(AliAnalysisTaskEventMixingBF, 1); // example of analysis\r
+};\r
+\r
+#endif\r
diff --git a/PWGCF/EBYE/BalanceFunctions/AliBalanceEventMixing.cxx b/PWGCF/EBYE/BalanceFunctions/AliBalanceEventMixing.cxx
new file mode 100644 (file)
index 0000000..bc3860e
--- /dev/null
@@ -0,0 +1,812 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+//           Balance Function class
+//   This is the class to deal with the Balance Function analysis
+//   Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
+//   Modified: Michael Weber, m.weber@cern.ch
+//-----------------------------------------------------------------
+
+
+//ROOT
+#include <Riostream.h>
+#include <TMath.h>
+#include <TAxis.h>
+#include <TH2D.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 "AliBalance.h"
+#include "AliBalanceEventMixing.h"
+
+ClassImp(AliBalanceEventMixing)
+
+//____________________________________________________________________//
+AliBalanceEventMixing::AliBalanceEventMixing() :
+  TObject(), 
+  bShuffle(kFALSE),
+  fAnalysisLevel("ESD"),
+  fAnalyzedEvents(0) ,
+  fCentralityId(0) ,
+  fCentStart(0.),
+  fCentStop(0.)
+{
+  // 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.;
+
+    fNn[i] = 0.0;
+    fNp[i] = 0.0;
+
+    for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
+      fNpp[i][j] = .0;
+      fNnn[i][j] = .0;
+      fNpn[i][j] = .0;
+      fNnp[i][j] = .0;
+      fB[i][j] = 0.0;
+      ferror[i][j] = 0.0;
+    }
+
+    fHistP[i]  = NULL;
+    fHistN[i]  = NULL;
+    fHistPP[i] = NULL;
+    fHistPN[i] = NULL;
+    fHistNP[i] = NULL;
+    fHistNN[i] = NULL;
+
+  }
+}
+
+
+//____________________________________________________________________//
+AliBalanceEventMixing::AliBalanceEventMixing(const AliBalanceEventMixing& balance):
+  TObject(balance), bShuffle(balance.bShuffle), 
+  fAnalysisLevel(balance.fAnalysisLevel),
+  fAnalyzedEvents(balance.fAnalyzedEvents), 
+  fCentralityId(balance.fCentralityId),
+  fCentStart(balance.fCentStart),
+  fCentStop(balance.fCentStop) {
+  //copy constructor
+  for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
+    fNn[i] = balance.fNn[i];
+    fNp[i] = balance.fNp[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 j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
+      fNpp[i][j] = .0;
+      fNnn[i][j] = .0;
+      fNpn[i][j] = .0;
+      fNnp[i][j] = .0;
+      fB[i][j] = 0.0;
+      ferror[i][j] = 0.0;
+    } 
+  }
+ }
+
+//____________________________________________________________________//
+AliBalanceEventMixing::~AliBalanceEventMixing() {
+  // 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 AliBalanceEventMixing::SetInterval(Int_t iAnalysisType,
+                            Double_t p1Start, Double_t p1Stop,
+                            Int_t ibins, Double_t p2Start, Double_t p2Stop) {
+  // Sets the analyzed interval. 
+  // Set the same Information for all analyses
+
+  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 AliBalanceEventMixing::InitHistograms() {
+  //Initialize the histograms
+  TString histName;
+  for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
+    histName = "fHistP"; histName += gBFAnalysisType[iAnalysisType]; 
+    if(bShuffle) histName.Append("_shuffle");
+    if(fCentralityId) histName += fCentralityId.Data();
+    fHistP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
+
+    histName = "fHistN"; histName += gBFAnalysisType[iAnalysisType]; 
+    if(bShuffle) histName.Append("_shuffle");
+    if(fCentralityId) histName += fCentralityId.Data();
+    fHistN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
+  
+    histName = "fHistPN"; histName += gBFAnalysisType[iAnalysisType]; 
+    if(bShuffle) histName.Append("_shuffle");
+    if(fCentralityId) histName += fCentralityId.Data();
+    fHistPN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
+    
+    histName = "fHistNP"; histName += gBFAnalysisType[iAnalysisType]; 
+    if(bShuffle) histName.Append("_shuffle");
+    if(fCentralityId) histName += fCentralityId.Data();
+    fHistNP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
+    
+    histName = "fHistPP"; histName += gBFAnalysisType[iAnalysisType]; 
+    if(bShuffle) histName.Append("_shuffle");
+    if(fCentralityId) histName += fCentralityId.Data();
+    fHistPP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
+    
+    histName = "fHistNN"; histName += gBFAnalysisType[iAnalysisType]; 
+    if(bShuffle) histName.Append("_shuffle");
+    if(fCentralityId) histName += fCentralityId.Data();
+    fHistNN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
+  }
+}
+
+//____________________________________________________________________//
+void AliBalanceEventMixing::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 AliBalanceEventMixing::CalculateBalance(Float_t fCentrality,vector<Double_t> **chargeVector,Int_t iMainTrack) {
+  // Calculates the balance function
+  // For the event mixing only for all combinations of the first track (main event) with all other tracks (mix event)
+  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("(AliBalanceEventMixing) Number of tracks: %d",gNtrack);
+
+  for(i = 0; i < gNtrack;i++){
+
+      // for event mixing: only store the track from the main event
+      if(iMainTrack > -1){ 
+       if(i>0) break;
+      }
+      
+      Short_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);
+      
+      //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] += 1.;
+             fHistP[iAnalysisType]->Fill(fCentrality,pseudorapidity);
+           }//charge > 0
+           if(charge < 0) {
+             fNn[iAnalysisType] += 1.;
+             fHistN[iAnalysisType]->Fill(fCentrality,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] += 1.;
+             fHistP[iAnalysisType]->Fill(fCentrality,phi);
+           }//charge > 0
+           if(charge < 0) {
+             fNn[iAnalysisType] += 1.;
+             fHistN[iAnalysisType]->Fill(fCentrality,phi);
+           }//charge < 0
+         }//p1 interval check
+       }//analysis type: phi
+       else {
+         if((rapidity >= fP1Start[iAnalysisType]) && (rapidity <= fP1Stop[iAnalysisType])) {
+           if(charge > 0) {
+             fNp[iAnalysisType] += 1.;
+             fHistP[iAnalysisType]->Fill(fCentrality,rapidity);
+           }//charge > 0
+           if(charge < 0) {
+             fNn[iAnalysisType] += 1.;
+             fHistN[iAnalysisType]->Fill(fCentrality,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.;
+
+  Short_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 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);
+
+     for(j = 0; j < i; j++) {
+
+      // for event mixing: only store the track from the main event
+      if(iMainTrack > -1){ 
+       if(j>0) break;
+      }
+
+      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);
+    
+       // filling the arrays
+
+       // RAPIDITY 
+        dy = TMath::Abs(rap1 - rap2);
+
+       // Eta
+       deta = TMath::Abs(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 = TMath::Abs(phi1 - phi2);
+       if(dphi>180) dphi = 360 - dphi;  //dphi should be between 0 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][iBin] += 1.;
+               fHistPP[kRapidity]->Fill(fCentrality,dy);
+             }
+             else if((charge1 < 0)&&(charge2 < 0)) {
+               fNnn[kRapidity][iBin] += 1.;
+               fHistNN[kRapidity]->Fill(fCentrality,dy);
+             }
+             else if((charge1 > 0)&&(charge2 < 0)) {
+               fNpn[kRapidity][iBin] += 1.;
+               fHistPN[kRapidity]->Fill(fCentrality,dy);
+             }
+             else if((charge1 < 0)&&(charge2 > 0)) {
+               fNpn[kRapidity][iBin] += 1.;
+                   fHistPN[kRapidity]->Fill(fCentrality,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][iBin] += 1.;
+               fHistPP[kEta]->Fill(fCentrality,deta);
+             }
+             if((charge1 < 0)&&(charge2 < 0)) {
+               fNnn[kEta][iBin] += 1.;
+                   fHistNN[kEta]->Fill(fCentrality,deta);
+             }
+             if((charge1 > 0)&&(charge2 < 0)) {
+               fNpn[kEta][iBin] += 1.;
+               fHistPN[kEta]->Fill(fCentrality,deta);
+             }
+             if((charge1 < 0)&&(charge2 > 0)) {
+               fNpn[kEta][iBin] += 1.;
+                   fHistPN[kEta]->Fill(fCentrality,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][iBin] += 1.;
+               fHistPP[kQlong]->Fill(fCentrality,qLong);
+             }
+             if((charge1 < 0)&&(charge2 < 0)) {
+               fNnn[kQlong][iBin] += 1.;
+               fHistNN[kQlong]->Fill(fCentrality,qLong);
+             }
+             if((charge1 > 0)&&(charge2 < 0)) {
+               fNpn[kQlong][iBin] += 1.;
+               fHistPN[kQlong]->Fill(fCentrality,qLong);
+             }
+             if((charge1 < 0)&&(charge2 > 0)) {
+               fNpn[kQlong][iBin] += 1.;
+               fHistPN[kQlong]->Fill(fCentrality,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][iBin] += 1.;
+               fHistPP[kQout]->Fill(fCentrality,qOut);
+                 }
+             if((charge1 < 0)&&(charge2 < 0)) {
+               fNnn[kQout][iBin] += 1.;
+               fHistNN[kQout]->Fill(fCentrality,qOut);
+             }
+             if((charge1 > 0)&&(charge2 < 0)) {
+               fNpn[kQout][iBin] += 1.;
+               fHistPN[kQout]->Fill(fCentrality,qOut);
+             }
+             if((charge1 < 0)&&(charge2 > 0)) {
+               fNpn[kQout][iBin] += 1.;
+               fHistPN[kQout]->Fill(fCentrality,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][iBin] += 1.;
+               fHistPP[kQside]->Fill(fCentrality,qSide);
+             }
+             if((charge1 < 0)&&(charge2 < 0)) {
+               fNnn[kQside][iBin] += 1.;
+               fHistNN[kQside]->Fill(fCentrality,qSide);
+             }
+             if((charge1 > 0)&&(charge2 < 0)) {
+               fNpn[kQside][iBin] += 1.;
+               fHistPN[kQside]->Fill(fCentrality,qSide);
+             }
+             if((charge1 < 0)&&(charge2 > 0)) {
+               fNpn[kQside][iBin] += 1.;
+               fHistPN[kQside]->Fill(fCentrality,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][iBin] += 1.;
+               fHistPP[kQinv]->Fill(fCentrality,qInv);
+             }
+             if((charge1 < 0)&&(charge2 < 0)) {
+               fNnn[kQinv][iBin] += 1.;
+               fHistNN[kQinv]->Fill(fCentrality,qInv);
+             }
+             if((charge1 > 0)&&(charge2 < 0)) {
+               fNpn[kQinv][iBin] += 1.;
+               fHistPN[kQinv]->Fill(fCentrality,qInv);
+             }
+             if((charge1 < 0)&&(charge2 > 0)) {
+               fNpn[kQinv][iBin] += 1.;
+               fHistPN[kQinv]->Fill(fCentrality,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][iBin] += 1.;
+               fHistPP[kPhi]->Fill(fCentrality,dphi);
+             }
+             if((charge1 < 0)&&(charge2 < 0)) {
+               fNnn[kPhi][iBin] += 1.;
+               fHistNN[kPhi]->Fill(fCentrality,dphi);
+             }
+             if((charge1 > 0)&&(charge2 < 0)) {
+               fNpn[kPhi][iBin] += 1.;
+               fHistPN[kPhi]->Fill(fCentrality,dphi);
+             }
+             if((charge1 < 0)&&(charge2 > 0)) {
+               fNpn[kPhi][iBin] += 1.;
+               fHistPN[kPhi]->Fill(fCentrality,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 AliBalanceEventMixing::GetBalance(Int_t iAnalysisType, Int_t p2) {
+  // Returns the value of the balance function in bin p2
+  fB[iAnalysisType][p2] = 0.5*(((fNpn[iAnalysisType][p2] - 2.*fNnn[iAnalysisType][p2])/fNn[iAnalysisType]) + ((fNpn[iAnalysisType][p2] - 2.*fNpp[iAnalysisType][p2])/fNp[iAnalysisType]))/fP2Step[iAnalysisType];
+  
+  return fB[iAnalysisType][p2];
+}
+    
+//____________________________________________________________________//
+Double_t AliBalanceEventMixing::GetError(Int_t iAnalysisType, 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][p2] = TMath::Sqrt(Double_t(fNpp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType]))
+                             + Double_t(fNnn[iAnalysisType][p2])/(Double_t(fNn[iAnalysisType])*Double_t(fNn[iAnalysisType]))
+                             + Double_t(fNpn[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType])) 
+                             + Double_t(fNnp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType]))
+                             //+ TMath::Power(fNpn[iAnalysisType][p2]-fNpp[iAnalysisType][p2],2)/TMath::Power(Double_t(fNp[iAnalysisType]),3)
+                             //+ TMath::Power(fNnp[iAnalysisType][p2]-fNnn[iAnalysisType][p2],2)/TMath::Power(Double_t(fNn[iAnalysisType]),3) 
+                              ) /fP2Step[iAnalysisType];*/
+
+  ferror[iAnalysisType][p2] = TMath::Sqrt( Double_t(fNpp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType])) + 
+                                          Double_t(fNnn[iAnalysisType][p2])/(Double_t(fNn[iAnalysisType])*Double_t(fNn[iAnalysisType])) + 
+                                          Double_t(fNpn[iAnalysisType][p2])*TMath::Power((0.5/Double_t(fNp[iAnalysisType]) + 0.5/Double_t(fNn[iAnalysisType])),2))/fP2Step[iAnalysisType];
+  
+  return ferror[iAnalysisType][p2];
+}
+//____________________________________________________________________//
+TGraphErrors *AliBalanceEventMixing::DrawBalance(Int_t iAnalysisType) {
+
+  // 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] == 0)||(fNn[iAnalysisType] == 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,i);
+    ber[i] = GetError(iAnalysisType,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 AliBalanceEventMixing::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: "<<gBFAnalysisType[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 *AliBalanceEventMixing::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow) {
+  //Returns the BF histogram, extracted from the 6 TH2D objects 
+  //(private members) of the AliBalanceEventMixing class.
+  //
+  // Acceptance correction: 
+  // - only for analysis type = kEta
+  // - only if etaWindow > 0 (default = -1.)
+  // - calculated as proposed by STAR 
+  //
+  TString gAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"};
+  TString histName = "gHistBalanceFunctionHistogram";
+  histName += gAnalysisType[iAnalysisType];
+
+  SetInterval(iAnalysisType, fHistP[iAnalysisType]->GetYaxis()->GetXmin(),
+             fHistP[iAnalysisType]->GetYaxis()->GetXmin(),
+             fHistPP[iAnalysisType]->GetNbinsY(),
+             fHistPP[iAnalysisType]->GetYaxis()->GetXmin(),
+             fHistPP[iAnalysisType]->GetYaxis()->GetXmax());
+
+  // determine the projection thresholds
+  Int_t binMinX, binMinY, binMinZ;
+  Int_t binMaxX, binMaxY, binMaxZ;
+
+  fHistPP[iAnalysisType]->GetBinXYZ(fHistPP[iAnalysisType]->FindBin(centrMin),binMinX,binMinY,binMinZ);
+  fHistPP[iAnalysisType]->GetBinXYZ(fHistPP[iAnalysisType]->FindBin(centrMax),binMaxX,binMaxY,binMaxZ);
+
+  TH1D *gHistBalanceFunctionHistogram = new TH1D(histName.Data(),"",fHistPP[iAnalysisType]->GetNbinsY(),fHistPP[iAnalysisType]->GetYaxis()->GetXmin(),fHistPP[iAnalysisType]->GetYaxis()->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]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistPN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
+  TH1D *hTemp2 = dynamic_cast<TH1D *>(fHistPN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f_copy",fHistPN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
+  TH1D *hTemp3 = dynamic_cast<TH1D *>(fHistNN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistNN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
+  TH1D *hTemp4 = dynamic_cast<TH1D *>(fHistPP[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistPP[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
+  TH1D *hTemp5 = dynamic_cast<TH1D *>(fHistN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
+  TH1D *hTemp6 = dynamic_cast<TH1D *>(fHistP[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistP[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
+
+  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]);
+  }
+
+  // do the acceptance correction (only for Eta and etaWindow > 0)
+  if(iAnalysisType == kEta && etaWindow > 0){
+    for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){
+      
+      Double_t notCorrected = gHistBalanceFunctionHistogram->GetBinContent(iBin+1);
+      Double_t corrected    = notCorrected / (1 - (gHistBalanceFunctionHistogram->GetBinCenter(iBin+1))/ etaWindow );
+      gHistBalanceFunctionHistogram->SetBinContent(iBin+1, corrected);
+      
+    }
+  }
+  
+  PrintResults(iAnalysisType,gHistBalanceFunctionHistogram);
+
+  return gHistBalanceFunctionHistogram;
+}
diff --git a/PWGCF/EBYE/BalanceFunctions/AliBalanceEventMixing.h b/PWGCF/EBYE/BalanceFunctions/AliBalanceEventMixing.h
new file mode 100644 (file)
index 0000000..5ff4481
--- /dev/null
@@ -0,0 +1,147 @@
+#ifndef ALIBALANCEEVENTMIXING_H
+#define ALIBALANCEEVENTMIXING_H
+/*  See cxx source for full Copyright notice */
+
+
+//-------------------------------------------------------------------------
+//                          Class AliBalanceEventMixing
+//   This is the class for the Balance Function analysis
+//
+//    Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
+//    Modified: Michael Weber, m.weber@cern.ch
+//-------------------------------------------------------------------------
+
+#include <TObject.h>
+#include "TString.h"
+
+#define ANALYSIS_TYPES 7
+#define MAXIMUM_NUMBER_OF_STEPS        1024
+
+class TGraphErrors;
+class TH1D;
+class TH2D;
+
+class AliBalanceEventMixing : public TObject {
+ public:
+  enum EAnalysisType {
+    kRapidity = 0, 
+    kEta      = 1, 
+    kQlong    = 2, 
+    kQout     = 3, 
+    kQside    = 4, 
+    kQinv     = 5, 
+    kPhi      = 6 
+  };
+
+  AliBalanceEventMixing();
+  AliBalanceEventMixing(const AliBalanceEventMixing& balance);
+  ~AliBalanceEventMixing();
+
+  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);
+  void SetCentralityInterval(Double_t cStart, Double_t cStop)  { fCentStart = cStart; fCentStop = cStop;};
+  void SetNp(Int_t analysisType, Double_t NpSet)   { fNp[analysisType] = NpSet; }
+  void SetNn(Int_t analysisType, Double_t NnSet)   { fNn[analysisType] = NnSet; }
+  void SetNpp(Int_t analysisType, Int_t ibin, Double_t NppSet) { if(ibin > -1 && ibin < MAXIMUM_NUMBER_OF_STEPS) fNpp[analysisType][ibin] = NppSet; }
+  void SetNpn(Int_t analysisType, Int_t ibin, Double_t NpnSet) { if(ibin > -1 && ibin < MAXIMUM_NUMBER_OF_STEPS) fNpn[analysisType][ibin] = NpnSet; }
+  void SetNnp(Int_t analysisType, Int_t ibin, Double_t NnpSet) { if(ibin > -1 && ibin < MAXIMUM_NUMBER_OF_STEPS) fNnp[analysisType][ibin] = NnpSet; }
+  void SetNnn(Int_t analysisType, Int_t ibin, Double_t NnnSet) { if(ibin > -1 && ibin < MAXIMUM_NUMBER_OF_STEPS) fNnn[analysisType][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) const { return 1.0*fNp[analysisType]; }
+  Double_t GetNn(Int_t analysisType) const { return 1.0*fNn[analysisType]; }
+  Double_t GetNnn(Int_t analysisType, Int_t p2) const { 
+    return 1.0*fNnn[analysisType][p2]; }
+  Double_t GetNpp(Int_t analysisType, Int_t p2) const { 
+    return 1.0*fNpp[analysisType][p2]; }
+  Double_t GetNpn(Int_t analysisType, Int_t p2) const { 
+    return 1.0*fNpn[analysisType][p2]; }  
+  Double_t GetNnp(Int_t analysisType, Int_t p2) const { 
+    return 1.0*fNnp[analysisType][p2]; }
+
+  void CalculateBalance(Float_t fCentrality, vector<Double_t> **chargeVector,Int_t iMainTrack = -1);
+  
+  Double_t GetBalance(Int_t a, Int_t p2);
+  Double_t GetError(Int_t a, Int_t p2);
+
+  TH2D *GetHistNp(Int_t iAnalysisType) { return fHistP[iAnalysisType];}
+  TH2D *GetHistNn(Int_t iAnalysisType) { return fHistN[iAnalysisType];}
+  TH2D *GetHistNpn(Int_t iAnalysisType) { return fHistPN[iAnalysisType];}
+  TH2D *GetHistNnp(Int_t iAnalysisType) { return fHistNP[iAnalysisType];}
+  TH2D *GetHistNpp(Int_t iAnalysisType) { return fHistPP[iAnalysisType];}
+  TH2D *GetHistNnn(Int_t iAnalysisType) { return fHistNN[iAnalysisType];}
+
+  void PrintAnalysisSettings();
+  TGraphErrors *DrawBalance(Int_t fAnalysisType);
+
+  void SetHistNp(Int_t iAnalysisType, TH2D *gHist) { 
+    fHistP[iAnalysisType] = gHist;}
+  void SetHistNn(Int_t iAnalysisType, TH2D *gHist) { 
+    fHistN[iAnalysisType] = gHist;}
+  void SetHistNpn(Int_t iAnalysisType, TH2D *gHist) { 
+    fHistPN[iAnalysisType] = gHist;}
+  void SetHistNnp(Int_t iAnalysisType, TH2D *gHist) { 
+    fHistNP[iAnalysisType] = gHist;}
+  void SetHistNpp(Int_t iAnalysisType, TH2D *gHist) { 
+    fHistPP[iAnalysisType] = gHist;}
+  void SetHistNnn(Int_t iAnalysisType, TH2D *gHist) { 
+    fHistNN[iAnalysisType] = gHist;}
+
+  TH1D *GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow = -1);
+  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_NUMBER_OF_STEPS]; //N(--)
+  Double_t fNpp[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //N(++)
+  Double_t fNpn[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //N(+-)
+  Double_t fNnp[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //N(-+)
+  Double_t fNn[ANALYSIS_TYPES], fNp[ANALYSIS_TYPES]; //number of pos./neg. inside the analyzed interval
+  
+  Double_t fB[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //BF matrix
+  Double_t ferror[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //error of the BF
+  
+  TH2D *fHistP[ANALYSIS_TYPES]; //N+
+  TH2D *fHistN[ANALYSIS_TYPES]; //N-
+  TH2D *fHistPN[ANALYSIS_TYPES]; //N+-
+  TH2D *fHistNP[ANALYSIS_TYPES]; //N-+
+  TH2D *fHistPP[ANALYSIS_TYPES]; //N++
+  TH2D *fHistNN[ANALYSIS_TYPES]; //N--
+
+  AliBalanceEventMixing & operator=(const AliBalanceEventMixing & ) {return *this;}
+
+  ClassDef(AliBalanceEventMixing, 3)
+};
+
+#endif
diff --git a/PWGCF/EBYE/macros/AddMixingHandler.C b/PWGCF/EBYE/macros/AddMixingHandler.C
new file mode 100644 (file)
index 0000000..816a21c
--- /dev/null
@@ -0,0 +1,39 @@
+#ifndef __CINT__//|
+#include <AliAnalysisManager.h>//|
+#include <AliMultiInputEventHandler.h>//|
+#include <EventMixing/EventMixing/AliMixEventPool.h>//|
+#include <EventMixing/EventMixing/AliMixEventCutObj.h>//|
+#include <EventMixing/EventMixing/AliMixInputEventHandler.h>//|
+#include <AliVEvent.h>//|
+#endif//|
+
+void AddMixingHandler(Double_t centMin = 70, Double_t centMax = 80, Double_t centStep = 2, AliMultiInputEventHandler* multiInputHandler, Bool_t useMC = kFALSE, Bool_t usePhysSel = kFALSE,TString opts = "")
+{
+
+   if (!multiInputHandler) return;
+
+   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+   const Int_t bufferSize = 1;
+   const Int_t mixNum = 5;
+   AliMixInputEventHandler *mixHandler = new AliMixInputEventHandler(bufferSize, mixNum);
+   mixHandler->SetInputHandlerForMixing(dynamic_cast<AliMultiInputEventHandler*>(mgr->GetInputEventHandler()));
+   AliMixEventPool *evPool = new AliMixEventPool();
+
+   //AliMixEventCutObj *multi = new AliMixEventCutObj(AliMixEventCutObj::kMultiplicity, 2, 10002, 10000);
+   AliMixEventCutObj *zvertex = new AliMixEventCutObj(AliMixEventCutObj::kZVertex, -10, 10, 5);
+
+   AliMixEventCutObj *centrality = new AliMixEventCutObj(AliMixEventCutObj::kCentrality, centMin, centMax, centStep, "V0M");
+
+   evPool->AddCut(centrality);
+   //evPool->AddCut(multi);
+   evPool->AddCut(zvertex);
+
+   // adds event pool (comment it and u will have default mixing)
+   mixHandler->SetEventPool(evPool);
+
+   // only use events with physics selection
+   if (usePhysSel) mixHandler->SelectCollisionCandidates(AliVEvent::kMB);
+
+   multiInputHandler->AddInputEventHandler(mixHandler);
+
+}
diff --git a/PWGCF/EBYE/macros/AddTaskBalanceEventMixing.C b/PWGCF/EBYE/macros/AddTaskBalanceEventMixing.C
new file mode 100644 (file)
index 0000000..a77fae2
--- /dev/null
@@ -0,0 +1,182 @@
+// 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
+AliAnalysisTaskEventMixingBF *AddTaskBalanceEventMixing(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
+\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("AddTaskEventMixingBF", "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("AddTaskEventMixingBF", "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 EventMixingBF configuration\r
+  //gROOT->LoadMacro("./configBalanceFunctionAnalysisEventMixing.C");\r
+  gROOT->LoadMacro("$ALICE_ROOT/PWGCF/EBYE/macros/configBalanceFunctionAnalysisEventMixing.C");\r
+  AliBalanceEventMixing *bf  = 0;  // Balance Function object\r
+  AliBalanceEventMixing *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("AddTaskEventMixingBF", "analysis type NOT known.");\r
+    return NULL;\r
+  }\r
+\r
+  // Create the task, add it to manager and configure it.\r
+  //===========================================================================\r
+  AliAnalysisTaskEventMixingBF *taskEventMixingBF = new AliAnalysisTaskEventMixingBF("TaskEventMixingBF");\r
+  taskEventMixingBF->SetAnalysisObject(bf);\r
+  if(gRunShuffling) taskEventMixingBF->SetShufflingObject(bfs);\r
+\r
+  taskEventMixingBF->SetCentralityPercentileRange(centrMin,centrMax);\r
+  if(analysisType == "ESD") {\r
+    AliESDtrackCuts *trackCuts = GetTrackCutsObject(ptMin,ptMax,etaMin,etaMax,maxTPCchi2,DCAxy,DCAz,minNClustersTPC);\r
+    taskEventMixingBF->SetAnalysisCutObject(trackCuts);\r
+    if(kUsePID) {\r
+      if(kUseBayesianPID)\r
+       taskEventMixingBF->SetUseBayesianPID(gMinAcceptedProbability);\r
+      else if(kUseNSigmaPID)\r
+       taskEventMixingBF->SetUseNSigmaPID(nSigmaMax);\r
+      taskEventMixingBF->SetParticleOfInterest(AliAnalysistaskEventMixingBF::kProton);\r
+      taskEventMixingBF->SetDetectorUsedForPID(AliAnalysisTaskEventMixingBF::kTOFpid);\r
+    }\r
+  }\r
+  else if(analysisType == "AOD") {\r
+    // pt and eta cut (pt_min, pt_max, eta_min, eta_max)\r
+    taskEventMixingBF->SetAODtrackCutBit(AODfilterBit);\r
+    taskEventMixingBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax);\r
+\r
+    // set extra DCA cuts (-1 no extra cut)\r
+    taskEventMixingBF->SetExtraDCACutsAOD(DCAxy,DCAz);\r
+\r
+    // set extra TPC chi2 / nr of clusters cut\r
+    taskEventMixingBF->SetExtraTPCCutsAOD(maxTPCchi2, minNClustersTPC);\r
+    \r
+  }\r
+  else if(analysisType == "MC") {\r
+    taskEventMixingBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax); \r
+  }\r
+\r
+  // offline trigger selection (AliVEvent.h)\r
+  // taskEventMixingBF->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) taskEventMixingBF->SelectCollisionCandidates(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral);\r
+  else                taskEventMixingBF->SelectCollisionCandidates(AliVEvent::kMB);\r
+\r
+  // centrality estimator (default = V0M)\r
+  taskEventMixingBF->SetCentralityEstimator(centralityEstimator);\r
+  \r
+  // vertex cut (x,y,z)\r
+  taskEventMixingBF->SetVertexDiamond(.3,.3,vertexZ);\r
+  \r
+\r
+\r
+  //bf->PrintAnalysisSettings();\r
+  mgr->AddTask(taskEventMixingBF);\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 *coutEventMixingBF = mgr->CreateContainer(Form("listEventMixingBF_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  if(gRunShuffling) AliAnalysisDataContainer *coutEventMixingBFS = mgr->CreateContainer(Form("listEventMixingBFShuffled_%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(taskEventMixingBF, 0, mgr->GetCommonInputContainer());\r
+  mgr->ConnectOutput(taskEventMixingBF, 1, coutQA);\r
+  mgr->ConnectOutput(taskEventMixingBF, 2, coutEventMixingBF);\r
+  if(gRunShuffling) mgr->ConnectOutput(taskEventMixingBF, 3, coutEventMixingBFS);\r
+  if(kUsePID && analysisType == "ESD") mgr->ConnectOutput(taskEventMixingBF, 4, coutQAPID);\r
+\r
+  return taskEventMixingBF;\r
+}\r
diff --git a/PWGCF/EBYE/macros/configBalanceFunctionAnalysisEventMixing.C b/PWGCF/EBYE/macros/configBalanceFunctionAnalysisEventMixing.C
new file mode 100644 (file)
index 0000000..c575724
--- /dev/null
@@ -0,0 +1,63 @@
+//__________________________________________________//\r
+AliBalanceEventMixing *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 AliBalanceEventMixing object and return it\r
+  AliBalanceEventMixing *gBalance = new AliBalanceEventMixing();\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(AliBalanceEventMixing::kRapidity,-0.8,0.8,64,0.0,1.6);  \r
+  //Eta\r
+  gBalance->SetInterval(AliBalanceEventMixing::kEta,-0.8,0.8,64,0.0,1.6);\r
+  //Qlong\r
+  gBalance->SetInterval(AliBalanceEventMixing::kQlong,-1,1,100,0.0,2.0);\r
+  //Qout\r
+  gBalance->SetInterval(AliBalanceEventMixing::kQout,-1,1,100,0.0,2.0);\r
+  //Qside\r
+  gBalance->SetInterval(AliBalanceEventMixing::kQside,-1,1,100,0.0,2.0);\r
+  //Qinv\r
+  gBalance->SetInterval(AliBalanceEventMixing::kQinv,-1,1,100,0.0,2.0);\r
+  //Phi\r
+  gBalance->SetInterval(AliBalanceEventMixing::kPhi,0.,360.,90,0.,180.0);\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
index 1c2ad3c..b3eedcd 100644 (file)
@@ -6,8 +6,10 @@
 
 #pragma link C++ class AliBalance+;
 #pragma link C++ class AliBalancePsi+;
+#pragma link C++ class AliBalanceEventMixing+;
 #pragma link C++ class AliAnalysisTaskBF+;
 #pragma link C++ class AliAnalysisTaskBFPsi+;
+#pragma link C++ class AliAnalysisTaskEventMixingBF+;
 
 #pragma link C++ class AliAnalysisTaskToyModel+;