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 AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const AliAnalysisTaskKMeans &res)
 150 : AliAnalysisTaskSE(res)
 151       ,fK(0)
 152       ,fNMin(0)
 153       ,fHists(0)
 154       ,fH1CEta(0)
 155       ,fH1CPhi(0)
 156       ,fH1CEtaR(0)
 157       ,fH2N1N2(0)
 158       ,fH1Pt(0)
 159       ,fH1PtC(0)
 160       ,fH1PtC1(0)
 161       ,fH1PtC2(0)
 162       ,fH1PtAS(0)
 163       ,fH1PtR(0)
 164       ,fH1SPt(0)
 165       ,fH1SPtC(0)
 166       ,fH1DPhi(0)
 167       ,fH1DR(0)
 168       ,fH1DRR(0)
 169       ,fH2DPhiEta(0)
 170       ,fH2DPhiEtaR(0)
 171       ,fH2DPhiEtaL(0)
 172       ,fH2DPhiEtaLR(0)
 173       ,fH2DPhiEtaC(0)
 174       ,fH2DPhiEtaCR(0)
 175       ,fH1Resp(0)
 176       ,fH1RespR(0)
 177       ,fH2Sigma(0)
 178       ,fH2SigmaR(0)
 179       ,fCuts(0)
 180 {
 181     // Dummy copy constructor
 182 }
 183
 184 AliAnalysisTaskKMeans& AliAnalysisTaskKMeans::operator=(const AliAnalysisTaskKMeans& /*trclass*/)
 185 {
 186     // Dummy assignment operator
 187     return *this;
 188 }
 189
 190
 191
 192 //________________________________________________________________________
 193 void AliAnalysisTaskKMeans::UserCreateOutputObjects()
 194 {
 195   // Create histograms
 196   // Called once
 197     fHists   = new TList();
 198     fH1CEta  = new TH1F("fH1CEta",   "eta distribution of clusters",        90, -1.5, 1.5);
 199     fH1CEtaR = new TH1F("fH1CEtaR",  "eta distribution of clusters",        90, -1.5, 1.5);
 200     fH1CPhi  = new TH1F("fH1CPhi",   "phi distribution of clusters",       157,  0.0, 2. * TMath::Pi());
 201     fH2N1N2  = new TH2F("fH2N1N2",   "multiplicity distribution",          50, 0., 50., 50, 0., 50.);
 202
 203     fH1Pt    = new TH1F("fH1Pt",     "pt distribution",50, 0., 10.);
 204     fH1PtC   = new TH1F("fH1PtC",    "pt distribution",50, 0., 10.);
 205     fH1PtC1  = new TH1F("fH1PtC1",   "pt distribution",50, 0., 10.);
 206     fH1PtC2  = new TH1F("fH1PtC2",   "pt distribution",50, 0., 10.);
 207     fH1PtAS  = new TH1F("fH1PtAS",   "pt distribution",50, 0., 10.);
 208     fH1PtR   = new TH1F("fH1PtR",    "pt distribution",50, 0., 10.);
 209
 210     fH1SPt    = new TH1F("fH1SPt",     "sum pt distribution",50, 0., 10.);
 211     fH1SPtC   = new TH1F("fH1SPtC",    "sum pt distribution",50, 0., 10.);
 212
 213     fH1DR         = new TH1F("fH1DR",     "dR distribution", 50, 0., 5.);
 214     fH1DPhi       = new TH1F("fH1DPhi",   "dPhi distribution", 31, 0., TMath::Pi());
 215     fH2DPhiEta    = new TH2F("fH2DPhiEta","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 216     fH2DPhiEtaR   = new TH2F("fH2DPhiEtaR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 217     fH2DPhiEtaL   = new TH2F("fH2DPhiEtaL","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 218     fH2DPhiEtaLR  = new TH2F("fH2DPhiEtaLR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 219     fH2DPhiEtaC   = new TH2F("fH2DPhiEtaC","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 220     fH2DPhiEtaCR  = new TH2F("fH2DPhiEtaCR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 221     fH1DRR        = new TH1F("fH1DRR",    "dR distribution", 50, 0., 5.);
 222     fH1Resp       = new TH1F("fH1Resp",   "Responsibility", 50, 0., 1.);
 223     fH1RespR      = new TH1F("fH1RespR",  "Responsibility", 50, 0., 1.);
 224     fH2Sigma      = new TH2F("fH2Sigma",  "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.);
 225     fH2SigmaR     = new TH2F("fH2SigmaR", "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.);
 226
 227
 228     fHDensity    = new TH1F("fHDensity",    "density distribution", 100, 0., 20.);
 229     fHCSize      = new TH1F("fHCSize",      "cluster size", 20, -0.5, 19.5);
 230
 231     fHNCluster   = new TH1F("fHNCluster",   "Number of clusters", 11, -0.5, 10.5);
 232     fHPtDensity  = new TH2F("fHPtDensity", "Pt vs density", 100, 0., 20., 50, 0., 10.);
 233
 234
 235     fHists->SetOwner();
 236
 237     fHists->Add(fH1CEta);
 238     fHists->Add(fH1CEtaR);
 239     fHists->Add(fH1CPhi);
 240     fHists->Add(fH2N1N2);
 241     fHists->Add(fH1Pt);
 242     fHists->Add(fH1PtR);
 243     fHists->Add(fH1PtC);
 244     fHists->Add(fH1PtC1);
 245     fHists->Add(fH1PtC2);
 246     fHists->Add(fH1PtAS);
 247     fHists->Add(fH1DR);
 248     fHists->Add(fH1SPtC);
 249     fHists->Add(fH1SPt);
 250     fHists->Add(fH1DPhi);
 251     fHists->Add(fH2DPhiEta);
 252     fHists->Add(fH2DPhiEtaR);
 253     fHists->Add(fH2DPhiEtaL);
 254     fHists->Add(fH2DPhiEtaLR);
 255     fHists->Add(fH2DPhiEtaC);
 256     fHists->Add(fH2DPhiEtaCR);
 257     fHists->Add(fH1DRR);
 258     fHists->Add(fH1RespR);
 259     fHists->Add(fH1Resp);
 260     fHists->Add(fH2Sigma);
 261     fHists->Add(fH2SigmaR);
 262     fHists->Add(fHCSize);
 263     fHists->Add(fHDensity);
 264     fHists->Add(fHNCluster);
 265     fHists->Add(fHPtDensity);
 266     //
 267     for (Int_t i = 0; i < 10; i++) {
 268       fA[i] = new AliKMeansResult(i+1);
 269       fB[i] = new AliKMeansResult(i+1);
 270     }
 271 }
 272
 273 //________________________________________________________________________
 274 void AliAnalysisTaskKMeans::UserExec(Option_t *)
 275 {
 276   // Main loop
 277   // Called for each event
 278     Double_t   phi [500];
 279     Double_t   phiR[500];
 280     Double_t    eta[500];
 281     Double_t   etaR[500];
 282     Double_t    pt [500];
 283
 284   if (!fInputEvent) {
 285     Printf("ERROR: fESD not available");
 286     return;
 287   }
 288
 289
 290   Int_t ic = 0;
 291
 292   //
 293   // Fill eta-phi positions
 294   //
 295
 296   Double_t ptmax = 0.;
 297   Int_t    icmax = -1;
 298
 299   for (Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++) {
 300       AliVParticle* track = fInputEvent->GetTrack(iTracks);
 301       if ((fCuts->AcceptTrack((AliESDtrack*)track)))
 302       {
 303           phi [ic] = track->Phi();
 304           phiR[ic] = 2. * TMath::Pi() * gRandom->Rndm();
 305           eta [ic] = track->Eta();
 306           etaR[ic] = 1.8 * gRandom->Rndm() - 0.9;
 307           pt  [ic] = track->Pt();
 308           if (pt[ic] > ptmax) {
 309               ptmax = pt[ic];
 310               icmax = ic;
 311           }
 312           ic++;
 313       }
 314   } //track loop
 315
 316   //
 317   // Cluster
 318   if (ic < fNMin) {
 319       PostData(1, fHists);
 320       return;
 321   }
 322   //
 323   Double_t rk0[10];
 324
 325   AliKMeansResult* res = 0;
 326
 327   for (Int_t i = 0; i < fK; i++) {
 328     res = fA[i];
 329     AliKMeansClustering::SoftKMeans2(i+1, ic, phi, eta, res->GetMx(),res->GetMy(), res->GetSigma2(), res->GetRk());
 330     res->Sort(ic, phi, eta);
 331     Int_t j = (res->GetInd())[0];
 332     rk0[i]  = (res->GetTarget())[j];
 333   }
 334
 335   Double_t* mPhi   = 0;
 336   Double_t* mEta   = 0;
 337   Double_t* sigma2 = 0;
 338   Float_t   rmax   = -1.;
 339   Int_t     imax   = 0;
 340
 341   for (Int_t i = 0; i < fK; i++) {
 342     if (rk0[i] > rmax)  {
 343         rmax = rk0[i];
 344         imax = i;
 345       }
 346   }
 347
 348
 349   res = fA[imax];
 350   mPhi    = res->GetMx();
 351   mEta    = res->GetMy();
 352   sigma2  = res->GetSigma2();
 353   Int_t nk = imax + 1;
 354   Int_t im = (res->GetInd())[0];
 355   fHDensity->Fill(rmax / TMath::Pi());
 356   fHCSize->Fill(rmax * sigma2[im]);
 357   fHNCluster->Fill(Float_t(nk));
 358
 359   Double_t dphic, detac;
 360
 361   if (rmax > 0.) {
 362     // Analysis
 363     //
 364     // Cluster Multiplicity
 365     Int_t mult[2] = {0, 0};
 366     //
 367     if (nk > 1) {
 368       dphic = DeltaPhi(mPhi[0], mPhi[1]);
 369       detac = TMath::Abs(mEta[0] - mEta[1]);
 370       fH2DPhiEtaC->Fill(dphic, detac);
 371     }
 372     //
 373     // Random cluster position
 374     Int_t ir = Int_t(Float_t(ic) * gRandom->Rndm());
 375
 376     Double_t crPhi = phi[ir];
 377     Double_t crEta = eta[ir];
 378     //
 379     Double_t sumPt  = 0;
 380     Double_t sumPtC = 0;
 381
 382     for (Int_t i = 0; i < 1; i++) {
 383       fH1CEta->Fill(mEta[im]);
 384       fH1CPhi->Fill(mPhi[im]);
 385       for (Int_t j = 0; j < ic; j++) {
 386         Double_t r    = DeltaR(mPhi[im], mEta[im], phi[j], eta[j]);
 387         Double_t dphi = DeltaPhi(mPhi[im], phi[j]);
 388         Double_t deta = mEta[im] - eta[j];
 389         Double_t rr   = DeltaR(crPhi, crEta, phi[j], eta[j]);
 390
 391         fH1DR->Fill(r);
 392         fH1DPhi->Fill(dphi);
 393         fH2DPhiEta->Fill(dphi, TMath::Abs(deta));
 394         if (j == icmax) fH2DPhiEtaL->Fill(dphi, TMath::Abs(deta));
 395
 396         if (r < 0.2) {
 397           fH1PtC2->Fill(pt[j]);
 398         }
 399         if (r < 0.3) {
 400           fH1PtC1->Fill(pt[j]);
 401           fHPtDensity->Fill(rmax/TMath::Pi(), pt[j]);
 402         }
 403         if (rr < 0.3) {
 404             fH1PtR->Fill(pt[j]);
 405         }
 406
 407         if (r < 0.4)  {
 408           sumPtC += pt[j];
 409           mult[i]++;
 410           fH1PtC->Fill(pt[j]);
 411         }
 412         if (r > 0.7 && dphi < (TMath::Pi() - 0.3)) {
 413           fH1Pt->Fill(pt[j]);
 414         }
 415
 416         if (r > 0.7 && r < 1.) {
 417           sumPt += pt[j];
 418         }
 419
 420         if (dphi > (TMath::Pi() - 0.3)) {
 421           fH1PtAS->Fill(pt[j]);
 422         }
 423       }
 424     }
 425
 426     fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
 427     fH1SPt ->Fill(sumPt);
 428     fH1SPtC->Fill(sumPtC);
 429   }
 430     //
 431     // Randomized phi
 432     //
 433     for (Int_t i = 0; i < fK; i++) {
 434       res = fB[i];
 435       AliKMeansClustering::SoftKMeans2(i+1, ic, phiR, etaR, res->GetMx(),res->GetMy(), res->GetSigma2(), res->GetRk());
 436       res->Sort(ic, phiR, etaR);
 437       //res->Sort();
 438       Int_t j = (res->GetInd())[0];
 439       rk0[i]  = (res->GetTarget())[j];
 440     }
 441
 442     rmax   = -1.;
 443     imax   = 0;
 444     for (Int_t i = 0; i < fK; i++) {
 445       if (rk0[i] > rmax) {
 446           rmax = rk0[i];
 447           imax = i;
 448         }
 449     }
 450     res = fB[imax];
 451     mPhi    = res->GetMx();
 452     mEta    = res->GetMy();
 453     sigma2  = res->GetSigma2();
 454     nk = imax + 1;
 455
 456     //  fH1RespR->Fill(rk[ind[0]]/(rk[0]+rk[1]));
 457     //
 458     // Cluster Multiplicity
 459     if (rmax > 0.) {
 460       for (Int_t i = 0; i < 1; i++) {
 461         Int_t im = (res->GetInd())[i];
 462         fH1CEtaR->Fill(mEta[im]);
 463       }
 464
 465       for (Int_t i = 0; i < 1; i++) {
 466         Int_t im = (res->GetInd())[i];
 467         for (Int_t j = 0; j < ic; j++) {
 468           Double_t dphi = DeltaPhi(mPhi[im], phiR[j]);
 469           Double_t deta = mEta[im] - etaR[j];
 470           Double_t r = DeltaR(mPhi[im], mEta[im], phiR[j], etaR[j]);
 471           fH1DRR->Fill(r);
 472           fH2DPhiEtaR->Fill(dphi, TMath::Abs(deta));
 473           if (j == icmax) fH2DPhiEtaLR->Fill(dphi, TMath::Abs(deta));
 474         }
 475       }
 476       if (nk > 1) {
 477         dphic = DeltaPhi(mPhi[0], mPhi[1]);
 478         detac = TMath::Abs(mEta[0] - mEta[1]);
 479         fH2DPhiEtaCR->Fill(dphic, detac);
 480       }
 481     }
 482     PostData(1, fHists);
 483 }
 484
 485 Double_t AliAnalysisTaskKMeans::DeltaPhi(Double_t phi1, Double_t phi2)
 486 {
 487     Double_t dphi = TMath::Abs(phi1 - phi2);
 488     if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
 489     return dphi;
 490 }
 491
 492 Double_t AliAnalysisTaskKMeans::DeltaR(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
 493 {
 494     Double_t dphi = DeltaPhi(phi1, phi2);
 495     Double_t deta = eta1 - eta2;
 496     return (TMath::Sqrt(dphi * dphi + deta * deta));
 497
 498 }
 499
 500 //________________________________________________________________________
 501 void AliAnalysisTaskKMeans::Terminate(Option_t *)
 502 {
 503 }
 504