new version of the LRC code (Christoph Mayer <Christoph.Mayer@cern.ch>)
authormiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 2 Dec 2013 13:13:33 +0000 (13:13 +0000)
committermiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 2 Dec 2013 13:13:33 +0000 (13:13 +0000)
PWGCF/Correlations/DPhi/AliAnalysisTaskLongRangeCorrelations.cxx
PWGCF/Correlations/DPhi/AliAnalysisTaskLongRangeCorrelations.h
PWGCF/Correlations/macros/AddTaskLongRangeCorrelations.C

index 4de6731..d0b297e 100644 (file)
-// -*- C++ -*-\r
-/**************************************************************************\r
- * Copyright(c) 2012,      ALICE Experiment at CERN, All rights reserved. *\r
- *                                                                        *\r
- * Author: The ALICE Off-line Project.                                    *\r
- * Contributors are mentioned in the code where appropriate.              *\r
- *                                                                        *\r
- * Permission to use, copy, modify and distribute this software and its   *\r
- * documentation strictly for non-commercial purposes is hereby granted   *\r
- * without fee, provided that the above copyright notice appears in all   *\r
- * copies and that both the copyright notice and this permission notice   *\r
- * appear in the supporting documentation. The authors make no claims     *\r
- * about the suitability of this software for any purpose. It is          *\r
- * provided "as is" without express or implied warranty.                  *\r
- **************************************************************************/\r
-// $Id: AliAnalysisTaskLongRangeCorrelations.cxx 359 2013-10-29 15:39:09Z cmayer $\r
-\r
-#include <numeric>\r
-#include <functional>\r
-#include <set>\r
-\r
-#include <TChain.h>\r
-#include <THashList.h>\r
-#include <THnSparse.h>\r
-#include <TH1.h>\r
-#include <TH2.h>\r
-#include <TH3.h>\r
-\r
-#include "AliEventPoolManager.h"\r
-#include "AliAnalysisManager.h"\r
-#include "AliAODEvent.h"\r
-#include "AliAODMCHeader.h"\r
-#include "AliAODMCParticle.h"\r
-#include "AliMCEvent.h"\r
-#include "AliInputEventHandler.h"\r
-#include "AliLog.h"\r
-#include "AliTHn.h"\r
-\r
-#include "AliMCEventHandler.h"\r
-#include "AliMCEvent.h"\r
-#include "AliStack.h"\r
-#include "AliAnalysisTaskLongRangeCorrelations.h"\r
-\r
-\r
-ClassImp(AliAnalysisTaskLongRangeCorrelations);\r
-ClassImp(LRCParticle);\r
-\r
-AliAnalysisTaskLongRangeCorrelations::AliAnalysisTaskLongRangeCorrelations(const char *name)\r
-  : AliAnalysisTaskSE(name)\r
-  , fOutputList(NULL)\r
-  , fVertexZ(NULL)\r
-  , fRunMixing(kFALSE)\r
-  , fPoolMgr(NULL)\r
-  , fMixingTracks(50000)\r
-  , fTrackFilter(128)\r
-  , fCentMin(0), fCentMax(20)\r
-  , fPtMin(0.2), fPtMax(1e10)\r
-  , fPhiMin(0.), fPhiMax(TMath::TwoPi())\r
-  , fMaxAbsVertexZ(10.)\r
-  , fSelectPrimaryMCParticles(0)\r
-  , fSelectPrimaryMCDataParticles(0)\r
-  , fNMin(-1)\r
-  , fNMax(-1)\r
-  , fnBinsCent( 220), fnBinsPt(400), fnBinsPhi(4),              fnBinsEta(120)\r
-  , fxMinCent( -5.0), fxMinPt( 0.0), fxMinPhi( 0.0),            fxMinEta(-1.5)\r
-  , fxMaxCent(105.0), fxMaxPt( 4.0), fxMaxPhi( TMath::TwoPi()), fxMaxEta( 1.5) {\r
-  DefineInput(0, TChain::Class());\r
-  DefineOutput(1, TList::Class());\r
-}\r
-\r
-AliAnalysisTaskLongRangeCorrelations::~AliAnalysisTaskLongRangeCorrelations() {\r
-  if (NULL != fOutputList) {\r
-    delete fOutputList;\r
-    fOutputList = NULL;\r
-  }\r
-}\r
-\r
-void AliAnalysisTaskLongRangeCorrelations::UserCreateOutputObjects() {\r
-  fOutputList = new THashList;\r
-  fOutputList->SetOwner(kTRUE);\r
-  fOutputList->SetName(GetOutputListName());\r
-\r
-  const Double_t vertexBins[] = {  // Zvtx bins\r
-    -10., -7., -5., -3., -1., 1., 3., 5., 7., 10.\r
-  };\r
-  const Int_t nVertexBins(sizeof(vertexBins)/sizeof(Double_t)-1);\r
-  fVertexZ = new TAxis(nVertexBins, vertexBins);\r
-  fVertexZ->SetName("VertexZAxis");\r
-  fOutputList->Add(fVertexZ);\r
-\r
-  fOutputList->Add(new TH1D("histVertexZ", ";vertex Z (cm)", nVertexBins, vertexBins));\r
-\r
-  // Event statistics\r
-  const char* eventStatLabels[] = { \r
-    "All Events",\r
-    "Physics Selection",\r
-    "Centrality Selection",\r
-    "Vertex Selection",\r
-    "Analyzed Events"\r
-  };\r
-  const size_t nEventStat(sizeof(eventStatLabels)/sizeof(const char*));\r
-  TH1* hStats(new TH1D("histEventStats", "histEventStats", nEventStat, -0.5, nEventStat-0.5));\r
-  for (size_t i(0); i<nEventStat; ++i)\r
-    hStats->GetXaxis()->SetBinLabel(1+i, eventStatLabels[i]);\r
-  fOutputList->Add(hStats);\r
-\r
-  // QA histograms\r
-  fOutputList->Add(new TH2D("histQACentrality", ";centrality V0M(%);selected;",\r
-                           fnBinsCent,fxMinCent,fxMaxCent, 2, -0.5, 1.5));\r
-  fOutputList->Add(new TH2D("histQAVertexZ", ";vertex-z (cm);selected;",\r
-                           800, -40., 40., 2, -0.5, 1.5));\r
-  fOutputList->Add(new TH1D("histQAMultiplicityBeforeCuts",\r
-                           ";tracks", 10000, 0, 10000));\r
-  fOutputList->Add(new TH1D("histQAMultiplicityAfterCuts",\r
-                           ";selected tracks", 10000, 0, 10000));\r
-  fOutputList->Add(new TH3D("histQACentPt", ";charge;centrality V0M(%);p_{T} (GeV/c)",\r
-                           2, -0.5, 1.5, fnBinsCent, fxMinCent, fxMaxCent, fnBinsPt, fxMinPt, fxMaxPt));\r
-  fOutputList->Add(new TH3D("histQAPhiEta", ";charge;#phi (rad);#eta",\r
-                           2, -0.5, 1.5, 200, 0.0, TMath::TwoPi(), 300, -1.5, 1.5));\r
-\r
-  // N(eta) distributions with different binnings\r
-  fOutputList->Add(new TH2D("histNEta_120", ";#eta;N", 120, -1.5, 1.5, 1000, -.5, 1000-.5));  // 0.025\r
-  fOutputList->Add(new TH2D("histNEta__60", ";#eta;N",  60, -1.5, 1.5, 1000, -.5, 1000-.5));  // 0.05\r
-  fOutputList->Add(new TH2D("histNEta__30", ";#eta;N",  30, -1.5, 1.5, 1000, -.5, 1000-.5));  // 0.1\r
-  fOutputList->Add(new TH2D("histNEta__15", ";#eta;N",  15, -1.5, 1.5, 1000, -.5, 1000-.5));  // 0.2\r
-\r
-  // Moments\r
-  fOutputList->Add(MakeHistPhiEta("histMoment1PhiEta_1"));\r
-  if (fRunMixing)\r
-    fOutputList->Add(MakeHistPhiEta("histMoment1PhiEta_2"));\r
-  fOutputList->Add(MakeHistPhiEtaPhiEta("histMoment2PhiEtaPhiEta_11"));\r
-  if (fRunMixing) {\r
-    fOutputList->Add(MakeHistPhiEtaPhiEta("histMoment2PhiEtaPhiEta_12"));\r
-    fOutputList->Add(MakeHistPhiEtaPhiEta("histMoment2PhiEtaPhiEta_22"));\r
-  }\r
-  \r
-  // add MC Histograms\r
-  const Int_t N(fOutputList->GetEntries());\r
-  for (Int_t i(0); i<N; ++i)\r
-    fOutputList->Add(fOutputList->At(i)->Clone(TString("MC_") + fOutputList->At(i)->GetName()));\r
-\r
-  if (fRunMixing)\r
-    SetupForMixing();\r
-\r
-  PostData(1, fOutputList);\r
-}\r
-\r
-void AliAnalysisTaskLongRangeCorrelations::UserExec(Option_t* ) {\r
-  AliAnalysisManager* pManager(AliAnalysisManager::GetAnalysisManager());\r
-  if (NULL == pManager) return;\r
-\r
-  AliInputEventHandler* pInputHandler(dynamic_cast<AliInputEventHandler*>(pManager->GetInputEventHandler()));\r
-  if (NULL == pInputHandler) return;\r
-\r
-  AliAODEvent* pAOD(dynamic_cast<AliAODEvent*>(InputEvent()));\r
-  if (NULL == pAOD) return;\r
-\r
-  AliAODHeader *pAODHeader = pAOD->GetHeader();\r
-  if (NULL == pAODHeader) return;\r
-\r
-  AliAODMCHeader *pAODMCHeader(dynamic_cast<AliAODMCHeader*>(pAOD->FindListObject(AliAODMCHeader::StdBranchName())));\r
-  const Bool_t isMC(NULL != pAODMCHeader);\r
-\r
-  TClonesArray *arrayMC(NULL);\r
-  if (isMC) {\r
-    arrayMC = dynamic_cast<TClonesArray*>(pAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName()));\r
-    if (NULL == arrayMC) return;\r
-  }\r
-\r
-  // --- event cuts MC ---\r
-  if (isMC) {\r
-    Fill("MC_histEventStats", 0.); // all events\r
-    Fill("MC_histEventStats", 1.); // events passing physics selection\r
-  }\r
-\r
-  // --- event cuts data ---\r
-  Fill("histEventStats", 0.); // all events\r
-\r
-  // physics selection\r
-  if (!pInputHandler->IsEventSelected()) return;\r
-  Fill("histEventStats", 1.); // events passing physics selection\r
-\r
-  // centrality selection\r
-  const Double_t centrality(pAODHeader->GetCentralityP()->GetCentralityPercentile("V0M"));\r
-  AliDebug(3, Form("centrality=%f", centrality));\r
-  const Bool_t centralitySelected(centrality > fCentMin && centrality < fCentMax);\r
-  Fill("histQACentrality", centrality, centralitySelected);\r
-  if (!centralitySelected) return;\r
-  Fill("histEventStats", 2.); // events passing centrality selection\r
-  if (isMC)\r
-    Fill("MC_histEventStats", 2.); // events passing centrality selection\r
-\r
-  // vertex selection -- data\r
-  const Int_t nVertex(pAOD->GetNumberOfVertices());\r
-  if (0 == nVertex) return;\r
-  const AliAODVertex* pVertex(pAOD->GetPrimaryVertex());\r
-  if (NULL == pVertex) return;\r
-  const Int_t nTracksPrimary(pVertex->GetNContributors());\r
-  if (nTracksPrimary < 1) return;\r
-\r
-  const Double_t zVertex(pVertex->GetZ());\r
-  const Bool_t vertexSelectedData(TMath::Abs(zVertex) < fMaxAbsVertexZ);\r
-\r
-  // vertex selection -- MC\r
-  Bool_t vertexSelectedMC(kTRUE);\r
-  if (isMC)\r
-    vertexSelectedMC = (TMath::Abs(pAODMCHeader->GetVtxZ()) < fMaxAbsVertexZ);\r
-\r
-  // combined vertex selection (data and MC)\r
-  const Bool_t vertexSelected(vertexSelectedData && vertexSelectedMC);\r
-  Fill("histQAVertexZ", zVertex, vertexSelected);\r
-  if (isMC)\r
-    Fill("MC_histQAVertexZ", pAODMCHeader->GetVtxZ(), vertexSelected);\r
-  if (!vertexSelected) return;\r
-  Fill("histEventStats",    3.); // events passing vertex selection\r
-  if (isMC)\r
-    Fill("MC_histEventStats", 3.); // events passing vertex selection\r
-\r
-  // ------------------\r
-  // event is accepted\r
-  // ------------------\r
-\r
-  TObjArray* tracksMain(GetAcceptedTracks(pAOD, arrayMC, centrality));\r
-  Fill("histQAMultiplicityBeforeCuts", pAOD->GetNumberOfTracks());\r
-  Fill("histQAMultiplicityAfterCuts", tracksMain->GetEntriesFast());\r
-  \r
-  if (fRunMixing) {  \r
-    AliEventPool* pEventPool(fPoolMgr->GetEventPool(centrality, pAOD->GetPrimaryVertex()->GetZ()));\r
-    if (NULL == pEventPool)\r
-      AliFatal(Form("No pool found for centrality = %f, zVtx = %f", centrality, pAOD->GetPrimaryVertex()->GetZ()));\r
-    \r
-//     pEventPool->PrintInfo();\r
-    if (pEventPool->IsReady()\r
-       || pEventPool->NTracksInPool()     > fMixingTracks/10\r
-       || pEventPool->GetCurrentNEvents() >= 5) {\r
-      Fill("histEventStats", 4.); // analyzed events\r
-      Fill("histVertexZ", zVertex);\r
-      const Int_t nMix(pEventPool->GetCurrentNEvents());\r
-      for (Int_t i(0); i<nMix; ++i) {\r
-       TObjArray* tracksMixed(pEventPool->GetEvent(i));\r
-       CalculateMoments("", tracksMain, tracksMixed, zVertex, 1./nMix);\r
-      }\r
-    }\r
-    // Update the Event pool\r
-    pEventPool->UpdatePool(tracksMain);\r
-  } else { // no mixing\r
-    Fill("histEventStats", 4.); // analyzed events\r
-    Fill("histVertexZ", zVertex);\r
-    CalculateMoments("", tracksMain, tracksMain, zVertex, 1.);\r
-    delete tracksMain;\r
-  }\r
-  \r
-  if (isMC) {\r
-    TObjArray* tracksMC(GetAcceptedTracks(arrayMC, centrality));\r
-    const Double_t x[2] = {\r
-      arrayMC->GetEntriesFast(),\r
-      tracksMC->GetEntriesFast()\r
-    };\r
-    Fill("MC_histQAMultiplicity", x);\r
-    Fill("MC_histEventStats", 4.); // analyzed MC events\r
-    Fill("MC_histVertexZ", pAOD->GetPrimaryVertex()->GetZ());\r
-    CalculateMoments("MC_", tracksMC, tracksMC, zVertex, 1.);\r
-    delete tracksMC;\r
-  }\r
-}\r
-\r
-void AliAnalysisTaskLongRangeCorrelations::Terminate(Option_t* ) {\r
-  //\r
-  fOutputList = dynamic_cast<TList*>(GetOutputData(1));\r
-  if (NULL == fOutputList) {\r
-    AliFatal("NULL == fOutputList");\r
-    return; // needed to avoid coverity warning\r
-  }\r
-}\r
-\r
-TString AliAnalysisTaskLongRangeCorrelations::GetOutputListName() const {\r
-  TString listName("listLRC");\r
-  listName += TString::Format("_%smix",         fRunMixing ? "" : "no");\r
-  listName += TString::Format("_trackFilt%d",   fTrackFilter);\r
-  listName += TString::Format("_cent%.0fT%.0f", fCentMin, fCentMax);\r
-  listName += TString::Format("_ptMin%.0fMeV",  1e3*fPtMin);\r
-  listName += TString::Format("_phi%.0fT%.0f",  TMath::RadToDeg()*fPhiMin, TMath::RadToDeg()*fPhiMax);\r
-  if ( 1 == fSelectPrimaryMCParticles)\r
-    listName += "_selPrimMC";\r
-  if (-1 == fSelectPrimaryMCParticles)\r
-    listName += "_selNonPrimMC";\r
-  if ( 1 == fSelectPrimaryMCDataParticles)\r
-    listName += "_selPrimMCData";\r
-  if (-1 == fSelectPrimaryMCDataParticles)\r
-    listName += "_selNonPrimMCData";\r
-  if (-1 != fNMin)\r
-    listName += TString::Format("_nMin%d", fNMin);\r
-  if (-1 != fNMax)\r
-    listName += TString::Format("_nMax%d", fNMax);\r
-  return listName;\r
-}\r
-\r
-void AliAnalysisTaskLongRangeCorrelations::SetupForMixing() {\r
-  const Int_t trackDepth(fMixingTracks);\r
-  const Int_t poolsize(1000); // Maximum number of events\r
-\r
-  Double_t centralityBins[] = { fCentMin, fCentMax };\r
-  const Int_t nCentralityBins(sizeof(centralityBins)/sizeof(Double_t) - 1);\r
-\r
-  fPoolMgr = new AliEventPoolManager(poolsize, trackDepth,\r
-                                    nCentralityBins, centralityBins,\r
-                                    fVertexZ->GetNbins()+1, \r
-                                    const_cast<Double_t*>(fVertexZ->GetXbins()->GetArray()));\r
-}\r
-\r
-THnSparse* AliAnalysisTaskLongRangeCorrelations::MakeHistSparsePhiEta(const char* name) const {\r
-  const Int_t nVertexZBins=fVertexZ->GetNbins();\r
-  const Int_t   nBinsM1[] = { fnBinsPhi, fnBinsEta, nVertexZBins     };\r
-  const Double_t xMinM1[] = {  fxMinPhi,  fxMinEta, 0.5              };\r
-  const Double_t xMaxM1[] = {  fxMaxPhi,  fxMaxEta, nVertexZBins+0.5 };\r
-  const TString title(TString(name)\r
-                     +";#phi;#eta;vertex Z (bin);");\r
-  return new THnSparseD(name, title.Data(), 3, nBinsM1, xMinM1, xMaxM1);\r
-}\r
-\r
-AliTHn* AliAnalysisTaskLongRangeCorrelations::MakeHistPhiEta(const char* name) const {\r
-  const Int_t nVertexZBins=fVertexZ->GetNbins();\r
-  const Int_t   nBinsM1[] = { fnBinsPhi, fnBinsEta, nVertexZBins     };\r
-  const Double_t xMinM1[] = {  fxMinPhi,  fxMinEta, 0.5              };\r
-  const Double_t xMaxM1[] = {  fxMaxPhi,  fxMaxEta, nVertexZBins+0.5 };\r
-  const TString title(TString(name)\r
-                     +";#phi;#eta;vertex Z (bin);");\r
-  AliTHn* h(new AliTHn(name, title.Data(), 1, 3, nBinsM1));\r
-  for (Int_t i=0; i<3; ++i)\r
-    h->SetBinLimits(i, xMinM1[i], xMaxM1[i]);\r
-  return h;\r
-}\r
-\r
-AliTHn* AliAnalysisTaskLongRangeCorrelations::MakeHistPhiEtaPhiEta(const char* name) const {\r
-  const Int_t nVertexZBins=fVertexZ->GetNbins();\r
-  const Int_t   nBinsM2[] = {  fnBinsPhi, fnBinsEta,  fnBinsPhi, fnBinsEta, nVertexZBins     };\r
-  const Double_t xMinM2[] = {  fxMinPhi,  fxMinEta,   fxMinPhi,  fxMinEta,  0.5              };\r
-  const Double_t xMaxM2[] = {  fxMaxPhi,  fxMaxEta,   fxMaxPhi,  fxMaxEta,  nVertexZBins+0.5 };\r
-  const TString title(TString(name)\r
-                     +";#phi_{1};#eta_{1}"\r
-                     +";#phi_{2};#eta_{2};vertex Z (bin);");\r
-  AliTHn* h(new AliTHn(name, title.Data(), 1, 5, nBinsM2));\r
-  for (Int_t i=0; i<5; ++i)\r
-    h->SetBinLimits(i, xMinM2[i], xMaxM2[i]);\r
-  return h;\r
-}\r
-\r
-TObjArray* AliAnalysisTaskLongRangeCorrelations::GetAcceptedTracks(AliAODEvent*  pAOD,\r
-                                                                  TClonesArray* arrayMC,\r
-                                                                  Double_t      centrality) {\r
-  TObjArray* tracks= new TObjArray;\r
-  tracks->SetOwner(kTRUE);\r
-\r
-  const Long64_t N(pAOD->GetNumberOfTracks());\r
-  AliDebug(5, Form("#tracks= %6lld %f", N, centrality));\r
-  for (Long64_t i(0); i<N; ++i) {\r
-    AliAODTrack* pAODTrack(dynamic_cast<AliAODTrack*>(pAOD->GetTrack(i)));\r
-    if (NULL == pAODTrack) continue;\r
-\r
-    // track filter selection\r
-    if (!pAODTrack->TestFilterBit(fTrackFilter)) continue;\r
-\r
-    // select only primary tracks\r
-    if (pAODTrack->GetType() != AliAODTrack::kPrimary) continue;\r
-\r
-    if (NULL != arrayMC) {\r
-      const Int_t label(pAODTrack->GetLabel());\r
-      AliAODMCParticle* mcParticle((label >= 0) \r
-                                  ? static_cast<AliAODMCParticle*>(arrayMC->At(label))\r
-                                  : NULL);\r
-      if (label >=0 && NULL == mcParticle)\r
-       AliFatal("MC particle not found");\r
-\r
-//       const Bool_t isPhysicalPrimary(mcParticle->IsPhysicalPrimary());\r
-//       Printf("mcParticle->IsPhysicalPrimary() %d", isPhysicalPrimary);\r
-\r
-      switch (fSelectPrimaryMCDataParticles) {\r
-      case -1:\r
-       if (label < 0) continue;\r
-       if (kTRUE  == mcParticle->IsPhysicalPrimary()) continue;\r
-       break;\r
-      case  0:\r
-       // NOP, take all tracks\r
-       break;\r
-      case  1:\r
-       if (label < 0) continue;\r
-       if (kFALSE == mcParticle->IsPhysicalPrimary()) continue;\r
-       break;\r
-      default:\r
-       AliFatal("fSelectPrimaryMCDataParticles != {-1,0,1}");\r
-      }            \r
-    }\r
-\r
-    // select only charged tracks\r
-    if (pAODTrack->Charge() == 0) continue;\r
-    \r
-    Fill("histQACentPt", pAODTrack->Charge()>0, centrality,       pAODTrack->Pt());\r
-    Fill("histQAPhiEta", pAODTrack->Charge()>0, pAODTrack->Phi(), pAODTrack->Eta());\r
-    if (pAODTrack->Phi() < fPhiMin || pAODTrack->Phi() > fPhiMax) continue;\r
-    if (pAODTrack->Pt()  < fPtMin  || pAODTrack->Pt()  > fPtMax)  continue;\r
-\r
-    tracks->Add(new LRCParticle(pAODTrack->Eta(), pAODTrack->Phi()));\r
-  } // next track\r
-  return tracks;\r
-}\r
-\r
-TObjArray* AliAnalysisTaskLongRangeCorrelations::GetAcceptedTracks(TClonesArray* tracksMC,\r
-                                                                  Double_t centrality) {\r
-  // for keeping track of MC labels\r
-  std::set<Int_t> labelSet;\r
-\r
-  TObjArray* tracks= new TObjArray;\r
-  tracks->SetOwner(kTRUE);\r
-\r
-  const Long64_t N(tracksMC->GetEntriesFast());\r
-  AliDebug(5, Form("#tracks= %6lld %f", N, centrality));\r
-  for (Long64_t i(0); i<N; ++i) {\r
-    AliAODMCParticle* pMCTrack(dynamic_cast<AliAODMCParticle*>(tracksMC->At(i)));\r
-    if (NULL == pMCTrack) continue;    \r
-\r
-    // no track filter selection for MC tracks    \r
-\r
-    if (labelSet.find(pMCTrack->Label()) != labelSet.end()) {\r
-      Printf("Duplicate Label= %3d", pMCTrack->Label());\r
-      continue;\r
-    }\r
-    labelSet.insert(pMCTrack->Label());    \r
-    \r
-//     Printf("isPrim = %d", pMCTrack->IsPhysicalPrimary());\r
-\r
-    switch (fSelectPrimaryMCParticles) {\r
-    case -1:\r
-      if (kTRUE == pMCTrack->IsPhysicalPrimary()) continue;\r
-      break;\r
-    case  0:\r
-      // NOP, take all MC tracks\r
-      break;\r
-    case  1:\r
-      if (kFALSE == pMCTrack->IsPhysicalPrimary()) continue;\r
-      break;\r
-    default:\r
-      AliFatal("fSelectPrimaryMCParticles != {-1,0,1}");\r
-    }\r
-\r
-    // select only charged tracks\r
-    if (pMCTrack->Charge() == 0) continue;\r
-\r
-    Fill("MC_histQACentPt", pMCTrack->Charge()>0, centrality,      pMCTrack->Pt());\r
-    Fill("MC_histQAPhiEta", pMCTrack->Charge()>0, pMCTrack->Phi(), pMCTrack->Eta());\r
-\r
-    if (pMCTrack->Phi() < fPhiMin || pMCTrack->Phi() > fPhiMax) continue;\r
-    if (pMCTrack->Pt()  < fPtMin  || pMCTrack->Pt()  > fPtMax)  continue;\r
-\r
-    tracks->Add(new LRCParticle(pMCTrack->Eta(), pMCTrack->Phi()));\r
-  } // next track\r
-  return tracks;\r
-}\r
-\r
-Bool_t AddTHnSparseToAliTHn(AliTHn* h, THnSparse* hs, Double_t weight) {\r
-  if (h->GetNVar() != hs->GetNdimensions())\r
-    return kFALSE;\r
-\r
-  const size_t nDim(hs->GetNdimensions());\r
-  \r
-  Int_t    *coord = new Int_t[nDim];\r
-  Double_t *x     = new Double_t[nDim];\r
-\r
-  const Long64_t N(hs->GetNbins());\r
-  for (Long64_t i(0); i<N; ++i) {\r
-    const Double_t n(hs->GetBinContent(i, coord));\r
-    for (size_t j=0; j<nDim; ++j) \r
-      x[j] = hs->GetAxis(j)->GetBinCenter(coord[j]);\r
-    h->Fill(x, 0, n*weight);\r
-  }\r
-\r
-  delete[] coord;\r
-  delete[] x;\r
-\r
-  return kTRUE;\r
-}\r
-\r
-void AliAnalysisTaskLongRangeCorrelations::CalculateMoments(TString prefix,\r
-                                                           TObjArray* tracks1,\r
-                                                           TObjArray* tracks2,\r
-                                                           Double_t vertexZ,\r
-                                                           Double_t weight) {\r
-  const Long64_t nc1(tracks1->GetEntriesFast());\r
-  if (fNMin != -1 && nc1 < fNMin) return;\r
-  if (fNMax != -1 && nc1 > fNMax) return;\r
-\r
-  const Long64_t nc2(tracks1->GetEntriesFast());\r
-  if (fNMin != -1 && nc2 < fNMin) return;\r
-  if (fNMax != -1 && nc2 > fNMax) return;\r
-\r
-  AliTHn* hN1(dynamic_cast<AliTHn*>(fOutputList->FindObject(prefix+"histMoment1PhiEta_1")));\r
-  if (NULL == hN1) return;\r
-\r
-  // <n_1>\r
-  THnSparse* hN1ForThisEvent(ComputeNForThisEvent(tracks1, "hN1", vertexZ));\r
-  AddTHnSparseToAliTHn(hN1,hN1ForThisEvent, weight); \r
-//   hN1->GetGrid(0)->GetGrid()->Add(hN1ForThisEvent, weight);\r
-\r
-  // n(eta) distributions\r
-  FillNEtaHist(prefix+"histNEta_300", hN1ForThisEvent, weight);\r
-  FillNEtaHist(prefix+"histNEta_120", hN1ForThisEvent, weight);\r
-  FillNEtaHist(prefix+"histNEta__30", hN1ForThisEvent, weight);\r
-  FillNEtaHist(prefix+"histNEta__15", hN1ForThisEvent, weight);\r
-\r
-  TObjArray* hNs(new TObjArray);\r
-\r
-  // <n_1 n_1>\r
-  hNs->AddAt(hN1ForThisEvent, 0);\r
-  hNs->AddAt(hN1ForThisEvent, 1);\r
-\r
-  ComputeNXForThisEvent(hNs, prefix+"histMoment2PhiEtaPhiEta_11", vertexZ, weight);\r
-\r
-  if (fRunMixing) {\r
-    AliTHn* hN2(dynamic_cast<AliTHn*>(fOutputList->FindObject(prefix+"histMoment1PhiEta_2")));\r
-    if (NULL == hN2) return;\r
-    // <n_2>\r
-    THnSparse* hN2ForThisEvent(ComputeNForThisEvent(tracks2, "hN2", vertexZ));\r
-    AddTHnSparseToAliTHn(hN2,hN2ForThisEvent, weight); \r
-//     hN2->GetGrid(0)->GetGrid()->Add(hN2ForThisEvent, weight);\r
-\r
-    // <n_1 n_2>\r
-    hNs->AddAt(hN2ForThisEvent, 1);\r
-    ComputeNXForThisEvent(hNs, prefix+"histMoment2PhiEtaPhiEta_12", vertexZ, weight);\r
-    \r
-    // <n_2 n_2>\r
-    hNs->AddAt(hN2ForThisEvent, 0);\r
-    ComputeNXForThisEvent(hNs, prefix+"histMoment2PhiEtaPhiEta_22", vertexZ, weight);\r
-\r
-    // clean up\r
-    delete hN2ForThisEvent;\r
-  }\r
-\r
-  // clean up\r
-  delete hNs;\r
-  delete hN1ForThisEvent;\r
-}\r
-\r
-void AliAnalysisTaskLongRangeCorrelations::FillNEtaHist(TString name,\r
-                                                       THnSparse* hs,\r
-                                                       Double_t weight) {\r
-\r
-  TH2* hSum(dynamic_cast<TH2*>(fOutputList->FindObject(name)));\r
-  if (NULL == hSum) return;\r
-\r
-  TH2* hPerEvent(dynamic_cast<TH2*>(hSum->Clone("hPerEvent")));\r
-  if (NULL == hPerEvent) return;\r
-  hPerEvent->Reset();  \r
-\r
-  TH1* h1PerEvent(hPerEvent->ProjectionX());\r
-\r
-  // fill h1PerEvent\r
-  const Long64_t N(hs->GetNbins());\r
-  for (Long64_t i(0); i<N; ++i) {\r
-    Int_t coord[2] = { 0, 0 };\r
-    const Double_t n(hs->GetBinContent(i, coord));\r
-    const Double_t eta(hs->GetAxis(1)->GetBinCenter(coord[1]));\r
-    h1PerEvent->Fill(eta, n);\r
-  }\r
\r
-  for (Int_t i(1); i<=h1PerEvent->GetNbinsX(); ++i)\r
-    hPerEvent->Fill(h1PerEvent->GetXaxis()->GetBinCenter(i),\r
-                   h1PerEvent->GetBinContent(i));\r
-\r
-  hSum->Add(hPerEvent, weight);\r
-\r
-  delete hPerEvent;\r
-  delete h1PerEvent;\r
-}\r
-\r
-THnSparse* AliAnalysisTaskLongRangeCorrelations::ComputeNForThisEvent(TObjArray* tracks,\r
-                                                                     const char* histName,\r
-                                                                     Double_t vertexZ) const {\r
-  const Double_t vertexZBin(fVertexZ->FindBin(vertexZ));\r
-  THnSparse* hN(MakeHistSparsePhiEta(histName));\r
-  const Long64_t nTracks(tracks->GetEntriesFast());\r
-  for (Long64_t i(0); i<nTracks; ++i) {\r
-    const LRCParticle* p(dynamic_cast<LRCParticle*>(tracks->At(i)));\r
-    if (NULL == p) continue;\r
-    const Double_t x[] = { p->Phi(), p->Eta(), vertexZBin };\r
-    hN->Fill(x);\r
-  }\r
-  return hN;\r
-}\r
-\r
-class MultiDimIterator {\r
-public:\r
-  MultiDimIterator(TObjArray* _fHs, Double_t vertexZBin)\r
-    : fHs(_fHs)\r
-    , fN(fHs->GetEntriesFast())\r
-    , fDims(fN,  0)\r
-    , fIdxs(fN,  0)\r
-    , fNs  (fN,  0)\r
-    , fX (2*fN+1, 0)\r
-    , fj(0) {\r
-    for (Long64_t i=0; i<fN; ++i) {\r
-      THnSparse* hNi(reinterpret_cast<THnSparse*>(fHs->At(i)));\r
-      if (NULL == hNi)\r
-       AliFatal("NULL == hNi");\r
-      fDims[i] = hNi->GetNbins();\r
-      AliDebug(3, Form("%lld %s %lld", i, hNi->GetName(), fDims[i]));\r
-      update(i);\r
-    }\r
-    fX[2*fN] = vertexZBin;\r
-  }\r
-\r
-  const char* ClassName() const { return "MultiDimIterator"; }\r
-  const Double_t* GetX() const { return &fX.front(); }\r
-  Double_t        GetN() const { // returns the product of all multiplicities\r
-    return std::accumulate(fNs.begin(), fNs.end(), Double_t(1), std::multiplies<Double_t>());\r
-  }\r
-  Bool_t end() const { return fj == fN; }\r
-\r
-  MultiDimIterator& operator++() {\r
-    Long64_t j(0);\r
-    for (; j<fN; ++j) {\r
-      ++fIdxs[j];\r
-      update(j);\r
-      if (fIdxs[j] < fDims[j]) break;\r
-      fIdxs[j] = 0;\r
-      update(j);\r
-    }\r
-    fj = j;\r
-    return *this;\r
-  }\r
-\r
-protected:\r
-  void update(Long64_t j) {\r
-    THnSparse* hs(reinterpret_cast<THnSparse*>(fHs->At(j)));\r
-    Int_t coord[] = { 0, 0 };\r
-    fNs[j] = hs->GetBinContent(fIdxs[j], coord);\r
-    for (Int_t k(0); k<2; ++k)\r
-      fX[2*j+k] = hs->GetAxis(k)->GetBinCenter(coord[k]);\r
-  }\r
-private:\r
-  MultiDimIterator(const MultiDimIterator&);\r
-  MultiDimIterator& operator=(const MultiDimIterator&);\r
-\r
-  TObjArray*            fHs;   // array of THnSparse histograms\r
-  const Long64_t        fN;    // number of histograms\r
-  std::vector<Long64_t> fDims; // number of filled bins for each THnSparse\r
-  std::vector<Long64_t> fIdxs; // indices\r
-  std::vector<Double_t> fNs;   // number of tracks\r
-  std::vector<Double_t> fX;    // coordinate\r
-  Long64_t              fj;    // state\r
-} ;\r
-\r
-void AliAnalysisTaskLongRangeCorrelations::ComputeNXForThisEvent(TObjArray* hNs,\r
-                                                                const char* histName,\r
-                                                                Double_t vertexZ,\r
-                                                                Double_t weight) {\r
-  if (NULL == fOutputList) return;\r
-\r
-  AliTHn* hs(dynamic_cast<AliTHn*>(fOutputList->FindObject(histName)));\r
-  if (hs == NULL) return;\r
-\r
-  for (MultiDimIterator mdi(hNs, fVertexZ->FindBin(vertexZ)); !mdi.end(); ++mdi)\r
-    hs->Fill(mdi.GetX(), 0, mdi.GetN()*weight);\r
-\r
-//   hs->FillParent();\r
-}\r
-\r
-void AliAnalysisTaskLongRangeCorrelations::Fill(const char* histName, Double_t x) {\r
-  if (NULL == fOutputList) return;\r
-  TH1* h = dynamic_cast<TH1*>(fOutputList->FindObject(histName));\r
-  if (h == NULL) return;\r
-  h->Fill(x);\r
-}\r
-void AliAnalysisTaskLongRangeCorrelations::Fill(const char* histName, Double_t x, Double_t y) {\r
-  if (NULL == fOutputList) return;\r
-  TH2* h = dynamic_cast<TH2*>(fOutputList->FindObject(histName));\r
-  if (h == NULL) return;\r
-  h->Fill(x, y);\r
-}\r
-void AliAnalysisTaskLongRangeCorrelations::Fill(const char* histName, Double_t x, Double_t y, Double_t z) {\r
-  if (NULL == fOutputList) return;\r
-  TH3* h = dynamic_cast<TH3*>(fOutputList->FindObject(histName));\r
-  if (h == NULL) return;\r
-  h->Fill(x, y, z);\r
-}\r
-void AliAnalysisTaskLongRangeCorrelations::Fill(const char* histName, const Double_t* x, Double_t w) {\r
-  if (NULL == fOutputList) return;\r
-  THnSparse* h = dynamic_cast<THnSparse*>(fOutputList->FindObject(histName));\r
-  if (h == NULL) return;\r
-  h->Fill(x, w);\r
-}\r
+// -*- C++ -*-
+/**************************************************************************
+ * Copyright(c) 2012,      ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+// $Id: AliAnalysisTaskLongRangeCorrelations.cxx 359 2013-10-29 15:39:09Z cmayer $
+
+#include <numeric>
+#include <functional>
+#include <set>
+
+#include <TChain.h>
+#include <THashList.h>
+#include <THnSparse.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TH3.h>
+
+#include "AliEventPoolManager.h"
+#include "AliAnalysisManager.h"
+#include "AliAODEvent.h"
+#include "AliAODMCHeader.h"
+#include "AliAODMCParticle.h"
+#include "AliMCEvent.h"
+#include "AliInputEventHandler.h"
+#include "AliLog.h"
+#include "AliTHn.h"
+
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+#include "AliAnalysisTaskLongRangeCorrelations.h"
+
+
+ClassImp(AliAnalysisTaskLongRangeCorrelations);
+ClassImp(LRCParticle);
+
+AliAnalysisTaskLongRangeCorrelations::AliAnalysisTaskLongRangeCorrelations(const char *name)
+  : AliAnalysisTaskSE(name)
+  , fOutputList(NULL)
+  , fVertexZ(NULL)
+  , fRunMixing(kFALSE)
+  , fPoolMgr(NULL)
+  , fMixingTracks(50000)
+  , fTrackFilter(128)
+  , fCentMin(0), fCentMax(20)
+  , fPtMin(0.2), fPtMax(1e10)
+  , fPhiMin(0.), fPhiMax(TMath::TwoPi())
+  , fMaxAbsVertexZ(10.)
+  , fSelectPrimaryMCParticles(0)
+  , fSelectPrimaryMCDataParticles(0)
+  , fNMin(-1)
+  , fNMax(-1)
+  , fnBinsCent( 220), fnBinsPt(400), fnBinsPhi(4),              fnBinsEta(120)
+  , fxMinCent( -5.0), fxMinPt( 0.0), fxMinPhi( 0.0),            fxMinEta(-1.5)
+  , fxMaxCent(105.0), fxMaxPt( 4.0), fxMaxPhi( TMath::TwoPi()), fxMaxEta( 1.5) {
+  DefineInput(0, TChain::Class());
+  DefineOutput(1, TList::Class());
+}
+
+AliAnalysisTaskLongRangeCorrelations::~AliAnalysisTaskLongRangeCorrelations() {
+  if (NULL != fOutputList) {
+    delete fOutputList;
+    fOutputList = NULL;
+  }
+}
+
+void AliAnalysisTaskLongRangeCorrelations::UserCreateOutputObjects() {
+  fOutputList = new THashList;
+  fOutputList->SetOwner(kTRUE);
+  fOutputList->SetName(GetOutputListName());
+
+  const Double_t vertexBins[] = {  // Zvtx bins
+    -10., -7., -5., -3., -1., 1., 3., 5., 7., 10.
+  };
+  const Int_t nVertexBins(sizeof(vertexBins)/sizeof(Double_t)-1);
+  fVertexZ = new TAxis(nVertexBins, vertexBins);
+  fVertexZ->SetName("VertexZAxis");
+  fOutputList->Add(fVertexZ);
+
+  fOutputList->Add(new TH1D("histVertexZ", ";vertex Z (cm)", nVertexBins, vertexBins));
+
+  // Event statistics
+  const char* eventStatLabels[] = { 
+    "All Events",
+    "Physics Selection",
+    "Centrality Selection",
+    "Vertex Selection",
+    "Analyzed Events"
+  };
+  const size_t nEventStat(sizeof(eventStatLabels)/sizeof(const char*));
+  TH1* hStats(new TH1D("histEventStats", "histEventStats", nEventStat, -0.5, nEventStat-0.5));
+  for (size_t i(0); i<nEventStat; ++i)
+    hStats->GetXaxis()->SetBinLabel(1+i, eventStatLabels[i]);
+  fOutputList->Add(hStats);
+
+  // QA histograms
+  fOutputList->Add(new TH2D("histQACentrality", ";centrality V0M(%);selected;",
+                           fnBinsCent,fxMinCent,fxMaxCent, 2, -0.5, 1.5));
+  fOutputList->Add(new TH2D("histQAVertexZ", ";vertex-z (cm);selected;",
+                           800, -40., 40., 2, -0.5, 1.5));
+  fOutputList->Add(new TH1D("histQAMultiplicityBeforeCuts",
+                           ";tracks", 10000, 0, 10000));
+  fOutputList->Add(new TH1D("histQAMultiplicityAfterCuts",
+                           ";selected tracks", 10000, 0, 10000));
+  fOutputList->Add(new TH3D("histQACentPt", ";charge;centrality V0M(%);p_{T} (GeV/c)",
+                           2, -0.5, 1.5, fnBinsCent, fxMinCent, fxMaxCent, fnBinsPt, fxMinPt, fxMaxPt));
+  fOutputList->Add(new TH3D("histQAPhiEta", ";charge;#phi (rad);#eta",
+                           2, -0.5, 1.5, 200, 0.0, TMath::TwoPi(), 300, -1.5, 1.5));
+
+  // N(eta) distributions with different binnings
+  fOutputList->Add(new TH2D("histNEta_120", ";#eta;N", 120, -1.5, 1.5, 1000, -.5, 1000-.5));  // 0.025
+  fOutputList->Add(new TH2D("histNEta__60", ";#eta;N",  60, -1.5, 1.5, 1000, -.5, 1000-.5));  // 0.05
+  fOutputList->Add(new TH2D("histNEta__30", ";#eta;N",  30, -1.5, 1.5, 1000, -.5, 1000-.5));  // 0.1
+  fOutputList->Add(new TH2D("histNEta__15", ";#eta;N",  15, -1.5, 1.5, 1000, -.5, 1000-.5));  // 0.2
+
+  // Moments
+  fOutputList->Add(MakeHistPhiEta("histMoment1PhiEta_1"));
+  if (fRunMixing)
+    fOutputList->Add(MakeHistPhiEta("histMoment1PhiEta_2"));
+  fOutputList->Add(MakeHistPhiEtaPhiEta("histMoment2PhiEtaPhiEta_11"));
+  if (fRunMixing) {
+    fOutputList->Add(MakeHistPhiEtaPhiEta("histMoment2PhiEtaPhiEta_12"));
+    fOutputList->Add(MakeHistPhiEtaPhiEta("histMoment2PhiEtaPhiEta_22"));
+  }
+  
+  // add MC Histograms
+  const Int_t N(fOutputList->GetEntries());
+  for (Int_t i(0); i<N; ++i)
+    fOutputList->Add(fOutputList->At(i)->Clone(TString("MC_") + fOutputList->At(i)->GetName()));
+
+  if (fRunMixing)
+    SetupForMixing();
+
+  PostData(1, fOutputList);
+}
+
+void AliAnalysisTaskLongRangeCorrelations::UserExec(Option_t* ) {
+  AliAnalysisManager* pManager(AliAnalysisManager::GetAnalysisManager());
+  if (NULL == pManager) return;
+
+  AliInputEventHandler* pInputHandler(dynamic_cast<AliInputEventHandler*>(pManager->GetInputEventHandler()));
+  if (NULL == pInputHandler) return;
+
+  AliAODEvent* pAOD(dynamic_cast<AliAODEvent*>(InputEvent()));
+  if (NULL == pAOD) return;
+
+  AliAODHeader *pAODHeader = pAOD->GetHeader();
+  if (NULL == pAODHeader) return;
+
+  AliAODMCHeader *pAODMCHeader(dynamic_cast<AliAODMCHeader*>(pAOD->FindListObject(AliAODMCHeader::StdBranchName())));
+  const Bool_t isMC(NULL != pAODMCHeader);
+
+  TClonesArray *arrayMC(NULL);
+  if (isMC) {
+    arrayMC = dynamic_cast<TClonesArray*>(pAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
+    if (NULL == arrayMC) return;
+  }
+
+  // --- event cuts MC ---
+  if (isMC) {
+    Fill("MC_histEventStats", 0.); // all events
+    Fill("MC_histEventStats", 1.); // events passing physics selection
+  }
+
+  // --- event cuts data ---
+  Fill("histEventStats", 0.); // all events
+
+  // physics selection
+  if (!pInputHandler->IsEventSelected()) return;
+  Fill("histEventStats", 1.); // events passing physics selection
+
+  // centrality selection
+  const Double_t centrality(pAODHeader->GetCentralityP()->GetCentralityPercentile("V0M"));
+  AliDebug(3, Form("centrality=%f", centrality));
+  const Bool_t centralitySelected(centrality > fCentMin && centrality < fCentMax);
+  Fill("histQACentrality", centrality, centralitySelected);
+  if (!centralitySelected) return;
+  Fill("histEventStats", 2.); // events passing centrality selection
+  if (isMC)
+    Fill("MC_histEventStats", 2.); // events passing centrality selection
+
+  // vertex selection -- data
+  const Int_t nVertex(pAOD->GetNumberOfVertices());
+  if (0 == nVertex) return;
+  const AliAODVertex* pVertex(pAOD->GetPrimaryVertex());
+  if (NULL == pVertex) return;
+  const Int_t nTracksPrimary(pVertex->GetNContributors());
+  if (nTracksPrimary < 1) return;
+
+  const Double_t zVertex(pVertex->GetZ());
+  const Bool_t vertexSelectedData(TMath::Abs(zVertex) < fMaxAbsVertexZ);
+
+  // vertex selection -- MC
+  Bool_t vertexSelectedMC(kTRUE);
+  if (isMC)
+    vertexSelectedMC = (TMath::Abs(pAODMCHeader->GetVtxZ()) < fMaxAbsVertexZ);
+
+  // combined vertex selection (data and MC)
+  const Bool_t vertexSelected(vertexSelectedData && vertexSelectedMC);
+  Fill("histQAVertexZ", zVertex, vertexSelected);
+  if (isMC)
+    Fill("MC_histQAVertexZ", pAODMCHeader->GetVtxZ(), vertexSelected);
+  if (!vertexSelected) return;
+  Fill("histEventStats",    3.); // events passing vertex selection
+  if (isMC)
+    Fill("MC_histEventStats", 3.); // events passing vertex selection
+
+  // ------------------
+  // event is accepted
+  // ------------------
+
+  TObjArray* tracksMain(GetAcceptedTracks(pAOD, arrayMC, centrality));
+  Fill("histQAMultiplicityBeforeCuts", pAOD->GetNumberOfTracks());
+  Fill("histQAMultiplicityAfterCuts", tracksMain->GetEntriesFast());
+  
+  if (fRunMixing) {  
+    AliEventPool* pEventPool(fPoolMgr->GetEventPool(centrality, pAOD->GetPrimaryVertex()->GetZ()));
+    if (NULL == pEventPool)
+      AliFatal(Form("No pool found for centrality = %f, zVtx = %f", centrality, pAOD->GetPrimaryVertex()->GetZ()));
+    
+//     pEventPool->PrintInfo();
+    if (pEventPool->IsReady()
+       || pEventPool->NTracksInPool()     > fMixingTracks/10
+       || pEventPool->GetCurrentNEvents() >= 5) {
+      Fill("histEventStats", 4.); // analyzed events
+      Fill("histVertexZ", zVertex);
+      const Int_t nMix(pEventPool->GetCurrentNEvents());
+      for (Int_t i(0); i<nMix; ++i) {
+       TObjArray* tracksMixed(pEventPool->GetEvent(i));
+       CalculateMoments("", tracksMain, tracksMixed, zVertex, 1./nMix);
+      }
+    }
+    // Update the Event pool
+    pEventPool->UpdatePool(tracksMain);
+  } else { // no mixing
+    Fill("histEventStats", 4.); // analyzed events
+    Fill("histVertexZ", zVertex);
+    CalculateMoments("", tracksMain, tracksMain, zVertex, 1.);
+    delete tracksMain;
+  }
+  
+  if (isMC) {
+    TObjArray* tracksMC(GetAcceptedTracks(arrayMC, centrality));
+    const Double_t x[2] = {
+      arrayMC->GetEntriesFast(),
+      tracksMC->GetEntriesFast()
+    };
+    Fill("MC_histQAMultiplicity", x);
+    Fill("MC_histEventStats", 4.); // analyzed MC events
+    Fill("MC_histVertexZ", pAOD->GetPrimaryVertex()->GetZ());
+    CalculateMoments("MC_", tracksMC, tracksMC, zVertex, 1.);
+    delete tracksMC;
+  }
+}
+
+void AliAnalysisTaskLongRangeCorrelations::Terminate(Option_t* ) {
+  //
+  fOutputList = dynamic_cast<TList*>(GetOutputData(1));
+  if (NULL == fOutputList) {
+    AliFatal("NULL == fOutputList");
+    return; // needed to avoid coverity warning
+  }
+}
+
+TString AliAnalysisTaskLongRangeCorrelations::GetOutputListName() const {
+  TString listName("listLRC");
+  listName += TString::Format("_%smix",         fRunMixing ? "" : "no");
+  listName += TString::Format("_trackFilt%d",   fTrackFilter);
+  listName += TString::Format("_cent%.0fT%.0f", fCentMin, fCentMax);
+  listName += TString::Format("_ptMin%.0fMeV",  1e3*fPtMin);
+  listName += TString::Format("_phi%.0fT%.0f",  TMath::RadToDeg()*fPhiMin, TMath::RadToDeg()*fPhiMax);
+  if ( 1 == fSelectPrimaryMCParticles)
+    listName += "_selPrimMC";
+  if (-1 == fSelectPrimaryMCParticles)
+    listName += "_selNonPrimMC";
+  if ( 1 == fSelectPrimaryMCDataParticles)
+    listName += "_selPrimMCData";
+  if (-1 == fSelectPrimaryMCDataParticles)
+    listName += "_selNonPrimMCData";
+  if (-1 != fNMin)
+    listName += TString::Format("_nMin%d", fNMin);
+  if (-1 != fNMax)
+    listName += TString::Format("_nMax%d", fNMax);
+  return listName;
+}
+
+void AliAnalysisTaskLongRangeCorrelations::SetupForMixing() {
+  const Int_t trackDepth(fMixingTracks);
+  const Int_t poolsize(1000); // Maximum number of events
+
+  Double_t centralityBins[] = { fCentMin, fCentMax };
+  const Int_t nCentralityBins(sizeof(centralityBins)/sizeof(Double_t) - 1);
+
+  fPoolMgr = new AliEventPoolManager(poolsize, trackDepth,
+                                    nCentralityBins, centralityBins,
+                                    fVertexZ->GetNbins()+1, 
+                                    const_cast<Double_t*>(fVertexZ->GetXbins()->GetArray()));
+}
+
+THnSparse* AliAnalysisTaskLongRangeCorrelations::MakeHistSparsePhiEta(const char* name) const {
+  const Int_t nVertexZBins=fVertexZ->GetNbins();
+  const Int_t   nBinsM1[] = { fnBinsPhi, fnBinsEta, nVertexZBins     };
+  const Double_t xMinM1[] = {  fxMinPhi,  fxMinEta, 0.5              };
+  const Double_t xMaxM1[] = {  fxMaxPhi,  fxMaxEta, nVertexZBins+0.5 };
+  const TString title(TString(name)
+                     +";#phi;#eta;vertex Z (bin);");
+  return new THnSparseD(name, title.Data(), 3, nBinsM1, xMinM1, xMaxM1);
+}
+
+AliTHn* AliAnalysisTaskLongRangeCorrelations::MakeHistPhiEta(const char* name) const {
+  const Int_t nVertexZBins=fVertexZ->GetNbins();
+  const Int_t   nBinsM1[] = { fnBinsPhi, fnBinsEta, nVertexZBins     };
+  const Double_t xMinM1[] = {  fxMinPhi,  fxMinEta, 0.5              };
+  const Double_t xMaxM1[] = {  fxMaxPhi,  fxMaxEta, nVertexZBins+0.5 };
+  const TString title(TString(name)
+                     +";#phi;#eta;vertex Z (bin);");
+  AliTHn* h(new AliTHn(name, title.Data(), 1, 3, nBinsM1));
+  for (Int_t i=0; i<3; ++i)
+    h->SetBinLimits(i, xMinM1[i], xMaxM1[i]);
+  return h;
+}
+
+AliTHn* AliAnalysisTaskLongRangeCorrelations::MakeHistPhiEtaPhiEta(const char* name) const {
+  const Int_t nVertexZBins=fVertexZ->GetNbins();
+  const Int_t   nBinsM2[] = {  fnBinsPhi, fnBinsEta,  fnBinsPhi, fnBinsEta, nVertexZBins     };
+  const Double_t xMinM2[] = {  fxMinPhi,  fxMinEta,   fxMinPhi,  fxMinEta,  0.5              };
+  const Double_t xMaxM2[] = {  fxMaxPhi,  fxMaxEta,   fxMaxPhi,  fxMaxEta,  nVertexZBins+0.5 };
+  const TString title(TString(name)
+                     +";#phi_{1};#eta_{1}"
+                     +";#phi_{2};#eta_{2};vertex Z (bin);");
+  AliTHn* h(new AliTHn(name, title.Data(), 1, 5, nBinsM2));
+  for (Int_t i=0; i<5; ++i)
+    h->SetBinLimits(i, xMinM2[i], xMaxM2[i]);
+  return h;
+}
+
+TObjArray* AliAnalysisTaskLongRangeCorrelations::GetAcceptedTracks(AliAODEvent*  pAOD,
+                                                                  TClonesArray* arrayMC,
+                                                                  Double_t      centrality) {
+  TObjArray* tracks= new TObjArray;
+  tracks->SetOwner(kTRUE);
+
+  const Long64_t N(pAOD->GetNumberOfTracks());
+  AliDebug(5, Form("#tracks= %6lld %f", N, centrality));
+  for (Long64_t i(0); i<N; ++i) {
+    AliAODTrack* pAODTrack(dynamic_cast<AliAODTrack*>(pAOD->GetTrack(i)));
+    if (NULL == pAODTrack) continue;
+
+    // track filter selection
+    if (!pAODTrack->TestFilterBit(fTrackFilter)) continue;
+
+    // select only primary tracks
+    if (pAODTrack->GetType() != AliAODTrack::kPrimary) continue;
+
+    if (NULL != arrayMC) {
+      const Int_t label(pAODTrack->GetLabel());
+      AliAODMCParticle* mcParticle((label >= 0) 
+                                  ? static_cast<AliAODMCParticle*>(arrayMC->At(label))
+                                  : NULL);
+      if (label >=0 && NULL == mcParticle)
+       AliFatal("MC particle not found");
+
+//       const Bool_t isPhysicalPrimary(mcParticle->IsPhysicalPrimary());
+//       Printf("mcParticle->IsPhysicalPrimary() %d", isPhysicalPrimary);
+
+      switch (fSelectPrimaryMCDataParticles) {
+      case -1:
+       if (label < 0) continue;
+       if (kTRUE  == mcParticle->IsPhysicalPrimary()) continue;
+       break;
+      case  0:
+       // NOP, take all tracks
+       break;
+      case  1:
+       if (label < 0) continue;
+       if (kFALSE == mcParticle->IsPhysicalPrimary()) continue;
+       break;
+      default:
+       AliFatal("fSelectPrimaryMCDataParticles != {-1,0,1}");
+      }            
+    }
+
+    // select only charged tracks
+    if (pAODTrack->Charge() == 0) continue;
+    
+    Fill("histQACentPt", pAODTrack->Charge()>0, centrality,       pAODTrack->Pt());
+    Fill("histQAPhiEta", pAODTrack->Charge()>0, pAODTrack->Phi(), pAODTrack->Eta());
+    if (pAODTrack->Phi() < fPhiMin || pAODTrack->Phi() > fPhiMax) continue;
+    if (pAODTrack->Pt()  < fPtMin  || pAODTrack->Pt()  > fPtMax)  continue;
+
+    tracks->Add(new LRCParticle(pAODTrack->Eta(), pAODTrack->Phi()));
+  } // next track
+  return tracks;
+}
+
+TObjArray* AliAnalysisTaskLongRangeCorrelations::GetAcceptedTracks(TClonesArray* tracksMC,
+                                                                  Double_t centrality) {
+  // for keeping track of MC labels
+  std::set<Int_t> labelSet;
+
+  TObjArray* tracks= new TObjArray;
+  tracks->SetOwner(kTRUE);
+
+  const Long64_t N(tracksMC->GetEntriesFast());
+  AliDebug(5, Form("#tracks= %6lld %f", N, centrality));
+  for (Long64_t i(0); i<N; ++i) {
+    AliAODMCParticle* pMCTrack(dynamic_cast<AliAODMCParticle*>(tracksMC->At(i)));
+    if (NULL == pMCTrack) continue;    
+
+    // no track filter selection for MC tracks    
+
+    if (labelSet.find(pMCTrack->Label()) != labelSet.end()) {
+      Printf("Duplicate Label= %3d", pMCTrack->Label());
+      continue;
+    }
+    labelSet.insert(pMCTrack->Label());    
+    
+//     Printf("isPrim = %d", pMCTrack->IsPhysicalPrimary());
+
+    switch (fSelectPrimaryMCParticles) {
+    case -1:
+      if (kTRUE == pMCTrack->IsPhysicalPrimary()) continue;
+      break;
+    case  0:
+      // NOP, take all MC tracks
+      break;
+    case  1:
+      if (kFALSE == pMCTrack->IsPhysicalPrimary()) continue;
+      break;
+    default:
+      AliFatal("fSelectPrimaryMCParticles != {-1,0,1}");
+    }
+
+    // select only charged tracks
+    if (pMCTrack->Charge() == 0) continue;
+
+    Fill("MC_histQACentPt", pMCTrack->Charge()>0, centrality,      pMCTrack->Pt());
+    Fill("MC_histQAPhiEta", pMCTrack->Charge()>0, pMCTrack->Phi(), pMCTrack->Eta());
+
+    if (pMCTrack->Phi() < fPhiMin || pMCTrack->Phi() > fPhiMax) continue;
+    if (pMCTrack->Pt()  < fPtMin  || pMCTrack->Pt()  > fPtMax)  continue;
+
+    tracks->Add(new LRCParticle(pMCTrack->Eta(), pMCTrack->Phi()));
+  } // next track
+  return tracks;
+}
+
+Bool_t AddTHnSparseToAliTHn(AliTHn* h, THnSparse* hs, Double_t weight) {
+  if (h->GetNVar() != hs->GetNdimensions())
+    return kFALSE;
+
+  const size_t nDim(hs->GetNdimensions());
+  
+  Int_t    *coord = new Int_t[nDim];
+  Double_t *x     = new Double_t[nDim];
+
+  const Long64_t N(hs->GetNbins());
+  for (Long64_t i(0); i<N; ++i) {
+    const Double_t n(hs->GetBinContent(i, coord));
+    for (size_t j=0; j<nDim; ++j) 
+      x[j] = hs->GetAxis(j)->GetBinCenter(coord[j]);
+    h->Fill(x, 0, n*weight);
+  }
+
+  delete[] coord;
+  delete[] x;
+
+  return kTRUE;
+}
+
+void AliAnalysisTaskLongRangeCorrelations::CalculateMoments(TString prefix,
+                                                           TObjArray* tracks1,
+                                                           TObjArray* tracks2,
+                                                           Double_t vertexZ,
+                                                           Double_t weight) {
+  const Long64_t nc1(tracks1->GetEntriesFast());
+  if (fNMin != -1 && nc1 < fNMin) return;
+  if (fNMax != -1 && nc1 > fNMax) return;
+
+  const Long64_t nc2(tracks1->GetEntriesFast());
+  if (fNMin != -1 && nc2 < fNMin) return;
+  if (fNMax != -1 && nc2 > fNMax) return;
+
+  AliTHn* hN1(dynamic_cast<AliTHn*>(fOutputList->FindObject(prefix+"histMoment1PhiEta_1")));
+  if (NULL == hN1) return;
+
+  // <n_1>
+  THnSparse* hN1ForThisEvent(ComputeNForThisEvent(tracks1, "hN1", vertexZ));
+  AddTHnSparseToAliTHn(hN1,hN1ForThisEvent, weight); 
+//   hN1->GetGrid(0)->GetGrid()->Add(hN1ForThisEvent, weight);
+
+  // n(eta) distributions
+  FillNEtaHist(prefix+"histNEta_300", hN1ForThisEvent, weight);
+  FillNEtaHist(prefix+"histNEta_120", hN1ForThisEvent, weight);
+  FillNEtaHist(prefix+"histNEta__30", hN1ForThisEvent, weight);
+  FillNEtaHist(prefix+"histNEta__15", hN1ForThisEvent, weight);
+
+  TObjArray* hNs(new TObjArray);
+
+  // <n_1 n_1>
+  hNs->AddAt(hN1ForThisEvent, 0);
+  hNs->AddAt(hN1ForThisEvent, 1);
+
+  ComputeNXForThisEvent(hNs, prefix+"histMoment2PhiEtaPhiEta_11", vertexZ, weight);
+
+  if (fRunMixing) {
+    AliTHn* hN2(dynamic_cast<AliTHn*>(fOutputList->FindObject(prefix+"histMoment1PhiEta_2")));
+    if (NULL == hN2) return;
+    // <n_2>
+    THnSparse* hN2ForThisEvent(ComputeNForThisEvent(tracks2, "hN2", vertexZ));
+    AddTHnSparseToAliTHn(hN2,hN2ForThisEvent, weight); 
+//     hN2->GetGrid(0)->GetGrid()->Add(hN2ForThisEvent, weight);
+
+    // <n_1 n_2>
+    hNs->AddAt(hN2ForThisEvent, 1);
+    ComputeNXForThisEvent(hNs, prefix+"histMoment2PhiEtaPhiEta_12", vertexZ, weight);
+    
+    // <n_2 n_2>
+    hNs->AddAt(hN2ForThisEvent, 0);
+    ComputeNXForThisEvent(hNs, prefix+"histMoment2PhiEtaPhiEta_22", vertexZ, weight);
+
+    // clean up
+    delete hN2ForThisEvent;
+  }
+
+  // clean up
+  delete hNs;
+  delete hN1ForThisEvent;
+}
+
+void AliAnalysisTaskLongRangeCorrelations::FillNEtaHist(TString name,
+                                                       THnSparse* hs,
+                                                       Double_t weight) {
+
+  TH2* hSum(dynamic_cast<TH2*>(fOutputList->FindObject(name)));
+  if (NULL == hSum) return;
+
+  TH2* hPerEvent(dynamic_cast<TH2*>(hSum->Clone("hPerEvent")));
+  if (NULL == hPerEvent) return;
+  hPerEvent->Reset();  
+
+  TH1* h1PerEvent(hPerEvent->ProjectionX());
+
+  // fill h1PerEvent
+  const Long64_t N(hs->GetNbins());
+  for (Long64_t i(0); i<N; ++i) {
+    Int_t coord[2] = { 0, 0 };
+    const Double_t n(hs->GetBinContent(i, coord));
+    const Double_t eta(hs->GetAxis(1)->GetBinCenter(coord[1]));
+    h1PerEvent->Fill(eta, n);
+  }
+  for (Int_t i(1); i<=h1PerEvent->GetNbinsX(); ++i)
+    hPerEvent->Fill(h1PerEvent->GetXaxis()->GetBinCenter(i),
+                   h1PerEvent->GetBinContent(i));
+
+  hSum->Add(hPerEvent, weight);
+
+  delete hPerEvent;
+  delete h1PerEvent;
+}
+
+THnSparse* AliAnalysisTaskLongRangeCorrelations::ComputeNForThisEvent(TObjArray* tracks,
+                                                                     const char* histName,
+                                                                     Double_t vertexZ) const {
+  const Double_t vertexZBin(fVertexZ->FindBin(vertexZ));
+  THnSparse* hN(MakeHistSparsePhiEta(histName));
+  const Long64_t nTracks(tracks->GetEntriesFast());
+  for (Long64_t i(0); i<nTracks; ++i) {
+    const LRCParticle* p(dynamic_cast<LRCParticle*>(tracks->At(i)));
+    if (NULL == p) continue;
+    const Double_t x[] = { p->Phi(), p->Eta(), vertexZBin };
+    hN->Fill(x);
+  }
+  return hN;
+}
+
+class MultiDimIterator {
+public:
+  MultiDimIterator(TObjArray* _fHs, Double_t vertexZBin)
+    : fHs(_fHs)
+    , fN(fHs->GetEntriesFast())
+    , fDims(fN,  0)
+    , fIdxs(fN,  0)
+    , fNs  (fN,  0)
+    , fX (2*fN+1, 0)
+    , fj(0) {
+    for (Long64_t i=0; i<fN; ++i) {
+      THnSparse* hNi(reinterpret_cast<THnSparse*>(fHs->At(i)));
+      if (NULL == hNi)
+       AliFatal("NULL == hNi");
+      fDims[i] = hNi->GetNbins();
+      AliDebug(3, Form("%lld %s %lld", i, hNi->GetName(), fDims[i]));
+      update(i);
+    }
+    fX[2*fN] = vertexZBin;
+  }
+
+  const char* ClassName() const { return "MultiDimIterator"; }
+  const Double_t* GetX() const { return &fX.front(); }
+  Double_t        GetN() const { // returns the product of all multiplicities
+    return std::accumulate(fNs.begin(), fNs.end(), Double_t(1), std::multiplies<Double_t>());
+  }
+  Bool_t end() const { return fj == fN; }
+
+  MultiDimIterator& operator++() {
+    Long64_t j(0);
+    for (; j<fN; ++j) {
+      ++fIdxs[j];
+      update(j);
+      if (fIdxs[j] < fDims[j]) break;
+      fIdxs[j] = 0;
+      update(j);
+    }
+    fj = j;
+    return *this;
+  }
+
+protected:
+  void update(Long64_t j) {
+    THnSparse* hs(reinterpret_cast<THnSparse*>(fHs->At(j)));
+    Int_t coord[] = { 0, 0 };
+    fNs[j] = hs->GetBinContent(fIdxs[j], coord);
+    for (Int_t k(0); k<2; ++k)
+      fX[2*j+k] = hs->GetAxis(k)->GetBinCenter(coord[k]);
+  }
+private:
+  MultiDimIterator(const MultiDimIterator&);
+  MultiDimIterator& operator=(const MultiDimIterator&);
+
+  TObjArray*            fHs;   // array of THnSparse histograms
+  const Long64_t        fN;    // number of histograms
+  std::vector<Long64_t> fDims; // number of filled bins for each THnSparse
+  std::vector<Long64_t> fIdxs; // indices
+  std::vector<Double_t> fNs;   // number of tracks
+  std::vector<Double_t> fX;    // coordinate
+  Long64_t              fj;    // state
+} ;
+
+void AliAnalysisTaskLongRangeCorrelations::ComputeNXForThisEvent(TObjArray* hNs,
+                                                                const char* histName,
+                                                                Double_t vertexZ,
+                                                                Double_t weight) {
+  if (NULL == fOutputList) return;
+
+  AliTHn* hs(dynamic_cast<AliTHn*>(fOutputList->FindObject(histName)));
+  if (hs == NULL) return;
+
+  for (MultiDimIterator mdi(hNs, fVertexZ->FindBin(vertexZ)); !mdi.end(); ++mdi)
+    hs->Fill(mdi.GetX(), 0, mdi.GetN()*weight);
+
+//   hs->FillParent();
+}
+
+void AliAnalysisTaskLongRangeCorrelations::Fill(const char* histName, Double_t x) {
+  if (NULL == fOutputList) return;
+  TH1* h = dynamic_cast<TH1*>(fOutputList->FindObject(histName));
+  if (h == NULL) return;
+  h->Fill(x);
+}
+void AliAnalysisTaskLongRangeCorrelations::Fill(const char* histName, Double_t x, Double_t y) {
+  if (NULL == fOutputList) return;
+  TH2* h = dynamic_cast<TH2*>(fOutputList->FindObject(histName));
+  if (h == NULL) return;
+  h->Fill(x, y);
+}
+void AliAnalysisTaskLongRangeCorrelations::Fill(const char* histName, Double_t x, Double_t y, Double_t z) {
+  if (NULL == fOutputList) return;
+  TH3* h = dynamic_cast<TH3*>(fOutputList->FindObject(histName));
+  if (h == NULL) return;
+  h->Fill(x, y, z);
+}
+void AliAnalysisTaskLongRangeCorrelations::Fill(const char* histName, const Double_t* x, Double_t w) {
+  if (NULL == fOutputList) return;
+  THnSparse* h = dynamic_cast<THnSparse*>(fOutputList->FindObject(histName));
+  if (h == NULL) return;
+  h->Fill(x, w);
+}
index c86e772..1920a84 100644 (file)
-// -*- C++ -*-\r
-// $Id: AliAnalysisTaskLongRangeCorrelations.h 341 2013-09-30 15:59:19Z cmayer $\r
-#ifndef _AliAnalysisTaskLongRangeCorrelations_H_\r
-#define _AliAnalysisTaskLongRangeCorrelations_H_\r
-\r
-/* Copyright(c) 2012, ALICE Experiment at CERN, All rights reserved.      *\r
- * See cxx source for full Copyright notice                               */\r
-\r
-////////////////////////////////////////////////////////////////////////\r
-//\r
-// Analysis class for Long-range Correlations\r
-//\r
-// Look for correlations in eta (and in phi)\r
-//\r
-// This class needs input AODs.\r
-// The output is a list of analysis-specific containers.\r
-//\r
-//    Author:\r
-//    Christoph Mayer\r
-// \r
-////////////////////////////////////////////////////////////////////////\r
-\r
-\r
-class TAxis;\r
-class TClonesArray;\r
-class TList;\r
-class TObjArray;\r
-\r
-class AliAODEvent;\r
-class AliAODHeader;\r
-class AliEventPoolManager;\r
-class AliTHn;\r
-\r
-#include <TObject.h>\r
-#include "AliAnalysisTaskSE.h"\r
-\r
-class AliAnalysisTaskLongRangeCorrelations : public AliAnalysisTaskSE { \r
-public: \r
-  AliAnalysisTaskLongRangeCorrelations(const char *name="AliAnalysisTaskLongRangeCorrelations");\r
-  virtual ~AliAnalysisTaskLongRangeCorrelations();\r
-\r
-  /*   virtual void NotifyRun(); */\r
-  virtual void UserCreateOutputObjects();\r
-  virtual void UserExec(Option_t *);\r
-  virtual void Terminate(Option_t *);\r
-\r
-  void SetRunMixing(Bool_t runMixing)      { fRunMixing    = runMixing; }\r
-  void SetMixingTracks(Int_t mixingTracks) { fMixingTracks = mixingTracks; }\r
-  void SetTrackFilter(Int_t trackFilter)   { fTrackFilter  = trackFilter; }\r
-\r
-  void SetCentralityRange(Double_t centMin, Double_t centMax) {\r
-    fCentMin = centMin; fCentMax = centMax;\r
-  }\r
-  void SetPtRange(Double_t ptMin, Double_t ptMax) {\r
-    fPtMin = ptMin; fPtMax = ptMax;\r
-  }\r
-  void SetPhiRange(Double_t phiMin, Double_t phiMax) {\r
-    fPhiMin = phiMin; fPhiMax = phiMax;\r
-  }\r
-  void SetMaxAbsVertexZ(Double_t maxAbsVertexZ) { fMaxAbsVertexZ = maxAbsVertexZ; }\r
-\r
-  void SetSelectPrimaryMCParticles(Int_t flagMC, Int_t flagMCData) {\r
-    fSelectPrimaryMCParticles     = flagMC;\r
-    fSelectPrimaryMCDataParticles = flagMCData;\r
-  }\r
-\r
-  void SetRangeN(Int_t nMin, Int_t nMax) {\r
-    fNMin = nMin;\r
-    fNMax = nMax;\r
-  }\r
-\r
-  TString GetOutputListName() const;\r
-\r
-protected:\r
-  // for now up to second moments:\r
-  // <n_1>(phi_1,eta_1)\r
-  // <n_2>(phi_2,eta_2) \r
-  // <n_1, n_2>(phi_1,eta_1;phi_2,eta_2)\r
-  void       CalculateMoments(TString, TObjArray*, TObjArray*, Double_t, Double_t); \r
-  void       ComputeNXForThisEvent(TObjArray*, const char*, Double_t, Double_t);\r
-  THnSparse* ComputeNForThisEvent(TObjArray*, const char*, Double_t) const;\r
-  void       FillNEtaHist(TString, THnSparse*, Double_t);\r
-\r
-  TObjArray* GetAcceptedTracks(AliAODEvent*, TClonesArray*, Double_t);\r
-  TObjArray* GetAcceptedTracks(TClonesArray*, Double_t);\r
-\r
-  // filling histograms by name\r
-  void Fill(const char*, Double_t);                           // TH1 weight=1\r
-  void Fill(const char*, Double_t, Double_t );                // TH2 weight=1\r
-  void Fill(const char*, Double_t, Double_t, Double_t);       // TH3 weight=1\r
-  void Fill(const char*, const Double_t*, Double_t weight=1); // THnSparse\r
-\r
-  void SetupForMixing();\r
-\r
-  THnSparse* MakeHistSparsePhiEta(const char* name) const;\r
-  AliTHn* MakeHistPhiEta(const char* name) const;\r
-  AliTHn* MakeHistPhiEtaPhiEta(const char* name) const; \r
-\r
-private:\r
-  TList*               fOutputList;         //! Output list\r
-  TAxis*               fVertexZ;            //! vertex z bins\r
-  Bool_t               fRunMixing;          //\r
-  AliEventPoolManager* fPoolMgr;            //! event pool manager  \r
-  Int_t                fMixingTracks;       //\r
-  Int_t                fTrackFilter;        //\r
-  Double_t             fCentMin, fCentMax;  // centrality range\r
-  Double_t             fPtMin, fPtMax;      // P_{T} range\r
-  Double_t             fPhiMin, fPhiMax;    // #phi range\r
-  Double_t             fMaxAbsVertexZ;      // max abs(zvertex)\r
-  Int_t                fSelectPrimaryMCParticles;     // 0: no, 1: yes, -1: only non-primary particles\r
-  Int_t                fSelectPrimaryMCDataParticles; // 0: no, 1: yes, -1: only non-primary particles\r
-  Int_t                fNMin;\r
-  Int_t                fNMax;\r
-  // histogram data\r
-  Int_t    fnBinsCent, fnBinsPt, fnBinsPhi, fnBinsEta;\r
-  Double_t fxMinCent,  fxMinPt,  fxMinPhi,  fxMinEta;\r
-  Double_t fxMaxCent,  fxMaxPt,  fxMaxPhi,  fxMaxEta;  \r
-\r
-  AliAnalysisTaskLongRangeCorrelations(const AliAnalysisTaskLongRangeCorrelations&); // not implemented\r
-  AliAnalysisTaskLongRangeCorrelations& operator=(const AliAnalysisTaskLongRangeCorrelations&); // not implemented\r
-  \r
-  ClassDef(AliAnalysisTaskLongRangeCorrelations, 1);\r
-} ; \r
-\r
-class LRCParticle : public TObject {\r
-public:\r
-  LRCParticle(Double_t eta=0, Double_t phi=0)\r
-    : fEta(eta), fPhi(phi) {}\r
-  virtual ~LRCParticle() {}\r
-\r
-  Double_t Eta() const { return fEta; }\r
-  Double_t Phi() const { return fPhi; }\r
-\r
-protected:\r
-private:\r
-  LRCParticle(const LRCParticle&);\r
-  LRCParticle& operator=(const LRCParticle&);\r
-\r
-  Double_t fEta;\r
-  Double_t fPhi;\r
-  ClassDef(LRCParticle, 1);\r
-} ;\r
-#endif // _AliAnalysisTaskLongRangeCorrelations_H_\r
-\r
+// -*- C++ -*-
+// $Id: AliAnalysisTaskLongRangeCorrelations.h 341 2013-09-30 15:59:19Z cmayer $
+#ifndef _AliAnalysisTaskLongRangeCorrelations_H_
+#define _AliAnalysisTaskLongRangeCorrelations_H_
+
+/* Copyright(c) 2012, ALICE Experiment at CERN, All rights reserved.      *
+ * See cxx source for full Copyright notice                               */
+
+////////////////////////////////////////////////////////////////////////
+//
+// Analysis class for Long-range Correlations
+//
+// Look for correlations in eta (and in phi)
+//
+// This class needs input AODs.
+// The output is a list of analysis-specific containers.
+//
+//    Author:
+//    Christoph Mayer
+// 
+////////////////////////////////////////////////////////////////////////
+
+
+class TAxis;
+class TClonesArray;
+class TList;
+class TObjArray;
+
+class AliAODEvent;
+class AliAODHeader;
+class AliEventPoolManager;
+class AliTHn;
+
+#include <TObject.h>
+#include "AliAnalysisTaskSE.h"
+
+class AliAnalysisTaskLongRangeCorrelations : public AliAnalysisTaskSE { 
+public: 
+  AliAnalysisTaskLongRangeCorrelations(const char *name="AliAnalysisTaskLongRangeCorrelations");
+  virtual ~AliAnalysisTaskLongRangeCorrelations();
+
+  /*   virtual void NotifyRun(); */
+  virtual void UserCreateOutputObjects();
+  virtual void UserExec(Option_t *);
+  virtual void Terminate(Option_t *);
+
+  void SetRunMixing(Bool_t runMixing)      { fRunMixing    = runMixing; }
+  void SetMixingTracks(Int_t mixingTracks) { fMixingTracks = mixingTracks; }
+  void SetTrackFilter(Int_t trackFilter)   { fTrackFilter  = trackFilter; }
+
+  void SetCentralityRange(Double_t centMin, Double_t centMax) {
+    fCentMin = centMin; fCentMax = centMax;
+  }
+  void SetPtRange(Double_t ptMin, Double_t ptMax) {
+    fPtMin = ptMin; fPtMax = ptMax;
+  }
+  void SetPhiRange(Double_t phiMin, Double_t phiMax) {
+    fPhiMin = phiMin; fPhiMax = phiMax;
+  }
+  void SetMaxAbsVertexZ(Double_t maxAbsVertexZ) { fMaxAbsVertexZ = maxAbsVertexZ; }
+
+  void SetSelectPrimaryMCParticles(Int_t flagMC, Int_t flagMCData) {
+    fSelectPrimaryMCParticles     = flagMC;
+    fSelectPrimaryMCDataParticles = flagMCData;
+  }
+
+  void SetRangeN(Int_t nMin, Int_t nMax) {
+    fNMin = nMin;
+    fNMax = nMax;
+  }
+
+  TString GetOutputListName() const;
+
+protected:
+  // for now up to second moments:
+  // <n_1>(phi_1,eta_1)
+  // <n_2>(phi_2,eta_2) 
+  // <n_1, n_2>(phi_1,eta_1;phi_2,eta_2)
+  void       CalculateMoments(TString, TObjArray*, TObjArray*, Double_t, Double_t); 
+  void       ComputeNXForThisEvent(TObjArray*, const char*, Double_t, Double_t);
+  THnSparse* ComputeNForThisEvent(TObjArray*, const char*, Double_t) const;
+  void       FillNEtaHist(TString, THnSparse*, Double_t);
+
+  TObjArray* GetAcceptedTracks(AliAODEvent*, TClonesArray*, Double_t);
+  TObjArray* GetAcceptedTracks(TClonesArray*, Double_t);
+
+  // filling histograms by name
+  void Fill(const char*, Double_t);                           // TH1 weight=1
+  void Fill(const char*, Double_t, Double_t );                // TH2 weight=1
+  void Fill(const char*, Double_t, Double_t, Double_t);       // TH3 weight=1
+  void Fill(const char*, const Double_t*, Double_t weight=1); // THnSparse
+
+  void SetupForMixing();
+
+  THnSparse* MakeHistSparsePhiEta(const char* name) const;
+  AliTHn* MakeHistPhiEta(const char* name) const;
+  AliTHn* MakeHistPhiEtaPhiEta(const char* name) const; 
+
+private:
+  TList*               fOutputList;         //! Output list
+  TAxis*               fVertexZ;            //! vertex z bins
+  Bool_t               fRunMixing;          //
+  AliEventPoolManager* fPoolMgr;            //! event pool manager  
+  Int_t                fMixingTracks;       //
+  Int_t                fTrackFilter;        //
+  Double_t             fCentMin, fCentMax;  // centrality range
+  Double_t             fPtMin, fPtMax;      // P_{T} range
+  Double_t             fPhiMin, fPhiMax;    // #phi range
+  Double_t             fMaxAbsVertexZ;      // max abs(zvertex)
+  Int_t                fSelectPrimaryMCParticles;     // 0: no, 1: yes, -1: only non-primary particles
+  Int_t                fSelectPrimaryMCDataParticles; // 0: no, 1: yes, -1: only non-primary particles
+  Int_t                fNMin;
+  Int_t                fNMax;
+  // histogram data
+  Int_t    fnBinsCent, fnBinsPt, fnBinsPhi, fnBinsEta;
+  Double_t fxMinCent,  fxMinPt,  fxMinPhi,  fxMinEta;
+  Double_t fxMaxCent,  fxMaxPt,  fxMaxPhi,  fxMaxEta;  
+
+  AliAnalysisTaskLongRangeCorrelations(const AliAnalysisTaskLongRangeCorrelations&); // not implemented
+  AliAnalysisTaskLongRangeCorrelations& operator=(const AliAnalysisTaskLongRangeCorrelations&); // not implemented
+  
+  ClassDef(AliAnalysisTaskLongRangeCorrelations, 1);
+} ; 
+
+class LRCParticle : public TObject {
+public:
+  LRCParticle(Double_t eta=0, Double_t phi=0)
+    : fEta(eta), fPhi(phi) {}
+  virtual ~LRCParticle() {}
+
+  Double_t Eta() const { return fEta; }
+  Double_t Phi() const { return fPhi; }
+
+protected:
+private:
+  LRCParticle(const LRCParticle&);
+  LRCParticle& operator=(const LRCParticle&);
+
+  Double_t fEta;
+  Double_t fPhi;
+  ClassDef(LRCParticle, 1);
+} ;
+#endif // _AliAnalysisTaskLongRangeCorrelations_H_
+
index 35034a2..8b0c17d 100644 (file)
@@ -1,58 +1,58 @@
-// -*- c++ -*-\r
-// $Id: AddTaskLongRangeCorrelations.C 341 2013-09-30 15:59:19Z cmayer $\r
-\r
-const Double_t centMin[] = {  0,   0,  10,  20,  30,  40,  50,  60,  70,  80 };\r
-const Double_t centMax[] = {  5,  10,  20,  30,  40,  50,  60,  70,  80, 100 };\r
-const Int_t    nMin[]    = { -1, 120,  70,  35,  20,   5,   0,   0,   0,  -1 };\r
-const Int_t    nMax[]    = { -1, 400, 275, 200, 145, 100,  68,  46,  30,  -1 };\r
-\r
-AliAnalysisTaskLongRangeCorrelations*\r
-AddTaskLongRangeCorrelations(Int_t  trackFilter  = 128, // TPC only\r
-                            Bool_t runMixing    = !kTRUE,\r
-                            Int_t  mixingTracks = 50000,\r
-                            Int_t selPrimMC  = 0, Int_t selPrimMCData = 0,\r
-                            Double_t ptMin   = 0.2, \r
-                            Double_t phiMin  = 0, Double_t phiMax  = TMath::TwoPi()) {\r
-\r
-  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
-  if (NULL == mgr) {\r
-    ::Error("AddTaskLongRangeCorrelations", "No analysis manager to connect to.");\r
-    return NULL;\r
-  }\r
-  if (NULL == mgr->GetInputEventHandler()) {\r
-    ::Error("AddTaskLongRangeCorrelations", "This task requires an input event handler");\r
-    return NULL;\r
-  }\r
-  TString type = mgr->GetInputEventHandler()->GetDataType();\r
-  if (type != "AOD") {\r
-    ::Error("AddTaskLongRangeCorrelations", "This task runs only on AOD data");\r
-    return NULL;\r
-  }\r
-\r
-  AliAnalysisTaskLongRangeCorrelations *taskLRC = NULL;\r
-  AliAnalysisDataContainer             *listLRC = NULL;\r
-  TString outputFileName = AliAnalysisManager::GetCommonFileName();\r
-  outputFileName += ":PWGCFEbyE.outputLongRangeCorrelations.root";\r
-\r
-  for (Int_t i=0; i<sizeof(centMin)/sizeof(Double_t); ++i) {\r
-    taskLRC = new AliAnalysisTaskLongRangeCorrelations("TaskLongRangeCorrelations");\r
-    taskLRC->SetRunMixing(runMixing);\r
-    taskLRC->SetMixingTracks(mixingTracks);\r
-    taskLRC->SetTrackFilter(trackFilter);\r
-    taskLRC->SetCentralityRange(centMin[i], centMax[i]);\r
-    taskLRC->SetPtRange(ptMin, 1e20);\r
-    taskLRC->SetPhiRange(phiMin, phiMax);\r
-    taskLRC->SelectCollisionCandidates(AliVEvent::kMB);\r
-    taskLRC->SetSelectPrimaryMCParticles(selPrimMC, selPrimMCData);\r
-    taskLRC->SetRangeN(nMin[i], nMax[i]);\r
-\r
-    listLRC = mgr->CreateContainer(taskLRC->GetOutputListName(), TList::Class(),\r
-                                  AliAnalysisManager::kOutputContainer,\r
-                                  outputFileName.Data());\r
-    mgr->AddTask(taskLRC);\r
-    mgr->ConnectInput(taskLRC,  0, mgr->GetCommonInputContainer());\r
-    mgr->ConnectOutput(taskLRC, 1, listLRC);\r
-  }\r
-\r
-  return taskLRC;\r
-}\r
+// -*- c++ -*-
+// $Id: AddTaskLongRangeCorrelations.C 341 2013-09-30 15:59:19Z cmayer $
+
+const Double_t centMin[] = {  0,   0,  10,  20,  30,  40,  50,  60,  70,  80 };
+const Double_t centMax[] = {  5,  10,  20,  30,  40,  50,  60,  70,  80, 100 };
+const Int_t    nMin[]    = { -1, 120,  70,  35,  20,   5,   0,   0,   0,  -1 };
+const Int_t    nMax[]    = { -1, 400, 275, 200, 145, 100,  68,  46,  30,  -1 };
+
+AliAnalysisTaskLongRangeCorrelations*
+AddTaskLongRangeCorrelations(Int_t  trackFilter  = 128, // TPC only
+                            Bool_t runMixing    = !kTRUE,
+                            Int_t  mixingTracks = 50000,
+                            Int_t selPrimMC  = 0, Int_t selPrimMCData = 0,
+                            Double_t ptMin   = 0.2, 
+                            Double_t phiMin  = 0, Double_t phiMax  = TMath::TwoPi()) {
+
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (NULL == mgr) {
+    ::Error("AddTaskLongRangeCorrelations", "No analysis manager to connect to.");
+    return NULL;
+  }
+  if (NULL == mgr->GetInputEventHandler()) {
+    ::Error("AddTaskLongRangeCorrelations", "This task requires an input event handler");
+    return NULL;
+  }
+  TString type = mgr->GetInputEventHandler()->GetDataType();
+  if (type != "AOD") {
+    ::Error("AddTaskLongRangeCorrelations", "This task runs only on AOD data");
+    return NULL;
+  }
+
+  AliAnalysisTaskLongRangeCorrelations *taskLRC = NULL;
+  AliAnalysisDataContainer             *listLRC = NULL;
+  TString outputFileName = AliAnalysisManager::GetCommonFileName();
+  outputFileName += ":PWGCFEbyE.outputLongRangeCorrelations.root";
+
+  for (Int_t i=0; i<sizeof(centMin)/sizeof(Double_t); ++i) {
+    taskLRC = new AliAnalysisTaskLongRangeCorrelations("TaskLongRangeCorrelations");
+    taskLRC->SetRunMixing(runMixing);
+    taskLRC->SetMixingTracks(mixingTracks);
+    taskLRC->SetTrackFilter(trackFilter);
+    taskLRC->SetCentralityRange(centMin[i], centMax[i]);
+    taskLRC->SetPtRange(ptMin, 1e20);
+    taskLRC->SetPhiRange(phiMin, phiMax);
+    taskLRC->SelectCollisionCandidates(AliVEvent::kMB);
+    taskLRC->SetSelectPrimaryMCParticles(selPrimMC, selPrimMCData);
+    taskLRC->SetRangeN(nMin[i], nMax[i]);
+
+    listLRC = mgr->CreateContainer(taskLRC->GetOutputListName(), TList::Class(),
+                                  AliAnalysisManager::kOutputContainer,
+                                  outputFileName.Data());
+    mgr->AddTask(taskLRC);
+    mgr->ConnectInput(taskLRC,  0, mgr->GetCommonInputContainer());
+    mgr->ConnectOutput(taskLRC, 1, listLRC);
+  }
+
+  return taskLRC;
+}