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