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