update DHC task (Constantin Loizides <cloizides@lbl.gov>)
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / FourierDecomposition / AliDhcTask.cxx
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"
7 #include "TFormula.h"
8 #include "TH1.h"
9 #include "TH2.h"
10 #include "TH3.h"
11 #include "TAxis.h"
12 #include "TProfile2D.h"
13 #include "TROOT.h"
14 #include "TTree.h"
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"
24 #include "AliESDMuonTrack.h"
25 #include "AliPool.h"
26 #include "AliVParticle.h"
27
28 using std::cout;
29 using std::endl;
30
31 ClassImp(AliDhcTask)
32
33
34 //________________________________________________________________________
35 AliDhcTask::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
52 //________________________________________________________________________
53 AliDhcTask::AliDhcTask(const char *name) 
54 : AliAnalysisTaskSE(name), fVerbosity(0), fEtaMax(1), fZVtxMax(10), fPtMin(0.25), fPtMax(15),
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),
61   fCentrality(99), fZVertex(99), fEsdTPCOnly(0), fPoolMgr(0),
62   fCentMethod("V0M"), fNBdeta(20), fNBdphi(36),
63   fBPtT(0x0), fBPtA(0x0), fBCent(0x0), fBZvtx(0x0),
64   fMixBCent(0x0), fMixBZvtx(0x0),
65   fHEffT(0x0), fHEffA(0x0)
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
76   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,Tracks "
77                "AOD:header,tracks,vertices,";
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);
91 }
92
93 //________________________________________________________________________
94 void AliDhcTask::UserCreateOutputObjects()
95 {
96   // Create histograms
97   // Called once (per slave on PROOF!)
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   }
106
107   fOutputList = new TList();
108   fOutputList->SetOwner(1);
109
110   fEsdTPCOnly = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
111   //fEsdTPCOnly->SetMinNClustersTPC(70);
112   fEsdTPCOnly->SetMinNCrossedRowsTPC(70);
113   fEsdTPCOnly->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
114
115   BookHistos();
116   InitEventMixer(); 
117   PostData(1, fOutputList);
118 }
119
120 //________________________________________________________________________
121 void AliDhcTask::BookHistos()
122 {
123   // Book histograms.
124
125   if (fVerbosity > 1) {
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     }
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   }
135
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   }
160
161   fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", nZvtx, zvtx, nCent, cent);
162   fOutputList->Add(fHEvt);
163   fHTrk = new TH2F("fHTrk", "Track-level variables", 
164                    100, 0, TMath::TwoPi(), 100, -fEtaMax, +fEtaMax);
165   fOutputList->Add(fHTrk);
166   
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);
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
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;
187   fHSs   = new TH2*[fNbins];
188   fHMs   = new TH2*[fNbins];
189   fHPts  = new TH1*[fNbins];
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");
194   fOutputList->Add(fIndex);
195
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) {
201           fHSs[count]  = 0;
202           fHMs[count]  = 0;
203           fHPts[count] = 0;
204           if (a>t) {
205             ++count;
206             continue;
207           }
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()),
221                                  fNBdphi,-0.5*TMath::Pi(),1.5*TMath::Pi(),fNBdeta,-2*fEtaMax,2*fEtaMax);
222           fHMs[count] = new TH2F(Form("hM%d",count), Form("Mixed %s;#Delta #varphi;#Delta #eta",title.Data()),
223                                  fNBdphi,-0.5*TMath::Pi(),1.5*TMath::Pi(),fNBdeta,-2*fEtaMax,2*fEtaMax);
224           fOutputList->Add(fHSs[count]);
225           fOutputList->Add(fHMs[count]);
226           Int_t tcount = (Int_t)fIndex->Eval(t,a,z,c);
227           if (fVerbosity>5)
228             cout << count << " " << tcount << ": " << title << endl;
229           if (count != tcount)
230             AliFatal(Form("Index mismatch: %d %d", count, tcount));
231           ++count;
232         }
233       }
234     }
235   }
236 }
237
238 //________________________________________________________________________
239 void 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
282 //________________________________________________________________________
283 void 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
290   // Centrality pools
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   }
297  
298   // Z-vertex pools
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   }
305
306   fPoolMgr = new AliEvtPoolManager();
307   fPoolMgr->SetTargetTrackDepth(fTrackDepth);
308   if (fVerbosity>4)
309     fPoolMgr->SetDebug(1);
310   fPoolMgr->InitEventPools(fPoolSize, nCentBins, centBins, nZvtxBins, zvtxbin);
311 }
312
313 //________________________________________________________________________
314 void AliDhcTask::UserExec(Option_t *) 
315 {
316   // Main loop, called for each event.
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.
325   MiniEvent* sTracks = new MiniEvent(0); // deleted by pool mgr.
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
341   Bool_t mcgen = 0;
342   if (fTracksName.Contains("Gen"))
343     mcgen = 1;
344
345   // Centrality, vertex, other event variables...
346   if (mcgen) {
347     fZVertex = 0;
348     TList *list = InputEvent()->GetList();
349     TClonesArray *tcaTracks = 0;
350     if (list) 
351       tcaTracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracksName));
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       }
367     }
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       }
380     }
381   }
382
383   // Fill Event histogram
384   fHEvt->Fill(fZVertex, fCentrality);
385   if (fCentrality > fBCent->GetXmax() || fCentrality < fBCent->GetXmin()) {
386     if (fVerbosity > 1)
387       AliInfo(Form("Event REJECTED (centrality out of range). fCentrality = %.1f", fCentrality));
388     return;
389   }
390
391   // Get array of selected tracks
392   if (dType == kESD) {
393     GetESDTracks(sTracks);
394   }
395   if (dType == kAOD) {
396     GetAODTracks(sTracks);
397   }
398
399   // Get pool containing tracks from other events like this one
400   AliEvtPool* pool = fPoolMgr->GetEventPool(fCentrality, fZVertex);
401   if (!pool) {
402     AliWarning(Form("No pool found. Centrality %f, ZVertex %f", 
403                     fCentrality, fZVertex));
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 {
423             Correlate(*evI, *evJ, kDiffEvt);
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);
433       Correlate(*sTracks, *bgTracks, kDiffEvt);
434     }
435   }
436
437   pool->UpdatePool(sTracks);
438   PostData(1, fOutputList);
439 }
440
441 //________________________________________________________________________
442 void AliDhcTask::GetESDTracks(MiniEvent* miniEvt)
443 {
444   // Loop twice: 1. Count sel. tracks. 2. Fill vector.
445
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       }
466       Bool_t trkOK = fEsdTPCOnly->AcceptTrack(esdtrack);
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       }
485
486       AliExternalTrackParam exParam;
487       Bool_t relate = newtrack->RelateToVertexTPC(vtxSPD,fESD->GetMagneticField(),kVeryBig,&exParam);
488       if (!relate) {
489         delete newtrack;
490         continue;
491       }
492
493       // set the constraint parameters to the track
494       newtrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
495
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++;
509     }
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     }
523   } else {
524     TList *list = InputEvent()->GetList();
525     TClonesArray *tcaTracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracksName));
526
527     if(!tcaTracks){
528       AliError("Ptr to tcaTracks zero");
529       return;
530     }
531
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     }
539
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));
551     }
552   }
553
554   if (fFillMuons) {
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);
559       if (muonTrack) {
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);
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;
574         miniEvt->push_back(AliMiniTrack(ptMu, etaMu, phiMu, signMu));
575       }
576     }
577   }
578 }
579
580 //________________________________________________________________________
581 void AliDhcTask::GetAODTracks(MiniEvent* miniEvt)
582 {
583   // Loop twice: 1. Count sel. tracks. 2. Fill vector.
584
585   if (fTracksName.IsNull()) {
586     Int_t nTrax = fAOD->GetNumberOfTracks();
587     Int_t nSelTrax = 0;
588
589     if (fVerbosity > 2)
590       AliInfo(Form("%d tracks in event",nTrax));
591
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++;
612     }
613
614     if (miniEvt)
615       miniEvt->reserve(nSelTrax);
616     else {
617       AliError("!miniEvt");
618       return;
619     }
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       }
628     
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   }
676
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     }
701   }
702 }
703
704 //________________________________________________________________________
705 Double_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 //________________________________________________________________________
723 Int_t AliDhcTask::Correlate(const MiniEvent &evt1, const MiniEvent &evt2, Int_t pairing)
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.
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
739   Int_t iMax = evt1.size();
740   Int_t jMax = evt2.size();
741
742   TH2  **hist = fHMs;
743   if (pairing == kSameEvt) {
744     hist = fHSs;
745     fHCent->Fill(fCentrality);
746     fHZvtx->Fill(fZVertex);
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;
754
755   Int_t ptindex = (Int_t)fIndex->Eval(1,1,zbin,cbin);
756
757   fHPtTrgEvt->Reset();
758   for (Int_t i=0; i<iMax; ++i) {
759     const AliMiniTrack &a(evt1.at(i));
760     Float_t pta  = a.Pt();
761     fHPtTrgEvt->Fill(pta);
762     if (pairing == kSameEvt) {
763       fHPts[ptindex]->Fill(pta);
764     }
765   }
766
767   Bool_t bCountTrg; // only count the trigger if an associated particle is found
768
769   for (Int_t i=0; i<iMax; ++i) {
770     // Trigger particles
771     const AliMiniTrack &a(evt1.at(i));
772
773     Float_t pta  = a.Pt();
774     Float_t etaa = a.Eta();
775     Float_t phia = a.Phi();
776     Int_t abin = fHPtTrg->FindBin(pta);
777     if (fHPtTrg->IsBinOverflow(abin) || fHPtTrg->IsBinUnderflow(abin))
778       continue;
779
780     if (etaa<fEtaTLo || etaa>fEtaTHi)
781       continue;
782
783     // efficiency weighting
784     Double_t effWtT = 1.0;
785     if (fHEffT) {
786       // trigger particle
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);
797     }
798     
799     if (pairing == kSameEvt) {
800       fHTrk->Fill(phia,etaa);
801       fHPtTrgNorm1S->Fill(pta,fCentrality,fZVertex,effWtT);
802     } else {
803       fHPtTrgNorm1M->Fill(pta,fCentrality,fZVertex,effWtT);
804     }
805     bCountTrg = kFALSE;
806
807     for (Int_t j=0; j<jMax; ++j) {
808       // Associated particles
809       if (pairing == kSameEvt && i==j)
810         continue;
811
812       const AliMiniTrack &b(evt2.at(j));
813       
814       Float_t ptb  = b.Pt();
815       Float_t etab = b.Eta();
816       Float_t phib = b.Phi();
817       if (pta < ptb)
818             continue;
819
820       Int_t bbin = fHPtAss->FindBin(ptb);
821       if (fHPtAss->IsBinOverflow(bbin) || fHPtAss->IsBinUnderflow(bbin))
822         continue;
823
824       if (etab<fEtaALo || etab>fEtaAHi)
825         continue;
826
827       Float_t dphi = DeltaPhi(phia, phib);
828       Float_t deta = etaa - etab;
829
830       Int_t index = globIndex+(abin-1)*nPtAssc+(bbin-1);
831       Double_t weight = 1.0;
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       }
841
842       // efficiency weighting
843       if (fHEffT) {
844         // trigger particle
845         weight *= effWtT;
846       }
847       if (fHEffA) {
848         // associated particle
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);
859       }
860       hist[index]->Fill(dphi,deta,weight);
861       bCountTrg = kTRUE;
862
863       if (pairing == kSameEvt) {
864         fHPtAss->Fill(ptb);
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);
872       }
873     }
874   }
875
876   return 1;
877 }
878
879 //________________________________________________________________________
880 void AliDhcTask::Terminate(Option_t *) 
881 {
882   // Draw result to the screen
883   // Called once at the end of the query
884 }
885
886 //________________________________________________________________________
887 Bool_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
915   if (name.CompareTo("TPCVertex")==0)
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 }
925
926 //________________________________________________________________________
927 Bool_t AliDhcTask::IsGoodMUONtrack(AliESDMuonTrack &track)
928 {
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 //________________________________________________________________________
947 Bool_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 //________________________________________________________________________
965 AliDhcTask::~AliDhcTask()
966 {
967   //Destructor
968   if (fPoolMgr) 
969     delete fPoolMgr;
970 }