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