/************************************************************************** * 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 *) { }