- Chosing optimal number of clusters
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 12 Feb 2010 17:08:43 +0000 (17:08 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 12 Feb 2010 17:08:43 +0000 (17:08 +0000)
JETAN/AliAnalysisTaskKMeans.cxx
JETAN/AliAnalysisTaskKMeans.h
JETAN/AliKMeansClustering.cxx
JETAN/AliKMeansClustering.h

index f27ed45..191aba1 100644 (file)
@@ -37,6 +37,7 @@
 #include "TProfile.h"
 #include "TMath.h"
 #include "TRandom.h"
+#include "TObjArray.h"
 
 #include "AliAnalysisTask.h"
 #include "AliAnalysisManager.h"
@@ -91,6 +92,10 @@ AliAnalysisTaskKMeans::AliAnalysisTaskKMeans()
     ,fH1RespR(0)
     ,fH2Sigma(0)
     ,fH2SigmaR(0)
+    ,fHDensity(0)
+    ,fHCSize(0)
+    ,fHNCluster(0)
+    ,fHPtDensity(0)
     ,fCuts(0)
 {
   //
@@ -129,6 +134,10 @@ AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name)
       ,fH1RespR(0)
       ,fH2Sigma(0)
       ,fH2SigmaR(0)
+      ,fHDensity(0)
+      ,fHCSize(0)
+      ,fHNCluster(0)
+      ,fHPtDensity(0)
       ,fCuts(0)
 {
   //
@@ -172,6 +181,14 @@ void AliAnalysisTaskKMeans::UserCreateOutputObjects()
     fH1RespR      = new TH1F("fH1RespR",  "Responsibility", 50, 0., 1.);
     fH2Sigma      = new TH2F("fH2Sigma",  "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.);
     fH2SigmaR     = new TH2F("fH2SigmaR", "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.);
+
+
+    fHDensity    = new TH1F("fHDensity",    "density distribution", 100, 0., 20.);
+    fHCSize      = new TH1F("fHCSize",      "cluster size", 20, -0.5, 19.5);
+
+    fHNCluster   = new TH1F("fHNCluster",   "Number of clusters", 11, -0.5, 10.5);
+    fHPtDensity  = new TH2F("fHPtDensity", "Pt vs density", 100, 0., 20., 50, 0., 10.);
+
  
     fHists->SetOwner();
 
@@ -199,10 +216,16 @@ void AliAnalysisTaskKMeans::UserCreateOutputObjects()
     fHists->Add(fH1RespR);
     fHists->Add(fH1Resp);    
     fHists->Add(fH2Sigma);    
-    fHists->Add(fH2SigmaR);    
-    
+    fHists->Add(fH2SigmaR);
+    fHists->Add(fHCSize);    
+    fHists->Add(fHDensity);
+    fHists->Add(fHNCluster);
+    fHists->Add(fHPtDensity);
     //
-    
+    for (Int_t i = 0; i < 10; i++) {
+      fA[i] = new AliKMeansResult(i+1);
+      fB[i] = new AliKMeansResult(i+1);
+    }
 }
 
 //________________________________________________________________________
@@ -212,16 +235,9 @@ void AliAnalysisTaskKMeans::UserExec(Option_t *)
   // Called for each event
     Double_t   phi [500];
     Double_t   phiR[500];
-    Double_t   eta[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];
-    Double_t sigma2[20];
-    Double_t sigmax2[20];
-    Double_t sigmay2[20];
-    Int_t       ind[20];
+    Double_t    pt [500];
 
   if (!fInputEvent) {
     Printf("ERROR: fESD not available");
@@ -230,17 +246,11 @@ void AliAnalysisTaskKMeans::UserExec(Option_t *)
   
 
   Int_t ic = 0;
+  
   //
   // Fill eta-phi positions
   //
-  /*
-  const AliMultiplicity *spdMult = fInputEvent->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;
   
@@ -248,11 +258,11 @@ void AliAnalysisTaskKMeans::UserExec(Option_t *)
       AliVParticle* track = fInputEvent->GetTrack(iTracks);
       if ((fCuts->AcceptTrack((AliESDtrack*)track))) 
       {
-         phi[ic]  = track->Phi();
+         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();
+         eta [ic] = track->Eta();
+         etaR[ic] = 1.8 * gRandom->Rndm() - 0.9;
+         pt  [ic] = track->Pt();
          if (pt[ic] > ptmax) {
              ptmax = pt[ic];
              icmax = ic;
@@ -268,30 +278,55 @@ void AliAnalysisTaskKMeans::UserExec(Option_t *)
       return;
   }
   //
-  Int_t ok;
+  Double_t* rk0 = new Double_t[10];
+
+  AliKMeansResult* res = 0;
+
+  for (Int_t i = 0; i < fK; i++) {
+    res = fA[i];
+    AliKMeansClustering::SoftKMeans2(i+1, ic, phi, eta, res->GetMx(),res->GetMy(), res->GetSigma2(), res->GetRk());
+    res->Sort(ic, phi, eta);
+    Int_t j = (res->GetInd())[0];
+    rk0[i]  = (res->GetTarget())[j];
+  }
+
+  Double_t* mPhi   = 0;
+  Double_t* mEta   = 0;
+  Double_t* sigma2 = 0;
+  Float_t   rmax   = -1.;
+  Int_t     imax   = 0;
+
+  for (Int_t i = 0; i < fK; i++) {
+    if (rk0[i] > rmax)  {
+       rmax = rk0[i];
+       imax = i;
+      }
+  }
+
+
+  res = fA[imax];
+  mPhi    = res->GetMx();
+  mEta    = res->GetMy();
+  sigma2  = res->GetSigma2();
+  Int_t nk = imax + 1;
+  Int_t im = (res->GetInd())[0];
+  fHDensity->Fill(rmax / TMath::Pi());
+  fHCSize->Fill(rmax * sigma2[im]);
+  fHNCluster->Fill(Float_t(nk));
+
   Double_t dphic, detac;
-  //
-  ok = AliKMeansClustering::SoftKMeans3(fK, ic, phi, eta, mPhi, mEta, sigmax2, sigmay2, rk);
-  if (ok) {
-    //
-    // Sort
-    for (Int_t i = 0; i < fK; i++) {
-      if (rk[i] > 1.) rk[i] /= TMath::Sqrt(sigmax2[i] * sigmay2[i]);
-      else rk[i] = 0.;
-    }    
 
-    TMath::Sort(fK, rk, ind);
-    fH1Resp->Fill(rk[ind[0]]/(rk[0]+rk[1]));
-    fH2Sigma->Fill(TMath::Sqrt(sigmax2[ind[0]]), TMath::Sqrt(sigmay2[ind[0]])); 
-    //
+  if (rmax > 0.) {
     // Analysis
     //
     // Cluster Multiplicity
     Int_t mult[2] = {0, 0};
     //
-    dphic = DeltaPhi(mPhi[0], mPhi[1]);
-    detac = TMath::Abs(mEta[0] - mEta[1]);
-    fH2DPhiEtaC->Fill(dphic, detac);  
+    if (nk > 1) {
+      dphic = DeltaPhi(mPhi[0], mPhi[1]);
+      detac = TMath::Abs(mEta[0] - mEta[1]);
+      fH2DPhiEtaC->Fill(dphic, detac);  
+    }
     //
     // Random cluster position
     Int_t ir = Int_t(Float_t(ic) * gRandom->Rndm());
@@ -301,14 +336,16 @@ void AliAnalysisTaskKMeans::UserExec(Option_t *)
     //
     Double_t sumPt  = 0;
     Double_t sumPtC = 0;
+
     for (Int_t i = 0; i < 1; i++) {
-      fH1CEta->Fill(mEta[ind[i]]);
-      fH1CPhi->Fill(mPhi[ind[i]]);      
+      fH1CEta->Fill(mEta[im]);
+      fH1CPhi->Fill(mPhi[im]);      
       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];
+       Double_t r    = DeltaR(mPhi[im], mEta[im], phi[j], eta[j]);
+       Double_t dphi = DeltaPhi(mPhi[im], phi[j]);
+       Double_t deta = mEta[im] - eta[j];
        Double_t rr   = DeltaR(crPhi, crEta, phi[j], eta[j]);
+
        fH1DR->Fill(r);
        fH1DPhi->Fill(dphi);
        fH2DPhiEta->Fill(dphi, TMath::Abs(deta));
@@ -317,29 +354,27 @@ void AliAnalysisTaskKMeans::UserExec(Option_t *)
        if (r < 0.2) {
          fH1PtC2->Fill(pt[j]);
        }
-       if (r < 0.3) 
-         {
-           fH1PtC1->Fill(pt[j]);
-         }
-       if (rr < 0.3) 
-         {
+       if (r < 0.3) {
+         fH1PtC1->Fill(pt[j]);
+         fHPtDensity->Fill(rmax/TMath::Pi(), pt[j]);
+       }
+       if (rr < 0.3) {
            fH1PtR->Fill(pt[j]);
-         }
+       }
 
-       if (r < 0.4) 
-         {
-           sumPtC += pt[j];
-           mult[i]++;
-           fH1PtC->Fill(pt[j]);
-         } 
-       if (r > 0.7) {
+       if (r < 0.4)  {
+         sumPtC += pt[j];
+         mult[i]++;
+         fH1PtC->Fill(pt[j]);
+       } 
+       if (r > 0.7 && dphi < (TMath::Pi() - 0.3)) {
          fH1Pt->Fill(pt[j]);
        }
        
        if (r > 0.7 && r < 1.) {
          sumPt += pt[j];
        }
-
+       
        if (dphi > (TMath::Pi() - 0.3)) {
          fH1PtAS->Fill(pt[j]);
        }
@@ -349,48 +384,60 @@ void AliAnalysisTaskKMeans::UserExec(Option_t *)
     fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
     fH1SPt ->Fill(sumPt);
     fH1SPtC->Fill(sumPtC);
-  } // ok
-
-  // Randomized phi
-  //
-  ok = AliKMeansClustering::SoftKMeans3(fK, ic, phiR, etaR, mPhi, mEta, sigmax2, sigmay2, rk);
-  if (ok) {
-    //
-    // Sort according to responsibility density
-    for (Int_t i = 0; i < fK; i++) {
-      if (rk[i] > 1.) rk[i] /= TMath::Sqrt(sigmax2[i] * sigmay2[i]);
-      else rk[i] = 0.;
-    }    
-
-    TMath::Sort(fK, rk, ind);
-    fH1RespR->Fill(rk[ind[0]]/(rk[0]+rk[1]));
-    fH2SigmaR->Fill(TMath::Sqrt(sigmax2[ind[0]]), TMath::Sqrt(sigmay2[ind[0]])); 
+  }
     //
-    // Analyse
+    // Randomized phi
     //
-    // Cluster Multiplicity
-    for (Int_t i = 0; i < 1; i++) {
-      fH1CEtaR->Fill(mEta[ind[i]]);
+    for (Int_t i = 0; i < fK; i++) {
+      res = fB[i];
+      AliKMeansClustering::SoftKMeans2(i+1, ic, phiR, etaR, res->GetMx(),res->GetMy(), res->GetSigma2(), res->GetRk());
+      res->Sort(ic, phiR, etaR);
+      //res->Sort();
+      Int_t j = (res->GetInd())[0];
+      rk0[i]  = (res->GetTarget())[j];
     }
     
-    for (Int_t i = 0; i < 1; i++) {
-      for (Int_t j = 0; j < ic; j++) {
-       Double_t dphi = DeltaPhi(mPhi[ind[i]], phiR[j]);
-       Double_t deta = mEta[ind[i]] - etaR[j];
-       Double_t r = DeltaR(mPhi[ind[i]], mEta[ind[i]], phiR[j], etaR[j]);
-       fH1DRR->Fill(r);
-       fH2DPhiEtaR->Fill(dphi, TMath::Abs(deta));
-       if (j == icmax) fH2DPhiEtaLR->Fill(dphi, TMath::Abs(deta));
+    rmax   = -1.;
+    imax   = 0;
+    for (Int_t i = 0; i < fK; i++) {
+      if (rk0[i] > rmax) {
+         rmax = rk0[i];
+         imax = i;
+       }
+    }
+    res = fB[imax];
+    mPhi    = res->GetMx();
+    mEta    = res->GetMy();
+    sigma2  = res->GetSigma2();
+    nk = imax + 1;
+    
+    //  fH1RespR->Fill(rk[ind[0]]/(rk[0]+rk[1]));
+    //
+    // Cluster Multiplicity
+    if (rmax > 0.) {
+      for (Int_t i = 0; i < 1; i++) {
+       Int_t im = (res->GetInd())[i];
+       fH1CEtaR->Fill(mEta[im]);
+      }
+      
+      for (Int_t i = 0; i < 1; i++) {
+       Int_t im = (res->GetInd())[i];
+       for (Int_t j = 0; j < ic; j++) {
+         Double_t dphi = DeltaPhi(mPhi[im], phiR[j]);
+         Double_t deta = mEta[im] - etaR[j];
+         Double_t r = DeltaR(mPhi[im], mEta[im], phiR[j], etaR[j]);
+         fH1DRR->Fill(r);
+         fH2DPhiEtaR->Fill(dphi, TMath::Abs(deta));
+         if (j == icmax) fH2DPhiEtaLR->Fill(dphi, TMath::Abs(deta));
+       }
+      }
+      if (nk > 1) {
+       dphic = DeltaPhi(mPhi[0], mPhi[1]);
+       detac = TMath::Abs(mEta[0] - mEta[1]);
+       fH2DPhiEtaCR->Fill(dphic, detac);  
       }
     }
-    dphic = DeltaPhi(mPhi[0], mPhi[1]);
-    detac = TMath::Abs(mEta[0] - mEta[1]);
-    fH2DPhiEtaCR->Fill(dphic, detac);  
-  } // ok
-  
-  //
-  // Post output data.
-  PostData(1, fHists);
+    PostData(1, fHists);
 }      
 
 Double_t AliAnalysisTaskKMeans::DeltaPhi(Double_t phi1, Double_t phi2)
index c1b6d59..a6dd6f1 100644 (file)
@@ -18,11 +18,10 @@ class TH1F;
 class TH2F;
 class TList;
 class TProfile;
-
 class AliESDEvent;
 class AliESDtrack;
 class AliESDtrackCuts;
-
+class AliKMeansResult;
 
 #include "AliAnalysisTaskSE.h"
 
@@ -43,32 +42,38 @@ class AliAnalysisTaskKMeans : public AliAnalysisTaskSE {
   // Others
   Int_t            fK;             // K                        
   Int_t            fNMin;          // Minimum multipicity                         
-  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*            fH1PtAS;        // away-side peak 
-  TH1F*            fH1PtR;         // away-side peak 
-  TH1F*            fH1SPt;         // sum pt
-  TH1F*            fH1SPtC;        // sum pt
-  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
-  TH2F*            fH2DPhiEtaLR;   // eta-phi of leading particle
-  TH2F*            fH2DPhiEtaC;    // eta-phi of Clusters
-  TH2F*            fH2DPhiEtaCR;   // eta-phi of Clusters
-  TH1F*            fH1Resp;        // responsibility
-  TH1F*            fH1RespR;       // responsibility
-  TH2F*            fH2Sigma;
-  TH2F*            fH2SigmaR;
+  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*            fH1PtAS;        //! away-side peak 
+  TH1F*            fH1PtR;         //! away-side peak 
+  TH1F*            fH1SPt;         //! sum pt
+  TH1F*            fH1SPtC;        //! sum pt
+  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
+  TH2F*            fH2DPhiEtaLR;   //! eta-phi of leading particle
+  TH2F*            fH2DPhiEtaC;    //! eta-phi of Clusters
+  TH2F*            fH2DPhiEtaCR;   //! eta-phi of Clusters
+  TH1F*            fH1Resp;        //! responsibility
+  TH1F*            fH1RespR;       //! responsibility
+  TH2F*            fH2Sigma;       //! sigma
+  TH2F*            fH2SigmaR;      //! sigma random
+  TH1F*            fHDensity;      //! Particle density
+  TH1F*            fHCSize;        //! Cluster Size
+  TH1F*            fHNCluster;     //! Number of clusters
+  TH2F*            fHPtDensity;    //! Pt vs density
+  AliKMeansResult* fA[10];         //!
+  AliKMeansResult* fB[10];         //!
   AliESDtrackCuts* fCuts;             // List of cuts
   ClassDef(AliAnalysisTaskKMeans, 1); // A k-means clustering analysis
 };
index 948c26a..9947da7 100755 (executable)
@@ -440,24 +440,62 @@ void AliKMeansClustering::OptimalInit(Int_t k, Int_t n, Double_t* x, Double_t* y
 }
 
 
-Int_t AliKMeansClustering::NTwoSigma(Int_t k, Int_t n, Double_t* x, Double_t* y, Double_t* mx, Double_t* my, 
-                                          Double_t* sigma2x, Double_t* sigma2y)
+ClassImp(AliKMeansResult)
+
+
+    
+AliKMeansResult::AliKMeansResult(Int_t k):
+    TObject(),
+    fK(k),
+    fMx    (new Double_t[k]),
+    fMy    (new Double_t[k]),
+    fSigma2(new Double_t[k]),
+    fRk    (new Double_t[k]),
+    fTarget(new Double_t[k]),
+    fInd   (new Int_t[k])
 {
-  //  
-  // Counting the number of data points within 2sigma of a cluster
-  //
-  Int_t n2sigma = 0;
-  for (Int_t j = 0; j < n; j++) {
-    Double_t nassign = 0;
-    for (Int_t i = 0; i < k; i++) {
-      Double_t dx = TMath::Abs(mx[i]-x[j]);
-      if (dx > TMath::Pi()) dx = 2. * TMath::Pi() - dx; 
-      Double_t dy = TMath::Abs(my[i]-y[j]);
-      if ((dx * dx / sigma2x[i] + dy * dy / sigma2y[i]) < 4.) nassign++;
-    } // clusters
-    if (nassign > 0) n2sigma++;
-  } // data points
-  return (n2sigma);
+// Constructor
+}
+
+AliKMeansResult::~AliKMeansResult()
+{
+// Destructor
+    delete[] fMx;
+    delete[] fMy;
+    delete[] fSigma2;
+    delete[] fRk;
+    delete[] fInd;
+    delete[] fTarget;
+}
+
+void AliKMeansResult::Sort()
+{
+// Sort clusters
+    for (Int_t i = 0; i < fK; i++) {
+       if (fRk[i] > 1.) fRk[i] /= fSigma2[i];
+       else fRk[i] = 0.;
+    }
+    
+    TMath::Sort(fK, fRk, fInd);
+    
 }
 
+void AliKMeansResult::Sort(Int_t n, Double_t* x, Double_t* y)
+{
+  // Build target array and sort
+  for (Int_t i = 0; i < fK; i++)
+    {
+      Int_t nc = 0;
+      for (Int_t j = 0; j < n; j++)
+       {
+         if (2. * AliKMeansClustering::d(fMx[i], fMy[i], x[j], y[j])  <  fSigma2[i]) nc++;
+       }
+      if (nc > 1) {
+       fTarget[i] = Double_t(nc) / fSigma2[i];
+      } else {
+       fTarget[i] = 0.;
+      }
+    }
 
+  TMath::Sort(fK, fTarget, fInd);
+}
index 56b50ed..502ca99 100755 (executable)
@@ -28,12 +28,37 @@ class AliKMeansClustering : public TObject
   static void  OptimalInit(Int_t k, Int_t n, Double_t* x, Double_t* y, Double_t* mx, Double_t* my);
   static void  SetBeta(Double_t beta) {fBeta = beta;}
   static Double_t d(Double_t mx, Double_t my, Double_t x, Double_t y);
-  static Int_t NTwoSigma(Int_t k, Int_t n, Double_t* x, Double_t* y, Double_t* mx, Double_t* my, 
-                               Double_t* sigma2x, Double_t* sigma2y);
 protected:
   static Double_t fBeta;
   
   ClassDef(AliKMeansClustering, 1)
 };
+
+class AliKMeansResult : public TObject
+{
+ public:
+  AliKMeansResult(Int_t k);
+  virtual ~AliKMeansResult();
+  Int_t      GetK()      const  {return fK;}
+  Double_t*  GetMx()     const  {return fMx;}
+  Double_t*  GetMy()     const  {return fMy;}
+  Double_t*  GetSigma2() const  {return fSigma2;}
+  Double_t*  GetRk()     const  {return fRk;}
+  Int_t*     GetInd()    const  {return fInd;}
+  Double_t*  GetTarget() const  {return fTarget;}
+  void       Sort();
+  void       Sort(Double_t* target);
+  void       Sort(Int_t n, Double_t* x, Double_t* y);  
+protected:
+  Int_t        fK;        //! Number of clusters
+  Double_t*    fMx;       //! Position x
+  Double_t*    fMy;       //! Position y
+  Double_t*    fSigma2;   //! Sigma2
+  Double_t*    fRk;       //! Responsibility
+  Double_t*    fTarget;   //! Target for sorting
+  Int_t*       fInd;      //! Index for sorting
  
+  ClassDef(AliKMeansResult, 1)
+};
+
 #endif