[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
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
e42693ad	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
70f2ce9d	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.);
42dc9410	207	fH1PtAS = new TH1F("fH1PtAS", "pt distribution",50, 0., 10.);
42dc9410	208	fH1PtR = new TH1F("fH1PtR", "pt distribution",50, 0., 10.);
70f2ce9d	209
a51f90df	210	fH1SPt = new TH1F("fH1SPt", "sum pt distribution",50, 0., 10.);
	211	fH1SPtC = new TH1F("fH1SPtC", "sum pt distribution",50, 0., 10.);
	212
5a9b5891	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.);
320bb308	222	fH1Resp = new TH1F("fH1Resp", "Responsibility", 50, 0., 1.);
320bb308	223	fH1RespR = new TH1F("fH1RespR", "Responsibility", 50, 0., 1.);
42dc9410	224	fH2Sigma = new TH2F("fH2Sigma", "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.);
42dc9410	225	fH2SigmaR = new TH2F("fH2SigmaR", "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.);
19c5a36a	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
42dc9410	234
70f2ce9d	235	fHists->SetOwner();
	236
	237	fHists->Add(fH1CEta);
	238	fHists->Add(fH1CEtaR);
	239	fHists->Add(fH1CPhi);
	240	fHists->Add(fH2N1N2);
	241	fHists->Add(fH1Pt);
42dc9410	242	fHists->Add(fH1PtR);
70f2ce9d	243	fHists->Add(fH1PtC);
	244	fHists->Add(fH1PtC1);
	245	fHists->Add(fH1PtC2);
42dc9410	246	fHists->Add(fH1PtAS);
70f2ce9d	247	fHists->Add(fH1DR);
a51f90df	248	fHists->Add(fH1SPtC);
a51f90df	249	fHists->Add(fH1SPt);
70f2ce9d	250	fHists->Add(fH1DPhi);
	251	fHists->Add(fH2DPhiEta);
	252	fHists->Add(fH2DPhiEtaR);
	253	fHists->Add(fH2DPhiEtaL);
5a9b5891	254	fHists->Add(fH2DPhiEtaLR);
a51f90df	255	fHists->Add(fH2DPhiEtaC);
a51f90df	256	fHists->Add(fH2DPhiEtaCR);
70f2ce9d	257	fHists->Add(fH1DRR);
320bb308	258	fHists->Add(fH1RespR);
320bb308	259	fHists->Add(fH1Resp);
42dc9410	260	fHists->Add(fH2Sigma);
19c5a36a	261	fHists->Add(fH2SigmaR);
	262	fHists->Add(fHCSize);
	263	fHists->Add(fHDensity);
	264	fHists->Add(fHNCluster);
	265	fHists->Add(fHPtDensity);
70f2ce9d	266	//
19c5a36a	267	for (Int_t i = 0; i < 10; i++) {
	268	fA[i] = new AliKMeansResult(i+1);
	269	fB[i] = new AliKMeansResult(i+1);
	270	}
70f2ce9d	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];
19c5a36a	280	Double_t eta[500];
70f2ce9d	281	Double_t etaR[500];
19c5a36a	282	Double_t pt [500];
70f2ce9d	283
	284	if (!fInputEvent) {
	285	Printf("ERROR: fESD not available");
	286	return;
	287	}
	288
70f2ce9d	289
70f2ce9d	290	Int_t ic = 0;
19c5a36a	291
70f2ce9d	292	//
	293	// Fill eta-phi positions
	294	//
19c5a36a	295
70f2ce9d	296	Double_t ptmax = 0.;
	297	Int_t icmax = -1;
	298
42dc9410	299	for (Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++) {
	300	AliVParticle* track = fInputEvent->GetTrack(iTracks);
	301	if ((fCuts->AcceptTrack((AliESDtrack*)track)))
70f2ce9d	302	{
19c5a36a	303	phi [ic] = track->Phi();
70f2ce9d	304	phiR[ic] = 2. * TMath::Pi() * gRandom->Rndm();
19c5a36a	305	eta [ic] = track->Eta();
	306	etaR[ic] = 1.8 * gRandom->Rndm() - 0.9;
	307	pt [ic] = track->Pt();
70f2ce9d	308	if (pt[ic] > ptmax) {
	309	ptmax = pt[ic];
	310	icmax = ic;
	311	}
	312	ic++;
	313	}
	314	} //track loop
	315
	316	//
	317	// Cluster
42dc9410	318	if (ic < fNMin) {
70f2ce9d	319	PostData(1, fHists);
	320	return;
	321	}
	322	//
19c5a36a	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
42dc9410	359	Double_t dphic, detac;
42dc9410	360
19c5a36a	361	if (rmax > 0.) {
42dc9410	362	// Analysis
	363	//
	364	// Cluster Multiplicity
	365	Int_t mult[2] = {0, 0};
	366	//
19c5a36a	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	}
42dc9410	372	//
	373	// Random cluster position
	374	Int_t ir = Int_t(Float_t(ic) * gRandom->Rndm());
a51f90df	375
42dc9410	376	Double_t crPhi = phi[ir];
	377	Double_t crEta = eta[ir];
	378	//
	379	Double_t sumPt = 0;
	380	Double_t sumPtC = 0;
19c5a36a	381
42dc9410	382	for (Int_t i = 0; i < 1; i++) {
19c5a36a	383	fH1CEta->Fill(mEta[im]);
19c5a36a	384	fH1CPhi->Fill(mPhi[im]);
70f2ce9d	385	for (Int_t j = 0; j < ic; j++) {
19c5a36a	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];
42dc9410	389	Double_t rr = DeltaR(crPhi, crEta, phi[j], eta[j]);
19c5a36a	390
42dc9410	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	}
19c5a36a	399	if (r < 0.3) {
	400	fH1PtC1->Fill(pt[j]);
	401	fHPtDensity->Fill(rmax/TMath::Pi(), pt[j]);
	402	}
	403	if (rr < 0.3) {
42dc9410	404	fH1PtR->Fill(pt[j]);
19c5a36a	405	}
42dc9410	406
19c5a36a	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)) {
42dc9410	413	fH1Pt->Fill(pt[j]);
	414	}
	415
	416	if (r > 0.7 && r < 1.) {
	417	sumPt += pt[j];
	418	}
19c5a36a	419
42dc9410	420	if (dphi > (TMath::Pi() - 0.3)) {
	421	fH1PtAS->Fill(pt[j]);
	422	}
70f2ce9d	423	}
42dc9410	424	}
70f2ce9d	425
42dc9410	426	fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
	427	fH1SPt ->Fill(sumPt);
	428	fH1SPtC->Fill(sumPtC);
19c5a36a	429	}
42dc9410	430	//
19c5a36a	431	// Randomized phi
42dc9410	432	//
19c5a36a	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];
42dc9410	440	}
42dc9410	441
19c5a36a	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);
70f2ce9d	480	}
42dc9410	481	}
19c5a36a	482	PostData(1, fHists);
70f2ce9d	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