]> 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
Changes to compile with Root6 on macosx64
[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 **************************************************************************/
e474b5d3 16// $Id: AliAnalysisTaskLongRangeCorrelations.cxx 407 2014-03-21 11:55:57Z cmayer $
43d6f3d4 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)
a742e0d7 64 , fDeltaEta(-1)
43d6f3d4 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
72AliAnalysisTaskLongRangeCorrelations::~AliAnalysisTaskLongRangeCorrelations() {
73 if (NULL != fOutputList) {
74 delete fOutputList;
75 fOutputList = NULL;
76 }
77}
78
79void 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
149void 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] = {
2942f542 257 static_cast<Double_t>(arrayMC->GetEntriesFast()),
258 static_cast<Double_t>(tracksMC->GetEntriesFast())
43d6f3d4 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
268void 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
277TString 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);
a742e0d7 296 if (fDeltaEta >= 0)
297 listName += TString::Format("_deltaEta%02.0f", 10*fDeltaEta);
43d6f3d4 298 return listName;
299}
300
301void 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
314THnSparse* 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
324AliTHn* 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
337AliTHn* 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
351TObjArray* 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
410TObjArray* 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
462Bool_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
485void AliAnalysisTaskLongRangeCorrelations::CalculateMoments(TString prefix,
486 TObjArray* tracks1,
487 TObjArray* tracks2,
488 Double_t vertexZ,
489 Double_t weight) {
a742e0d7 490 THnSparse* hN1ForThisEvent(ComputeNForThisEvent(tracks1, "hN1", vertexZ));
491
e474b5d3 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);
a742e0d7 505
e474b5d3 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 }
43d6f3d4 512
513 AliTHn* hN1(dynamic_cast<AliTHn*>(fOutputList->FindObject(prefix+"histMoment1PhiEta_1")));
514 if (NULL == hN1) return;
515
516 // <n_1>
43d6f3d4 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
559void 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
591THnSparse* 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
606class MultiDimIterator {
607public:
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
647protected:
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 }
655private:
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
668void 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
683void 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}
689void 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}
695void 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}
701void 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}