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