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