[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 "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)
    ,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)
      ,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.);
 
    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);    
    
    //
    
}

//________________________________________________________________________
void AliAnalysisTaskKMeans::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event
    Double_t   phi [500];
    Double_t   phiR[500];
    Double_t   eta[500];
    Double_t   etaR[500];
    Double_t   pt [500];
    Double_t   mPhi[20];
    Double_t   mEta[20];
    Double_t     rk[20];
    Double_t sigma2[20];
    Double_t sigmax2[20];
    Double_t sigmay2[20];
    Int_t       ind[20];

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

  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;
  
  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] = 2. * gRandom->Rndm() - 1.;
	  pt[ic]   = track->Pt();
	  if (pt[ic] > ptmax) {
	      ptmax = pt[ic];
	      icmax = ic;
	  }
	  ic++;
      }
  } //track loop 

  //
  // Cluster
  if (ic < fNMin) {
      PostData(1, fHists);
      return;
  }
  //
  Int_t ok;
  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]])); 
    //
    // 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);  
    //
    // 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[ind[i]]);
      fH1CPhi->Fill(mPhi[ind[i]]);      
      for (Int_t j = 0; j < ic; j++) {
	Double_t r    = DeltaR(mPhi[ind[i]], mEta[ind[i]], phi[j], eta[j]);
	Double_t dphi = DeltaPhi(mPhi[ind[i]], phi[j]);
	Double_t deta = mEta[ind[i]] - eta[j];
	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]);
	  }
	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) {
	  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);
  } // 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
    //
    // Cluster Multiplicity
    for (Int_t i = 0; i < 1; i++) {
      fH1CEtaR->Fill(mEta[ind[i]]);
    }
    
    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));
      }
    }
    dphic = DeltaPhi(mPhi[0], mPhi[1]);
    detac = TMath::Abs(mEta[0] - mEta[1]);
    fH2DPhiEtaCR->Fill(dphic, detac);  
  } // ok
  
  //
  // Post output data.
  PostData(1, fHists);
}      

Double_t AliAnalysisTaskKMeans::DeltaPhi(Double_t phi1, Double_t phi2)
{
    Double_t dphi = TMath::Abs(phi1 - phi2);
    if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
    return dphi;
}

Double_t AliAnalysisTaskKMeans::DeltaR(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
{
    Double_t dphi = DeltaPhi(phi1, phi2);
    Double_t deta = eta1 - eta2;
    return (TMath::Sqrt(dphi * dphi + deta * deta));
    
}

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