JETAN/AliAnalysisTaskKMeans.cxx

   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$ */
  17 //-------------------------------------------------------------------------
  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
  24 //-------------------------------------------------------------------------
  25
  26
  27
  28 #include "TChain.h"
  29 #include "TTree.h"
  30 #include "TH1F.h"
  31 #include "TH2F.h"
  32 #include "TProfile.h"
  33
  34 #include "TList.h"
  35 #include "TParticle.h"
  36 #include "TParticlePDG.h"
  37 #include "TProfile.h"
  38 #include "TMath.h"
  39 #include "TRandom.h"
  40
  41 #include "AliAnalysisTask.h"
  42 #include "AliAnalysisManager.h"
  43
  44 #include "AliESDEvent.h"
  45 #include "AliStack.h"
  46 #include "AliESDVertex.h"
  47 #include "AliESDInputHandler.h"
  48 #include "AliESDtrackCuts.h"
  49 #include "AliMultiplicity.h"
  50
  51 #include "AliMCParticle.h"
  52 #include "AliMCEvent.h"
  53 #include "AliAnalysisTaskKMeans.h"
  54 #include "AliTrackReference.h"
  55 #include "AliStack.h"
  56 #include "AliHeader.h"
  57 #include "AliKMeansClustering.h"
  58
  59
  60
  61 ClassImp(AliAnalysisTaskKMeans)
  62
  63 AliAnalysisTaskKMeans::AliAnalysisTaskKMeans()
  64     : AliAnalysisTaskSE()
  65     ,fK(0)
  66     ,fHists(0)
  67     ,fH1CEta(0)
  68     ,fH1CPhi(0)
  69     ,fH1CEtaR(0)
  70     ,fH2N1N2(0)
  71     ,fH1Pt(0)
  72     ,fH1PtC(0)
  73     ,fH1PtC1(0)
  74     ,fH1PtC2(0)
  75     ,fH1SPt(0)
  76     ,fH1SPtC(0)
  77     ,fH1DPhi(0)
  78     ,fH1DR(0)
  79     ,fH1DRR(0)
  80     ,fH2DPhiEta(0)
  81     ,fH2DPhiEtaR(0)
  82     ,fH2DPhiEtaL(0)
  83     ,fH2DPhiEtaLR(0)
  84     ,fH2DPhiEtaC(0)
  85     ,fH2DPhiEtaCR(0)
  86     ,fH1Resp(0)
  87     ,fH1RespR(0)
  88     ,fCuts(0)
  89 {
  90   //
  91   // Constructor
  92   //
  93 }
  94
  95 //________________________________________________________________________
  96 AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name)
  97     : AliAnalysisTaskSE(name)
  98       ,fK(0)
  99       ,fHists(0)
 100       ,fH1CEta(0)
 101       ,fH1CPhi(0)
 102       ,fH1CEtaR(0)
 103       ,fH2N1N2(0)
 104       ,fH1Pt(0)
 105       ,fH1PtC(0)
 106       ,fH1PtC1(0)
 107       ,fH1PtC2(0)
 108       ,fH1SPt(0)
 109       ,fH1SPtC(0)
 110       ,fH1DPhi(0)
 111       ,fH1DR(0)
 112       ,fH1DRR(0)
 113       ,fH2DPhiEta(0)
 114       ,fH2DPhiEtaR(0)
 115       ,fH2DPhiEtaL(0)
 116       ,fH2DPhiEtaLR(0)
 117       ,fH2DPhiEtaC(0)
 118       ,fH2DPhiEtaCR(0)
 119       ,fH1Resp(0)
 120       ,fH1RespR(0)
 121       ,fCuts(0)
 122 {
 123   //
 124   // Constructor
 125   //
 126   DefineOutput(1,  TList::Class());
 127 }
 128
 129
 130 //________________________________________________________________________
 131 void AliAnalysisTaskKMeans::UserCreateOutputObjects()
 132 {
 133   // Create histograms
 134   // Called once
 135     fHists   = new TList();
 136     fH1CEta  = new TH1F("fH1CEta",   "eta distribution of clusters",        90, -1.5, 1.5);
 137     fH1CEtaR = new TH1F("fH1CEtaR",  "eta distribution of clusters",        90, -1.5, 1.5);
 138     fH1CPhi  = new TH1F("fH1CPhi",   "phi distribution of clusters",       157,  0.0, 2. * TMath::Pi());
 139     fH2N1N2  = new TH2F("fH2N1N2",   "multiplicity distribution",          50, 0., 50., 50, 0., 50.);
 140
 141     fH1Pt    = new TH1F("fH1Pt",     "pt distribution",50, 0., 10.);
 142     fH1PtC   = new TH1F("fH1PtC",    "pt distribution",50, 0., 10.);
 143     fH1PtC1  = new TH1F("fH1PtC1",   "pt distribution",50, 0., 10.);
 144     fH1PtC2  = new TH1F("fH1PtC2",   "pt distribution",50, 0., 10.);
 145
 146     fH1SPt    = new TH1F("fH1SPt",     "sum pt distribution",50, 0., 10.);
 147     fH1SPtC   = new TH1F("fH1SPtC",    "sum pt distribution",50, 0., 10.);
 148
 149     fH1DR         = new TH1F("fH1DR",     "dR distribution", 50, 0., 5.);
 150     fH1DPhi       = new TH1F("fH1DPhi",   "dPhi distribution", 31, 0., TMath::Pi());
 151     fH2DPhiEta    = new TH2F("fH2DPhiEta","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 152     fH2DPhiEtaR   = new TH2F("fH2DPhiEtaR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 153     fH2DPhiEtaL   = new TH2F("fH2DPhiEtaL","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 154     fH2DPhiEtaLR  = new TH2F("fH2DPhiEtaLR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 155     fH2DPhiEtaC   = new TH2F("fH2DPhiEtaC","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 156     fH2DPhiEtaCR  = new TH2F("fH2DPhiEtaCR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 157     fH1DRR        = new TH1F("fH1DRR",    "dR distribution", 50, 0., 5.);
 158     fH1Resp       = new TH1F("fH1Resp",   "Responsibility", 50, 0., 1.);
 159     fH1RespR      = new TH1F("fH1RespR",  "Responsibility", 50, 0., 1.);
 160     fHists->SetOwner();
 161
 162     fHists->Add(fH1CEta);
 163     fHists->Add(fH1CEtaR);
 164     fHists->Add(fH1CPhi);
 165     fHists->Add(fH2N1N2);
 166     fHists->Add(fH1Pt);
 167     fHists->Add(fH1PtC);
 168     fHists->Add(fH1PtC1);
 169     fHists->Add(fH1PtC2);
 170     fHists->Add(fH1DR);
 171     fHists->Add(fH1SPtC);
 172     fHists->Add(fH1SPt);
 173     fHists->Add(fH1DPhi);
 174     fHists->Add(fH2DPhiEta);
 175     fHists->Add(fH2DPhiEtaR);
 176     fHists->Add(fH2DPhiEtaL);
 177     fHists->Add(fH2DPhiEtaLR);
 178     fHists->Add(fH2DPhiEtaC);
 179     fHists->Add(fH2DPhiEtaCR);
 180     fHists->Add(fH1DRR);
 181     fHists->Add(fH1RespR);
 182     fHists->Add(fH1Resp);
 183     //
 184     AliKMeansClustering::SetBeta(4.);
 185
 186 }
 187
 188 //________________________________________________________________________
 189 void AliAnalysisTaskKMeans::UserExec(Option_t *)
 190 {
 191   // Main loop
 192   // Called for each event
 193     Double_t   phi [500];
 194     Double_t   phiR[500];
 195     Double_t   eta[500];
 196     Double_t   etaR[500];
 197     Double_t   pt [500];
 198     Double_t   mPhi[20];
 199     Double_t   mEta[20];
 200     Double_t     rk[20];
 201     Int_t       ind[20];
 202
 203   if (!fInputEvent) {
 204     Printf("ERROR: fESD not available");
 205     return;
 206   }
 207
 208   AliESDEvent* esdE = (AliESDEvent*) fInputEvent;
 209
 210   Int_t ic = 0;
 211   //
 212   // Fill eta-phi positions
 213   //
 214   /*
 215   const AliMultiplicity *spdMult = esdE->GetMultiplicity();
 216   for (Int_t i = 0; i < spdMult->GetNumberOfTracklets(); ++i) {
 217           phi[ic] = spdMult->GetPhi(i);
 218           eta[ic] = spdMult->GetEta(i);
 219           ic++;
 220   }
 221   */
 222   Double_t ptmax = 0.;
 223   Int_t    icmax = -1;
 224
 225   for (Int_t iTracks = 0; iTracks < esdE->GetNumberOfTracks(); iTracks++) {
 226       AliESDtrack* track = esdE->GetTrack(iTracks);
 227       if ((fCuts->AcceptTrack(track)))
 228       {
 229           phi[ic]  = track->Phi();
 230           phiR[ic] = 2. * TMath::Pi() * gRandom->Rndm();
 231           eta[ic]  = track->Eta();
 232           etaR[ic] = 2. * gRandom->Rndm() - 1.;
 233           pt[ic]   = track->Pt();
 234           if (pt[ic] > ptmax) {
 235               ptmax = pt[ic];
 236               icmax = ic;
 237           }
 238           ic++;
 239       }
 240   } //track loop
 241
 242   //
 243   // Cluster
 244   if (ic < 10) {
 245       PostData(1, fHists);
 246       return;
 247   }
 248   //
 249   AliKMeansClustering::SoftKMeans(fK, ic, phi, eta, mPhi, mEta, rk);
 250   //
 251   // Sort
 252   TMath::Sort(fK, rk, ind);
 253   fH1Resp->Fill(rk[ind[0]]/(rk[0]+rk[1]));
 254   //
 255   // Analyse
 256   //
 257   // Cluster Multiplicity
 258   Int_t mult[2] = {0, 0};
 259   //
 260   Double_t dphic = DeltaPhi(mPhi[0], mPhi[1]);
 261   Double_t detac = TMath::Abs(mEta[0] - mEta[1]);
 262   fH2DPhiEtaC->Fill(dphic, detac);
 263
 264   //
 265   Double_t sumPt  = 0;
 266   Double_t sumPtC = 0;
 267   for (Int_t i = 0; i < 1; i++) {
 268       fH1CEta->Fill(mEta[ind[i]]);
 269       fH1CPhi->Fill(mPhi[ind[i]]);
 270       for (Int_t j = 0; j < ic; j++) {
 271           Double_t r    = DeltaR(mPhi[ind[i]], mEta[ind[i]], phi[j], eta[j]);
 272           Double_t dphi = DeltaPhi(mPhi[ind[i]], phi[j]);
 273           Double_t deta = mEta[ind[i]] - eta[j];
 274
 275           fH1DR->Fill(r);
 276           fH1DPhi->Fill(dphi);
 277           fH2DPhiEta->Fill(dphi, TMath::Abs(deta));
 278           if (j == icmax) fH2DPhiEtaL->Fill(dphi, TMath::Abs(deta));
 279
 280           if (r < 0.2) {
 281               fH1PtC2->Fill(pt[j]);
 282           }
 283           if (r < 0.3)
 284           {
 285               fH1PtC1->Fill(pt[j]);
 286           }
 287           if (r < 0.4)
 288           {
 289             sumPtC += pt[j];
 290             mult[i]++;
 291             fH1PtC->Fill(pt[j]);
 292           }
 293           if (r > 0.7) {
 294               fH1Pt->Fill(pt[j]);
 295           }
 296
 297           if (r > 0.7 && r < 1.) {
 298             sumPt += pt[j];
 299           }
 300       }
 301   }
 302
 303   fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
 304   fH1SPt ->Fill(sumPt);
 305   fH1SPtC->Fill(sumPtC);
 306
 307
 308   // Randomized phi
 309   //
 310   AliKMeansClustering::SoftKMeans(fK, ic, phiR, etaR, mPhi, mEta, rk);
 311   //
 312   // Sort
 313   TMath::Sort(fK, rk, ind);
 314   fH1RespR->Fill(rk[ind[0]]/(rk[0]+rk[1]));
 315   //
 316   // Analyse
 317   //
 318   // Cluster Multiplicity
 319   for (Int_t i = 0; i < 1; i++) {
 320       fH1CEtaR->Fill(mEta[ind[i]]);
 321   }
 322
 323   for (Int_t i = 0; i < 1; i++) {
 324       for (Int_t j = 0; j < ic; j++) {
 325           Double_t dphi = DeltaPhi(mPhi[ind[i]], phiR[j]);
 326           Double_t deta = mEta[ind[i]] - etaR[j];
 327           Double_t r = DeltaR(mPhi[ind[i]], mEta[ind[i]], phiR[j], etaR[j]);
 328           fH1DRR->Fill(r);
 329           fH2DPhiEtaR->Fill(dphi, TMath::Abs(deta));
 330           if (j == icmax) fH2DPhiEtaLR->Fill(dphi, TMath::Abs(deta));
 331       }
 332   }
 333   dphic = DeltaPhi(mPhi[0], mPhi[1]);
 334   detac = TMath::Abs(mEta[0] - mEta[1]);
 335   fH2DPhiEtaCR->Fill(dphic, detac);
 336
 337
 338   //
 339   // Post output data.
 340   PostData(1, fHists);
 341 }
 342
 343 Double_t AliAnalysisTaskKMeans::DeltaPhi(Double_t phi1, Double_t phi2)
 344 {
 345     Double_t dphi = TMath::Abs(phi1 - phi2);
 346     if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
 347     return dphi;
 348 }
 349
 350 Double_t AliAnalysisTaskKMeans::DeltaR(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
 351 {
 352     Double_t dphi = DeltaPhi(phi1, phi2);
 353     Double_t deta = eta1 - eta2;
 354     return (TMath::Sqrt(dphi * dphi + deta * deta));
 355
 356 }
 357
 358 //________________________________________________________________________
 359 void AliAnalysisTaskKMeans::Terminate(Option_t *)
 360 {
 361 }
 362