New task for k-means clustering
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 21 Jan 2010 11:24:00 +0000 (11:24 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 21 Jan 2010 11:24:00 +0000 (11:24 +0000)
JETAN/AliAnalysisTaskKMeans.cxx [new file with mode: 0644]
JETAN/AliAnalysisTaskKMeans.h [new file with mode: 0644]
JETAN/AliKMeansClustering.cxx [new file with mode: 0755]
JETAN/AliKMeansClustering.h [new file with mode: 0755]
JETAN/JETANLinkDef.h
JETAN/libJETAN.pkg
PWG4/macros/AddTaskKMeans.C [new file with mode: 0644]

diff --git a/JETAN/AliAnalysisTaskKMeans.cxx b/JETAN/AliAnalysisTaskKMeans.cxx
new file mode 100644 (file)
index 0000000..87b9897
--- /dev/null
@@ -0,0 +1,313 @@
+/**************************************************************************
+ * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+//-------------------------------------------------------------------------
+//     Analysis Task that uses the Soft K-Means Algorithm to find clusters in
+//     the eta-phi space of Minimum Bias. No pt information is used for the clustering.
+//     
+//
+//     Author: Andreas Morsch (CERN)
+//     andreas.morsch@cern.ch
+//-------------------------------------------------------------------------
+
+
+
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TProfile.h"
+
+#include "TList.h"
+#include "TParticle.h"
+#include "TParticlePDG.h"
+#include "TProfile.h"
+#include "TMath.h"
+#include "TRandom.h"
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+
+#include "AliESDEvent.h"
+#include "AliStack.h"
+#include "AliESDVertex.h"
+#include "AliESDInputHandler.h"
+#include "AliESDtrackCuts.h"
+#include "AliMultiplicity.h"
+
+#include "AliMCParticle.h"
+#include "AliMCEvent.h"
+#include "AliAnalysisTaskKMeans.h"
+#include "AliTrackReference.h"
+#include "AliStack.h"
+#include "AliHeader.h"
+#include "AliKMeansClustering.h"
+
+
+
+ClassImp(AliAnalysisTaskKMeans)
+
+AliAnalysisTaskKMeans::AliAnalysisTaskKMeans() 
+    : AliAnalysisTaskSE() 
+    ,fHists(0)
+    ,fH1CEta(0)
+    ,fH1CPhi(0)
+    ,fH1CEtaR(0)
+    ,fH2N1N2(0)
+    ,fH1Pt(0)
+    ,fH1PtC(0)
+    ,fH1PtC1(0)
+    ,fH1PtC2(0)
+    ,fH1DPhi(0)
+    ,fH1DR(0)
+    ,fH1DRR(0)
+    ,fH2DPhiEta(0)
+    ,fH2DPhiEtaR(0)
+    ,fH2DPhiEtaL(0)
+    ,fCuts(0)
+{
+  //
+  // Constructor
+  //
+}
+
+//________________________________________________________________________
+AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name) 
+    : AliAnalysisTaskSE(name) 
+      ,fHists(0)
+      ,fH1CEta(0)
+      ,fH1CPhi(0)
+      ,fH1CEtaR(0)
+      ,fH2N1N2(0)
+      ,fH1Pt(0)
+      ,fH1PtC(0)
+      ,fH1PtC1(0)
+      ,fH1PtC2(0)
+      ,fH1DPhi(0)
+      ,fH1DR(0)
+      ,fH1DRR(0)
+      ,fH2DPhiEta(0)
+      ,fH2DPhiEtaR(0)
+      ,fH2DPhiEtaL(0)
+      ,fCuts(0)
+{
+  //
+  // Constructor
+  //
+  DefineOutput(1,  TList::Class());
+}
+
+
+//________________________________________________________________________
+void AliAnalysisTaskKMeans::UserCreateOutputObjects()
+{
+  // Create histograms
+  // Called once
+    fHists   = new TList();
+    fH1CEta  = new TH1F("fH1CEta",   "eta distribution of clusters",        90, -1.5, 1.5);
+    fH1CEtaR = new TH1F("fH1CEtaR",  "eta distribution of clusters",        90, -1.5, 1.5);
+    fH1CPhi  = new TH1F("fH1CPhi",   "phi distribution of clusters",       157,  0.0, 2. * TMath::Pi());
+    fH2N1N2  = new TH2F("fH2N1N2",   "multiplicity distribution",          50, 0., 50., 50, 0., 50.);
+    
+    fH1Pt    = new TH1F("fH1Pt",     "pt distribution",50, 0., 10.);
+    fH1PtC   = new TH1F("fH1PtC",    "pt distribution",50, 0., 10.);
+    fH1PtC1  = new TH1F("fH1PtC1",   "pt distribution",50, 0., 10.);
+    fH1PtC2  = new TH1F("fH1PtC2",   "pt distribution",50, 0., 10.);
+
+    fH1DR       = new TH1F("fH1DR",     "dR distribution", 50, 0., 5.);
+    fH1DPhi     = new TH1F("fH1DPhi",   "dPhi distribution", 31, 0., TMath::Pi());
+    fH2DPhiEta  = new TH2F("fH2DPhiEta","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
+    fH2DPhiEtaR = new TH2F("fH2DPhiEtaR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
+    fH2DPhiEtaL = new TH2F("fH2DPhiEtaL","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
+    fH1DRR      = new TH1F("fH1DRR",    "dR distribution", 50, 0., 5.);
+    
+    fHists->SetOwner();
+
+    fHists->Add(fH1CEta);
+    fHists->Add(fH1CEtaR);
+    fHists->Add(fH1CPhi);
+    fHists->Add(fH2N1N2);
+    fHists->Add(fH1Pt);
+    fHists->Add(fH1PtC);
+    fHists->Add(fH1PtC1);
+    fHists->Add(fH1PtC2);
+    fHists->Add(fH1DR);
+    fHists->Add(fH1DPhi);
+    fHists->Add(fH2DPhiEta);
+    fHists->Add(fH2DPhiEtaR);
+    fHists->Add(fH2DPhiEtaL);
+    fHists->Add(fH1DRR);
+    
+    //
+    AliKMeansClustering::SetBeta(20.);
+    
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskKMeans::UserExec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
+    Double_t   phi [500];
+    Double_t   phiR[500];
+    Double_t   eta[500];
+    Double_t   etaR[500];
+    Double_t   pt [500];
+    Double_t   mPhi[20];
+    Double_t   mEta[20];
+    Double_t     rk[20];
+    Int_t       ind[20];
+
+  if (!fInputEvent) {
+    Printf("ERROR: fESD not available");
+    return;
+  }
+  
+  AliESDEvent* esdE = (AliESDEvent*) fInputEvent;
+
+  Int_t ic = 0;
+  //
+  // Fill eta-phi positions
+  //
+  /*
+  const AliMultiplicity *spdMult = esdE->GetMultiplicity();
+  for (Int_t i = 0; i < spdMult->GetNumberOfTracklets(); ++i) {
+         phi[ic] = spdMult->GetPhi(i);
+         eta[ic] = spdMult->GetEta(i);
+         ic++;
+  }
+  */
+  Double_t ptmax = 0.;
+  Int_t    icmax = -1;
+  
+  for (Int_t iTracks = 0; iTracks < esdE->GetNumberOfTracks(); iTracks++) {
+      AliESDtrack* track = esdE->GetTrack(iTracks);
+      if ((fCuts->AcceptTrack(track))) 
+      {
+         phi[ic]  = track->Phi();
+         phiR[ic] = 2. * TMath::Pi() * gRandom->Rndm();
+         eta[ic]  = track->Eta();
+         etaR[ic] = 2. * gRandom->Rndm() - 1.;
+         pt[ic]   = track->Pt();
+         if (pt[ic] > ptmax) {
+             ptmax = pt[ic];
+             icmax = ic;
+         }
+         ic++;
+      }
+  } //track loop 
+
+  //
+  // Cluster
+  if (ic < 10) {
+      PostData(1, fHists);
+      return;
+  }
+  //
+  AliKMeansClustering::SoftKMeans(2, ic, phi, eta, mPhi, mEta, rk);
+  //
+  // Sort
+  TMath::Sort(2, rk, ind);
+  //
+  // Analyse
+  //
+  // Cluster Multiplicity
+  Int_t mult[2] = {0, 0};
+  
+  for (Int_t i = 0; i < 1; i++) {
+      fH1CEta->Fill(mEta[ind[i]]);
+      fH1CPhi->Fill(mPhi[ind[i]]);      
+      for (Int_t j = 0; j < ic; j++) {
+         Double_t r    = DeltaR(mPhi[ind[i]], mEta[ind[i]], phi[j], eta[j]);
+         Double_t dphi = DeltaPhi(mPhi[ind[i]], phi[j]);
+         Double_t deta = mEta[ind[i]] - eta[j];
+         
+         fH1DR->Fill(r);
+         fH1DPhi->Fill(dphi);
+         fH2DPhiEta->Fill(dphi, TMath::Abs(deta));
+         if (j == icmax) fH2DPhiEtaL->Fill(dphi, TMath::Abs(deta));
+
+         if (r < 0.2) {
+             fH1PtC2->Fill(pt[j]);
+         }
+         if (r < 0.3) 
+         {
+             fH1PtC1->Fill(pt[j]);
+         }
+         if (r < 0.4) 
+         {
+             mult[i]++;
+             fH1PtC->Fill(pt[j]);
+         } 
+         if (r > 0.7) {
+             fH1Pt->Fill(pt[j]);
+         }
+      }
+  }
+
+  fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
+
+
+
+  // Randomized phi
+  //
+  AliKMeansClustering::SoftKMeans(2, ic, phiR, etaR, mPhi, mEta, rk);
+  //
+  // Sort
+  TMath::Sort(2, rk, ind);
+  //
+  // Analyse
+  //
+  // Cluster Multiplicity
+  for (Int_t i = 0; i < 1; i++) {
+      fH1CEtaR->Fill(mEta[ind[i]]);
+  }
+  for (Int_t i = 0; i < 2; i++) {
+      for (Int_t j = 0; j < ic; j++) {
+         Double_t dphi = DeltaPhi(mPhi[ind[i]], phi[j]);
+         Double_t deta = mEta[ind[i]] - eta[j];
+         Double_t r = DeltaR(mPhi[ind[i]], mEta[ind[i]], phi[j], eta[j]);
+         fH1DRR->Fill(r);
+         fH2DPhiEtaR->Fill(dphi, TMath::Abs(deta));
+      }
+  }
+
+  
+  //
+  // Post output data.
+  PostData(1, fHists);
+}      
+
+Double_t AliAnalysisTaskKMeans::DeltaPhi(Double_t phi1, Double_t phi2)
+{
+    Double_t dphi = TMath::Abs(phi1 - phi2);
+    if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
+    return dphi;
+}
+
+Double_t AliAnalysisTaskKMeans::DeltaR(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
+{
+    Double_t dphi = DeltaPhi(phi1, phi2);
+    Double_t deta = eta1 - eta2;
+    return (TMath::Sqrt(dphi * dphi + deta * deta));
+    
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskKMeans::Terminate(Option_t *) 
+{
+}  
+
diff --git a/JETAN/AliAnalysisTaskKMeans.h b/JETAN/AliAnalysisTaskKMeans.h
new file mode 100644 (file)
index 0000000..1127c1e
--- /dev/null
@@ -0,0 +1,63 @@
+#ifndef AliAnalysisTaskKMeans_cxx
+#define AliAnalysisTaskKMeans_cxx
+ /* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+//-------------------------------------------------------------------------
+//     Analysis Task that uses the Soft K-Means Algorithm to find clusters in
+//     the eta-phi space of Minimum Bias. No pt information is used for the clustering.
+//     
+//
+//     Author: Andreas Morsch (CERN)
+//     andreas.morsch@cern.ch
+//-------------------------------------------------------------------------
+
+class TH1F;
+class TH2F;
+class TList;
+class TProfile;
+
+class AliESDEvent;
+class AliESDtrack;
+class AliESDtrackCuts;
+
+
+#include "AliAnalysisTaskSE.h"
+
+class AliAnalysisTaskKMeans : public AliAnalysisTaskSE {
+ public:
+    AliAnalysisTaskKMeans();
+  AliAnalysisTaskKMeans(const char *name);
+  virtual ~AliAnalysisTaskKMeans() {}
+  virtual void     UserCreateOutputObjects();
+  virtual void     UserExec(Option_t *option);
+  virtual void     Terminate(Option_t *);
+  virtual void     SetCuts(AliESDtrackCuts* cuts) {fCuts = cuts;}
+  virtual Double_t DeltaPhi(Double_t phi1, Double_t phi2);
+  virtual Double_t DeltaR(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2);
+  
+ private:
+  // Others
+  TList*           fHists;         // Histograms
+  TH1F*            fH1CEta;        // Eta distribution of clusters
+  TH1F*            fH1CPhi;        // Phi distribution of clusters  
+  TH1F*            fH1CEtaR;       // Eta distribution of clusters for rndm evnt
+  TH2F*            fH2N1N2;        // Cluster sizes 
+  TH1F*            fH1Pt;          // pt outside clusters
+  TH1F*            fH1PtC;         // pt outside clusters
+  TH1F*            fH1PtC1;        // pt dr > 0.4
+  TH1F*            fH1PtC2;        // pt dr > 0.2 
+  TH1F*            fH1DPhi;        // Dphi wr to cluster
+  TH1F*            fH1DR;          // DR   wr to cluster
+  TH1F*            fH1DRR;         // DR   wr to cluster from rndm events   
+  TH2F*            fH2DPhiEta;     // eta-phi wr to cluster
+  TH2F*            fH2DPhiEtaR;    // eta-phi wr to cluster for rndm events 
+  TH2F*            fH2DPhiEtaL;    // eta-phi of leading particle
+  
+  AliESDtrackCuts* fCuts;             // List of cuts
+  ClassDef(AliAnalysisTaskKMeans, 1); // A k-means clustering analysis
+};
+
+#endif
diff --git a/JETAN/AliKMeansClustering.cxx b/JETAN/AliKMeansClustering.cxx
new file mode 100755 (executable)
index 0000000..478a07c
--- /dev/null
@@ -0,0 +1,129 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+// Implemenatation of the K-Means Clustering Algorithm
+// See http://en.wikipedia.org/wiki/K-means_clustering and references therein.
+//
+// This particular implementation is the so called Soft K-means algorithm.
+// It has been modified to work on the cylindrical topology in eta-phi space.
+//
+// Author: Andreas Morsch (CERN)
+// andreas.morsch@cern.ch
+#include "AliKMeansClustering.h"
+#include <TMath.h>
+#include <TRandom.h>
+
+ClassImp(AliKMeansClustering)
+
+Double_t AliKMeansClustering::fBeta = 10.;
+void AliKMeansClustering::SoftKMeans(Int_t k, Int_t n, Double_t* x, Double_t* y, Double_t* mx, Double_t* my , Double_t* rk )
+{
+    //
+    // The soft K-means algorithm
+    //
+    Int_t i,j;
+    //
+    // (1) Initialisation of the k means
+
+    for (i = 0; i < k; i++) {
+       mx[i] = 2. * TMath::Pi() * gRandom->Rndm();
+       my[i] = -1. + 2. * gRandom->Rndm();
+    }
+
+    //
+    // (2a) The responsibilities
+    Double_t** r = new Double_t*[n]; // responsibilities
+    for (j = 0; j < n; j++) {r[j] = new Double_t[k];}
+    //
+    // (2b) Normalisation
+    Double_t* nr = new Double_t[n];
+
+    // (3) Iterations
+    Int_t nit = 0;
+    
+    while(1) {
+       nit++;
+      //
+      // Assignment step
+      //
+      for (j = 0; j < n; j++) {
+       nr[j] = 0.;
+       for (i = 0; i < k; i++) {
+         r[j][i] = TMath::Exp(- fBeta * d(mx[i], my[i], x[j], y[j]));
+         nr[j] += r[j][i];
+       } // mean i
+      } // data point j
+       
+      for (j = 0; j < n; j++) {
+       for (i = 0; i < k; i++) {
+         r[j][i] /=  nr[j];
+       } // mean i
+      } // data point j
+      
+       //
+       // Update step
+      Double_t di = 0;
+      
+      for (i = 0; i < k; i++) {
+         Double_t oldx = mx[i];
+         Double_t oldy = my[i];
+         
+         mx[i] = x[0];
+         my[i] = y[0];
+         rk[i] = r[0][i];
+       
+       for (j = 1; j < n; j++) {
+           Double_t xx =  x[j];
+//
+// Here we have to take into acount the cylinder topology where phi is defined mod 2xpi
+// If two coordinates are separated by more than pi in phi one has to be shifted by +/- 2 pi
+
+           Double_t dx = mx[i] - x[j];
+           if (dx >  TMath::Pi()) xx += 2. * TMath::Pi();
+           if (dx < -TMath::Pi()) xx -= 2. * TMath::Pi();
+           mx[i] = mx[i] * rk[i] + r[j][i] * xx;
+           my[i] = my[i] * rk[i] + r[j][i] * y[j];
+           rk[i] += r[j][i];
+           mx[i] /= rk[i];
+           my[i] /= rk[i];         
+           if (mx[i] > 2. * TMath::Pi()) mx[i] -= 2. * TMath::Pi();
+           if (mx[i] < 0.              ) mx[i] += 2. * TMath::Pi();
+       } // Data
+       di += d(mx[i], my[i], oldx, oldy);
+      } // means 
+       //
+       // ending condition
+      if (di < 1.e-8 || nit > 1000) break;
+    } // while
+
+// Clean-up    
+    delete[] nr;
+    for (j = 0; j < n; j++) delete[] r[j];
+    delete[] r;
+}
+
+Double_t AliKMeansClustering::d(Double_t mx, Double_t my, Double_t x, Double_t y)
+{
+    //
+    // Distance definition 
+    // Quasi - Euclidian on the eta-phi cylinder
+    
+    Double_t dx = TMath::Abs(mx-x);
+    if (dx > TMath::Pi()) dx = 2. * TMath::Pi() - dx;
+    
+    return (dx * dx + (my - y) * (my - y));
+}
diff --git a/JETAN/AliKMeansClustering.h b/JETAN/AliKMeansClustering.h
new file mode 100755 (executable)
index 0000000..86fd92d
--- /dev/null
@@ -0,0 +1,32 @@
+#ifndef ALIKMEANSCLUSTERING_H
+#define ALIKMEANSCLUSTERING_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// Implemenatation of the K-Means Clustering Algorithm
+// http://en.wikipedia.org/wiki/K-means_clustering
+// This particular implementation is the so called Soft K-means algorithm.
+// It has been modified to work on the cylindrical topology in eta-phi space.
+//
+// Author: Andreas Morsch (CERN)
+// andreas.morsch@cern.ch
+
+#include <TObject.h>
+class AliKMeansClustering : public TObject
+{
+ public:
+  AliKMeansClustering()          {}
+  virtual ~AliKMeansClustering() {}
+  
+  static void SoftKMeans(Int_t k, Int_t n, Double_t* x, Double_t* y, Double_t* mx, Double_t* my , Double_t* rk );
+  static void SetBeta(Double_t beta) {fBeta = beta;}
+  static Double_t d(Double_t mx, Double_t my, Double_t x, Double_t y);
+protected:
+  static Double_t fBeta;
+  
+  ClassDef(AliKMeansClustering, 1)
+};
+#endif
index 6387535..c86c228 100644 (file)
@@ -43,4 +43,6 @@
 #pragma        link C++ class AliJetHistos+;
 #pragma        link C++ class AliAnalysisTaskDiJets+;
 #pragma        link C++ class AliEventShape+;
+#pragma        link C++ class AliKMeansClustering+;
+#pragma        link C++ class AliAnalysisTaskKMeans+;
 #endif
index 05f4de4..8fde5ba 100644 (file)
@@ -26,7 +26,9 @@ SRCS =        AliJetHeader.cxx \
        AliJetMCReaderHeader.cxx \
        AliJetHistos.cxx \
        AliAnalysisTaskDiJets.cxx \
-       AliEventShape.cxx
+       AliEventShape.cxx \
+       AliKMeansClustering.cxx \
+       AliAnalysisTaskKMeans.cxx
 
 HDRS:= $(SRCS:.cxx=.h) 
 DHDR= JETANLinkDef.h
diff --git a/PWG4/macros/AddTaskKMeans.C b/PWG4/macros/AddTaskKMeans.C
new file mode 100644 (file)
index 0000000..8f4d34e
--- /dev/null
@@ -0,0 +1,46 @@
+AliAnalysisTaskKMeans *AddTaskKMeans()\r
+{\r
+// Creates a dijet task, configures it and adds it to the analysis manager.\r
+\r
+   // Get the pointer to the existing analysis manager via the static access method.\r
+   //==============================================================================\r
+   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
+   if (!mgr) {\r
+      ::Error("AddTaskJets", "No analysis manager to connect to.");\r
+      return NULL;\r
+   }\r
+\r
+   // Check the analysis type using the event handlers connected to the analysis manager.\r
+   //==============================================================================\r
+   if (!mgr->GetInputEventHandler()) {\r
+      ::Error("AddTaskKMeans", "This task requires an input event handler");\r
+      return NULL;\r
+   }\r
+\r
+   // Create the task and configure it.\r
+   //===========================================================================\r
+\r
+   AliAnalysisTaskKMeans *taskKMeans = new AliAnalysisTaskKMeans("K-Means Analysis");\r
+   taskKMeans->SetDebugLevel(0);\r
+   AliESDtrackCuts* esdTrackCutsL = new AliESDtrackCuts("AliESDtrackCuts", "Standard");\r
+   esdTrackCutsL->SetMinNClustersTPC(50);\r
+   esdTrackCutsL->SetRequireTPCRefit(kTRUE);\r
+   esdTrackCutsL->SetRequireITSRefit(kTRUE);\r
+   esdTrackCutsL->SetMaxDCAToVertexXY(3.);\r
+   esdTrackCutsL->SetMaxDCAToVertexZ(3.);\r
+   esdTrackCutsL->SetAcceptKinkDaughters(kFALSE);\r
+   taskKMeans->SetCuts(esdTrackCutsL);\r
+   mgr->AddTask(taskKMeans);\r
+\r
+   AliAnalysisDataContainer* cout_kmeans = mgr->CreateContainer("KMeans", TList::Class(),AliAnalysisManager::kOutputContainer,\r
+     Form("%s:PWG4_KMeans", AliAnalysisManager::GetCommonFileName()));\r
+\r
+   // Create ONLY the output containers for the data produced by the task.\r
+   // Get and connect other common input/output containers via the manager as below\r
+   //==============================================================================\r
+   mgr->ConnectInput  (taskKMeans, 0, mgr->GetCommonInputContainer());\r
+   mgr->ConnectOutput (taskKMeans, 0, mgr->GetCommonOutputContainer());\r
+   mgr->ConnectOutput (taskKMeans, 1, cout_kmeans);\r
+\r
+   return taskKMeans;\r
+}\r