]> 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
Restoring the functionallity of the PAR files: partial fix
[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
109 // Constructor
110 //
111 for (Int_t i=0; i< 10; i++) {
112 fA[i] = 0;
113 fB[i] = 0;
114 }
115}
116
117//________________________________________________________________________
118AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name)
119 : AliAnalysisTaskSE(name)
120 ,fK(0)
121 ,fNMin(0)
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)
131 ,fH1PtAS(0)
132 ,fH1PtR(0)
133 ,fH1SPt(0)
134 ,fH1SPtC(0)
135 ,fH1DPhi(0)
136 ,fH1DR(0)
137 ,fH1DRR(0)
138 ,fH2DPhiEta(0)
139 ,fH2DPhiEtaR(0)
140 ,fH2DPhiEtaL(0)
141 ,fH2DPhiEtaLR(0)
142 ,fH2DPhiEtaC(0)
143 ,fH2DPhiEtaCR(0)
144 ,fH1Resp(0)
145 ,fH1RespR(0)
146 ,fH2Sigma(0)
147 ,fH2SigmaR(0)
148 ,fHDensity(0)
149 ,fHCSize(0)
150 ,fHNCluster(0)
151 ,fHPtDensity(0)
152 ,fHDPhi(0)
153 ,fH2EtaPhi(0)
154 ,fH2REtaPhi(0)
155 ,fCuts(0)
156
157{
158 //
159 // Constructor
160 //
161 for (Int_t i=0; i< 10; i++) {
162 fA[i] = 0;
163 fB[i] = 0;
164 }
165 DefineOutput(1, TList::Class());
166}
167
168AliAnalysisTaskKMeans::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)
198 ,fHDensity(0)
199 ,fHCSize(0)
200 ,fHNCluster(0)
201 ,fHPtDensity(0)
202 ,fHDPhi(0)
203 ,fH2EtaPhi(0)
204 ,fH2REtaPhi(0)
205 ,fCuts(0)
206{
207 // Dummy copy constructor
208 for (Int_t i=0; i< 10; i++) {
209 fA[i] = 0;
210 fB[i] = 0;
211 }
212}
213
214AliAnalysisTaskKMeans& AliAnalysisTaskKMeans::operator=(const AliAnalysisTaskKMeans& /*trclass*/)
215{
216 // Dummy assignment operator
217 return *this;
218}
219
220
221
222//________________________________________________________________________
223void 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.);
237 fH1PtAS = new TH1F("fH1PtAS", "pt distribution",50, 0., 10.);
238 fH1PtR = new TH1F("fH1PtR", "pt distribution",50, 0., 10.);
239
240 fH1SPt = new TH1F("fH1SPt", "sum pt distribution",50, 0., 10.);
241 fH1SPtC = new TH1F("fH1SPtC", "sum pt distribution",50, 0., 10.);
242
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.);
252 fH1Resp = new TH1F("fH1Resp", "Responsibility", 50, 0., 1.);
253 fH1RespR = new TH1F("fH1RespR", "Responsibility", 50, 0., 1.);
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.);
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.);
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());
265 fHists->SetOwner();
266
267 fHists->Add(fH1CEta);
268 fHists->Add(fH1CEtaR);
269 fHists->Add(fH1CPhi);
270 fHists->Add(fH2N1N2);
271 fHists->Add(fH1Pt);
272 fHists->Add(fH1PtR);
273 fHists->Add(fH1PtC);
274 fHists->Add(fH1PtC1);
275 fHists->Add(fH1PtC2);
276 fHists->Add(fH1PtAS);
277 fHists->Add(fH1DR);
278 fHists->Add(fH1SPtC);
279 fHists->Add(fH1SPt);
280 fHists->Add(fH1DPhi);
281 fHists->Add(fH2DPhiEta);
282 fHists->Add(fH2DPhiEtaR);
283 fHists->Add(fH2DPhiEtaL);
284 fHists->Add(fH2DPhiEtaLR);
285 fHists->Add(fH2DPhiEtaC);
286 fHists->Add(fH2DPhiEtaCR);
287 fHists->Add(fH1DRR);
288 fHists->Add(fH1RespR);
289 fHists->Add(fH1Resp);
290 fHists->Add(fH2Sigma);
291 fHists->Add(fH2SigmaR);
292 fHists->Add(fHCSize);
293 fHists->Add(fHDensity);
294 fHists->Add(fHNCluster);
295 fHists->Add(fHPtDensity);
296 fHists->Add(fHDPhi);
297 fHists->Add(fH2EtaPhi);
298 //
299 for (Int_t i = 0; i < 10; i++) {
300 fA[i] = new AliKMeansResult(i+1);
301 fB[i] = new AliKMeansResult(i+1);
302 }
303}
304
305//________________________________________________________________________
306void AliAnalysisTaskKMeans::UserExec(Option_t *)
307{
308 // Main loop
309 // Called for each event
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};
315
316 if (!fInputEvent) {
317 Printf("ERROR: fESD not available");
318 return;
319 }
320
321
322 Int_t ic = 0;
323
324 //
325 // Fill eta-phi positions
326 //
327
328 Double_t ptmax = 0.;
329 Int_t icmax = -1;
330
331 for (Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++) {
332 AliVParticle* track = fInputEvent->GetTrack(iTracks);
333 if ((fCuts->AcceptTrack((AliESDtrack*)track)))
334 {
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
351 if (pt[ic] > ptmax) {
352 ptmax = pt[ic];
353 icmax = ic;
354 }
355
356 fH2EtaPhi->Fill(eta[ic], phi[ic]);
357 ic++;
358 }
359 } //track loop
360
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
368 //
369 // Cluster
370 if (ic < fNMin) {
371 PostData(1, fHists);
372 return;
373 }
374
375 //
376 Double_t rk0[10];
377 AliKMeansResult* res = 0;
378 AliKMeansResult best(10);
379 Float_t rmaxG = -1.;
380
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) {
391 rmax = rk0[i];
392 imax = i;
393 }
394 }
395 if (rmax > rmaxG) {
396 rmaxG = rmax;
397 best.CopyResults(fA[imax]);
398 }
399 }
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]);
409 fHNCluster->Fill(Float_t(nk));
410
411 Double_t dphic, detac;
412
413 if (rmaxG > 0. && TMath::Abs(etaC) < 0.4) {
414 // Analysis
415 //
416 // Cluster Multiplicity
417 Int_t mult[2] = {0, 0};
418 //
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 }
424 //
425 // Random cluster position
426 Int_t ir = Int_t(Float_t(ic) * gRandom->Rndm());
427
428 Double_t crPhi = phi[ir];
429 Double_t crEta = eta[ir];
430 //
431 Double_t sumPt = 0;
432 Double_t sumPtC = 0;
433
434 for (Int_t i = 0; i < 1; i++) {
435 fH1CEta->Fill(mEta[im]);
436 fH1CPhi->Fill(mPhi[im]);
437 for (Int_t j = 0; j < ic; j++) {
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];
441 Double_t rr = DeltaR(crPhi, crEta, phi[j], eta[j]);
442
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 }
451 if (r < 0.3) {
452 fH1PtC1->Fill(pt[j]);
453 fHPtDensity->Fill(rmaxG/TMath::Pi(), pt[j]);
454 }
455 if (rr < 0.3) {
456 fH1PtR->Fill(pt[j]);
457 }
458
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)) {
465 fH1Pt->Fill(pt[j]);
466 }
467
468 if (r > 0.7 && r < 1.) {
469 sumPt += pt[j];
470 }
471
472 if (dphi > (TMath::Pi() - 0.3)) {
473 fH1PtAS->Fill(pt[j], 1.);
474 }
475 }
476 }
477
478 fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
479 fH1SPt ->Fill(sumPt);
480 fH1SPtC->Fill(sumPtC);
481 }
482 //
483 // Randomized phi
484 //
485 rmaxG = -1.;
486 for (Int_t k = 0; k < 20; k++) {
487 Float_t rmax = -1.;
488 Int_t imax = 0;
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);
493 Int_t j = (res->GetInd())[0];
494 rk0[i] = (res->GetTarget())[j];
495 if (rk0[i] > rmax) {
496 rmax = rk0[i];
497 imax = i;
498 }
499 }
500 if (rmax > rmaxG) {
501 rmaxG = rmax;
502 best.CopyResults(fB[imax]);
503 }
504 }
505
506 mPhi = best.GetMx();
507 mEta = best.GetMy();
508 nk = best.GetK();
509 im = (best.GetInd())[0];
510 etaC = mEta[im];
511
512 //
513 // Cluster Multiplicity
514 if (rmaxG > 0. && TMath::Abs(etaC) < 0.4) {
515 for (Int_t i = 0; i < 1; i++) {
516 im = (best.GetInd())[i];
517 fH1CEtaR->Fill(mEta[im]);
518 }
519
520 for (Int_t i = 0; i < 1; i++) {
521 im = (best.GetInd())[i];
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);
535 }
536 }
537 PostData(1, fHists);
538}
539
540Double_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
547Double_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//________________________________________________________________________
556void AliAnalysisTaskKMeans::Terminate(Option_t *)
557{
558}
559