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