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