[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());
}

AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const AliAnalysisTaskKMeans &res)
: AliAnalysisTaskSE(res) 
      ,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)
      ,fCuts(0)
{
    // Dummy copy constructor
}

AliAnalysisTaskKMeans& AliAnalysisTaskKMeans::operator=(const AliAnalysisTaskKMeans& /*trclass*/)
{
    // Dummy assignment operator
    return *this;
}


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