PWG4/UserTasks/DiHadronCorrelations/AliDhcTask_opt.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 "THnSparse.h"
  12 #include "TROOT.h"
  13 #include "TTree.h"
  14 #include "AliAODEvent.h"
  15 #include "AliAODInputHandler.h"
  16 #include "AliAnalysisManager.h"
  17 #include "AliAnalysisTask.h"
  18 #include "AliCentrality.h"
  19 #include "AliDhcTask_opt.h"
  20 #include "AliESDEvent.h"
  21 #include "AliESDInputHandler.h"
  22 #include "AliESDtrackCuts.h"
  23 #include "KiddiePoolClasses.h"
  24 #include "AliVParticle.h"
  25
  26 ClassImp(AliDhcTask)
  27
  28 //________________________________________________________________________
  29 AliDhcTask::AliDhcTask(const char *name)
  30 : AliAnalysisTaskSE(name), fVerbosity(0), fEtaMax(1), fZVtxMax(10), fPtMin(0.25), fPtMax(15),
  31   fESD(0), fAOD(0), fOutputList(0), fHistPt(0), fHEvt(0), fHTrk(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 KiddiePoolManager();
 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 = 0; // Vector of selected MiniTracks.
 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     sTracks = GetESDTrax();
 262   }
 263   if (dType == kAOD) {
 264     sTracks = GetAODTrax();
 265   }
 266
 267   // Get pool containing tracks from other events like this one
 268   KiddiePool* pool = fPoolMgr->GetEventPool(fCentrality, fZVertex);
 269   if (!pool) {
 270     AliWarning(Form("No pool found. Centrality %f, ZVertex %f", fCentrality, fZVertex));
 271     return;
 272   }
 273
 274   if (!pool->IsReady()) {
 275     pool->UpdatePool(sTracks);
 276     return;
 277   }
 278
 279   if (pool->IsFirstReady()) {
 280     // recover events missed before the pool is ready
 281     Int_t nEvs = pool->GetCurrentNEvents();
 282     if (nEvs>1) {
 283       for (Int_t i=0; i<nEvs; ++i) {
 284         MiniEvent* evI = pool->GetEvent(i);
 285         for (Int_t j=0; j<nEvs; ++j) {
 286           MiniEvent* evJ = pool->GetEvent(j);
 287           if (i==j) {
 288             Correlate(*evI, *evJ, kSameEvt);
 289           } else {
 290             Correlate(*evI, *evJ, kDiffEvt, 1.0);
 291           }
 292         }
 293       }
 294     }
 295   } else { /* standard case: same event, then mix*/
 296     Correlate(*sTracks, *sTracks, kSameEvt);
 297     Int_t nMix = pool->GetCurrentNEvents();
 298     for (Int_t jMix=0; jMix<nMix; ++jMix) {
 299       MiniEvent* bgTracks = pool->GetEvent(jMix);
 300       Correlate(*sTracks, *bgTracks, kDiffEvt, 1.0);
 301     }
 302   }
 303
 304   pool->UpdatePool(sTracks);
 305   PostData(1, fOutputList);
 306   return;
 307 }
 308
 309 //________________________________________________________________________
 310 MiniEvent* AliDhcTask::GetESDTrax() const
 311 {
 312   // Loop twice: 1. Count sel. tracks. 2. Fill vector.
 313
 314   const AliESDVertex *vtxSPD = fESD->GetPrimaryVertexSPD();
 315   if (!vtxSPD)
 316     return 0;
 317
 318   Int_t nTrax = fESD->GetNumberOfTracks();
 319   if (fVerbosity > 2)
 320     AliInfo(Form("%d tracks in event",nTrax));
 321
 322   // Loop 1.
 323   Int_t nSelTrax = 0;
 324   TObjArray arr(nTrax);
 325   arr.SetOwner(1);
 326
 327   for (Int_t i = 0; i < nTrax; ++i) {
 328     AliESDtrack* esdtrack = fESD->GetTrack(i);
 329     if (!esdtrack) {
 330       AliError(Form("Couldn't get ESD track %d\n", i));
 331       continue;
 332     }
 333     Bool_t trkOK = fEsdTrackCutsTPCOnly->AcceptTrack(esdtrack);
 334     if (!trkOK)
 335       continue;
 336     Double_t pt = esdtrack->Pt();
 337     Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
 338     if (!ptOK)
 339       continue;
 340     Double_t eta = esdtrack->Eta();
 341     if (TMath::Abs(eta) > fEtaMax)
 342       continue;
 343
 344     // create a tpc only track
 345     AliESDtrack *newtrack = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 346     if(!newtrack)
 347       continue;
 348     if (newtrack->Pt()<=0) {
 349       delete newtrack;
 350       continue;
 351     }
 352
 353     AliExternalTrackParam exParam;
 354     Bool_t relate = newtrack->RelateToVertexTPC(vtxSPD,fESD->GetMagneticField(),kVeryBig,&exParam);
 355     if (!relate) {
 356       delete newtrack;
 357       continue;
 358     }
 359
 360     // set the constraint parameters to the track
 361     newtrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 362
 363     pt = newtrack->Pt();
 364     ptOK = pt >= fPtMin && pt < fPtMax;
 365     if (!ptOK) {
 366       delete newtrack;
 367       continue;
 368     }
 369     eta  = esdtrack->Eta();
 370     if (TMath::Abs(eta) > fEtaMax) {
 371       delete newtrack;
 372       continue;
 373     }
 374     arr.Add(newtrack);
 375     nSelTrax++;
 376   }
 377
 378   MiniEvent* miniEvt = new MiniEvent(0);
 379   miniEvt->reserve(nSelTrax);
 380
 381   // Loop 2.
 382   for (Int_t i = 0; i < nSelTrax; ++i) {
 383     AliESDtrack* esdtrack = static_cast<AliESDtrack*>(arr.At(i));
 384     if (!esdtrack) {
 385       AliError(Form("Couldn't get ESD track %d\n", i));
 386       continue;
 387     }
 388     Double_t pt = esdtrack->Pt();
 389     Double_t eta  = esdtrack->Eta();
 390     Double_t phi  = esdtrack->Phi();
 391     Int_t    sign = esdtrack->Charge() > 0 ? 1 : -1;
 392     miniEvt->push_back(MiniTrack(pt, eta, phi, sign));
 393   }
 394   return miniEvt;
 395 }
 396
 397 //________________________________________________________________________
 398 MiniEvent* AliDhcTask::GetAODTrax() const
 399 {
 400   // Loop twice: 1. Count sel. tracks. 2. Fill vector.
 401
 402   Int_t nTrax = fAOD->GetNumberOfTracks();
 403   Int_t nSelTrax = 0;
 404
 405   if (fVerbosity > 2)
 406     AliInfo(Form("%d tracks in event",nTrax));
 407
 408   // Loop 1.
 409   for (Int_t i = 0; i < nTrax; ++i) {
 410     AliAODTrack* aodtrack = fAOD->GetTrack(i);
 411     if (!aodtrack) {
 412       AliError(Form("Couldn't get AOD track %d\n", i));
 413       continue;
 414     }
 415     // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
 416     UInt_t tpcOnly = 1 << 7;
 417     Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly);
 418     if (!trkOK)
 419       continue;
 420     Double_t pt = aodtrack->Pt();
 421     Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
 422     if (!ptOK)
 423       continue;
 424     Double_t eta = aodtrack->Eta();
 425     if (TMath::Abs(eta) > fEtaMax)
 426       continue;
 427     nSelTrax++;
 428   }
 429
 430   MiniEvent* miniEvt = new MiniEvent(0);
 431   miniEvt->reserve(nSelTrax);
 432
 433   // Loop 2.
 434   for (Int_t i = 0; i < nTrax; ++i) {
 435     AliAODTrack* aodtrack = fAOD->GetTrack(i);
 436     if (!aodtrack) {
 437       AliError(Form("Couldn't get AOD track %d\n", i));
 438       continue;
 439     }
 440
 441     // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C
 442     UInt_t tpcOnly = 1 << 7;
 443     Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly);
 444     if (!trkOK)
 445       continue;
 446     Double_t pt = aodtrack->Pt();
 447     Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
 448     if (!ptOK)
 449       continue;
 450     Double_t eta  = aodtrack->Eta();
 451     if (TMath::Abs(eta) > fEtaMax)
 452       continue;
 453
 454     Double_t phi  = aodtrack->Phi();
 455     Int_t    sign = aodtrack->Charge() > 0 ? 1 : -1;
 456     miniEvt->push_back(MiniTrack(pt, eta, phi, sign));
 457   }
 458   return miniEvt;
 459 }
 460
 461 //________________________________________________________________________
 462 Double_t AliDhcTask::DeltaPhi(Double_t phia, Double_t phib,
 463                               Double_t rangeMin, Double_t rangeMax) const
 464 {
 465   Double_t dphi = -999;
 466   Double_t pi = TMath::Pi();
 467
 468   if (phia < 0)         phia += 2*pi;
 469   else if (phia > 2*pi) phia -= 2*pi;
 470   if (phib < 0)         phib += 2*pi;
 471   else if (phib > 2*pi) phib -= 2*pi;
 472   dphi = phib - phia;
 473   if (dphi < rangeMin)      dphi += 2*pi;
 474   else if (dphi > rangeMax) dphi -= 2*pi;
 475
 476   return dphi;
 477 }
 478
 479 //________________________________________________________________________
 480 Int_t AliDhcTask::Correlate(const MiniEvent &evt1, const MiniEvent &evt2,
 481                             Int_t pairing, Double_t /*weight*/)
 482 {
 483   // Triggered angular correlations. If pairing is kSameEvt, particles
 484   // within evt1 are correlated. If kDiffEvt, correlate triggers from
 485   // evt1 with partners from evt2.
 486
 487   Int_t cbin = fHCent->FindBin(fCentrality);
 488   if (fHCent->IsBinOverflow(cbin) ||
 489       fHCent->IsBinUnderflow(cbin))
 490     return 0;
 491
 492   Int_t zbin = fHZvtx->FindBin(fZVertex);
 493   if (fHZvtx->IsBinOverflow(zbin) ||
 494       fHZvtx->IsBinUnderflow(zbin))
 495     return 0;
 496
 497   Int_t iMax = evt1.size();
 498   Int_t jMax = evt2.size();
 499
 500   TH2  **hist = fHMs;
 501   if (pairing == kSameEvt) {
 502     hist = fHSs;
 503     fHCent->AddBinContent(cbin);
 504     fHZvtx->AddBinContent(zbin);
 505   }
 506
 507   Int_t nZvtx = fHZvtx->GetNbinsX();
 508   Int_t nPtTrig = fHPtTrg->GetNbinsX();
 509   Int_t nPtAssc = fHPtAss->GetNbinsX();
 510
 511   Int_t globIndex = (cbin-1)*nZvtx*nPtTrig*nPtAssc+(zbin-1)*nPtTrig*nPtAssc;
 512
 513   for (Int_t i=0; i<iMax; ++i) {
 514
 515     // Trigger particles
 516     const MiniTrack &a(evt1.at(i));
 517
 518     Float_t pta  = a.Pt();
 519     Int_t abin = fHPtTrg->FindBin(pta);
 520     if (fHPtTrg->IsBinOverflow(abin) ||
 521         fHPtTrg->IsBinUnderflow(abin))
 522       continue;
 523
 524     if (pairing == kSameEvt) {
 525       fHistPt->Fill(pta);
 526       fHTrk->Fill(a.Phi(),a.Eta());
 527       fHPtTrg->AddBinContent(abin);
 528     }
 529
 530     for (Int_t j=0; j<jMax; ++j) {
 531       // Associated particles
 532       if (pairing == kSameEvt && i==j)
 533         continue;
 534
 535       const MiniTrack &b(evt2.at(j));
 536
 537       Float_t ptb  = b.Pt();
 538       if (pta < ptb)
 539         continue;
 540
 541       Int_t bbin = fHPtTrg->FindBin(ptb);
 542       if (fHPtAss->IsBinOverflow(bbin) ||
 543           fHPtAss->IsBinUnderflow(bbin))
 544         continue;
 545
 546       Float_t dphi = DeltaPhi(a.Phi(), b.Phi());
 547       Float_t deta = a.Eta() - b.Eta();
 548
 549       Int_t index = globIndex+(abin-1)*nPtAssc+(bbin-1);
 550       hist[index]->Fill(dphi,deta);
 551
 552       if (pairing == kSameEvt) {
 553         fHPtAss->AddBinContent(bbin);
 554         fMeanPtTrg[cbin-1]->Fill(pta,ptb,pta);
 555         fMeanPtAss[cbin-1]->Fill(pta,ptb,ptb);
 556         fMean2PtTrg[cbin-1]->Fill(pta,ptb,pta*pta);
 557         fMean2PtAss[cbin-1]->Fill(pta,ptb,ptb*ptb);
 558       }
 559     }
 560   }
 561
 562   return 0;
 563 }
 564
 565 //________________________________________________________________________
 566 void AliDhcTask::Terminate(Option_t *)
 567 {
 568   // Draw result to the screen
 569   // Called once at the end of the query
 570
 571   delete fPoolMgr;
 572
 573   fHCent->SetEntries(fHCent->Integral());
 574   fHZvtx->SetEntries(fHZvtx->Integral());
 575   fHPtTrg->SetEntries(fHPtTrg->Integral());
 576   fHPtAss->SetEntries(fHPtAss->Integral());
 577
 578   if (gROOT->IsBatch())
 579     return;
 580
 581   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 582   if (!fOutputList) {
 583     AliError("Output list not available");
 584     return;
 585   }
 586
 587   fHistPt = dynamic_cast<TH1F*> (fOutputList->At(0));
 588   if (!fHistPt) {
 589     AliError("ERROR: fHistPt not available\n");
 590     return;
 591   }
 592
 593   TCanvas *c1 = new TCanvas("AliDhcTask","Pt",10,10,510,510);
 594   c1->cd(1)->SetLogy();
 595   fHistPt->DrawCopy("E");
 596 }
 597
 598 //________________________________________________________________________
 599 Bool_t AliDhcTask::VertexOk(TObject* obj) const
 600 {
 601   // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
 602
 603   Int_t nContributors  = 0;
 604   Double_t zVertex     = 999;
 605   TString name("");
 606
 607   if (obj->InheritsFrom("AliESDEvent")) {
 608     AliESDEvent* esdevt = (AliESDEvent*) obj;
 609     const AliESDVertex* vtx = esdevt->GetPrimaryVertex();
 610     if (!vtx)
 611       return 0;
 612     nContributors = vtx->GetNContributors();
 613     zVertex       = vtx->GetZ();
 614     name          = vtx->GetName();
 615   }
 616   else if (obj->InheritsFrom("AliAODEvent")) {
 617     AliAODEvent* aodevt = (AliAODEvent*) obj;
 618     if (aodevt->GetNumberOfVertices() < 1)
 619       return 0;
 620     const AliAODVertex* vtx = aodevt->GetPrimaryVertex();
 621     nContributors = vtx->GetNContributors();
 622     zVertex       = vtx->GetZ();
 623     name          = vtx->GetName();
 624   }
 625
 626   // Reject if TPC-only vertex
 627   if (name.CompareTo("TPCVertex")==0)
 628     return kFALSE;
 629
 630   // Check # contributors and range...
 631   if( nContributors < 1 || TMath::Abs(zVertex) > fZVtxMax ) {
 632     return kFALSE;
 633   }
 634
 635   return kTRUE;
 636 }