PWGCF/Correlations/DPhi/AliAnalysisTaskLongRangeCorrelations.cxx

   1 // -*- C++ -*-
   2 /**************************************************************************
   3  * Copyright(c) 2012,      ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16 // $Id: AliAnalysisTaskLongRangeCorrelations.cxx 283 2013-03-26 19:07:58Z cmayer $
  17
  18 #include <numeric>
  19 #include <functional>
  20 #include <set>
  21
  22 #include <TChain.h>
  23 #include <THashList.h>
  24 #include <THnSparse.h>
  25 #include <TH1.h>
  26 #include <TH2.h>
  27 #include <TH3.h>
  28
  29 #include "AliEventPoolManager.h"
  30 #include "AliAnalysisManager.h"
  31 #include "AliAODEvent.h"
  32 #include "AliAODMCHeader.h"
  33 #include "AliAODMCParticle.h"
  34 #include "AliMCEvent.h"
  35 #include "AliInputEventHandler.h"
  36 #include "AliLog.h"
  37
  38 #include "AliAnalysisTaskLongRangeCorrelations.h"
  39
  40
  41 ClassImp(AliAnalysisTaskLongRangeCorrelations);
  42 ClassImp(LRCParticle);
  43
  44 AliAnalysisTaskLongRangeCorrelations::AliAnalysisTaskLongRangeCorrelations(const char *name)
  45   : AliAnalysisTaskSE(name)
  46   , fOutputList(NULL)
  47   , fVertexZ(NULL)
  48   , fRunMixing(kFALSE)
  49   , fPoolMgr(NULL)
  50   , fMixingTracks(50000)
  51   , fTrackFilter(128)
  52   , fCentMin(0), fCentMax(20)
  53   , fPtMin(0.2), fPtMax(1e10)
  54   , fPhiMin(0.), fPhiMax(TMath::TwoPi())
  55   , fMaxAbsVertexZ(10.)
  56   , fnBinsCent( 220), fnBinsPt(400), fnBinsPhi(4),              fnBinsEta(120)
  57   , fxMinCent( -5.0), fxMinPt( 0.0), fxMinPhi( 0.0),            fxMinEta(-1.5)
  58   , fxMaxCent(105.0), fxMaxPt( 4.0), fxMaxPhi( TMath::TwoPi()), fxMaxEta( 1.5) {
  59   DefineInput(0, TChain::Class());
  60   DefineOutput(1, TList::Class());
  61 }
  62
  63 AliAnalysisTaskLongRangeCorrelations::~AliAnalysisTaskLongRangeCorrelations() {
  64   if (NULL != fOutputList) {
  65     delete fOutputList;
  66     fOutputList = NULL;
  67   }
  68 }
  69
  70 void AliAnalysisTaskLongRangeCorrelations::UserCreateOutputObjects() {
  71   fOutputList = new THashList;
  72   fOutputList->SetOwner(kTRUE);
  73   fOutputList->SetName(GetOutputListName());
  74
  75   const Double_t vertexBins[] = {  // Zvtx bins
  76     -10., -7., -5., -3., -1., 1., 3., 5., 7., 10.
  77   };
  78   fVertexZ = new TAxis(sizeof(vertexBins)/sizeof(Double_t)-1, vertexBins);
  79   fVertexZ->SetName("VertexZAxis");
  80   fOutputList->Add(fVertexZ);
  81
  82   // Event statistics
  83   const char* eventStatLabels[] = {
  84     "All Events",
  85     "Physics Selection",
  86     "Centrality Selection",
  87     "Vertex Selection",
  88     "Analyzed Events"
  89   };
  90   const size_t nEventStat(sizeof(eventStatLabels)/sizeof(const char*));
  91   TH1* hStats(new TH1D("histEventStats", "histEventStats", nEventStat, -0.5, nEventStat-0.5));
  92   for (size_t i(0); i<nEventStat; ++i)
  93     hStats->GetXaxis()->SetBinLabel(1+i, eventStatLabels[i]);
  94   fOutputList->Add(hStats);
  95
  96   // QA histograms
  97   fOutputList->Add(new TH2D("histQACentrality", ";centrality V0M(%);selected;",
  98                             fnBinsCent,fxMinCent,fxMaxCent, 2, -0.5, 1.5));
  99   fOutputList->Add(new TH2D("histQAVertexZ", ";vertex-z (cm);selected;",
 100                             800, -40., 40., 2, -0.5, 1.5));
 101   fOutputList->Add(new TH1D("histQAMultiplicityBeforeCuts",
 102                             ";tracks", 10000, 0, 10000));
 103   fOutputList->Add(new TH1D("histQAMultiplicityAfterCuts",
 104                             ";selected tracks", 10000, 0, 10000));
 105   fOutputList->Add(new TH3D("histQACentPt", ";charge;centrality V0M(%);p_{T} (GeV/c)",
 106                             2, -0.5, 1.5, fnBinsCent, fxMinCent, fxMaxCent, fnBinsPt, fxMinPt, fxMaxPt));
 107   fOutputList->Add(new TH3D("histQAPhiEta", ";charge;#phi (rad);#eta",
 108                             2, -0.5, 1.5, 200, 0.0, TMath::TwoPi(), 300, -1.5, 1.5));
 109
 110   // N(eta) distributions with different binnings
 111   fOutputList->Add(new TH2D("histNEta_300", ";#eta;N", 300, -1.5, 1.5, 1000, 0., 1000.));  // 0.01
 112   fOutputList->Add(new TH2D("histNEta_120", ";#eta;N", 120, -1.5, 1.5, 1000, 0., 1000.));  // 0.025
 113   fOutputList->Add(new TH2D("histNEta__30", ";#eta;N",  30, -1.5, 1.5, 1000, 0., 1000.));  // 0.1
 114   fOutputList->Add(new TH2D("histNEta__15", ";#eta;N",  15, -1.5, 1.5, 1000, 0., 1000.));  // 0.2
 115
 116   // Moments
 117   fOutputList->Add(MakeHistSparsePhiEta("histMoment1PhiEta_1"));
 118   fOutputList->Add(MakeHistSparsePhiEta("histMoment1PhiEta_2"));
 119   fOutputList->Add(MakeHistSparsePhiEtaPhiEta("histMoment2PhiEtaPhiEta_11"));
 120   fOutputList->Add(MakeHistSparsePhiEtaPhiEta("histMoment2PhiEtaPhiEta_12"));
 121   fOutputList->Add(MakeHistSparsePhiEtaPhiEta("histMoment2PhiEtaPhiEta_22"));
 122
 123   // add MC Histograms
 124   const Int_t N(fOutputList->GetEntries());
 125   for (Int_t i(0); i<N; ++i)
 126     fOutputList->Add(fOutputList->At(i)->Clone(TString("MC_") + fOutputList->At(i)->GetName()));
 127
 128   if (fRunMixing)
 129     SetupForMixing();
 130
 131   PostData(1, fOutputList);
 132 }
 133
 134 void AliAnalysisTaskLongRangeCorrelations::UserExec(Option_t* ) {
 135   AliAnalysisManager* pManager(AliAnalysisManager::GetAnalysisManager());
 136   if (NULL == pManager) return;
 137
 138   AliInputEventHandler* pInputHandler(dynamic_cast<AliInputEventHandler*>(pManager->GetInputEventHandler()));
 139   if (NULL == pInputHandler) return;
 140
 141   AliAODEvent* pAOD(dynamic_cast<AliAODEvent*>(InputEvent()));
 142   if (NULL == pAOD) return;
 143
 144   AliAODHeader *pAODHeader = pAOD->GetHeader();
 145   if (NULL == pAODHeader) return;
 146
 147   AliAODMCHeader *pAODMCHeader(dynamic_cast<AliAODMCHeader*>(pAOD->FindListObject(AliAODMCHeader::StdBranchName())));
 148   const Bool_t isMC(NULL != pAODMCHeader);
 149
 150   TClonesArray *arrayMC(NULL);
 151   if (isMC) {
 152     arrayMC = dynamic_cast<TClonesArray*>(pAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
 153     if (NULL == arrayMC) return;
 154   }
 155
 156   // --- event cuts MC ---
 157   if (isMC) {
 158     Fill("MC_histEventStats", 0.); // all events
 159     Fill("MC_histEventStats", 1.); // events passing physics selection
 160     Fill("MC_histEventStats", 2.); // events passing centrality selection
 161     // vertex cut in MC
 162     const Bool_t vertexSelected(TMath::Abs(pAODMCHeader->GetVtxZ()) < fMaxAbsVertexZ);
 163     Fill("MC_histQAVertexZ", pAODMCHeader->GetVtxZ(), vertexSelected);
 164     if (!vertexSelected) return;
 165     Fill("MC_histEventStats", 3.); // events passing vertex selection
 166   }
 167
 168   // --- event cuts data ---
 169   Fill("histEventStats", 0.); // all events
 170
 171   if (!pInputHandler->IsEventSelected()) return;
 172   Fill("histEventStats", 1.); // events passing physics selection
 173
 174   const Double_t centrality(pAODHeader->GetCentralityP()->GetCentralityPercentile("V0M"));
 175   AliDebug(3, Form("centrality=%f", centrality));
 176   const Bool_t centralitySelected(centrality > fCentMin && centrality < fCentMax);
 177   Fill("histQACentrality", centrality, centralitySelected);
 178   if (!centralitySelected) return;
 179   Fill("histEventStats", 2.); // events passing centrality selection
 180
 181   // vertex selection
 182   const Int_t nVertex(pAOD->GetNumberOfVertices());
 183   if (0 == nVertex) return;
 184   const AliAODVertex* pVertex(pAOD->GetPrimaryVertex());
 185   if (NULL == pVertex) return;
 186   const Int_t nTracksPrimary(pVertex->GetNContributors());
 187   if (nTracksPrimary < 1) return;
 188
 189   const Double_t zVertex(pVertex->GetZ());
 190   const Bool_t vertexSelected(TMath::Abs(zVertex) < fMaxAbsVertexZ);
 191   Fill("histQAVertexZ", zVertex, vertexSelected);
 192   if (!vertexSelected) return;
 193   Fill("histEventStats", 3.); // events passing vertex selection
 194
 195   // event is accepted
 196   TObjArray* tracksMain(GetAcceptedTracks(pAOD, centrality));
 197   Fill("histQAMultiplicityBeforeCuts", pAOD->GetNumberOfTracks());
 198   Fill("histQAMultiplicityAfterCuts", tracksMain->GetEntriesFast());
 199
 200   if (fRunMixing) {
 201     AliEventPool* pEventPool(fPoolMgr->GetEventPool(centrality, pAOD->GetPrimaryVertex()->GetZ()));
 202     if (NULL == pEventPool)
 203       AliFatal(Form("No pool found for centrality = %f, zVtx = %f", centrality, pAOD->GetPrimaryVertex()->GetZ()));
 204
 205 //     pEventPool->PrintInfo();
 206     if (pEventPool->IsReady()
 207         || pEventPool->NTracksInPool()     > fMixingTracks/10
 208         || pEventPool->GetCurrentNEvents() >= 5) {
 209       Fill("histEventStats", 4.); // analyzed events
 210       const Int_t nMix(pEventPool->GetCurrentNEvents());
 211       for (Int_t i(0); i<nMix; ++i) {
 212         TObjArray* tracksMixed(pEventPool->GetEvent(i));
 213         CalculateMoments("", tracksMain, tracksMixed, zVertex, 1./nMix);
 214       }
 215     }
 216     // Update the Event pool
 217     pEventPool->UpdatePool(tracksMain);
 218   } else { // no mixing
 219     Fill("histEventStats", 4.); // analyzed events
 220     CalculateMoments("", tracksMain, tracksMain, zVertex, 1.);
 221     delete tracksMain;
 222   }
 223
 224   if (isMC) {
 225     TObjArray* tracksMC(GetAcceptedTracks(arrayMC, centrality));
 226     const Double_t x[2] = {
 227       arrayMC->GetEntriesFast(),
 228       tracksMC->GetEntriesFast()
 229     };
 230     Fill("MC_histQAMultiplicity", x);
 231     Fill("MC_histEventStats", 4.); // analyzed MC events
 232     CalculateMoments("MC_", tracksMC, tracksMC, zVertex, 1.);
 233     delete tracksMC;
 234   }
 235 }
 236
 237 void AliAnalysisTaskLongRangeCorrelations::Terminate(Option_t* ) {
 238   //
 239 }
 240
 241 TString AliAnalysisTaskLongRangeCorrelations::GetOutputListName() const {
 242   TString listName("listLRC");
 243   listName += TString::Format("_%smix",         fRunMixing ? "" : "no");
 244   listName += TString::Format("_trackFilt%d",   fTrackFilter);
 245   listName += TString::Format("_cent%.0fT%.0f", fCentMin, fCentMax);
 246   listName += TString::Format("_ptMin%.0fMeV",  1e3*fPtMin);
 247   listName += TString::Format("_phi%.0fT%.0f",  TMath::RadToDeg()*fPhiMin, TMath::RadToDeg()*fPhiMax);
 248   return listName;
 249 }
 250
 251 void AliAnalysisTaskLongRangeCorrelations::SetupForMixing() {
 252   const Int_t trackDepth(fMixingTracks);
 253   const Int_t poolsize(1000); // Maximum number of events
 254
 255   Double_t centralityBins[] = { fCentMin, fCentMax };
 256   const Int_t nCentralityBins(sizeof(centralityBins)/sizeof(Double_t) - 1);
 257
 258   fPoolMgr = new AliEventPoolManager(poolsize, trackDepth,
 259                                      nCentralityBins, centralityBins,
 260                                      fVertexZ->GetNbins()+1,
 261                                      const_cast<Double_t*>(fVertexZ->GetXbins()->GetArray()));
 262 }
 263
 264 THnSparse* AliAnalysisTaskLongRangeCorrelations::MakeHistSparsePhiEta(const char* name) const {
 265   const Int_t nVertexZBins=fVertexZ->GetNbins();
 266   const Int_t   nBinsM1[] = { fnBinsPhi, fnBinsEta, nVertexZBins     };
 267   const Double_t xMinM1[] = {  fxMinPhi,  fxMinEta, 0.5              };
 268   const Double_t xMaxM1[] = {  fxMaxPhi,  fxMaxEta, nVertexZBins+0.5 };
 269   const TString title(TString(name)
 270                       +";#phi;#eta;vertex Z (bin);");
 271   return new THnSparseD(name, title.Data(), 3, nBinsM1, xMinM1, xMaxM1);
 272 }
 273 THnSparse* AliAnalysisTaskLongRangeCorrelations::MakeHistSparsePhiEtaPhiEta(const char* name) const {
 274   const Int_t nVertexZBins=fVertexZ->GetNbins();
 275   const Int_t   nBinsM2[] = {  fnBinsPhi, fnBinsEta,  fnBinsPhi, fnBinsEta, nVertexZBins     };
 276   const Double_t xMinM2[] = {  fxMinPhi,  fxMinEta,   fxMinPhi,  fxMinEta,  0.5              };
 277   const Double_t xMaxM2[] = {  fxMaxPhi,  fxMaxEta,   fxMaxPhi,  fxMaxEta,  nVertexZBins+0.5 };
 278   const TString title(TString(name)
 279                       +";#phi_{1};#eta_{1}"
 280                       +";#phi_{2};#eta_{2};vertex Z (bin);");
 281   return new THnSparseD(name, title.Data(), 5, nBinsM2, xMinM2, xMaxM2);
 282 }
 283
 284 TObjArray* AliAnalysisTaskLongRangeCorrelations::GetAcceptedTracks(AliAODEvent* pAOD,
 285                                                                    Double_t centrality) {
 286   TObjArray* tracks= new TObjArray;
 287   tracks->SetOwner(kTRUE);
 288
 289   const Long64_t N(pAOD->GetNumberOfTracks());
 290   AliDebug(5, Form("#tracks= %6lld %f", N, centrality));
 291   for (Long64_t i(0); i<N; ++i) {
 292     AliAODTrack* pAODTrack(dynamic_cast<AliAODTrack*>(pAOD->GetTrack(i)));
 293     if (NULL == pAODTrack) continue;
 294
 295     // track filter selection
 296     if (!pAODTrack->TestFilterBit(fTrackFilter)) continue;
 297
 298     // select only primary tracks
 299     if (pAODTrack->GetType() != AliAODTrack::kPrimary) continue;
 300
 301     // select only charged tracks
 302     if (pAODTrack->Charge() == 0) continue;
 303
 304     Fill("histQACentPt", pAODTrack->Charge()>0, centrality,       pAODTrack->Pt());
 305     Fill("histQAPhiEta", pAODTrack->Charge()>0, pAODTrack->Phi(), pAODTrack->Eta());
 306     if (pAODTrack->Phi() < fPhiMin || pAODTrack->Phi() > fPhiMax) continue;
 307     if (pAODTrack->Pt()  < fPtMin  || pAODTrack->Pt()  > fPtMax)  continue;
 308
 309     tracks->Add(new LRCParticle(pAODTrack->Eta(), pAODTrack->Phi()));
 310   } // next track
 311   return tracks;
 312 }
 313
 314 TObjArray* AliAnalysisTaskLongRangeCorrelations::GetAcceptedTracks(TClonesArray* tracksMC,
 315                                                                    Double_t centrality) {
 316   // for keeping track of MC labels
 317   std::set<Int_t> labelSet;
 318
 319   TObjArray* tracks= new TObjArray;
 320   tracks->SetOwner(kTRUE);
 321
 322   const Long64_t N(tracksMC->GetEntriesFast());
 323   AliDebug(5, Form("#tracks= %6lld %f", N, centrality));
 324   for (Long64_t i(0); i<N; ++i) {
 325     AliAODMCParticle* pMCTrack(dynamic_cast<AliAODMCParticle*>(tracksMC->At(i)));
 326     if (NULL == pMCTrack) continue;
 327
 328     // no track filter selection for MC tracks
 329
 330     if (labelSet.find(pMCTrack->Label()) != labelSet.end()) {
 331       Printf("Duplicate Label= %3d", pMCTrack->Label());
 332       continue;
 333     }
 334     labelSet.insert(pMCTrack->Label());
 335
 336     // select only primary tracks
 337     if (kFALSE == pMCTrack->IsPhysicalPrimary()) continue;
 338
 339     // select only charged tracks
 340     if (pMCTrack->Charge() == 0) continue;
 341
 342     Fill("MC_histQACentPt", pMCTrack->Charge()>0, centrality,      pMCTrack->Pt());
 343     Fill("MC_histQAPhiEta", pMCTrack->Charge()>0, pMCTrack->Phi(), pMCTrack->Eta());
 344
 345     if (pMCTrack->Phi() < fPhiMin || pMCTrack->Phi() > fPhiMax) continue;
 346     if (pMCTrack->Pt()  < fPtMin  || pMCTrack->Pt()  > fPtMax)  continue;
 347
 348     tracks->Add(new LRCParticle(pMCTrack->Eta(), pMCTrack->Phi()));
 349   } // next track
 350   return tracks;
 351 }
 352
 353 void AliAnalysisTaskLongRangeCorrelations::CalculateMoments(TString prefix,
 354                                                             TObjArray* tracks1,
 355                                                             TObjArray* tracks2,
 356                                                             Double_t vertexZ,
 357                                                             Double_t weight) {
 358   THnSparse* hN1(dynamic_cast<THnSparse*>(fOutputList->FindObject(prefix+"histMoment1PhiEta_1")));
 359   THnSparse* hN2(dynamic_cast<THnSparse*>(fOutputList->FindObject(prefix+"histMoment1PhiEta_2")));
 360   if (NULL == hN1) return;
 361   if (NULL == hN2) return;
 362
 363   // <n_1>
 364   THnSparse* hN1ForThisEvent(ComputeNForThisEvent(tracks1, "hN1", vertexZ));
 365   hN1->Add(hN1ForThisEvent, weight);
 366
 367   // <n_2>
 368   THnSparse* hN2ForThisEvent(ComputeNForThisEvent(tracks2, "hN2", vertexZ));
 369   hN2->Add(hN2ForThisEvent, weight);
 370
 371   // n(eta) distributions
 372   FillNEtaHist(prefix+"histNEta_300", hN1ForThisEvent, weight);
 373   FillNEtaHist(prefix+"histNEta_120", hN1ForThisEvent, weight);
 374   FillNEtaHist(prefix+"histNEta__30", hN1ForThisEvent, weight);
 375   FillNEtaHist(prefix+"histNEta__15", hN1ForThisEvent, weight);
 376
 377   TObjArray* hNs(new TObjArray);
 378
 379   // <n_1 n_1>
 380   hNs->AddAt(hN1ForThisEvent, 0);
 381   hNs->AddAt(hN1ForThisEvent, 1);
 382   ComputeNXForThisEvent(hNs, prefix+"histMoment2PhiEtaPhiEta_11", vertexZ, weight);
 383
 384   // <n_1 n_2>
 385   hNs->AddAt(hN2ForThisEvent, 1);
 386   ComputeNXForThisEvent(hNs, prefix+"histMoment2PhiEtaPhiEta_12", vertexZ, weight);
 387
 388   // <n_2 n_2>
 389   hNs->AddAt(hN2ForThisEvent, 0);
 390   ComputeNXForThisEvent(hNs, prefix+"histMoment2PhiEtaPhiEta_22", vertexZ, weight);
 391
 392   // clean up
 393   delete hNs;
 394   delete hN1ForThisEvent;
 395   delete hN2ForThisEvent;
 396 }
 397
 398 void AliAnalysisTaskLongRangeCorrelations::FillNEtaHist(TString name,
 399                                                         THnSparse* hs,
 400                                                         Double_t weight) {
 401
 402   TH2* hSum(dynamic_cast<TH2*>(fOutputList->FindObject(name)));
 403   if (NULL == hSum) return;
 404
 405   TH2* hPerEvent(dynamic_cast<TH2*>(hSum->Clone("hPerEvent")));
 406   if (NULL == hPerEvent) return;
 407   hPerEvent->Reset();
 408
 409   // fill hPerEvent
 410   const Long64_t N(hs->GetNbins());
 411   for (Long64_t i(0); i<N; ++i) {
 412     Int_t coord[2] = { 0, 0 };
 413     const Double_t n(hs->GetBinContent(i, coord));
 414     const Double_t eta(hs->GetAxis(1)->GetBinCenter(coord[1]));
 415     hPerEvent->Fill(eta, n);
 416   }
 417
 418   // add zero counts for those eta bins with zero tracks
 419   TH1* h(hPerEvent->ProjectionX());
 420   for (Int_t i(1); i<=h->GetNbinsX(); ++i)
 421     hPerEvent->SetBinContent(i,1, Double_t(h->GetBinContent(i) == 0));
 422
 423   hSum->Add(hPerEvent, weight);
 424
 425   delete h;
 426   delete hPerEvent;
 427 }
 428
 429 THnSparse* AliAnalysisTaskLongRangeCorrelations::ComputeNForThisEvent(TObjArray* tracks,
 430                                                                       const char* histName,
 431                                                                       Double_t vertexZ) const {
 432   const Double_t vertexZBin(fVertexZ->FindBin(vertexZ));
 433   THnSparse* hN(MakeHistSparsePhiEta(histName));
 434   const Long64_t nTracks(tracks->GetEntriesFast());
 435   for (Long64_t i(0); i<nTracks; ++i) {
 436     const LRCParticle* p(dynamic_cast<LRCParticle*>(tracks->At(i)));
 437     if (NULL == p) continue;
 438     const Double_t x[] = { p->Phi(), p->Eta(), vertexZBin };
 439     hN->Fill(x);
 440   }
 441   return hN;
 442 }
 443
 444 class MultiDimIterator {
 445 public:
 446   MultiDimIterator(TObjArray* _fHs, Double_t vertexZBin)
 447     : fHs(_fHs)
 448     , fN(fHs->GetEntriesFast())
 449     , fDims(fN,  0)
 450     , fIdxs(fN,  0)
 451     , fNs  (fN,  0)
 452     , fX (2*fN+1, 0)
 453     , fj(0) {
 454     for (Long64_t i=0; i<fN; ++i) {
 455       THnSparse* hNi(reinterpret_cast<THnSparse*>(fHs->At(i)));
 456       if (NULL == hNi)
 457         AliFatal("NULL == hNi");
 458       fDims[i] = hNi->GetNbins();
 459       AliDebug(3, Form("%lld %s %lld", i, hNi->GetName(), fDims[i]));
 460       update(i);
 461     }
 462     fX[2*fN] = vertexZBin;
 463   }
 464
 465   const char* ClassName() const { return "MultiDimIterator"; }
 466   const Double_t* GetX() const { return &fX.front(); }
 467   Double_t        GetN() const { // returns the product of all multiplicities
 468     return std::accumulate(fNs.begin(), fNs.end(), Double_t(1), std::multiplies<Double_t>());
 469   }
 470   Bool_t end() const { return fj == fN; }
 471
 472   MultiDimIterator& operator++() {
 473     Long64_t j(0);
 474     for (; j<fN; ++j) {
 475       ++fIdxs[j];
 476       update(j);
 477       if (fIdxs[j] < fDims[j]) break;
 478       fIdxs[j] = 0;
 479       update(j);
 480     }
 481     fj = j;
 482     return *this;
 483   }
 484
 485 protected:
 486   void update(Long64_t j) {
 487     THnSparse* hs(reinterpret_cast<THnSparse*>(fHs->At(j)));
 488     Int_t coord[] = { 0, 0 };
 489     fNs[j] = hs->GetBinContent(fIdxs[j], coord);
 490     for (Int_t k(0); k<2; ++k)
 491       fX[2*j+k] = hs->GetAxis(k)->GetBinCenter(coord[k]);
 492   }
 493 private:
 494   MultiDimIterator(const MultiDimIterator&);
 495   MultiDimIterator& operator=(const MultiDimIterator&);
 496
 497   TObjArray*            fHs;   // array of THnSparse histograms
 498   const Long64_t        fN;    // number of histograms
 499   std::vector<Long64_t> fDims; // number of filled bins for each THnSparse
 500   std::vector<Long64_t> fIdxs; // indices
 501   std::vector<Double_t> fNs;   // number of tracks
 502   std::vector<Double_t> fX;    // coordinate
 503   Long64_t              fj;    // state
 504 } ;
 505
 506 void AliAnalysisTaskLongRangeCorrelations::ComputeNXForThisEvent(TObjArray* hNs,
 507                                                                  const char* histName,
 508                                                                  Double_t vertexZ,
 509                                                                  Double_t weight) {
 510   if (NULL == fOutputList) return;
 511
 512   THnSparse* hs(dynamic_cast<THnSparse*>(fOutputList->FindObject(histName)));
 513   if (hs == NULL) return;
 514
 515   for (MultiDimIterator mdi(hNs, fVertexZ->FindBin(vertexZ)); !mdi.end(); ++mdi)
 516     hs->Fill(mdi.GetX(), mdi.GetN()*weight);
 517 }
 518
 519 void AliAnalysisTaskLongRangeCorrelations::Fill(const char* histName, Double_t x) {
 520   if (NULL == fOutputList) return;
 521   TH1* h = dynamic_cast<TH1*>(fOutputList->FindObject(histName));
 522   if (h == NULL) return;
 523   h->Fill(x);
 524 }
 525 void AliAnalysisTaskLongRangeCorrelations::Fill(const char* histName, Double_t x, Double_t y) {
 526   if (NULL == fOutputList) return;
 527   TH2* h = dynamic_cast<TH2*>(fOutputList->FindObject(histName));
 528   if (h == NULL) return;
 529   h->Fill(x, y);
 530 }
 531 void AliAnalysisTaskLongRangeCorrelations::Fill(const char* histName, Double_t x, Double_t y, Double_t z) {
 532   if (NULL == fOutputList) return;
 533   TH3* h = dynamic_cast<TH3*>(fOutputList->FindObject(histName));
 534   if (h == NULL) return;
 535   h->Fill(x, y, z);
 536 }
 537 void AliAnalysisTaskLongRangeCorrelations::Fill(const char* histName, const Double_t* x, Double_t w) {
 538   if (NULL == fOutputList) return;
 539   THnSparse* h = dynamic_cast<THnSparse*>(fOutputList->FindObject(histName));
 540   if (h == NULL) return;
 541   h->Fill(x, w);
 542 }