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