Triggered Balance Function analysis (multidimensions)
authormiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 30 Apr 2012 12:43:18 +0000 (12:43 +0000)
committermiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 30 Apr 2012 12:43:18 +0000 (12:43 +0000)
PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskTriggeredBF.cxx [new file with mode: 0755]
PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskTriggeredBF.h [new file with mode: 0755]
PWGCF/EBYE/BalanceFunctions/AliBalanceEventMixing.h
PWGCF/EBYE/BalanceFunctions/AliBalanceTriggered.cxx [new file with mode: 0644]
PWGCF/EBYE/BalanceFunctions/AliBalanceTriggered.h [new file with mode: 0644]
PWGCF/EBYE/macros/AddTaskBalanceTriggered.C [new file with mode: 0644]

diff --git a/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskTriggeredBF.cxx b/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskTriggeredBF.cxx
new file mode 100755 (executable)
index 0000000..e781e7c
--- /dev/null
@@ -0,0 +1,460 @@
+#include "TChain.h"\r
+#include "TList.h"\r
+#include "TCanvas.h"\r
+#include "TLorentzVector.h"\r
+#include "TGraphErrors.h"\r
+#include "TH1F.h"\r
+#include "TH2F.h"\r
+#include "TArrayF.h"\r
+#include "TF1.h"\r
+#include "TRandom.h"\r
+\r
+#include "AliLog.h"\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+#include "AliAnalysisManager.h"\r
+\r
+#include "AliESDVertex.h"\r
+#include "AliESDEvent.h"\r
+#include "AliESDInputHandler.h"\r
+#include "AliAODEvent.h"\r
+#include "AliAODTrack.h"\r
+#include "AliAODInputHandler.h"\r
+#include "AliGenEventHeader.h"\r
+#include "AliGenHijingEventHeader.h"\r
+#include "AliMCEventHandler.h"\r
+#include "AliMCEvent.h"\r
+#include "AliMixInputEventHandler.h"\r
+#include "AliStack.h"\r
+\r
+#include "TH2D.h"    \r
+#include "AliTHn.h"              \r
+\r
+#include "AliAnalysisTaskTriggeredBF.h"\r
+#include "AliBalanceTriggered.h"\r
+\r
+\r
+// Analysis task for the TriggeredBF code\r
+// Authors: Panos.Christakoglou@nikhef.nl, m.weber@cern.ch\r
+\r
+ClassImp(AliAnalysisTaskTriggeredBF)\r
+\r
+//________________________________________________________________________\r
+AliAnalysisTaskTriggeredBF::AliAnalysisTaskTriggeredBF(const char *name) \r
+: AliAnalysisTaskSE(name), \r
+  fBalance(0),\r
+  fRunShuffling(kFALSE),\r
+  fShuffledBalance(0),\r
+  fList(0),\r
+  fListTriggeredBF(0),\r
+  fListTriggeredBFS(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
+  fCentralityEstimator("V0M"),\r
+  fUseCentrality(kFALSE),\r
+  fCentralityPercentileMin(0.), \r
+  fCentralityPercentileMax(5.),\r
+  fImpactParameterMin(0.),\r
+  fImpactParameterMax(20.),\r
+  fUseMultiplicity(kFALSE),\r
+  fNumberOfAcceptedTracksMin(0),\r
+  fNumberOfAcceptedTracksMax(10000),\r
+  fHistNumberOfAcceptedTracks(0),\r
+  fUseOfflineTrigger(kFALSE),\r
+  fVxMax(0.3),\r
+  fVyMax(0.3),\r
+  fVzMax(10.),\r
+  nAODtrackCutBit(128),\r
+  fPtMin(0.3),\r
+  fPtMax(1.5),\r
+  fEtaMin(-0.8),\r
+  fEtaMax(-0.8),\r
+  fDCAxyCut(-1),\r
+  fDCAzCut(-1),\r
+  fTPCchi2Cut(-1),\r
+  fNClustersTPCCut(-1)\r
+{\r
+  // Constructor\r
+  // Define input and output slots here\r
+  // Input slot #0 works with a TChain\r
+  DefineInput(0, TChain::Class());\r
+  // Output slot #0 writes into a TH1 container\r
+  DefineOutput(1, TList::Class());\r
+  DefineOutput(2, TList::Class());\r
+  DefineOutput(3, TList::Class());\r
+}\r
+\r
+//________________________________________________________________________\r
+AliAnalysisTaskTriggeredBF::~AliAnalysisTaskTriggeredBF() {\r
+\r
+  // Destructor\r
+\r
+}\r
+\r
+//________________________________________________________________________\r
+void AliAnalysisTaskTriggeredBF::UserCreateOutputObjects() {\r
+  // Create histograms\r
+  // Called once\r
+  if(!fBalance) {\r
+    fBalance = new AliBalanceTriggered();\r
+    fBalance->SetAnalysisLevel("AOD");\r
+  }\r
+  if(fRunShuffling) {\r
+    if(!fShuffledBalance) {\r
+      fShuffledBalance = new AliBalanceTriggered();\r
+      fShuffledBalance->SetAnalysisLevel("AOD");\r
+    }\r
+  }\r
+\r
+  //QA list\r
+  fList = new TList();\r
+  fList->SetName("listQA");\r
+  fList->SetOwner();\r
+\r
+  //Balance Function list\r
+  fListTriggeredBF = new TList();\r
+  fListTriggeredBF->SetName("listTriggeredBF");\r
+  fListTriggeredBF->SetOwner();\r
+\r
+  if(fRunShuffling) {\r
+    fListTriggeredBFS = new TList();\r
+    fListTriggeredBFS->SetName("listTriggeredBFShuffled");\r
+    fListTriggeredBFS->SetOwner();\r
+  }\r
+  \r
+  \r
+  //Event stats.\r
+  TString gCutName[4] = {"Total","Offline trigger",\r
+                         "Vertex","Analyzed"};\r
+  fHistEventStats = new TH1F("fHistEventStats",\r
+                             "Event statistics;;N_{events}",\r
+                             4,0.5,4.5);\r
+  for(Int_t i = 1; i <= 4; i++)\r
+    fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
+  fList->Add(fHistEventStats);\r
+  \r
+  TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r
+  fHistCentStats = new TH2F("fHistCentStats",\r
+                           "Centrality statistics;;Cent percentile",\r
+                           9,-0.5,8.5,220,-5,105);\r
+  for(Int_t i = 1; i <= 9; i++)\r
+    fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
+  fList->Add(fHistCentStats);\r
+  \r
+  fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);\r
+  fList->Add(fHistTriggerStats);\r
+  \r
+  fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);\r
+  fList->Add(fHistTrackStats);\r
+\r
+  fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);\r
+  fList->Add(fHistNumberOfAcceptedTracks);\r
+\r
+  // Vertex distributions\r
+  fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);\r
+  fList->Add(fHistVx);\r
+  fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);\r
+  fList->Add(fHistVy);\r
+  fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);\r
+  fList->Add(fHistVz);\r
+\r
+  // QA histograms\r
+  fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);\r
+  fList->Add(fHistClus);\r
+  fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);\r
+  fList->Add(fHistChi2);\r
+  fHistDCA  = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); \r
+  fList->Add(fHistDCA);\r
+  fHistPt   = new TH1F("fHistPt","p_{T} distribution",200,0,10);\r
+  fList->Add(fHistPt);\r
+  fHistEta  = new TH1F("fHistEta","#eta distribution",200,-2,2);\r
+  fList->Add(fHistEta);\r
+  fHistPhi  = new TH1F("fHistPhi","#phi distribution",200,-20,380);\r
+  fList->Add(fHistPhi);\r
+  fHistPhiBefore  = new TH1F("fHistPhiBefore","#phi distribution",200,0.,2*TMath::Pi());\r
+  fList->Add(fHistPhiBefore);\r
+  fHistPhiAfter  = new TH1F("fHistPhiAfter","#phi distribution",200,0.,2*TMath::Pi());\r
+  fList->Add(fHistPhiAfter);\r
+  fHistV0M  = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);\r
+  fList->Add(fHistV0M);\r
+  TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};\r
+  fHistRefTracks  = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);\r
+  for(Int_t i = 1; i <= 6; i++)\r
+    fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
+  fList->Add(fHistRefTracks);\r
+\r
+\r
+\r
+  // Balance function histograms\r
+  // Initialize histograms if not done yet\r
+  if(!fBalance->GetHistNp()){\r
+    AliWarning("Histograms not yet initialized! --> Will be done now");\r
+    AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
+    fBalance->InitHistograms();\r
+  }\r
+\r
+  if(fRunShuffling) {\r
+    if(!fShuffledBalance->GetHistNp()) {\r
+      AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
+      AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
+      fShuffledBalance->InitHistograms();\r
+    }\r
+  }\r
+\r
+  fListTriggeredBF->Add(fBalance->GetHistNp());\r
+  fListTriggeredBF->Add(fBalance->GetHistNn());\r
+  fListTriggeredBF->Add(fBalance->GetHistNpn());\r
+  fListTriggeredBF->Add(fBalance->GetHistNnn());\r
+  fListTriggeredBF->Add(fBalance->GetHistNpp());\r
+  fListTriggeredBF->Add(fBalance->GetHistNnp());\r
+  \r
+  if(fRunShuffling) {\r
+    fListTriggeredBFS->Add(fShuffledBalance->GetHistNp());\r
+    fListTriggeredBFS->Add(fShuffledBalance->GetHistNn());\r
+    fListTriggeredBFS->Add(fShuffledBalance->GetHistNpn());\r
+    fListTriggeredBFS->Add(fShuffledBalance->GetHistNnn());\r
+    fListTriggeredBFS->Add(fShuffledBalance->GetHistNpp());\r
+    fListTriggeredBFS->Add(fShuffledBalance->GetHistNnp());\r
+  }  \r
+\r
+\r
\r
+\r
+  // Post output data.\r
+  PostData(1, fList);\r
+  PostData(2, fListTriggeredBF);\r
+  if(fRunShuffling) PostData(3, fListTriggeredBFS);\r
+}\r
+\r
+//________________________________________________________________________\r
+void AliAnalysisTaskTriggeredBF::UserExec(Option_t *) {\r
+  // Main loop\r
+  // Called for each event\r
+\r
+  TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
+\r
+  Float_t fCentrality           = 0.;\r
+  \r
+  // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)\r
+  vector<Double_t> *chargeVector[9];          // original charge\r
+  for(Int_t i = 0; i < 9; i++){\r
+    chargeVector[i]        = new vector<Double_t>;\r
+  }\r
+  \r
+  Double_t v_charge;\r
+  Double_t v_y;\r
+  Double_t v_eta;\r
+  Double_t v_phi;\r
+  Double_t v_p[3];\r
+  Double_t v_pt;\r
+  Double_t v_E;\r
+\r
+  // -------------------------------------------------------------                  \r
+  // AOD analysis (vertex and track cuts also here!!!!)\r
+  if(gAnalysisLevel == "AOD") {\r
+    AliAODEvent* aodEventMain = dynamic_cast<AliAODEvent*>(InputEvent()); \r
+    if(!aodEventMain) {\r
+      AliError("aodEventMain not available");\r
+      return;\r
+    }\r
+    \r
+    \r
+    AliAODHeader *aodHeaderMain = aodEventMain->GetHeader();\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
+    \r
+    if(fUseOfflineTrigger)\r
+      isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
+    \r
+    if(isSelectedMain) {\r
+      fHistEventStats->Fill(2); //triggered events\r
+      \r
+      //Centrality stuff (centrality in AOD header)\r
+      if(fUseCentrality) {\r
+       fCentrality = aodHeaderMain->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
+       \r
+       // 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
+      \r
+      if(vertexMain) {\r
+       Double32_t fCovMain[6];\r
+       vertexMain->GetCovarianceMatrix(fCovMain);\r
+       \r
+       if(vertexMain->GetNContributors() > 0) {\r
+         if(fCovMain[5] != 0) {\r
+           fHistEventStats->Fill(3); //events with a proper vertex\r
+           if(TMath::Abs(vertexMain->GetX()) < fVxMax) {\r
+             if(TMath::Abs(vertexMain->GetY()) < fVyMax) {\r
+               if(TMath::Abs(vertexMain->GetZ()) < fVzMax) {\r
+                 fHistEventStats->Fill(4); //analyzed events\r
+                 fHistVx->Fill(vertexMain->GetX());\r
+                 fHistVy->Fill(vertexMain->GetY());\r
+                 fHistVz->Fill(vertexMain->GetZ());\r
+                 \r
+                 // Loop over tracks in main event\r
+                 for (Int_t iTracksMain = 0; iTracksMain < aodEventMain->GetNumberOfTracks(); iTracksMain++) {\r
+                   AliAODTrack* aodTrackMain = dynamic_cast<AliAODTrack *>(aodEventMain->GetTrack(iTracksMain));\r
+                   if (!aodTrackMain) {\r
+                     AliError(Form("Could not receive track %d", iTracksMain));\r
+                     continue;\r
+                   }\r
+                   \r
+                   // AOD track cuts\r
+                   \r
+                   // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
+                   // take only TPC only tracks \r
+                   fHistTrackStats->Fill(aodTrackMain->GetFilterMap());\r
+                   if(!aodTrackMain->TestFilterBit(nAODtrackCutBit)) continue;\r
+                   \r
+                   v_charge = aodTrackMain->Charge();\r
+                   v_y      = aodTrackMain->Y();\r
+                   v_eta    = aodTrackMain->Eta();\r
+                   v_phi    = aodTrackMain->Phi() * TMath::RadToDeg();\r
+                   v_E      = aodTrackMain->E();\r
+                   v_pt     = aodTrackMain->Pt();\r
+                   aodTrackMain->PxPyPz(v_p);\r
+                   \r
+                   Float_t DCAxy = aodTrackMain->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
+                   Float_t DCAz  = aodTrackMain->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
+                   \r
+                   \r
+                   // Kinematics cuts from ESD track cuts\r
+                   if( v_pt < fPtMin || v_pt > fPtMax)      continue;\r
+                   if( v_eta < fEtaMin || v_eta > fEtaMax)  continue;\r
+                   \r
+                   // Extra DCA cuts (for systematic studies [!= -1])\r
+                   if( fDCAxyCut != -1 && fDCAzCut != -1){\r
+                     if(TMath::Sqrt((DCAxy*DCAxy)/(fDCAxyCut*fDCAxyCut)+(DCAz*DCAz)/(fDCAzCut*fDCAzCut)) > 1 ){\r
+                       continue;  // 2D cut\r
+                     }\r
+                   }\r
+                   \r
+                   // Extra TPC cuts (for systematic studies [!= -1])\r
+                   if( fTPCchi2Cut != -1 && aodTrackMain->Chi2perNDF() > fTPCchi2Cut){\r
+                     continue;\r
+                   }\r
+                   if( fNClustersTPCCut != -1 && aodTrackMain->GetTPCNcls() < fNClustersTPCCut){\r
+                     continue;\r
+                   }\r
+                   \r
+                   // fill QA histograms\r
+                   fHistClus->Fill(aodTrackMain->GetITSNcls(),aodTrackMain->GetTPCNcls());\r
+                   fHistDCA->Fill(DCAz,DCAxy);\r
+                   fHistChi2->Fill(aodTrackMain->Chi2perNDF());\r
+                   fHistPt->Fill(v_pt);\r
+                   fHistEta->Fill(v_eta);\r
+                   fHistPhi->Fill(v_phi);\r
+                   \r
+                   // fill charge vector\r
+                   chargeVector[0]->push_back(v_charge);\r
+                   chargeVector[1]->push_back(v_y);\r
+                   chargeVector[2]->push_back(v_eta);\r
+                   chargeVector[3]->push_back(v_phi);\r
+                   chargeVector[4]->push_back(v_p[0]);\r
+                   chargeVector[5]->push_back(v_p[1]);\r
+                   chargeVector[6]->push_back(v_p[2]);\r
+                   chargeVector[7]->push_back(v_pt);\r
+                   chargeVector[8]->push_back(v_E);\r
+                   \r
+                 } //track loop\r
+                 \r
+                   // calculate balance function\r
+                 fBalance->FillBalance(fCentrality,chargeVector);\r
+                 // clean charge vector afterwards\r
+                 for(Int_t i = 0; i < 9; i++){                \r
+                   chargeVector[i]->clear();\r
+                 }\r
+\r
+               }//Vz cut\r
+             }//Vy cut\r
+           }//Vx cut\r
+         }//proper vertex resolution\r
+       }//proper number of contributors\r
+      }//vertex object valid\r
+    }//triggered event \r
+  }//AOD analysis\r
+  else{\r
+    AliError("Triggered Balance Function analysis only for AODs!");\r
+  }\r
+}     \r
+\r
+//________________________________________________________________________\r
+void  AliAnalysisTaskTriggeredBF::FinishTaskOutput(){\r
+\r
+  if (!fBalance) {\r
+    AliError("fBalance not available");\r
+    return;\r
+  }  \r
+  if(fRunShuffling) {\r
+    if (!fShuffledBalance) {\r
+      AliError("fShuffledBalance not available");\r
+      return;\r
+    }\r
+  }\r
+\r
+}\r
+\r
+//________________________________________________________________________\r
+void AliAnalysisTaskTriggeredBF::Terminate(Option_t *) {\r
+  // Called once at the end of the query\r
+\r
+  // not implemented ...\r
+\r
+}\r
+\r
+void AliAnalysisTaskTriggeredBF::UserExecMix(Option_t *)\r
+{\r
+\r
+  // not yet done for event mixing!\r
+  return;\r
+\r
+}\r
+\r
diff --git a/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskTriggeredBF.h b/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskTriggeredBF.h
new file mode 100755 (executable)
index 0000000..3da0a2b
--- /dev/null
@@ -0,0 +1,159 @@
+#ifndef ALIANALYSISTASKTRIGGEREDBF_CXX\r
+#define ALIANALYSISTASKTRIGGEREDBF_CXX\r
+\r
+// Analysis task for the TriggeredBF 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 AliBalanceTriggered;\r
+class AliESDtrackCuts;\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+#include "AliBalanceTriggered.h"\r
+\r
\r
+\r
+class AliAnalysisTaskTriggeredBF : public AliAnalysisTaskSE {\r
+ public:\r
+  AliAnalysisTaskTriggeredBF(const char *name = "AliAnalysisTaskTriggeredBF");\r
+  virtual ~AliAnalysisTaskTriggeredBF(); \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(AliBalanceTriggered *const analysis) {\r
+    fBalance         = analysis;\r
+    }\r
+  void SetShufflingObject(AliBalanceTriggered *const analysisShuffled) {\r
+    fRunShuffling = kTRUE;\r
+    fShuffledBalance = analysisShuffled;\r
+  }\r
\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
+  //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
+\r
+ private:\r
+  AliBalanceTriggered *fBalance; //TriggeredBF object\r
+  Bool_t fRunShuffling;//run shuffling or not\r
+  AliBalanceTriggered *fShuffledBalance; //TriggeredBF object (shuffled)\r
+  TList *fList; //fList object\r
+  TList *fListTriggeredBF; //fList object\r
+  TList *fListTriggeredBFS; //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
\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
\r
+  AliAnalysisTaskTriggeredBF(const AliAnalysisTaskTriggeredBF&); // not implemented\r
+  AliAnalysisTaskTriggeredBF& operator=(const AliAnalysisTaskTriggeredBF&); // not implemented\r
+  \r
+  ClassDef(AliAnalysisTaskTriggeredBF, 1); // example of analysis\r
+};\r
+\r
+#endif\r
index 5ff4481..dcd080a 100644 (file)
@@ -141,7 +141,7 @@ class AliBalanceEventMixing : public TObject {
 
   AliBalanceEventMixing & operator=(const AliBalanceEventMixing & ) {return *this;}
 
-  ClassDef(AliBalanceEventMixing, 3)
+  ClassDef(AliBalanceEventMixing, 1)
 };
 
 #endif
diff --git a/PWGCF/EBYE/BalanceFunctions/AliBalanceTriggered.cxx b/PWGCF/EBYE/BalanceFunctions/AliBalanceTriggered.cxx
new file mode 100644 (file)
index 0000000..7c296e2
--- /dev/null
@@ -0,0 +1,455 @@
+/**************************************************************************
+ * 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 <TH1D.h>
+
+
+#include <AliTHn.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 "AliBalanceTriggered.h"
+
+ClassImp(AliBalanceTriggered)
+
+//____________________________________________________________________//
+AliBalanceTriggered::AliBalanceTriggered() :
+  TObject(), 
+  bShuffle(kFALSE),
+  fAnalysisLevel("AOD"),
+  fHistP(0x0),
+  fHistN(0x0),
+  fHistPN(0x0),
+  fHistNP(0x0),
+  fHistPP(0x0),
+  fHistNN(0x0)
+{
+  // Default constructor
+}
+
+
+//____________________________________________________________________//
+AliBalanceTriggered::AliBalanceTriggered(const AliBalanceTriggered& balance):
+  TObject(balance), bShuffle(balance.bShuffle), 
+  fAnalysisLevel(balance.fAnalysisLevel),
+  fHistP(balance.fHistP),
+  fHistN(balance.fHistN),
+  fHistPN(balance.fHistPN),
+  fHistNP(balance.fHistNP),
+  fHistPP(balance.fHistPP),
+  fHistNN(balance.fHistNN)
+{
+  //copy constructor
+
+}
+
+//____________________________________________________________________//
+AliBalanceTriggered::~AliBalanceTriggered() {
+  // Destructor
+  delete fHistP;
+  delete fHistN;
+  delete fHistPN;
+  delete fHistNP;
+  delete fHistPP;
+  delete fHistNN;
+  
+
+}
+
+//____________________________________________________________________//
+void AliBalanceTriggered::InitHistograms() {
+
+  //Initialize the histograms
+
+  TString title    = "";      // histogram title
+  Int_t anaSteps   = 1;       // analysis steps
+
+  // single particle histograms
+  Int_t iBinSingle[nTrackVarsSingle];         // binning for track variables
+  Double_t* dBinsSingle[nTrackVarsSingle];    // bins for track variables  
+  TString axisTitleSingle[nTrackVarsSingle];  // axis titles for track variables
+  
+  // two particle histograms
+  Int_t iBinPair[nTrackVarsPair];         // binning for track variables
+  Double_t* dBinsPair[nTrackVarsPair];    // bins for track variables  
+  TString axisTitlePair[nTrackVarsPair];  // axis titles for track variables
+
+   
+  //-----------------------------------------------------------
+  // histogram settings (hard coded!)
+  //-----------------------------------------------------------
+
+  // eta
+  const Int_t kNEtaBins = 40;
+  Double_t etaBins[kNEtaBins+1];
+  for(Int_t i = 0; i < kNEtaBins+1; i++){
+    etaBins[i] = -1.0 + i * 0.05;
+  } 
+  iBinSingle[0]       = kNEtaBins;
+  dBinsSingle[0]      = etaBins;
+  axisTitleSingle[0]  = "#eta"; 
+  
+  // delta eta
+  const Int_t kNDeltaEtaBins = 40;
+  Double_t deltaEtaBins[kNDeltaEtaBins+1];
+  for(Int_t i = 0; i < kNDeltaEtaBins+1; i++){
+    deltaEtaBins[i] = -2.0 + i * 0.1;
+  } 
+  iBinPair[0]       = kNDeltaEtaBins;
+  dBinsPair[0]      = deltaEtaBins;
+  axisTitlePair[0]  = "#Delta #eta"; 
+
+  // phi
+  const Int_t kNPhiBins = 72;
+  Double_t phiBins[kNPhiBins+1];
+  for(Int_t i = 0; i < kNPhiBins+1; i++){
+    phiBins[i] = 0.0 + i * 5.;
+  } 
+  iBinSingle[1]       = kNPhiBins;
+  dBinsSingle[1]      = phiBins;
+  axisTitleSingle[1]  = "#phi (#circ)"; 
+  
+  // delta phi
+  const Int_t kNDeltaPhiBins = 72;
+  Double_t deltaPhiBins[kNDeltaPhiBins+1];
+  for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){
+    deltaPhiBins[i] = -180.0 + i * 5.;
+  } 
+  iBinPair[1]       = kNDeltaPhiBins;
+  dBinsPair[1]      = deltaPhiBins;
+  axisTitlePair[1]  = "#Delta #phi (#circ)"; 
+
+  // pT
+  const Int_t kNPtBins = 22;
+  Double_t pTBins[kNPtBins+1] = {0.15, 0.2, 0.3, 0.4, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 15.0, 20.0};
+  iBinSingle[2]       = kNPtBins;
+  dBinsSingle[2]      = pTBins;
+  axisTitleSingle[2]  = "p_{T,trig} (GeV/c)";
+  iBinPair[2]       = kNPtBins;
+  dBinsPair[2]      = pTBins;
+  axisTitlePair[2]  = "p_{T} (GeV/c)";
+  iBinPair[3]       = kNPtBins;
+  dBinsPair[3]      = pTBins;
+  axisTitlePair[3]  = "p_{T,trig} (GeV/c)"; 
+
+  // centrality
+  const Int_t kNCentBins = 9;
+  Double_t centBins[kNCentBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
+
+  iBinSingle[3]       = kNCentBins;
+  dBinsSingle[3]      = centBins;
+  axisTitleSingle[3]  = "centrality"; 
+
+  iBinPair[4]       = kNCentBins;
+  dBinsPair[4]      = centBins;
+  axisTitlePair[4]  = "centrality";
+
+  //-----------------------------------------------------------
+  
+
+
+  // User settings (depending on the analysis settings [only one at the moment])
+  if(fAnalysisLevel == "AOD"){
+    title = "hdEtaVsdPhi";
+  }
+
+  //-----------------------------------------------------------
+  //-----------------------------------------------------------
+  // Histogram creation
+
+  // histogram for negative particles
+  fHistN = new AliTHn(Form("fHistN"), Form("%s_N",title.Data()), anaSteps, nTrackVarsSingle, iBinSingle);
+  for (Int_t j=0; j<nTrackVarsSingle; j++)
+    {
+      fHistN->SetBinLimits(j, dBinsSingle[j]);
+      fHistN->SetVarTitle(j, axisTitleSingle[j]);
+    }
+
+  // histogram for positive particles
+  fHistP = new AliTHn(Form("fHistP"), Form("%s_P",title.Data()), anaSteps, nTrackVarsSingle, iBinSingle);
+  for (Int_t j=0; j<nTrackVarsSingle; j++)
+    {
+      fHistP->SetBinLimits(j, dBinsSingle[j]);
+      fHistP->SetVarTitle(j, axisTitleSingle[j]);
+    }
+
+  // histogram for +- pairs
+  fHistPN = new AliTHn(Form("fHistPN"), Form("%s_PN",title.Data()), anaSteps, nTrackVarsPair, iBinPair);
+  for (Int_t j=0; j<nTrackVarsPair; j++)
+    {
+      fHistPN->SetBinLimits(j, dBinsPair[j]);
+      fHistPN->SetVarTitle(j, axisTitlePair[j]);
+    }
+
+  // histogram for -+ pairs
+  fHistNP = new AliTHn(Form("fHistNP"), Form("%s_NP",title.Data()), anaSteps, nTrackVarsPair, iBinPair);
+  for (Int_t j=0; j<nTrackVarsPair; j++)
+    {
+      fHistNP->SetBinLimits(j, dBinsPair[j]);
+      fHistNP->SetVarTitle(j, axisTitlePair[j]);
+    }
+
+  // histogram for ++ pairs
+  fHistPP = new AliTHn(Form("fHistPP"), Form("%s_PP",title.Data()), anaSteps, nTrackVarsPair, iBinPair);
+  for (Int_t j=0; j<nTrackVarsPair; j++)
+    {
+      fHistPP->SetBinLimits(j, dBinsPair[j]);
+      fHistPP->SetVarTitle(j, axisTitlePair[j]);
+    }
+
+  // histogram for -- pairs
+  fHistNN = new AliTHn(Form("fHistNN"), Form("%s_NN",title.Data()), anaSteps, nTrackVarsPair, iBinPair);
+  for (Int_t j=0; j<nTrackVarsPair; j++)
+    {
+      fHistNN->SetBinLimits(j, dBinsPair[j]);
+      fHistNN->SetVarTitle(j, axisTitlePair[j]);
+    }
+
+  //-----------------------------------------------------------
+  //-----------------------------------------------------------
+
+
+}
+
+//____________________________________________________________________//
+void AliBalanceTriggered::FillBalance(Float_t fCentrality,vector<Double_t> **chargeVector) {
+  // Calculates the balance function
+
+  // Initialize histograms if not done yet
+  if(!fHistPN){
+    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();
+  Double_t trackVarsSingle[nTrackVarsSingle];
+  Double_t trackVarsPair[nTrackVarsPair];
+
+  // 1st particle loop
+  for(Int_t i = 0; i < gNtrack;i++){
+
+    Short_t  charge = (Short_t) chargeVector[0]->at(i);
+    trackVarsSingle[0]    =  chargeVector[2]->at(i);  //eta
+    trackVarsSingle[1]    =  chargeVector[3]->at(i);  //phi
+    trackVarsSingle[2]    =  chargeVector[7]->at(i);  //pt trigger
+    trackVarsSingle[3]    =  fCentrality;             //centrality (really as variable here????)
+
+    //fill single particle histograms
+    if(charge > 0)  fHistP->Fill(trackVarsSingle,0,1.); 
+    else            fHistN->Fill(trackVarsSingle,0,1.); 
+
+    // 2nd particle loop
+    for(Int_t j = 0; j < i; j++) {
+
+      // need check for single particle region!!!???
+      //
+      
+      Short_t charge2 = (Short_t) chargeVector[0]->at(j);
+      trackVarsPair[0]    =  chargeVector[2]->at(i) - chargeVector[2]->at(j) ;  //delta eta
+      trackVarsPair[1]    =  chargeVector[3]->at(i) - chargeVector[3]->at(j);  //delta phi
+      trackVarsPair[2]    =  chargeVector[7]->at(j);  //pt
+      trackVarsPair[3]    =  chargeVector[7]->at(i);  //pt trigger
+      trackVarsPair[4]    =  fCentrality;             //centrality (really as variable here????)
+
+    if( charge > 0 && charge2 < 0)  fHistPN->Fill(trackVarsPair,0,1.); 
+    else if( charge < 0 && charge2 > 0)  fHistNP->Fill(trackVarsPair,0,1.); 
+    else if( charge > 0 && charge2 > 0)  fHistPP->Fill(trackVarsPair,0,1.); 
+    else if( charge < 0 && charge2 < 0)  fHistNN->Fill(trackVarsPair,0,1.); 
+    else AliWarning("Wrong charge combination!");
+
+    }//end of 2nd particle loop
+  }//end of 1st particle loop
+}  
+
+
+TH1D* AliBalanceTriggered::GetBalanceFunctionHistogram1D(Int_t var, Double_t pTMin, Double_t pTMax, Double_t centrMin, Double_t centrMax){
+
+  // check which variable should be analyzed
+  // 0 = Delta eta
+  // 1 = Delta phi
+
+  if( var < 0 || var > 1){
+    AliError("Only Variable 0 (= Delta eta) or 1 (= Delta phi) allowed");
+    return NULL;
+  }
+
+
+  // Choose region to analyze 
+  // for Single Histograms (P,N):       2 = pT,trigger; 3 = centrality
+  // for Pair Histograms (PN,NP,NN,PP): 2 = pT; 3 = pT,trigger; 4 = centrality
+
+  // pT trigger
+  fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMin,pTMax); 
+  fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMin,pTMax); 
+  fHistPN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMin,pTMax); 
+  fHistNP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMin,pTMax); 
+  fHistPP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMin,pTMax); 
+  fHistNN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMin,pTMax); 
+
+  // pT
+  fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax); 
+  fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax); 
+  fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax); 
+  fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax); 
+
+  // centrality
+  fHistP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(centrMin,centrMax); 
+  fHistN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(centrMin,centrMax); 
+  fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax); 
+  fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax); 
+  fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax); 
+  fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax); 
+  
+
+  // Project into the wanted space (1st: analysis step, 2nd: axis)
+  TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,var);
+  TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,var);
+  TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,var);
+  TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,var);
+  TH1D* hTemp5 = (TH1D*)fHistP->Project(0,var);
+  TH1D* hTemp6 = (TH1D*)fHistN->Project(0,var);
+
+  TH1D* gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
+  gHistBalanceFunctionHistogram->Reset();
+  
+  // Calculate BF
+  if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)) {
+    hTemp1->Sumw2();
+    hTemp2->Sumw2();
+    hTemp3->Sumw2();
+    hTemp4->Sumw2();
+    hTemp1->Add(hTemp3,-1.);
+    hTemp1->Scale(1./hTemp5->GetEntries());
+    hTemp2->Add(hTemp4,-1.);
+    hTemp2->Scale(1./hTemp6->GetEntries());
+    gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
+    gHistBalanceFunctionHistogram->Scale(0.5/1.);
+  }
+
+  return gHistBalanceFunctionHistogram;
+}
+
+TH1D* AliBalanceTriggered::GetHistogram1D(Int_t histo, Int_t var, Double_t pTMin, Double_t pTMax, Double_t centrMin, Double_t centrMax){
+
+  // check which variable should be analyzed
+  //
+  // pair histograms:
+  // 0 = Delta eta 
+  // 1 = Delta phi
+  // 2 = pT, trigger
+  // 3 = centrality
+  //
+  // pair histograms:
+  // 0 = Delta eta 
+  // 1 = Delta phi
+  // 2 = pT
+  // 3 = pT, trigger
+  // 4 = centrality
+
+  if(histo < 0 || histo > 5){
+    AliError("Only 6 histograms available: 0(P), 1(N), 2(PN), 3(NP), 4(PP), 5(NN)");
+    return NULL;
+  }
+
+  if( histo > 1 && (var < 0 || var > 5)){
+    AliError("Only Variable 0 to 4 allowed for pair histograms (histo > 1)");
+    return NULL;
+  }
+  if( histo < 2 && (var < 0 || var > 4)){
+    AliError("Only Variable 0 to 3 allowed for single histograms (histo < 2)");
+    return NULL;
+  }
+
+  // get the histogram
+  AliTHn *gTHn = NULL;
+  switch(histo){
+  case 0:
+    gTHn = fHistP;
+    break;
+
+  case 1:
+    gTHn = fHistN;
+    break;
+
+  case 2:
+    gTHn = fHistPN;
+    break;
+
+  case 3:
+    gTHn = fHistNP;
+    break;
+
+  case 4:
+    gTHn = fHistPP;
+    break;
+
+  case 5:
+    gTHn = fHistNN;
+    break;
+
+  default:
+    break;
+
+  }
+
+  if(!gTHn){
+    AliError(Form("AliTHn number %d = NULL",histo));
+    return NULL;
+  }
+
+  // Choose region to analyze 
+  // for Single Histograms (P,N):       2 = pT,trigger; 3 = centrality
+  // for Pair Histograms (PN,NP,NN,PP): 2 = pT; 3 = pT,trigger; 4 = centrality
+
+  // pT trigger
+  gTHn->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMin,pTMax); 
+  // pT
+  if(histo > 1) gTHn->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax); 
+  // centrality
+  if(histo < 2) gTHn->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(centrMin,centrMax); 
+  else          gTHn->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax); 
+
+  // Project into the wanted space (1st: analysis step, 2nd: axis)
+  TH1D* gHisto = (TH1D*)gTHn->Project(0,var);
+
+  return gHisto;
+}
diff --git a/PWGCF/EBYE/BalanceFunctions/AliBalanceTriggered.h b/PWGCF/EBYE/BalanceFunctions/AliBalanceTriggered.h
new file mode 100644 (file)
index 0000000..7b58eab
--- /dev/null
@@ -0,0 +1,84 @@
+#ifndef ALIBALANCETRIGGERED_H
+#define ALIBALANCETRIGGERED_H
+/*  See cxx source for full Copyright notice */
+
+
+//-------------------------------------------------------------------------
+//                          Class AliBalanceTriggered
+//   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"
+
+class TGraphErrors;
+class TH1D;
+class AliTHn;
+
+const Int_t nTrackVarsSingle = 4;       // track variables in histogram (eta, phi, pTtrig, centrality)
+const Int_t nTrackVarsPair   = 5;       // track variables in histogram (dEta, dPhi, pT, pTtrig, centrality)
+
+class AliBalanceTriggered : public TObject {
+ public:
+
+  AliBalanceTriggered();
+  AliBalanceTriggered(const AliBalanceTriggered& balance);
+  ~AliBalanceTriggered();
+  
+  // analysis setters
+  void SetAnalysisLevel(const char* analysisLevel) {
+    fAnalysisLevel = analysisLevel;}
+  void SetShuffle(Bool_t shuffle) {bShuffle = shuffle;}
+
+  // analysis getters
+  const char* GetAnalysisLevel() {return fAnalysisLevel.Data();}
+  const Bool_t GetShuffle() {return bShuffle;}
+
+  // initialize histograms
+  void InitHistograms(void);
+
+  // histogram setters
+  void SetHistNp(AliTHn *gHist)  { fHistP  = gHist;}
+  void SetHistNn(AliTHn *gHist)  { fHistN  = gHist;}
+  void SetHistNpn(AliTHn *gHist) { fHistPN = gHist;}
+  void SetHistNnp(AliTHn *gHist) { fHistNP = gHist;}
+  void SetHistNpp(AliTHn *gHist) { fHistPP = gHist;}
+  void SetHistNnn(AliTHn *gHist) { fHistNN = gHist;}
+  // histogram getters
+  AliTHn *GetHistNp() { return fHistP;}
+  AliTHn *GetHistNn() { return fHistN;}
+  AliTHn *GetHistNpn() { return fHistPN;}
+  AliTHn *GetHistNnp() { return fHistNP;}
+  AliTHn *GetHistNpp() { return fHistPP;}
+  AliTHn *GetHistNnn() { return fHistNN;}
+
+  // Fill balance function histograms
+  void FillBalance(Float_t fCentrality,vector<Double_t> **chargeVector);
+  // Get the balance function histogram 
+  TH1D *GetBalanceFunctionHistogram1D(Int_t var, Double_t pTMin, Double_t pTMax, Double_t centrMin, Double_t centrMax);
+
+  // Get 1D histogram
+  TH1D* GetHistogram1D(Int_t histo, Int_t var, Double_t pTMin, Double_t pTMax, Double_t centrMin, Double_t centrMax);
+
+ private:
+  Bool_t bShuffle; //shuffled balance function object
+  TString fAnalysisLevel; // ESD, AOD or MC
+
+  AliTHn *fHistP;  //N+
+  AliTHn *fHistN;  //N-
+  AliTHn *fHistPN; //N+-
+  AliTHn *fHistNP; //N-+
+  AliTHn *fHistPP; //N++
+  AliTHn *fHistNN; //N--
+
+  AliBalanceTriggered & operator=(const AliBalanceTriggered & ) {return *this;}
+
+  ClassDef(AliBalanceTriggered, 1)
+};
+
+#endif
diff --git a/PWGCF/EBYE/macros/AddTaskBalanceTriggered.C b/PWGCF/EBYE/macros/AddTaskBalanceTriggered.C
new file mode 100644 (file)
index 0000000..69d9f5e
--- /dev/null
@@ -0,0 +1,153 @@
+// now in options\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
+AliAnalysisTaskTriggeredBF *AddTaskBalanceTriggered(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.15,\r
+                                                Double_t ptMax=20,\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
+  TString outputFileName(fileNameBase);\r
+  outputFileName.Append(".root");\r
+\r
+  //===========================================================================\r
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
+  if (!mgr) {\r
+    ::Error("AddTaskTriggeredBF", "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("AddTaskTriggeredBF", "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
+  // setup the balance function objects\r
+  AliBalanceTriggered *bf  = 0;  // Balance Function object\r
+  AliBalanceTriggered *bfs = 0;  // shuffled Balance function object\r
+  \r
+  if (analysisType=="AOD"){\r
+    \r
+    bf = new AliBalanceTriggered();\r
+    bf->SetAnalysisLevel(analysisType);\r
+    bf->SetShuffle(kFALSE);\r
+    bf->InitHistograms();\r
+    \r
+    if(gRunShuffling){\r
+      bfs = new AliBalanceTriggered();\r
+      bfs->SetAnalysisLevel(analysisType);\r
+      bfs->SetShuffle(kTRUE);\r
+      bfs->InitHistograms();\r
+    }\r
+  }\r
+  else{\r
+    ::Error("AddTaskTriggeredBF", "analysis type NOT supported.");\r
+    return NULL;\r
+  }\r
+\r
+  // Create the task, add it to manager and configure it.\r
+  //===========================================================================\r
+  AliAnalysisTaskTriggeredBF *taskTriggeredBF = new AliAnalysisTaskTriggeredBF("TaskTriggeredBF");\r
+  taskTriggeredBF->SetAnalysisObject(bf);\r
+  if(gRunShuffling) taskTriggeredBF->SetShufflingObject(bfs);\r
+\r
+  taskTriggeredBF->SetCentralityPercentileRange(centrMin,centrMax);\r
+  if(analysisType == "AOD") {\r
+    // pt and eta cut (pt_min, pt_max, eta_min, eta_max)\r
+    taskTriggeredBF->SetAODtrackCutBit(AODfilterBit);\r
+    taskTriggeredBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax);\r
+    \r
+    // set extra DCA cuts (-1 no extra cut)\r
+    taskTriggeredBF->SetExtraDCACutsAOD(DCAxy,DCAz);\r
+    \r
+    // set extra TPC chi2 / nr of clusters cut\r
+    taskTriggeredBF->SetExtraTPCCutsAOD(maxTPCchi2, minNClustersTPC);\r
+    \r
+  }\r
+\r
+  // offline trigger selection (AliVEvent.h)\r
+  // taskTriggeredBF->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) taskTriggeredBF->SelectCollisionCandidates(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral);\r
+  else                taskTriggeredBF->SelectCollisionCandidates(AliVEvent::kMB);\r
+\r
+  // centrality estimator (default = V0M)\r
+  taskTriggeredBF->SetCentralityEstimator(centralityEstimator);\r
+  \r
+  // vertex cut (x,y,z)\r
+  taskTriggeredBF->SetVertexDiamond(.3,.3,vertexZ);\r
+  \r
+  //bf->PrintAnalysisSettings();\r
+  mgr->AddTask(taskTriggeredBF);\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 *coutTriggeredBF = mgr->CreateContainer(Form("listTriggeredBF_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  if(gRunShuffling) AliAnalysisDataContainer *coutTriggeredBFS = mgr->CreateContainer(Form("listTriggeredBFShuffled_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+\r
+  mgr->ConnectInput(taskTriggeredBF, 0, mgr->GetCommonInputContainer());\r
+  mgr->ConnectOutput(taskTriggeredBF, 1, coutQA);\r
+  mgr->ConnectOutput(taskTriggeredBF, 2, coutTriggeredBF);\r
+  if(gRunShuffling) mgr->ConnectOutput(taskTriggeredBF, 3, coutTriggeredBFS);\r
+\r
+  return taskTriggeredBF;\r
+}\r