update DHC task (Constantin Loizides <cloizides@lbl.gov>)
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / FourierDecomposition / AliDhcTask.cxx
CommitLineData
48a61f36 1// Dihadron correlations task - simple task to read ESD or AOD input,
2// calculate same- and mixed-event correlations, and fill THnSparse
3// output. -A. Adare, Apr 2011
4
5#include "TCanvas.h"
6#include "TChain.h"
8433ff41 7#include "TFormula.h"
48a61f36 8#include "TH1.h"
9#include "TH2.h"
beb13b6c 10#include "TH3.h"
d688e049 11#include "TAxis.h"
8433ff41 12#include "TProfile2D.h"
48a61f36 13#include "TROOT.h"
14#include "TTree.h"
48a61f36 15#include "AliAODEvent.h"
16#include "AliAODInputHandler.h"
17#include "AliAnalysisManager.h"
18#include "AliAnalysisTask.h"
19#include "AliCentrality.h"
20#include "AliDhcTask.h"
21#include "AliESDEvent.h"
22#include "AliESDInputHandler.h"
23#include "AliESDtrackCuts.h"
beb13b6c 24#include "AliESDMuonTrack.h"
8433ff41 25#include "AliPool.h"
48a61f36 26#include "AliVParticle.h"
27
c64cb1f6 28using std::cout;
29using std::endl;
30
48a61f36 31ClassImp(AliDhcTask)
32
beb13b6c 33
34//________________________________________________________________________
35AliDhcTask::AliDhcTask()
36: AliAnalysisTaskSE(), fVerbosity(0), fEtaMax(1), fZVtxMax(10), fPtMin(0.25), fPtMax(15),
37 fTrackDepth(1000), fPoolSize(200), fTracksName(), fDoWeights(kFALSE), fFillMuons(kFALSE),
38 fEtaTLo(-1.0), fEtaTHi(1.0), fEtaALo(-1.0), fEtaAHi(1.0),
39 fESD(0x0), fAOD(0x0), fOutputList(0x0), fHEvt(0x0), fHTrk(0x0),
40 fHPtAss(0x0), fHPtTrg(0x0), fHPtTrgEvt(0x0),
41 fHPtTrgNorm1S(0x0), fHPtTrgNorm1M(0x0), fHPtTrgNorm2S(0x0), fHPtTrgNorm2M(0x0),
42 fHCent(0x0), fHZvtx(0x0), fNbins(0), fHSs(0x0), fHMs(0x0), fHPts(0x0), fIndex(0x0),
43 fCentrality(99), fZVertex(99), fEsdTPCOnly(0), fPoolMgr(0),
44 fCentMethod("V0M"), fNBdeta(20), fNBdphi(36),
45 fBPtT(0x0), fBPtA(0x0), fBCent(0x0), fBZvtx(0x0),
46 fMixBCent(0x0), fMixBZvtx(0x0),
47 fHEffT(0x0), fHEffA(0x0)
48{
49
50}
51
48a61f36 52//________________________________________________________________________
53AliDhcTask::AliDhcTask(const char *name)
42036329 54: AliAnalysisTaskSE(name), fVerbosity(0), fEtaMax(1), fZVtxMax(10), fPtMin(0.25), fPtMax(15),
beb13b6c 55 fTrackDepth(1000), fPoolSize(200), fTracksName(), fDoWeights(kFALSE), fFillMuons(kFALSE),
56 fEtaTLo(-1.0), fEtaTHi(1.0), fEtaALo(-1.0), fEtaAHi(1.0),
57 fESD(0x0), fAOD(0x0), fOutputList(0x0), fHEvt(0x0), fHTrk(0x0),
58 fHPtAss(0x0), fHPtTrg(0x0), fHPtTrgEvt(0x0),
59 fHPtTrgNorm1S(0x0), fHPtTrgNorm1M(0x0), fHPtTrgNorm2S(0x0), fHPtTrgNorm2M(0x0),
60 fHCent(0x0), fHZvtx(0x0), fNbins(0), fHSs(0x0), fHMs(0x0), fHPts(0x0), fIndex(0x0),
8a9d3e12 61 fCentrality(99), fZVertex(99), fEsdTPCOnly(0), fPoolMgr(0),
d688e049 62 fCentMethod("V0M"), fNBdeta(20), fNBdphi(36),
63 fBPtT(0x0), fBPtA(0x0), fBCent(0x0), fBZvtx(0x0),
beb13b6c 64 fMixBCent(0x0), fMixBZvtx(0x0),
65 fHEffT(0x0), fHEffA(0x0)
48a61f36 66{
67 // Constructor
68
69 // Define input and output slots here
70 // Input slot #0 works with a TChain
71 DefineInput(0, TChain::Class());
72 // Output slot #0 id reserved by the base class for AOD
73 // Output slot #1 writes into a TH1 container
74 DefineOutput(1, TList::Class());
75
48a61f36 76 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,Tracks "
77 "AOD:header,tracks,vertices,";
d688e049 78
79 Double_t ptt[4] = {0.25, 1.0, 2.0, 15.0};
80 fBPtT = new TAxis(3,ptt);
81 Double_t pta[4] = {0.25, 1.0, 2.0, 15.0};
82 fBPtA = new TAxis(3,pta);
83 Double_t cent[2] = {-100.0, 100.0};
84 fBCent = new TAxis(1,cent);
85 Double_t zvtx[2] = {-10, 10};
86 fBZvtx = new TAxis(1,zvtx);
87 Double_t centmix[2] = {-100.0, 100.0};
88 fMixBCent = new TAxis(1,centmix);
89 Double_t zvtxmix[9] = {-10,-6,-4,-2,0,2,4,6,10};
90 fMixBZvtx = new TAxis(8,zvtxmix);
48a61f36 91}
92
93//________________________________________________________________________
94void AliDhcTask::UserCreateOutputObjects()
95{
96 // Create histograms
97 // Called once (per slave on PROOF!)
beb13b6c 98 if (fVerbosity > 1) {
99 AliInfo("Initialize Dhc Task");
100 AliInfo(Form(" efficiency correct triggers? %d", fHEffT!=0));
101 AliInfo(Form(" efficiency correct associates? %d", fHEffA!=0));
102 AliInfo(Form(" fill muons? %d", fFillMuons));
103 AliInfo(Form(" trigger eta range %f .. %f", fEtaTLo, fEtaTHi));
104 AliInfo(Form(" associate eta range %f .. %f", fEtaALo, fEtaAHi));
105 }
48a61f36 106
48a61f36 107 fOutputList = new TList();
108 fOutputList->SetOwner(1);
109
8a9d3e12 110 fEsdTPCOnly = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
111 //fEsdTPCOnly->SetMinNClustersTPC(70);
112 fEsdTPCOnly->SetMinNCrossedRowsTPC(70);
113 fEsdTPCOnly->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
48a61f36 114
8433ff41 115 BookHistos();
116 InitEventMixer();
48a61f36 117 PostData(1, fOutputList);
118}
119
120//________________________________________________________________________
121void AliDhcTask::BookHistos()
122{
8433ff41 123 // Book histograms.
124
d688e049 125 if (fVerbosity > 1) {
8a9d3e12 126 AliInfo(Form("Number of pt(t) bins: %d", fBPtT->GetNbins()));
127 for (Int_t i=1; i<=fBPtT->GetNbins(); i++) {
128 AliInfo(Form("pt(t) bin %d, %f to %f", i, fBPtT->GetBinLowEdge(i), fBPtT->GetBinUpEdge(i)));
129 }
d688e049 130 AliInfo(Form("Number of pt(a) bins: %d", fBPtA->GetNbins()));
131 for (Int_t i=1; i<=fBPtA->GetNbins(); i++) {
132 AliInfo(Form("pt(a) bin %d, %f to %f", i, fBPtA->GetBinLowEdge(i), fBPtA->GetBinUpEdge(i)));
133 }
134 }
48a61f36 135
d688e049 136 Int_t nPtAssc=fBPtA->GetNbins();
137 Int_t nPtTrig=fBPtT->GetNbins();
138 Int_t nCent=fBCent->GetNbins();
139 Int_t nZvtx=fBZvtx->GetNbins();
140 Double_t ptt[nPtTrig+1];
141 ptt[0] = fBPtT->GetBinLowEdge(1);
142 for (Int_t i=1; i<=nPtTrig; i++) {
143 ptt[i] = fBPtT->GetBinUpEdge(i);
144 }
145 Double_t pta[nPtAssc+1];
146 pta[0] = fBPtA->GetBinLowEdge(1);
147 for (Int_t i=1; i<=nPtAssc; i++) {
148 pta[i] = fBPtA->GetBinUpEdge(i);
149 }
150 Double_t cent[nCent+1];
151 cent[0] = fBCent->GetBinLowEdge(1);
152 for (Int_t i=1; i<=nCent; i++) {
153 cent[i] = fBCent->GetBinUpEdge(i);
154 }
155 Double_t zvtx[nZvtx+1];
156 zvtx[0] = fBZvtx->GetBinLowEdge(1);
157 for (Int_t i=1; i<=nZvtx; i++) {
158 zvtx[i] = fBZvtx->GetBinUpEdge(i);
159 }
48a61f36 160
8a9d3e12 161 fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", nZvtx, zvtx, nCent, cent);
8433ff41 162 fOutputList->Add(fHEvt);
48a61f36 163 fHTrk = new TH2F("fHTrk", "Track-level variables",
8433ff41 164 100, 0, TMath::TwoPi(), 100, -fEtaMax, +fEtaMax);
165 fOutputList->Add(fHTrk);
166
8433ff41 167 fHPtAss = new TH1F("fHPtAss","PtAssoc;P_{T} (GeV/c) [GeV/c]",nPtAssc,pta);
168 fOutputList->Add(fHPtAss);
169 fHPtTrg = new TH1F("fHPtTrg","PtTrig;P_{T} (GeV/c) [GeV/c]",nPtTrig,ptt);
170 fOutputList->Add(fHPtTrg);
beb13b6c 171 fHPtTrgEvt = new TH1F("fHPtTrgEvt","PtTrig;P_{T} (GeV/c) [GeV/c]",nPtTrig,ptt);
172 fHPtTrgNorm1S = new TH3F("fHPtTrgNorm1S","PtTrig;P_{T} (GeV/c) [GeV/c];centrality;z_{vtx}",nPtTrig,ptt,nCent,cent,nZvtx,zvtx);
173 fOutputList->Add(fHPtTrgNorm1S);
174 fHPtTrgNorm1M = (TH3*) fHPtTrgNorm1S->Clone("fHPtTrgNorm1M");
175 fOutputList->Add(fHPtTrgNorm1M);
176 fHPtTrgNorm2S = (TH3*) fHPtTrgNorm1S->Clone("fHPtTrgNorm2S");
177 fOutputList->Add(fHPtTrgNorm2S);
178 fHPtTrgNorm2M = (TH3*) fHPtTrgNorm1S->Clone("fHPtTrgNorm2M");
179 fOutputList->Add(fHPtTrgNorm2M);
180
8433ff41 181 fHCent = new TH1F("fHCent","Cent;bins",nCent,cent);
182 fOutputList->Add(fHCent);
183 fHZvtx = new TH1F("fHZvtx","Zvertex;bins",nZvtx,zvtx);
184 fOutputList->Add(fHZvtx);
185
186 fNbins = nPtTrig*nPtAssc*nCent*nZvtx;
8a9d3e12 187 fHSs = new TH2*[fNbins];
188 fHMs = new TH2*[fNbins];
189 fHPts = new TH1*[fNbins];
8433ff41 190
191 fIndex = new TFormula("GlobIndex","(t-1)*[0]*[1]*[2]+(z-1)*[0]*[1]+(x-1)*[0]+(y-1)+0*[4]");
192 fIndex->SetParameters(nPtTrig,nPtAssc,nZvtx,nCent);
193 fIndex->SetParNames("NTrigBins","NAssocBins", "NZvertexBins", "NCentBins");
beb13b6c 194 fOutputList->Add(fIndex);
8433ff41 195
8433ff41 196 Int_t count = 0;
197 for (Int_t c=1; c<=nCent; ++c) {
198 for (Int_t z=1; z<=nZvtx; ++z) {
199 for (Int_t t=1; t<=nPtTrig; ++t) {
200 for (Int_t a=1; a<=nPtAssc; ++a) {
8a9d3e12 201 fHSs[count] = 0;
202 fHMs[count] = 0;
203 fHPts[count] = 0;
8433ff41 204 if (a>t) {
205 ++count;
206 continue;
207 }
8a9d3e12 208 if (t==1 && a==1) {
209 TString title(Form("cen=%d (%.1f to %.1f), zVtx=%d (%.1f to %.1f)",
210 c, fBCent->GetBinLowEdge(c), fBCent->GetBinUpEdge(c),
211 z, fBZvtx->GetBinLowEdge(z), fBZvtx->GetBinUpEdge(z)));
212 fHPts[count] = new TH1F(Form("hPt%d",count), Form("Ptdist %s;p_{T} (GeV/c)",title.Data()), 200,0,20);
213 fOutputList->Add(fHPts[count]);
214 }
215 TString title(Form("cen=%d (%.1f to %.1f), zVtx=%d (%.1f to %.1f), trig=%d (%.1f to %.1f), assoc=%d (%.1f to %.1f)",
216 c, fBCent->GetBinLowEdge(c), fBCent->GetBinUpEdge(c),
217 z, fBZvtx->GetBinLowEdge(z), fBZvtx->GetBinUpEdge(z),
218 t, fBPtT->GetBinLowEdge(t), fBPtT->GetBinUpEdge(t),
219 a, fBPtA->GetBinLowEdge(a), fBPtA->GetBinUpEdge(a)));
220 fHSs[count] = new TH2F(Form("hS%d",count), Form("Signal %s;#Delta #varphi;#Delta #eta",title.Data()),
d688e049 221 fNBdphi,-0.5*TMath::Pi(),1.5*TMath::Pi(),fNBdeta,-2*fEtaMax,2*fEtaMax);
8a9d3e12 222 fHMs[count] = new TH2F(Form("hM%d",count), Form("Mixed %s;#Delta #varphi;#Delta #eta",title.Data()),
d688e049 223 fNBdphi,-0.5*TMath::Pi(),1.5*TMath::Pi(),fNBdeta,-2*fEtaMax,2*fEtaMax);
8433ff41 224 fOutputList->Add(fHSs[count]);
225 fOutputList->Add(fHMs[count]);
8a9d3e12 226 Int_t tcount = (Int_t)fIndex->Eval(t,a,z,c);
8433ff41 227 if (fVerbosity>5)
8a9d3e12 228 cout << count << " " << tcount << ": " << title << endl;
229 if (count != tcount)
230 AliFatal(Form("Index mismatch: %d %d", count, tcount));
8433ff41 231 ++count;
232 }
233 }
234 }
235 }
48a61f36 236}
237
beb13b6c 238//________________________________________________________________________
239void AliDhcTask::SetAnaMode(Int_t iAna=kHH)
240{
241 if (iAna==kHH) {
242 SetFillMuons(kFALSE);
243 fEtaTLo = -1.0;
244 fEtaTHi = 1.0;
245 fEtaALo = -1.0;
246 fEtaAHi = 1.0;
247 } else if (iAna==kMuH) {
248 SetFillMuons(kTRUE);
249 fEtaTLo = -5.0;
250 fEtaTHi = -2.0;
251 fEtaALo = -1.0;
252 fEtaAHi = 1.0;
253 } else if (iAna==kHMu) {
254 SetFillMuons(kTRUE);
255 fEtaTLo = -1.0;
256 fEtaTHi = 1.0;
257 fEtaALo = -5.0;
258 fEtaAHi = -2.0;
259 } else if (iAna==kMuMu) {
260 SetFillMuons(kTRUE);
261 fEtaTLo = -5.0;
262 fEtaTHi = -2.0;
263 fEtaALo = -5.0;
264 fEtaAHi = -2.0;
265 } else if (iAna==kPSide) {
266 SetFillMuons(kFALSE);
267 fEtaTLo = -1.0;
268 fEtaTHi = -0.465;
269 fEtaALo = -1.0;
270 fEtaAHi = -0.465;
271 } else if (iAna==kASide) {
272 SetFillMuons(kFALSE);
273 fEtaTLo = -0.465;
274 fEtaTHi = 1.0;
275 fEtaALo = -0.465;
276 fEtaAHi = 1.0;
277 } else {
278 AliInfo(Form("Unrecognized analysis option: %d", iAna));
279 }
280}
281
48a61f36 282//________________________________________________________________________
283void AliDhcTask::InitEventMixer()
284{
285 // The effective pool size in events is set by trackDepth, so more
286 // low-mult events are required to maintain the threshold than
287 // high-mult events. Centrality pools are indep. of data histogram
288 // binning, no need to match.
289
48a61f36 290 // Centrality pools
d688e049 291 Int_t nCentBins=fMixBCent->GetNbins();
292 Double_t centBins[nCentBins+1];
293 centBins[0] = fMixBCent->GetBinLowEdge(1);
294 for (Int_t i=1; i<=nCentBins; i++) {
295 centBins[i] = fMixBCent->GetBinUpEdge(i);
296 }
48a61f36 297
298 // Z-vertex pools
d688e049 299 Int_t nZvtxBins=fMixBZvtx->GetNbins();
300 Double_t zvtxbin[nZvtxBins+1];
301 zvtxbin[0] = fMixBZvtx->GetBinLowEdge(1);
302 for (Int_t i=1; i<=nZvtxBins; i++) {
303 zvtxbin[i] = fMixBZvtx->GetBinUpEdge(i);
304 }
8433ff41 305
306 fPoolMgr = new AliEvtPoolManager();
d688e049 307 fPoolMgr->SetTargetTrackDepth(fTrackDepth);
8433ff41 308 if (fVerbosity>4)
309 fPoolMgr->SetDebug(1);
d688e049 310 fPoolMgr->InitEventPools(fPoolSize, nCentBins, centBins, nZvtxBins, zvtxbin);
48a61f36 311}
312
313//________________________________________________________________________
314void AliDhcTask::UserExec(Option_t *)
315{
316 // Main loop, called for each event.
48a61f36 317 static int iEvent = -1; ++iEvent;
318
319 if (fVerbosity>2) {
320 if (iEvent % 10 == 0)
321 cout << iEvent << endl;
322 }
323
324 Int_t dType = -1; // Will be set to kESD or kAOD.
f6701c6e 325 MiniEvent* sTracks = new MiniEvent(0); // deleted by pool mgr.
48a61f36 326
327 LoadBranches();
328
329 // Get event pointers, check for signs of life
330 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
331 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
332 if (fESD)
333 dType = kESD;
334 else if (fAOD)
335 dType = kAOD;
336 else {
337 AliError("Neither fESD nor fAOD available");
338 return;
339 }
340
d688e049 341 Bool_t mcgen = 0;
342 if (fTracksName.Contains("Gen"))
343 mcgen = 1;
344
48a61f36 345 // Centrality, vertex, other event variables...
d688e049 346 if (mcgen) {
347 fZVertex = 0;
348 TList *list = InputEvent()->GetList();
beb13b6c 349 TClonesArray *tcaTracks = 0;
350 if (list)
351 tcaTracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracksName));
d688e049 352 if (tcaTracks)
353 fCentrality = tcaTracks->GetEntries();
354 } else {
355 if (dType == kESD) {
356 if (!VertexOk(fESD)) {
357 if (fVerbosity > 1)
358 AliInfo(Form("Event REJECTED (ESD vertex not OK). z = %.1f", fZVertex));
359 return;
360 }
361 const AliESDVertex* vertex = fESD->GetPrimaryVertex();
362 fZVertex = vertex->GetZ();
363 if(fESD->GetCentrality()) {
364 fCentrality =
365 fESD->GetCentrality()->GetCentralityPercentile(fCentMethod);
366 }
48a61f36 367 }
d688e049 368 if (dType == kAOD) {
369 const AliAODVertex* vertex = fAOD->GetPrimaryVertex();
370 fZVertex = vertex->GetZ();
371 if (!VertexOk(fAOD)) {
372 if (fVerbosity > 1)
373 Info("Exec()", "Event REJECTED (AOD vertex not OK). z = %.1f", fZVertex);
374 return;
375 }
376 const AliCentrality *aodCent = fAOD->GetHeader()->GetCentralityP();
377 if (aodCent) {
378 fCentrality = aodCent->GetCentralityPercentile(fCentMethod);
379 }
48a61f36 380 }
381 }
d688e049 382
48a61f36 383 // Fill Event histogram
384 fHEvt->Fill(fZVertex, fCentrality);
d688e049 385 if (fCentrality > fBCent->GetXmax() || fCentrality < fBCent->GetXmin()) {
9370528f 386 if (fVerbosity > 1)
387 AliInfo(Form("Event REJECTED (centrality out of range). fCentrality = %.1f", fCentrality));
48a61f36 388 return;
389 }
390
391 // Get array of selected tracks
392 if (dType == kESD) {
f6701c6e 393 GetESDTracks(sTracks);
48a61f36 394 }
395 if (dType == kAOD) {
f6701c6e 396 GetAODTracks(sTracks);
48a61f36 397 }
398
399 // Get pool containing tracks from other events like this one
8433ff41 400 AliEvtPool* pool = fPoolMgr->GetEventPool(fCentrality, fZVertex);
48a61f36 401 if (!pool) {
c71ff6f4 402 AliWarning(Form("No pool found. Centrality %f, ZVertex %f",
403 fCentrality, fZVertex));
48a61f36 404 return;
405 }
406
407 if (!pool->IsReady()) {
408 pool->UpdatePool(sTracks);
409 return;
410 }
411
412 if (pool->IsFirstReady()) {
413 // recover events missed before the pool is ready
414 Int_t nEvs = pool->GetCurrentNEvents();
415 if (nEvs>1) {
416 for (Int_t i=0; i<nEvs; ++i) {
417 MiniEvent* evI = pool->GetEvent(i);
418 for (Int_t j=0; j<nEvs; ++j) {
419 MiniEvent* evJ = pool->GetEvent(j);
420 if (i==j) {
421 Correlate(*evI, *evJ, kSameEvt);
422 } else {
d688e049 423 Correlate(*evI, *evJ, kDiffEvt);
48a61f36 424 }
425 }
426 }
427 }
428 } else { /* standard case: same event, then mix*/
429 Correlate(*sTracks, *sTracks, kSameEvt);
430 Int_t nMix = pool->GetCurrentNEvents();
431 for (Int_t jMix=0; jMix<nMix; ++jMix) {
432 MiniEvent* bgTracks = pool->GetEvent(jMix);
d688e049 433 Correlate(*sTracks, *bgTracks, kDiffEvt);
48a61f36 434 }
435 }
436
48a61f36 437 pool->UpdatePool(sTracks);
438 PostData(1, fOutputList);
48a61f36 439}
440
441//________________________________________________________________________
f6701c6e 442void AliDhcTask::GetESDTracks(MiniEvent* miniEvt)
48a61f36 443{
444 // Loop twice: 1. Count sel. tracks. 2. Fill vector.
445
d688e049 446 if (fTracksName.IsNull()) {
447 const AliESDVertex *vtxSPD = fESD->GetPrimaryVertexSPD();
448 if (!vtxSPD)
449 return;
450
451 Int_t nTrax = fESD->GetNumberOfTracks();
452 if (fVerbosity > 2)
453 AliInfo(Form("%d tracks in event",nTrax));
454
455 // Loop 1.
456 Int_t nSelTrax = 0;
457 TObjArray arr(nTrax);
458 arr.SetOwner(1);
459
460 for (Int_t i = 0; i < nTrax; ++i) {
461 AliESDtrack* esdtrack = fESD->GetTrack(i);
462 if (!esdtrack) {
463 AliError(Form("Couldn't get ESD track %d\n", i));
464 continue;
465 }
8a9d3e12 466 Bool_t trkOK = fEsdTPCOnly->AcceptTrack(esdtrack);
d688e049 467 if (!trkOK)
468 continue;
469 Double_t pt = esdtrack->Pt();
470 Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
471 if (!ptOK)
472 continue;
473 Double_t eta = esdtrack->Eta();
474 if (TMath::Abs(eta) > fEtaMax)
475 continue;
476
477 // create a tpc only track
478 AliESDtrack *newtrack = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
479 if(!newtrack)
480 continue;
481 if (newtrack->Pt()<=0) {
482 delete newtrack;
483 continue;
484 }
8433ff41 485
d688e049 486 AliExternalTrackParam exParam;
487 Bool_t relate = newtrack->RelateToVertexTPC(vtxSPD,fESD->GetMagneticField(),kVeryBig,&exParam);
488 if (!relate) {
489 delete newtrack;
490 continue;
491 }
8433ff41 492
d688e049 493 // set the constraint parameters to the track
494 newtrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
8433ff41 495
d688e049 496 pt = newtrack->Pt();
497 ptOK = pt >= fPtMin && pt < fPtMax;
498 if (!ptOK) {
499 delete newtrack;
500 continue;
501 }
502 eta = esdtrack->Eta();
503 if (TMath::Abs(eta) > fEtaMax) {
504 delete newtrack;
505 continue;
506 }
507 arr.Add(newtrack);
508 nSelTrax++;
8433ff41 509 }
beb13b6c 510
511 for(Int_t itrack = 0; itrack < nSelTrax; itrack++) {
512 AliVTrack *esdtrack = static_cast<AliESDtrack*>(arr.At(itrack));
513 if(!esdtrack) {
514 AliError(Form("ERROR: Could not retrieve esdtrack %d",itrack));
515 continue;
516 }
517 Double_t pt = esdtrack->Pt();
518 Double_t eta = esdtrack->Eta();
519 Double_t phi = esdtrack->Phi();
520 Int_t sign = esdtrack->Charge() > 0 ? 1 : -1;
521 miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
522 }
b673a083 523 } else {
524 TList *list = InputEvent()->GetList();
525 TClonesArray *tcaTracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracksName));
f018e7c8 526
b673a083 527 if(!tcaTracks){
528 AliError("Ptr to tcaTracks zero");
529 return;
530 }
f018e7c8 531
b673a083 532 const Int_t ntracks = tcaTracks->GetEntries();
533 if (miniEvt)
534 miniEvt->reserve(ntracks);
535 else {
536 AliError("Ptr to miniEvt zero");
537 return;
538 }
48a61f36 539
b673a083 540 for (Int_t itrack = 0; itrack < ntracks; itrack++) {
541 AliVTrack *vtrack = static_cast<AliVTrack*>(tcaTracks->At(itrack));
542 if (!vtrack) {
543 AliError(Form("ERROR: Could not retrieve track %d",itrack));
544 continue;
545 }
546 Double_t pt = vtrack->Pt();
547 Double_t eta = vtrack->Eta();
548 Double_t phi = vtrack->Phi();
549 Int_t sign = vtrack->Charge() > 0 ? 1 : -1;
550 miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
48a61f36 551 }
48a61f36 552 }
beb13b6c 553
554 if (fFillMuons) {
beb13b6c 555 // count good muons
556 Int_t nGoodMuons = 0;
557 for (Int_t iMu = 0; iMu<fESD->GetNumberOfMuonTracks(); iMu++) {
558 AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iMu);
b673a083 559 if (muonTrack) {
beb13b6c 560 if (IsGoodMUONtrack(*muonTrack)) nGoodMuons++;
561 }
562 }
563 miniEvt->reserve(miniEvt->size()+nGoodMuons);
564 // fill them into the mini event
565 for (Int_t iMu = 0; iMu<fESD->GetNumberOfMuonTracks(); iMu++) {
566 AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iMu);
b673a083 567 if (muonTrack) {
568 if (!IsGoodMUONtrack(*muonTrack))
569 continue;
570 Double_t ptMu = muonTrack->Pt();
571 Double_t etaMu = muonTrack->Eta();
572 Double_t phiMu = muonTrack->Phi();
573 Double_t signMu = muonTrack->Charge() > 0 ? 1 : -1;
beb13b6c 574 miniEvt->push_back(AliMiniTrack(ptMu, etaMu, phiMu, signMu));
575 }
576 }
577 }
48a61f36 578}
579
580//________________________________________________________________________
f6701c6e 581void AliDhcTask::GetAODTracks(MiniEvent* miniEvt)
48a61f36 582{
583 // Loop twice: 1. Count sel. tracks. 2. Fill vector.
584
b673a083 585 if (fTracksName.IsNull()) {
586 Int_t nTrax = fAOD->GetNumberOfTracks();
587 Int_t nSelTrax = 0;
48a61f36 588
b673a083 589 if (fVerbosity > 2)
590 AliInfo(Form("%d tracks in event",nTrax));
48a61f36 591
b673a083 592 // Loop 1.
593 for (Int_t i = 0; i < nTrax; ++i) {
594 AliAODTrack* aodtrack = fAOD->GetTrack(i);
595 if (!aodtrack) {
596 AliError(Form("Couldn't get AOD track %d\n", i));
597 continue;
598 }
599 // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
600 UInt_t tpcOnly = 1 << 7;
601 Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly);
602 if (!trkOK)
603 continue;
604 Double_t pt = aodtrack->Pt();
605 Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
606 if (!ptOK)
607 continue;
608 Double_t eta = aodtrack->Eta();
609 if (TMath::Abs(eta) > fEtaMax)
610 continue;
611 nSelTrax++;
48a61f36 612 }
48a61f36 613
b673a083 614 if (miniEvt)
615 miniEvt->reserve(nSelTrax);
616 else {
617 AliError("!miniEvt");
618 return;
48a61f36 619 }
b673a083 620
621 // Loop 2.
622 for (Int_t i = 0; i < nTrax; ++i) {
623 AliAODTrack* aodtrack = fAOD->GetTrack(i);
624 if (!aodtrack) {
625 AliError(Form("Couldn't get AOD track %d\n", i));
626 continue;
627 }
48a61f36 628
b673a083 629 // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
630 UInt_t tpcOnly = 1 << 7;
631 Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly);
632 if (!trkOK)
633 continue;
634 Double_t pt = aodtrack->Pt();
635 Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
636 if (!ptOK)
637 continue;
638 Double_t eta = aodtrack->Eta();
639 if (TMath::Abs(eta) > fEtaMax)
640 continue;
641
642 Double_t phi = aodtrack->Phi();
643 Int_t sign = aodtrack->Charge() > 0 ? 1 : -1;
644 miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
645 }
646 } else {
647 TList *list = InputEvent()->GetList();
648 TClonesArray *tcaTracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracksName));
649
650 if (!tcaTracks){
651 AliError("Ptr to tcaTracks zero");
652 return;
653 }
654
655 const Int_t ntracks = tcaTracks->GetEntries();
656 if (miniEvt)
657 miniEvt->reserve(ntracks);
658 else {
659 AliError("Ptr to miniEvt zero");
660 return;
661 }
662
663 for (Int_t itrack = 0; itrack < ntracks; itrack++) {
664 AliVTrack *vtrack = static_cast<AliVTrack*>(tcaTracks->At(itrack));
665 if (!vtrack) {
666 AliError(Form("ERROR: Could not retrieve vtrack %d",itrack));
667 continue;
668 }
669 Double_t pt = vtrack->Pt();
670 Double_t eta = vtrack->Eta();
671 Double_t phi = vtrack->Phi();
672 Int_t sign = vtrack->Charge() > 0 ? 1 : -1;
673 miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
674 }
675 }
48a61f36 676
b673a083 677 if (fFillMuons) {
678 // count good muons
679 Int_t nGoodMuons = 0;
680 for (Int_t iMu = 0; iMu<fAOD->GetNumberOfTracks(); iMu++) {
681 AliAODTrack* muonTrack = fAOD->GetTrack(iMu);
682 if(muonTrack) {
683 if (IsGoodMUONtrack(*muonTrack))
684 nGoodMuons++;
685 }
686 }
687 miniEvt->reserve(miniEvt->size()+nGoodMuons);
688 // fill them into the mini event
689 for (Int_t iMu = 0; iMu<fAOD->GetNumberOfTracks(); iMu++) {
690 AliAODTrack* muonTrack = fAOD->GetTrack(iMu);
691 if (muonTrack) {
692 if (!IsGoodMUONtrack(*muonTrack))
693 continue;
694 Double_t ptMu = muonTrack->Pt();
695 Double_t etaMu = muonTrack->Eta();
696 Double_t phiMu = muonTrack->Phi();
697 Double_t signMu = muonTrack->Charge() > 0 ? 1 : -1;
698 miniEvt->push_back(AliMiniTrack(ptMu, etaMu, phiMu, signMu));
699 }
700 }
48a61f36 701 }
48a61f36 702}
703
704//________________________________________________________________________
705Double_t AliDhcTask::DeltaPhi(Double_t phia, Double_t phib,
706 Double_t rangeMin, Double_t rangeMax) const
707{
708 Double_t dphi = -999;
709 Double_t pi = TMath::Pi();
710
711 if (phia < 0) phia += 2*pi;
712 else if (phia > 2*pi) phia -= 2*pi;
713 if (phib < 0) phib += 2*pi;
714 else if (phib > 2*pi) phib -= 2*pi;
715 dphi = phib - phia;
716 if (dphi < rangeMin) dphi += 2*pi;
717 else if (dphi > rangeMax) dphi -= 2*pi;
718
719 return dphi;
720}
721
722//________________________________________________________________________
d688e049 723Int_t AliDhcTask::Correlate(const MiniEvent &evt1, const MiniEvent &evt2, Int_t pairing)
48a61f36 724{
725 // Triggered angular correlations. If pairing is kSameEvt, particles
726 // within evt1 are correlated. If kDiffEvt, correlate triggers from
727 // evt1 with partners from evt2.
8433ff41 728
729 Int_t cbin = fHCent->FindBin(fCentrality);
730 if (fHCent->IsBinOverflow(cbin) ||
731 fHCent->IsBinUnderflow(cbin))
732 return 0;
733
734 Int_t zbin = fHZvtx->FindBin(fZVertex);
735 if (fHZvtx->IsBinOverflow(zbin) ||
736 fHZvtx->IsBinUnderflow(zbin))
737 return 0;
738
48a61f36 739 Int_t iMax = evt1.size();
740 Int_t jMax = evt2.size();
741
8433ff41 742 TH2 **hist = fHMs;
743 if (pairing == kSameEvt) {
744 hist = fHSs;
b673a083 745 fHCent->Fill(fCentrality);
746 fHZvtx->Fill(fZVertex);
8433ff41 747 }
748
749 Int_t nZvtx = fHZvtx->GetNbinsX();
750 Int_t nPtTrig = fHPtTrg->GetNbinsX();
751 Int_t nPtAssc = fHPtAss->GetNbinsX();
752
753 Int_t globIndex = (cbin-1)*nZvtx*nPtTrig*nPtAssc+(zbin-1)*nPtTrig*nPtAssc;
48a61f36 754
8a9d3e12 755 Int_t ptindex = (Int_t)fIndex->Eval(1,1,zbin,cbin);
756
beb13b6c 757 fHPtTrgEvt->Reset();
383b3bca 758 for (Int_t i=0; i<iMax; ++i) {
759 const AliMiniTrack &a(evt1.at(i));
760 Float_t pta = a.Pt();
beb13b6c 761 fHPtTrgEvt->Fill(pta);
383b3bca 762 if (pairing == kSameEvt) {
8a9d3e12 763 fHPts[ptindex]->Fill(pta);
d688e049 764 }
765 }
73a051c1 766
beb13b6c 767 Bool_t bCountTrg; // only count the trigger if an associated particle is found
48a61f36 768
beb13b6c 769 for (Int_t i=0; i<iMax; ++i) {
48a61f36 770 // Trigger particles
8433ff41 771 const AliMiniTrack &a(evt1.at(i));
772
773 Float_t pta = a.Pt();
beb13b6c 774 Float_t etaa = a.Eta();
7e990b20 775 Float_t phia = a.Phi();
8433ff41 776 Int_t abin = fHPtTrg->FindBin(pta);
beb13b6c 777 if (fHPtTrg->IsBinOverflow(abin) || fHPtTrg->IsBinUnderflow(abin))
778 continue;
779
780 if (etaa<fEtaTLo || etaa>fEtaTHi)
8433ff41 781 continue;
48a61f36 782
beb13b6c 783 // efficiency weighting
beb13b6c 784 Double_t effWtT = 1.0;
785 if (fHEffT) {
786 // trigger particle
7e990b20 787 const Int_t nEffDimT = fHEffT->GetNdimensions();
788 Int_t effBinT[nEffDimT];
789 effBinT[0] = fHEffT->GetAxis(0)->FindBin(etaa);
790 effBinT[1] = fHEffT->GetAxis(1)->FindBin(pta);
791 effBinT[2] = fHEffT->GetAxis(2)->FindBin(fCentrality);
792 effBinT[3] = fHEffT->GetAxis(3)->FindBin(fZVertex);
793 if (nEffDimT>4) {
794 effBinT[4] = fHEffT->GetAxis(4)->FindBin(phia);
795 }
796 effWtT = fHEffT->GetBinContent(effBinT);
beb13b6c 797 }
798
48a61f36 799 if (pairing == kSameEvt) {
7e990b20 800 fHTrk->Fill(phia,etaa);
beb13b6c 801 fHPtTrgNorm1S->Fill(pta,fCentrality,fZVertex,effWtT);
802 } else {
803 fHPtTrgNorm1M->Fill(pta,fCentrality,fZVertex,effWtT);
48a61f36 804 }
beb13b6c 805 bCountTrg = kFALSE;
48a61f36 806
8433ff41 807 for (Int_t j=0; j<jMax; ++j) {
48a61f36 808 // Associated particles
809 if (pairing == kSameEvt && i==j)
beb13b6c 810 continue;
48a61f36 811
8433ff41 812 const AliMiniTrack &b(evt2.at(j));
48a61f36 813
48a61f36 814 Float_t ptb = b.Pt();
beb13b6c 815 Float_t etab = b.Eta();
7e990b20 816 Float_t phib = b.Phi();
beb13b6c 817 if (pta < ptb)
818 continue;
48a61f36 819
73a051c1 820 Int_t bbin = fHPtAss->FindBin(ptb);
beb13b6c 821 if (fHPtAss->IsBinOverflow(bbin) || fHPtAss->IsBinUnderflow(bbin))
822 continue;
823
824 if (etab<fEtaALo || etab>fEtaAHi)
825 continue;
8433ff41 826
7e990b20 827 Float_t dphi = DeltaPhi(phia, phib);
beb13b6c 828 Float_t deta = etaa - etab;
48a61f36 829
8433ff41 830 Int_t index = globIndex+(abin-1)*nPtAssc+(bbin-1);
8a9d3e12 831 Double_t weight = 1.0;
beb13b6c 832 // number of trigger weight event-by-event (CMS method)
833 if (fDoWeights) {
834 Double_t nTrgs = fHPtTrgEvt->GetBinContent(abin);
835 if (nTrgs==0.0) {
836 AliError(Form("Trying to work with number of triggers weight = %f",nTrgs));
837 return 0;
838 }
839 weight *= 1./nTrgs;
840 }
8a9d3e12 841
beb13b6c 842 // efficiency weighting
843 if (fHEffT) {
844 // trigger particle
845 weight *= effWtT;
73a051c1 846 }
beb13b6c 847 if (fHEffA) {
848 // associated particle
7e990b20 849 const Int_t nEffDimA = fHEffA->GetNdimensions();
850 Int_t effBinA[nEffDimA];
851 effBinA[0] = fHEffA->GetAxis(0)->FindBin(etab);
852 effBinA[1] = fHEffA->GetAxis(1)->FindBin(ptb);
853 effBinA[2] = fHEffA->GetAxis(2)->FindBin(fCentrality);
854 effBinA[3] = fHEffA->GetAxis(3)->FindBin(fZVertex);
855 if (nEffDimA>4) {
856 effBinA[4] = fHEffA->GetAxis(4)->FindBin(phib);
857 }
858 weight *= fHEffA->GetBinContent(effBinA);
beb13b6c 859 }
beb13b6c 860 hist[index]->Fill(dphi,deta,weight);
861 bCountTrg = kTRUE;
8433ff41 862
863 if (pairing == kSameEvt) {
b673a083 864 fHPtAss->Fill(ptb);
beb13b6c 865 }
866 }
867 if (bCountTrg) {
868 if (pairing == kSameEvt) {
869 fHPtTrgNorm2S->Fill(pta,fCentrality,fZVertex,effWtT);
870 } else {
871 fHPtTrgNorm2M->Fill(pta,fCentrality,fZVertex,effWtT);
8433ff41 872 }
48a61f36 873 }
874 }
875
beb13b6c 876 return 1;
48a61f36 877}
878
879//________________________________________________________________________
880void AliDhcTask::Terminate(Option_t *)
881{
882 // Draw result to the screen
883 // Called once at the end of the query
48a61f36 884}
885
886//________________________________________________________________________
887Bool_t AliDhcTask::VertexOk(TObject* obj) const
888{
889 // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
890
891 Int_t nContributors = 0;
892 Double_t zVertex = 999;
893 TString name("");
894
895 if (obj->InheritsFrom("AliESDEvent")) {
896 AliESDEvent* esdevt = (AliESDEvent*) obj;
897 const AliESDVertex* vtx = esdevt->GetPrimaryVertex();
898 if (!vtx)
899 return 0;
900 nContributors = vtx->GetNContributors();
901 zVertex = vtx->GetZ();
902 name = vtx->GetName();
903 }
904 else if (obj->InheritsFrom("AliAODEvent")) {
905 AliAODEvent* aodevt = (AliAODEvent*) obj;
906 if (aodevt->GetNumberOfVertices() < 1)
907 return 0;
908 const AliAODVertex* vtx = aodevt->GetPrimaryVertex();
909 nContributors = vtx->GetNContributors();
910 zVertex = vtx->GetZ();
911 name = vtx->GetName();
912 }
913
914 // Reject if TPC-only vertex
8433ff41 915 if (name.CompareTo("TPCVertex")==0)
48a61f36 916 return kFALSE;
917
918 // Check # contributors and range...
919 if( nContributors < 1 || TMath::Abs(zVertex) > fZVtxMax ) {
920 return kFALSE;
921 }
922
923 return kTRUE;
924}
beb13b6c 925
926//________________________________________________________________________
927Bool_t AliDhcTask::IsGoodMUONtrack(AliESDMuonTrack &track)
928{
b673a083 929 // Applying track cuts for MUON tracks
930
931 if (!track.ContainTrackerData())
932 return kFALSE;
933 if (!track.ContainTriggerData())
934 return kFALSE;
935 Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
936 if ((thetaTrackAbsEnd < 2.) || (thetaTrackAbsEnd > 10.))
937 return kFALSE;
938 Double_t eta = track.Eta();
939 if ((eta < -4.) || (eta > -2.5))
940 return kFALSE;
941 if (track.GetMatchTrigger() < 0.5)
942 return kFALSE;
943 return kTRUE;
944}
945
946//________________________________________________________________________
947Bool_t AliDhcTask::IsGoodMUONtrack(AliAODTrack &track)
948{
949 // Applying track cuts for MUON tracks
950
951 if (!track.IsMuonTrack())
952 return kFALSE;
953 Double_t dThetaAbs = TMath::ATan(track.GetRAtAbsorberEnd()/505.)* TMath::RadToDeg();
954 if ((dThetaAbs<2.) || (dThetaAbs>10.))
955 return kFALSE;
956 Double_t dEta = track.Eta();
957 if ((dEta<-4.) || (dEta>2.5))
958 return kFALSE;
959 if (track.GetMatchTrigger()<0.5)
960 return kFALSE;
961 return kTRUE;
962}
963
964//________________________________________________________________________
965AliDhcTask::~AliDhcTask()
966{
967 //Destructor
968 if (fPoolMgr)
969 delete fPoolMgr;
beb13b6c 970}