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     AliInfo(Form("Event REJECTED (centrality out of range). fCentrality = %.1f", fCentrality));
 387     return;
 388   }
 389
 390   // Get array of selected tracks
 391   if (dType == kESD) {
 392     GetESDTracks(sTracks);
 393   }
 394   if (dType == kAOD) {
 395     GetAODTracks(sTracks);
 396   }
 397
 398   // Get pool containing tracks from other events like this one
 399   AliEvtPool* pool = fPoolMgr->GetEventPool(fCentrality, fZVertex);
 400   if (!pool) {
 401     AliWarning(Form("No pool found. Centrality %f, ZVertex %f",
 402                     fCentrality, fZVertex));
 403     return;
 404   }
 405
 406   if (!pool->IsReady()) {
 407     pool->UpdatePool(sTracks);
 408     return;
 409   }
 410
 411   if (pool->IsFirstReady()) {
 412     // recover events missed before the pool is ready
 413     Int_t nEvs = pool->GetCurrentNEvents();
 414     if (nEvs>1) {
 415       for (Int_t i=0; i<nEvs; ++i) {
 416         MiniEvent* evI = pool->GetEvent(i);
 417         for (Int_t j=0; j<nEvs; ++j) {
 418           MiniEvent* evJ = pool->GetEvent(j);
 419           if (i==j) {
 420             Correlate(*evI, *evJ, kSameEvt);
 421           } else {
 422             Correlate(*evI, *evJ, kDiffEvt);
 423           }
 424         }
 425       }
 426     }
 427   } else { /* standard case: same event, then mix*/
 428     Correlate(*sTracks, *sTracks, kSameEvt);
 429     Int_t nMix = pool->GetCurrentNEvents();
 430     for (Int_t jMix=0; jMix<nMix; ++jMix) {
 431       MiniEvent* bgTracks = pool->GetEvent(jMix);
 432       Correlate(*sTracks, *bgTracks, kDiffEvt);
 433     }
 434   }
 435
 436   pool->UpdatePool(sTracks);
 437   PostData(1, fOutputList);
 438 }
 439
 440 //________________________________________________________________________
 441 void AliDhcTask::GetESDTracks(MiniEvent* miniEvt)
 442 {
 443   // Loop twice: 1. Count sel. tracks. 2. Fill vector.
 444
 445   if (fTracksName.IsNull()) {
 446     const AliESDVertex *vtxSPD = fESD->GetPrimaryVertexSPD();
 447     if (!vtxSPD)
 448       return;
 449
 450     Int_t nTrax = fESD->GetNumberOfTracks();
 451     if (fVerbosity > 2)
 452       AliInfo(Form("%d tracks in event",nTrax));
 453
 454     // Loop 1.
 455     Int_t nSelTrax = 0;
 456     TObjArray arr(nTrax);
 457     arr.SetOwner(1);
 458
 459     for (Int_t i = 0; i < nTrax; ++i) {
 460       AliESDtrack* esdtrack = fESD->GetTrack(i);
 461       if (!esdtrack) {
 462         AliError(Form("Couldn't get ESD track %d\n", i));
 463         continue;
 464       }
 465       Bool_t trkOK = fEsdTPCOnly->AcceptTrack(esdtrack);
 466       if (!trkOK)
 467         continue;
 468       Double_t pt = esdtrack->Pt();
 469       Bool_t ptOK = pt >= fPtMin && pt < fPtMax;
 470       if (!ptOK)
 471         continue;
 472       Double_t eta = esdtrack->Eta();
 473       if (TMath::Abs(eta) > fEtaMax)
 474         continue;
 475
 476       // create a tpc only track
 477       AliESDtrack *newtrack = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 478       if(!newtrack)
 479         continue;
 480       if (newtrack->Pt()<=0) {
 481         delete newtrack;
 482         continue;
 483       }
 484
 485       AliExternalTrackParam exParam;
 486       Bool_t relate = newtrack->RelateToVertexTPC(vtxSPD,fESD->GetMagneticField(),kVeryBig,&exParam);
 487       if (!relate) {
 488         delete newtrack;
 489         continue;
 490       }
 491
 492       // set the constraint parameters to the track
 493       newtrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 494
 495       pt = newtrack->Pt();
 496       ptOK = pt >= fPtMin && pt < fPtMax;
 497       if (!ptOK) {
 498         delete newtrack;
 499         continue;
 500       }
 501       eta  = esdtrack->Eta();
 502       if (TMath::Abs(eta) > fEtaMax) {
 503         delete newtrack;
 504         continue;
 505       }
 506       arr.Add(newtrack);
 507       nSelTrax++;
 508     }
 509
 510     for(Int_t itrack = 0; itrack < nSelTrax; itrack++) {
 511       AliVTrack *esdtrack = static_cast<AliESDtrack*>(arr.At(itrack));
 512       if(!esdtrack) {
 513         AliError(Form("ERROR: Could not retrieve esdtrack %d",itrack));
 514         continue;
 515       }
 516       Double_t pt = esdtrack->Pt();
 517       Double_t eta  = esdtrack->Eta();
 518       Double_t phi  = esdtrack->Phi();
 519       Int_t    sign = esdtrack->Charge() > 0 ? 1 : -1;
 520       miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
 521     }
 522     return;
 523   }
 524
 525   TList *list = InputEvent()->GetList();
 526   TClonesArray *tcaTracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracksName));
 527
 528   if(!tcaTracks){
 529     AliError("Ptr to tcaTracks zero");
 530     return;
 531   }
 532
 533   const Int_t ntracks = tcaTracks->GetEntries();
 534   if (miniEvt)
 535     miniEvt->reserve(ntracks);
 536   else {
 537     AliError("Ptr to miniEvt zero");
 538     return;
 539   }
 540
 541   for(Int_t itrack = 0; itrack < ntracks; itrack++) {
 542     AliVTrack *esdtrack = static_cast<AliESDtrack*>(tcaTracks->At(itrack));
 543     if(!esdtrack) {
 544       AliError(Form("ERROR: Could not retrieve esdtrack %d",itrack));
 545       continue;
 546     }
 547     Double_t pt = esdtrack->Pt();
 548     Double_t eta  = esdtrack->Eta();
 549     Double_t phi  = esdtrack->Phi();
 550     Int_t    sign = esdtrack->Charge() > 0 ? 1 : -1;
 551     miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign));
 552   }
 553
 554   if (fFillMuons) {
 555     Double_t ptMu, etaMu, phiMu;
 556     Int_t signMu;
 557     // count good muons
 558     Int_t nGoodMuons = 0;
 559     for (Int_t iMu = 0; iMu<fESD->GetNumberOfMuonTracks(); iMu++) {
 560       AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iMu);
 561       if(muonTrack) {
 562           if (IsGoodMUONtrack(*muonTrack)) nGoodMuons++;
 563       }
 564     }
 565     miniEvt->reserve(miniEvt->size()+nGoodMuons);
 566     // fill them into the mini event
 567     for (Int_t iMu = 0; iMu<fESD->GetNumberOfMuonTracks(); iMu++) {
 568       AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iMu);
 569       if(muonTrack) {
 570         if(!IsGoodMUONtrack(*muonTrack)) continue;
 571         ptMu   = muonTrack->Pt();
 572         etaMu  = muonTrack->Eta();
 573         phiMu  = muonTrack->Phi();
 574         signMu = muonTrack->Charge() > 0 ? 1 : -1;
 575         miniEvt->push_back(AliMiniTrack(ptMu, etaMu, phiMu, signMu));
 576       }
 577     }
 578   }
 579 }
 580
 581 //________________________________________________________________________
 582 void AliDhcTask::GetAODTracks(MiniEvent* miniEvt)
 583 {
 584   // Loop twice: 1. Count sel. tracks. 2. Fill vector.
 585
 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 }
 647
 648 //________________________________________________________________________
 649 Double_t AliDhcTask::DeltaPhi(Double_t phia, Double_t phib,
 650                               Double_t rangeMin, Double_t rangeMax) const
 651 {
 652   Double_t dphi = -999;
 653   Double_t pi = TMath::Pi();
 654
 655   if (phia < 0)         phia += 2*pi;
 656   else if (phia > 2*pi) phia -= 2*pi;
 657   if (phib < 0)         phib += 2*pi;
 658   else if (phib > 2*pi) phib -= 2*pi;
 659   dphi = phib - phia;
 660   if (dphi < rangeMin)      dphi += 2*pi;
 661   else if (dphi > rangeMax) dphi -= 2*pi;
 662
 663   return dphi;
 664 }
 665
 666 //________________________________________________________________________
 667 Int_t AliDhcTask::Correlate(const MiniEvent &evt1, const MiniEvent &evt2, Int_t pairing)
 668 {
 669   // Triggered angular correlations. If pairing is kSameEvt, particles
 670   // within evt1 are correlated. If kDiffEvt, correlate triggers from
 671   // evt1 with partners from evt2.
 672
 673   Int_t cbin = fHCent->FindBin(fCentrality);
 674   if (fHCent->IsBinOverflow(cbin) ||
 675       fHCent->IsBinUnderflow(cbin))
 676     return 0;
 677
 678   Int_t zbin = fHZvtx->FindBin(fZVertex);
 679   if (fHZvtx->IsBinOverflow(zbin) ||
 680       fHZvtx->IsBinUnderflow(zbin))
 681     return 0;
 682
 683   Int_t iMax = evt1.size();
 684   Int_t jMax = evt2.size();
 685
 686   TH2  **hist = fHMs;
 687   if (pairing == kSameEvt) {
 688     hist = fHSs;
 689     fHCent->AddBinContent(cbin);
 690     fHZvtx->AddBinContent(zbin);
 691   }
 692
 693   Int_t nZvtx = fHZvtx->GetNbinsX();
 694   Int_t nPtTrig = fHPtTrg->GetNbinsX();
 695   Int_t nPtAssc = fHPtAss->GetNbinsX();
 696
 697   Int_t globIndex = (cbin-1)*nZvtx*nPtTrig*nPtAssc+(zbin-1)*nPtTrig*nPtAssc;
 698
 699   Int_t ptindex = (Int_t)fIndex->Eval(1,1,zbin,cbin);
 700
 701   fHPtTrgEvt->Reset();
 702   for (Int_t i=0; i<iMax; ++i) {
 703     const AliMiniTrack &a(evt1.at(i));
 704     Float_t pta  = a.Pt();
 705     fHPtTrgEvt->Fill(pta);
 706     if (pairing == kSameEvt) {
 707       fHPts[ptindex]->Fill(pta);
 708     }
 709   }
 710
 711   Bool_t bCountTrg; // only count the trigger if an associated particle is found
 712
 713   for (Int_t i=0; i<iMax; ++i) {
 714     // Trigger particles
 715     const AliMiniTrack &a(evt1.at(i));
 716
 717     Float_t pta  = a.Pt();
 718     Float_t etaa = a.Eta();
 719     Float_t phia = a.Phi();
 720     Int_t abin = fHPtTrg->FindBin(pta);
 721     if (fHPtTrg->IsBinOverflow(abin) || fHPtTrg->IsBinUnderflow(abin))
 722       continue;
 723
 724     if (etaa<fEtaTLo || etaa>fEtaTHi)
 725       continue;
 726
 727     // efficiency weighting
 728     Double_t effWtT = 1.0;
 729     if (fHEffT) {
 730       // trigger particle
 731       const Int_t nEffDimT = fHEffT->GetNdimensions();
 732       Int_t effBinT[nEffDimT];
 733       effBinT[0] = fHEffT->GetAxis(0)->FindBin(etaa);
 734       effBinT[1] = fHEffT->GetAxis(1)->FindBin(pta);
 735       effBinT[2] = fHEffT->GetAxis(2)->FindBin(fCentrality);
 736       effBinT[3] = fHEffT->GetAxis(3)->FindBin(fZVertex);
 737       if (nEffDimT>4) {
 738         effBinT[4] = fHEffT->GetAxis(4)->FindBin(phia);
 739       }
 740       effWtT = fHEffT->GetBinContent(effBinT);
 741     }
 742
 743     if (pairing == kSameEvt) {
 744       fHTrk->Fill(phia,etaa);
 745       fHPtTrgNorm1S->Fill(pta,fCentrality,fZVertex,effWtT);
 746     } else {
 747       fHPtTrgNorm1M->Fill(pta,fCentrality,fZVertex,effWtT);
 748     }
 749     bCountTrg = kFALSE;
 750
 751     for (Int_t j=0; j<jMax; ++j) {
 752       // Associated particles
 753       if (pairing == kSameEvt && i==j)
 754         continue;
 755
 756       const AliMiniTrack &b(evt2.at(j));
 757
 758       Float_t ptb  = b.Pt();
 759       Float_t etab = b.Eta();
 760       Float_t phib = b.Phi();
 761       if (pta < ptb)
 762             continue;
 763
 764       Int_t bbin = fHPtAss->FindBin(ptb);
 765       if (fHPtAss->IsBinOverflow(bbin) || fHPtAss->IsBinUnderflow(bbin))
 766         continue;
 767
 768       if (etab<fEtaALo || etab>fEtaAHi)
 769         continue;
 770
 771       Float_t dphi = DeltaPhi(phia, phib);
 772       Float_t deta = etaa - etab;
 773
 774       Int_t index = globIndex+(abin-1)*nPtAssc+(bbin-1);
 775       Double_t weight = 1.0;
 776       // number of trigger weight event-by-event (CMS method)
 777       if (fDoWeights) {
 778         Double_t nTrgs = fHPtTrgEvt->GetBinContent(abin);
 779         if (nTrgs==0.0) {
 780           AliError(Form("Trying to work with number of triggers weight = %f",nTrgs));
 781           return 0;
 782         }
 783         weight *= 1./nTrgs;
 784       }
 785
 786       // efficiency weighting
 787       if (fHEffT) {
 788         // trigger particle
 789         weight *= effWtT;
 790       }
 791       if (fHEffA) {
 792         // associated particle
 793         const Int_t nEffDimA = fHEffA->GetNdimensions();
 794         Int_t effBinA[nEffDimA];
 795         effBinA[0] = fHEffA->GetAxis(0)->FindBin(etab);
 796         effBinA[1] = fHEffA->GetAxis(1)->FindBin(ptb);
 797         effBinA[2] = fHEffA->GetAxis(2)->FindBin(fCentrality);
 798         effBinA[3] = fHEffA->GetAxis(3)->FindBin(fZVertex);
 799         if (nEffDimA>4) {
 800           effBinA[4] = fHEffA->GetAxis(4)->FindBin(phib);
 801         }
 802         weight *= fHEffA->GetBinContent(effBinA);
 803       }
 804       hist[index]->Fill(dphi,deta,weight);
 805       bCountTrg = kTRUE;
 806
 807       if (pairing == kSameEvt) {
 808         fHPtAss->AddBinContent(bbin);
 809       }
 810     }
 811     if (bCountTrg) {
 812       if (pairing == kSameEvt) {
 813         fHPtTrgNorm2S->Fill(pta,fCentrality,fZVertex,effWtT);
 814       } else {
 815         fHPtTrgNorm2M->Fill(pta,fCentrality,fZVertex,effWtT);
 816       }
 817     }
 818   }
 819
 820   return 1;
 821 }
 822
 823 //________________________________________________________________________
 824 void AliDhcTask::Terminate(Option_t *)
 825 {
 826   // Draw result to the screen
 827   // Called once at the end of the query
 828
 829   delete fPoolMgr;
 830   fHCent->SetEntries(fHCent->Integral());
 831   fHZvtx->SetEntries(fHZvtx->Integral());
 832   fHPtTrg->SetEntries(fHPtTrg->Integral());
 833   fHPtAss->SetEntries(fHPtAss->Integral());
 834 }
 835
 836 //________________________________________________________________________
 837 Bool_t AliDhcTask::VertexOk(TObject* obj) const
 838 {
 839   // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
 840
 841   Int_t nContributors  = 0;
 842   Double_t zVertex     = 999;
 843   TString name("");
 844
 845   if (obj->InheritsFrom("AliESDEvent")) {
 846     AliESDEvent* esdevt = (AliESDEvent*) obj;
 847     const AliESDVertex* vtx = esdevt->GetPrimaryVertex();
 848     if (!vtx)
 849       return 0;
 850     nContributors = vtx->GetNContributors();
 851     zVertex       = vtx->GetZ();
 852     name          = vtx->GetName();
 853   }
 854   else if (obj->InheritsFrom("AliAODEvent")) {
 855     AliAODEvent* aodevt = (AliAODEvent*) obj;
 856     if (aodevt->GetNumberOfVertices() < 1)
 857       return 0;
 858     const AliAODVertex* vtx = aodevt->GetPrimaryVertex();
 859     nContributors = vtx->GetNContributors();
 860     zVertex       = vtx->GetZ();
 861     name          = vtx->GetName();
 862   }
 863
 864   // Reject if TPC-only vertex
 865   if (name.CompareTo("TPCVertex")==0)
 866     return kFALSE;
 867
 868   // Check # contributors and range...
 869   if( nContributors < 1 || TMath::Abs(zVertex) > fZVtxMax ) {
 870     return kFALSE;
 871   }
 872
 873   return kTRUE;
 874 }
 875
 876 //________________________________________________________________________
 877 Bool_t AliDhcTask::IsGoodMUONtrack(AliESDMuonTrack &track)
 878 {
 879     // Applying track cuts for MUON tracks
 880     if(!track.ContainTrackerData()) return kFALSE;
 881     if(!track.ContainTriggerData()) return kFALSE;
 882     Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
 883     Double_t eta = track.Eta();
 884     if(!(thetaTrackAbsEnd >= 2. && thetaTrackAbsEnd <= 10. && eta >= -4. && eta <= -2.5)) return kFALSE;
 885     if(track.GetMatchTrigger() <= 0) return kFALSE;
 886     // if(track.Pt() <= 1.0) return kFALSE;
 887     //  if(track.GetNormalizedChi2() >= 4.0) return kFALSE;
 888
 889     return kTRUE;
 890 }