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
  41 #include "AliAnalysisTask.h"
  42 #include "AliAnalysisManager.h"
  43
  44 #include "AliESDEvent.h"
  45 #include "AliStack.h"
  46 #include "AliESDVertex.h"
  47 #include "AliESDInputHandler.h"
  48 #include "AliESDtrackCuts.h"
  49 #include "AliMultiplicity.h"
  50
  51 #include "AliMCParticle.h"
  52 #include "AliMCEvent.h"
  53 #include "AliAnalysisTaskKMeans.h"
  54 #include "AliTrackReference.h"
  55 #include "AliStack.h"
  56 #include "AliHeader.h"
  57 #include "AliKMeansClustering.h"
  58
  59
  60
  61 ClassImp(AliAnalysisTaskKMeans)
  62
  63 AliAnalysisTaskKMeans::AliAnalysisTaskKMeans()
  64     : AliAnalysisTaskSE()
  65     ,fHists(0)
  66     ,fH1CEta(0)
  67     ,fH1CPhi(0)
  68     ,fH1CEtaR(0)
  69     ,fH2N1N2(0)
  70     ,fH1Pt(0)
  71     ,fH1PtC(0)
  72     ,fH1PtC1(0)
  73     ,fH1PtC2(0)
  74     ,fH1SPt(0)
  75     ,fH1SPtC(0)
  76     ,fH1DPhi(0)
  77     ,fH1DR(0)
  78     ,fH1DRR(0)
  79     ,fH2DPhiEta(0)
  80     ,fH2DPhiEtaR(0)
  81     ,fH2DPhiEtaL(0)
  82     ,fH2DPhiEtaLR(0)
  83     ,fH2DPhiEtaC(0)
  84     ,fH2DPhiEtaCR(0)
  85     ,fCuts(0)
  86 {
  87   //
  88   // Constructor
  89   //
  90 }
  91
  92 //________________________________________________________________________
  93 AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name)
  94     : AliAnalysisTaskSE(name)
  95       ,fHists(0)
  96       ,fH1CEta(0)
  97       ,fH1CPhi(0)
  98       ,fH1CEtaR(0)
  99       ,fH2N1N2(0)
 100       ,fH1Pt(0)
 101       ,fH1PtC(0)
 102       ,fH1PtC1(0)
 103       ,fH1PtC2(0)
 104       ,fH1SPt(0)
 105       ,fH1SPtC(0)
 106       ,fH1DPhi(0)
 107       ,fH1DR(0)
 108       ,fH1DRR(0)
 109       ,fH2DPhiEta(0)
 110       ,fH2DPhiEtaR(0)
 111       ,fH2DPhiEtaL(0)
 112       ,fH2DPhiEtaLR(0)
 113       ,fH2DPhiEtaC(0)
 114       ,fH2DPhiEtaCR(0)
 115       ,fCuts(0)
 116 {
 117   //
 118   // Constructor
 119   //
 120   DefineOutput(1,  TList::Class());
 121 }
 122
 123
 124 //________________________________________________________________________
 125 void AliAnalysisTaskKMeans::UserCreateOutputObjects()
 126 {
 127   // Create histograms
 128   // Called once
 129     fHists   = new TList();
 130     fH1CEta  = new TH1F("fH1CEta",   "eta distribution of clusters",        90, -1.5, 1.5);
 131     fH1CEtaR = new TH1F("fH1CEtaR",  "eta distribution of clusters",        90, -1.5, 1.5);
 132     fH1CPhi  = new TH1F("fH1CPhi",   "phi distribution of clusters",       157,  0.0, 2. * TMath::Pi());
 133     fH2N1N2  = new TH2F("fH2N1N2",   "multiplicity distribution",          50, 0., 50., 50, 0., 50.);
 134
 135     fH1Pt    = new TH1F("fH1Pt",     "pt distribution",50, 0., 10.);
 136     fH1PtC   = new TH1F("fH1PtC",    "pt distribution",50, 0., 10.);
 137     fH1PtC1  = new TH1F("fH1PtC1",   "pt distribution",50, 0., 10.);
 138     fH1PtC2  = new TH1F("fH1PtC2",   "pt distribution",50, 0., 10.);
 139
 140     fH1SPt    = new TH1F("fH1SPt",     "sum pt distribution",50, 0., 10.);
 141     fH1SPtC   = new TH1F("fH1SPtC",    "sum pt distribution",50, 0., 10.);
 142
 143     fH1DR         = new TH1F("fH1DR",     "dR distribution", 50, 0., 5.);
 144     fH1DPhi       = new TH1F("fH1DPhi",   "dPhi distribution", 31, 0., TMath::Pi());
 145     fH2DPhiEta    = new TH2F("fH2DPhiEta","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 146     fH2DPhiEtaR   = new TH2F("fH2DPhiEtaR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 147     fH2DPhiEtaL   = new TH2F("fH2DPhiEtaL","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 148     fH2DPhiEtaLR  = new TH2F("fH2DPhiEtaLR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 149     fH2DPhiEtaC   = new TH2F("fH2DPhiEtaC","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 150     fH2DPhiEtaCR  = new TH2F("fH2DPhiEtaCR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
 151     fH1DRR        = new TH1F("fH1DRR",    "dR distribution", 50, 0., 5.);
 152
 153     fHists->SetOwner();
 154
 155     fHists->Add(fH1CEta);
 156     fHists->Add(fH1CEtaR);
 157     fHists->Add(fH1CPhi);
 158     fHists->Add(fH2N1N2);
 159     fHists->Add(fH1Pt);
 160     fHists->Add(fH1PtC);
 161     fHists->Add(fH1PtC1);
 162     fHists->Add(fH1PtC2);
 163     fHists->Add(fH1DR);
 164     fHists->Add(fH1SPtC);
 165     fHists->Add(fH1SPt);
 166     fHists->Add(fH1DPhi);
 167     fHists->Add(fH2DPhiEta);
 168     fHists->Add(fH2DPhiEtaR);
 169     fHists->Add(fH2DPhiEtaL);
 170     fHists->Add(fH2DPhiEtaLR);
 171     fHists->Add(fH2DPhiEtaC);
 172     fHists->Add(fH2DPhiEtaCR);
 173     fHists->Add(fH1DRR);
 174
 175     //
 176     AliKMeansClustering::SetBeta(20.);
 177
 178 }
 179
 180 //________________________________________________________________________
 181 void AliAnalysisTaskKMeans::UserExec(Option_t *)
 182 {
 183   // Main loop
 184   // Called for each event
 185     Double_t   phi [500];
 186     Double_t   phiR[500];
 187     Double_t   eta[500];
 188     Double_t   etaR[500];
 189     Double_t   pt [500];
 190     Double_t   mPhi[20];
 191     Double_t   mEta[20];
 192     Double_t     rk[20];
 193     Int_t       ind[20];
 194
 195   if (!fInputEvent) {
 196     Printf("ERROR: fESD not available");
 197     return;
 198   }
 199
 200   AliESDEvent* esdE = (AliESDEvent*) fInputEvent;
 201
 202   Int_t ic = 0;
 203   //
 204   // Fill eta-phi positions
 205   //
 206   /*
 207   const AliMultiplicity *spdMult = esdE->GetMultiplicity();
 208   for (Int_t i = 0; i < spdMult->GetNumberOfTracklets(); ++i) {
 209           phi[ic] = spdMult->GetPhi(i);
 210           eta[ic] = spdMult->GetEta(i);
 211           ic++;
 212   }
 213   */
 214   Double_t ptmax = 0.;
 215   Int_t    icmax = -1;
 216
 217   for (Int_t iTracks = 0; iTracks < esdE->GetNumberOfTracks(); iTracks++) {
 218       AliESDtrack* track = esdE->GetTrack(iTracks);
 219       if ((fCuts->AcceptTrack(track)))
 220       {
 221           phi[ic]  = track->Phi();
 222           phiR[ic] = 2. * TMath::Pi() * gRandom->Rndm();
 223           eta[ic]  = track->Eta();
 224           etaR[ic] = 2. * gRandom->Rndm() - 1.;
 225           pt[ic]   = track->Pt();
 226           if (pt[ic] > ptmax) {
 227               ptmax = pt[ic];
 228               icmax = ic;
 229           }
 230           ic++;
 231       }
 232   } //track loop
 233
 234   //
 235   // Cluster
 236   if (ic < 10) {
 237       PostData(1, fHists);
 238       return;
 239   }
 240   //
 241   AliKMeansClustering::SoftKMeans(2, ic, phi, eta, mPhi, mEta, rk);
 242   //
 243   // Sort
 244   TMath::Sort(2, rk, ind);
 245   //
 246   // Analyse
 247   //
 248   // Cluster Multiplicity
 249   Int_t mult[2] = {0, 0};
 250   //
 251   Double_t dphic = DeltaPhi(mPhi[0], mPhi[1]);
 252   Double_t detac = TMath::Abs(mEta[0] - mEta[1]);
 253   fH2DPhiEtaC->Fill(dphic, detac);
 254
 255   //
 256   Double_t sumPt  = 0;
 257   Double_t sumPtC = 0;
 258   for (Int_t i = 0; i < 1; i++) {
 259       fH1CEta->Fill(mEta[ind[i]]);
 260       fH1CPhi->Fill(mPhi[ind[i]]);
 261       for (Int_t j = 0; j < ic; j++) {
 262           Double_t r    = DeltaR(mPhi[ind[i]], mEta[ind[i]], phi[j], eta[j]);
 263           Double_t dphi = DeltaPhi(mPhi[ind[i]], phi[j]);
 264           Double_t deta = mEta[ind[i]] - eta[j];
 265
 266           fH1DR->Fill(r);
 267           fH1DPhi->Fill(dphi);
 268           fH2DPhiEta->Fill(dphi, TMath::Abs(deta));
 269           if (j == icmax) fH2DPhiEtaL->Fill(dphi, TMath::Abs(deta));
 270
 271           if (r < 0.2) {
 272               fH1PtC2->Fill(pt[j]);
 273           }
 274           if (r < 0.3)
 275           {
 276               fH1PtC1->Fill(pt[j]);
 277           }
 278           if (r < 0.4)
 279           {
 280             sumPtC += pt[j];
 281             mult[i]++;
 282             fH1PtC->Fill(pt[j]);
 283           }
 284           if (r > 0.7) {
 285               fH1Pt->Fill(pt[j]);
 286           }
 287
 288           if (r > 0.7 && r < 1.) {
 289             sumPt += pt[j];
 290           }
 291       }
 292   }
 293
 294   fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
 295   fH1SPt ->Fill(sumPt);
 296   fH1SPtC->Fill(sumPtC);
 297
 298
 299   // Randomized phi
 300   //
 301   AliKMeansClustering::SoftKMeans(2, ic, phiR, etaR, mPhi, mEta, rk);
 302   //
 303   // Sort
 304   TMath::Sort(2, rk, ind);
 305   //
 306   // Analyse
 307   //
 308   // Cluster Multiplicity
 309   for (Int_t i = 0; i < 1; i++) {
 310       fH1CEtaR->Fill(mEta[ind[i]]);
 311   }
 312
 313   for (Int_t i = 0; i < 1; i++) {
 314       for (Int_t j = 0; j < ic; j++) {
 315           Double_t dphi = DeltaPhi(mPhi[ind[i]], phiR[j]);
 316           Double_t deta = mEta[ind[i]] - etaR[j];
 317           Double_t r = DeltaR(mPhi[ind[i]], mEta[ind[i]], phiR[j], etaR[j]);
 318           fH1DRR->Fill(r);
 319           fH2DPhiEtaR->Fill(dphi, TMath::Abs(deta));
 320           if (j == icmax) fH2DPhiEtaLR->Fill(dphi, TMath::Abs(deta));
 321       }
 322   }
 323   dphic = DeltaPhi(mPhi[0], mPhi[1]);
 324   detac = TMath::Abs(mEta[0] - mEta[1]);
 325   fH2DPhiEtaCR->Fill(dphic, detac);
 326
 327
 328   //
 329   // Post output data.
 330   PostData(1, fHists);
 331 }
 332
 333 Double_t AliAnalysisTaskKMeans::DeltaPhi(Double_t phi1, Double_t phi2)
 334 {
 335     Double_t dphi = TMath::Abs(phi1 - phi2);
 336     if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
 337     return dphi;
 338 }
 339
 340 Double_t AliAnalysisTaskKMeans::DeltaR(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
 341 {
 342     Double_t dphi = DeltaPhi(phi1, phi2);
 343     Double_t deta = eta1 - eta2;
 344     return (TMath::Sqrt(dphi * dphi + deta * deta));
 345
 346 }
 347
 348 //________________________________________________________________________
 349 void AliAnalysisTaskKMeans::Terminate(Option_t *)
 350 {
 351 }
 352