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