Update Balance Function correction Maps (Alis Rodriguez Manso <alisrm@nikhef.nl>)
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskKMeans.cxx
index 24162f9..c47df6f 100644 (file)
  **************************************************************************/
 
 /* $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 "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"
@@ -61,7 +65,9 @@
 ClassImp(AliAnalysisTaskKMeans)
 
 AliAnalysisTaskKMeans::AliAnalysisTaskKMeans() 
-    : AliAnalysisTaskSE() 
+    : AliAnalysisTaskSE()
+    ,fK(0)
+    ,fNMin(0)
     ,fHists(0)
     ,fH1CEta(0)
     ,fH1CPhi(0)
@@ -71,6 +77,8 @@ AliAnalysisTaskKMeans::AliAnalysisTaskKMeans()
     ,fH1PtC(0)
     ,fH1PtC1(0)
     ,fH1PtC2(0)
+    ,fH1PtAS(0)
+    ,fH1PtR(0)
     ,fH1SPt(0)
     ,fH1SPtC(0)
     ,fH1DPhi(0)
@@ -82,16 +90,35 @@ AliAnalysisTaskKMeans::AliAnalysisTaskKMeans()
     ,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)
@@ -101,6 +128,8 @@ AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name)
       ,fH1PtC(0)
       ,fH1PtC1(0)
       ,fH1PtC2(0)
+      ,fH1PtAS(0)
+      ,fH1PtR(0)
       ,fH1SPt(0)
       ,fH1SPtC(0)
       ,fH1DPhi(0)
@@ -112,14 +141,83 @@ AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name)
       ,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()
@@ -136,6 +234,8 @@ 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.);
 
     fH1SPt    = new TH1F("fH1SPt",     "sum pt distribution",50, 0., 10.);
     fH1SPtC   = new TH1F("fH1SPtC",    "sum pt distribution",50, 0., 10.);
@@ -149,7 +249,19 @@ void AliAnalysisTaskKMeans::UserCreateOutputObjects()
     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);
@@ -157,9 +269,11 @@ 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);
@@ -171,10 +285,21 @@ void AliAnalysisTaskKMeans::UserCreateOutputObjects()
     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);
     //
-    AliKMeansClustering::SetBeta(20.);
-    
+    for (Int_t i = 0; i < 10; i++) {
+      fA[i] = new AliKMeansResult(i+1);
+      fB[i] = new AliKMeansResult(i+1);
+    }
 }
 
 //________________________________________________________________________
@@ -182,152 +307,234 @@ 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];
-    Double_t   mPhi[20];
-    Double_t   mEta[20];
-    Double_t     rk[20];
-    Int_t       ind[20];
+    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;
   }
   
-  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();
-         phiR[ic] = 2. * TMath::Pi() * gRandom->Rndm();
-         eta[ic]  = track->Eta();
-         etaR[ic] = 2. * gRandom->Rndm() - 1.;
-         pt[ic]   = track->Pt();
+       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 < 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};
-  //
-  Double_t dphic = DeltaPhi(mPhi[0], mPhi[1]);
-  Double_t detac = TMath::Abs(mEta[0] - mEta[1]);
-  fH2DPhiEtaC->Fill(dphic, detac);  
 
   //
-  Double_t sumPt  = 0;
-  Double_t sumPtC = 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) 
-         {
-           sumPtC += pt[j];
-           mult[i]++;
-           fH1PtC->Fill(pt[j]);
-         } 
-         if (r > 0.7) {
-             fH1Pt->Fill(pt[j]);
-         }
-
-         if (r > 0.7 && r < 1.) {
-           sumPt += pt[j];
-         }
+  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];
 
-  fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
-  fH1SPt ->Fill(sumPt);
-  fH1SPtC->Fill(sumPtC);
-  
+  fHDensity->Fill(rmaxG / TMath::Pi());
+  fHCSize->Fill(2.28 * rmaxG * sigma2[im]);
+  fHNCluster->Fill(Float_t(nk));
 
-  // 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]]);
-  }
+  Double_t dphic, detac;
 
-  for (Int_t i = 0; i < 1; i++) {
+  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 dphi = DeltaPhi(mPhi[ind[i]], phiR[j]);
-         Double_t deta = mEta[ind[i]] - etaR[j];
-         Double_t r = DeltaR(mPhi[ind[i]], mEta[ind[i]], phiR[j], etaR[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));
+       }
       }
-  }
-  dphic = DeltaPhi(mPhi[0], mPhi[1]);
-  detac = TMath::Abs(mEta[0] - mEta[1]);
-  fH2DPhiEtaCR->Fill(dphic, detac);  
-
-  
-  //
-  // 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)