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