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 #include "TObjArray.h"
  41
  42 #include "AliAnalysisTask.h"
  43 #include "AliAnalysisManager.h"
  44
  45 #include "AliVEvent.h"
  46 #include "AliESDEvent.h"
  47 #include "AliStack.h"
  48 #include "AliESDVertex.h"
  49 #include "AliESDInputHandler.h"
  50 #include "AliESDtrackCuts.h"
  51 #include "AliMultiplicity.h"
  52
  53 #include "AliMCParticle.h"
  54 #include "AliMCEvent.h"
  55 #include "AliAnalysisTaskKMeans.h"
  56 #include "AliTrackReference.h"
  57 #include "AliStack.h"
  58 #include "AliHeader.h"
  59 #include "AliKMeansClustering.h"
  60
  61
  62
  63 ClassImp(AliAnalysisTaskKMeans)
  64
  65 AliAnalysisTaskKMeans::AliAnalysisTaskKMeans()
  66     : AliAnalysisTaskSE()
  67     ,fK(0)
  68     ,fNMin(0)
  69     ,fHists(0)
  70     ,fH1CEta(0)
  71     ,fH1CPhi(0)
  72     ,fH1CEtaR(0)
  73     ,fH2N1N2(0)
  74     ,fH1Pt(0)
  75     ,fH1PtC(0)
  76     ,fH1PtC1(0)
  77     ,fH1PtC2(0)
  78     ,fH1PtAS(0)
  79     ,fH1PtR(0)
  80     ,fH1SPt(0)
  81     ,fH1SPtC(0)
  82     ,fH1DPhi(0)
  83     ,fH1DR(0)
  84     ,fH1DRR(0)
  85     ,fH2DPhiEta(0)
  86     ,fH2DPhiEtaR(0)
  87     ,fH2DPhiEtaL(0)
  88     ,fH2DPhiEtaLR(0)
  89     ,fH2DPhiEtaC(0)
  90     ,fH2DPhiEtaCR(0)
  91     ,fH1Resp(0)
  92     ,fH1RespR(0)
  93     ,fH2Sigma(0)
  94     ,fH2SigmaR(0)
  95     ,fHDensity(0)
  96     ,fHCSize(0)
  97     ,fHNCluster(0)
  98     ,fHPtDensity(0)
  99     ,fCuts(0)
 100 {
 101   //
 102   // Constructor
 103   //
 104 }
 105
 106 //________________________________________________________________________
 107 AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name)
 108     : AliAnalysisTaskSE(name)
 109       ,fK(0)
 110       ,fNMin(0)
 111       ,fHists(0)
 112       ,fH1CEta(0)
 113       ,fH1CPhi(0)
 114       ,fH1CEtaR(0)
 115       ,fH2N1N2(0)
 116       ,fH1Pt(0)
 117       ,fH1PtC(0)
 118       ,fH1PtC1(0)
 119       ,fH1PtC2(0)
 120       ,fH1PtAS(0)
 121       ,fH1PtR(0)
 122       ,fH1SPt(0)
 123       ,fH1SPtC(0)
 124       ,fH1DPhi(0)
 125       ,fH1DR(0)
 126       ,fH1DRR(0)
 127       ,fH2DPhiEta(0)
 128       ,fH2DPhiEtaR(0)
 129       ,fH2DPhiEtaL(0)
 130       ,fH2DPhiEtaLR(0)
 131       ,fH2DPhiEtaC(0)
 132       ,fH2DPhiEtaCR(0)
 133       ,fH1Resp(0)
 134       ,fH1RespR(0)
 135       ,fH2Sigma(0)
 136       ,fH2SigmaR(0)
 137       ,fHDensity(0)
 138       ,fHCSize(0)
 139       ,fHNCluster(0)
 140       ,fHPtDensity(0)
 141       ,fCuts(0)
 142 {
 143   //
 144   // Constructor
 145   //
 146   DefineOutput(1,  TList::Class());
 147 }
 148
 149
 150 //________________________________________________________________________
 151 void AliAnalysisTaskKMeans::UserCreateOutputObjects()
 152 {
 153   // Create histograms
 154   // Called once
 155     fHists   = new TList();
 156     fH1CEta  = new TH1F("fH1CEta",   "eta distribution of clusters",        90, -1.5, 1.5);
 157     fH1CEtaR = new TH1F("fH1CEtaR",  "eta distribution of clusters",        90, -1.5, 1.5);
 158     fH1CPhi  = new TH1F("fH1CPhi",   "phi distribution of clusters",       157,  0.0, 2. * TMath::Pi());
 159     fH2N1N2  = new TH2F("fH2N1N2",   "multiplicity distribution",          50, 0., 50., 50, 0., 50.);
 160
 161     fH1Pt    = new TH1F("fH1Pt",     "pt distribution",50, 0., 10.);
 162     fH1PtC   = new TH1F("fH1PtC",    "pt distribution",50, 0., 10.);
 163     fH1PtC1  = new TH1F("fH1PtC1",   "pt distribution",50, 0., 10.);
 164     fH1PtC2  = new TH1F("fH1PtC2",   "pt distribution",50, 0., 10.);
 165     fH1PtAS  = new TH1F("fH1PtAS",   "pt distribution",50, 0., 10.);
 166     fH1PtR   = new TH1F("fH1PtR",    "pt distribution",50, 0., 10.);
 167
 168     fH1SPt    = new TH1F("fH1SPt",     "sum pt distribution",50, 0., 10.);
 169     fH1SPtC   = new TH1F("fH1SPtC",    "sum pt distribution",50, 0., 10.);
 170
 171     fH1DR         = new TH1F("fH1DR",     "dR distribution", 50, 0., 5.);
 172     fH1DPhi       = new TH1F("fH1DPhi",   "dPhi distribution", 31, 0., TMath::Pi());
 173     fH2DPhiEta    = new TH2F("fH2DPhiEta","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 174     fH2DPhiEtaR   = new TH2F("fH2DPhiEtaR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 175     fH2DPhiEtaL   = new TH2F("fH2DPhiEtaL","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 176     fH2DPhiEtaLR  = new TH2F("fH2DPhiEtaLR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 177     fH2DPhiEtaC   = new TH2F("fH2DPhiEtaC","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 178     fH2DPhiEtaCR  = new TH2F("fH2DPhiEtaCR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 179     fH1DRR        = new TH1F("fH1DRR",    "dR distribution", 50, 0., 5.);
 180     fH1Resp       = new TH1F("fH1Resp",   "Responsibility", 50, 0., 1.);
 181     fH1RespR      = new TH1F("fH1RespR",  "Responsibility", 50, 0., 1.);
 182     fH2Sigma      = new TH2F("fH2Sigma",  "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.);
 183     fH2SigmaR     = new TH2F("fH2SigmaR", "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.);
 184
 185
 186     fHDensity    = new TH1F("fHDensity",    "density distribution", 100, 0., 20.);
 187     fHCSize      = new TH1F("fHCSize",      "cluster size", 20, -0.5, 19.5);
 188
 189     fHNCluster   = new TH1F("fHNCluster",   "Number of clusters", 11, -0.5, 10.5);
 190     fHPtDensity  = new TH2F("fHPtDensity", "Pt vs density", 100, 0., 20., 50, 0., 10.);
 191
 192
 193     fHists->SetOwner();
 194
 195     fHists->Add(fH1CEta);
 196     fHists->Add(fH1CEtaR);
 197     fHists->Add(fH1CPhi);
 198     fHists->Add(fH2N1N2);
 199     fHists->Add(fH1Pt);
 200     fHists->Add(fH1PtR);
 201     fHists->Add(fH1PtC);
 202     fHists->Add(fH1PtC1);
 203     fHists->Add(fH1PtC2);
 204     fHists->Add(fH1PtAS);
 205     fHists->Add(fH1DR);
 206     fHists->Add(fH1SPtC);
 207     fHists->Add(fH1SPt);
 208     fHists->Add(fH1DPhi);
 209     fHists->Add(fH2DPhiEta);
 210     fHists->Add(fH2DPhiEtaR);
 211     fHists->Add(fH2DPhiEtaL);
 212     fHists->Add(fH2DPhiEtaLR);
 213     fHists->Add(fH2DPhiEtaC);
 214     fHists->Add(fH2DPhiEtaCR);
 215     fHists->Add(fH1DRR);
 216     fHists->Add(fH1RespR);
 217     fHists->Add(fH1Resp);
 218     fHists->Add(fH2Sigma);
 219     fHists->Add(fH2SigmaR);
 220     fHists->Add(fHCSize);
 221     fHists->Add(fHDensity);
 222     fHists->Add(fHNCluster);
 223     fHists->Add(fHPtDensity);
 224     //
 225     for (Int_t i = 0; i < 10; i++) {
 226       fA[i] = new AliKMeansResult(i+1);
 227       fB[i] = new AliKMeansResult(i+1);
 228     }
 229 }
 230
 231 //________________________________________________________________________
 232 void AliAnalysisTaskKMeans::UserExec(Option_t *)
 233 {
 234   // Main loop
 235   // Called for each event
 236     Double_t   phi [500];
 237     Double_t   phiR[500];
 238     Double_t    eta[500];
 239     Double_t   etaR[500];
 240     Double_t    pt [500];
 241
 242   if (!fInputEvent) {
 243     Printf("ERROR: fESD not available");
 244     return;
 245   }
 246
 247
 248   Int_t ic = 0;
 249
 250   //
 251   // Fill eta-phi positions
 252   //
 253
 254   Double_t ptmax = 0.;
 255   Int_t    icmax = -1;
 256
 257   for (Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++) {
 258       AliVParticle* track = fInputEvent->GetTrack(iTracks);
 259       if ((fCuts->AcceptTrack((AliESDtrack*)track)))
 260       {
 261           phi [ic] = track->Phi();
 262           phiR[ic] = 2. * TMath::Pi() * gRandom->Rndm();
 263           eta [ic] = track->Eta();
 264           etaR[ic] = 1.8 * gRandom->Rndm() - 0.9;
 265           pt  [ic] = track->Pt();
 266           if (pt[ic] > ptmax) {
 267               ptmax = pt[ic];
 268               icmax = ic;
 269           }
 270           ic++;
 271       }
 272   } //track loop
 273
 274   //
 275   // Cluster
 276   if (ic < fNMin) {
 277       PostData(1, fHists);
 278       return;
 279   }
 280   //
 281   Double_t* rk0 = new Double_t[10];
 282
 283   AliKMeansResult* res = 0;
 284
 285   for (Int_t i = 0; i < fK; i++) {
 286     res = fA[i];
 287     AliKMeansClustering::SoftKMeans2(i+1, ic, phi, eta, res->GetMx(),res->GetMy(), res->GetSigma2(), res->GetRk());
 288     res->Sort(ic, phi, eta);
 289     Int_t j = (res->GetInd())[0];
 290     rk0[i]  = (res->GetTarget())[j];
 291   }
 292
 293   Double_t* mPhi   = 0;
 294   Double_t* mEta   = 0;
 295   Double_t* sigma2 = 0;
 296   Float_t   rmax   = -1.;
 297   Int_t     imax   = 0;
 298
 299   for (Int_t i = 0; i < fK; i++) {
 300     if (rk0[i] > rmax)  {
 301         rmax = rk0[i];
 302         imax = i;
 303       }
 304   }
 305
 306
 307   res = fA[imax];
 308   mPhi    = res->GetMx();
 309   mEta    = res->GetMy();
 310   sigma2  = res->GetSigma2();
 311   Int_t nk = imax + 1;
 312   Int_t im = (res->GetInd())[0];
 313   fHDensity->Fill(rmax / TMath::Pi());
 314   fHCSize->Fill(rmax * sigma2[im]);
 315   fHNCluster->Fill(Float_t(nk));
 316
 317   Double_t dphic, detac;
 318
 319   if (rmax > 0.) {
 320     // Analysis
 321     //
 322     // Cluster Multiplicity
 323     Int_t mult[2] = {0, 0};
 324     //
 325     if (nk > 1) {
 326       dphic = DeltaPhi(mPhi[0], mPhi[1]);
 327       detac = TMath::Abs(mEta[0] - mEta[1]);
 328       fH2DPhiEtaC->Fill(dphic, detac);
 329     }
 330     //
 331     // Random cluster position
 332     Int_t ir = Int_t(Float_t(ic) * gRandom->Rndm());
 333
 334     Double_t crPhi = phi[ir];
 335     Double_t crEta = eta[ir];
 336     //
 337     Double_t sumPt  = 0;
 338     Double_t sumPtC = 0;
 339
 340     for (Int_t i = 0; i < 1; i++) {
 341       fH1CEta->Fill(mEta[im]);
 342       fH1CPhi->Fill(mPhi[im]);
 343       for (Int_t j = 0; j < ic; j++) {
 344         Double_t r    = DeltaR(mPhi[im], mEta[im], phi[j], eta[j]);
 345         Double_t dphi = DeltaPhi(mPhi[im], phi[j]);
 346         Double_t deta = mEta[im] - eta[j];
 347         Double_t rr   = DeltaR(crPhi, crEta, phi[j], eta[j]);
 348
 349         fH1DR->Fill(r);
 350         fH1DPhi->Fill(dphi);
 351         fH2DPhiEta->Fill(dphi, TMath::Abs(deta));
 352         if (j == icmax) fH2DPhiEtaL->Fill(dphi, TMath::Abs(deta));
 353
 354         if (r < 0.2) {
 355           fH1PtC2->Fill(pt[j]);
 356         }
 357         if (r < 0.3) {
 358           fH1PtC1->Fill(pt[j]);
 359           fHPtDensity->Fill(rmax/TMath::Pi(), pt[j]);
 360         }
 361         if (rr < 0.3) {
 362             fH1PtR->Fill(pt[j]);
 363         }
 364
 365         if (r < 0.4)  {
 366           sumPtC += pt[j];
 367           mult[i]++;
 368           fH1PtC->Fill(pt[j]);
 369         }
 370         if (r > 0.7 && dphi < (TMath::Pi() - 0.3)) {
 371           fH1Pt->Fill(pt[j]);
 372         }
 373
 374         if (r > 0.7 && r < 1.) {
 375           sumPt += pt[j];
 376         }
 377
 378         if (dphi > (TMath::Pi() - 0.3)) {
 379           fH1PtAS->Fill(pt[j]);
 380         }
 381       }
 382     }
 383
 384     fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
 385     fH1SPt ->Fill(sumPt);
 386     fH1SPtC->Fill(sumPtC);
 387   }
 388     //
 389     // Randomized phi
 390     //
 391     for (Int_t i = 0; i < fK; i++) {
 392       res = fB[i];
 393       AliKMeansClustering::SoftKMeans2(i+1, ic, phiR, etaR, res->GetMx(),res->GetMy(), res->GetSigma2(), res->GetRk());
 394       res->Sort(ic, phiR, etaR);
 395       //res->Sort();
 396       Int_t j = (res->GetInd())[0];
 397       rk0[i]  = (res->GetTarget())[j];
 398     }
 399
 400     rmax   = -1.;
 401     imax   = 0;
 402     for (Int_t i = 0; i < fK; i++) {
 403       if (rk0[i] > rmax) {
 404           rmax = rk0[i];
 405           imax = i;
 406         }
 407     }
 408     res = fB[imax];
 409     mPhi    = res->GetMx();
 410     mEta    = res->GetMy();
 411     sigma2  = res->GetSigma2();
 412     nk = imax + 1;
 413
 414     //  fH1RespR->Fill(rk[ind[0]]/(rk[0]+rk[1]));
 415     //
 416     // Cluster Multiplicity
 417     if (rmax > 0.) {
 418       for (Int_t i = 0; i < 1; i++) {
 419         Int_t im = (res->GetInd())[i];
 420         fH1CEtaR->Fill(mEta[im]);
 421       }
 422
 423       for (Int_t i = 0; i < 1; i++) {
 424         Int_t im = (res->GetInd())[i];
 425         for (Int_t j = 0; j < ic; j++) {
 426           Double_t dphi = DeltaPhi(mPhi[im], phiR[j]);
 427           Double_t deta = mEta[im] - etaR[j];
 428           Double_t r = DeltaR(mPhi[im], mEta[im], phiR[j], etaR[j]);
 429           fH1DRR->Fill(r);
 430           fH2DPhiEtaR->Fill(dphi, TMath::Abs(deta));
 431           if (j == icmax) fH2DPhiEtaLR->Fill(dphi, TMath::Abs(deta));
 432         }
 433       }
 434       if (nk > 1) {
 435         dphic = DeltaPhi(mPhi[0], mPhi[1]);
 436         detac = TMath::Abs(mEta[0] - mEta[1]);
 437         fH2DPhiEtaCR->Fill(dphic, detac);
 438       }
 439     }
 440     PostData(1, fHists);
 441 }
 442
 443 Double_t AliAnalysisTaskKMeans::DeltaPhi(Double_t phi1, Double_t phi2)
 444 {
 445     Double_t dphi = TMath::Abs(phi1 - phi2);
 446     if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
 447     return dphi;
 448 }
 449
 450 Double_t AliAnalysisTaskKMeans::DeltaR(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
 451 {
 452     Double_t dphi = DeltaPhi(phi1, phi2);
 453     Double_t deta = eta1 - eta2;
 454     return (TMath::Sqrt(dphi * dphi + deta * deta));
 455
 456 }
 457
 458 //________________________________________________________________________
 459 void AliAnalysisTaskKMeans::Terminate(Option_t *)
 460 {
 461 }
 462