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