]>
Commit | Line | Data |
---|---|---|
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 | ||
5 | #include "TCanvas.h" | |
6 | #include "TChain.h" | |
8433ff41 | 7 | #include "TFormula.h" |
48a61f36 | 8 | #include "TH1.h" |
9 | #include "TH2.h" | |
d688e049 | 10 | #include "TAxis.h" |
8433ff41 | 11 | #include "TProfile2D.h" |
48a61f36 | 12 | #include "TROOT.h" |
13 | #include "TTree.h" | |
48a61f36 | 14 | #include "AliAODEvent.h" |
15 | #include "AliAODInputHandler.h" | |
16 | #include "AliAnalysisManager.h" | |
17 | #include "AliAnalysisTask.h" | |
18 | #include "AliCentrality.h" | |
19 | #include "AliDhcTask.h" | |
20 | #include "AliESDEvent.h" | |
21 | #include "AliESDInputHandler.h" | |
22 | #include "AliESDtrackCuts.h" | |
8433ff41 | 23 | #include "AliPool.h" |
48a61f36 | 24 | #include "AliVParticle.h" |
25 | ||
c64cb1f6 | 26 | using std::cout; |
27 | using std::endl; | |
28 | ||
48a61f36 | 29 | ClassImp(AliDhcTask) |
30 | ||
31 | //________________________________________________________________________ | |
32 | AliDhcTask::AliDhcTask(const char *name) | |
42036329 | 33 | : AliAnalysisTaskSE(name), fVerbosity(0), fEtaMax(1), fZVtxMax(10), fPtMin(0.25), fPtMax(15), |
d688e049 | 34 | fTrackDepth(1000), fPoolSize(200), fTracksName(), fDoWeights(0), |
42036329 | 35 | fESD(0), fAOD(0), fOutputList(0), fHistPt(0), fHEvt(0), fHTrk(0), fHPtAss(0), fHPtTrg(0), |
36 | fHCent(0), fHZvtx(0), fNbins(0), fHSs(0), fHMs(0), | |
37 | fIndex(0), fMeanPtTrg(0), fMeanPtAss(0), fMean2PtTrg(0), fMean2PtAss(0), | |
d688e049 | 38 | fCentrality(99), fZVertex(99), fEsdTrackCutsTPCOnly(0), fPoolMgr(0), |
39 | fCentMethod("V0M"), fNBdeta(20), fNBdphi(36), | |
40 | fBPtT(0x0), fBPtA(0x0), fBCent(0x0), fBZvtx(0x0), | |
41 | fMixBCent(0x0), fMixBZvtx(0x0) | |
48a61f36 | 42 | { |
43 | // Constructor | |
44 | ||
45 | // Define input and output slots here | |
46 | // Input slot #0 works with a TChain | |
47 | DefineInput(0, TChain::Class()); | |
48 | // Output slot #0 id reserved by the base class for AOD | |
49 | // Output slot #1 writes into a TH1 container | |
50 | DefineOutput(1, TList::Class()); | |
51 | ||
48a61f36 | 52 | fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,Tracks " |
53 | "AOD:header,tracks,vertices,"; | |
d688e049 | 54 | |
55 | Double_t ptt[4] = {0.25, 1.0, 2.0, 15.0}; | |
56 | fBPtT = new TAxis(3,ptt); | |
57 | Double_t pta[4] = {0.25, 1.0, 2.0, 15.0}; | |
58 | fBPtA = new TAxis(3,pta); | |
59 | Double_t cent[2] = {-100.0, 100.0}; | |
60 | fBCent = new TAxis(1,cent); | |
61 | Double_t zvtx[2] = {-10, 10}; | |
62 | fBZvtx = new TAxis(1,zvtx); | |
63 | Double_t centmix[2] = {-100.0, 100.0}; | |
64 | fMixBCent = new TAxis(1,centmix); | |
65 | Double_t zvtxmix[9] = {-10,-6,-4,-2,0,2,4,6,10}; | |
66 | fMixBZvtx = new TAxis(8,zvtxmix); | |
48a61f36 | 67 | } |
68 | ||
69 | //________________________________________________________________________ | |
70 | void AliDhcTask::UserCreateOutputObjects() | |
71 | { | |
72 | // Create histograms | |
73 | // Called once (per slave on PROOF!) | |
74 | ||
48a61f36 | 75 | fOutputList = new TList(); |
76 | fOutputList->SetOwner(1); | |
77 | ||
8433ff41 | 78 | fEsdTrackCutsTPCOnly = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts(); |
79 | //fEsdTrackCutsTPCOnly->SetMinNClustersTPC(70); | |
80 | fEsdTrackCutsTPCOnly->SetMinNCrossedRowsTPC(70); | |
81 | fEsdTrackCutsTPCOnly->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8); | |
48a61f36 | 82 | |
8433ff41 | 83 | BookHistos(); |
84 | InitEventMixer(); | |
48a61f36 | 85 | PostData(1, fOutputList); |
86 | } | |
87 | ||
88 | //________________________________________________________________________ | |
89 | void AliDhcTask::BookHistos() | |
90 | { | |
8433ff41 | 91 | // Book histograms. |
92 | ||
d688e049 | 93 | if (fVerbosity > 1) { |
94 | AliInfo(Form("Number of pt(a) bins: %d", fBPtA->GetNbins())); | |
95 | for (Int_t i=1; i<=fBPtA->GetNbins(); i++) { | |
96 | AliInfo(Form("pt(a) bin %d, %f to %f", i, fBPtA->GetBinLowEdge(i), fBPtA->GetBinUpEdge(i))); | |
97 | } | |
98 | } | |
48a61f36 | 99 | |
d688e049 | 100 | Int_t nPtAssc=fBPtA->GetNbins(); |
101 | Int_t nPtTrig=fBPtT->GetNbins(); | |
102 | Int_t nCent=fBCent->GetNbins(); | |
103 | Int_t nZvtx=fBZvtx->GetNbins(); | |
104 | Double_t ptt[nPtTrig+1]; | |
105 | ptt[0] = fBPtT->GetBinLowEdge(1); | |
106 | for (Int_t i=1; i<=nPtTrig; i++) { | |
107 | ptt[i] = fBPtT->GetBinUpEdge(i); | |
108 | } | |
109 | Double_t pta[nPtAssc+1]; | |
110 | pta[0] = fBPtA->GetBinLowEdge(1); | |
111 | for (Int_t i=1; i<=nPtAssc; i++) { | |
112 | pta[i] = fBPtA->GetBinUpEdge(i); | |
113 | } | |
114 | Double_t cent[nCent+1]; | |
115 | cent[0] = fBCent->GetBinLowEdge(1); | |
116 | for (Int_t i=1; i<=nCent; i++) { | |
117 | cent[i] = fBCent->GetBinUpEdge(i); | |
118 | } | |
119 | Double_t zvtx[nZvtx+1]; | |
120 | zvtx[0] = fBZvtx->GetBinLowEdge(1); | |
121 | for (Int_t i=1; i<=nZvtx; i++) { | |
122 | zvtx[i] = fBZvtx->GetBinUpEdge(i); | |
123 | } | |
48a61f36 | 124 | |
125 | // Event histo | |
8433ff41 | 126 | fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", 30, -15, 15, 101, 0, 101); |
127 | fOutputList->Add(fHEvt); | |
48a61f36 | 128 | // Track histo |
129 | fHTrk = new TH2F("fHTrk", "Track-level variables", | |
8433ff41 | 130 | 100, 0, TMath::TwoPi(), 100, -fEtaMax, +fEtaMax); |
131 | fOutputList->Add(fHTrk); | |
132 | ||
48a61f36 | 133 | // Left over from the tutorial :) |
8433ff41 | 134 | fHistPt = new TH1F("fHistPt", "P_{T} distribution", 200, 0., fPtMax); |
48a61f36 | 135 | fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
136 | fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); | |
137 | fHistPt->SetMarkerStyle(kFullCircle); | |
8433ff41 | 138 | fOutputList->Add(fHistPt); |
139 | ||
140 | fHPtAss = new TH1F("fHPtAss","PtAssoc;P_{T} (GeV/c) [GeV/c]",nPtAssc,pta); | |
141 | fOutputList->Add(fHPtAss); | |
142 | fHPtTrg = new TH1F("fHPtTrg","PtTrig;P_{T} (GeV/c) [GeV/c]",nPtTrig,ptt); | |
143 | fOutputList->Add(fHPtTrg); | |
144 | fHCent = new TH1F("fHCent","Cent;bins",nCent,cent); | |
145 | fOutputList->Add(fHCent); | |
146 | fHZvtx = new TH1F("fHZvtx","Zvertex;bins",nZvtx,zvtx); | |
147 | fOutputList->Add(fHZvtx); | |
148 | ||
149 | fNbins = nPtTrig*nPtAssc*nCent*nZvtx; | |
150 | fHSs = new TH2*[fNbins]; | |
151 | fHMs = new TH2*[fNbins]; | |
152 | ||
153 | fIndex = new TFormula("GlobIndex","(t-1)*[0]*[1]*[2]+(z-1)*[0]*[1]+(x-1)*[0]+(y-1)+0*[4]"); | |
154 | fIndex->SetParameters(nPtTrig,nPtAssc,nZvtx,nCent); | |
155 | fIndex->SetParNames("NTrigBins","NAssocBins", "NZvertexBins", "NCentBins"); | |
156 | //fOutputList->Add(fIndex); | |
157 | ||
158 | fMeanPtTrg = new TProfile2D*[nPtTrig*nPtAssc]; | |
159 | fMeanPtAss = new TProfile2D*[nPtTrig*nPtAssc]; | |
160 | fMean2PtTrg = new TProfile2D*[nPtTrig*nPtAssc]; | |
161 | fMean2PtAss = new TProfile2D*[nPtTrig*nPtAssc]; | |
162 | for (Int_t c=1; c<=nCent; ++c) { | |
163 | TString title(Form("cen=%d (%.1f)",c,fHCent->GetBinCenter(c))); | |
164 | fMeanPtTrg[c-1] = new TProfile2D(Form("hMeanPtTrgCen%d",c),title,nPtTrig,ptt,nPtAssc,pta); | |
165 | fMeanPtAss[c-1] = new TProfile2D(Form("hMeanPtAssCen%d",c),title,nPtTrig,ptt,nPtAssc,pta); | |
166 | fMean2PtTrg[c-1] = new TProfile2D(Form("hMean2PtTrgCen%d",c),title,nPtTrig,ptt,nPtAssc,pta); | |
167 | fMean2PtAss[c-1] = new TProfile2D(Form("hMean2PtAssCen%d",c),title,nPtTrig,ptt,nPtAssc,pta); | |
168 | fOutputList->Add(fMeanPtTrg[c-1]); | |
169 | fOutputList->Add(fMeanPtAss[c-1]); | |
170 | fOutputList->Add(fMean2PtTrg[c-1]); | |
171 | fOutputList->Add(fMean2PtAss[c-1]); | |
172 | } | |
173 | ||
174 | Int_t count = 0; | |
175 | for (Int_t c=1; c<=nCent; ++c) { | |
176 | for (Int_t z=1; z<=nZvtx; ++z) { | |
177 | for (Int_t t=1; t<=nPtTrig; ++t) { | |
178 | for (Int_t a=1; a<=nPtAssc; ++a) { | |
179 | fHSs[count] = 0; | |
180 | fHMs[count] = 0; | |
181 | if (a>t) { | |
182 | ++count; | |
183 | continue; | |
184 | } | |
185 | TString title(Form("cen=%d (%.1f), zVtx=%d (%.1f), trig=%d (%.1f), assc=%d (%.1f)", | |
186 | c,fHCent->GetBinCenter(c), z,fHZvtx->GetBinCenter(z), | |
187 | t,fHPtTrg->GetBinCenter(t),a, fHPtAss->GetBinCenter(a))); | |
188 | fHSs[count] = new TH2F(Form("hS%d",count), Form("Signal %s",title.Data()), | |
d688e049 | 189 | fNBdphi,-0.5*TMath::Pi(),1.5*TMath::Pi(),fNBdeta,-2*fEtaMax,2*fEtaMax); |
8433ff41 | 190 | fHMs[count] = new TH2F(Form("hM%d",count), Form("Signal %s",title.Data()), |
d688e049 | 191 | fNBdphi,-0.5*TMath::Pi(),1.5*TMath::Pi(),fNBdeta,-2*fEtaMax,2*fEtaMax); |
8433ff41 | 192 | fOutputList->Add(fHSs[count]); |
193 | fOutputList->Add(fHMs[count]); | |
194 | if (fVerbosity>5) | |
195 | cout << count << " " << fIndex->Eval(t,a,z,c) << ": " << title << endl; | |
196 | ++count; | |
197 | } | |
198 | } | |
199 | } | |
200 | } | |
48a61f36 | 201 | } |
202 | ||
203 | //________________________________________________________________________ | |
204 | void AliDhcTask::InitEventMixer() | |
205 | { | |
206 | // The effective pool size in events is set by trackDepth, so more | |
207 | // low-mult events are required to maintain the threshold than | |
208 | // high-mult events. Centrality pools are indep. of data histogram | |
209 | // binning, no need to match. | |
210 | ||
48a61f36 | 211 | // Centrality pools |
d688e049 | 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 | } | |
48a61f36 | 218 | |
219 | // Z-vertex pools | |
d688e049 | 220 | Int_t nZvtxBins=fMixBZvtx->GetNbins(); |
221 | Double_t zvtxbin[nZvtxBins+1]; | |
222 | zvtxbin[0] = fMixBZvtx->GetBinLowEdge(1); | |
223 | for (Int_t i=1; i<=nZvtxBins; i++) { | |
224 | zvtxbin[i] = fMixBZvtx->GetBinUpEdge(i); | |
225 | } | |
8433ff41 | 226 | |
227 | fPoolMgr = new AliEvtPoolManager(); | |
d688e049 | 228 | fPoolMgr->SetTargetTrackDepth(fTrackDepth); |
8433ff41 | 229 | if (fVerbosity>4) |
230 | fPoolMgr->SetDebug(1); | |
d688e049 | 231 | fPoolMgr->InitEventPools(fPoolSize, nCentBins, centBins, nZvtxBins, zvtxbin); |
48a61f36 | 232 | } |
233 | ||
234 | //________________________________________________________________________ | |
235 | void AliDhcTask::UserExec(Option_t *) | |
236 | { | |
237 | // Main loop, called for each event. | |
238 | ||
239 | static int iEvent = -1; ++iEvent; | |
240 | ||
241 | if (fVerbosity>2) { | |
242 | if (iEvent % 10 == 0) | |
243 | cout << iEvent << endl; | |
244 | } | |
245 | ||
246 | Int_t dType = -1; // Will be set to kESD or kAOD. | |
f6701c6e | 247 | MiniEvent* sTracks = new MiniEvent(0); // deleted by pool mgr. |
48a61f36 | 248 | |
249 | LoadBranches(); | |
250 | ||
251 | // Get event pointers, check for signs of life | |
252 | fESD = dynamic_cast<AliESDEvent*>(InputEvent()); | |
253 | fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); | |
254 | if (fESD) | |
255 | dType = kESD; | |
256 | else if (fAOD) | |
257 | dType = kAOD; | |
258 | else { | |
259 | AliError("Neither fESD nor fAOD available"); | |
260 | return; | |
261 | } | |
262 | ||
d688e049 | 263 | Bool_t mcgen = 0; |
264 | if (fTracksName.Contains("Gen")) | |
265 | mcgen = 1; | |
266 | ||
48a61f36 | 267 | // Centrality, vertex, other event variables... |
d688e049 | 268 | if (mcgen) { |
269 | fZVertex = 0; | |
270 | TList *list = InputEvent()->GetList(); | |
271 | TClonesArray *tcaTracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracksName)); | |
272 | if (tcaTracks) | |
273 | fCentrality = tcaTracks->GetEntries(); | |
274 | } else { | |
275 | if (dType == kESD) { | |
276 | if (!VertexOk(fESD)) { | |
277 | if (fVerbosity > 1) | |
278 | AliInfo(Form("Event REJECTED (ESD vertex not OK). z = %.1f", fZVertex)); | |
279 | return; | |
280 | } | |
281 | const AliESDVertex* vertex = fESD->GetPrimaryVertex(); | |
282 | fZVertex = vertex->GetZ(); | |
283 | if(fESD->GetCentrality()) { | |
284 | fCentrality = | |
285 | fESD->GetCentrality()->GetCentralityPercentile(fCentMethod); | |
286 | } | |
48a61f36 | 287 | } |
d688e049 | 288 | if (dType == kAOD) { |
289 | const AliAODVertex* vertex = fAOD->GetPrimaryVertex(); | |
290 | fZVertex = vertex->GetZ(); | |
291 | if (!VertexOk(fAOD)) { | |
292 | if (fVerbosity > 1) | |
293 | Info("Exec()", "Event REJECTED (AOD vertex not OK). z = %.1f", fZVertex); | |
294 | return; | |
295 | } | |
296 | const AliCentrality *aodCent = fAOD->GetHeader()->GetCentralityP(); | |
297 | if (aodCent) { | |
298 | fCentrality = aodCent->GetCentralityPercentile(fCentMethod); | |
299 | } | |
48a61f36 | 300 | } |
301 | } | |
d688e049 | 302 | |
48a61f36 | 303 | // Fill Event histogram |
304 | fHEvt->Fill(fZVertex, fCentrality); | |
305 | ||
d688e049 | 306 | if (fCentrality > fBCent->GetXmax() || fCentrality < fBCent->GetXmin()) { |
307 | AliInfo(Form("Event REJECTED (centrality out of range). fCentrality = %.1f", fCentrality)); | |
48a61f36 | 308 | return; |
309 | } | |
310 | ||
311 | // Get array of selected tracks | |
312 | if (dType == kESD) { | |
f6701c6e | 313 | GetESDTracks(sTracks); |
48a61f36 | 314 | } |
315 | if (dType == kAOD) { | |
f6701c6e | 316 | GetAODTracks(sTracks); |
48a61f36 | 317 | } |
318 | ||
319 | // Get pool containing tracks from other events like this one | |
8433ff41 | 320 | AliEvtPool* pool = fPoolMgr->GetEventPool(fCentrality, fZVertex); |
48a61f36 | 321 | if (!pool) { |
c71ff6f4 | 322 | AliWarning(Form("No pool found. Centrality %f, ZVertex %f", |
323 | fCentrality, fZVertex)); | |
48a61f36 | 324 | return; |
325 | } | |
326 | ||
327 | if (!pool->IsReady()) { | |
328 | pool->UpdatePool(sTracks); | |
329 | return; | |
330 | } | |
331 | ||
332 | if (pool->IsFirstReady()) { | |
333 | // recover events missed before the pool is ready | |
334 | Int_t nEvs = pool->GetCurrentNEvents(); | |
335 | if (nEvs>1) { | |
336 | for (Int_t i=0; i<nEvs; ++i) { | |
337 | MiniEvent* evI = pool->GetEvent(i); | |
338 | for (Int_t j=0; j<nEvs; ++j) { | |
339 | MiniEvent* evJ = pool->GetEvent(j); | |
340 | if (i==j) { | |
341 | Correlate(*evI, *evJ, kSameEvt); | |
342 | } else { | |
d688e049 | 343 | Correlate(*evI, *evJ, kDiffEvt); |
48a61f36 | 344 | } |
345 | } | |
346 | } | |
347 | } | |
348 | } else { /* standard case: same event, then mix*/ | |
349 | Correlate(*sTracks, *sTracks, kSameEvt); | |
350 | Int_t nMix = pool->GetCurrentNEvents(); | |
351 | for (Int_t jMix=0; jMix<nMix; ++jMix) { | |
352 | MiniEvent* bgTracks = pool->GetEvent(jMix); | |
d688e049 | 353 | Correlate(*sTracks, *bgTracks, kDiffEvt); |
48a61f36 | 354 | } |
355 | } | |
356 | ||
48a61f36 | 357 | pool->UpdatePool(sTracks); |
358 | PostData(1, fOutputList); | |
48a61f36 | 359 | } |
360 | ||
361 | //________________________________________________________________________ | |
f6701c6e | 362 | void AliDhcTask::GetESDTracks(MiniEvent* miniEvt) |
48a61f36 | 363 | { |
364 | // Loop twice: 1. Count sel. tracks. 2. Fill vector. | |
365 | ||
d688e049 | 366 | if (fTracksName.IsNull()) { |
367 | const AliESDVertex *vtxSPD = fESD->GetPrimaryVertexSPD(); | |
368 | if (!vtxSPD) | |
369 | return; | |
370 | ||
371 | Int_t nTrax = fESD->GetNumberOfTracks(); | |
372 | if (fVerbosity > 2) | |
373 | AliInfo(Form("%d tracks in event",nTrax)); | |
374 | ||
375 | // Loop 1. | |
376 | Int_t nSelTrax = 0; | |
377 | TObjArray arr(nTrax); | |
378 | arr.SetOwner(1); | |
379 | ||
380 | for (Int_t i = 0; i < nTrax; ++i) { | |
381 | AliESDtrack* esdtrack = fESD->GetTrack(i); | |
382 | if (!esdtrack) { | |
383 | AliError(Form("Couldn't get ESD track %d\n", i)); | |
384 | continue; | |
385 | } | |
386 | Bool_t trkOK = fEsdTrackCutsTPCOnly->AcceptTrack(esdtrack); | |
387 | if (!trkOK) | |
388 | continue; | |
389 | Double_t pt = esdtrack->Pt(); | |
390 | Bool_t ptOK = pt >= fPtMin && pt < fPtMax; | |
391 | if (!ptOK) | |
392 | continue; | |
393 | Double_t eta = esdtrack->Eta(); | |
394 | if (TMath::Abs(eta) > fEtaMax) | |
395 | continue; | |
396 | ||
397 | // create a tpc only track | |
398 | AliESDtrack *newtrack = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID()); | |
399 | if(!newtrack) | |
400 | continue; | |
401 | if (newtrack->Pt()<=0) { | |
402 | delete newtrack; | |
403 | continue; | |
404 | } | |
8433ff41 | 405 | |
d688e049 | 406 | AliExternalTrackParam exParam; |
407 | Bool_t relate = newtrack->RelateToVertexTPC(vtxSPD,fESD->GetMagneticField(),kVeryBig,&exParam); | |
408 | if (!relate) { | |
409 | delete newtrack; | |
410 | continue; | |
411 | } | |
8433ff41 | 412 | |
d688e049 | 413 | // set the constraint parameters to the track |
414 | newtrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance()); | |
8433ff41 | 415 | |
d688e049 | 416 | pt = newtrack->Pt(); |
417 | ptOK = pt >= fPtMin && pt < fPtMax; | |
418 | if (!ptOK) { | |
419 | delete newtrack; | |
420 | continue; | |
421 | } | |
422 | eta = esdtrack->Eta(); | |
423 | if (TMath::Abs(eta) > fEtaMax) { | |
424 | delete newtrack; | |
425 | continue; | |
426 | } | |
427 | arr.Add(newtrack); | |
428 | nSelTrax++; | |
8433ff41 | 429 | } |
d688e049 | 430 | return; |
48a61f36 | 431 | } |
432 | ||
d688e049 | 433 | TList *list = InputEvent()->GetList(); |
434 | TClonesArray *tcaTracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracksName)); | |
435 | const Int_t ntracks = tcaTracks->GetEntries(); | |
f6701c6e | 436 | if (miniEvt) |
d688e049 | 437 | miniEvt->reserve(ntracks); |
f6701c6e | 438 | else { |
d688e049 | 439 | AliError("Ptr to miniEvt zero"); |
f6701c6e | 440 | return; |
441 | } | |
48a61f36 | 442 | |
d688e049 | 443 | for(Int_t itrack = 0; itrack < ntracks; itrack++) { |
444 | AliVTrack *esdtrack = static_cast<AliESDtrack*>(tcaTracks->At(itrack)); | |
445 | if(!esdtrack) { | |
446 | AliError(Form("ERROR: Could not retrieve esdtrack %d",itrack)); | |
48a61f36 | 447 | continue; |
448 | } | |
48a61f36 | 449 | Double_t pt = esdtrack->Pt(); |
48a61f36 | 450 | Double_t eta = esdtrack->Eta(); |
48a61f36 | 451 | Double_t phi = esdtrack->Phi(); |
452 | Int_t sign = esdtrack->Charge() > 0 ? 1 : -1; | |
8433ff41 | 453 | miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign)); |
48a61f36 | 454 | } |
48a61f36 | 455 | } |
456 | ||
457 | //________________________________________________________________________ | |
f6701c6e | 458 | void AliDhcTask::GetAODTracks(MiniEvent* miniEvt) |
48a61f36 | 459 | { |
460 | // Loop twice: 1. Count sel. tracks. 2. Fill vector. | |
461 | ||
462 | Int_t nTrax = fAOD->GetNumberOfTracks(); | |
463 | Int_t nSelTrax = 0; | |
464 | ||
465 | if (fVerbosity > 2) | |
466 | AliInfo(Form("%d tracks in event",nTrax)); | |
467 | ||
468 | // Loop 1. | |
469 | for (Int_t i = 0; i < nTrax; ++i) { | |
470 | AliAODTrack* aodtrack = fAOD->GetTrack(i); | |
471 | if (!aodtrack) { | |
472 | AliError(Form("Couldn't get AOD track %d\n", i)); | |
473 | continue; | |
474 | } | |
475 | // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C | |
476 | UInt_t tpcOnly = 1 << 7; | |
477 | Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly); | |
478 | if (!trkOK) | |
479 | continue; | |
480 | Double_t pt = aodtrack->Pt(); | |
481 | Bool_t ptOK = pt >= fPtMin && pt < fPtMax; | |
482 | if (!ptOK) | |
483 | continue; | |
484 | Double_t eta = aodtrack->Eta(); | |
8433ff41 | 485 | if (TMath::Abs(eta) > fEtaMax) |
48a61f36 | 486 | continue; |
487 | nSelTrax++; | |
488 | } | |
489 | ||
f6701c6e | 490 | if (miniEvt) |
491 | miniEvt->reserve(nSelTrax); | |
492 | else { | |
493 | AliError("!miniEvt"); | |
494 | return; | |
495 | } | |
496 | ||
48a61f36 | 497 | // Loop 2. |
498 | for (Int_t i = 0; i < nTrax; ++i) { | |
499 | AliAODTrack* aodtrack = fAOD->GetTrack(i); | |
500 | if (!aodtrack) { | |
501 | AliError(Form("Couldn't get AOD track %d\n", i)); | |
502 | continue; | |
503 | } | |
504 | ||
505 | // See $ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C | |
506 | UInt_t tpcOnly = 1 << 7; | |
507 | Bool_t trkOK = aodtrack->TestFilterBit(tpcOnly); | |
508 | if (!trkOK) | |
509 | continue; | |
510 | Double_t pt = aodtrack->Pt(); | |
511 | Bool_t ptOK = pt >= fPtMin && pt < fPtMax; | |
512 | if (!ptOK) | |
513 | continue; | |
514 | Double_t eta = aodtrack->Eta(); | |
8433ff41 | 515 | if (TMath::Abs(eta) > fEtaMax) |
48a61f36 | 516 | continue; |
517 | ||
518 | Double_t phi = aodtrack->Phi(); | |
519 | Int_t sign = aodtrack->Charge() > 0 ? 1 : -1; | |
8433ff41 | 520 | miniEvt->push_back(AliMiniTrack(pt, eta, phi, sign)); |
48a61f36 | 521 | } |
48a61f36 | 522 | } |
523 | ||
524 | //________________________________________________________________________ | |
525 | Double_t AliDhcTask::DeltaPhi(Double_t phia, Double_t phib, | |
526 | Double_t rangeMin, Double_t rangeMax) const | |
527 | { | |
528 | Double_t dphi = -999; | |
529 | Double_t pi = TMath::Pi(); | |
530 | ||
531 | if (phia < 0) phia += 2*pi; | |
532 | else if (phia > 2*pi) phia -= 2*pi; | |
533 | if (phib < 0) phib += 2*pi; | |
534 | else if (phib > 2*pi) phib -= 2*pi; | |
535 | dphi = phib - phia; | |
536 | if (dphi < rangeMin) dphi += 2*pi; | |
537 | else if (dphi > rangeMax) dphi -= 2*pi; | |
538 | ||
539 | return dphi; | |
540 | } | |
541 | ||
542 | //________________________________________________________________________ | |
d688e049 | 543 | Int_t AliDhcTask::Correlate(const MiniEvent &evt1, const MiniEvent &evt2, Int_t pairing) |
48a61f36 | 544 | { |
545 | // Triggered angular correlations. If pairing is kSameEvt, particles | |
546 | // within evt1 are correlated. If kDiffEvt, correlate triggers from | |
547 | // evt1 with partners from evt2. | |
8433ff41 | 548 | |
549 | Int_t cbin = fHCent->FindBin(fCentrality); | |
550 | if (fHCent->IsBinOverflow(cbin) || | |
551 | fHCent->IsBinUnderflow(cbin)) | |
552 | return 0; | |
553 | ||
554 | Int_t zbin = fHZvtx->FindBin(fZVertex); | |
555 | if (fHZvtx->IsBinOverflow(zbin) || | |
556 | fHZvtx->IsBinUnderflow(zbin)) | |
557 | return 0; | |
558 | ||
48a61f36 | 559 | Int_t iMax = evt1.size(); |
560 | Int_t jMax = evt2.size(); | |
561 | ||
8433ff41 | 562 | TH2 **hist = fHMs; |
563 | if (pairing == kSameEvt) { | |
564 | hist = fHSs; | |
565 | fHCent->AddBinContent(cbin); | |
566 | fHZvtx->AddBinContent(zbin); | |
567 | } | |
568 | ||
569 | Int_t nZvtx = fHZvtx->GetNbinsX(); | |
570 | Int_t nPtTrig = fHPtTrg->GetNbinsX(); | |
571 | Int_t nPtAssc = fHPtAss->GetNbinsX(); | |
572 | ||
573 | Int_t globIndex = (cbin-1)*nZvtx*nPtTrig*nPtAssc+(zbin-1)*nPtTrig*nPtAssc; | |
48a61f36 | 574 | |
d688e049 | 575 | Double_t weight = 1; |
576 | if (fDoWeights) { // Count trigger particles in this event | |
577 | for (Int_t i=0; i<iMax; ++i) { | |
578 | const AliMiniTrack &a(evt1.at(i)); | |
579 | Float_t pta = a.Pt(); | |
580 | Int_t abin = fHPtTrg->FindBin(pta); | |
581 | if (fHPtTrg->IsBinOverflow(abin) || | |
582 | fHPtTrg->IsBinUnderflow(abin)) | |
583 | continue; | |
584 | ++weight; | |
585 | } | |
586 | } | |
587 | ||
48a61f36 | 588 | for (Int_t i=0; i<iMax; ++i) { |
589 | ||
590 | // Trigger particles | |
8433ff41 | 591 | const AliMiniTrack &a(evt1.at(i)); |
592 | ||
593 | Float_t pta = a.Pt(); | |
594 | Int_t abin = fHPtTrg->FindBin(pta); | |
595 | if (fHPtTrg->IsBinOverflow(abin) || | |
596 | fHPtTrg->IsBinUnderflow(abin)) | |
597 | continue; | |
48a61f36 | 598 | |
599 | if (pairing == kSameEvt) { | |
8433ff41 | 600 | fHistPt->Fill(pta); |
48a61f36 | 601 | fHTrk->Fill(a.Phi(),a.Eta()); |
8433ff41 | 602 | fHPtTrg->AddBinContent(abin); |
48a61f36 | 603 | } |
604 | ||
8433ff41 | 605 | for (Int_t j=0; j<jMax; ++j) { |
48a61f36 | 606 | // Associated particles |
607 | if (pairing == kSameEvt && i==j) | |
608 | continue; | |
609 | ||
8433ff41 | 610 | const AliMiniTrack &b(evt2.at(j)); |
48a61f36 | 611 | |
48a61f36 | 612 | Float_t ptb = b.Pt(); |
613 | if (pta < ptb) | |
614 | continue; | |
615 | ||
8433ff41 | 616 | Int_t bbin = fHPtTrg->FindBin(ptb); |
617 | if (fHPtAss->IsBinOverflow(bbin) || | |
618 | fHPtAss->IsBinUnderflow(bbin)) | |
619 | continue; | |
620 | ||
48a61f36 | 621 | Float_t dphi = DeltaPhi(a.Phi(), b.Phi()); |
622 | Float_t deta = a.Eta() - b.Eta(); | |
623 | ||
8433ff41 | 624 | Int_t index = globIndex+(abin-1)*nPtAssc+(bbin-1); |
d688e049 | 625 | hist[index]->Fill(dphi,deta,1./weight); |
8433ff41 | 626 | |
627 | if (pairing == kSameEvt) { | |
628 | fHPtAss->AddBinContent(bbin); | |
629 | fMeanPtTrg[cbin-1]->Fill(pta,ptb,pta); | |
630 | fMeanPtAss[cbin-1]->Fill(pta,ptb,ptb); | |
631 | fMean2PtTrg[cbin-1]->Fill(pta,ptb,pta*pta); | |
632 | fMean2PtAss[cbin-1]->Fill(pta,ptb,ptb*ptb); | |
633 | } | |
48a61f36 | 634 | } |
635 | } | |
636 | ||
637 | return 0; | |
638 | } | |
639 | ||
640 | //________________________________________________________________________ | |
641 | void AliDhcTask::Terminate(Option_t *) | |
642 | { | |
643 | // Draw result to the screen | |
644 | // Called once at the end of the query | |
645 | ||
8433ff41 | 646 | delete fPoolMgr; |
647 | ||
648 | fHCent->SetEntries(fHCent->Integral()); | |
649 | fHZvtx->SetEntries(fHZvtx->Integral()); | |
650 | fHPtTrg->SetEntries(fHPtTrg->Integral()); | |
651 | fHPtAss->SetEntries(fHPtAss->Integral()); | |
652 | ||
48a61f36 | 653 | if (gROOT->IsBatch()) |
654 | return; | |
655 | ||
656 | fOutputList = dynamic_cast<TList*> (GetOutputData(1)); | |
657 | if (!fOutputList) { | |
658 | AliError("Output list not available"); | |
659 | return; | |
660 | } | |
661 | ||
c71ff6f4 | 662 | fHistPt = dynamic_cast<TH1F*> (fOutputList->FindObject("fHistPt")); |
48a61f36 | 663 | if (!fHistPt) { |
664 | AliError("ERROR: fHistPt not available\n"); | |
665 | return; | |
666 | } | |
48a61f36 | 667 | TCanvas *c1 = new TCanvas("AliDhcTask","Pt",10,10,510,510); |
668 | c1->cd(1)->SetLogy(); | |
669 | fHistPt->DrawCopy("E"); | |
670 | } | |
671 | ||
672 | //________________________________________________________________________ | |
673 | Bool_t AliDhcTask::VertexOk(TObject* obj) const | |
674 | { | |
675 | // Modified from AliAnalyseLeadingTrackUE::VertexSelection() | |
676 | ||
677 | Int_t nContributors = 0; | |
678 | Double_t zVertex = 999; | |
679 | TString name(""); | |
680 | ||
681 | if (obj->InheritsFrom("AliESDEvent")) { | |
682 | AliESDEvent* esdevt = (AliESDEvent*) obj; | |
683 | const AliESDVertex* vtx = esdevt->GetPrimaryVertex(); | |
684 | if (!vtx) | |
685 | return 0; | |
686 | nContributors = vtx->GetNContributors(); | |
687 | zVertex = vtx->GetZ(); | |
688 | name = vtx->GetName(); | |
689 | } | |
690 | else if (obj->InheritsFrom("AliAODEvent")) { | |
691 | AliAODEvent* aodevt = (AliAODEvent*) obj; | |
692 | if (aodevt->GetNumberOfVertices() < 1) | |
693 | return 0; | |
694 | const AliAODVertex* vtx = aodevt->GetPrimaryVertex(); | |
695 | nContributors = vtx->GetNContributors(); | |
696 | zVertex = vtx->GetZ(); | |
697 | name = vtx->GetName(); | |
698 | } | |
699 | ||
700 | // Reject if TPC-only vertex | |
8433ff41 | 701 | if (name.CompareTo("TPCVertex")==0) |
48a61f36 | 702 | return kFALSE; |
703 | ||
704 | // Check # contributors and range... | |
705 | if( nContributors < 1 || TMath::Abs(zVertex) > fZVtxMax ) { | |
706 | return kFALSE; | |
707 | } | |
708 | ||
709 | return kTRUE; | |
710 | } |