- Chosing optimal number of clusters
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskKMeans.cxx
CommitLineData
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"
19c5a36a 40#include "TObjArray.h"
70f2ce9d 41
42#include "AliAnalysisTask.h"
43#include "AliAnalysisManager.h"
44
42dc9410 45#include "AliVEvent.h"
70f2ce9d 46#include "AliESDEvent.h"
47#include "AliStack.h"
48#include "AliESDVertex.h"
49#include "AliESDInputHandler.h"
50#include "AliESDtrackCuts.h"
51#include "AliMultiplicity.h"
52
53#include "AliMCParticle.h"
54#include "AliMCEvent.h"
55#include "AliAnalysisTaskKMeans.h"
56#include "AliTrackReference.h"
57#include "AliStack.h"
58#include "AliHeader.h"
59#include "AliKMeansClustering.h"
60
61
62
63ClassImp(AliAnalysisTaskKMeans)
64
65AliAnalysisTaskKMeans::AliAnalysisTaskKMeans()
320bb308 66 : AliAnalysisTaskSE()
67 ,fK(0)
42dc9410 68 ,fNMin(0)
70f2ce9d 69 ,fHists(0)
70 ,fH1CEta(0)
71 ,fH1CPhi(0)
72 ,fH1CEtaR(0)
73 ,fH2N1N2(0)
74 ,fH1Pt(0)
75 ,fH1PtC(0)
76 ,fH1PtC1(0)
77 ,fH1PtC2(0)
42dc9410 78 ,fH1PtAS(0)
79 ,fH1PtR(0)
a51f90df 80 ,fH1SPt(0)
81 ,fH1SPtC(0)
70f2ce9d 82 ,fH1DPhi(0)
83 ,fH1DR(0)
84 ,fH1DRR(0)
85 ,fH2DPhiEta(0)
86 ,fH2DPhiEtaR(0)
87 ,fH2DPhiEtaL(0)
5a9b5891 88 ,fH2DPhiEtaLR(0)
a51f90df 89 ,fH2DPhiEtaC(0)
90 ,fH2DPhiEtaCR(0)
320bb308 91 ,fH1Resp(0)
92 ,fH1RespR(0)
42dc9410 93 ,fH2Sigma(0)
94 ,fH2SigmaR(0)
19c5a36a 95 ,fHDensity(0)
96 ,fHCSize(0)
97 ,fHNCluster(0)
98 ,fHPtDensity(0)
70f2ce9d 99 ,fCuts(0)
100{
101 //
102 // Constructor
103 //
104}
105
106//________________________________________________________________________
107AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name)
108 : AliAnalysisTaskSE(name)
320bb308 109 ,fK(0)
42dc9410 110 ,fNMin(0)
70f2ce9d 111 ,fHists(0)
112 ,fH1CEta(0)
113 ,fH1CPhi(0)
114 ,fH1CEtaR(0)
115 ,fH2N1N2(0)
116 ,fH1Pt(0)
117 ,fH1PtC(0)
118 ,fH1PtC1(0)
119 ,fH1PtC2(0)
42dc9410 120 ,fH1PtAS(0)
121 ,fH1PtR(0)
a51f90df 122 ,fH1SPt(0)
123 ,fH1SPtC(0)
70f2ce9d 124 ,fH1DPhi(0)
125 ,fH1DR(0)
126 ,fH1DRR(0)
127 ,fH2DPhiEta(0)
128 ,fH2DPhiEtaR(0)
129 ,fH2DPhiEtaL(0)
5a9b5891 130 ,fH2DPhiEtaLR(0)
a51f90df 131 ,fH2DPhiEtaC(0)
132 ,fH2DPhiEtaCR(0)
320bb308 133 ,fH1Resp(0)
134 ,fH1RespR(0)
42dc9410 135 ,fH2Sigma(0)
136 ,fH2SigmaR(0)
19c5a36a 137 ,fHDensity(0)
138 ,fHCSize(0)
139 ,fHNCluster(0)
140 ,fHPtDensity(0)
70f2ce9d 141 ,fCuts(0)
142{
143 //
144 // Constructor
145 //
146 DefineOutput(1, TList::Class());
147}
148
149
150//________________________________________________________________________
151void AliAnalysisTaskKMeans::UserCreateOutputObjects()
152{
153 // Create histograms
154 // Called once
155 fHists = new TList();
156 fH1CEta = new TH1F("fH1CEta", "eta distribution of clusters", 90, -1.5, 1.5);
157 fH1CEtaR = new TH1F("fH1CEtaR", "eta distribution of clusters", 90, -1.5, 1.5);
158 fH1CPhi = new TH1F("fH1CPhi", "phi distribution of clusters", 157, 0.0, 2. * TMath::Pi());
159 fH2N1N2 = new TH2F("fH2N1N2", "multiplicity distribution", 50, 0., 50., 50, 0., 50.);
160
161 fH1Pt = new TH1F("fH1Pt", "pt distribution",50, 0., 10.);
162 fH1PtC = new TH1F("fH1PtC", "pt distribution",50, 0., 10.);
163 fH1PtC1 = new TH1F("fH1PtC1", "pt distribution",50, 0., 10.);
164 fH1PtC2 = new TH1F("fH1PtC2", "pt distribution",50, 0., 10.);
42dc9410 165 fH1PtAS = new TH1F("fH1PtAS", "pt distribution",50, 0., 10.);
166 fH1PtR = new TH1F("fH1PtR", "pt distribution",50, 0., 10.);
70f2ce9d 167
a51f90df 168 fH1SPt = new TH1F("fH1SPt", "sum pt distribution",50, 0., 10.);
169 fH1SPtC = new TH1F("fH1SPtC", "sum pt distribution",50, 0., 10.);
170
5a9b5891 171 fH1DR = new TH1F("fH1DR", "dR distribution", 50, 0., 5.);
172 fH1DPhi = new TH1F("fH1DPhi", "dPhi distribution", 31, 0., TMath::Pi());
173 fH2DPhiEta = new TH2F("fH2DPhiEta","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
174 fH2DPhiEtaR = new TH2F("fH2DPhiEtaR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
175 fH2DPhiEtaL = new TH2F("fH2DPhiEtaL","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
176 fH2DPhiEtaLR = new TH2F("fH2DPhiEtaLR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
177 fH2DPhiEtaC = new TH2F("fH2DPhiEtaC","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
178 fH2DPhiEtaCR = new TH2F("fH2DPhiEtaCR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
179 fH1DRR = new TH1F("fH1DRR", "dR distribution", 50, 0., 5.);
320bb308 180 fH1Resp = new TH1F("fH1Resp", "Responsibility", 50, 0., 1.);
181 fH1RespR = new TH1F("fH1RespR", "Responsibility", 50, 0., 1.);
42dc9410 182 fH2Sigma = new TH2F("fH2Sigma", "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.);
183 fH2SigmaR = new TH2F("fH2SigmaR", "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.);
19c5a36a 184
185
186 fHDensity = new TH1F("fHDensity", "density distribution", 100, 0., 20.);
187 fHCSize = new TH1F("fHCSize", "cluster size", 20, -0.5, 19.5);
188
189 fHNCluster = new TH1F("fHNCluster", "Number of clusters", 11, -0.5, 10.5);
190 fHPtDensity = new TH2F("fHPtDensity", "Pt vs density", 100, 0., 20., 50, 0., 10.);
191
42dc9410 192
70f2ce9d 193 fHists->SetOwner();
194
195 fHists->Add(fH1CEta);
196 fHists->Add(fH1CEtaR);
197 fHists->Add(fH1CPhi);
198 fHists->Add(fH2N1N2);
199 fHists->Add(fH1Pt);
42dc9410 200 fHists->Add(fH1PtR);
70f2ce9d 201 fHists->Add(fH1PtC);
202 fHists->Add(fH1PtC1);
203 fHists->Add(fH1PtC2);
42dc9410 204 fHists->Add(fH1PtAS);
70f2ce9d 205 fHists->Add(fH1DR);
a51f90df 206 fHists->Add(fH1SPtC);
207 fHists->Add(fH1SPt);
70f2ce9d 208 fHists->Add(fH1DPhi);
209 fHists->Add(fH2DPhiEta);
210 fHists->Add(fH2DPhiEtaR);
211 fHists->Add(fH2DPhiEtaL);
5a9b5891 212 fHists->Add(fH2DPhiEtaLR);
a51f90df 213 fHists->Add(fH2DPhiEtaC);
214 fHists->Add(fH2DPhiEtaCR);
70f2ce9d 215 fHists->Add(fH1DRR);
320bb308 216 fHists->Add(fH1RespR);
217 fHists->Add(fH1Resp);
42dc9410 218 fHists->Add(fH2Sigma);
19c5a36a 219 fHists->Add(fH2SigmaR);
220 fHists->Add(fHCSize);
221 fHists->Add(fHDensity);
222 fHists->Add(fHNCluster);
223 fHists->Add(fHPtDensity);
70f2ce9d 224 //
19c5a36a 225 for (Int_t i = 0; i < 10; i++) {
226 fA[i] = new AliKMeansResult(i+1);
227 fB[i] = new AliKMeansResult(i+1);
228 }
70f2ce9d 229}
230
231//________________________________________________________________________
232void AliAnalysisTaskKMeans::UserExec(Option_t *)
233{
234 // Main loop
235 // Called for each event
236 Double_t phi [500];
237 Double_t phiR[500];
19c5a36a 238 Double_t eta[500];
70f2ce9d 239 Double_t etaR[500];
19c5a36a 240 Double_t pt [500];
70f2ce9d 241
242 if (!fInputEvent) {
243 Printf("ERROR: fESD not available");
244 return;
245 }
246
70f2ce9d 247
248 Int_t ic = 0;
19c5a36a 249
70f2ce9d 250 //
251 // Fill eta-phi positions
252 //
19c5a36a 253
70f2ce9d 254 Double_t ptmax = 0.;
255 Int_t icmax = -1;
256
42dc9410 257 for (Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++) {
258 AliVParticle* track = fInputEvent->GetTrack(iTracks);
259 if ((fCuts->AcceptTrack((AliESDtrack*)track)))
70f2ce9d 260 {
19c5a36a 261 phi [ic] = track->Phi();
70f2ce9d 262 phiR[ic] = 2. * TMath::Pi() * gRandom->Rndm();
19c5a36a 263 eta [ic] = track->Eta();
264 etaR[ic] = 1.8 * gRandom->Rndm() - 0.9;
265 pt [ic] = track->Pt();
70f2ce9d 266 if (pt[ic] > ptmax) {
267 ptmax = pt[ic];
268 icmax = ic;
269 }
270 ic++;
271 }
272 } //track loop
273
274 //
275 // Cluster
42dc9410 276 if (ic < fNMin) {
70f2ce9d 277 PostData(1, fHists);
278 return;
279 }
280 //
19c5a36a 281 Double_t* rk0 = new Double_t[10];
282
283 AliKMeansResult* res = 0;
284
285 for (Int_t i = 0; i < fK; i++) {
286 res = fA[i];
287 AliKMeansClustering::SoftKMeans2(i+1, ic, phi, eta, res->GetMx(),res->GetMy(), res->GetSigma2(), res->GetRk());
288 res->Sort(ic, phi, eta);
289 Int_t j = (res->GetInd())[0];
290 rk0[i] = (res->GetTarget())[j];
291 }
292
293 Double_t* mPhi = 0;
294 Double_t* mEta = 0;
295 Double_t* sigma2 = 0;
296 Float_t rmax = -1.;
297 Int_t imax = 0;
298
299 for (Int_t i = 0; i < fK; i++) {
300 if (rk0[i] > rmax) {
301 rmax = rk0[i];
302 imax = i;
303 }
304 }
305
306
307 res = fA[imax];
308 mPhi = res->GetMx();
309 mEta = res->GetMy();
310 sigma2 = res->GetSigma2();
311 Int_t nk = imax + 1;
312 Int_t im = (res->GetInd())[0];
313 fHDensity->Fill(rmax / TMath::Pi());
314 fHCSize->Fill(rmax * sigma2[im]);
315 fHNCluster->Fill(Float_t(nk));
316
42dc9410 317 Double_t dphic, detac;
42dc9410 318
19c5a36a 319 if (rmax > 0.) {
42dc9410 320 // Analysis
321 //
322 // Cluster Multiplicity
323 Int_t mult[2] = {0, 0};
324 //
19c5a36a 325 if (nk > 1) {
326 dphic = DeltaPhi(mPhi[0], mPhi[1]);
327 detac = TMath::Abs(mEta[0] - mEta[1]);
328 fH2DPhiEtaC->Fill(dphic, detac);
329 }
42dc9410 330 //
331 // Random cluster position
332 Int_t ir = Int_t(Float_t(ic) * gRandom->Rndm());
a51f90df 333
42dc9410 334 Double_t crPhi = phi[ir];
335 Double_t crEta = eta[ir];
336 //
337 Double_t sumPt = 0;
338 Double_t sumPtC = 0;
19c5a36a 339
42dc9410 340 for (Int_t i = 0; i < 1; i++) {
19c5a36a 341 fH1CEta->Fill(mEta[im]);
342 fH1CPhi->Fill(mPhi[im]);
70f2ce9d 343 for (Int_t j = 0; j < ic; j++) {
19c5a36a 344 Double_t r = DeltaR(mPhi[im], mEta[im], phi[j], eta[j]);
345 Double_t dphi = DeltaPhi(mPhi[im], phi[j]);
346 Double_t deta = mEta[im] - eta[j];
42dc9410 347 Double_t rr = DeltaR(crPhi, crEta, phi[j], eta[j]);
19c5a36a 348
42dc9410 349 fH1DR->Fill(r);
350 fH1DPhi->Fill(dphi);
351 fH2DPhiEta->Fill(dphi, TMath::Abs(deta));
352 if (j == icmax) fH2DPhiEtaL->Fill(dphi, TMath::Abs(deta));
353
354 if (r < 0.2) {
355 fH1PtC2->Fill(pt[j]);
356 }
19c5a36a 357 if (r < 0.3) {
358 fH1PtC1->Fill(pt[j]);
359 fHPtDensity->Fill(rmax/TMath::Pi(), pt[j]);
360 }
361 if (rr < 0.3) {
42dc9410 362 fH1PtR->Fill(pt[j]);
19c5a36a 363 }
42dc9410 364
19c5a36a 365 if (r < 0.4) {
366 sumPtC += pt[j];
367 mult[i]++;
368 fH1PtC->Fill(pt[j]);
369 }
370 if (r > 0.7 && dphi < (TMath::Pi() - 0.3)) {
42dc9410 371 fH1Pt->Fill(pt[j]);
372 }
373
374 if (r > 0.7 && r < 1.) {
375 sumPt += pt[j];
376 }
19c5a36a 377
42dc9410 378 if (dphi > (TMath::Pi() - 0.3)) {
379 fH1PtAS->Fill(pt[j]);
380 }
70f2ce9d 381 }
42dc9410 382 }
70f2ce9d 383
42dc9410 384 fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
385 fH1SPt ->Fill(sumPt);
386 fH1SPtC->Fill(sumPtC);
19c5a36a 387 }
42dc9410 388 //
19c5a36a 389 // Randomized phi
42dc9410 390 //
19c5a36a 391 for (Int_t i = 0; i < fK; i++) {
392 res = fB[i];
393 AliKMeansClustering::SoftKMeans2(i+1, ic, phiR, etaR, res->GetMx(),res->GetMy(), res->GetSigma2(), res->GetRk());
394 res->Sort(ic, phiR, etaR);
395 //res->Sort();
396 Int_t j = (res->GetInd())[0];
397 rk0[i] = (res->GetTarget())[j];
42dc9410 398 }
399
19c5a36a 400 rmax = -1.;
401 imax = 0;
402 for (Int_t i = 0; i < fK; i++) {
403 if (rk0[i] > rmax) {
404 rmax = rk0[i];
405 imax = i;
406 }
407 }
408 res = fB[imax];
409 mPhi = res->GetMx();
410 mEta = res->GetMy();
411 sigma2 = res->GetSigma2();
412 nk = imax + 1;
413
414 // fH1RespR->Fill(rk[ind[0]]/(rk[0]+rk[1]));
415 //
416 // Cluster Multiplicity
417 if (rmax > 0.) {
418 for (Int_t i = 0; i < 1; i++) {
419 Int_t im = (res->GetInd())[i];
420 fH1CEtaR->Fill(mEta[im]);
421 }
422
423 for (Int_t i = 0; i < 1; i++) {
424 Int_t im = (res->GetInd())[i];
425 for (Int_t j = 0; j < ic; j++) {
426 Double_t dphi = DeltaPhi(mPhi[im], phiR[j]);
427 Double_t deta = mEta[im] - etaR[j];
428 Double_t r = DeltaR(mPhi[im], mEta[im], phiR[j], etaR[j]);
429 fH1DRR->Fill(r);
430 fH2DPhiEtaR->Fill(dphi, TMath::Abs(deta));
431 if (j == icmax) fH2DPhiEtaLR->Fill(dphi, TMath::Abs(deta));
432 }
433 }
434 if (nk > 1) {
435 dphic = DeltaPhi(mPhi[0], mPhi[1]);
436 detac = TMath::Abs(mEta[0] - mEta[1]);
437 fH2DPhiEtaCR->Fill(dphic, detac);
70f2ce9d 438 }
42dc9410 439 }
19c5a36a 440 PostData(1, fHists);
70f2ce9d 441}
442
443Double_t AliAnalysisTaskKMeans::DeltaPhi(Double_t phi1, Double_t phi2)
444{
445 Double_t dphi = TMath::Abs(phi1 - phi2);
446 if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
447 return dphi;
448}
449
450Double_t AliAnalysisTaskKMeans::DeltaR(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
451{
452 Double_t dphi = DeltaPhi(phi1, phi2);
453 Double_t deta = eta1 - eta2;
454 return (TMath::Sqrt(dphi * dphi + deta * deta));
455
456}
457
458//________________________________________________________________________
459void AliAnalysisTaskKMeans::Terminate(Option_t *)
460{
461}
462