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