[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
  //
}

//________________________________________________________________________
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
  //
  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
}

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