#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"
ClassImp(AliAnalysisTaskKMeans)
AliAnalysisTaskKMeans::AliAnalysisTaskKMeans()
- : AliAnalysisTaskSE()
+ : AliAnalysisTaskSE()
+ ,fK(0)
+ ,fNMin(0)
,fHists(0)
,fH1CEta(0)
,fH1CPhi(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)
{
//
//________________________________________________________________________
AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name)
: AliAnalysisTaskSE(name)
+ ,fK(0)
+ ,fNMin(0)
,fHists(0)
,fH1CEta(0)
,fH1CPhi(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)
{
//
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.);
- 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.);
- fH1DRR = new TH1F("fH1DRR", "dR distribution", 50, 0., 5.);
-
+ 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(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);
//
- AliKMeansClustering::SetBeta(20.);
-
+ for (Int_t i = 0; i < 10; i++) {
+ fA[i] = new AliKMeansResult(i+1);
+ fB[i] = new AliKMeansResult(i+1);
+ }
}
//________________________________________________________________________
// Called for each event
Double_t phi [500];
Double_t phiR[500];
- Double_t eta[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];
- Int_t ind[20];
+ Double_t pt [500];
if (!fInputEvent) {
Printf("ERROR: fESD not available");
return;
}
- AliESDEvent* esdE = (AliESDEvent*) fInputEvent;
Int_t ic = 0;
+
//
// Fill eta-phi positions
//
- /*
- const AliMultiplicity *spdMult = esdE->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 < esdE->GetNumberOfTracks(); iTracks++) {
- AliESDtrack* track = esdE->GetTrack(iTracks);
- if ((fCuts->AcceptTrack(track)))
+ for (Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++) {
+ AliVParticle* track = fInputEvent->GetTrack(iTracks);
+ if ((fCuts->AcceptTrack((AliESDtrack*)track)))
{
- phi[ic] = track->Phi();
+ 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();
+ 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;
//
// Cluster
- if (ic < 10) {
+ if (ic < fNMin) {
PostData(1, fHists);
return;
}
//
- AliKMeansClustering::SoftKMeans(2, ic, phi, eta, mPhi, mEta, rk);
- //
- // Sort
- TMath::Sort(2, rk, ind);
- //
- // Analyse
- //
- // Cluster Multiplicity
- Int_t mult[2] = {0, 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];
-
- 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 (r < 0.4)
- {
- mult[i]++;
- fH1PtC->Fill(pt[j]);
- }
- if (r > 0.7) {
- fH1Pt->Fill(pt[j]);
- }
+ 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;
}
}
- fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
+ 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;
- // Randomized phi
- //
- AliKMeansClustering::SoftKMeans(2, ic, phiR, etaR, mPhi, mEta, rk);
- //
- // Sort
- TMath::Sort(2, rk, ind);
- //
- // Analyse
- //
- // Cluster Multiplicity
- for (Int_t i = 0; i < 1; i++) {
- fH1CEtaR->Fill(mEta[ind[i]]);
- }
- for (Int_t i = 0; i < 2; i++) {
+ 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 dphi = DeltaPhi(mPhi[ind[i]], phi[j]);
- Double_t deta = mEta[ind[i]] - eta[j];
- Double_t r = DeltaR(mPhi[ind[i]], mEta[ind[i]], phi[j], eta[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));
+ }
}
- }
-
-
- //
- // Post output data.
- PostData(1, fHists);
+ 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)
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff