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