478033480563282ee33e6b130ff35bf199255711
[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 "TProfile2D.h"
11 #include "TROOT.h"
12 #include "TTree.h"
13 #include "AliAODEvent.h"
14 #include "AliAODInputHandler.h"
15 #include "AliAnalysisManager.h"
16 #include "AliAnalysisTask.h"
17 #include "AliCentrality.h"
18 #include "AliDhcTask.h"
19 #include "AliESDEvent.h"
20 #include "AliESDInputHandler.h"
21 #include "AliESDtrackCuts.h"
22 #include "AliPool.h"
23 #include "AliVParticle.h"
24
25 ClassImp(AliDhcTask)
26
27 //________________________________________________________________________
28 AliDhcTask::AliDhcTask(const char *name) 
29 : AliAnalysisTaskSE(name), fVerbosity(0), fEtaMax(1), fZVtxMax(10), fPtMin(0.25), fPtMax(15),
30   fESD(0), fAOD(0), fOutputList(0), fHistPt(0), fHEvt(0), fHTrk(0), fHPtAss(0), fHPtTrg(0),
31   fHCent(0), fHZvtx(0), fNbins(0), fHSs(0), fHMs(0),
32   fIndex(0), fMeanPtTrg(0), fMeanPtAss(0), fMean2PtTrg(0), fMean2PtAss(0),
33   fCentrality(99), fZVertex(99), fEsdTrackCutsTPCOnly(0), fPoolMgr(0)
34 {
35   // Constructor
36
37   // Define input and output slots here
38   // Input slot #0 works with a TChain
39   DefineInput(0, TChain::Class());
40   // Output slot #0 id reserved by the base class for AOD
41   // Output slot #1 writes into a TH1 container
42   DefineOutput(1, TList::Class());
43
44   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,Tracks "
45                "AOD:header,tracks,vertices,";
46 }
47
48 //________________________________________________________________________
49 void AliDhcTask::UserCreateOutputObjects()
50 {
51   // Create histograms
52   // Called once (per slave on PROOF!)
53
54   fOutputList = new TList();
55   fOutputList->SetOwner(1);
56
57   fEsdTrackCutsTPCOnly = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
58   //fEsdTrackCutsTPCOnly->SetMinNClustersTPC(70);
59   fEsdTrackCutsTPCOnly->SetMinNCrossedRowsTPC(70);
60   fEsdTrackCutsTPCOnly->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
61
62   BookHistos();
63   InitEventMixer(); 
64   PostData(1, fOutputList);
65 }
66
67 //________________________________________________________________________
68 void AliDhcTask::BookHistos()
69 {
70   // Book histograms.
71
72   Int_t nDeta=20, nPtAssc=12, nPtTrig=12, nCent=12, nDphi=36, nZvtx=8;
73
74   Double_t ptt[]  = {0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 8.0, 15};
75   Double_t pta[]  = {0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 8.0, 15};
76   Double_t cent[] = {0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 90};
77   Double_t zvtx[] = {-10, -6, -4, -2, 0, 2, 4, 6, 10};
78
79   // Event histo
80   fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", 30, -15, 15, 101, 0, 101);
81   fOutputList->Add(fHEvt);
82   // Track histo
83   fHTrk = new TH2F("fHTrk", "Track-level variables", 
84                    100, 0, TMath::TwoPi(), 100, -fEtaMax, +fEtaMax);
85   fOutputList->Add(fHTrk);
86   
87   // Left over from the tutorial :)
88   fHistPt = new TH1F("fHistPt", "P_{T} distribution", 200, 0., fPtMax);
89   fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
90   fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
91   fHistPt->SetMarkerStyle(kFullCircle);
92   fOutputList->Add(fHistPt);
93
94   fHPtAss = new TH1F("fHPtAss","PtAssoc;P_{T} (GeV/c) [GeV/c]",nPtAssc,pta);
95   fOutputList->Add(fHPtAss);
96   fHPtTrg = new TH1F("fHPtTrg","PtTrig;P_{T} (GeV/c) [GeV/c]",nPtTrig,ptt);
97   fOutputList->Add(fHPtTrg);
98   fHCent = new TH1F("fHCent","Cent;bins",nCent,cent);
99   fOutputList->Add(fHCent);
100   fHZvtx = new TH1F("fHZvtx","Zvertex;bins",nZvtx,zvtx);
101   fOutputList->Add(fHZvtx);
102
103   fNbins = nPtTrig*nPtAssc*nCent*nZvtx;
104   fHSs = new TH2*[fNbins];
105   fHMs = new TH2*[fNbins];
106
107   fIndex = new TFormula("GlobIndex","(t-1)*[0]*[1]*[2]+(z-1)*[0]*[1]+(x-1)*[0]+(y-1)+0*[4]");
108   fIndex->SetParameters(nPtTrig,nPtAssc,nZvtx,nCent);
109   fIndex->SetParNames("NTrigBins","NAssocBins", "NZvertexBins", "NCentBins");
110   //fOutputList->Add(fIndex);
111
112   fMeanPtTrg = new TProfile2D*[nPtTrig*nPtAssc];
113   fMeanPtAss = new TProfile2D*[nPtTrig*nPtAssc];
114   fMean2PtTrg = new TProfile2D*[nPtTrig*nPtAssc];
115   fMean2PtAss = new TProfile2D*[nPtTrig*nPtAssc];
116   for (Int_t c=1; c<=nCent; ++c) {
117     TString title(Form("cen=%d (%.1f)",c,fHCent->GetBinCenter(c)));
118     fMeanPtTrg[c-1]  = new TProfile2D(Form("hMeanPtTrgCen%d",c),title,nPtTrig,ptt,nPtAssc,pta);
119     fMeanPtAss[c-1]  = new TProfile2D(Form("hMeanPtAssCen%d",c),title,nPtTrig,ptt,nPtAssc,pta);
120     fMean2PtTrg[c-1] = new TProfile2D(Form("hMean2PtTrgCen%d",c),title,nPtTrig,ptt,nPtAssc,pta);
121     fMean2PtAss[c-1] = new TProfile2D(Form("hMean2PtAssCen%d",c),title,nPtTrig,ptt,nPtAssc,pta);
122     fOutputList->Add(fMeanPtTrg[c-1]);
123     fOutputList->Add(fMeanPtAss[c-1]);
124     fOutputList->Add(fMean2PtTrg[c-1]);
125     fOutputList->Add(fMean2PtAss[c-1]);
126   }
127
128   Int_t count = 0;
129   for (Int_t c=1; c<=nCent; ++c) {
130     for (Int_t z=1; z<=nZvtx; ++z) {
131       for (Int_t t=1; t<=nPtTrig; ++t) {
132         for (Int_t a=1; a<=nPtAssc; ++a) {
133           fHSs[count] = 0;
134           fHMs[count] = 0;
135           if (a>t) {
136             ++count;
137             continue;
138           }
139           TString title(Form("cen=%d (%.1f), zVtx=%d (%.1f), trig=%d (%.1f), assc=%d (%.1f)",
140                              c,fHCent->GetBinCenter(c), z,fHZvtx->GetBinCenter(z),
141                              t,fHPtTrg->GetBinCenter(t),a, fHPtAss->GetBinCenter(a)));
142           fHSs[count] = new TH2F(Form("hS%d",count), Form("Signal %s",title.Data()),
143                                  nDphi,-0.5*TMath::Pi(),1.5*TMath::Pi(),nDeta,-2*fEtaMax,2*fEtaMax);
144           fHMs[count] = new TH2F(Form("hM%d",count), Form("Signal %s",title.Data()),
145                                  nDphi,-0.5*TMath::Pi(),1.5*TMath::Pi(),nDeta,-2*fEtaMax,2*fEtaMax);
146           fOutputList->Add(fHSs[count]);
147           fOutputList->Add(fHMs[count]);
148           if (fVerbosity>5)
149             cout << count << " " << fIndex->Eval(t,a,z,c) << ": " << title << endl;
150           ++count;
151         }
152       }
153     }
154   }
155
156   return;
157 }
158
159 //________________________________________________________________________
160 void AliDhcTask::InitEventMixer()
161 {
162   // The effective pool size in events is set by trackDepth, so more
163   // low-mult events are required to maintain the threshold than
164   // high-mult events. Centrality pools are indep. of data histogram
165   // binning, no need to match.
166
167   Int_t trackDepth = 1000;   // # tracks to fill pool
168   Int_t poolsize   = 200;    // Maximum number of events
169
170   // Centrality pools
171   Int_t nCentBins  = 12;
172   Double_t centBins[] = {0,1,2,3,4,5,10,20,30,40,50,60,90.1};
173  
174   //Int_t nCentBins  = 1;
175   //Double_t centBins[] = {-1,100.1};
176  
177   // Z-vertex pools
178   Int_t nZvtxBins  = 8;
179   Double_t zvtxbin[] = {-10,-6,-4,-2,0,2,4,6,10};
180
181   fPoolMgr = new AliEvtPoolManager();
182   fPoolMgr->SetTargetTrackDepth(trackDepth);
183   if (fVerbosity>4)
184     fPoolMgr->SetDebug(1);
185   fPoolMgr->InitEventPools(poolsize, nCentBins, centBins, nZvtxBins, zvtxbin);
186   
187   return;
188 }
189
190 //________________________________________________________________________
191 void AliDhcTask::UserExec(Option_t *) 
192 {
193   // Main loop, called for each event.
194
195   static int iEvent = -1; ++iEvent;
196
197   if (fVerbosity>2) {
198     if (iEvent % 10 == 0) 
199       cout << iEvent << endl;
200   }
201
202   Int_t dType = -1;       // Will be set to kESD or kAOD.
203   MiniEvent* sTracks = new MiniEvent(0); // deleted by pool mgr.
204   Double_t centCL1 = -1;
205
206   LoadBranches();
207
208   // Get event pointers, check for signs of life
209   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
210   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
211   if (fESD)
212     dType = kESD;
213   else if (fAOD)
214     dType = kAOD;
215   else {
216     AliError("Neither fESD nor fAOD available");
217     return;
218   }
219
220   // Centrality, vertex, other event variables...
221   if (dType == kESD) {
222     if (!VertexOk(fESD)) {
223       if (fVerbosity > 1)
224         AliInfo(Form("Event REJECTED. z = %.1f", fZVertex));
225       return;
226     }
227     const AliESDVertex* vertex = fESD->GetPrimaryVertex();
228     fZVertex = vertex->GetZ();
229     if(fESD->GetCentrality()) {
230       fCentrality = 
231         fESD->GetCentrality()->GetCentralityPercentile("V0M");
232       centCL1 =
233         fESD->GetCentrality()->GetCentralityPercentile("CL1");
234     }
235   }
236   if (dType == kAOD) {
237     const AliAODVertex* vertex = fAOD->GetPrimaryVertex();
238     fZVertex = vertex->GetZ();
239     if (!VertexOk(fAOD)) {
240       if (fVerbosity > 1)
241         Info("Exec()", "Event REJECTED. z = %.1f", fZVertex);
242       return;
243     }
244     const AliCentrality *aodCent = fAOD->GetHeader()->GetCentralityP();
245     if (aodCent) {
246       fCentrality = aodCent->GetCentralityPercentile("V0M");
247       centCL1     = aodCent->GetCentralityPercentile("CL1");
248     }
249   }
250   
251   // Fill Event histogram
252   fHEvt->Fill(fZVertex, fCentrality);
253
254   if (fCentrality > 90. || fCentrality < 0) {
255     AliInfo(Form("Event REJECTED. fCentrality = %.1f", fCentrality));
256     return;
257   }
258
259   // Get array of selected tracks
260   if (dType == kESD) {
261     GetESDTracks(sTracks);
262   }
263   if (dType == kAOD) {
264     GetAODTracks(sTracks);
265   }
266
267   // Get pool containing tracks from other events like this one
268   AliEvtPool* pool = fPoolMgr->GetEventPool(fCentrality, fZVertex);
269   if (!pool) {
270     AliWarning(Form("No pool found. Centrality %f, ZVertex %f", 
271                     fCentrality, fZVertex));
272     return;
273   }
274
275   if (!pool->IsReady()) {
276     pool->UpdatePool(sTracks);
277     return;
278   }
279
280   if (pool->IsFirstReady()) {
281     // recover events missed before the pool is ready
282     Int_t nEvs = pool->GetCurrentNEvents();
283     if (nEvs>1) {
284       for (Int_t i=0; i<nEvs; ++i) {
285         MiniEvent* evI = pool->GetEvent(i);
286         for (Int_t j=0; j<nEvs; ++j) {
287           MiniEvent* evJ = pool->GetEvent(j);
288           if (i==j) {
289             Correlate(*evI, *evJ, kSameEvt);
290           } else {
291             Correlate(*evI, *evJ, kDiffEvt, 1.0);
292           }
293         }
294       }
295     }
296   } else { /* standard case: same event, then mix*/
297     Correlate(*sTracks, *sTracks, kSameEvt);  
298     Int_t nMix = pool->GetCurrentNEvents();
299     for (Int_t jMix=0; jMix<nMix; ++jMix) {
300       MiniEvent* bgTracks = pool->GetEvent(jMix);
301       Correlate(*sTracks, *bgTracks, kDiffEvt, 1.0);
302     }
303   }
304
305   pool->UpdatePool(sTracks);
306   PostData(1, fOutputList);
307   return;
308 }
309
310 //________________________________________________________________________
311 void AliDhcTask::GetESDTracks(MiniEvent* miniEvt)
312 {
313   // Loop twice: 1. Count sel. tracks. 2. Fill vector.
314
315   const AliESDVertex *vtxSPD = fESD->GetPrimaryVertexSPD();
316   if (!vtxSPD)
317     return;
318
319   Int_t nTrax = fESD->GetNumberOfTracks();
320   if (fVerbosity > 2)
321     AliInfo(Form("%d tracks in event",nTrax));
322
323   // Loop 1.
324   Int_t nSelTrax = 0;
325   TObjArray arr(nTrax);
326   arr.SetOwner(1);
327
328   for (Int_t i = 0; i < nTrax; ++i) {
329     AliESDtrack* esdtrack = fESD->GetTrack(i);
330     if (!esdtrack) {
331       AliError(Form("Couldn't get ESD track %d\n", i));
332       continue;
333     }
334     Bool_t trkOK = fEsdTrackCutsTPCOnly->AcceptTrack(esdtrack);
335     if (!trkOK)
336       continue;
337     Double_t pt = esdtrack->Pt();
338     Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
339     if (!ptOK)
340       continue;
341     Double_t eta = esdtrack->Eta();
342     if (TMath::Abs(eta) > fEtaMax)
343       continue;
344
345     // create a tpc only track
346     AliESDtrack *newtrack = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
347     if(!newtrack)
348       continue;
349     if (newtrack->Pt()<=0) {
350       delete newtrack;
351       continue;
352     }
353
354     AliExternalTrackParam exParam;
355     Bool_t relate = newtrack->RelateToVertexTPC(vtxSPD,fESD->GetMagneticField(),kVeryBig,&exParam);
356     if (!relate) {
357       delete newtrack;
358       continue;
359     }
360
361     // set the constraint parameters to the track
362     newtrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
363
364     pt = newtrack->Pt();
365     ptOK = pt >= fPtMin && pt < fPtMax;
366     if (!ptOK) {
367       delete newtrack;
368       continue;
369     }
370     eta  = esdtrack->Eta();
371     if (TMath::Abs(eta) > fEtaMax) {
372       delete newtrack;
373       continue;
374     }
375     arr.Add(newtrack);
376     nSelTrax++;
377   }
378
379   if (miniEvt)
380     miniEvt->reserve(nSelTrax);
381   else {
382     AliError("!miniEvt");
383     return;
384   }
385
386   // Loop 2.
387   for (Int_t i = 0; i < nSelTrax; ++i) {
388     AliESDtrack* esdtrack = static_cast<AliESDtrack*>(arr.At(i));
389     if (!esdtrack) {
390       AliError(Form("Couldn't get ESD track %d\n", i));
391       continue;
392     }
393     Double_t pt = esdtrack->Pt();
394     Double_t eta  = esdtrack->Eta();
395     Double_t phi  = esdtrack->Phi();
396     Int_t    sign = esdtrack->Charge() > 0 ? 1 : -1;
397     miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
398   }
399   //  return miniEvt;
400 }
401
402 //________________________________________________________________________
403 void AliDhcTask::GetAODTracks(MiniEvent* miniEvt)
404 {
405   // Loop twice: 1. Count sel. tracks. 2. Fill vector.
406
407   Int_t nTrax = fAOD->GetNumberOfTracks();
408   Int_t nSelTrax = 0;
409
410   if (fVerbosity > 2)
411     AliInfo(Form("%d tracks in event",nTrax));
412
413   // Loop 1.
414   for (Int_t i = 0; i < nTrax; ++i) {
415     AliAODTrack* aodtrack = fAOD->GetTrack(i);
416     if (!aodtrack) {
417       AliError(Form("Couldn't get AOD track %d\n", i));
418       continue;
419     }
420     // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
421     UInt_t tpcOnly = 1 << 7;
422     Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly);
423     if (!trkOK)
424       continue;
425     Double_t pt = aodtrack->Pt();
426     Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
427     if (!ptOK)
428       continue;
429     Double_t eta = aodtrack->Eta();
430     if (TMath::Abs(eta) > fEtaMax)
431       continue;
432     nSelTrax++;
433   }
434
435   if (miniEvt)
436     miniEvt->reserve(nSelTrax);
437   else {
438     AliError("!miniEvt");
439     return;
440   }
441   
442   // Loop 2.  
443   for (Int_t i = 0; i < nTrax; ++i) {
444     AliAODTrack* aodtrack = fAOD->GetTrack(i);
445     if (!aodtrack) {
446       AliError(Form("Couldn't get AOD track %d\n", i));
447       continue;
448     }
449     
450     // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
451     UInt_t tpcOnly = 1 << 7;
452     Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly);
453     if (!trkOK)
454       continue;
455     Double_t pt = aodtrack->Pt();
456     Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
457     if (!ptOK)
458       continue;
459     Double_t eta  = aodtrack->Eta();
460     if (TMath::Abs(eta) > fEtaMax)
461       continue;
462
463     Double_t phi  = aodtrack->Phi();
464     Int_t    sign = aodtrack->Charge() > 0 ? 1 : -1;
465     miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
466   }
467   //  return miniEvt;
468 }
469
470 //________________________________________________________________________
471 Double_t AliDhcTask::DeltaPhi(Double_t phia, Double_t phib,
472                               Double_t rangeMin, Double_t rangeMax) const
473 {
474   Double_t dphi = -999;
475   Double_t pi = TMath::Pi();
476   
477   if (phia < 0)         phia += 2*pi;
478   else if (phia > 2*pi) phia -= 2*pi;
479   if (phib < 0)         phib += 2*pi;
480   else if (phib > 2*pi) phib -= 2*pi;
481   dphi = phib - phia;
482   if (dphi < rangeMin)      dphi += 2*pi;
483   else if (dphi > rangeMax) dphi -= 2*pi;
484   
485   return dphi;
486 }
487
488 //________________________________________________________________________
489 Int_t AliDhcTask::Correlate(const MiniEvent &evt1, const MiniEvent &evt2, 
490                             Int_t pairing, Double_t /*weight*/)
491 {
492   // Triggered angular correlations. If pairing is kSameEvt, particles
493   // within evt1 are correlated. If kDiffEvt, correlate triggers from
494   // evt1 with partners from evt2.
495
496   Int_t cbin = fHCent->FindBin(fCentrality);
497   if (fHCent->IsBinOverflow(cbin) ||
498       fHCent->IsBinUnderflow(cbin))
499     return 0;
500
501   Int_t zbin = fHZvtx->FindBin(fZVertex);
502   if (fHZvtx->IsBinOverflow(zbin) ||
503       fHZvtx->IsBinUnderflow(zbin))
504     return 0;
505
506   Int_t iMax = evt1.size();
507   Int_t jMax = evt2.size();
508
509   TH2  **hist = fHMs;
510   if (pairing == kSameEvt) {
511     hist = fHSs;
512     fHCent->AddBinContent(cbin);
513     fHZvtx->AddBinContent(zbin);
514   }
515
516   Int_t nZvtx = fHZvtx->GetNbinsX();
517   Int_t nPtTrig = fHPtTrg->GetNbinsX();
518   Int_t nPtAssc = fHPtAss->GetNbinsX();
519
520   Int_t globIndex = (cbin-1)*nZvtx*nPtTrig*nPtAssc+(zbin-1)*nPtTrig*nPtAssc;
521
522   for (Int_t i=0; i<iMax; ++i) {
523
524     // Trigger particles
525     const AliMiniTrack &a(evt1.at(i));
526
527     Float_t pta  = a.Pt();
528     Int_t abin = fHPtTrg->FindBin(pta);
529     if (fHPtTrg->IsBinOverflow(abin) ||
530         fHPtTrg->IsBinUnderflow(abin))
531       continue;
532
533     if (pairing == kSameEvt) {
534       fHistPt->Fill(pta);
535       fHTrk->Fill(a.Phi(),a.Eta());
536       fHPtTrg->AddBinContent(abin);
537     }
538
539     for (Int_t j=0; j<jMax; ++j) {
540       // Associated particles
541       if (pairing == kSameEvt && i==j)
542         continue;
543
544       const AliMiniTrack &b(evt2.at(j));
545       
546       Float_t ptb  = b.Pt();
547       if (pta < ptb) 
548         continue;
549
550       Int_t bbin = fHPtTrg->FindBin(ptb);
551       if (fHPtAss->IsBinOverflow(bbin) ||
552           fHPtAss->IsBinUnderflow(bbin))
553         continue;
554
555       Float_t dphi = DeltaPhi(a.Phi(), b.Phi());
556       Float_t deta = a.Eta() - b.Eta();
557
558       Int_t index = globIndex+(abin-1)*nPtAssc+(bbin-1);
559       hist[index]->Fill(dphi,deta);
560
561       if (pairing == kSameEvt) {
562         fHPtAss->AddBinContent(bbin);
563         fMeanPtTrg[cbin-1]->Fill(pta,ptb,pta);
564         fMeanPtAss[cbin-1]->Fill(pta,ptb,ptb);
565         fMean2PtTrg[cbin-1]->Fill(pta,ptb,pta*pta);
566         fMean2PtAss[cbin-1]->Fill(pta,ptb,ptb*ptb);
567       }
568     }
569   }
570
571   return 0;
572 }
573
574 //________________________________________________________________________
575 void AliDhcTask::Terminate(Option_t *) 
576 {
577   // Draw result to the screen
578   // Called once at the end of the query
579
580   delete fPoolMgr;
581
582   fHCent->SetEntries(fHCent->Integral());
583   fHZvtx->SetEntries(fHZvtx->Integral());
584   fHPtTrg->SetEntries(fHPtTrg->Integral());
585   fHPtAss->SetEntries(fHPtAss->Integral());
586
587   if (gROOT->IsBatch())
588     return;
589
590   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
591   if (!fOutputList) {
592     AliError("Output list not available");
593     return;
594   }
595   
596   fHistPt = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistPt"));
597   if (!fHistPt) {
598     AliError("ERROR: fHistPt not available\n");
599     return;
600   }
601   TCanvas *c1 = new TCanvas("AliDhcTask","Pt",10,10,510,510);
602   c1->cd(1)->SetLogy();
603   fHistPt->DrawCopy("E");
604 }
605
606 //________________________________________________________________________
607 Bool_t AliDhcTask::VertexOk(TObject* obj) const
608 {
609   // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
610   
611   Int_t nContributors  = 0;
612   Double_t zVertex     = 999;
613   TString name("");
614   
615   if (obj->InheritsFrom("AliESDEvent")) {
616     AliESDEvent* esdevt = (AliESDEvent*) obj;
617     const AliESDVertex* vtx = esdevt->GetPrimaryVertex();
618     if (!vtx)
619       return 0;
620     nContributors = vtx->GetNContributors();
621     zVertex       = vtx->GetZ();
622     name          = vtx->GetName();
623   }
624   else if (obj->InheritsFrom("AliAODEvent")) { 
625     AliAODEvent* aodevt = (AliAODEvent*) obj;
626     if (aodevt->GetNumberOfVertices() < 1)
627       return 0;
628     const AliAODVertex* vtx = aodevt->GetPrimaryVertex();
629     nContributors = vtx->GetNContributors();
630     zVertex       = vtx->GetZ();
631     name          = vtx->GetName();
632   }
633   
634   // Reject if TPC-only vertex
635   if (name.CompareTo("TPCVertex")==0)
636     return kFALSE;
637   
638   // Check # contributors and range...
639   if( nContributors < 1 || TMath::Abs(zVertex) > fZVtxMax ) {
640     return kFALSE;
641   }
642   
643   return kTRUE;
644 }