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