1 /**************************************************************************
2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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.
22 // Author: Andreas Morsch (CERN)
23 // andreas.morsch@cern.ch
24 //---------------------------------------------------------------------------------------
36 #include "TParticle.h"
37 #include "TParticlePDG.h"
41 #include "TObjArray.h"
43 #include "AliAnalysisTask.h"
44 #include "AliAnalysisManager.h"
46 #include "AliVEvent.h"
47 #include "AliESDEvent.h"
48 #include "AliExternalTrackParam.h"
50 #include "AliESDVertex.h"
51 #include "AliESDInputHandler.h"
52 #include "AliESDtrackCuts.h"
53 #include "AliMultiplicity.h"
55 #include "AliMCParticle.h"
56 #include "AliMCEvent.h"
57 #include "AliAnalysisTaskKMeans.h"
58 #include "AliTrackReference.h"
60 #include "AliHeader.h"
61 #include "AliKMeansClustering.h"
65 ClassImp(AliAnalysisTaskKMeans)
67 AliAnalysisTaskKMeans::AliAnalysisTaskKMeans()
111 for (Int_t i=0; i< 10; i++) {
117 //________________________________________________________________________
118 AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name)
119 : AliAnalysisTaskSE(name)
161 for (Int_t i=0; i< 10; i++) {
165 DefineOutput(1, TList::Class());
168 AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const AliAnalysisTaskKMeans &res)
169 : AliAnalysisTaskSE(res)
207 // Dummy copy constructor
208 for (Int_t i=0; i< 10; i++) {
214 AliAnalysisTaskKMeans& AliAnalysisTaskKMeans::operator=(const AliAnalysisTaskKMeans& /*trclass*/)
216 // Dummy assignment operator
222 //________________________________________________________________________
223 void AliAnalysisTaskKMeans::UserCreateOutputObjects()
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.);
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.);
240 fH1SPt = new TH1F("fH1SPt", "sum pt distribution",50, 0., 10.);
241 fH1SPtC = new TH1F("fH1SPtC", "sum pt distribution",50, 0., 10.);
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.);
258 fHDensity = new TH1F("fHDensity", "density distribution", 100, 0., 20.);
259 fHCSize = new TH1F("fHCSize", "cluster size", 20, -0.5, 19.5);
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());
267 fHists->Add(fH1CEta);
268 fHists->Add(fH1CEtaR);
269 fHists->Add(fH1CPhi);
270 fHists->Add(fH2N1N2);
274 fHists->Add(fH1PtC1);
275 fHists->Add(fH1PtC2);
276 fHists->Add(fH1PtAS);
278 fHists->Add(fH1SPtC);
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);
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);
297 fHists->Add(fH2EtaPhi);
299 for (Int_t i = 0; i < 10; i++) {
300 fA[i] = new AliKMeansResult(i+1);
301 fB[i] = new AliKMeansResult(i+1);
305 //________________________________________________________________________
306 void AliAnalysisTaskKMeans::UserExec(Option_t *)
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};
317 Printf("ERROR: fESD not available");
325 // Fill eta-phi positions
331 for (Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++) {
332 AliVParticle* track = fInputEvent->GetTrack(iTracks);
333 if ((fCuts->AcceptTrack((AliESDtrack*)track)))
335 const AliExternalTrackParam * tpcT = ((AliESDtrack*) track)->GetTPCInnerParam();
337 if (TMath::Abs(tpcT->Eta()) > 0.9) continue;
339 phi [ic] = tpcT->Phi();
340 eta [ic] = tpcT->Eta();
341 pt [ic] = tpcT->Pt();
344 fH2REtaPhi->GetRandom2(etaR[ic], phiR[ic]);
346 phiR[ic] = 2. * TMath::Pi() * gRandom->Rndm();
347 etaR[ic] = 1.8 * gRandom->Rndm() - 0.9;
351 if (pt[ic] > ptmax) {
356 fH2EtaPhi->Fill(eta[ic], phi[ic]);
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]);
377 AliKMeansResult* res = 0;
378 AliKMeansResult best(10);
381 for (Int_t k = 0; k < 20; k++) {
384 for (Int_t i = 0; i < fK; 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];
397 best.CopyResults(fA[imax]);
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];
407 fHDensity->Fill(rmaxG / TMath::Pi());
408 fHCSize->Fill(2.28 * rmaxG * sigma2[im]);
409 fHNCluster->Fill(Float_t(nk));
411 Double_t dphic, detac;
413 if (rmaxG > 0. && TMath::Abs(etaC) < 0.4) {
416 // Cluster Multiplicity
417 Int_t mult[2] = {0, 0};
420 dphic = DeltaPhi(mPhi[0], mPhi[1]);
421 detac = TMath::Abs(mEta[0] - mEta[1]);
422 fH2DPhiEtaC->Fill(dphic, detac);
425 // Random cluster position
426 Int_t ir = Int_t(Float_t(ic) * gRandom->Rndm());
428 Double_t crPhi = phi[ir];
429 Double_t crEta = eta[ir];
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]);
445 fH2DPhiEta->Fill(dphi, TMath::Abs(deta));
446 if (j == icmax) fH2DPhiEtaL->Fill(dphi, TMath::Abs(deta));
449 fH1PtC2->Fill(pt[j]);
452 fH1PtC1->Fill(pt[j]);
453 fHPtDensity->Fill(rmaxG/TMath::Pi(), pt[j]);
464 if (r > 0.7 && dphi < (TMath::Pi() - 0.3)) {
468 if (r > 0.7 && r < 1.) {
472 if (dphi > (TMath::Pi() - 0.3)) {
473 fH1PtAS->Fill(pt[j], 1.);
478 fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
479 fH1SPt ->Fill(sumPt);
480 fH1SPtC->Fill(sumPtC);
486 for (Int_t k = 0; k < 20; k++) {
489 for (Int_t i = 0; i < fK; 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];
502 best.CopyResults(fB[imax]);
509 im = (best.GetInd())[0];
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]);
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]);
527 fH2DPhiEtaR->Fill(dphi, TMath::Abs(deta));
528 if (j == icmax) fH2DPhiEtaLR->Fill(dphi, TMath::Abs(deta));
532 dphic = DeltaPhi(mPhi[0], mPhi[1]);
533 detac = TMath::Abs(mEta[0] - mEta[1]);
534 fH2DPhiEtaCR->Fill(dphic, detac);
540 Double_t AliAnalysisTaskKMeans::DeltaPhi(Double_t phi1, Double_t phi2)
542 Double_t dphi = TMath::Abs(phi1 - phi2);
543 if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
547 Double_t AliAnalysisTaskKMeans::DeltaR(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
549 Double_t dphi = DeltaPhi(phi1, phi2);
550 Double_t deta = eta1 - eta2;
551 return (TMath::Sqrt(dphi * dphi + deta * deta));
555 //________________________________________________________________________
556 void AliAnalysisTaskKMeans::Terminate(Option_t *)