]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - JETAN/AliAnalysisTaskKMeans.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Provide place where T0 group will store resolutions inside ESD data,
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskKMeans.cxx
... / ...
CommitLineData
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 "TFile.h"
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"
41#include "TObjArray.h"
42
43#include "AliAnalysisTask.h"
44#include "AliAnalysisManager.h"
45
46#include "AliVEvent.h"
47#include "AliESDEvent.h"
48#include "AliExternalTrackParam.h"
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()
68 : AliAnalysisTaskSE()
69 ,fK(0)
70 ,fNMin(0)
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)
80 ,fH1PtAS(0)
81 ,fH1PtR(0)
82 ,fH1SPt(0)
83 ,fH1SPtC(0)
84 ,fH1DPhi(0)
85 ,fH1DR(0)
86 ,fH1DRR(0)
87 ,fH2DPhiEta(0)
88 ,fH2DPhiEtaR(0)
89 ,fH2DPhiEtaL(0)
90 ,fH2DPhiEtaLR(0)
91 ,fH2DPhiEtaC(0)
92 ,fH2DPhiEtaCR(0)
93 ,fH1Resp(0)
94 ,fH1RespR(0)
95 ,fH2Sigma(0)
96 ,fH2SigmaR(0)
97 ,fHDensity(0)
98 ,fHCSize(0)
99 ,fHNCluster(0)
100 ,fHPtDensity(0)
101 ,fHDPhi(0)
102 ,fH2EtaPhi(0)
103 ,fH2REtaPhi(0)
104 ,fCuts(0)
105
106{
107 //
108 // Constructor
109 //
110}
111
112//________________________________________________________________________
113AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name)
114 : AliAnalysisTaskSE(name)
115 ,fK(0)
116 ,fNMin(0)
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)
126 ,fH1PtAS(0)
127 ,fH1PtR(0)
128 ,fH1SPt(0)
129 ,fH1SPtC(0)
130 ,fH1DPhi(0)
131 ,fH1DR(0)
132 ,fH1DRR(0)
133 ,fH2DPhiEta(0)
134 ,fH2DPhiEtaR(0)
135 ,fH2DPhiEtaL(0)
136 ,fH2DPhiEtaLR(0)
137 ,fH2DPhiEtaC(0)
138 ,fH2DPhiEtaCR(0)
139 ,fH1Resp(0)
140 ,fH1RespR(0)
141 ,fH2Sigma(0)
142 ,fH2SigmaR(0)
143 ,fHDensity(0)
144 ,fHCSize(0)
145 ,fHNCluster(0)
146 ,fHPtDensity(0)
147 ,fHDPhi(0)
148 ,fH2EtaPhi(0)
149 ,fH2REtaPhi(0)
150 ,fCuts(0)
151
152{
153 //
154 // Constructor
155 //
156 DefineOutput(1, TList::Class());
157}
158
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)
189 ,fHDensity(0)
190 ,fHCSize(0)
191 ,fHNCluster(0)
192 ,fHPtDensity(0)
193 ,fHDPhi(0)
194 ,fH2EtaPhi(0)
195 ,fH2REtaPhi(0)
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
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.);
224 fH1PtAS = new TH1F("fH1PtAS", "pt distribution",50, 0., 10.);
225 fH1PtR = new TH1F("fH1PtR", "pt distribution",50, 0., 10.);
226
227 fH1SPt = new TH1F("fH1SPt", "sum pt distribution",50, 0., 10.);
228 fH1SPtC = new TH1F("fH1SPtC", "sum pt distribution",50, 0., 10.);
229
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.);
239 fH1Resp = new TH1F("fH1Resp", "Responsibility", 50, 0., 1.);
240 fH1RespR = new TH1F("fH1RespR", "Responsibility", 50, 0., 1.);
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.);
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.);
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());
252 fHists->SetOwner();
253
254 fHists->Add(fH1CEta);
255 fHists->Add(fH1CEtaR);
256 fHists->Add(fH1CPhi);
257 fHists->Add(fH2N1N2);
258 fHists->Add(fH1Pt);
259 fHists->Add(fH1PtR);
260 fHists->Add(fH1PtC);
261 fHists->Add(fH1PtC1);
262 fHists->Add(fH1PtC2);
263 fHists->Add(fH1PtAS);
264 fHists->Add(fH1DR);
265 fHists->Add(fH1SPtC);
266 fHists->Add(fH1SPt);
267 fHists->Add(fH1DPhi);
268 fHists->Add(fH2DPhiEta);
269 fHists->Add(fH2DPhiEtaR);
270 fHists->Add(fH2DPhiEtaL);
271 fHists->Add(fH2DPhiEtaLR);
272 fHists->Add(fH2DPhiEtaC);
273 fHists->Add(fH2DPhiEtaCR);
274 fHists->Add(fH1DRR);
275 fHists->Add(fH1RespR);
276 fHists->Add(fH1Resp);
277 fHists->Add(fH2Sigma);
278 fHists->Add(fH2SigmaR);
279 fHists->Add(fHCSize);
280 fHists->Add(fHDensity);
281 fHists->Add(fHNCluster);
282 fHists->Add(fHPtDensity);
283 fHists->Add(fHDPhi);
284 fHists->Add(fH2EtaPhi);
285 //
286 for (Int_t i = 0; i < 10; i++) {
287 fA[i] = new AliKMeansResult(i+1);
288 fB[i] = new AliKMeansResult(i+1);
289 }
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];
299 Double_t eta[500];
300 Double_t etaR[500];
301 Double_t pt [500];
302
303 if (!fInputEvent) {
304 Printf("ERROR: fESD not available");
305 return;
306 }
307
308
309 Int_t ic = 0;
310
311 //
312 // Fill eta-phi positions
313 //
314
315 Double_t ptmax = 0.;
316 Int_t icmax = -1;
317
318 for (Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++) {
319 AliVParticle* track = fInputEvent->GetTrack(iTracks);
320 if ((fCuts->AcceptTrack((AliESDtrack*)track)))
321 {
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
338 if (pt[ic] > ptmax) {
339 ptmax = pt[ic];
340 icmax = ic;
341 }
342
343 fH2EtaPhi->Fill(eta[ic], phi[ic]);
344 ic++;
345 }
346 } //track loop
347
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
355 //
356 // Cluster
357 if (ic < fNMin) {
358 PostData(1, fHists);
359 return;
360 }
361
362 //
363 Double_t rk0[10];
364 AliKMeansResult* res = 0;
365 AliKMeansResult best(10);
366 Float_t rmaxG = -1.;
367
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) {
378 rmax = rk0[i];
379 imax = i;
380 }
381 }
382 if (rmax > rmaxG) {
383 rmaxG = rmax;
384 best.CopyResults(fA[imax]);
385 }
386 }
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]);
396 fHNCluster->Fill(Float_t(nk));
397
398 Double_t dphic, detac;
399
400 if (rmaxG > 0. && TMath::Abs(etaC) < 0.4) {
401 // Analysis
402 //
403 // Cluster Multiplicity
404 Int_t mult[2] = {0, 0};
405 //
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 }
411 //
412 // Random cluster position
413 Int_t ir = Int_t(Float_t(ic) * gRandom->Rndm());
414
415 Double_t crPhi = phi[ir];
416 Double_t crEta = eta[ir];
417 //
418 Double_t sumPt = 0;
419 Double_t sumPtC = 0;
420
421 for (Int_t i = 0; i < 1; i++) {
422 fH1CEta->Fill(mEta[im]);
423 fH1CPhi->Fill(mPhi[im]);
424 for (Int_t j = 0; j < ic; j++) {
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];
428 Double_t rr = DeltaR(crPhi, crEta, phi[j], eta[j]);
429
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 }
438 if (r < 0.3) {
439 fH1PtC1->Fill(pt[j]);
440 fHPtDensity->Fill(rmaxG/TMath::Pi(), pt[j]);
441 }
442 if (rr < 0.3) {
443 fH1PtR->Fill(pt[j]);
444 }
445
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)) {
452 fH1Pt->Fill(pt[j]);
453 }
454
455 if (r > 0.7 && r < 1.) {
456 sumPt += pt[j];
457 }
458
459 if (dphi > (TMath::Pi() - 0.3)) {
460 fH1PtAS->Fill(pt[j], 1.);
461 }
462 }
463 }
464
465 fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
466 fH1SPt ->Fill(sumPt);
467 fH1SPtC->Fill(sumPtC);
468 }
469 //
470 // Randomized phi
471 //
472 rmaxG = -1.;
473 for (Int_t k = 0; k < 20; k++) {
474 Float_t rmax = -1.;
475 Int_t imax = 0;
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);
480 Int_t j = (res->GetInd())[0];
481 rk0[i] = (res->GetTarget())[j];
482 if (rk0[i] > rmax) {
483 rmax = rk0[i];
484 imax = i;
485 }
486 }
487 if (rmax > rmaxG) {
488 rmaxG = rmax;
489 best.CopyResults(fB[imax]);
490 }
491 }
492
493 mPhi = best.GetMx();
494 mEta = best.GetMy();
495 nk = best.GetK();
496 im = (best.GetInd())[0];
497 etaC = mEta[im];
498
499 //
500 // Cluster Multiplicity
501 if (rmaxG > 0. && TMath::Abs(etaC) < 0.4) {
502 for (Int_t i = 0; i < 1; i++) {
503 im = (best.GetInd())[i];
504 fH1CEtaR->Fill(mEta[im]);
505 }
506
507 for (Int_t i = 0; i < 1; i++) {
508 im = (best.GetInd())[i];
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);
522 }
523 }
524 PostData(1, fHists);
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