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