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