]>
Commit | Line | Data |
---|---|---|
70f2ce9d | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id$ */ | |
17 | //------------------------------------------------------------------------- | |
18 | // Analysis Task that uses the Soft K-Means Algorithm to find clusters in | |
19 | // the eta-phi space of Minimum Bias. No pt information is used for the clustering. | |
20 | // | |
21 | // | |
22 | // Author: Andreas Morsch (CERN) | |
23 | // andreas.morsch@cern.ch | |
24 | //------------------------------------------------------------------------- | |
25 | ||
26 | ||
27 | ||
28 | #include "TChain.h" | |
29 | #include "TTree.h" | |
30 | #include "TH1F.h" | |
31 | #include "TH2F.h" | |
32 | #include "TProfile.h" | |
33 | ||
34 | #include "TList.h" | |
35 | #include "TParticle.h" | |
36 | #include "TParticlePDG.h" | |
37 | #include "TProfile.h" | |
38 | #include "TMath.h" | |
39 | #include "TRandom.h" | |
40 | ||
41 | #include "AliAnalysisTask.h" | |
42 | #include "AliAnalysisManager.h" | |
43 | ||
42dc9410 | 44 | #include "AliVEvent.h" |
70f2ce9d | 45 | #include "AliESDEvent.h" |
46 | #include "AliStack.h" | |
47 | #include "AliESDVertex.h" | |
48 | #include "AliESDInputHandler.h" | |
49 | #include "AliESDtrackCuts.h" | |
50 | #include "AliMultiplicity.h" | |
51 | ||
52 | #include "AliMCParticle.h" | |
53 | #include "AliMCEvent.h" | |
54 | #include "AliAnalysisTaskKMeans.h" | |
55 | #include "AliTrackReference.h" | |
56 | #include "AliStack.h" | |
57 | #include "AliHeader.h" | |
58 | #include "AliKMeansClustering.h" | |
59 | ||
60 | ||
61 | ||
62 | ClassImp(AliAnalysisTaskKMeans) | |
63 | ||
64 | AliAnalysisTaskKMeans::AliAnalysisTaskKMeans() | |
320bb308 | 65 | : AliAnalysisTaskSE() |
66 | ,fK(0) | |
42dc9410 | 67 | ,fNMin(0) |
70f2ce9d | 68 | ,fHists(0) |
69 | ,fH1CEta(0) | |
70 | ,fH1CPhi(0) | |
71 | ,fH1CEtaR(0) | |
72 | ,fH2N1N2(0) | |
73 | ,fH1Pt(0) | |
74 | ,fH1PtC(0) | |
75 | ,fH1PtC1(0) | |
76 | ,fH1PtC2(0) | |
42dc9410 | 77 | ,fH1PtAS(0) |
78 | ,fH1PtR(0) | |
a51f90df | 79 | ,fH1SPt(0) |
80 | ,fH1SPtC(0) | |
70f2ce9d | 81 | ,fH1DPhi(0) |
82 | ,fH1DR(0) | |
83 | ,fH1DRR(0) | |
84 | ,fH2DPhiEta(0) | |
85 | ,fH2DPhiEtaR(0) | |
86 | ,fH2DPhiEtaL(0) | |
5a9b5891 | 87 | ,fH2DPhiEtaLR(0) |
a51f90df | 88 | ,fH2DPhiEtaC(0) |
89 | ,fH2DPhiEtaCR(0) | |
320bb308 | 90 | ,fH1Resp(0) |
91 | ,fH1RespR(0) | |
42dc9410 | 92 | ,fH2Sigma(0) |
93 | ,fH2SigmaR(0) | |
70f2ce9d | 94 | ,fCuts(0) |
95 | { | |
96 | // | |
97 | // Constructor | |
98 | // | |
99 | } | |
100 | ||
101 | //________________________________________________________________________ | |
102 | AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name) | |
103 | : AliAnalysisTaskSE(name) | |
320bb308 | 104 | ,fK(0) |
42dc9410 | 105 | ,fNMin(0) |
70f2ce9d | 106 | ,fHists(0) |
107 | ,fH1CEta(0) | |
108 | ,fH1CPhi(0) | |
109 | ,fH1CEtaR(0) | |
110 | ,fH2N1N2(0) | |
111 | ,fH1Pt(0) | |
112 | ,fH1PtC(0) | |
113 | ,fH1PtC1(0) | |
114 | ,fH1PtC2(0) | |
42dc9410 | 115 | ,fH1PtAS(0) |
116 | ,fH1PtR(0) | |
a51f90df | 117 | ,fH1SPt(0) |
118 | ,fH1SPtC(0) | |
70f2ce9d | 119 | ,fH1DPhi(0) |
120 | ,fH1DR(0) | |
121 | ,fH1DRR(0) | |
122 | ,fH2DPhiEta(0) | |
123 | ,fH2DPhiEtaR(0) | |
124 | ,fH2DPhiEtaL(0) | |
5a9b5891 | 125 | ,fH2DPhiEtaLR(0) |
a51f90df | 126 | ,fH2DPhiEtaC(0) |
127 | ,fH2DPhiEtaCR(0) | |
320bb308 | 128 | ,fH1Resp(0) |
129 | ,fH1RespR(0) | |
42dc9410 | 130 | ,fH2Sigma(0) |
131 | ,fH2SigmaR(0) | |
70f2ce9d | 132 | ,fCuts(0) |
133 | { | |
134 | // | |
135 | // Constructor | |
136 | // | |
137 | DefineOutput(1, TList::Class()); | |
138 | } | |
139 | ||
140 | ||
141 | //________________________________________________________________________ | |
142 | void AliAnalysisTaskKMeans::UserCreateOutputObjects() | |
143 | { | |
144 | // Create histograms | |
145 | // Called once | |
146 | fHists = new TList(); | |
147 | fH1CEta = new TH1F("fH1CEta", "eta distribution of clusters", 90, -1.5, 1.5); | |
148 | fH1CEtaR = new TH1F("fH1CEtaR", "eta distribution of clusters", 90, -1.5, 1.5); | |
149 | fH1CPhi = new TH1F("fH1CPhi", "phi distribution of clusters", 157, 0.0, 2. * TMath::Pi()); | |
150 | fH2N1N2 = new TH2F("fH2N1N2", "multiplicity distribution", 50, 0., 50., 50, 0., 50.); | |
151 | ||
152 | fH1Pt = new TH1F("fH1Pt", "pt distribution",50, 0., 10.); | |
153 | fH1PtC = new TH1F("fH1PtC", "pt distribution",50, 0., 10.); | |
154 | fH1PtC1 = new TH1F("fH1PtC1", "pt distribution",50, 0., 10.); | |
155 | fH1PtC2 = new TH1F("fH1PtC2", "pt distribution",50, 0., 10.); | |
42dc9410 | 156 | fH1PtAS = new TH1F("fH1PtAS", "pt distribution",50, 0., 10.); |
157 | fH1PtR = new TH1F("fH1PtR", "pt distribution",50, 0., 10.); | |
70f2ce9d | 158 | |
a51f90df | 159 | fH1SPt = new TH1F("fH1SPt", "sum pt distribution",50, 0., 10.); |
160 | fH1SPtC = new TH1F("fH1SPtC", "sum pt distribution",50, 0., 10.); | |
161 | ||
5a9b5891 | 162 | fH1DR = new TH1F("fH1DR", "dR distribution", 50, 0., 5.); |
163 | fH1DPhi = new TH1F("fH1DPhi", "dPhi distribution", 31, 0., TMath::Pi()); | |
164 | fH2DPhiEta = new TH2F("fH2DPhiEta","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.); | |
165 | fH2DPhiEtaR = new TH2F("fH2DPhiEtaR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.); | |
166 | fH2DPhiEtaL = new TH2F("fH2DPhiEtaL","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.); | |
167 | fH2DPhiEtaLR = new TH2F("fH2DPhiEtaLR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.); | |
168 | fH2DPhiEtaC = new TH2F("fH2DPhiEtaC","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.); | |
169 | fH2DPhiEtaCR = new TH2F("fH2DPhiEtaCR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.); | |
170 | fH1DRR = new TH1F("fH1DRR", "dR distribution", 50, 0., 5.); | |
320bb308 | 171 | fH1Resp = new TH1F("fH1Resp", "Responsibility", 50, 0., 1.); |
172 | fH1RespR = new TH1F("fH1RespR", "Responsibility", 50, 0., 1.); | |
42dc9410 | 173 | fH2Sigma = new TH2F("fH2Sigma", "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.); |
174 | fH2SigmaR = new TH2F("fH2SigmaR", "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.); | |
175 | ||
70f2ce9d | 176 | fHists->SetOwner(); |
177 | ||
178 | fHists->Add(fH1CEta); | |
179 | fHists->Add(fH1CEtaR); | |
180 | fHists->Add(fH1CPhi); | |
181 | fHists->Add(fH2N1N2); | |
182 | fHists->Add(fH1Pt); | |
42dc9410 | 183 | fHists->Add(fH1PtR); |
70f2ce9d | 184 | fHists->Add(fH1PtC); |
185 | fHists->Add(fH1PtC1); | |
186 | fHists->Add(fH1PtC2); | |
42dc9410 | 187 | fHists->Add(fH1PtAS); |
70f2ce9d | 188 | fHists->Add(fH1DR); |
a51f90df | 189 | fHists->Add(fH1SPtC); |
190 | fHists->Add(fH1SPt); | |
70f2ce9d | 191 | fHists->Add(fH1DPhi); |
192 | fHists->Add(fH2DPhiEta); | |
193 | fHists->Add(fH2DPhiEtaR); | |
194 | fHists->Add(fH2DPhiEtaL); | |
5a9b5891 | 195 | fHists->Add(fH2DPhiEtaLR); |
a51f90df | 196 | fHists->Add(fH2DPhiEtaC); |
197 | fHists->Add(fH2DPhiEtaCR); | |
70f2ce9d | 198 | fHists->Add(fH1DRR); |
320bb308 | 199 | fHists->Add(fH1RespR); |
200 | fHists->Add(fH1Resp); | |
42dc9410 | 201 | fHists->Add(fH2Sigma); |
202 | fHists->Add(fH2SigmaR); | |
203 | ||
70f2ce9d | 204 | // |
70f2ce9d | 205 | |
206 | } | |
207 | ||
208 | //________________________________________________________________________ | |
209 | void AliAnalysisTaskKMeans::UserExec(Option_t *) | |
210 | { | |
211 | // Main loop | |
212 | // Called for each event | |
213 | Double_t phi [500]; | |
214 | Double_t phiR[500]; | |
215 | Double_t eta[500]; | |
216 | Double_t etaR[500]; | |
217 | Double_t pt [500]; | |
218 | Double_t mPhi[20]; | |
219 | Double_t mEta[20]; | |
220 | Double_t rk[20]; | |
42dc9410 | 221 | Double_t sigma2[20]; |
222 | Double_t sigmax2[20]; | |
223 | Double_t sigmay2[20]; | |
70f2ce9d | 224 | Int_t ind[20]; |
225 | ||
226 | if (!fInputEvent) { | |
227 | Printf("ERROR: fESD not available"); | |
228 | return; | |
229 | } | |
230 | ||
70f2ce9d | 231 | |
232 | Int_t ic = 0; | |
233 | // | |
234 | // Fill eta-phi positions | |
235 | // | |
236 | /* | |
42dc9410 | 237 | const AliMultiplicity *spdMult = fInputEvent->GetMultiplicity(); |
70f2ce9d | 238 | for (Int_t i = 0; i < spdMult->GetNumberOfTracklets(); ++i) { |
239 | phi[ic] = spdMult->GetPhi(i); | |
240 | eta[ic] = spdMult->GetEta(i); | |
241 | ic++; | |
242 | } | |
243 | */ | |
244 | Double_t ptmax = 0.; | |
245 | Int_t icmax = -1; | |
246 | ||
42dc9410 | 247 | for (Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++) { |
248 | AliVParticle* track = fInputEvent->GetTrack(iTracks); | |
249 | if ((fCuts->AcceptTrack((AliESDtrack*)track))) | |
70f2ce9d | 250 | { |
251 | phi[ic] = track->Phi(); | |
252 | phiR[ic] = 2. * TMath::Pi() * gRandom->Rndm(); | |
253 | eta[ic] = track->Eta(); | |
254 | etaR[ic] = 2. * gRandom->Rndm() - 1.; | |
255 | pt[ic] = track->Pt(); | |
256 | if (pt[ic] > ptmax) { | |
257 | ptmax = pt[ic]; | |
258 | icmax = ic; | |
259 | } | |
260 | ic++; | |
261 | } | |
262 | } //track loop | |
263 | ||
264 | // | |
265 | // Cluster | |
42dc9410 | 266 | if (ic < fNMin) { |
70f2ce9d | 267 | PostData(1, fHists); |
268 | return; | |
269 | } | |
270 | // | |
42dc9410 | 271 | Int_t ok; |
272 | Double_t dphic, detac; | |
a51f90df | 273 | // |
42dc9410 | 274 | ok = AliKMeansClustering::SoftKMeans3(fK, ic, phi, eta, mPhi, mEta, sigmax2, sigmay2, rk); |
275 | if (ok) { | |
276 | // | |
277 | // Sort | |
278 | for (Int_t i = 0; i < fK; i++) { | |
279 | if (rk[i] > 1.) rk[i] /= TMath::Sqrt(sigmax2[i] * sigmay2[i]); | |
280 | else rk[i] = 0.; | |
281 | } | |
282 | ||
283 | TMath::Sort(fK, rk, ind); | |
284 | fH1Resp->Fill(rk[ind[0]]/(rk[0]+rk[1])); | |
285 | fH2Sigma->Fill(TMath::Sqrt(sigmax2[ind[0]]), TMath::Sqrt(sigmay2[ind[0]])); | |
286 | // | |
287 | // Analysis | |
288 | // | |
289 | // Cluster Multiplicity | |
290 | Int_t mult[2] = {0, 0}; | |
291 | // | |
292 | dphic = DeltaPhi(mPhi[0], mPhi[1]); | |
293 | detac = TMath::Abs(mEta[0] - mEta[1]); | |
294 | fH2DPhiEtaC->Fill(dphic, detac); | |
295 | // | |
296 | // Random cluster position | |
297 | Int_t ir = Int_t(Float_t(ic) * gRandom->Rndm()); | |
a51f90df | 298 | |
42dc9410 | 299 | Double_t crPhi = phi[ir]; |
300 | Double_t crEta = eta[ir]; | |
301 | // | |
302 | Double_t sumPt = 0; | |
303 | Double_t sumPtC = 0; | |
304 | for (Int_t i = 0; i < 1; i++) { | |
70f2ce9d | 305 | fH1CEta->Fill(mEta[ind[i]]); |
306 | fH1CPhi->Fill(mPhi[ind[i]]); | |
307 | for (Int_t j = 0; j < ic; j++) { | |
42dc9410 | 308 | Double_t r = DeltaR(mPhi[ind[i]], mEta[ind[i]], phi[j], eta[j]); |
309 | Double_t dphi = DeltaPhi(mPhi[ind[i]], phi[j]); | |
310 | Double_t deta = mEta[ind[i]] - eta[j]; | |
311 | Double_t rr = DeltaR(crPhi, crEta, phi[j], eta[j]); | |
312 | fH1DR->Fill(r); | |
313 | fH1DPhi->Fill(dphi); | |
314 | fH2DPhiEta->Fill(dphi, TMath::Abs(deta)); | |
315 | if (j == icmax) fH2DPhiEtaL->Fill(dphi, TMath::Abs(deta)); | |
316 | ||
317 | if (r < 0.2) { | |
318 | fH1PtC2->Fill(pt[j]); | |
319 | } | |
320 | if (r < 0.3) | |
321 | { | |
322 | fH1PtC1->Fill(pt[j]); | |
70f2ce9d | 323 | } |
42dc9410 | 324 | if (rr < 0.3) |
70f2ce9d | 325 | { |
42dc9410 | 326 | fH1PtR->Fill(pt[j]); |
70f2ce9d | 327 | } |
42dc9410 | 328 | |
329 | if (r < 0.4) | |
70f2ce9d | 330 | { |
a51f90df | 331 | sumPtC += pt[j]; |
332 | mult[i]++; | |
333 | fH1PtC->Fill(pt[j]); | |
70f2ce9d | 334 | } |
42dc9410 | 335 | if (r > 0.7) { |
336 | fH1Pt->Fill(pt[j]); | |
337 | } | |
338 | ||
339 | if (r > 0.7 && r < 1.) { | |
340 | sumPt += pt[j]; | |
341 | } | |
342 | ||
343 | if (dphi > (TMath::Pi() - 0.3)) { | |
344 | fH1PtAS->Fill(pt[j]); | |
345 | } | |
70f2ce9d | 346 | } |
42dc9410 | 347 | } |
70f2ce9d | 348 | |
42dc9410 | 349 | fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1])); |
350 | fH1SPt ->Fill(sumPt); | |
351 | fH1SPtC->Fill(sumPtC); | |
352 | } // ok | |
70f2ce9d | 353 | |
354 | // Randomized phi | |
355 | // | |
42dc9410 | 356 | ok = AliKMeansClustering::SoftKMeans3(fK, ic, phiR, etaR, mPhi, mEta, sigmax2, sigmay2, rk); |
357 | if (ok) { | |
358 | // | |
359 | // Sort according to responsibility density | |
360 | for (Int_t i = 0; i < fK; i++) { | |
361 | if (rk[i] > 1.) rk[i] /= TMath::Sqrt(sigmax2[i] * sigmay2[i]); | |
362 | else rk[i] = 0.; | |
363 | } | |
364 | ||
365 | TMath::Sort(fK, rk, ind); | |
366 | fH1RespR->Fill(rk[ind[0]]/(rk[0]+rk[1])); | |
367 | fH2SigmaR->Fill(TMath::Sqrt(sigmax2[ind[0]]), TMath::Sqrt(sigmay2[ind[0]])); | |
368 | // | |
369 | // Analyse | |
370 | // | |
371 | // Cluster Multiplicity | |
372 | for (Int_t i = 0; i < 1; i++) { | |
70f2ce9d | 373 | fH1CEtaR->Fill(mEta[ind[i]]); |
42dc9410 | 374 | } |
375 | ||
376 | for (Int_t i = 0; i < 1; i++) { | |
70f2ce9d | 377 | for (Int_t j = 0; j < ic; j++) { |
42dc9410 | 378 | Double_t dphi = DeltaPhi(mPhi[ind[i]], phiR[j]); |
379 | Double_t deta = mEta[ind[i]] - etaR[j]; | |
380 | Double_t r = DeltaR(mPhi[ind[i]], mEta[ind[i]], phiR[j], etaR[j]); | |
381 | fH1DRR->Fill(r); | |
382 | fH2DPhiEtaR->Fill(dphi, TMath::Abs(deta)); | |
383 | if (j == icmax) fH2DPhiEtaLR->Fill(dphi, TMath::Abs(deta)); | |
70f2ce9d | 384 | } |
42dc9410 | 385 | } |
386 | dphic = DeltaPhi(mPhi[0], mPhi[1]); | |
387 | detac = TMath::Abs(mEta[0] - mEta[1]); | |
388 | fH2DPhiEtaCR->Fill(dphic, detac); | |
389 | } // ok | |
70f2ce9d | 390 | |
391 | // | |
392 | // Post output data. | |
393 | PostData(1, fHists); | |
394 | } | |
395 | ||
396 | Double_t AliAnalysisTaskKMeans::DeltaPhi(Double_t phi1, Double_t phi2) | |
397 | { | |
398 | Double_t dphi = TMath::Abs(phi1 - phi2); | |
399 | if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi; | |
400 | return dphi; | |
401 | } | |
402 | ||
403 | Double_t AliAnalysisTaskKMeans::DeltaR(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2) | |
404 | { | |
405 | Double_t dphi = DeltaPhi(phi1, phi2); | |
406 | Double_t deta = eta1 - eta2; | |
407 | return (TMath::Sqrt(dphi * dphi + deta * deta)); | |
408 | ||
409 | } | |
410 | ||
411 | //________________________________________________________________________ | |
412 | void AliAnalysisTaskKMeans::Terminate(Option_t *) | |
413 | { | |
414 | } | |
415 |