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