Code as used for final data
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 24 Jun 2011 08:06:25 +0000 (08:06 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 24 Jun 2011 08:06:25 +0000 (08:06 +0000)
PWG4/UserTasks/DiHadronCorrelations/AddTask_cloizide_Dhc.C [new file with mode: 0644]
PWG4/UserTasks/DiHadronCorrelations/AliDhcTask.cxx [new file with mode: 0644]
PWG4/UserTasks/DiHadronCorrelations/AliDhcTask.h [new file with mode: 0644]
PWG4/UserTasks/DiHadronCorrelations/AliDhcTask_opt.cxx [new file with mode: 0644]
PWG4/UserTasks/DiHadronCorrelations/AliDhcTask_opt.h [new file with mode: 0644]
PWG4/UserTasks/DiHadronCorrelations/CheckOutput.C [new file with mode: 0644]
PWG4/UserTasks/DiHadronCorrelations/KiddiePoolClasses.cxx [new file with mode: 0644]
PWG4/UserTasks/DiHadronCorrelations/KiddiePoolClasses.h [new file with mode: 0644]
PWG4/UserTasks/DiHadronCorrelations/RunDhcTask.C [new file with mode: 0644]
PWG4/UserTasks/DiHadronCorrelations/rootlogon.C [new file with mode: 0644]

diff --git a/PWG4/UserTasks/DiHadronCorrelations/AddTask_cloizide_Dhc.C b/PWG4/UserTasks/DiHadronCorrelations/AddTask_cloizide_Dhc.C
new file mode 100644 (file)
index 0000000..ade131e
--- /dev/null
@@ -0,0 +1,33 @@
+AliAnalysisTask *AddTask_cloizide_Dhc() 
+{
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    Error("AddTask_cloizide_Dhc", "No analysis manager found.");
+    return 0;
+  }
+
+  TString list=gSystem->Getenv("LIST");
+  Bool_t isPass1 = list.Contains("pass1");
+  Bool_t isPass2 = list.Contains("pass2");
+
+  if (isPass1)
+    return 0;
+
+  gROOT->ProcessLine(".L KiddiePoolClasses.cxx+");
+  gROOT->ProcessLine(".L AliDhcTask_opt.cxx+");
+
+  AliDhcTask *dhcTask = new AliDhcTask("Task_cloizide_Dhc");
+  dhcTask->SelectCollisionCandidates(AliVEvent::kMB);
+  dhcTask->SetVerbosity(10);
+  mgr->AddTask(dhcTask);
+
+  AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+  mgr->ConnectInput(dhcTask,0,cinput);
+  AliAnalysisDataContainer *co1 = mgr->CreateContainer("Cont_cloizide_DhcAna",
+                                                       TList::Class(), 
+                                                       AliAnalysisManager::kOutputContainer,
+                                                       "cloizide_DhcAna.root");
+  mgr->ConnectOutput(dhcTask,1,co1);
+
+  return 0;
+}
diff --git a/PWG4/UserTasks/DiHadronCorrelations/AliDhcTask.cxx b/PWG4/UserTasks/DiHadronCorrelations/AliDhcTask.cxx
new file mode 100644 (file)
index 0000000..8a9769b
--- /dev/null
@@ -0,0 +1,538 @@
+// Dihadron correlations task - simple task to read ESD or AOD input,
+// calculate same- and mixed-event correlations, and fill THnSparse
+// output. -A. Adare, Apr 2011
+
+#include "TCanvas.h"
+#include "TChain.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "THnSparse.h"
+#include "TROOT.h"
+#include "TTree.h"
+
+#include "AliAODEvent.h"
+#include "AliAODInputHandler.h"
+#include "AliAnalysisManager.h"
+#include "AliAnalysisTask.h"
+#include "AliCentrality.h"
+#include "AliDhcTask.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliESDtrackCuts.h"
+#include "KiddiePoolClasses.h"
+#include "AliVParticle.h"
+
+ClassImp(AliDhcTask)
+
+//________________________________________________________________________
+AliDhcTask::AliDhcTask(const char *name) 
+: AliAnalysisTaskSE(name), fVerbosity(0), fESD(0), fAOD(0), fOutputList(0), 
+  fHistPt(0), fHEvt(0), fHTrk(0), fHS(0), fHM(0), fPoolMgr(0), 
+  fCentrality(99), fZVertex(99), fZVtxMax(10), fPtMin(0), fPtMax(100), 
+  fInputHandler(0), fEsdTrackCutsTPCOnly(0)
+{
+  // Constructor
+
+  // Define input and output slots here
+  // Input slot #0 works with a TChain
+  DefineInput(0, TChain::Class());
+  // Output slot #0 id reserved by the base class for AOD
+  // Output slot #1 writes into a TH1 container
+  DefineOutput(1, TList::Class());
+
+  // Initialize cut variables
+  fZVtxMax = 10;    // cm
+  fPtMin   = 0.50;  // GeV/c
+  fPtMax   = 15.;
+
+  /*
+  // Not used for anything now...
+  fInputHandler = (AliInputEventHandler*)
+    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
+  if (!fInputHandler) {
+    Error("AliDhcTask()", "Did not get input handler");
+  }
+  */
+
+  fEsdTrackCutsTPCOnly = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
+  fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,Tracks "
+               "AOD:header,tracks,vertices,";
+}
+
+//________________________________________________________________________
+void AliDhcTask::UserCreateOutputObjects()
+{
+  // Create histograms
+  // Called once (per slave on PROOF!)
+
+  BookHistos();
+  InitEventMixer(); 
+
+  fOutputList = new TList();
+  fOutputList->SetOwner(1);
+
+  fOutputList->Add(fHistPt);
+  fOutputList->Add(fHS);
+  fOutputList->Add(fHM);
+  fOutputList->Add(fHEvt);
+  fOutputList->Add(fHTrk);
+
+  PostData(1, fOutputList);
+}
+
+//________________________________________________________________________
+void AliDhcTask::BookHistos()
+{
+  // Setup for THnSparse correlation histos.
+  // Important! Order bins[] according to ePairHistAxes.
+
+  const Int_t ndims = 6;
+  Int_t nDeta=22, nPtAssc=11, nPtTrig=11, nCent=8, nDphi=36, nZvtx=8;
+  Int_t bins[ndims] = {nDeta, nPtAssc, nPtTrig, nCent, nDphi, nZvtx };
+  Double_t xmin[ndims] = { -2.2, 0.5, 0.5, 0,  -0.5*TMath::Pi(), -10 };
+  Double_t xmax[ndims] = { +2.2, 15., 15., 90, +1.5*TMath::Pi(), +10 };
+  fHS = new THnSparseF("fHS", "Same evt Correlations", 
+                      ndims, bins, xmin, xmax);
+
+  // Override the uniform binning for some axes. Keep xmin and xmax
+  // limits the same, but allow variable bin widths inside.
+  Double_t ptt[]  = {0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 6.0, 8.0, 10, 15};
+  Double_t pta[]  = {0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 6.0, 8.0, 10, 15};
+  Double_t cent[] = {0, 2, 10, 20, 30, 40, 50, 60, 90};
+  Double_t zvtx[] = {-10, -6, -4, -2, 0, 2, 4, 6, 10};
+  fHS->SetBinEdges(kPtTrig, ptt);
+  fHS->SetBinEdges(kPtAssc, pta);
+  fHS->SetBinEdges(kCent, cent);
+  fHS->SetBinEdges(kZvtx, zvtx); // Match InitEventMixer - no point in going finer
+
+  fHM = (THnSparse*) fHS->Clone("fHM");
+  fHM->SetTitle("Mixed evt Correlations");
+
+  // Event histo
+  fHEvt = new TH2F("fHEvt", "Event-level variables", 30, -15, 15, 101, 0, 101);
+  // Track histo
+  fHTrk = new TH2F("fHTrk", "Track-level variables", 
+                  100, 0, TMath::TwoPi(), 100, -2, +2);
+
+  // Left over from the tutorial :)
+  fHistPt = new TH1F("fHistPt", "P_{T} distribution", 200, 0., 20.);
+  fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
+  fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
+  fHistPt->SetMarkerStyle(kFullCircle);
+  
+  return;
+}
+
+//________________________________________________________________________
+void AliDhcTask::InitEventMixer()
+{
+  // The effective pool size in events is set by trackDepth, so more
+  // low-mult events are required to maintain the threshold than
+  // high-mult events. Centrality pools are indep. of data histogram
+  // binning, no need to match.
+
+  Int_t trackDepth = 1000;   // # tracks to fill pool
+  Int_t poolsize   = 200;    // Maximum number of events
+
+  // Centrality pools
+  Int_t nCentBins  = 9;
+  Double_t centBins[] = {0,1,2,5,10,20,30,40,60,90.1};
+  //Int_t nCentBins  = 1;
+  //Double_t centBins[] = {-1,100.1};
+  // Z-vertex pools
+  Int_t nZvtxBins  = 8;
+  Double_t zvtxbin[] = {-10, -6, -4, -2, 0, 2, 4, 6, 10};
+  
+  fPoolMgr = new KiddiePoolManager(poolsize, trackDepth, nCentBins, 
+                                  centBins, nZvtxBins, zvtxbin);
+  return;
+}
+
+//________________________________________________________________________
+void AliDhcTask::UserExec(Option_t *) 
+{
+  // Main loop, called for each event.
+
+  static int iEvent = -1; ++iEvent;
+
+  if (fVerbosity>2) {
+    if (iEvent % 10 == 0) 
+      cout << iEvent << endl;
+  }
+
+  Int_t dType = -1;       // Will be set to kESD or kAOD.
+  MiniEvent* sTracks = 0; // Vector of selected MiniTracks.
+  Double_t centCL1 = -1;
+
+  LoadBranches();
+
+  // Get event pointers, check for signs of life
+  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
+  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
+  if (fESD)
+    dType = kESD;
+  else if (fAOD)
+    dType = kAOD;
+  else {
+    AliError("Neither fESD nor fAOD available");
+    return;
+  }
+
+  // Centrality, vertex, other event variables...
+  if (dType == kESD) {
+    if (!VertexOk(fESD)) {
+      if (fVerbosity > 1)
+       AliInfo(Form("Event REJECTED. z = %.1f", fZVertex));
+      return;
+    }
+    const AliESDVertex* vertex = fESD->GetPrimaryVertex();
+    fZVertex = vertex->GetZ();
+    if(fESD->GetCentrality()) {
+      fCentrality = 
+       fESD->GetCentrality()->GetCentralityPercentile("V0M");
+      centCL1 =
+       fESD->GetCentrality()->GetCentralityPercentile("CL1");
+    }
+  }
+  if (dType == kAOD) {
+    const AliAODVertex* vertex = fAOD->GetPrimaryVertex();
+    fZVertex = vertex->GetZ();
+    if (!VertexOk(fAOD)) {
+      if (fVerbosity > 1)
+       Info("Exec()", "Event REJECTED. z = %.1f", fZVertex);
+      return;
+    }
+    const AliCentrality *aodCent = fAOD->GetHeader()->GetCentralityP();
+    if (aodCent) {
+      fCentrality = aodCent->GetCentralityPercentile("V0M");
+      centCL1     = aodCent->GetCentralityPercentile("CL1");
+    }
+  }
+  
+  // Fill Event histogram
+  fHEvt->Fill(fZVertex, fCentrality);
+
+  if (fCentrality > 90. || fCentrality < 0) {
+    AliInfo(Form("Event REJECTED. fCentrality = %.1f", fCentrality));
+    return;
+  }
+
+  // Get array of selected tracks
+  if (dType == kESD) {
+    sTracks = GetESDTrax();
+  }
+  if (dType == kAOD) {
+    sTracks = GetAODTrax();
+  }
+
+  // Get pool containing tracks from other events like this one
+  KiddiePool* pool = fPoolMgr->GetEventPool(fCentrality, fZVertex);
+  if (!pool) {
+    AliFatal(Form("No pool found. Centrality %f, ZVertex %f", fCentrality, fZVertex));
+    return;
+  }
+
+  if (!pool->IsReady()) {
+    pool->UpdatePool(sTracks);
+    return;
+  }
+
+  if (pool->IsFirstReady()) {
+    // recover events missed before the pool is ready
+    Int_t nEvs = pool->GetCurrentNEvents();
+    if (nEvs>1) {
+      for (Int_t i=0; i<nEvs; ++i) {
+       MiniEvent* evI = pool->GetEvent(i);
+       for (Int_t j=0; j<nEvs; ++j) {
+         MiniEvent* evJ = pool->GetEvent(j);
+         if (i==j) {
+           Correlate(*evI, *evJ, kSameEvt);
+         } else {
+           Correlate(*evI, *evJ, kDiffEvt, 1.0);
+         }
+       }
+      }
+    }
+  } else { /* standard case: same event, then mix*/
+    Correlate(*sTracks, *sTracks, kSameEvt);  
+    Int_t nMix = pool->GetCurrentNEvents();
+    for (Int_t jMix=0; jMix<nMix; ++jMix) {
+      MiniEvent* bgTracks = pool->GetEvent(jMix);
+      Correlate(*sTracks, *bgTracks, kDiffEvt, 1.0);
+    }
+  }
+
+  if (fVerbosity>4) {
+    cout << "Output of SAME  THnSparse: " << fHS->GetSparseFractionBins() << " " << fHS->GetSparseFractionMem() << endl; 
+    cout << "Output of MIXED THnSparse: " << fHM->GetSparseFractionBins() << " " << fHM->GetSparseFractionMem() << endl; 
+  }
+
+  pool->UpdatePool(sTracks);
+  PostData(1, fOutputList);
+  return;
+}
+
+//________________________________________________________________________
+MiniEvent* AliDhcTask::GetESDTrax() const
+{
+  // Loop twice: 1. Count sel. tracks. 2. Fill vector.
+
+  Int_t nTrax = fESD->GetNumberOfTracks();
+  Int_t nSelTrax = 0;
+
+  if (fVerbosity > 2)
+    AliInfo(Form("%d tracks in event",nTrax));
+
+  // Loop 1.
+  for (Int_t i = 0; i < nTrax; ++i) {
+    AliESDtrack* esdtrack = fESD->GetTrack(i);
+    if (!esdtrack) {
+      AliError(Form("Couldn't get ESD track %d\n", i));
+      continue;
+    }
+    Bool_t trkOK = fEsdTrackCutsTPCOnly->AcceptTrack(esdtrack);
+    if (!trkOK)
+      continue;
+    Double_t pt = esdtrack->Pt();
+    Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
+    if (!ptOK)
+      continue;
+    Double_t eta = esdtrack->Eta();
+    if (TMath::Abs(eta) > 1.0)
+      continue;
+    nSelTrax++;
+  }
+
+  MiniEvent* miniEvt = new MiniEvent(0);
+  miniEvt->reserve(nSelTrax);
+
+  // Loop 2.  
+  for (Int_t i = 0; i < nTrax; ++i) {
+    AliESDtrack* esdtrack = fESD->GetTrack(i);
+    if (!esdtrack) {
+      AliError(Form("Couldn't get ESD track %d\n", i));
+      continue;
+    }
+    Bool_t trkOK = fEsdTrackCutsTPCOnly->AcceptTrack(esdtrack);
+    if (!trkOK)
+      continue;
+    Double_t pt = esdtrack->Pt();
+    Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
+    if (!ptOK)
+      continue;
+    Double_t eta  = esdtrack->Eta();
+    if (TMath::Abs(eta) > 1.0)
+      continue;
+
+    Double_t phi  = esdtrack->Phi();
+    Int_t    sign = esdtrack->Charge() > 0 ? 1 : -1;
+    miniEvt->push_back(MiniTrack(pt, eta, phi, sign));
+  }
+  return miniEvt;
+}
+
+//________________________________________________________________________
+MiniEvent* AliDhcTask::GetAODTrax() const
+{
+  // Loop twice: 1. Count sel. tracks. 2. Fill vector.
+
+  Int_t nTrax = fAOD->GetNumberOfTracks();
+  Int_t nSelTrax = 0;
+
+  if (fVerbosity > 2)
+    AliInfo(Form("%d tracks in event",nTrax));
+
+  // Loop 1.
+  for (Int_t i = 0; i < nTrax; ++i) {
+    AliAODTrack* aodtrack = fAOD->GetTrack(i);
+    if (!aodtrack) {
+      AliError(Form("Couldn't get AOD track %d\n", i));
+      continue;
+    }
+    // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
+    UInt_t tpcOnly = 1 << 7;
+    Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly);
+    if (!trkOK)
+      continue;
+    Double_t pt = aodtrack->Pt();
+    Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
+    if (!ptOK)
+      continue;
+    Double_t eta = aodtrack->Eta();
+    if (TMath::Abs(eta) > 1.0)
+      continue;
+    nSelTrax++;
+  }
+
+  MiniEvent* miniEvt = new MiniEvent(0);
+  miniEvt->reserve(nSelTrax);
+
+  // Loop 2.  
+  for (Int_t i = 0; i < nTrax; ++i) {
+    AliAODTrack* aodtrack = fAOD->GetTrack(i);
+    if (!aodtrack) {
+      AliError(Form("Couldn't get AOD track %d\n", i));
+      continue;
+    }
+    
+    // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
+    UInt_t tpcOnly = 1 << 7;
+    Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly);
+    if (!trkOK)
+      continue;
+    Double_t pt = aodtrack->Pt();
+    Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
+    if (!ptOK)
+      continue;
+    Double_t eta  = aodtrack->Eta();
+    if (TMath::Abs(eta) > 1.0)
+      continue;
+
+    Double_t phi  = aodtrack->Phi();
+    Int_t    sign = aodtrack->Charge() > 0 ? 1 : -1;
+    miniEvt->push_back(MiniTrack(pt, eta, phi, sign));
+  }
+  return miniEvt;
+}
+
+//________________________________________________________________________
+Double_t AliDhcTask::DeltaPhi(Double_t phia, Double_t phib,
+                             Double_t rangeMin, Double_t rangeMax) const
+{
+  Double_t dphi = -999;
+  Double_t pi = TMath::Pi();
+  
+  if (phia < 0)         phia += 2*pi;
+  else if (phia > 2*pi) phia -= 2*pi;
+  if (phib < 0)         phib += 2*pi;
+  else if (phib > 2*pi) phib -= 2*pi;
+  dphi = phib - phia;
+  if (dphi < rangeMin)      dphi += 2*pi;
+  else if (dphi > rangeMax) dphi -= 2*pi;
+  
+  return dphi;
+}
+
+//________________________________________________________________________
+Int_t AliDhcTask::Correlate(const MiniEvent &evt1, const MiniEvent &evt2, 
+                           Int_t pairing, Double_t weight)
+{
+  // Triggered angular correlations. If pairing is kSameEvt, particles
+  // within evt1 are correlated. If kDiffEvt, correlate triggers from
+  // evt1 with partners from evt2.
+  
+  Int_t iMax = evt1.size();
+  Int_t jMax = evt2.size();
+
+  THnSparse *hist = fHM;
+  if (pairing == kSameEvt)
+    hist = fHS;
+
+  for (Int_t i=0; i<iMax; ++i) {
+
+    // Trigger particles
+    const MiniTrack &a(evt1.at(i));
+
+    if (pairing == kSameEvt) {
+      fHistPt->Fill(a.Pt());
+      fHTrk->Fill(a.Phi(),a.Eta());
+    }
+
+    for (int j=0; j<jMax; ++j) {
+      // Associated particles
+      if (pairing == kSameEvt && i==j)
+       continue;
+
+      const MiniTrack &b(evt2.at(j));
+      
+      Float_t pta  = a.Pt();
+      Float_t ptb  = b.Pt();
+      if (pta < ptb) 
+       continue;
+
+      Float_t dphi = DeltaPhi(a.Phi(), b.Phi());
+      Float_t deta = a.Eta() - b.Eta();
+
+      Double_t x[6]; // Match ndims in fHS
+      x[kDeta]       = deta;
+      x[kPtAssc]     = ptb;
+      x[kPtTrig]     = pta;
+      x[kCent]       = fCentrality;
+      x[kDphi]       = dphi;
+      x[kZvtx]       = fZVertex;
+      //x[kChargeComb] = a.Sign() * b.Sign(); 
+      
+      hist->Fill(x, weight);
+    }
+  }
+
+  return 0;
+}
+
+//________________________________________________________________________
+void AliDhcTask::Terminate(Option_t *) 
+{
+  // Draw result to the screen
+  // Called once at the end of the query
+
+  if (gROOT->IsBatch())
+    return;
+
+  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
+  if (!fOutputList) {
+    AliError("Output list not available");
+    return;
+  }
+  
+  fHistPt = dynamic_cast<TH1F*> (fOutputList->At(0));
+  if (!fHistPt) {
+    AliError("ERROR: fHistPt not available\n");
+    return;
+  }
+   
+  TCanvas *c1 = new TCanvas("AliDhcTask","Pt",10,10,510,510);
+  c1->cd(1)->SetLogy();
+  fHistPt->DrawCopy("E");
+}
+
+//________________________________________________________________________
+Bool_t AliDhcTask::VertexOk(TObject* obj) const
+{
+  // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
+  
+  Int_t nContributors  = 0;
+  Double_t zVertex     = 999;
+  TString name("");
+  
+  if (obj->InheritsFrom("AliESDEvent")) {
+    AliESDEvent* esdevt = (AliESDEvent*) obj;
+    const AliESDVertex* vtx = esdevt->GetPrimaryVertex();
+    if (!vtx)
+      return 0;
+    nContributors = vtx->GetNContributors();
+    zVertex       = vtx->GetZ();
+    name          = vtx->GetName();
+  }
+  else if (obj->InheritsFrom("AliAODEvent")) { 
+    AliAODEvent* aodevt = (AliAODEvent*) obj;
+    if (aodevt->GetNumberOfVertices() < 1)
+      return 0;
+    const AliAODVertex* vtx = aodevt->GetPrimaryVertex();
+    nContributors = vtx->GetNContributors();
+    zVertex       = vtx->GetZ();
+    name          = vtx->GetName();
+  }
+  
+  // Reject if TPC-only vertex
+  if (name.CompareTo("TPCVertex"))
+    return kFALSE;
+  
+  // Check # contributors and range...
+  if( nContributors < 1 || TMath::Abs(zVertex) > fZVtxMax ) {
+    return kFALSE;
+  }
+  
+  return kTRUE;
+}
diff --git a/PWG4/UserTasks/DiHadronCorrelations/AliDhcTask.h b/PWG4/UserTasks/DiHadronCorrelations/AliDhcTask.h
new file mode 100644 (file)
index 0000000..4b654ce
--- /dev/null
@@ -0,0 +1,81 @@
+// Dihadron correlations task - simple task to read ESD or AOD input,
+// calculate same- and mixed-event correlations, and fill THnSparse
+// output. -A. Adare, Apr 2011
+
+
+#ifndef AliDhcTask_cxx
+#define AliDhcTask_cxx
+
+class TH1;
+class TH2;
+class THnSparse;
+class TObject;
+class TObjArray;
+class AliESDEvent;
+class AliAODEvent;
+class AliESDtrackCuts;
+class KiddiePoolManager;
+
+#include "AliAnalysisTaskSE.h"
+#include "KiddiePoolClasses.h"
+
+class AliDhcTask : public AliAnalysisTaskSE {
+ public:
+  AliDhcTask() : 
+    AliAnalysisTaskSE(), fVerbosity(0), fESD(0), fAOD(0), fOutputList(0), 
+    fHistPt(0), fHEvt(0), fHTrk(0), fHS(0), fHM(0), fPoolMgr(0), fCentrality(99), 
+    fZVertex(99), fZVtxMax(10), fPtMin(0), fPtMax(100), fInputHandler(0), 
+    fEsdTrackCutsTPCOnly(0) {}
+  AliDhcTask(const char *name);
+  virtual ~AliDhcTask() {}
+  
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+  virtual void   Terminate(Option_t *);
+  void SetVerbosity(const Int_t v) { fVerbosity = v; }
+
+ protected:
+  enum ePairHistAxes  {kDeta, kPtAssc, kPtTrig, kCent, kDphi,
+                      kZvtx, kChargeComb};
+  enum eEventHistAxes {kZvtxEvt, kCentV0M, kCentCL1};
+  enum eTrackHistAxes {kPhiTrk, kEtaTrk};
+  enum ePairingScheme {kSameEvt, kDiffEvt};
+  enum eDataType {kESD, kAOD};
+
+  void BookHistos();
+  void InitEventMixer();
+  MiniEvent* GetESDTrax() const;
+  MiniEvent* GetAODTrax() const;
+  Bool_t VertexOk(TObject* obj) const;
+  Double_t DeltaPhi(Double_t phia, Double_t phib, 
+                   Double_t rangeMin = -TMath::Pi()/2, 
+                   Double_t rangeMax = 3*TMath::Pi()/2) const;
+  Int_t Correlate(const MiniEvent &arr1, const MiniEvent &arr2, 
+                 Int_t pairing = kSameEvt, Double_t weight = 1.);
+
+ private:
+  Int_t        fVerbosity;       // 0 = silence
+  AliESDEvent *fESD;             //! ESD object
+  AliAODEvent *fAOD;             //! AOD object
+  TList       *fOutputList;      //! Output list
+  TH1F        *fHistPt;          //! Pt spectrum
+  TH2         *fHEvt;            //! Cent, vtx, etc.
+  TH2         *fHTrk;            //! Phi, Eta, etc.
+  THnSparse   *fHS;              //! Same-evt correlations
+  THnSparse   *fHM;              //! Diff-evt correlations
+  KiddiePoolManager* fPoolMgr;   //! Event mixer
+  Double_t    fCentrality;       //! V0M for now
+  Double_t    fZVertex;          //! Of current event
+  Double_t    fZVtxMax;          //! Max |z| cut (cm)
+  Double_t    fPtMin;            //! Min pt cut
+  Double_t    fPtMax;            //! Max pt cut
+  AliInputEventHandler*  fInputHandler;
+  AliESDtrackCuts* fEsdTrackCutsTPCOnly;
+
+  AliDhcTask(const AliDhcTask&);            // not implemented
+  AliDhcTask &operator=(const AliDhcTask&); // not implemented
+
+  ClassDef(AliDhcTask, 1);
+};
+
+#endif
diff --git a/PWG4/UserTasks/DiHadronCorrelations/AliDhcTask_opt.cxx b/PWG4/UserTasks/DiHadronCorrelations/AliDhcTask_opt.cxx
new file mode 100644 (file)
index 0000000..62d6bdd
--- /dev/null
@@ -0,0 +1,636 @@
+// Dihadron correlations task - simple task to read ESD or AOD input,
+// calculate same- and mixed-event correlations, and fill THnSparse
+// output. -A. Adare, Apr 2011
+
+#include "TCanvas.h"
+#include "TChain.h"
+#include "TFormula.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TProfile2D.h"
+#include "THnSparse.h"
+#include "TROOT.h"
+#include "TTree.h"
+#include "AliAODEvent.h"
+#include "AliAODInputHandler.h"
+#include "AliAnalysisManager.h"
+#include "AliAnalysisTask.h"
+#include "AliCentrality.h"
+#include "AliDhcTask_opt.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliESDtrackCuts.h"
+#include "KiddiePoolClasses.h"
+#include "AliVParticle.h"
+
+ClassImp(AliDhcTask)
+
+//________________________________________________________________________
+AliDhcTask::AliDhcTask(const char *name) 
+: AliAnalysisTaskSE(name), fVerbosity(0), fEtaMax(1), fZVtxMax(10), fPtMin(0.25), fPtMax(15), 
+  fESD(0), fAOD(0), fOutputList(0), fHistPt(0), fHEvt(0), fHTrk(0), fHSs(0), fHMs(0), 
+  fIndex(0), fMeanPtTrg(0), fMeanPtAss(0), fMean2PtTrg(0), fMean2PtAss(0),
+  fCentrality(99), fZVertex(99), fEsdTrackCutsTPCOnly(0), fPoolMgr(0)
+{
+  // Constructor
+
+  // Define input and output slots here
+  // Input slot #0 works with a TChain
+  DefineInput(0, TChain::Class());
+  // Output slot #0 id reserved by the base class for AOD
+  // Output slot #1 writes into a TH1 container
+  DefineOutput(1, TList::Class());
+
+  fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,Tracks "
+               "AOD:header,tracks,vertices,";
+}
+
+//________________________________________________________________________
+void AliDhcTask::UserCreateOutputObjects()
+{
+  // Create histograms
+  // Called once (per slave on PROOF!)
+
+  fOutputList = new TList();
+  fOutputList->SetOwner(1);
+
+  fEsdTrackCutsTPCOnly = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
+  //fEsdTrackCutsTPCOnly->SetMinNClustersTPC(70);
+  fEsdTrackCutsTPCOnly->SetMinNCrossedRowsTPC(70);
+  fEsdTrackCutsTPCOnly->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
+
+  BookHistos();
+  InitEventMixer(); 
+  PostData(1, fOutputList);
+}
+
+//________________________________________________________________________
+void AliDhcTask::BookHistos()
+{
+  // Book histograms.
+
+  Int_t nDeta=20, nPtAssc=12, nPtTrig=12, nCent=12, nDphi=36, nZvtx=8;
+
+  Double_t ptt[]  = {0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 8.0, 15};
+  Double_t pta[]  = {0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 8.0, 15};
+  Double_t cent[] = {0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 90};
+  Double_t zvtx[] = {-10, -6, -4, -2, 0, 2, 4, 6, 10};
+
+  // Event histo
+  fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", 30, -15, 15, 101, 0, 101);
+  fOutputList->Add(fHEvt);
+  // Track histo
+  fHTrk = new TH2F("fHTrk", "Track-level variables", 
+                  100, 0, TMath::TwoPi(), 100, -fEtaMax, +fEtaMax);
+  fOutputList->Add(fHTrk);
+  
+  // Left over from the tutorial :)
+  fHistPt = new TH1F("fHistPt", "P_{T} distribution", 200, 0., fPtMax);
+  fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
+  fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
+  fHistPt->SetMarkerStyle(kFullCircle);
+  fOutputList->Add(fHistPt);
+
+  fHPtAss = new TH1F("fHPtAss","PtAssoc;P_{T} (GeV/c) [GeV/c]",nPtAssc,pta);
+  fOutputList->Add(fHPtAss);
+  fHPtTrg = new TH1F("fHPtTrg","PtTrig;P_{T} (GeV/c) [GeV/c]",nPtTrig,ptt);
+  fOutputList->Add(fHPtTrg);
+  fHCent = new TH1F("fHCent","Cent;bins",nCent,cent);
+  fOutputList->Add(fHCent);
+  fHZvtx = new TH1F("fHZvtx","Zvertex;bins",nZvtx,zvtx);
+  fOutputList->Add(fHZvtx);
+
+  fNbins = nPtTrig*nPtAssc*nCent*nZvtx;
+  fHSs = new TH2*[fNbins];
+  fHMs = new TH2*[fNbins];
+
+  fIndex = new TFormula("GlobIndex","(t-1)*[0]*[1]*[2]+(z-1)*[0]*[1]+(x-1)*[0]+(y-1)+0*[4]");
+  fIndex->SetParameters(nPtTrig,nPtAssc,nZvtx,nCent);
+  fIndex->SetParNames("NTrigBins","NAssocBins", "NZvertexBins", "NCentBins");
+  //fOutputList->Add(fIndex);
+
+  fMeanPtTrg = new TProfile2D*[nPtTrig*nPtAssc];
+  fMeanPtAss = new TProfile2D*[nPtTrig*nPtAssc];
+  fMean2PtTrg = new TProfile2D*[nPtTrig*nPtAssc];
+  fMean2PtAss = new TProfile2D*[nPtTrig*nPtAssc];
+  for (Int_t c=1; c<=nCent; ++c) {
+    TString title(Form("cen=%d (%.1f)",c,fHCent->GetBinCenter(c)));
+    fMeanPtTrg[c-1]  = new TProfile2D(Form("hMeanPtTrgCen%d",c),title,nPtTrig,ptt,nPtAssc,pta);
+    fMeanPtAss[c-1]  = new TProfile2D(Form("hMeanPtAssCen%d",c),title,nPtTrig,ptt,nPtAssc,pta);
+    fMean2PtTrg[c-1] = new TProfile2D(Form("hMean2PtTrgCen%d",c),title,nPtTrig,ptt,nPtAssc,pta);
+    fMean2PtAss[c-1] = new TProfile2D(Form("hMean2PtAssCen%d",c),title,nPtTrig,ptt,nPtAssc,pta);
+    fOutputList->Add(fMeanPtTrg[c-1]);
+    fOutputList->Add(fMeanPtAss[c-1]);
+    fOutputList->Add(fMean2PtTrg[c-1]);
+    fOutputList->Add(fMean2PtAss[c-1]);
+  }
+
+  Int_t count = 0;
+  for (Int_t c=1; c<=nCent; ++c) {
+    for (Int_t z=1; z<=nZvtx; ++z) {
+      for (Int_t t=1; t<=nPtTrig; ++t) {
+       for (Int_t a=1; a<=nPtAssc; ++a) {
+         fHSs[count] = 0;
+         fHMs[count] = 0;
+         if (a>t) {
+           ++count;
+           continue;
+         }
+         TString title(Form("cen=%d (%.1f), zVtx=%d (%.1f), trig=%d (%.1f), assc=%d (%.1f)",
+                            c,fHCent->GetBinCenter(c), z,fHZvtx->GetBinCenter(z),
+                            t,fHPtTrg->GetBinCenter(t),a, fHPtAss->GetBinCenter(a)));
+         fHSs[count] = new TH2F(Form("hS%d",count), Form("Signal %s",title.Data()),
+                                nDphi,-0.5*TMath::Pi(),1.5*TMath::Pi(),nDeta,-2*fEtaMax,2*fEtaMax);
+         fHMs[count] = new TH2F(Form("hM%d",count), Form("Signal %s",title.Data()),
+                                nDphi,-0.5*TMath::Pi(),1.5*TMath::Pi(),nDeta,-2*fEtaMax,2*fEtaMax);
+         fOutputList->Add(fHSs[count]);
+         fOutputList->Add(fHMs[count]);
+         if (fVerbosity>5)
+           cout << count << " " << fIndex->Eval(t,a,z,c) << ": " << title << endl;
+         ++count;
+       }
+      }
+    }
+  }
+
+  return;
+}
+
+//________________________________________________________________________
+void AliDhcTask::InitEventMixer()
+{
+  // The effective pool size in events is set by trackDepth, so more
+  // low-mult events are required to maintain the threshold than
+  // high-mult events. Centrality pools are indep. of data histogram
+  // binning, no need to match.
+
+  Int_t trackDepth = 1000;   // # tracks to fill pool
+  Int_t poolsize   = 200;    // Maximum number of events
+
+  // Centrality pools
+  Int_t nCentBins  = 12;
+  Double_t centBins[] = {0,1,2,3,4,5,10,20,30,40,50,60,90.1};
+  //Int_t nCentBins  = 1;
+  //Double_t centBins[] = {-1,100.1};
+  // Z-vertex pools
+  Int_t nZvtxBins  = 8;
+  Double_t zvtxbin[] = {-10,-6,-4,-2,0,2,4,6,10};
+
+  fPoolMgr = new KiddiePoolManager();
+  fPoolMgr->SetTargetTrackDepth(trackDepth);
+  if (fVerbosity>4)
+    fPoolMgr->SetDebug(1);
+  fPoolMgr->InitEventPools(poolsize, nCentBins, centBins, nZvtxBins, zvtxbin);
+  
+  return;
+}
+
+//________________________________________________________________________
+void AliDhcTask::UserExec(Option_t *) 
+{
+  // Main loop, called for each event.
+
+  static int iEvent = -1; ++iEvent;
+
+  if (fVerbosity>2) {
+    if (iEvent % 10 == 0) 
+      cout << iEvent << endl;
+  }
+
+  Int_t dType = -1;       // Will be set to kESD or kAOD.
+  MiniEvent* sTracks = 0; // Vector of selected MiniTracks.
+  Double_t centCL1 = -1;
+
+  LoadBranches();
+
+  // Get event pointers, check for signs of life
+  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
+  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
+  if (fESD)
+    dType = kESD;
+  else if (fAOD)
+    dType = kAOD;
+  else {
+    AliError("Neither fESD nor fAOD available");
+    return;
+  }
+
+  // Centrality, vertex, other event variables...
+  if (dType == kESD) {
+    if (!VertexOk(fESD)) {
+      if (fVerbosity > 1)
+       AliInfo(Form("Event REJECTED. z = %.1f", fZVertex));
+      return;
+    }
+    const AliESDVertex* vertex = fESD->GetPrimaryVertex();
+    fZVertex = vertex->GetZ();
+    if(fESD->GetCentrality()) {
+      fCentrality = 
+       fESD->GetCentrality()->GetCentralityPercentile("V0M");
+      centCL1 =
+       fESD->GetCentrality()->GetCentralityPercentile("CL1");
+    }
+  }
+  if (dType == kAOD) {
+    const AliAODVertex* vertex = fAOD->GetPrimaryVertex();
+    fZVertex = vertex->GetZ();
+    if (!VertexOk(fAOD)) {
+      if (fVerbosity > 1)
+       Info("Exec()", "Event REJECTED. z = %.1f", fZVertex);
+      return;
+    }
+    const AliCentrality *aodCent = fAOD->GetHeader()->GetCentralityP();
+    if (aodCent) {
+      fCentrality = aodCent->GetCentralityPercentile("V0M");
+      centCL1     = aodCent->GetCentralityPercentile("CL1");
+    }
+  }
+  
+  // Fill Event histogram
+  fHEvt->Fill(fZVertex, fCentrality);
+
+  if (fCentrality > 90. || fCentrality < 0) {
+    AliInfo(Form("Event REJECTED. fCentrality = %.1f", fCentrality));
+    return;
+  }
+
+  // Get array of selected tracks
+  if (dType == kESD) {
+    sTracks = GetESDTrax();
+  }
+  if (dType == kAOD) {
+    sTracks = GetAODTrax();
+  }
+
+  // Get pool containing tracks from other events like this one
+  KiddiePool* pool = fPoolMgr->GetEventPool(fCentrality, fZVertex);
+  if (!pool) {
+    AliWarning(Form("No pool found. Centrality %f, ZVertex %f", fCentrality, fZVertex));
+    return;
+  }
+
+  if (!pool->IsReady()) {
+    pool->UpdatePool(sTracks);
+    return;
+  }
+
+  if (pool->IsFirstReady()) {
+    // recover events missed before the pool is ready
+    Int_t nEvs = pool->GetCurrentNEvents();
+    if (nEvs>1) {
+      for (Int_t i=0; i<nEvs; ++i) {
+       MiniEvent* evI = pool->GetEvent(i);
+       for (Int_t j=0; j<nEvs; ++j) {
+         MiniEvent* evJ = pool->GetEvent(j);
+         if (i==j) {
+           Correlate(*evI, *evJ, kSameEvt);
+         } else {
+           Correlate(*evI, *evJ, kDiffEvt, 1.0);
+         }
+       }
+      }
+    }
+  } else { /* standard case: same event, then mix*/
+    Correlate(*sTracks, *sTracks, kSameEvt);  
+    Int_t nMix = pool->GetCurrentNEvents();
+    for (Int_t jMix=0; jMix<nMix; ++jMix) {
+      MiniEvent* bgTracks = pool->GetEvent(jMix);
+      Correlate(*sTracks, *bgTracks, kDiffEvt, 1.0);
+    }
+  }
+
+  pool->UpdatePool(sTracks);
+  PostData(1, fOutputList);
+  return;
+}
+
+//________________________________________________________________________
+MiniEvent* AliDhcTask::GetESDTrax() const
+{
+  // Loop twice: 1. Count sel. tracks. 2. Fill vector.
+
+  const AliESDVertex *vtxSPD = fESD->GetPrimaryVertexSPD();
+  if (!vtxSPD)
+    return 0;
+
+  Int_t nTrax = fESD->GetNumberOfTracks();
+  if (fVerbosity > 2)
+    AliInfo(Form("%d tracks in event",nTrax));
+
+  // Loop 1.
+  Int_t nSelTrax = 0;
+  TObjArray arr(nTrax);
+  arr.SetOwner(1);
+
+  for (Int_t i = 0; i < nTrax; ++i) {
+    AliESDtrack* esdtrack = fESD->GetTrack(i);
+    if (!esdtrack) {
+      AliError(Form("Couldn't get ESD track %d\n", i));
+      continue;
+    }
+    Bool_t trkOK = fEsdTrackCutsTPCOnly->AcceptTrack(esdtrack);
+    if (!trkOK)
+      continue;
+    Double_t pt = esdtrack->Pt();
+    Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
+    if (!ptOK)
+      continue;
+    Double_t eta = esdtrack->Eta();
+    if (TMath::Abs(eta) > fEtaMax)
+      continue;
+
+    // create a tpc only track
+    AliESDtrack *newtrack = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
+    if(!newtrack)
+      continue;
+    if (newtrack->Pt()<=0) {
+      delete newtrack;
+      continue;
+    }
+
+    AliExternalTrackParam exParam;
+    Bool_t relate = newtrack->RelateToVertexTPC(vtxSPD,fESD->GetMagneticField(),kVeryBig,&exParam);
+    if (!relate) {
+      delete newtrack;
+      continue;
+    }
+
+    // set the constraint parameters to the track
+    newtrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
+
+    pt = newtrack->Pt();
+    ptOK = pt >= fPtMin && pt < fPtMax;
+    if (!ptOK) {
+      delete newtrack;
+      continue;
+    }
+    eta  = esdtrack->Eta();
+    if (TMath::Abs(eta) > fEtaMax) {
+      delete newtrack;
+      continue;
+    }
+    arr.Add(newtrack);
+    nSelTrax++;
+  }
+
+  MiniEvent* miniEvt = new MiniEvent(0);
+  miniEvt->reserve(nSelTrax);
+
+  // Loop 2.
+  for (Int_t i = 0; i < nSelTrax; ++i) {
+    AliESDtrack* esdtrack = static_cast<AliESDtrack*>(arr.At(i));
+    if (!esdtrack) {
+      AliError(Form("Couldn't get ESD track %d\n", i));
+      continue;
+    }
+    Double_t pt = esdtrack->Pt();
+    Double_t eta  = esdtrack->Eta();
+    Double_t phi  = esdtrack->Phi();
+    Int_t    sign = esdtrack->Charge() > 0 ? 1 : -1;
+    miniEvt->push_back(MiniTrack(pt, eta, phi, sign));
+  }
+  return miniEvt;
+}
+
+//________________________________________________________________________
+MiniEvent* AliDhcTask::GetAODTrax() const
+{
+  // Loop twice: 1. Count sel. tracks. 2. Fill vector.
+
+  Int_t nTrax = fAOD->GetNumberOfTracks();
+  Int_t nSelTrax = 0;
+
+  if (fVerbosity > 2)
+    AliInfo(Form("%d tracks in event",nTrax));
+
+  // Loop 1.
+  for (Int_t i = 0; i < nTrax; ++i) {
+    AliAODTrack* aodtrack = fAOD->GetTrack(i);
+    if (!aodtrack) {
+      AliError(Form("Couldn't get AOD track %d\n", i));
+      continue;
+    }
+    // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
+    UInt_t tpcOnly = 1 << 7;
+    Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly);
+    if (!trkOK)
+      continue;
+    Double_t pt = aodtrack->Pt();
+    Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
+    if (!ptOK)
+      continue;
+    Double_t eta = aodtrack->Eta();
+    if (TMath::Abs(eta) > fEtaMax)
+      continue;
+    nSelTrax++;
+  }
+
+  MiniEvent* miniEvt = new MiniEvent(0);
+  miniEvt->reserve(nSelTrax);
+
+  // Loop 2.  
+  for (Int_t i = 0; i < nTrax; ++i) {
+    AliAODTrack* aodtrack = fAOD->GetTrack(i);
+    if (!aodtrack) {
+      AliError(Form("Couldn't get AOD track %d\n", i));
+      continue;
+    }
+    
+    // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
+    UInt_t tpcOnly = 1 << 7;
+    Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly);
+    if (!trkOK)
+      continue;
+    Double_t pt = aodtrack->Pt();
+    Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
+    if (!ptOK)
+      continue;
+    Double_t eta  = aodtrack->Eta();
+    if (TMath::Abs(eta) > fEtaMax)
+      continue;
+
+    Double_t phi  = aodtrack->Phi();
+    Int_t    sign = aodtrack->Charge() > 0 ? 1 : -1;
+    miniEvt->push_back(MiniTrack(pt, eta, phi, sign));
+  }
+  return miniEvt;
+}
+
+//________________________________________________________________________
+Double_t AliDhcTask::DeltaPhi(Double_t phia, Double_t phib,
+                             Double_t rangeMin, Double_t rangeMax) const
+{
+  Double_t dphi = -999;
+  Double_t pi = TMath::Pi();
+  
+  if (phia < 0)         phia += 2*pi;
+  else if (phia > 2*pi) phia -= 2*pi;
+  if (phib < 0)         phib += 2*pi;
+  else if (phib > 2*pi) phib -= 2*pi;
+  dphi = phib - phia;
+  if (dphi < rangeMin)      dphi += 2*pi;
+  else if (dphi > rangeMax) dphi -= 2*pi;
+  
+  return dphi;
+}
+
+//________________________________________________________________________
+Int_t AliDhcTask::Correlate(const MiniEvent &evt1, const MiniEvent &evt2, 
+                           Int_t pairing, Double_t /*weight*/)
+{
+  // Triggered angular correlations. If pairing is kSameEvt, particles
+  // within evt1 are correlated. If kDiffEvt, correlate triggers from
+  // evt1 with partners from evt2.
+
+  Int_t cbin = fHCent->FindBin(fCentrality);
+  if (fHCent->IsBinOverflow(cbin) ||
+      fHCent->IsBinUnderflow(cbin))
+    return 0;
+
+  Int_t zbin = fHZvtx->FindBin(fZVertex);
+  if (fHZvtx->IsBinOverflow(zbin) ||
+      fHZvtx->IsBinUnderflow(zbin))
+    return 0;
+
+  Int_t iMax = evt1.size();
+  Int_t jMax = evt2.size();
+
+  TH2  **hist = fHMs;
+  if (pairing == kSameEvt) {
+    hist = fHSs;
+    fHCent->AddBinContent(cbin);
+    fHZvtx->AddBinContent(zbin);
+  }
+
+  Int_t nZvtx = fHZvtx->GetNbinsX();
+  Int_t nPtTrig = fHPtTrg->GetNbinsX();
+  Int_t nPtAssc = fHPtAss->GetNbinsX();
+
+  Int_t globIndex = (cbin-1)*nZvtx*nPtTrig*nPtAssc+(zbin-1)*nPtTrig*nPtAssc;
+
+  for (Int_t i=0; i<iMax; ++i) {
+
+    // Trigger particles
+    const MiniTrack &a(evt1.at(i));
+
+    Float_t pta  = a.Pt();
+    Int_t abin = fHPtTrg->FindBin(pta);
+    if (fHPtTrg->IsBinOverflow(abin) ||
+       fHPtTrg->IsBinUnderflow(abin))
+      continue;
+
+    if (pairing == kSameEvt) {
+      fHistPt->Fill(pta);
+      fHTrk->Fill(a.Phi(),a.Eta());
+      fHPtTrg->AddBinContent(abin);
+    }
+
+    for (Int_t j=0; j<jMax; ++j) {
+      // Associated particles
+      if (pairing == kSameEvt && i==j)
+       continue;
+
+      const MiniTrack &b(evt2.at(j));
+      
+      Float_t ptb  = b.Pt();
+      if (pta < ptb) 
+       continue;
+
+      Int_t bbin = fHPtTrg->FindBin(ptb);
+      if (fHPtAss->IsBinOverflow(bbin) ||
+         fHPtAss->IsBinUnderflow(bbin))
+       continue;
+
+      Float_t dphi = DeltaPhi(a.Phi(), b.Phi());
+      Float_t deta = a.Eta() - b.Eta();
+
+      Int_t index = globIndex+(abin-1)*nPtAssc+(bbin-1);
+      hist[index]->Fill(dphi,deta);
+
+      if (pairing == kSameEvt) {
+       fHPtAss->AddBinContent(bbin);
+       fMeanPtTrg[cbin-1]->Fill(pta,ptb,pta);
+       fMeanPtAss[cbin-1]->Fill(pta,ptb,ptb);
+       fMean2PtTrg[cbin-1]->Fill(pta,ptb,pta*pta);
+       fMean2PtAss[cbin-1]->Fill(pta,ptb,ptb*ptb);
+      }
+    }
+  }
+
+  return 0;
+}
+
+//________________________________________________________________________
+void AliDhcTask::Terminate(Option_t *) 
+{
+  // Draw result to the screen
+  // Called once at the end of the query
+
+  delete fPoolMgr;
+
+  fHCent->SetEntries(fHCent->Integral());
+  fHZvtx->SetEntries(fHZvtx->Integral());
+  fHPtTrg->SetEntries(fHPtTrg->Integral());
+  fHPtAss->SetEntries(fHPtAss->Integral());
+
+  if (gROOT->IsBatch())
+    return;
+
+  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
+  if (!fOutputList) {
+    AliError("Output list not available");
+    return;
+  }
+  
+  fHistPt = dynamic_cast<TH1F*> (fOutputList->At(0));
+  if (!fHistPt) {
+    AliError("ERROR: fHistPt not available\n");
+    return;
+  }
+   
+  TCanvas *c1 = new TCanvas("AliDhcTask","Pt",10,10,510,510);
+  c1->cd(1)->SetLogy();
+  fHistPt->DrawCopy("E");
+}
+
+//________________________________________________________________________
+Bool_t AliDhcTask::VertexOk(TObject* obj) const
+{
+  // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
+  
+  Int_t nContributors  = 0;
+  Double_t zVertex     = 999;
+  TString name("");
+  
+  if (obj->InheritsFrom("AliESDEvent")) {
+    AliESDEvent* esdevt = (AliESDEvent*) obj;
+    const AliESDVertex* vtx = esdevt->GetPrimaryVertex();
+    if (!vtx)
+      return 0;
+    nContributors = vtx->GetNContributors();
+    zVertex       = vtx->GetZ();
+    name          = vtx->GetName();
+  }
+  else if (obj->InheritsFrom("AliAODEvent")) { 
+    AliAODEvent* aodevt = (AliAODEvent*) obj;
+    if (aodevt->GetNumberOfVertices() < 1)
+      return 0;
+    const AliAODVertex* vtx = aodevt->GetPrimaryVertex();
+    nContributors = vtx->GetNContributors();
+    zVertex       = vtx->GetZ();
+    name          = vtx->GetName();
+  }
+  
+  // Reject if TPC-only vertex
+  if (name.CompareTo("TPCVertex")==0)
+    return kFALSE;
+  
+  // Check # contributors and range...
+  if( nContributors < 1 || TMath::Abs(zVertex) > fZVtxMax ) {
+    return kFALSE;
+  }
+  
+  return kTRUE;
+}
diff --git a/PWG4/UserTasks/DiHadronCorrelations/AliDhcTask_opt.h b/PWG4/UserTasks/DiHadronCorrelations/AliDhcTask_opt.h
new file mode 100644 (file)
index 0000000..2ce87ec
--- /dev/null
@@ -0,0 +1,92 @@
+// Dihadron correlations task - simple task to read ESD or AOD input,
+// calculate same- and mixed-event correlations, and fill THnSparse
+// output. -A. Adare, Apr 2011
+
+
+#ifndef AliDhcTask_cxx
+#define AliDhcTask_cxx
+
+class TFormula;
+class TH1;
+class TH2;
+class TObjArray;
+class TObject;
+class TProfile2D;
+class AliAODEvent;
+class AliESDEvent;
+class AliESDtrackCuts;
+class KiddiePoolManager;
+
+#include "AliAnalysisTaskSE.h"
+#include "KiddiePoolClasses.h"
+
+class AliDhcTask : public AliAnalysisTaskSE {
+ public:
+  AliDhcTask() : 
+    AliAnalysisTaskSE(), fVerbosity(0), fEtaMax(1), fZVtxMax(10), fPtMin(0.25), fPtMax(15),
+    fESD(0), fAOD(0), fOutputList(0), fHistPt(0), fHEvt(0), fHTrk(0), fHPtAss(0), 
+    fHPtTrg(0), fHCent(0), fHZvtx(0), fNbins(0), fHSs(0), fHMs(0), fIndex(0), 
+    fCentrality(99), fZVertex(99), fEsdTrackCutsTPCOnly(0), fPoolMgr(0) {}
+
+  AliDhcTask(const char *name);
+  virtual ~AliDhcTask() {}
+  
+  void         SetVerbosity(Int_t v)                  { fVerbosity = v;         }
+  void         SetPtRange(Double_t min, Double_t max) { fPtMin=min; fPtMax=max; }
+  void         SetEtaMax(Double_t eta)                { fEtaMax = eta;          }
+  void         SetZvtx(Double_t zvtx)                 { fZVtxMax = zvtx;        }
+
+ protected:
+  enum ePairingScheme {kSameEvt, kDiffEvt};
+  enum eDataType      {kESD, kAOD};
+
+  void         BookHistos();
+  void         InitEventMixer();
+  MiniEvent*   GetESDTrax() const;
+  MiniEvent*   GetAODTrax() const;
+  Bool_t       VertexOk(TObject* obj) const;
+  Double_t     DeltaPhi(Double_t phia, Double_t phib, 
+                       Double_t rangeMin = -TMath::Pi()/2, 
+                       Double_t rangeMax = 3*TMath::Pi()/2) const;
+  Int_t        Correlate(const MiniEvent &arr1, const MiniEvent &arr2, 
+                        Int_t pairing = kSameEvt, Double_t weight = 1.);
+  void         UserCreateOutputObjects();
+  void         UserExec(Option_t *option);
+  void         Terminate(Option_t *);
+
+ private:
+  Int_t        fVerbosity;       //  0 = silence
+  Double_t     fEtaMax;          //  Max |eta| cut (cm)
+  Double_t     fZVtxMax;         //  Max |z| cut (cm)
+  Double_t     fPtMin;           //  Min pt cut
+  Double_t     fPtMax;           //  Max pt cut
+  AliESDEvent *fESD;             //! ESD object
+  AliAODEvent *fAOD;             //! AOD object
+  TList       *fOutputList;      //! Output list
+  TH1F        *fHistPt;          //! Pt spectrum
+  TH2         *fHEvt;            //! Cent, vtx, etc.
+  TH2         *fHTrk;            //! Phi, Eta, etc.
+  TH1         *fHPtAss;          //! Pt ass 
+  TH1         *fHPtTrg;          //! Pt trg
+  TH1         *fHCent;           //! Centrality
+  TH1         *fHZvtx;           //! Zvertex
+  Int_t        fNbins;           //! Number of histogram bins
+  TH2        **fHSs;             //! Same-evt correlations
+  TH2        **fHMs;             //! Diff-evt correlations
+  TFormula    *fIndex;           //! Index for histograms
+  TProfile2D **fMeanPtTrg;       //! Mean pt trig 
+  TProfile2D **fMeanPtAss;       //! Mean pt ass
+  TProfile2D **fMean2PtTrg;      //! RMS pt trig 
+  TProfile2D **fMean2PtAss;      //! RMS pt ass
+  Double_t     fCentrality;      //! V0M for now
+  Double_t     fZVertex;         //! Of current event
+  AliESDtrackCuts   *fEsdTrackCutsTPCOnly; //! Track cuts
+  KiddiePoolManager *fPoolMgr;             //! Event mixer
+
+  AliDhcTask(const AliDhcTask&);            // not implemented
+  AliDhcTask &operator=(const AliDhcTask&); // not implemented
+
+  ClassDef(AliDhcTask, 2);
+};
+
+#endif
diff --git a/PWG4/UserTasks/DiHadronCorrelations/CheckOutput.C b/PWG4/UserTasks/DiHadronCorrelations/CheckOutput.C
new file mode 100644 (file)
index 0000000..268b60d
--- /dev/null
@@ -0,0 +1,143 @@
+// "$MYDATA/gsi/gsi.root"
+
+const char* fileName = "$MYDATA/gsi02/gsi02a.root";
+const char* listName = "Cont_cloizide_DhcAna";
+
+// const char* fileName = "dhctask_output_etamax1.root";
+// const char* listName = "dhclist";
+
+// From AliDhcTask.h...keep current!
+const Int_t nPairAxes = 6;
+const Int_t nEvtAxes  = 3;
+const Int_t nTrkAxes  = 2;
+//const Int_t nExtra    = 3;
+
+TString sCorrAxes[] = {"#Delta#eta", "assoc. p_{T}", "trig. p_{T}", 
+                      "Centrality", "#Delta#phi (rad)",
+                      "z-vertex (cm)", "Charge Comb."};
+TString sEvtAxes[] = {"z-vertex", "V0M Centrality (%)", "CL1 Centrality"};
+TString sTrkAxes[] = {"phi", "eta"};
+
+enum ePairHistAxes  {kDeta, kPtAssc, kPtTrig, kCent, kDphi,
+                    kZvtx};
+enum eEventHistAxes {kZvtxEvt, kCentV0M, kCentCL1};
+enum eTrackHistAxes {kPhiTrk, kEtaTrk};
+
+TObjArray* hists = new TObjArray();
+
+//gROOT->LoadMacro("../../../common/ProjectionUtils.C+");
+void CheckOutput()
+{
+  gStyle->SetOptTitle(0);
+
+  TFile* inFile = new TFile(fileName, "read");
+  TList* dhcList = inFile->Get(listName);
+  if (!dhcList) {
+    Error("","No List %s in %s", fileName, listName);
+    return;
+  }
+
+  THnSparse* hs = dhcList->FindObject("fHS");
+  THnSparse* hm = dhcList->FindObject("fHM");
+
+  if (!hs) {
+    Error("","No hs in list %s", listName);
+    return;
+  }
+
+  TH2F* hEvt = dhcList->FindObject("fHEvt");
+  TH2F* hTrk = dhcList->FindObject("fHTrk");
+  hists->Add(hEvt);
+  hists->Add(hTrk);
+  TH1D* heta = hTrk->ProjectionY();
+  hists->Add(heta);
+
+  TH1D* hzvtx = hEvt->ProjectionX("hzvtx", 1, 90);
+  hists->Add(hzvtx);
+  TH1D* hcent = hEvt->ProjectionY("hcent");
+  hists->Add(hcent);
+  // Correlation projections
+
+  // hs->GetAxis(kZvtx)->SetRangeUser(-1.9, 1.9);
+  // hm->GetAxis(kZvtx)->SetRangeUser(-1.9, 1.9);
+
+  // hs->GetAxis(kCent)->SetRangeUser(0, 4.99);
+  // hm->GetAxis(kCent)->SetRangeUser(0, 4.99);
+
+  TH2D* sig = hs->Projection(kDeta,kDphi); sig->SetName("sig");
+  TH2D* bkg = hm->Projection(kDeta,kDphi); sig->SetName("bkg");
+  sig->Sumw2(); bkg->Sumw2();
+  sig->Scale(1./sig->Integral()); 
+  bkg->Scale(1./bkg->Integral()); 
+  sig->Divide(bkg);
+  sig->GetYaxis()->SetRangeUser(-1.999, 1.999); 
+  TH2D* hcz = hs->Projection(kZvtx, kCent);
+  hcz->SetName("zvtx_vs_cent");
+
+  // 1D
+  TH1D* hzr = 0;
+
+  TH1D *hsVar[nPairAxes], *hmVar[nPairAxes], *evtVar[nEvtAxes], 
+    *trkVar[nTrkAxes];
+
+  for (int i=0; i<nPairAxes; i++) {
+
+    hsVar[i] = hs->Projection(i);
+    hsVar[i]->SetFillColor(kRed-7);
+    hsVar[i]->Sumw2();
+    hsVar[i]->SetTitle(Form("title;%s;same-evt. pairs", sCorrAxes[i].Data()));
+    hmVar[i] = hm->Projection(i);
+    hmVar[i]->SetFillColor(kAzure-9);
+    hmVar[i]->Sumw2();
+    hmVar[i]->SetTitle(Form("title;%s;mixed-evt. pairs", sCorrAxes[i].Data()));
+    hists->Add(hsVar[i]);
+    hists->Add(hmVar[i]);
+    if(i==kZvtx) {
+      hzr = (TH1D*) hsVar[i]->Clone();
+      hzr->SetFillColor(kYellow);
+      hzr->Divide(hmVar[i]);
+      hzr->SetTitle("track count ratio;z-vertex (cm);same/mixed evt track multiplicity");
+      hists->Add(hzr);
+    }
+  }
+  hists->Add(sig);
+
+  /*
+  for (int i=0; i<nEvtAxes; i++) {
+    evtVar[i] = hEvt->Projection(i);
+    evtVar[i]->SetFillColor(kGray);
+    evtVar[i]->SetTitle(Form("title;%s;events", sEvtAxes[i].Data()));
+    hists->Add(evtVar[i]);
+    if (i==1) cout << evtVar[i]->Integral(1, 89) << endl;
+  }
+  for (int i=0; i<nTrkAxes; i++) {
+    trkVar[i] = hTrk->Projection(i);
+    trkVar[i]->SetTitle(Form("title;%s;tracks", sTrkAxes[i].Data()));
+    hists->Add(trkVar[i]);
+  }
+
+  //hists->Add(hcz);
+  */
+
+  for (int n=0; n<hists->GetEntries(); n++) {
+    TCanvas* c = new TCanvas(Form("c%d",n), Form("c%d",n), 
+                      10*n, 10*n, 700, 500);
+    TObject* obj = hists->At(n);
+    if (obj->InheritsFrom("TH2")) {
+      TH2* h2 = (TH2*) obj;
+      h2->Draw("surf1");
+    }
+    else if (obj->InheritsFrom("TH1")) {
+      TH1* h = (TH1*) obj;
+      h->SetLineWidth(2);
+      h->SetFillStyle(1001); // solid
+      if(h->GetFillColor()==kNone) {
+       h->SetFillColor(kGreen-9);
+      }
+      h->Draw("histf");
+      h->Draw("histepsame");
+    }
+  }
+
+  return;
+}
diff --git a/PWG4/UserTasks/DiHadronCorrelations/KiddiePoolClasses.cxx b/PWG4/UserTasks/DiHadronCorrelations/KiddiePoolClasses.cxx
new file mode 100644 (file)
index 0000000..ad521a3
--- /dev/null
@@ -0,0 +1,306 @@
+#include "KiddiePoolClasses.h"
+
+// ===================================================================
+//                             KiddiePool
+// ===================================================================
+
+ClassImp(KiddiePool)
+
+KiddiePool::~KiddiePool()
+{
+  while (!fEvents.empty()) {
+    MiniEvent* fe= fEvents.front();
+    delete fe;
+    fe = 0;
+    fEvents.pop_front();
+  }
+}
+
+void 
+KiddiePool::PrintInfo() const
+{
+  cout << Form("%20s: %d events", "Pool capacity", fMixDepth) << endl;
+  cout << Form("%20s: %d events, %d tracks", "Current size", 
+              GetCurrentNEvents(), NTracksInPool()) << endl;
+  cout << Form("%20s: %.1f to %.1f", "Sub-event mult.", fMultMin, fMultMax) << endl;
+  cout << Form("%20s: %.1f to %.1f", "Z-vtx range", fZvtxMin, fZvtxMax) << endl;
+
+  return;
+}
+
+Bool_t 
+KiddiePool::EventMatchesBin(Int_t mult, Double_t zvtx) const
+{
+  return EventMatchesBin((Double_t) mult, zvtx);
+}
+
+Bool_t 
+KiddiePool::EventMatchesBin(Double_t mult, Double_t zvtx) const
+{
+  // Lower bin limit included; upper limit excluded.
+
+  Bool_t multOK = (mult >= fMultMin && mult < fMultMax);
+  Bool_t zvtxOK = (zvtx >= fZvtxMin && zvtx < fZvtxMax);
+  return (multOK && zvtxOK);
+}
+
+Int_t 
+KiddiePool::NTracksInPool() const
+{
+  // Number of tracks for this cent, zvtx bin; possibly includes many events.
+
+  Int_t ntrk=0;
+  for (Int_t i=0; i<(Int_t)fEvents.size(); ++i) {
+    ntrk += fNTracksInEvent.at(i);
+  }
+  return ntrk;
+}
+
+Int_t
+KiddiePool::SetEventMultRange(Int_t multMin, Int_t multMax)
+{
+  fMultMin = (Double_t)multMin;
+  fMultMax = (Double_t)multMax;
+  return 0;
+}
+
+Int_t
+KiddiePool::SetEventMultRange(Double_t multMin, Double_t multMax)
+{
+  fMultMin = multMin;
+  fMultMax = multMax;
+  return 0;
+}
+
+Int_t
+KiddiePool::SetEventZvtxRange(Double_t zvtxMin, Double_t zvtxMax)
+{
+  fZvtxMin = zvtxMin;
+  fZvtxMax = zvtxMax;
+  return 0;
+}
+
+Int_t
+KiddiePool::GlobalEventIndex(Int_t j) const
+{
+  // Index returned from passing local pool event index.
+
+  if (j < 0 || j >= (Int_t)fEventIndex.size()) {
+    cout << "ERROR in KiddiePool::GlobalEventIndex(): "
+        << " Invalid index " << j << endl;
+    return -99;
+  }
+  return fEventIndex.at(j);
+}
+
+Int_t
+KiddiePool::UpdatePool(MiniEvent* miniEvt)
+{
+  // A rolling buffer (a double-ended queue) is updated by removing
+  // the oldest event, and appending the newest.
+
+  static Int_t iEvent = -1; 
+  iEvent++;
+
+  Int_t mult = miniEvt->size(); // # tracks in this (mini) event
+  Int_t nTrk = NTracksInPool();
+
+  if (nTrk < fTargetTrackDepth && ((nTrk + mult) >= fTargetTrackDepth)) 
+    fNTimes++;
+
+  // remove 0th element before appending this event
+  Bool_t removeFirstEvent = 0;
+  if (nTrk>fTargetTrackDepth) {
+    Int_t nTrksFirstEvent= fNTracksInEvent.front();
+    Int_t diff = nTrk - nTrksFirstEvent + mult;
+    if (diff>fTargetTrackDepth)
+      removeFirstEvent = 1;
+  }
+  if (removeFirstEvent) {
+    MiniEvent* oldestEvent = fEvents.front();
+    delete oldestEvent;
+    fEvents.pop_front();         // remove first track array 
+    fNTracksInEvent.pop_front(); // remove first int
+    fEventIndex.pop_front();
+  }
+
+  fNTracksInEvent.push_back(mult);
+  fEvents.push_back(miniEvt);
+  fEventIndex.push_back(iEvent);
+
+  if (fNTimes==1) {
+    fFirstFilled = kTRUE;
+    if (KiddiePool::fDebug) {
+      cout << "\nPool " << MultBinIndex() << ", " << ZvtxBinIndex() 
+           << " ready at event "<< iEvent;
+      PrintInfo();
+      cout << endl;
+    }
+    fNTimes++; // See this message exactly once/pool
+  } else {
+    fFirstFilled = kFALSE;
+  }
+
+  fWasUpdated = true;
+
+  if (KiddiePool::fDebug) {
+    cout << " Event " << fEventIndex.back();
+    cout << " PoolDepth = " << GetCurrentNEvents(); 
+    cout << " NTracksInCurrentEvent = " << " " << mult << endl;
+  }
+
+  return fEvents.size();
+}
+
+MiniEvent* 
+KiddiePool::GetEvent(Int_t i) const
+{
+  if (i<0 || i>=(Int_t)fEvents.size()) {
+    cout << "KiddiePool::GetEvent(" 
+        << i << "): Invalid index" << endl;
+    return 0x0;
+  }
+
+  MiniEvent* ev = fEvents.at(i);
+  return ev;
+}
+
+Int_t
+KiddiePool::NTracksInEvent(Int_t iEvent) const
+{
+  // Return number of tracks in iEvent, which is the local pool index.
+
+  Int_t n = -1;
+  Int_t curEvent = fEventIndex.back();
+  Int_t offset = curEvent - iEvent;
+  Int_t pos = fEventIndex.size() - offset - 1;
+
+  if (offset==0)
+    n = fNTracksInEvent.back();
+  else if (offset < 0 || iEvent < 0) {
+    n = 0;
+  }
+  else if (offset > 0 && offset <= (int)fEventIndex.size()) {
+    n = fNTracksInEvent.at(pos);
+  }
+  else
+    cout << "Event info no longer in memory" << endl;
+  return n;
+}
+
+// ===================================================================
+//                          KiddiePoolManager
+// ===================================================================
+
+
+ClassImp(KiddiePoolManager)
+
+KiddiePoolManager::KiddiePoolManager(Int_t depth,     Int_t minNTracks,
+                                    Int_t nMultBins, Double_t *multbins,
+                                    Int_t nZvtxBins, Double_t *zvtxbins) 
+: fDebug(0), fNMultBins(0), fNZvtxBins(0), fEvPool(0), fTargetTrackDepth(minNTracks) 
+{
+  // Constructor.
+
+  InitEventPools(depth, nMultBins, multbins, nZvtxBins, zvtxbins);
+  cout << "KiddiePoolManager initialized." << endl;
+}
+
+KiddiePoolManager::~KiddiePoolManager()
+{
+  for (Int_t iM=0; iM<fNMultBins; iM++) {
+    for (Int_t iZ=0; iZ<fNZvtxBins; iZ++) {
+      KiddiePool* pool = fEvPool.at(iM).at(iZ);
+      delete pool;
+    }
+  }
+}
+
+Int_t KiddiePoolManager::InitEventPools(Int_t depth, 
+                                       Int_t nMultBins, Double_t *multbin, 
+                                       Int_t nZvtxBins, Double_t *zvtxbin)
+{
+  // Assign KiddiePoolManager members.
+
+  fNMultBins = nMultBins;
+  fNZvtxBins = nZvtxBins;
+
+  for (Int_t iM=0; iM<fNMultBins; iM++) {
+    std::vector<KiddiePool*> evp;
+    for (Int_t iZ=0; iZ<fNZvtxBins; iZ++) {
+      evp.push_back(new KiddiePool(depth, 
+                                  multbin[iM], multbin[iM+1], 
+                                  zvtxbin[iZ], zvtxbin[iZ+1] ));
+    }
+    fEvPool.push_back(evp);
+  }
+  
+  for (Int_t iM=0; iM<nMultBins; iM++) {
+    for (Int_t iZ=0; iZ<nZvtxBins; iZ++) {
+      fEvPool.at(iM).at(iZ)->SetMultBinIndex(iM);
+      fEvPool.at(iM).at(iZ)->SetZvtxBinIndex(iZ);
+      fEvPool.at(iM).at(iZ)->SetTargetTrackDepth(fTargetTrackDepth);
+      fEvPool.at(iM).at(iZ)->SetDebug(fDebug);
+    }
+  }
+  
+  if (fDebug) {
+    cout << "fEvPool outer size: " << fEvPool.size() << endl;
+    for (Int_t iM=0; iM<nMultBins; iM++) {
+      for (Int_t iZ=0; iZ<nZvtxBins; iZ++) {
+       if(fEvPool.at(iM).at(iZ)) {
+         cout << "multiplicity bin: " << iM;
+         cout << ", z-vertex bin: " << iZ;
+         fEvPool.at(iM).at(iZ)->PrintInfo();
+       }
+      }
+    }
+  }
+  
+  return fEvPool.size();
+}
+
+KiddiePool*
+KiddiePoolManager::GetEventPool(Int_t iMult, Int_t iZvtx) const
+{
+  if (iMult < 0 || iMult >= fNMultBins) 
+    return 0x0;
+  if (iZvtx < 0 || iZvtx >= fNZvtxBins) 
+    return 0x0;
+  return fEvPool.at(iMult).at(iZvtx);
+}
+
+KiddiePool*
+KiddiePoolManager::GetEventPool(Int_t centVal, Double_t zVtxVal) const
+{
+  return GetEventPool((Double_t)centVal, zVtxVal);
+}
+
+KiddiePool*
+KiddiePoolManager::GetEventPool(Double_t centVal, Double_t zVtxVal) const
+{
+  // Return appropriate pool for this centrality and z-vertex value.
+
+  for (Int_t iM=0; iM<fNMultBins; iM++) {
+    for (Int_t iZ=0; iZ<fNZvtxBins; iZ++) {
+      KiddiePool* pool = GetEventPool(iM, iZ);
+      if (pool->EventMatchesBin(centVal, zVtxVal))
+       return pool;
+    }
+  }
+  return 0x0;
+}
+
+Int_t
+KiddiePoolManager::UpdatePools(MiniEvent* miniEvt)
+{
+  // Call UpdatePool for all bins.
+
+  for (Int_t iM=0; iM<fNMultBins; iM++) {
+    for (Int_t iZ=0; iZ<fNZvtxBins; iZ++) {
+      if (fEvPool.at(iM).at(iZ)->UpdatePool(miniEvt) > -1)
+        break;
+    }
+  }
+  return 0;
+}
diff --git a/PWG4/UserTasks/DiHadronCorrelations/KiddiePoolClasses.h b/PWG4/UserTasks/DiHadronCorrelations/KiddiePoolClasses.h
new file mode 100644 (file)
index 0000000..3de4ac8
--- /dev/null
@@ -0,0 +1,168 @@
+#ifndef KiddiePoolClasses_h
+#define KiddiePoolClasses_h
+
+#include <vector>
+#include <deque>
+#include <Rtypes.h>
+#include <Riostream.h>
+#include <TObject.h>
+#include <TFile.h>
+#include <TMath.h>
+#include <TSystem.h>
+
+class MiniTrack;
+
+// Low-memory event mixing classes:
+// - KiddiePoolManager: Maintain 2D vector (z,cent) of KiddiePools.
+// - KiddiePool: collections of sub-events in z and mult/cent bins.
+// - MiniTrack: super-lightweight track "struct" class.
+//
+// Stores a buffer of tracks that updates continuously. The track type
+// contained by the pools can be anything inheriting from
+// TObject. Pools are updated based on maintaining a minimum fixed
+// number of tracks. Multiplicity/centrality and z-vertex bins must be
+// passed in at initialization. For example of implementation, see
+// $ALICE_ROOT/PWG4/JetTasks/AliAnalysisTaskPhiCorrelations.cxx
+//
+// Authors: A. Adare and C. Loizides
+
+typedef std::vector< MiniTrack > MiniEvent;
+
+class MiniTrack
+{
+ public:
+  MiniTrack() : fPt(0), fEta(0), fPhi(0), fSign(0) {};
+  MiniTrack(Double_t pt, Double_t eta, Double_t phi, Int_t sign) :
+   fPt(pt), fEta(eta), fPhi(phi), fSign(sign>=0?1:0) {};
+  Float_t Eta()  const { return fEta; }
+  Float_t Phi()  const { return fPhi; }
+  Float_t Pt()   const { return fPt;  }
+  Int_t   Sign() const { if (fSign) return 1; return -1; }
+
+ protected: 
+  Float_t fPt;
+  Float_t fEta;
+  Float_t fPhi;
+  Bool_t  fSign;
+};
+
+class KiddiePool : public TObject
+{
+ public:
+ KiddiePool(Int_t d) : 
+   fEvents(0),
+   fNTracksInEvent(0),
+   fEventIndex(0),
+   fMixDepth(d), 
+   fMultMin(-999), 
+   fMultMax(+999), 
+   fZvtxMin(-999), 
+   fZvtxMax(+999), 
+   fWasUpdated(0), 
+   fMultBinIndex(0), 
+    fZvtxBinIndex(0), 
+   fDebug(0), 
+   fTargetTrackDepth(0),
+   fFirstFilled(0),
+   fNTimes(0) {;}
+  
+ KiddiePool(Int_t d, Double_t multMin, Double_t multMax, 
+           Double_t zvtxMin, Double_t zvtxMax) : 
+   fEvents(0),
+   fNTracksInEvent(0),
+   fEventIndex(0),
+   fMixDepth(d), 
+   fMultMin(multMin), 
+   fMultMax(multMax), 
+   fZvtxMin(zvtxMin),
+   fZvtxMax(zvtxMax),
+   fWasUpdated(0),
+    fMultBinIndex(0), 
+   fZvtxBinIndex(0),
+   fDebug(0),
+   fTargetTrackDepth(0),
+   fFirstFilled(0),
+   fNTimes(0) {;}
+  ~KiddiePool();
+  
+  Bool_t      EventMatchesBin(Int_t mult,    Double_t zvtx) const;
+  Bool_t      EventMatchesBin(Double_t mult, Double_t zvtx) const;
+  Bool_t      IsReady()                    const { return NTracksInPool() >= fTargetTrackDepth; }
+  Bool_t      IsFirstReady()               const { return fFirstFilled;   }
+  Int_t       GetNTimes()                  const { return fNTimes;        }
+  Int_t       GetCurrentNEvents()          const { return fEvents.size(); }
+  Int_t       GlobalEventIndex(Int_t j)    const;
+  MiniEvent*  GetEvent(Int_t i)            const;
+  Int_t       MultBinIndex()               const { return fMultBinIndex; }
+  Int_t       NTracksInEvent(Int_t iEvent) const;
+  Int_t       NTracksInCurrentEvent()      const { if (fNTracksInEvent.empty()) return 0;  
+                                                   return fNTracksInEvent.back(); }
+  void        PrintInfo()                  const;
+  Int_t       NTracksInPool()              const;
+  Bool_t      WasUpdated()                 const { return fWasUpdated; }
+  Int_t       ZvtxBinIndex()               const { return fZvtxBinIndex; }
+  void        SetDebug(Bool_t b)                 { fDebug = b; }
+  void        SetTargetTrackDepth(Int_t d)       { fTargetTrackDepth = d; }
+  Int_t       SetEventMultRange(Int_t    multMin, Int_t multMax);
+  Int_t       SetEventMultRange(Double_t multMin, Double_t multMax);
+  Int_t       SetEventZvtxRange(Double_t zvtxMin, Double_t zvtxMax);
+  void        SetMultBinIndex(Int_t iM)          { fMultBinIndex = iM; }
+  void        SetZvtxBinIndex(Int_t iZ)          { fZvtxBinIndex = iZ; }
+  Int_t       UpdatePool(MiniEvent* miniEvt);
+  
+protected:
+  deque< MiniEvent* >   fEvents;              //Pool of events storing Minitracks
+  deque<int>            fNTracksInEvent;      //Tracks in event
+  deque<int>            fEventIndex;          //Original event index
+  Int_t                 fMixDepth;            //Number of evts. to mix with
+  Double_t              fMultMin, fMultMax;   //Track multiplicity bin range
+  Double_t              fZvtxMin, fZvtxMax;   //Event z-vertex bin range
+  Bool_t                fWasUpdated;          //Evt. succesfully passed selection?
+  Int_t                 fMultBinIndex;        //Multiplicity bin
+  Int_t                 fZvtxBinIndex;        //Zvertex bin
+  Int_t                 fDebug;               //If 1 then debug on
+  Int_t                 fTargetTrackDepth;    //Number of tracks, once full
+  Bool_t                fFirstFilled;         //Init to false
+  Int_t                 fNTimes;              //Number of times init. condition reached
+
+  ClassDef(KiddiePool,1) // Event pool class
+};
+
+class KiddiePoolManager : public TObject
+{
+public:
+  KiddiePoolManager() : 
+    fDebug(0),
+    fNMultBins(0), 
+    fNZvtxBins(0),
+    fEvPool(0),
+    fTargetTrackDepth(0) {;}
+  KiddiePoolManager(Int_t maxEvts, Int_t minNTracks,
+                   Int_t nMultBins, Double_t *multbins,
+                   Int_t nZvtxBins, Double_t *zvtxbins);
+
+  ~KiddiePoolManager();
+  
+  // First uses bin indices, second uses the variables themselves.
+  KiddiePool *GetEventPool(Int_t iMult, Int_t iZvtx) const;
+  KiddiePool *GetEventPool(Int_t    centVal, Double_t zvtxVal) const;
+  KiddiePool *GetEventPool(Double_t centVal, Double_t zvtxVal) const;
+
+  Int_t       InitEventPools(Int_t depth, 
+                            Int_t nmultbins, Double_t *multbins, 
+                            Int_t nzvtxbins, Double_t *zvtxbins);
+  void        SetTargetTrackDepth(Int_t d) { fTargetTrackDepth = d;}
+  Int_t       UpdatePools(MiniEvent*);
+  void        SetDebug(Bool_t b) { fDebug = b; }
+
+protected:
+  Int_t      fDebug;                         // If 1 then debug on
+  Int_t      fNMultBins;                     // number mult bins
+  Int_t      fNZvtxBins;                     // number vertex bins
+  vector< vector<KiddiePool*> > fEvPool;     // pool in bins of [fMultBin][fZvtxBin]
+  Int_t      fTargetTrackDepth;              // Required track size, same for all pools.
+
+  ClassDef(KiddiePoolManager,1)
+};
+#endif
diff --git a/PWG4/UserTasks/DiHadronCorrelations/RunDhcTask.C b/PWG4/UserTasks/DiHadronCorrelations/RunDhcTask.C
new file mode 100644 (file)
index 0000000..f4d9126
--- /dev/null
@@ -0,0 +1,185 @@
+// Macro to run AliDhcTask - A. Adare, Apr 2011
+// Choose dataType = "esd" or "aod".
+// Choose runMode  = "local", "proof", "grid", or "terminate".
+// Choose proofCluster = "skaf.saske.sk" or "alice-caf.cern.ch".
+// When running on proof: use root, not aliroot.
+
+Int_t verbosity      = 0;
+Long64_t nEvents     = 200;
+TString runMode      = "local";
+TString dataType     = "aod";
+TString localDir     = "/Users/adare/esd/alice/data/2010/LHC10h/000139107/ESDs/pass2/";
+TString esdDir       = localDir + TString("10000139107001.120/");
+TString aodDir       = localDir + TString("AOD049/0008/");
+TString proofDataset = "/alice/data/LHC10h_000137848_p2_AOD049";
+//TString proofCluster = "aadare@alice-caf.cern.ch";
+TString proofCluster = "aadare@skaf.saske.sk";
+TString alirootVer   = "VO_ALICE@AliRoot::v4-21-21-AN";
+TString libsBase     = "Tree:Geom:VMC:STEERBase:ESD:AOD";
+//TString libsExtra    = "ANALYSIS:ANALYSISalice:CORRFW:JETAN:PWG4JetTasks";
+TString libsExtra    = "ANALYSIS:ANALYSISalice";
+TList* libList = new TList();
+TChain* chain = 0;
+
+void RunDhcTask()
+{
+  // Creates the large file memstat_<pid>.root.
+  // Call TMemStat::Show() on the promp to analyze it.
+  if (0)
+    TMemStat mm("gnubuiltin");
+
+  if (runMode=="local")
+    LocalSetup();
+  else if (runMode=="proof") {
+    ProofSetup();
+  }
+
+  AddDhcTask();
+  
+  /// ---  All set up now...run Analysis!  ---
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr->InitAnalysis()) 
+    return;
+  if (verbosity > 0)
+    mgr->PrintStatus();
+  
+  if (runMode=="local") {
+    mgr->StartAnalysis(runMode.Data(), chain, nEvents);
+  }
+  else if (runMode=="proof") {
+    mgr->StartAnalysis(runMode.Data(), proofDataset.Data(), nEvents);
+  }
+  return;
+}
+
+void LocalSetup()
+{
+  // Load libraries 
+  TObjArray* libs = libsExtra.Tokenize(":");
+  for (Int_t i=0; i<libs->GetEntries(); i++) {
+    const char* name = libs->At(i)->GetName();
+    Int_t ret = gSystem->Load(Form("lib%s", name));
+    if (ret > 0)
+      Warning("LocalSetup()", "lib%s already loaded", name);
+    if (ret < 0)
+      Warning("LocalSetup()", "lib%s failed to load", name);
+  }
+  
+  // Add headers
+  gSystem->AddIncludePath("-I$ALICE_ROOT/include");
+  
+  // Load macros/tasks
+  if (dataType == "esd") {
+    gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C");
+  }
+  gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+  gROOT->LoadMacro("KiddiePoolClasses.cxx+g");  
+  gROOT->LoadMacro("AliDhcTask.cxx+g");  
+  
+  // Setup chain
+  TString chainName = dataType.Contains("aod") ? "aodTree" : "esdTree";
+  chain = new TChain(chainName.Data());
+  if (dataType.Contains("aod"))
+    chain->Add(Form("%sAliAOD.root", aodDir.Data()));
+  if (dataType.Contains("esd"))
+    chain->Add(Form("%sAliESDs.root", esdDir.Data()));
+  return;
+}
+
+void ProofSetup()
+{
+  gEnv->SetValue("XSec.GSI.DelegProxy", "2");
+      
+  // Get PROOF pointer
+  TProof *proof = gProof;
+  if (proof) {
+    TProof::Reset(proofCluster.Data()); // add kTRUE as a 2nd arg for hard reset
+    //    proof->Reset("");
+    cerr << "Trying to reset proof" << endl;
+  } else {
+    proof = TProof::Open(proofCluster.Data(), "workers=60");
+  }
+  if (!proof) {
+    cerr << "Connection to " << proofCluster.Data()
+        << " failed!" << endl;
+    return;
+  }
+
+  if (0) {
+    proof->ShowPackages();
+    proof->ShowDataSets();
+  }
+  
+  // Load libraries   
+  int  ret = 0;
+  libList->Add(new TNamed("ALIROOT_EXTRA_LIBS", libsBase.Data()));
+  libList->Add(new TNamed("ALIROOT_EXTRA_LIBS", libsExtra.Data()));
+
+  ret = proof->EnablePackage(alirootVer.Data(), libList);
+  
+  proof->EnablePackage(alirootVer.Data(), 0);
+  if (ret) {
+    Error("ProofSetup()", "Failed to load all libs.");
+    return;
+  }
+
+  // Load macros/tasks
+  if (dataType == "esd") {
+    proof->Load("AddTaskCentrality.C");
+  }
+  proof->Load("AddTaskPhysicsSelection.C");
+  proof->Load("KiddiePoolClasses.cxx+g");  
+  proof->Load("AliDhcTask.cxx+g", 0);
+  return;
+}
+
+void AddDhcTask()
+{
+  // Need the following macros loaded beforehand:
+  // "$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C" (ESD only)
+  // "$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C"
+  // "AliDhcTask.cxx+g"
+
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr)
+    mgr = new AliAnalysisManager("DhcAnalysis");
+  AliESDInputHandler* esdH = new AliESDInputHandler();
+  AliAODInputHandler* aodH = new AliAODInputHandler();
+  
+  if (dataType.Contains("aod"))
+    mgr->SetInputEventHandler(aodH);
+  if (dataType.Contains("esd"))
+    mgr->SetInputEventHandler(esdH);
+  
+  if (dataType == "esd") {
+    AliCentralitySelectionTask *centTask = AddTaskCentrality();
+    centTask->SelectCollisionCandidates(AliVEvent::kAny);
+    centTask->SetPass(2);
+  }
+  
+  AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection();
+  
+  AliDhcTask *dhcTask = new AliDhcTask("DhcTask");
+  dhcTask->SelectCollisionCandidates(AliVEvent::kMB);
+  dhcTask->SetVerbosity(verbosity);
+
+  // Add task(s)
+  mgr->AddTask(dhcTask);
+  
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+  AliAnalysisDataContainer *coutput = 
+    mgr->CreateContainer("dhclist", TList::Class(),   
+                        AliAnalysisManager::kOutputContainer, 
+                        "dhctask_output.root");
+  
+  // Connect input/output
+  mgr->ConnectInput(dhcTask, 0, cinput);
+  mgr->ConnectOutput(dhcTask, 1, coutput);
+  
+  // Enable debug printouts
+  mgr->SetDebugLevel(verbosity);
+
+  return;
+}
+
diff --git a/PWG4/UserTasks/DiHadronCorrelations/rootlogon.C b/PWG4/UserTasks/DiHadronCorrelations/rootlogon.C
new file mode 100644 (file)
index 0000000..35adf6d
--- /dev/null
@@ -0,0 +1,22 @@
+{
+  cout << "Loading " << gSystem->Getenv("PWD") 
+       << "/rootlogon.C" << endl;
+  cout << "Using ROOT version " 
+       << gROOT->GetVersion() << endl;
+
+  gROOT->SetStyle("Plain");
+  gStyle->SetPalette(1);
+
+  gSystem->Load("libTree.so");
+  if (gSystem->Getenv("TMPDIR")) 
+    gSystem->SetBuildDir(gSystem->Getenv("TMPDIR"));
+  if (1) {
+      delete gRandom;
+      gRandom = new TRandom3(0);
+   }
+
+  // Include directories
+  // gSystem->AddIncludePath("-I\"common\" ");
+  // gROOT->LoadMacro("common/Utils.C+g");
+
+}