[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 "TFile.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 "AliExternalTrackParam.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)
    ,fHDPhi(0)
    ,fH2EtaPhi(0)
    ,fH2REtaPhi(0)
    ,fCuts(0)

{
  //

  // Constructor
  //
    for (Int_t i=0; i< 10; i++) {
	fA[i] = 0;
	fB[i] = 0;
    }
}

//________________________________________________________________________
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)
      ,fHDPhi(0)
      ,fH2EtaPhi(0)
      ,fH2REtaPhi(0)
      ,fCuts(0)

{
  //
  // Constructor
  //
    for (Int_t i=0; i< 10; i++) {
	fA[i] = 0;
	fB[i] = 0;
    }
  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)
      ,fHDensity(0)
      ,fHCSize(0)  
      ,fHNCluster(0)
      ,fHPtDensity(0)
      ,fHDPhi(0)
      ,fH2EtaPhi(0)
      ,fH2REtaPhi(0)
      ,fCuts(0)
{
    // Dummy copy constructor
    for (Int_t i=0; i< 10; i++) {
	fA[i] = 0;
	fB[i] = 0;
    }
}

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.);
    fHDPhi       = new TH1F("fHDPhi",   "phi correlation", 100, 0., 0.5);
    fH2EtaPhi    = new TH2F("fH2EtaPhi", "Eta Phi", 200, -1., 1., 628, 0., 2. * TMath::Pi());  
    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);
    fHists->Add(fHDPhi);
    fHists->Add(fH2EtaPhi);
    //
    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] = {0};
    Double_t   phiR[500] = {0};
    Double_t    eta[500] = {0};
    Double_t   etaR[500] = {0};
    Double_t    pt [500] = {0};

  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))) 
      {
	const AliExternalTrackParam * tpcT = ((AliESDtrack*) track)->GetTPCInnerParam();
	if (!tpcT) continue;
	if (TMath::Abs(tpcT->Eta()) > 0.9) continue;

	  phi [ic] = tpcT->Phi();
	  eta [ic] = tpcT->Eta();
	  pt  [ic] = tpcT->Pt();

	  if (fH2REtaPhi) {
	    fH2REtaPhi->GetRandom2(etaR[ic], phiR[ic]);
	  } else {
	    phiR[ic] = 2. * TMath::Pi() * gRandom->Rndm();
	    etaR[ic] = 1.8 * gRandom->Rndm() - 0.9;
	  }
	  

	  if (pt[ic] > ptmax) {
	      ptmax = pt[ic];
	      icmax = ic;
	  }

	  fH2EtaPhi->Fill(eta[ic], phi[ic]);
	  ic++;
      }
  } //track loop 

  for (Int_t i = 0; i < ic; i++) {
    for (Int_t j = i+1; j < ic; j++) {
      Double_t dphi = TMath::Abs(phi[i] - phi[j]);
      fHDPhi->Fill(dphi);
    }
  }

  //
  // Cluster
  if (ic < fNMin) {
      PostData(1, fHists);
      return;
  }

  //
  Double_t rk0[10];
  AliKMeansResult* res = 0;
  AliKMeansResult best(10);
  Float_t   rmaxG = -1.;

  for (Int_t k = 0; k < 20; k++) {
    Float_t   rmax   = -1.;
    Int_t     imax   = 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];
      if (rk0[i] > rmax)  {
	rmax = rk0[i];
	imax = i;
      }
    }
    if (rmax > rmaxG) {
      rmaxG = rmax;
      best.CopyResults(fA[imax]);
    }
  }
  Double_t* mPhi     = best.GetMx();
  Double_t* mEta     = best.GetMy();
  Double_t* sigma2   = best.GetSigma2();
  Int_t     nk       = best.GetK();
  Int_t     im       = (best.GetInd())[0];
  Double_t  etaC     = mEta[im];

  fHDensity->Fill(rmaxG / TMath::Pi());
  fHCSize->Fill(2.28 * rmaxG * sigma2[im]);
  fHNCluster->Fill(Float_t(nk));

  Double_t dphic, detac;

  if (rmaxG > 0. && TMath::Abs(etaC) < 0.4) {
    // 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(rmaxG/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], 1.);
	}
      }
    }

    fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
    fH1SPt ->Fill(sumPt);
    fH1SPtC->Fill(sumPtC);
  }
    //
    // Randomized phi
    //
  rmaxG = -1.;
  for (Int_t k = 0; k < 20; k++) {
    Float_t rmax   = -1.;
    Int_t   imax   =  0;
    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);
      Int_t j = (res->GetInd())[0];
      rk0[i]  = (res->GetTarget())[j];
      if (rk0[i] > rmax) {
	rmax = rk0[i];
	imax = i;
      }
    }
    if (rmax > rmaxG) {
      rmaxG = rmax;
      best.CopyResults(fB[imax]);
    }    
  }
  
    mPhi    = best.GetMx();
    mEta    = best.GetMy();
    nk      = best.GetK();
    im      = (best.GetInd())[0];
    etaC    = mEta[im];    

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