[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskKMeans.cxx

/**************************************************************************
 * 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 "TObjArray.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliVEvent.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()
    ,fK(0)
    ,fNMin(0)
    ,fHists(0)
    ,fH1CEta(0)
    ,fH1CPhi(0)
    ,fH1CEtaR(0)
    ,fH2N1N2(0)
    ,fH1Pt(0)
    ,fH1PtC(0)
    ,fH1PtC1(0)
    ,fH1PtC2(0)
    ,fH1PtAS(0)
    ,fH1PtR(0)
    ,fH1SPt(0)
    ,fH1SPtC(0)
    ,fH1DPhi(0)
    ,fH1DR(0)
    ,fH1DRR(0)
    ,fH2DPhiEta(0)
    ,fH2DPhiEtaR(0)
    ,fH2DPhiEtaL(0)
    ,fH2DPhiEtaLR(0)
    ,fH2DPhiEtaC(0)
    ,fH2DPhiEtaCR(0)
    ,fH1Resp(0)
    ,fH1RespR(0)
    ,fH2Sigma(0)
    ,fH2SigmaR(0)
    ,fHDensity(0)
    ,fHCSize(0)
    ,fHNCluster(0)
    ,fHPtDensity(0)
    ,fCuts(0)
{
  //
  // Constructor
  //
}

//________________________________________________________________________
AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name) 
    : AliAnalysisTaskSE(name) 
      ,fK(0)
      ,fNMin(0)
      ,fHists(0)
      ,fH1CEta(0)
      ,fH1CPhi(0)
      ,fH1CEtaR(0)
      ,fH2N1N2(0)
      ,fH1Pt(0)
      ,fH1PtC(0)
      ,fH1PtC1(0)
      ,fH1PtC2(0)
      ,fH1PtAS(0)
      ,fH1PtR(0)
      ,fH1SPt(0)
      ,fH1SPtC(0)
      ,fH1DPhi(0)
      ,fH1DR(0)
      ,fH1DRR(0)
      ,fH2DPhiEta(0)
      ,fH2DPhiEtaR(0)
      ,fH2DPhiEtaL(0)
      ,fH2DPhiEtaLR(0)
      ,fH2DPhiEtaC(0)
      ,fH2DPhiEtaCR(0)
      ,fH1Resp(0)
      ,fH1RespR(0)
      ,fH2Sigma(0)
      ,fH2SigmaR(0)
      ,fHDensity(0)
      ,fHCSize(0)
      ,fHNCluster(0)
      ,fHPtDensity(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.);
    fH1PtAS  = new TH1F("fH1PtAS",   "pt distribution",50, 0., 10.);
    fH1PtR   = new TH1F("fH1PtR",    "pt distribution",50, 0., 10.);

    fH1SPt    = new TH1F("fH1SPt",     "sum pt distribution",50, 0., 10.);
    fH1SPtC   = new TH1F("fH1SPtC",    "sum 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.);
    fH2DPhiEtaLR  = new TH2F("fH2DPhiEtaLR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
    fH2DPhiEtaC   = new TH2F("fH2DPhiEtaC","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
    fH2DPhiEtaCR  = new TH2F("fH2DPhiEtaCR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
    fH1DRR        = new TH1F("fH1DRR",    "dR distribution", 50, 0., 5.);
    fH1Resp       = new TH1F("fH1Resp",   "Responsibility", 50, 0., 1.);
    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();

    fHists->Add(fH1CEta);
    fHists->Add(fH1CEtaR);
    fHists->Add(fH1CPhi);
    fHists->Add(fH2N1N2);
    fHists->Add(fH1Pt);
    fHists->Add(fH1PtR);
    fHists->Add(fH1PtC);
    fHists->Add(fH1PtC1);
    fHists->Add(fH1PtC2);
    fHists->Add(fH1PtAS);
    fHists->Add(fH1DR);
    fHists->Add(fH1SPtC);
    fHists->Add(fH1SPt);
    fHists->Add(fH1DPhi);
    fHists->Add(fH2DPhiEta);
    fHists->Add(fH2DPhiEtaR);
    fHists->Add(fH2DPhiEtaL);
    fHists->Add(fH2DPhiEtaLR);
    fHists->Add(fH2DPhiEtaC);
    fHists->Add(fH2DPhiEtaCR);
    fHists->Add(fH1DRR);
    fHists->Add(fH1RespR);
    fHists->Add(fH1Resp);    
    fHists->Add(fH2Sigma);    
    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);
    }
}

//________________________________________________________________________
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];

  if (!fInputEvent) {
    Printf("ERROR: fESD not available");
    return;
  }
  

  Int_t ic = 0;
  
  //
  // Fill eta-phi positions
  //

  Double_t ptmax = 0.;
  Int_t    icmax = -1;
  
  for (Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++) {
      AliVParticle* track = fInputEvent->GetTrack(iTracks);
      if ((fCuts->AcceptTrack((AliESDtrack*)track))) 
      {
	  phi [ic] = track->Phi();
	  phiR[ic] = 2. * TMath::Pi() * gRandom->Rndm();
	  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;
	  }
	  ic++;
      }
  } //track loop 

  //
  // Cluster
  if (ic < fNMin) {
      PostData(1, fHists);
      return;
  }
  //
  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;

  if (rmax > 0.) {
    // Analysis
    //
    // Cluster Multiplicity
    Int_t mult[2] = {0, 0};
    //
    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());

    Double_t crPhi = phi[ir];
    Double_t crEta = eta[ir];
    //
    Double_t sumPt  = 0;
    Double_t sumPtC = 0;

    for (Int_t i = 0; i < 1; i++) {
      fH1CEta->Fill(mEta[im]);
      fH1CPhi->Fill(mPhi[im]);      
      for (Int_t j = 0; j < ic; 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));
	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]);
	  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 && 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]);
	}
      }
    }

    fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
    fH1SPt ->Fill(sumPt);
    fH1SPtC->Fill(sumPtC);
  }
    //
    // Randomized phi
    //
    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];
    }
    
    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);  
      }
    }
    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 *) 
{
}
Commit	Line	Data
70f2ce9d	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 "TTree.h"
	30	#include "TH1F.h"
	31	#include "TH2F.h"
	32	#include "TProfile.h"
	33
	34	#include "TList.h"
	35	#include "TParticle.h"
	36	#include "TParticlePDG.h"
	37	#include "TProfile.h"
	38	#include "TMath.h"
	39	#include "TRandom.h"
19c5a36a	40	#include "TObjArray.h"
70f2ce9d	41
	42	#include "AliAnalysisTask.h"
	43	#include "AliAnalysisManager.h"
	44
42dc9410	45	#include "AliVEvent.h"
70f2ce9d	46	#include "AliESDEvent.h"
	47	#include "AliStack.h"
	48	#include "AliESDVertex.h"
	49	#include "AliESDInputHandler.h"
	50	#include "AliESDtrackCuts.h"
	51	#include "AliMultiplicity.h"
	52
	53	#include "AliMCParticle.h"
	54	#include "AliMCEvent.h"
	55	#include "AliAnalysisTaskKMeans.h"
	56	#include "AliTrackReference.h"
	57	#include "AliStack.h"
	58	#include "AliHeader.h"
	59	#include "AliKMeansClustering.h"
	60
	61
	62
	63	ClassImp(AliAnalysisTaskKMeans)
	64
	65	AliAnalysisTaskKMeans::AliAnalysisTaskKMeans()
320bb308	66	: AliAnalysisTaskSE()
320bb308	67	,fK(0)
42dc9410	68	,fNMin(0)
70f2ce9d	69	,fHists(0)
	70	,fH1CEta(0)
	71	,fH1CPhi(0)
	72	,fH1CEtaR(0)
	73	,fH2N1N2(0)
	74	,fH1Pt(0)
	75	,fH1PtC(0)
	76	,fH1PtC1(0)
	77	,fH1PtC2(0)
42dc9410	78	,fH1PtAS(0)
42dc9410	79	,fH1PtR(0)
a51f90df	80	,fH1SPt(0)
a51f90df	81	,fH1SPtC(0)
70f2ce9d	82	,fH1DPhi(0)
	83	,fH1DR(0)
	84	,fH1DRR(0)
	85	,fH2DPhiEta(0)
	86	,fH2DPhiEtaR(0)
	87	,fH2DPhiEtaL(0)
5a9b5891	88	,fH2DPhiEtaLR(0)
a51f90df	89	,fH2DPhiEtaC(0)
a51f90df	90	,fH2DPhiEtaCR(0)
320bb308	91	,fH1Resp(0)
320bb308	92	,fH1RespR(0)
42dc9410	93	,fH2Sigma(0)
42dc9410	94	,fH2SigmaR(0)
19c5a36a	95	,fHDensity(0)
	96	,fHCSize(0)
	97	,fHNCluster(0)
	98	,fHPtDensity(0)
70f2ce9d	99	,fCuts(0)
	100	{
	101	//
	102	// Constructor
	103	//
	104	}
	105
	106	//________________________________________________________________________
	107	AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name)
	108	: AliAnalysisTaskSE(name)
320bb308	109	,fK(0)
42dc9410	110	,fNMin(0)
70f2ce9d	111	,fHists(0)
	112	,fH1CEta(0)
	113	,fH1CPhi(0)
	114	,fH1CEtaR(0)
	115	,fH2N1N2(0)
	116	,fH1Pt(0)
	117	,fH1PtC(0)
	118	,fH1PtC1(0)
	119	,fH1PtC2(0)
42dc9410	120	,fH1PtAS(0)
42dc9410	121	,fH1PtR(0)
a51f90df	122	,fH1SPt(0)
a51f90df	123	,fH1SPtC(0)
70f2ce9d	124	,fH1DPhi(0)
	125	,fH1DR(0)
	126	,fH1DRR(0)
	127	,fH2DPhiEta(0)
	128	,fH2DPhiEtaR(0)
	129	,fH2DPhiEtaL(0)
5a9b5891	130	,fH2DPhiEtaLR(0)
a51f90df	131	,fH2DPhiEtaC(0)
a51f90df	132	,fH2DPhiEtaCR(0)
320bb308	133	,fH1Resp(0)
320bb308	134	,fH1RespR(0)
42dc9410	135	,fH2Sigma(0)
42dc9410	136	,fH2SigmaR(0)
19c5a36a	137	,fHDensity(0)
	138	,fHCSize(0)
	139	,fHNCluster(0)
	140	,fHPtDensity(0)
70f2ce9d	141	,fCuts(0)
	142	{
	143	//
	144	// Constructor
	145	//
	146	DefineOutput(1, TList::Class());
	147	}
	148
	149
	150	//________________________________________________________________________
	151	void AliAnalysisTaskKMeans::UserCreateOutputObjects()
	152	{
	153	// Create histograms
	154	// Called once
	155	fHists = new TList();
	156	fH1CEta = new TH1F("fH1CEta", "eta distribution of clusters", 90, -1.5, 1.5);
	157	fH1CEtaR = new TH1F("fH1CEtaR", "eta distribution of clusters", 90, -1.5, 1.5);
	158	fH1CPhi = new TH1F("fH1CPhi", "phi distribution of clusters", 157, 0.0, 2. * TMath::Pi());
	159	fH2N1N2 = new TH2F("fH2N1N2", "multiplicity distribution", 50, 0., 50., 50, 0., 50.);
	160
	161	fH1Pt = new TH1F("fH1Pt", "pt distribution",50, 0., 10.);
	162	fH1PtC = new TH1F("fH1PtC", "pt distribution",50, 0., 10.);
	163	fH1PtC1 = new TH1F("fH1PtC1", "pt distribution",50, 0., 10.);
	164	fH1PtC2 = new TH1F("fH1PtC2", "pt distribution",50, 0., 10.);
42dc9410	165	fH1PtAS = new TH1F("fH1PtAS", "pt distribution",50, 0., 10.);
42dc9410	166	fH1PtR = new TH1F("fH1PtR", "pt distribution",50, 0., 10.);
70f2ce9d	167
a51f90df	168	fH1SPt = new TH1F("fH1SPt", "sum pt distribution",50, 0., 10.);
	169	fH1SPtC = new TH1F("fH1SPtC", "sum pt distribution",50, 0., 10.);
	170
5a9b5891	171	fH1DR = new TH1F("fH1DR", "dR distribution", 50, 0., 5.);
	172	fH1DPhi = new TH1F("fH1DPhi", "dPhi distribution", 31, 0., TMath::Pi());
	173	fH2DPhiEta = new TH2F("fH2DPhiEta","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
	174	fH2DPhiEtaR = new TH2F("fH2DPhiEtaR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
	175	fH2DPhiEtaL = new TH2F("fH2DPhiEtaL","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
	176	fH2DPhiEtaLR = new TH2F("fH2DPhiEtaLR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
	177	fH2DPhiEtaC = new TH2F("fH2DPhiEtaC","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
	178	fH2DPhiEtaCR = new TH2F("fH2DPhiEtaCR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
	179	fH1DRR = new TH1F("fH1DRR", "dR distribution", 50, 0., 5.);
320bb308	180	fH1Resp = new TH1F("fH1Resp", "Responsibility", 50, 0., 1.);
320bb308	181	fH1RespR = new TH1F("fH1RespR", "Responsibility", 50, 0., 1.);
42dc9410	182	fH2Sigma = new TH2F("fH2Sigma", "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.);
42dc9410	183	fH2SigmaR = new TH2F("fH2SigmaR", "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.);
19c5a36a	184
	185
	186	fHDensity = new TH1F("fHDensity", "density distribution", 100, 0., 20.);
	187	fHCSize = new TH1F("fHCSize", "cluster size", 20, -0.5, 19.5);
	188
	189	fHNCluster = new TH1F("fHNCluster", "Number of clusters", 11, -0.5, 10.5);
	190	fHPtDensity = new TH2F("fHPtDensity", "Pt vs density", 100, 0., 20., 50, 0., 10.);
	191
42dc9410	192
70f2ce9d	193	fHists->SetOwner();
	194
	195	fHists->Add(fH1CEta);
	196	fHists->Add(fH1CEtaR);
	197	fHists->Add(fH1CPhi);
	198	fHists->Add(fH2N1N2);
	199	fHists->Add(fH1Pt);
42dc9410	200	fHists->Add(fH1PtR);
70f2ce9d	201	fHists->Add(fH1PtC);
	202	fHists->Add(fH1PtC1);
	203	fHists->Add(fH1PtC2);
42dc9410	204	fHists->Add(fH1PtAS);
70f2ce9d	205	fHists->Add(fH1DR);
a51f90df	206	fHists->Add(fH1SPtC);
a51f90df	207	fHists->Add(fH1SPt);
70f2ce9d	208	fHists->Add(fH1DPhi);
	209	fHists->Add(fH2DPhiEta);
	210	fHists->Add(fH2DPhiEtaR);
	211	fHists->Add(fH2DPhiEtaL);
5a9b5891	212	fHists->Add(fH2DPhiEtaLR);
a51f90df	213	fHists->Add(fH2DPhiEtaC);
a51f90df	214	fHists->Add(fH2DPhiEtaCR);
70f2ce9d	215	fHists->Add(fH1DRR);
320bb308	216	fHists->Add(fH1RespR);
320bb308	217	fHists->Add(fH1Resp);
42dc9410	218	fHists->Add(fH2Sigma);
19c5a36a	219	fHists->Add(fH2SigmaR);
	220	fHists->Add(fHCSize);
	221	fHists->Add(fHDensity);
	222	fHists->Add(fHNCluster);
	223	fHists->Add(fHPtDensity);
70f2ce9d	224	//
19c5a36a	225	for (Int_t i = 0; i < 10; i++) {
	226	fA[i] = new AliKMeansResult(i+1);
	227	fB[i] = new AliKMeansResult(i+1);
	228	}
70f2ce9d	229	}
	230
	231	//________________________________________________________________________
	232	void AliAnalysisTaskKMeans::UserExec(Option_t *)
	233	{
	234	// Main loop
	235	// Called for each event
	236	Double_t phi [500];
	237	Double_t phiR[500];
19c5a36a	238	Double_t eta[500];
70f2ce9d	239	Double_t etaR[500];
19c5a36a	240	Double_t pt [500];
70f2ce9d	241
	242	if (!fInputEvent) {
	243	Printf("ERROR: fESD not available");
	244	return;
	245	}
	246
70f2ce9d	247
70f2ce9d	248	Int_t ic = 0;
19c5a36a	249
70f2ce9d	250	//
	251	// Fill eta-phi positions
	252	//
19c5a36a	253
70f2ce9d	254	Double_t ptmax = 0.;
	255	Int_t icmax = -1;
	256
42dc9410	257	for (Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++) {
	258	AliVParticle* track = fInputEvent->GetTrack(iTracks);
	259	if ((fCuts->AcceptTrack((AliESDtrack*)track)))
70f2ce9d	260	{
19c5a36a	261	phi [ic] = track->Phi();
70f2ce9d	262	phiR[ic] = 2. * TMath::Pi() * gRandom->Rndm();
19c5a36a	263	eta [ic] = track->Eta();
	264	etaR[ic] = 1.8 * gRandom->Rndm() - 0.9;
	265	pt [ic] = track->Pt();
70f2ce9d	266	if (pt[ic] > ptmax) {
	267	ptmax = pt[ic];
	268	icmax = ic;
	269	}
	270	ic++;
	271	}
	272	} //track loop
	273
	274	//
	275	// Cluster
42dc9410	276	if (ic < fNMin) {
70f2ce9d	277	PostData(1, fHists);
	278	return;
	279	}
	280	//
19c5a36a	281	Double_t* rk0 = new Double_t[10];
	282
	283	AliKMeansResult* res = 0;
	284
	285	for (Int_t i = 0; i < fK; i++) {
	286	res = fA[i];
	287	AliKMeansClustering::SoftKMeans2(i+1, ic, phi, eta, res->GetMx(),res->GetMy(), res->GetSigma2(), res->GetRk());
	288	res->Sort(ic, phi, eta);
	289	Int_t j = (res->GetInd())[0];
	290	rk0[i] = (res->GetTarget())[j];
	291	}
	292
	293	Double_t* mPhi = 0;
	294	Double_t* mEta = 0;
	295	Double_t* sigma2 = 0;
	296	Float_t rmax = -1.;
	297	Int_t imax = 0;
	298
	299	for (Int_t i = 0; i < fK; i++) {
	300	if (rk0[i] > rmax) {
	301	rmax = rk0[i];
	302	imax = i;
	303	}
	304	}
	305
	306
	307	res = fA[imax];
	308	mPhi = res->GetMx();
	309	mEta = res->GetMy();
	310	sigma2 = res->GetSigma2();
	311	Int_t nk = imax + 1;
	312	Int_t im = (res->GetInd())[0];
	313	fHDensity->Fill(rmax / TMath::Pi());
	314	fHCSize->Fill(rmax * sigma2[im]);
	315	fHNCluster->Fill(Float_t(nk));
	316
42dc9410	317	Double_t dphic, detac;
42dc9410	318
19c5a36a	319	if (rmax > 0.) {
42dc9410	320	// Analysis
	321	//
	322	// Cluster Multiplicity
	323	Int_t mult[2] = {0, 0};
	324	//
19c5a36a	325	if (nk > 1) {
	326	dphic = DeltaPhi(mPhi[0], mPhi[1]);
	327	detac = TMath::Abs(mEta[0] - mEta[1]);
	328	fH2DPhiEtaC->Fill(dphic, detac);
	329	}
42dc9410	330	//
	331	// Random cluster position
	332	Int_t ir = Int_t(Float_t(ic) * gRandom->Rndm());
a51f90df	333
42dc9410	334	Double_t crPhi = phi[ir];
	335	Double_t crEta = eta[ir];
	336	//
	337	Double_t sumPt = 0;
	338	Double_t sumPtC = 0;
19c5a36a	339
42dc9410	340	for (Int_t i = 0; i < 1; i++) {
19c5a36a	341	fH1CEta->Fill(mEta[im]);
19c5a36a	342	fH1CPhi->Fill(mPhi[im]);
70f2ce9d	343	for (Int_t j = 0; j < ic; j++) {
19c5a36a	344	Double_t r = DeltaR(mPhi[im], mEta[im], phi[j], eta[j]);
	345	Double_t dphi = DeltaPhi(mPhi[im], phi[j]);
	346	Double_t deta = mEta[im] - eta[j];
42dc9410	347	Double_t rr = DeltaR(crPhi, crEta, phi[j], eta[j]);
19c5a36a	348
42dc9410	349	fH1DR->Fill(r);
	350	fH1DPhi->Fill(dphi);
	351	fH2DPhiEta->Fill(dphi, TMath::Abs(deta));
	352	if (j == icmax) fH2DPhiEtaL->Fill(dphi, TMath::Abs(deta));
	353
	354	if (r < 0.2) {
	355	fH1PtC2->Fill(pt[j]);
	356	}
19c5a36a	357	if (r < 0.3) {
	358	fH1PtC1->Fill(pt[j]);
	359	fHPtDensity->Fill(rmax/TMath::Pi(), pt[j]);
	360	}
	361	if (rr < 0.3) {
42dc9410	362	fH1PtR->Fill(pt[j]);
19c5a36a	363	}
42dc9410	364
19c5a36a	365	if (r < 0.4) {
	366	sumPtC += pt[j];
	367	mult[i]++;
	368	fH1PtC->Fill(pt[j]);
	369	}
	370	if (r > 0.7 && dphi < (TMath::Pi() - 0.3)) {
42dc9410	371	fH1Pt->Fill(pt[j]);
	372	}
	373
	374	if (r > 0.7 && r < 1.) {
	375	sumPt += pt[j];
	376	}
19c5a36a	377
42dc9410	378	if (dphi > (TMath::Pi() - 0.3)) {
	379	fH1PtAS->Fill(pt[j]);
	380	}
70f2ce9d	381	}
42dc9410	382	}
70f2ce9d	383
42dc9410	384	fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
	385	fH1SPt ->Fill(sumPt);
	386	fH1SPtC->Fill(sumPtC);
19c5a36a	387	}
42dc9410	388	//
19c5a36a	389	// Randomized phi
42dc9410	390	//
19c5a36a	391	for (Int_t i = 0; i < fK; i++) {
	392	res = fB[i];
	393	AliKMeansClustering::SoftKMeans2(i+1, ic, phiR, etaR, res->GetMx(),res->GetMy(), res->GetSigma2(), res->GetRk());
	394	res->Sort(ic, phiR, etaR);
	395	//res->Sort();
	396	Int_t j = (res->GetInd())[0];
	397	rk0[i] = (res->GetTarget())[j];
42dc9410	398	}
42dc9410	399
19c5a36a	400	rmax = -1.;
	401	imax = 0;
	402	for (Int_t i = 0; i < fK; i++) {
	403	if (rk0[i] > rmax) {
	404	rmax = rk0[i];
	405	imax = i;
	406	}
	407	}
	408	res = fB[imax];
	409	mPhi = res->GetMx();
	410	mEta = res->GetMy();
	411	sigma2 = res->GetSigma2();
	412	nk = imax + 1;
	413
	414	// fH1RespR->Fill(rk[ind[0]]/(rk[0]+rk[1]));
	415	//
	416	// Cluster Multiplicity
	417	if (rmax > 0.) {
	418	for (Int_t i = 0; i < 1; i++) {
	419	Int_t im = (res->GetInd())[i];
	420	fH1CEtaR->Fill(mEta[im]);
	421	}
	422
	423	for (Int_t i = 0; i < 1; i++) {
	424	Int_t im = (res->GetInd())[i];
	425	for (Int_t j = 0; j < ic; j++) {
	426	Double_t dphi = DeltaPhi(mPhi[im], phiR[j]);
	427	Double_t deta = mEta[im] - etaR[j];
	428	Double_t r = DeltaR(mPhi[im], mEta[im], phiR[j], etaR[j]);
	429	fH1DRR->Fill(r);
	430	fH2DPhiEtaR->Fill(dphi, TMath::Abs(deta));
	431	if (j == icmax) fH2DPhiEtaLR->Fill(dphi, TMath::Abs(deta));
	432	}
	433	}
	434	if (nk > 1) {
	435	dphic = DeltaPhi(mPhi[0], mPhi[1]);
	436	detac = TMath::Abs(mEta[0] - mEta[1]);
	437	fH2DPhiEtaCR->Fill(dphic, detac);
70f2ce9d	438	}
42dc9410	439	}
19c5a36a	440	PostData(1, fHists);
70f2ce9d	441	}
	442
	443	Double_t AliAnalysisTaskKMeans::DeltaPhi(Double_t phi1, Double_t phi2)
	444	{
	445	Double_t dphi = TMath::Abs(phi1 - phi2);
	446	if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
	447	return dphi;
	448	}
	449
	450	Double_t AliAnalysisTaskKMeans::DeltaR(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
	451	{
	452	Double_t dphi = DeltaPhi(phi1, phi2);
	453	Double_t deta = eta1 - eta2;
	454	return (TMath::Sqrt(dphi * dphi + deta * deta));
	455
	456	}
	457
	458	//________________________________________________________________________
	459	void AliAnalysisTaskKMeans::Terminate(Option_t *)
	460	{
	461	}
	462