X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliAnalysisTaskKMeans.cxx;h=191aba1974eeebe6be3bb786b4b05a71e68e04ef;hb=b06a50a5f78aa1afba6e6c3f5e569754ded23417;hp=87b9897046d755e0167bbcf08f2fba9304422e3f;hpb=70f2ce9da81205a0ae2c5f86e67a4c2e41ec5516;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliAnalysisTaskKMeans.cxx b/JETAN/AliAnalysisTaskKMeans.cxx index 87b9897046d..191aba1974e 100644 --- a/JETAN/AliAnalysisTaskKMeans.cxx +++ b/JETAN/AliAnalysisTaskKMeans.cxx @@ -37,10 +37,12 @@ #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" @@ -61,7 +63,9 @@ ClassImp(AliAnalysisTaskKMeans) AliAnalysisTaskKMeans::AliAnalysisTaskKMeans() - : AliAnalysisTaskSE() + : AliAnalysisTaskSE() + ,fK(0) + ,fNMin(0) ,fHists(0) ,fH1CEta(0) ,fH1CPhi(0) @@ -71,12 +75,27 @@ AliAnalysisTaskKMeans::AliAnalysisTaskKMeans() ,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) { // @@ -87,6 +106,8 @@ AliAnalysisTaskKMeans::AliAnalysisTaskKMeans() //________________________________________________________________________ AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name) : AliAnalysisTaskSE(name) + ,fK(0) + ,fNMin(0) ,fHists(0) ,fH1CEta(0) ,fH1CPhi(0) @@ -96,12 +117,27 @@ AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name) ,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) { // @@ -126,14 +162,34 @@ void AliAnalysisTaskKMeans::UserCreateOutputObjects() 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); @@ -141,19 +197,35 @@ void AliAnalysisTaskKMeans::UserCreateOutputObjects() 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); + } } //________________________________________________________________________ @@ -163,45 +235,34 @@ void AliAnalysisTaskKMeans::UserExec(Option_t *) // 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; @@ -212,83 +273,171 @@ void AliAnalysisTaskKMeans::UserExec(Option_t *) // // 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)