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 359 2013-10-29 15:39:09Z 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 #include "AliTHn.h"
  38
  39 #include "AliMCEventHandler.h"
  40 #include "AliMCEvent.h"
  41 #include "AliStack.h"
  42 #include "AliAnalysisTaskLongRangeCorrelations.h"
  43
  44
  45 ClassImp(AliAnalysisTaskLongRangeCorrelations);
  46 ClassImp(LRCParticle);
  47
  48 AliAnalysisTaskLongRangeCorrelations::AliAnalysisTaskLongRangeCorrelations(const char *name)
  49   : AliAnalysisTaskSE(name)
  50   , fOutputList(NULL)
  51   , fVertexZ(NULL)
  52   , fRunMixing(kFALSE)
  53   , fPoolMgr(NULL)
  54   , fMixingTracks(50000)
  55   , fTrackFilter(128)
  56   , fCentMin(0), fCentMax(20)
  57   , fPtMin(0.2), fPtMax(1e10)
  58   , fPhiMin(0.), fPhiMax(TMath::TwoPi())
  59   , fMaxAbsVertexZ(10.)
  60   , fSelectPrimaryMCParticles(0)
  61   , fSelectPrimaryMCDataParticles(0)
  62   , fNMin(-1)
  63   , fNMax(-1)
  64   , fnBinsCent( 220), fnBinsPt(400), fnBinsPhi(4),              fnBinsEta(120)
  65   , fxMinCent( -5.0), fxMinPt( 0.0), fxMinPhi( 0.0),            fxMinEta(-1.5)
  66   , fxMaxCent(105.0), fxMaxPt( 4.0), fxMaxPhi( TMath::TwoPi()), fxMaxEta( 1.5) {
  67   DefineInput(0, TChain::Class());
  68   DefineOutput(1, TList::Class());
  69 }
  70
  71 AliAnalysisTaskLongRangeCorrelations::~AliAnalysisTaskLongRangeCorrelations() {
  72   if (NULL != fOutputList) {
  73     delete fOutputList;
  74     fOutputList = NULL;
  75   }
  76 }
  77
  78 void AliAnalysisTaskLongRangeCorrelations::UserCreateOutputObjects() {
  79   fOutputList = new THashList;
  80   fOutputList->SetOwner(kTRUE);
  81   fOutputList->SetName(GetOutputListName());
  82
  83   const Double_t vertexBins[] = {  // Zvtx bins
  84     -10., -7., -5., -3., -1., 1., 3., 5., 7., 10.
  85   };
  86   const Int_t nVertexBins(sizeof(vertexBins)/sizeof(Double_t)-1);
  87   fVertexZ = new TAxis(nVertexBins, vertexBins);
  88   fVertexZ->SetName("VertexZAxis");
  89   fOutputList->Add(fVertexZ);
  90
  91   fOutputList->Add(new TH1D("histVertexZ", ";vertex Z (cm)", nVertexBins, vertexBins));
  92
  93   // Event statistics
  94   const char* eventStatLabels[] = {
  95     "All Events",
  96     "Physics Selection",
  97     "Centrality Selection",
  98     "Vertex Selection",
  99     "Analyzed Events"
 100   };
 101   const size_t nEventStat(sizeof(eventStatLabels)/sizeof(const char*));
 102   TH1* hStats(new TH1D("histEventStats", "histEventStats", nEventStat, -0.5, nEventStat-0.5));
 103   for (size_t i(0); i<nEventStat; ++i)
 104     hStats->GetXaxis()->SetBinLabel(1+i, eventStatLabels[i]);
 105   fOutputList->Add(hStats);
 106
 107   // QA histograms
 108   fOutputList->Add(new TH2D("histQACentrality", ";centrality V0M(%);selected;",
 109                             fnBinsCent,fxMinCent,fxMaxCent, 2, -0.5, 1.5));
 110   fOutputList->Add(new TH2D("histQAVertexZ", ";vertex-z (cm);selected;",
 111                             800, -40., 40., 2, -0.5, 1.5));
 112   fOutputList->Add(new TH1D("histQAMultiplicityBeforeCuts",
 113                             ";tracks", 10000, 0, 10000));
 114   fOutputList->Add(new TH1D("histQAMultiplicityAfterCuts",
 115                             ";selected tracks", 10000, 0, 10000));
 116   fOutputList->Add(new TH3D("histQACentPt", ";charge;centrality V0M(%);p_{T} (GeV/c)",
 117                             2, -0.5, 1.5, fnBinsCent, fxMinCent, fxMaxCent, fnBinsPt, fxMinPt, fxMaxPt));
 118   fOutputList->Add(new TH3D("histQAPhiEta", ";charge;#phi (rad);#eta",
 119                             2, -0.5, 1.5, 200, 0.0, TMath::TwoPi(), 300, -1.5, 1.5));
 120
 121   // N(eta) distributions with different binnings
 122   fOutputList->Add(new TH2D("histNEta_120", ";#eta;N", 120, -1.5, 1.5, 1000, -.5, 1000-.5));  // 0.025
 123   fOutputList->Add(new TH2D("histNEta__60", ";#eta;N",  60, -1.5, 1.5, 1000, -.5, 1000-.5));  // 0.05
 124   fOutputList->Add(new TH2D("histNEta__30", ";#eta;N",  30, -1.5, 1.5, 1000, -.5, 1000-.5));  // 0.1
 125   fOutputList->Add(new TH2D("histNEta__15", ";#eta;N",  15, -1.5, 1.5, 1000, -.5, 1000-.5));  // 0.2
 126
 127   // Moments
 128   fOutputList->Add(MakeHistPhiEta("histMoment1PhiEta_1"));
 129   if (fRunMixing)
 130     fOutputList->Add(MakeHistPhiEta("histMoment1PhiEta_2"));
 131   fOutputList->Add(MakeHistPhiEtaPhiEta("histMoment2PhiEtaPhiEta_11"));
 132   if (fRunMixing) {
 133     fOutputList->Add(MakeHistPhiEtaPhiEta("histMoment2PhiEtaPhiEta_12"));
 134     fOutputList->Add(MakeHistPhiEtaPhiEta("histMoment2PhiEtaPhiEta_22"));
 135   }
 136
 137   // add MC Histograms
 138   const Int_t N(fOutputList->GetEntries());
 139   for (Int_t i(0); i<N; ++i)
 140     fOutputList->Add(fOutputList->At(i)->Clone(TString("MC_") + fOutputList->At(i)->GetName()));
 141
 142   if (fRunMixing)
 143     SetupForMixing();
 144
 145   PostData(1, fOutputList);
 146 }
 147
 148 void AliAnalysisTaskLongRangeCorrelations::UserExec(Option_t* ) {
 149   AliAnalysisManager* pManager(AliAnalysisManager::GetAnalysisManager());
 150   if (NULL == pManager) return;
 151
 152   AliInputEventHandler* pInputHandler(dynamic_cast<AliInputEventHandler*>(pManager->GetInputEventHandler()));
 153   if (NULL == pInputHandler) return;
 154
 155   AliAODEvent* pAOD(dynamic_cast<AliAODEvent*>(InputEvent()));
 156   if (NULL == pAOD) return;
 157
 158   AliAODHeader *pAODHeader = pAOD->GetHeader();
 159   if (NULL == pAODHeader) return;
 160
 161   AliAODMCHeader *pAODMCHeader(dynamic_cast<AliAODMCHeader*>(pAOD->FindListObject(AliAODMCHeader::StdBranchName())));
 162   const Bool_t isMC(NULL != pAODMCHeader);
 163
 164   TClonesArray *arrayMC(NULL);
 165   if (isMC) {
 166     arrayMC = dynamic_cast<TClonesArray*>(pAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
 167     if (NULL == arrayMC) return;
 168   }
 169
 170   // --- event cuts MC ---
 171   if (isMC) {
 172     Fill("MC_histEventStats", 0.); // all events
 173     Fill("MC_histEventStats", 1.); // events passing physics selection
 174   }
 175
 176   // --- event cuts data ---
 177   Fill("histEventStats", 0.); // all events
 178
 179   // physics selection
 180   if (!pInputHandler->IsEventSelected()) return;
 181   Fill("histEventStats", 1.); // events passing physics selection
 182
 183   // centrality selection
 184   const Double_t centrality(pAODHeader->GetCentralityP()->GetCentralityPercentile("V0M"));
 185   AliDebug(3, Form("centrality=%f", centrality));
 186   const Bool_t centralitySelected(centrality > fCentMin && centrality < fCentMax);
 187   Fill("histQACentrality", centrality, centralitySelected);
 188   if (!centralitySelected) return;
 189   Fill("histEventStats", 2.); // events passing centrality selection
 190   if (isMC)
 191     Fill("MC_histEventStats", 2.); // events passing centrality selection
 192
 193   // vertex selection -- data
 194   const Int_t nVertex(pAOD->GetNumberOfVertices());
 195   if (0 == nVertex) return;
 196   const AliAODVertex* pVertex(pAOD->GetPrimaryVertex());
 197   if (NULL == pVertex) return;
 198   const Int_t nTracksPrimary(pVertex->GetNContributors());
 199   if (nTracksPrimary < 1) return;
 200
 201   const Double_t zVertex(pVertex->GetZ());
 202   const Bool_t vertexSelectedData(TMath::Abs(zVertex) < fMaxAbsVertexZ);
 203
 204   // vertex selection -- MC
 205   Bool_t vertexSelectedMC(kTRUE);
 206   if (isMC)
 207     vertexSelectedMC = (TMath::Abs(pAODMCHeader->GetVtxZ()) < fMaxAbsVertexZ);
 208
 209   // combined vertex selection (data and MC)
 210   const Bool_t vertexSelected(vertexSelectedData && vertexSelectedMC);
 211   Fill("histQAVertexZ", zVertex, vertexSelected);
 212   if (isMC)
 213     Fill("MC_histQAVertexZ", pAODMCHeader->GetVtxZ(), vertexSelected);
 214   if (!vertexSelected) return;
 215   Fill("histEventStats",    3.); // events passing vertex selection
 216   if (isMC)
 217     Fill("MC_histEventStats", 3.); // events passing vertex selection
 218
 219   // ------------------
 220   // event is accepted
 221   // ------------------
 222
 223   TObjArray* tracksMain(GetAcceptedTracks(pAOD, arrayMC, centrality));
 224   Fill("histQAMultiplicityBeforeCuts", pAOD->GetNumberOfTracks());
 225   Fill("histQAMultiplicityAfterCuts", tracksMain->GetEntriesFast());
 226
 227   if (fRunMixing) {
 228     AliEventPool* pEventPool(fPoolMgr->GetEventPool(centrality, pAOD->GetPrimaryVertex()->GetZ()));
 229     if (NULL == pEventPool)
 230       AliFatal(Form("No pool found for centrality = %f, zVtx = %f", centrality, pAOD->GetPrimaryVertex()->GetZ()));
 231
 232 //     pEventPool->PrintInfo();
 233     if (pEventPool->IsReady()
 234         || pEventPool->NTracksInPool()     > fMixingTracks/10
 235         || pEventPool->GetCurrentNEvents() >= 5) {
 236       Fill("histEventStats", 4.); // analyzed events
 237       Fill("histVertexZ", zVertex);
 238       const Int_t nMix(pEventPool->GetCurrentNEvents());
 239       for (Int_t i(0); i<nMix; ++i) {
 240         TObjArray* tracksMixed(pEventPool->GetEvent(i));
 241         CalculateMoments("", tracksMain, tracksMixed, zVertex, 1./nMix);
 242       }
 243     }
 244     // Update the Event pool
 245     pEventPool->UpdatePool(tracksMain);
 246   } else { // no mixing
 247     Fill("histEventStats", 4.); // analyzed events
 248     Fill("histVertexZ", zVertex);
 249     CalculateMoments("", tracksMain, tracksMain, zVertex, 1.);
 250     delete tracksMain;
 251   }
 252
 253   if (isMC) {
 254     TObjArray* tracksMC(GetAcceptedTracks(arrayMC, centrality));
 255     const Double_t x[2] = {
 256       arrayMC->GetEntriesFast(),
 257       tracksMC->GetEntriesFast()
 258     };
 259     Fill("MC_histQAMultiplicity", x);
 260     Fill("MC_histEventStats", 4.); // analyzed MC events
 261     Fill("MC_histVertexZ", pAOD->GetPrimaryVertex()->GetZ());
 262     CalculateMoments("MC_", tracksMC, tracksMC, zVertex, 1.);
 263     delete tracksMC;
 264   }
 265 }
 266
 267 void AliAnalysisTaskLongRangeCorrelations::Terminate(Option_t* ) {
 268   //
 269   fOutputList = dynamic_cast<TList*>(GetOutputData(1));
 270   if (NULL == fOutputList) {
 271     AliFatal("NULL == fOutputList");
 272     return; // needed to avoid coverity warning
 273   }
 274 }
 275
 276 TString AliAnalysisTaskLongRangeCorrelations::GetOutputListName() const {
 277   TString listName("listLRC");
 278   listName += TString::Format("_%smix",         fRunMixing ? "" : "no");
 279   listName += TString::Format("_trackFilt%d",   fTrackFilter);
 280   listName += TString::Format("_cent%.0fT%.0f", fCentMin, fCentMax);
 281   listName += TString::Format("_ptMin%.0fMeV",  1e3*fPtMin);
 282   listName += TString::Format("_phi%.0fT%.0f",  TMath::RadToDeg()*fPhiMin, TMath::RadToDeg()*fPhiMax);
 283   if ( 1 == fSelectPrimaryMCParticles)
 284     listName += "_selPrimMC";
 285   if (-1 == fSelectPrimaryMCParticles)
 286     listName += "_selNonPrimMC";
 287   if ( 1 == fSelectPrimaryMCDataParticles)
 288     listName += "_selPrimMCData";
 289   if (-1 == fSelectPrimaryMCDataParticles)
 290     listName += "_selNonPrimMCData";
 291   if (-1 != fNMin)
 292     listName += TString::Format("_nMin%d", fNMin);
 293   if (-1 != fNMax)
 294     listName += TString::Format("_nMax%d", fNMax);
 295   return listName;
 296 }
 297
 298 void AliAnalysisTaskLongRangeCorrelations::SetupForMixing() {
 299   const Int_t trackDepth(fMixingTracks);
 300   const Int_t poolsize(1000); // Maximum number of events
 301
 302   Double_t centralityBins[] = { fCentMin, fCentMax };
 303   const Int_t nCentralityBins(sizeof(centralityBins)/sizeof(Double_t) - 1);
 304
 305   fPoolMgr = new AliEventPoolManager(poolsize, trackDepth,
 306                                      nCentralityBins, centralityBins,
 307                                      fVertexZ->GetNbins()+1,
 308                                      const_cast<Double_t*>(fVertexZ->GetXbins()->GetArray()));
 309 }
 310
 311 THnSparse* AliAnalysisTaskLongRangeCorrelations::MakeHistSparsePhiEta(const char* name) const {
 312   const Int_t nVertexZBins=fVertexZ->GetNbins();
 313   const Int_t   nBinsM1[] = { fnBinsPhi, fnBinsEta, nVertexZBins     };
 314   const Double_t xMinM1[] = {  fxMinPhi,  fxMinEta, 0.5              };
 315   const Double_t xMaxM1[] = {  fxMaxPhi,  fxMaxEta, nVertexZBins+0.5 };
 316   const TString title(TString(name)
 317                       +";#phi;#eta;vertex Z (bin);");
 318   return new THnSparseD(name, title.Data(), 3, nBinsM1, xMinM1, xMaxM1);
 319 }
 320
 321 AliTHn* AliAnalysisTaskLongRangeCorrelations::MakeHistPhiEta(const char* name) const {
 322   const Int_t nVertexZBins=fVertexZ->GetNbins();
 323   const Int_t   nBinsM1[] = { fnBinsPhi, fnBinsEta, nVertexZBins     };
 324   const Double_t xMinM1[] = {  fxMinPhi,  fxMinEta, 0.5              };
 325   const Double_t xMaxM1[] = {  fxMaxPhi,  fxMaxEta, nVertexZBins+0.5 };
 326   const TString title(TString(name)
 327                       +";#phi;#eta;vertex Z (bin);");
 328   AliTHn* h(new AliTHn(name, title.Data(), 1, 3, nBinsM1));
 329   for (Int_t i=0; i<3; ++i)
 330     h->SetBinLimits(i, xMinM1[i], xMaxM1[i]);
 331   return h;
 332 }
 333
 334 AliTHn* AliAnalysisTaskLongRangeCorrelations::MakeHistPhiEtaPhiEta(const char* name) const {
 335   const Int_t nVertexZBins=fVertexZ->GetNbins();
 336   const Int_t   nBinsM2[] = {  fnBinsPhi, fnBinsEta,  fnBinsPhi, fnBinsEta, nVertexZBins     };
 337   const Double_t xMinM2[] = {  fxMinPhi,  fxMinEta,   fxMinPhi,  fxMinEta,  0.5              };
 338   const Double_t xMaxM2[] = {  fxMaxPhi,  fxMaxEta,   fxMaxPhi,  fxMaxEta,  nVertexZBins+0.5 };
 339   const TString title(TString(name)
 340                       +";#phi_{1};#eta_{1}"
 341                       +";#phi_{2};#eta_{2};vertex Z (bin);");
 342   AliTHn* h(new AliTHn(name, title.Data(), 1, 5, nBinsM2));
 343   for (Int_t i=0; i<5; ++i)
 344     h->SetBinLimits(i, xMinM2[i], xMaxM2[i]);
 345   return h;
 346 }
 347
 348 TObjArray* AliAnalysisTaskLongRangeCorrelations::GetAcceptedTracks(AliAODEvent*  pAOD,
 349                                                                    TClonesArray* arrayMC,
 350                                                                    Double_t      centrality) {
 351   TObjArray* tracks= new TObjArray;
 352   tracks->SetOwner(kTRUE);
 353
 354   const Long64_t N(pAOD->GetNumberOfTracks());
 355   AliDebug(5, Form("#tracks= %6lld %f", N, centrality));
 356   for (Long64_t i(0); i<N; ++i) {
 357     AliAODTrack* pAODTrack(dynamic_cast<AliAODTrack*>(pAOD->GetTrack(i)));
 358     if (NULL == pAODTrack) continue;
 359
 360     // track filter selection
 361     if (!pAODTrack->TestFilterBit(fTrackFilter)) continue;
 362
 363     // select only primary tracks
 364     if (pAODTrack->GetType() != AliAODTrack::kPrimary) continue;
 365
 366     if (NULL != arrayMC) {
 367       const Int_t label(pAODTrack->GetLabel());
 368       AliAODMCParticle* mcParticle((label >= 0)
 369                                    ? static_cast<AliAODMCParticle*>(arrayMC->At(label))
 370                                    : NULL);
 371       if (label >=0 && NULL == mcParticle)
 372         AliFatal("MC particle not found");
 373
 374 //       const Bool_t isPhysicalPrimary(mcParticle->IsPhysicalPrimary());
 375 //       Printf("mcParticle->IsPhysicalPrimary() %d", isPhysicalPrimary);
 376
 377       switch (fSelectPrimaryMCDataParticles) {
 378       case -1:
 379         if (label < 0) continue;
 380         if (kTRUE  == mcParticle->IsPhysicalPrimary()) continue;
 381         break;
 382       case  0:
 383         // NOP, take all tracks
 384         break;
 385       case  1:
 386         if (label < 0) continue;
 387         if (kFALSE == mcParticle->IsPhysicalPrimary()) continue;
 388         break;
 389       default:
 390         AliFatal("fSelectPrimaryMCDataParticles != {-1,0,1}");
 391       }
 392     }
 393
 394     // select only charged tracks
 395     if (pAODTrack->Charge() == 0) continue;
 396
 397     Fill("histQACentPt", pAODTrack->Charge()>0, centrality,       pAODTrack->Pt());
 398     Fill("histQAPhiEta", pAODTrack->Charge()>0, pAODTrack->Phi(), pAODTrack->Eta());
 399     if (pAODTrack->Phi() < fPhiMin || pAODTrack->Phi() > fPhiMax) continue;
 400     if (pAODTrack->Pt()  < fPtMin  || pAODTrack->Pt()  > fPtMax)  continue;
 401
 402     tracks->Add(new LRCParticle(pAODTrack->Eta(), pAODTrack->Phi()));
 403   } // next track
 404   return tracks;
 405 }
 406
 407 TObjArray* AliAnalysisTaskLongRangeCorrelations::GetAcceptedTracks(TClonesArray* tracksMC,
 408                                                                    Double_t centrality) {
 409   // for keeping track of MC labels
 410   std::set<Int_t> labelSet;
 411
 412   TObjArray* tracks= new TObjArray;
 413   tracks->SetOwner(kTRUE);
 414
 415   const Long64_t N(tracksMC->GetEntriesFast());
 416   AliDebug(5, Form("#tracks= %6lld %f", N, centrality));
 417   for (Long64_t i(0); i<N; ++i) {
 418     AliAODMCParticle* pMCTrack(dynamic_cast<AliAODMCParticle*>(tracksMC->At(i)));
 419     if (NULL == pMCTrack) continue;
 420
 421     // no track filter selection for MC tracks
 422
 423     if (labelSet.find(pMCTrack->Label()) != labelSet.end()) {
 424       Printf("Duplicate Label= %3d", pMCTrack->Label());
 425       continue;
 426     }
 427     labelSet.insert(pMCTrack->Label());
 428
 429 //     Printf("isPrim = %d", pMCTrack->IsPhysicalPrimary());
 430
 431     switch (fSelectPrimaryMCParticles) {
 432     case -1:
 433       if (kTRUE == pMCTrack->IsPhysicalPrimary()) continue;
 434       break;
 435     case  0:
 436       // NOP, take all MC tracks
 437       break;
 438     case  1:
 439       if (kFALSE == pMCTrack->IsPhysicalPrimary()) continue;
 440       break;
 441     default:
 442       AliFatal("fSelectPrimaryMCParticles != {-1,0,1}");
 443     }
 444
 445     // select only charged tracks
 446     if (pMCTrack->Charge() == 0) continue;
 447
 448     Fill("MC_histQACentPt", pMCTrack->Charge()>0, centrality,      pMCTrack->Pt());
 449     Fill("MC_histQAPhiEta", pMCTrack->Charge()>0, pMCTrack->Phi(), pMCTrack->Eta());
 450
 451     if (pMCTrack->Phi() < fPhiMin || pMCTrack->Phi() > fPhiMax) continue;
 452     if (pMCTrack->Pt()  < fPtMin  || pMCTrack->Pt()  > fPtMax)  continue;
 453
 454     tracks->Add(new LRCParticle(pMCTrack->Eta(), pMCTrack->Phi()));
 455   } // next track
 456   return tracks;
 457 }
 458
 459 Bool_t AddTHnSparseToAliTHn(AliTHn* h, THnSparse* hs, Double_t weight) {
 460   if (h->GetNVar() != hs->GetNdimensions())
 461     return kFALSE;
 462
 463   const size_t nDim(hs->GetNdimensions());
 464
 465   Int_t    *coord = new Int_t[nDim];
 466   Double_t *x     = new Double_t[nDim];
 467
 468   const Long64_t N(hs->GetNbins());
 469   for (Long64_t i(0); i<N; ++i) {
 470     const Double_t n(hs->GetBinContent(i, coord));
 471     for (size_t j=0; j<nDim; ++j)
 472       x[j] = hs->GetAxis(j)->GetBinCenter(coord[j]);
 473     h->Fill(x, 0, n*weight);
 474   }
 475
 476   delete[] coord;
 477   delete[] x;
 478
 479   return kTRUE;
 480 }
 481
 482 void AliAnalysisTaskLongRangeCorrelations::CalculateMoments(TString prefix,
 483                                                             TObjArray* tracks1,
 484                                                             TObjArray* tracks2,
 485                                                             Double_t vertexZ,
 486                                                             Double_t weight) {
 487   const Long64_t nc1(tracks1->GetEntriesFast());
 488   if (fNMin != -1 && nc1 < fNMin) return;
 489   if (fNMax != -1 && nc1 > fNMax) return;
 490
 491   const Long64_t nc2(tracks1->GetEntriesFast());
 492   if (fNMin != -1 && nc2 < fNMin) return;
 493   if (fNMax != -1 && nc2 > fNMax) return;
 494
 495   AliTHn* hN1(dynamic_cast<AliTHn*>(fOutputList->FindObject(prefix+"histMoment1PhiEta_1")));
 496   if (NULL == hN1) return;
 497
 498   // <n_1>
 499   THnSparse* hN1ForThisEvent(ComputeNForThisEvent(tracks1, "hN1", vertexZ));
 500   AddTHnSparseToAliTHn(hN1,hN1ForThisEvent, weight);
 501 //   hN1->GetGrid(0)->GetGrid()->Add(hN1ForThisEvent, weight);
 502
 503   // n(eta) distributions
 504   FillNEtaHist(prefix+"histNEta_300", hN1ForThisEvent, weight);
 505   FillNEtaHist(prefix+"histNEta_120", hN1ForThisEvent, weight);
 506   FillNEtaHist(prefix+"histNEta__30", hN1ForThisEvent, weight);
 507   FillNEtaHist(prefix+"histNEta__15", hN1ForThisEvent, weight);
 508
 509   TObjArray* hNs(new TObjArray);
 510
 511   // <n_1 n_1>
 512   hNs->AddAt(hN1ForThisEvent, 0);
 513   hNs->AddAt(hN1ForThisEvent, 1);
 514
 515   ComputeNXForThisEvent(hNs, prefix+"histMoment2PhiEtaPhiEta_11", vertexZ, weight);
 516
 517   if (fRunMixing) {
 518     AliTHn* hN2(dynamic_cast<AliTHn*>(fOutputList->FindObject(prefix+"histMoment1PhiEta_2")));
 519     if (NULL == hN2) return;
 520     // <n_2>
 521     THnSparse* hN2ForThisEvent(ComputeNForThisEvent(tracks2, "hN2", vertexZ));
 522     AddTHnSparseToAliTHn(hN2,hN2ForThisEvent, weight);
 523 //     hN2->GetGrid(0)->GetGrid()->Add(hN2ForThisEvent, weight);
 524
 525     // <n_1 n_2>
 526     hNs->AddAt(hN2ForThisEvent, 1);
 527     ComputeNXForThisEvent(hNs, prefix+"histMoment2PhiEtaPhiEta_12", vertexZ, weight);
 528
 529     // <n_2 n_2>
 530     hNs->AddAt(hN2ForThisEvent, 0);
 531     ComputeNXForThisEvent(hNs, prefix+"histMoment2PhiEtaPhiEta_22", vertexZ, weight);
 532
 533     // clean up
 534     delete hN2ForThisEvent;
 535   }
 536
 537   // clean up
 538   delete hNs;
 539   delete hN1ForThisEvent;
 540 }
 541
 542 void AliAnalysisTaskLongRangeCorrelations::FillNEtaHist(TString name,
 543                                                         THnSparse* hs,
 544                                                         Double_t weight) {
 545
 546   TH2* hSum(dynamic_cast<TH2*>(fOutputList->FindObject(name)));
 547   if (NULL == hSum) return;
 548
 549   TH2* hPerEvent(dynamic_cast<TH2*>(hSum->Clone("hPerEvent")));
 550   if (NULL == hPerEvent) return;
 551   hPerEvent->Reset();
 552
 553   TH1* h1PerEvent(hPerEvent->ProjectionX());
 554
 555   // fill h1PerEvent
 556   const Long64_t N(hs->GetNbins());
 557   for (Long64_t i(0); i<N; ++i) {
 558     Int_t coord[2] = { 0, 0 };
 559     const Double_t n(hs->GetBinContent(i, coord));
 560     const Double_t eta(hs->GetAxis(1)->GetBinCenter(coord[1]));
 561     h1PerEvent->Fill(eta, n);
 562   }
 563
 564   for (Int_t i(1); i<=h1PerEvent->GetNbinsX(); ++i)
 565     hPerEvent->Fill(h1PerEvent->GetXaxis()->GetBinCenter(i),
 566                     h1PerEvent->GetBinContent(i));
 567
 568   hSum->Add(hPerEvent, weight);
 569
 570   delete hPerEvent;
 571   delete h1PerEvent;
 572 }
 573
 574 THnSparse* AliAnalysisTaskLongRangeCorrelations::ComputeNForThisEvent(TObjArray* tracks,
 575                                                                       const char* histName,
 576                                                                       Double_t vertexZ) const {
 577   const Double_t vertexZBin(fVertexZ->FindBin(vertexZ));
 578   THnSparse* hN(MakeHistSparsePhiEta(histName));
 579   const Long64_t nTracks(tracks->GetEntriesFast());
 580   for (Long64_t i(0); i<nTracks; ++i) {
 581     const LRCParticle* p(dynamic_cast<LRCParticle*>(tracks->At(i)));
 582     if (NULL == p) continue;
 583     const Double_t x[] = { p->Phi(), p->Eta(), vertexZBin };
 584     hN->Fill(x);
 585   }
 586   return hN;
 587 }
 588
 589 class MultiDimIterator {
 590 public:
 591   MultiDimIterator(TObjArray* _fHs, Double_t vertexZBin)
 592     : fHs(_fHs)
 593     , fN(fHs->GetEntriesFast())
 594     , fDims(fN,  0)
 595     , fIdxs(fN,  0)
 596     , fNs  (fN,  0)
 597     , fX (2*fN+1, 0)
 598     , fj(0) {
 599     for (Long64_t i=0; i<fN; ++i) {
 600       THnSparse* hNi(reinterpret_cast<THnSparse*>(fHs->At(i)));
 601       if (NULL == hNi)
 602         AliFatal("NULL == hNi");
 603       fDims[i] = hNi->GetNbins();
 604       AliDebug(3, Form("%lld %s %lld", i, hNi->GetName(), fDims[i]));
 605       update(i);
 606     }
 607     fX[2*fN] = vertexZBin;
 608   }
 609
 610   const char* ClassName() const { return "MultiDimIterator"; }
 611   const Double_t* GetX() const { return &fX.front(); }
 612   Double_t        GetN() const { // returns the product of all multiplicities
 613     return std::accumulate(fNs.begin(), fNs.end(), Double_t(1), std::multiplies<Double_t>());
 614   }
 615   Bool_t end() const { return fj == fN; }
 616
 617   MultiDimIterator& operator++() {
 618     Long64_t j(0);
 619     for (; j<fN; ++j) {
 620       ++fIdxs[j];
 621       update(j);
 622       if (fIdxs[j] < fDims[j]) break;
 623       fIdxs[j] = 0;
 624       update(j);
 625     }
 626     fj = j;
 627     return *this;
 628   }
 629
 630 protected:
 631   void update(Long64_t j) {
 632     THnSparse* hs(reinterpret_cast<THnSparse*>(fHs->At(j)));
 633     Int_t coord[] = { 0, 0 };
 634     fNs[j] = hs->GetBinContent(fIdxs[j], coord);
 635     for (Int_t k(0); k<2; ++k)
 636       fX[2*j+k] = hs->GetAxis(k)->GetBinCenter(coord[k]);
 637   }
 638 private:
 639   MultiDimIterator(const MultiDimIterator&);
 640   MultiDimIterator& operator=(const MultiDimIterator&);
 641
 642   TObjArray*            fHs;   // array of THnSparse histograms
 643   const Long64_t        fN;    // number of histograms
 644   std::vector<Long64_t> fDims; // number of filled bins for each THnSparse
 645   std::vector<Long64_t> fIdxs; // indices
 646   std::vector<Double_t> fNs;   // number of tracks
 647   std::vector<Double_t> fX;    // coordinate
 648   Long64_t              fj;    // state
 649 } ;
 650
 651 void AliAnalysisTaskLongRangeCorrelations::ComputeNXForThisEvent(TObjArray* hNs,
 652                                                                  const char* histName,
 653                                                                  Double_t vertexZ,
 654                                                                  Double_t weight) {
 655   if (NULL == fOutputList) return;
 656
 657   AliTHn* hs(dynamic_cast<AliTHn*>(fOutputList->FindObject(histName)));
 658   if (hs == NULL) return;
 659
 660   for (MultiDimIterator mdi(hNs, fVertexZ->FindBin(vertexZ)); !mdi.end(); ++mdi)
 661     hs->Fill(mdi.GetX(), 0, mdi.GetN()*weight);
 662
 663 //   hs->FillParent();
 664 }
 665
 666 void AliAnalysisTaskLongRangeCorrelations::Fill(const char* histName, Double_t x) {
 667   if (NULL == fOutputList) return;
 668   TH1* h = dynamic_cast<TH1*>(fOutputList->FindObject(histName));
 669   if (h == NULL) return;
 670   h->Fill(x);
 671 }
 672 void AliAnalysisTaskLongRangeCorrelations::Fill(const char* histName, Double_t x, Double_t y) {
 673   if (NULL == fOutputList) return;
 674   TH2* h = dynamic_cast<TH2*>(fOutputList->FindObject(histName));
 675   if (h == NULL) return;
 676   h->Fill(x, y);
 677 }
 678 void AliAnalysisTaskLongRangeCorrelations::Fill(const char* histName, Double_t x, Double_t y, Double_t z) {
 679   if (NULL == fOutputList) return;
 680   TH3* h = dynamic_cast<TH3*>(fOutputList->FindObject(histName));
 681   if (h == NULL) return;
 682   h->Fill(x, y, z);
 683 }
 684 void AliAnalysisTaskLongRangeCorrelations::Fill(const char* histName, const Double_t* x, Double_t w) {
 685   if (NULL == fOutputList) return;
 686   THnSparse* h = dynamic_cast<THnSparse*>(fOutputList->FindObject(histName));
 687   if (h == NULL) return;
 688   h->Fill(x, w);
 689 }