]>
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$ */ | |
b31787ca | 17 | //--------------------------------------------------------------------------------------- |
70f2ce9d | 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 | |
b31787ca | 24 | //--------------------------------------------------------------------------------------- |
70f2ce9d | 25 | |
26 | ||
27 | ||
28 | #include "TChain.h" | |
b31787ca | 29 | #include "TFile.h" |
70f2ce9d | 30 | #include "TTree.h" |
31 | #include "TH1F.h" | |
32 | #include "TH2F.h" | |
33 | #include "TProfile.h" | |
34 | ||
35 | #include "TList.h" | |
36 | #include "TParticle.h" | |
37 | #include "TParticlePDG.h" | |
38 | #include "TProfile.h" | |
39 | #include "TMath.h" | |
40 | #include "TRandom.h" | |
19c5a36a | 41 | #include "TObjArray.h" |
70f2ce9d | 42 | |
43 | #include "AliAnalysisTask.h" | |
44 | #include "AliAnalysisManager.h" | |
45 | ||
42dc9410 | 46 | #include "AliVEvent.h" |
70f2ce9d | 47 | #include "AliESDEvent.h" |
b31787ca | 48 | #include "AliExternalTrackParam.h" |
70f2ce9d | 49 | #include "AliStack.h" |
50 | #include "AliESDVertex.h" | |
51 | #include "AliESDInputHandler.h" | |
52 | #include "AliESDtrackCuts.h" | |
53 | #include "AliMultiplicity.h" | |
54 | ||
55 | #include "AliMCParticle.h" | |
56 | #include "AliMCEvent.h" | |
57 | #include "AliAnalysisTaskKMeans.h" | |
58 | #include "AliTrackReference.h" | |
59 | #include "AliStack.h" | |
60 | #include "AliHeader.h" | |
61 | #include "AliKMeansClustering.h" | |
62 | ||
63 | ||
64 | ||
65 | ClassImp(AliAnalysisTaskKMeans) | |
66 | ||
67 | AliAnalysisTaskKMeans::AliAnalysisTaskKMeans() | |
320bb308 | 68 | : AliAnalysisTaskSE() |
69 | ,fK(0) | |
42dc9410 | 70 | ,fNMin(0) |
70f2ce9d | 71 | ,fHists(0) |
72 | ,fH1CEta(0) | |
73 | ,fH1CPhi(0) | |
74 | ,fH1CEtaR(0) | |
75 | ,fH2N1N2(0) | |
76 | ,fH1Pt(0) | |
77 | ,fH1PtC(0) | |
78 | ,fH1PtC1(0) | |
79 | ,fH1PtC2(0) | |
42dc9410 | 80 | ,fH1PtAS(0) |
81 | ,fH1PtR(0) | |
a51f90df | 82 | ,fH1SPt(0) |
83 | ,fH1SPtC(0) | |
70f2ce9d | 84 | ,fH1DPhi(0) |
85 | ,fH1DR(0) | |
86 | ,fH1DRR(0) | |
87 | ,fH2DPhiEta(0) | |
88 | ,fH2DPhiEtaR(0) | |
89 | ,fH2DPhiEtaL(0) | |
5a9b5891 | 90 | ,fH2DPhiEtaLR(0) |
a51f90df | 91 | ,fH2DPhiEtaC(0) |
92 | ,fH2DPhiEtaCR(0) | |
320bb308 | 93 | ,fH1Resp(0) |
94 | ,fH1RespR(0) | |
42dc9410 | 95 | ,fH2Sigma(0) |
96 | ,fH2SigmaR(0) | |
19c5a36a | 97 | ,fHDensity(0) |
98 | ,fHCSize(0) | |
99 | ,fHNCluster(0) | |
100 | ,fHPtDensity(0) | |
b31787ca | 101 | ,fHDPhi(0) |
102 | ,fH2EtaPhi(0) | |
103 | ,fH2REtaPhi(0) | |
70f2ce9d | 104 | ,fCuts(0) |
f2b8ec84 | 105 | |
70f2ce9d | 106 | { |
107 | // | |
f34058c0 | 108 | |
70f2ce9d | 109 | // Constructor |
110 | // | |
f34058c0 | 111 | for (Int_t i=0; i< 10; i++) { |
112 | fA[i] = 0; | |
113 | fB[i] = 0; | |
114 | } | |
70f2ce9d | 115 | } |
116 | ||
117 | //________________________________________________________________________ | |
118 | AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name) | |
119 | : AliAnalysisTaskSE(name) | |
320bb308 | 120 | ,fK(0) |
42dc9410 | 121 | ,fNMin(0) |
70f2ce9d | 122 | ,fHists(0) |
123 | ,fH1CEta(0) | |
124 | ,fH1CPhi(0) | |
125 | ,fH1CEtaR(0) | |
126 | ,fH2N1N2(0) | |
127 | ,fH1Pt(0) | |
128 | ,fH1PtC(0) | |
129 | ,fH1PtC1(0) | |
130 | ,fH1PtC2(0) | |
42dc9410 | 131 | ,fH1PtAS(0) |
132 | ,fH1PtR(0) | |
a51f90df | 133 | ,fH1SPt(0) |
134 | ,fH1SPtC(0) | |
70f2ce9d | 135 | ,fH1DPhi(0) |
136 | ,fH1DR(0) | |
137 | ,fH1DRR(0) | |
138 | ,fH2DPhiEta(0) | |
139 | ,fH2DPhiEtaR(0) | |
140 | ,fH2DPhiEtaL(0) | |
5a9b5891 | 141 | ,fH2DPhiEtaLR(0) |
a51f90df | 142 | ,fH2DPhiEtaC(0) |
143 | ,fH2DPhiEtaCR(0) | |
320bb308 | 144 | ,fH1Resp(0) |
145 | ,fH1RespR(0) | |
42dc9410 | 146 | ,fH2Sigma(0) |
147 | ,fH2SigmaR(0) | |
19c5a36a | 148 | ,fHDensity(0) |
149 | ,fHCSize(0) | |
150 | ,fHNCluster(0) | |
151 | ,fHPtDensity(0) | |
b31787ca | 152 | ,fHDPhi(0) |
153 | ,fH2EtaPhi(0) | |
154 | ,fH2REtaPhi(0) | |
70f2ce9d | 155 | ,fCuts(0) |
f2b8ec84 | 156 | |
70f2ce9d | 157 | { |
158 | // | |
159 | // Constructor | |
160 | // | |
f34058c0 | 161 | for (Int_t i=0; i< 10; i++) { |
162 | fA[i] = 0; | |
163 | fB[i] = 0; | |
164 | } | |
70f2ce9d | 165 | DefineOutput(1, TList::Class()); |
166 | } | |
167 | ||
e42693ad | 168 | AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const AliAnalysisTaskKMeans &res) |
169 | : AliAnalysisTaskSE(res) | |
170 | ,fK(0) | |
171 | ,fNMin(0) | |
172 | ,fHists(0) | |
173 | ,fH1CEta(0) | |
174 | ,fH1CPhi(0) | |
175 | ,fH1CEtaR(0) | |
176 | ,fH2N1N2(0) | |
177 | ,fH1Pt(0) | |
178 | ,fH1PtC(0) | |
179 | ,fH1PtC1(0) | |
180 | ,fH1PtC2(0) | |
181 | ,fH1PtAS(0) | |
182 | ,fH1PtR(0) | |
183 | ,fH1SPt(0) | |
184 | ,fH1SPtC(0) | |
185 | ,fH1DPhi(0) | |
186 | ,fH1DR(0) | |
187 | ,fH1DRR(0) | |
188 | ,fH2DPhiEta(0) | |
189 | ,fH2DPhiEtaR(0) | |
190 | ,fH2DPhiEtaL(0) | |
191 | ,fH2DPhiEtaLR(0) | |
192 | ,fH2DPhiEtaC(0) | |
193 | ,fH2DPhiEtaCR(0) | |
194 | ,fH1Resp(0) | |
195 | ,fH1RespR(0) | |
196 | ,fH2Sigma(0) | |
197 | ,fH2SigmaR(0) | |
f2b8ec84 | 198 | ,fHDensity(0) |
199 | ,fHCSize(0) | |
200 | ,fHNCluster(0) | |
201 | ,fHPtDensity(0) | |
b31787ca | 202 | ,fHDPhi(0) |
203 | ,fH2EtaPhi(0) | |
204 | ,fH2REtaPhi(0) | |
e42693ad | 205 | ,fCuts(0) |
206 | { | |
207 | // Dummy copy constructor | |
f34058c0 | 208 | for (Int_t i=0; i< 10; i++) { |
209 | fA[i] = 0; | |
210 | fB[i] = 0; | |
211 | } | |
e42693ad | 212 | } |
213 | ||
214 | AliAnalysisTaskKMeans& AliAnalysisTaskKMeans::operator=(const AliAnalysisTaskKMeans& /*trclass*/) | |
215 | { | |
216 | // Dummy assignment operator | |
217 | return *this; | |
218 | } | |
219 | ||
220 | ||
70f2ce9d | 221 | |
222 | //________________________________________________________________________ | |
223 | void AliAnalysisTaskKMeans::UserCreateOutputObjects() | |
224 | { | |
225 | // Create histograms | |
226 | // Called once | |
227 | fHists = new TList(); | |
228 | fH1CEta = new TH1F("fH1CEta", "eta distribution of clusters", 90, -1.5, 1.5); | |
229 | fH1CEtaR = new TH1F("fH1CEtaR", "eta distribution of clusters", 90, -1.5, 1.5); | |
230 | fH1CPhi = new TH1F("fH1CPhi", "phi distribution of clusters", 157, 0.0, 2. * TMath::Pi()); | |
231 | fH2N1N2 = new TH2F("fH2N1N2", "multiplicity distribution", 50, 0., 50., 50, 0., 50.); | |
232 | ||
233 | fH1Pt = new TH1F("fH1Pt", "pt distribution",50, 0., 10.); | |
234 | fH1PtC = new TH1F("fH1PtC", "pt distribution",50, 0., 10.); | |
235 | fH1PtC1 = new TH1F("fH1PtC1", "pt distribution",50, 0., 10.); | |
236 | fH1PtC2 = new TH1F("fH1PtC2", "pt distribution",50, 0., 10.); | |
42dc9410 | 237 | fH1PtAS = new TH1F("fH1PtAS", "pt distribution",50, 0., 10.); |
238 | fH1PtR = new TH1F("fH1PtR", "pt distribution",50, 0., 10.); | |
70f2ce9d | 239 | |
a51f90df | 240 | fH1SPt = new TH1F("fH1SPt", "sum pt distribution",50, 0., 10.); |
241 | fH1SPtC = new TH1F("fH1SPtC", "sum pt distribution",50, 0., 10.); | |
242 | ||
5a9b5891 | 243 | fH1DR = new TH1F("fH1DR", "dR distribution", 50, 0., 5.); |
244 | fH1DPhi = new TH1F("fH1DPhi", "dPhi distribution", 31, 0., TMath::Pi()); | |
245 | fH2DPhiEta = new TH2F("fH2DPhiEta","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.); | |
246 | fH2DPhiEtaR = new TH2F("fH2DPhiEtaR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.); | |
247 | fH2DPhiEtaL = new TH2F("fH2DPhiEtaL","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.); | |
248 | fH2DPhiEtaLR = new TH2F("fH2DPhiEtaLR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.); | |
249 | fH2DPhiEtaC = new TH2F("fH2DPhiEtaC","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.); | |
250 | fH2DPhiEtaCR = new TH2F("fH2DPhiEtaCR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.); | |
251 | fH1DRR = new TH1F("fH1DRR", "dR distribution", 50, 0., 5.); | |
320bb308 | 252 | fH1Resp = new TH1F("fH1Resp", "Responsibility", 50, 0., 1.); |
253 | fH1RespR = new TH1F("fH1RespR", "Responsibility", 50, 0., 1.); | |
42dc9410 | 254 | fH2Sigma = new TH2F("fH2Sigma", "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.); |
255 | fH2SigmaR = new TH2F("fH2SigmaR", "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.); | |
19c5a36a | 256 | |
257 | ||
258 | fHDensity = new TH1F("fHDensity", "density distribution", 100, 0., 20.); | |
259 | fHCSize = new TH1F("fHCSize", "cluster size", 20, -0.5, 19.5); | |
260 | ||
261 | fHNCluster = new TH1F("fHNCluster", "Number of clusters", 11, -0.5, 10.5); | |
262 | fHPtDensity = new TH2F("fHPtDensity", "Pt vs density", 100, 0., 20., 50, 0., 10.); | |
b31787ca | 263 | fHDPhi = new TH1F("fHDPhi", "phi correlation", 100, 0., 0.5); |
264 | fH2EtaPhi = new TH2F("fH2EtaPhi", "Eta Phi", 200, -1., 1., 628, 0., 2. * TMath::Pi()); | |
70f2ce9d | 265 | fHists->SetOwner(); |
266 | ||
267 | fHists->Add(fH1CEta); | |
268 | fHists->Add(fH1CEtaR); | |
269 | fHists->Add(fH1CPhi); | |
270 | fHists->Add(fH2N1N2); | |
271 | fHists->Add(fH1Pt); | |
42dc9410 | 272 | fHists->Add(fH1PtR); |
70f2ce9d | 273 | fHists->Add(fH1PtC); |
274 | fHists->Add(fH1PtC1); | |
275 | fHists->Add(fH1PtC2); | |
42dc9410 | 276 | fHists->Add(fH1PtAS); |
70f2ce9d | 277 | fHists->Add(fH1DR); |
a51f90df | 278 | fHists->Add(fH1SPtC); |
279 | fHists->Add(fH1SPt); | |
70f2ce9d | 280 | fHists->Add(fH1DPhi); |
281 | fHists->Add(fH2DPhiEta); | |
282 | fHists->Add(fH2DPhiEtaR); | |
283 | fHists->Add(fH2DPhiEtaL); | |
5a9b5891 | 284 | fHists->Add(fH2DPhiEtaLR); |
a51f90df | 285 | fHists->Add(fH2DPhiEtaC); |
286 | fHists->Add(fH2DPhiEtaCR); | |
70f2ce9d | 287 | fHists->Add(fH1DRR); |
320bb308 | 288 | fHists->Add(fH1RespR); |
289 | fHists->Add(fH1Resp); | |
42dc9410 | 290 | fHists->Add(fH2Sigma); |
19c5a36a | 291 | fHists->Add(fH2SigmaR); |
292 | fHists->Add(fHCSize); | |
293 | fHists->Add(fHDensity); | |
294 | fHists->Add(fHNCluster); | |
295 | fHists->Add(fHPtDensity); | |
b31787ca | 296 | fHists->Add(fHDPhi); |
297 | fHists->Add(fH2EtaPhi); | |
70f2ce9d | 298 | // |
19c5a36a | 299 | for (Int_t i = 0; i < 10; i++) { |
300 | fA[i] = new AliKMeansResult(i+1); | |
301 | fB[i] = new AliKMeansResult(i+1); | |
302 | } | |
70f2ce9d | 303 | } |
304 | ||
305 | //________________________________________________________________________ | |
306 | void AliAnalysisTaskKMeans::UserExec(Option_t *) | |
307 | { | |
308 | // Main loop | |
309 | // Called for each event | |
f34058c0 | 310 | Double_t phi [500] = {0}; |
311 | Double_t phiR[500] = {0}; | |
312 | Double_t eta[500] = {0}; | |
313 | Double_t etaR[500] = {0}; | |
314 | Double_t pt [500] = {0}; | |
70f2ce9d | 315 | |
316 | if (!fInputEvent) { | |
317 | Printf("ERROR: fESD not available"); | |
318 | return; | |
319 | } | |
320 | ||
70f2ce9d | 321 | |
322 | Int_t ic = 0; | |
19c5a36a | 323 | |
70f2ce9d | 324 | // |
325 | // Fill eta-phi positions | |
326 | // | |
19c5a36a | 327 | |
70f2ce9d | 328 | Double_t ptmax = 0.; |
329 | Int_t icmax = -1; | |
330 | ||
42dc9410 | 331 | for (Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++) { |
332 | AliVParticle* track = fInputEvent->GetTrack(iTracks); | |
333 | if ((fCuts->AcceptTrack((AliESDtrack*)track))) | |
70f2ce9d | 334 | { |
b31787ca | 335 | const AliExternalTrackParam * tpcT = ((AliESDtrack*) track)->GetTPCInnerParam(); |
336 | if (!tpcT) continue; | |
337 | if (TMath::Abs(tpcT->Eta()) > 0.9) continue; | |
338 | ||
339 | phi [ic] = tpcT->Phi(); | |
340 | eta [ic] = tpcT->Eta(); | |
341 | pt [ic] = tpcT->Pt(); | |
342 | ||
343 | if (fH2REtaPhi) { | |
344 | fH2REtaPhi->GetRandom2(etaR[ic], phiR[ic]); | |
345 | } else { | |
346 | phiR[ic] = 2. * TMath::Pi() * gRandom->Rndm(); | |
347 | etaR[ic] = 1.8 * gRandom->Rndm() - 0.9; | |
348 | } | |
349 | ||
350 | ||
70f2ce9d | 351 | if (pt[ic] > ptmax) { |
352 | ptmax = pt[ic]; | |
353 | icmax = ic; | |
354 | } | |
b31787ca | 355 | |
356 | fH2EtaPhi->Fill(eta[ic], phi[ic]); | |
70f2ce9d | 357 | ic++; |
358 | } | |
359 | } //track loop | |
360 | ||
b31787ca | 361 | for (Int_t i = 0; i < ic; i++) { |
362 | for (Int_t j = i+1; j < ic; j++) { | |
363 | Double_t dphi = TMath::Abs(phi[i] - phi[j]); | |
364 | fHDPhi->Fill(dphi); | |
365 | } | |
366 | } | |
367 | ||
70f2ce9d | 368 | // |
369 | // Cluster | |
42dc9410 | 370 | if (ic < fNMin) { |
70f2ce9d | 371 | PostData(1, fHists); |
372 | return; | |
373 | } | |
b31787ca | 374 | |
70f2ce9d | 375 | // |
43760a73 | 376 | Double_t rk0[10]; |
19c5a36a | 377 | AliKMeansResult* res = 0; |
b31787ca | 378 | AliKMeansResult best(10); |
379 | Float_t rmaxG = -1.; | |
19c5a36a | 380 | |
b31787ca | 381 | for (Int_t k = 0; k < 20; k++) { |
382 | Float_t rmax = -1.; | |
383 | Int_t imax = 0; | |
384 | for (Int_t i = 0; i < fK; i++) { | |
385 | res = fA[i]; | |
386 | AliKMeansClustering::SoftKMeans2(i+1, ic, phi, eta, res->GetMx(),res->GetMy(), res->GetSigma2(), res->GetRk()); | |
387 | res->Sort(ic, phi, eta); | |
388 | Int_t j = (res->GetInd())[0]; | |
389 | rk0[i] = (res->GetTarget())[j]; | |
390 | if (rk0[i] > rmax) { | |
19c5a36a | 391 | rmax = rk0[i]; |
392 | imax = i; | |
393 | } | |
b31787ca | 394 | } |
395 | if (rmax > rmaxG) { | |
396 | rmaxG = rmax; | |
397 | best.CopyResults(fA[imax]); | |
398 | } | |
19c5a36a | 399 | } |
b31787ca | 400 | Double_t* mPhi = best.GetMx(); |
401 | Double_t* mEta = best.GetMy(); | |
402 | Double_t* sigma2 = best.GetSigma2(); | |
403 | Int_t nk = best.GetK(); | |
404 | Int_t im = (best.GetInd())[0]; | |
405 | Double_t etaC = mEta[im]; | |
406 | ||
407 | fHDensity->Fill(rmaxG / TMath::Pi()); | |
408 | fHCSize->Fill(2.28 * rmaxG * sigma2[im]); | |
19c5a36a | 409 | fHNCluster->Fill(Float_t(nk)); |
410 | ||
42dc9410 | 411 | Double_t dphic, detac; |
42dc9410 | 412 | |
b31787ca | 413 | if (rmaxG > 0. && TMath::Abs(etaC) < 0.4) { |
42dc9410 | 414 | // Analysis |
415 | // | |
416 | // Cluster Multiplicity | |
417 | Int_t mult[2] = {0, 0}; | |
418 | // | |
19c5a36a | 419 | if (nk > 1) { |
420 | dphic = DeltaPhi(mPhi[0], mPhi[1]); | |
421 | detac = TMath::Abs(mEta[0] - mEta[1]); | |
422 | fH2DPhiEtaC->Fill(dphic, detac); | |
423 | } | |
42dc9410 | 424 | // |
425 | // Random cluster position | |
426 | Int_t ir = Int_t(Float_t(ic) * gRandom->Rndm()); | |
a51f90df | 427 | |
42dc9410 | 428 | Double_t crPhi = phi[ir]; |
429 | Double_t crEta = eta[ir]; | |
430 | // | |
431 | Double_t sumPt = 0; | |
432 | Double_t sumPtC = 0; | |
19c5a36a | 433 | |
42dc9410 | 434 | for (Int_t i = 0; i < 1; i++) { |
19c5a36a | 435 | fH1CEta->Fill(mEta[im]); |
436 | fH1CPhi->Fill(mPhi[im]); | |
70f2ce9d | 437 | for (Int_t j = 0; j < ic; j++) { |
19c5a36a | 438 | Double_t r = DeltaR(mPhi[im], mEta[im], phi[j], eta[j]); |
439 | Double_t dphi = DeltaPhi(mPhi[im], phi[j]); | |
440 | Double_t deta = mEta[im] - eta[j]; | |
42dc9410 | 441 | Double_t rr = DeltaR(crPhi, crEta, phi[j], eta[j]); |
19c5a36a | 442 | |
42dc9410 | 443 | fH1DR->Fill(r); |
444 | fH1DPhi->Fill(dphi); | |
445 | fH2DPhiEta->Fill(dphi, TMath::Abs(deta)); | |
446 | if (j == icmax) fH2DPhiEtaL->Fill(dphi, TMath::Abs(deta)); | |
447 | ||
448 | if (r < 0.2) { | |
449 | fH1PtC2->Fill(pt[j]); | |
450 | } | |
19c5a36a | 451 | if (r < 0.3) { |
452 | fH1PtC1->Fill(pt[j]); | |
b31787ca | 453 | fHPtDensity->Fill(rmaxG/TMath::Pi(), pt[j]); |
19c5a36a | 454 | } |
455 | if (rr < 0.3) { | |
42dc9410 | 456 | fH1PtR->Fill(pt[j]); |
19c5a36a | 457 | } |
42dc9410 | 458 | |
19c5a36a | 459 | if (r < 0.4) { |
460 | sumPtC += pt[j]; | |
461 | mult[i]++; | |
462 | fH1PtC->Fill(pt[j]); | |
463 | } | |
464 | if (r > 0.7 && dphi < (TMath::Pi() - 0.3)) { | |
42dc9410 | 465 | fH1Pt->Fill(pt[j]); |
466 | } | |
467 | ||
468 | if (r > 0.7 && r < 1.) { | |
469 | sumPt += pt[j]; | |
470 | } | |
19c5a36a | 471 | |
42dc9410 | 472 | if (dphi > (TMath::Pi() - 0.3)) { |
b31787ca | 473 | fH1PtAS->Fill(pt[j], 1.); |
42dc9410 | 474 | } |
70f2ce9d | 475 | } |
42dc9410 | 476 | } |
70f2ce9d | 477 | |
42dc9410 | 478 | fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1])); |
479 | fH1SPt ->Fill(sumPt); | |
480 | fH1SPtC->Fill(sumPtC); | |
19c5a36a | 481 | } |
42dc9410 | 482 | // |
19c5a36a | 483 | // Randomized phi |
42dc9410 | 484 | // |
b31787ca | 485 | rmaxG = -1.; |
486 | for (Int_t k = 0; k < 20; k++) { | |
487 | Float_t rmax = -1.; | |
488 | Int_t imax = 0; | |
19c5a36a | 489 | for (Int_t i = 0; i < fK; i++) { |
490 | res = fB[i]; | |
491 | AliKMeansClustering::SoftKMeans2(i+1, ic, phiR, etaR, res->GetMx(),res->GetMy(), res->GetSigma2(), res->GetRk()); | |
492 | res->Sort(ic, phiR, etaR); | |
19c5a36a | 493 | Int_t j = (res->GetInd())[0]; |
494 | rk0[i] = (res->GetTarget())[j]; | |
19c5a36a | 495 | if (rk0[i] > rmax) { |
b31787ca | 496 | rmax = rk0[i]; |
497 | imax = i; | |
498 | } | |
19c5a36a | 499 | } |
b31787ca | 500 | if (rmax > rmaxG) { |
501 | rmaxG = rmax; | |
502 | best.CopyResults(fB[imax]); | |
503 | } | |
504 | } | |
505 | ||
506 | mPhi = best.GetMx(); | |
507 | mEta = best.GetMy(); | |
b31787ca | 508 | nk = best.GetK(); |
509 | im = (best.GetInd())[0]; | |
510 | etaC = mEta[im]; | |
511 | ||
19c5a36a | 512 | // |
513 | // Cluster Multiplicity | |
b31787ca | 514 | if (rmaxG > 0. && TMath::Abs(etaC) < 0.4) { |
19c5a36a | 515 | for (Int_t i = 0; i < 1; i++) { |
c87f0b8c | 516 | im = (best.GetInd())[i]; |
19c5a36a | 517 | fH1CEtaR->Fill(mEta[im]); |
518 | } | |
519 | ||
520 | for (Int_t i = 0; i < 1; i++) { | |
c87f0b8c | 521 | im = (best.GetInd())[i]; |
19c5a36a | 522 | for (Int_t j = 0; j < ic; j++) { |
523 | Double_t dphi = DeltaPhi(mPhi[im], phiR[j]); | |
524 | Double_t deta = mEta[im] - etaR[j]; | |
525 | Double_t r = DeltaR(mPhi[im], mEta[im], phiR[j], etaR[j]); | |
526 | fH1DRR->Fill(r); | |
527 | fH2DPhiEtaR->Fill(dphi, TMath::Abs(deta)); | |
528 | if (j == icmax) fH2DPhiEtaLR->Fill(dphi, TMath::Abs(deta)); | |
529 | } | |
530 | } | |
531 | if (nk > 1) { | |
532 | dphic = DeltaPhi(mPhi[0], mPhi[1]); | |
533 | detac = TMath::Abs(mEta[0] - mEta[1]); | |
534 | fH2DPhiEtaCR->Fill(dphic, detac); | |
70f2ce9d | 535 | } |
42dc9410 | 536 | } |
19c5a36a | 537 | PostData(1, fHists); |
70f2ce9d | 538 | } |
539 | ||
540 | Double_t AliAnalysisTaskKMeans::DeltaPhi(Double_t phi1, Double_t phi2) | |
541 | { | |
542 | Double_t dphi = TMath::Abs(phi1 - phi2); | |
543 | if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi; | |
544 | return dphi; | |
545 | } | |
546 | ||
547 | Double_t AliAnalysisTaskKMeans::DeltaR(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2) | |
548 | { | |
549 | Double_t dphi = DeltaPhi(phi1, phi2); | |
550 | Double_t deta = eta1 - eta2; | |
551 | return (TMath::Sqrt(dphi * dphi + deta * deta)); | |
552 | ||
553 | } | |
554 | ||
555 | //________________________________________________________________________ | |
556 | void AliAnalysisTaskKMeans::Terminate(Option_t *) | |
557 | { | |
558 | } | |
559 |