]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/Correlations/DPhi/AliAnalysisTaskLongRangeCorrelations.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Corrected end-of-line behavior
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / AliAnalysisTaskLongRangeCorrelations.cxx
CommitLineData
43d6f3d4 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
45ClassImp(AliAnalysisTaskLongRangeCorrelations);
46ClassImp(LRCParticle);
47
48AliAnalysisTaskLongRangeCorrelations::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
71AliAnalysisTaskLongRangeCorrelations::~AliAnalysisTaskLongRangeCorrelations() {
72 if (NULL != fOutputList) {
73 delete fOutputList;
74 fOutputList = NULL;
75 }
76}
77
78void 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
148void 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
267void 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
276TString 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
298void 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
311THnSparse* 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
321AliTHn* 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
334AliTHn* 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
348TObjArray* 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
407TObjArray* 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
459Bool_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
482void 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
542void 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
574THnSparse* 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
589class MultiDimIterator {
590public:
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
630protected:
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 }
638private:
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
651void 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
666void 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}
672void 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}
678void 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}
684void 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}