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