c47df6f0f5cdd85940b8519c1bab542b29a1a81c
[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 "TFile.h"
30 #include "TTree.h"
31 #include "TH1F.h"
32 #include "TH2F.h"
33 #include "TProfile.h"
34
35 #include "TList.h"
36 #include "TParticle.h"
37 #include "TParticlePDG.h"
38 #include "TProfile.h"
39 #include "TMath.h"
40 #include "TRandom.h"
41 #include "TObjArray.h"
42
43 #include "AliAnalysisTask.h"
44 #include "AliAnalysisManager.h"
45
46 #include "AliVEvent.h"
47 #include "AliESDEvent.h"
48 #include "AliExternalTrackParam.h"
49 #include "AliStack.h"
50 #include "AliESDVertex.h"
51 #include "AliESDInputHandler.h"
52 #include "AliESDtrackCuts.h"
53 #include "AliMultiplicity.h"
54
55 #include "AliMCParticle.h"
56 #include "AliMCEvent.h"
57 #include "AliAnalysisTaskKMeans.h"
58 #include "AliTrackReference.h"
59 #include "AliStack.h"
60 #include "AliHeader.h"
61 #include "AliKMeansClustering.h"
62
63
64
65 ClassImp(AliAnalysisTaskKMeans)
66
67 AliAnalysisTaskKMeans::AliAnalysisTaskKMeans() 
68     : AliAnalysisTaskSE()
69     ,fK(0)
70     ,fNMin(0)
71     ,fHists(0)
72     ,fH1CEta(0)
73     ,fH1CPhi(0)
74     ,fH1CEtaR(0)
75     ,fH2N1N2(0)
76     ,fH1Pt(0)
77     ,fH1PtC(0)
78     ,fH1PtC1(0)
79     ,fH1PtC2(0)
80     ,fH1PtAS(0)
81     ,fH1PtR(0)
82     ,fH1SPt(0)
83     ,fH1SPtC(0)
84     ,fH1DPhi(0)
85     ,fH1DR(0)
86     ,fH1DRR(0)
87     ,fH2DPhiEta(0)
88     ,fH2DPhiEtaR(0)
89     ,fH2DPhiEtaL(0)
90     ,fH2DPhiEtaLR(0)
91     ,fH2DPhiEtaC(0)
92     ,fH2DPhiEtaCR(0)
93     ,fH1Resp(0)
94     ,fH1RespR(0)
95     ,fH2Sigma(0)
96     ,fH2SigmaR(0)
97     ,fHDensity(0)
98     ,fHCSize(0)
99     ,fHNCluster(0)
100     ,fHPtDensity(0)
101     ,fHDPhi(0)
102     ,fH2EtaPhi(0)
103     ,fH2REtaPhi(0)
104     ,fCuts(0)
105
106 {
107   //
108
109   // Constructor
110   //
111     for (Int_t i=0; i< 10; i++) {
112         fA[i] = 0;
113         fB[i] = 0;
114     }
115 }
116
117 //________________________________________________________________________
118 AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const char *name) 
119     : AliAnalysisTaskSE(name) 
120       ,fK(0)
121       ,fNMin(0)
122       ,fHists(0)
123       ,fH1CEta(0)
124       ,fH1CPhi(0)
125       ,fH1CEtaR(0)
126       ,fH2N1N2(0)
127       ,fH1Pt(0)
128       ,fH1PtC(0)
129       ,fH1PtC1(0)
130       ,fH1PtC2(0)
131       ,fH1PtAS(0)
132       ,fH1PtR(0)
133       ,fH1SPt(0)
134       ,fH1SPtC(0)
135       ,fH1DPhi(0)
136       ,fH1DR(0)
137       ,fH1DRR(0)
138       ,fH2DPhiEta(0)
139       ,fH2DPhiEtaR(0)
140       ,fH2DPhiEtaL(0)
141       ,fH2DPhiEtaLR(0)
142       ,fH2DPhiEtaC(0)
143       ,fH2DPhiEtaCR(0)
144       ,fH1Resp(0)
145       ,fH1RespR(0)
146       ,fH2Sigma(0)
147       ,fH2SigmaR(0)
148       ,fHDensity(0)
149       ,fHCSize(0)
150       ,fHNCluster(0)
151       ,fHPtDensity(0)
152       ,fHDPhi(0)
153       ,fH2EtaPhi(0)
154       ,fH2REtaPhi(0)
155       ,fCuts(0)
156
157 {
158   //
159   // Constructor
160   //
161     for (Int_t i=0; i< 10; i++) {
162         fA[i] = 0;
163         fB[i] = 0;
164     }
165   DefineOutput(1,  TList::Class());
166 }
167
168 AliAnalysisTaskKMeans::AliAnalysisTaskKMeans(const AliAnalysisTaskKMeans &res)
169 : AliAnalysisTaskSE(res) 
170       ,fK(0)
171       ,fNMin(0)
172       ,fHists(0)
173       ,fH1CEta(0)
174       ,fH1CPhi(0)
175       ,fH1CEtaR(0)
176       ,fH2N1N2(0)
177       ,fH1Pt(0)
178       ,fH1PtC(0)
179       ,fH1PtC1(0)
180       ,fH1PtC2(0)
181       ,fH1PtAS(0)
182       ,fH1PtR(0)
183       ,fH1SPt(0)
184       ,fH1SPtC(0)
185       ,fH1DPhi(0)
186       ,fH1DR(0)
187       ,fH1DRR(0)
188       ,fH2DPhiEta(0)
189       ,fH2DPhiEtaR(0)
190       ,fH2DPhiEtaL(0)
191       ,fH2DPhiEtaLR(0)
192       ,fH2DPhiEtaC(0)
193       ,fH2DPhiEtaCR(0)
194       ,fH1Resp(0)
195       ,fH1RespR(0)
196       ,fH2Sigma(0)
197       ,fH2SigmaR(0)
198       ,fHDensity(0)
199       ,fHCSize(0)  
200       ,fHNCluster(0)
201       ,fHPtDensity(0)
202       ,fHDPhi(0)
203       ,fH2EtaPhi(0)
204       ,fH2REtaPhi(0)
205       ,fCuts(0)
206 {
207     // Dummy copy constructor
208     for (Int_t i=0; i< 10; i++) {
209         fA[i] = 0;
210         fB[i] = 0;
211     }
212 }
213
214 AliAnalysisTaskKMeans& AliAnalysisTaskKMeans::operator=(const AliAnalysisTaskKMeans& /*trclass*/)
215 {
216     // Dummy assignment operator
217     return *this;
218 }
219
220
221
222 //________________________________________________________________________
223 void AliAnalysisTaskKMeans::UserCreateOutputObjects()
224 {
225   // Create histograms
226   // Called once
227     fHists   = new TList();
228     fH1CEta  = new TH1F("fH1CEta",   "eta distribution of clusters",        90, -1.5, 1.5);
229     fH1CEtaR = new TH1F("fH1CEtaR",  "eta distribution of clusters",        90, -1.5, 1.5);
230     fH1CPhi  = new TH1F("fH1CPhi",   "phi distribution of clusters",       157,  0.0, 2. * TMath::Pi());
231     fH2N1N2  = new TH2F("fH2N1N2",   "multiplicity distribution",          50, 0., 50., 50, 0., 50.);
232     
233     fH1Pt    = new TH1F("fH1Pt",     "pt distribution",50, 0., 10.);
234     fH1PtC   = new TH1F("fH1PtC",    "pt distribution",50, 0., 10.);
235     fH1PtC1  = new TH1F("fH1PtC1",   "pt distribution",50, 0., 10.);
236     fH1PtC2  = new TH1F("fH1PtC2",   "pt distribution",50, 0., 10.);
237     fH1PtAS  = new TH1F("fH1PtAS",   "pt distribution",50, 0., 10.);
238     fH1PtR   = new TH1F("fH1PtR",    "pt distribution",50, 0., 10.);
239
240     fH1SPt    = new TH1F("fH1SPt",     "sum pt distribution",50, 0., 10.);
241     fH1SPtC   = new TH1F("fH1SPtC",    "sum pt distribution",50, 0., 10.);
242
243     fH1DR         = new TH1F("fH1DR",     "dR distribution", 50, 0., 5.);
244     fH1DPhi       = new TH1F("fH1DPhi",   "dPhi distribution", 31, 0., TMath::Pi());
245     fH2DPhiEta    = new TH2F("fH2DPhiEta","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
246     fH2DPhiEtaR   = new TH2F("fH2DPhiEtaR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
247     fH2DPhiEtaL   = new TH2F("fH2DPhiEtaL","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
248     fH2DPhiEtaLR  = new TH2F("fH2DPhiEtaLR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
249     fH2DPhiEtaC   = new TH2F("fH2DPhiEtaC","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
250     fH2DPhiEtaCR  = new TH2F("fH2DPhiEtaCR","eta phi distribution", 31, 0., TMath::Pi(), 20, 0., 2.);
251     fH1DRR        = new TH1F("fH1DRR",    "dR distribution", 50, 0., 5.);
252     fH1Resp       = new TH1F("fH1Resp",   "Responsibility", 50, 0., 1.);
253     fH1RespR      = new TH1F("fH1RespR",  "Responsibility", 50, 0., 1.);
254     fH2Sigma      = new TH2F("fH2Sigma",  "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.);
255     fH2SigmaR     = new TH2F("fH2SigmaR", "Sigma", 31, 0., TMath::Pi(), 20, 0., 2.);
256
257
258     fHDensity    = new TH1F("fHDensity",    "density distribution", 100, 0., 20.);
259     fHCSize      = new TH1F("fHCSize",      "cluster size", 20, -0.5, 19.5);
260
261     fHNCluster   = new TH1F("fHNCluster",   "Number of clusters", 11, -0.5, 10.5);
262     fHPtDensity  = new TH2F("fHPtDensity", "Pt vs density", 100, 0., 20., 50, 0., 10.);
263     fHDPhi       = new TH1F("fHDPhi",   "phi correlation", 100, 0., 0.5);
264     fH2EtaPhi    = new TH2F("fH2EtaPhi", "Eta Phi", 200, -1., 1., 628, 0., 2. * TMath::Pi());  
265     fHists->SetOwner();
266
267     fHists->Add(fH1CEta);
268     fHists->Add(fH1CEtaR);
269     fHists->Add(fH1CPhi);
270     fHists->Add(fH2N1N2);
271     fHists->Add(fH1Pt);
272     fHists->Add(fH1PtR);
273     fHists->Add(fH1PtC);
274     fHists->Add(fH1PtC1);
275     fHists->Add(fH1PtC2);
276     fHists->Add(fH1PtAS);
277     fHists->Add(fH1DR);
278     fHists->Add(fH1SPtC);
279     fHists->Add(fH1SPt);
280     fHists->Add(fH1DPhi);
281     fHists->Add(fH2DPhiEta);
282     fHists->Add(fH2DPhiEtaR);
283     fHists->Add(fH2DPhiEtaL);
284     fHists->Add(fH2DPhiEtaLR);
285     fHists->Add(fH2DPhiEtaC);
286     fHists->Add(fH2DPhiEtaCR);
287     fHists->Add(fH1DRR);
288     fHists->Add(fH1RespR);
289     fHists->Add(fH1Resp);    
290     fHists->Add(fH2Sigma);    
291     fHists->Add(fH2SigmaR);
292     fHists->Add(fHCSize);    
293     fHists->Add(fHDensity);
294     fHists->Add(fHNCluster);
295     fHists->Add(fHPtDensity);
296     fHists->Add(fHDPhi);
297     fHists->Add(fH2EtaPhi);
298     //
299     for (Int_t i = 0; i < 10; i++) {
300       fA[i] = new AliKMeansResult(i+1);
301       fB[i] = new AliKMeansResult(i+1);
302     }
303 }
304
305 //________________________________________________________________________
306 void AliAnalysisTaskKMeans::UserExec(Option_t *) 
307 {
308   // Main loop
309   // Called for each event
310     Double_t   phi [500] = {0};
311     Double_t   phiR[500] = {0};
312     Double_t    eta[500] = {0};
313     Double_t   etaR[500] = {0};
314     Double_t    pt [500] = {0};
315
316   if (!fInputEvent) {
317     Printf("ERROR: fESD not available");
318     return;
319   }
320   
321
322   Int_t ic = 0;
323   
324   //
325   // Fill eta-phi positions
326   //
327
328   Double_t ptmax = 0.;
329   Int_t    icmax = -1;
330   
331   for (Int_t iTracks = 0; iTracks < fInputEvent->GetNumberOfTracks(); iTracks++) {
332       AliVParticle* track = fInputEvent->GetTrack(iTracks);
333       if ((fCuts->AcceptTrack((AliESDtrack*)track))) 
334       {
335         const AliExternalTrackParam * tpcT = ((AliESDtrack*) track)->GetTPCInnerParam();
336         if (!tpcT) continue;
337         if (TMath::Abs(tpcT->Eta()) > 0.9) continue;
338
339           phi [ic] = tpcT->Phi();
340           eta [ic] = tpcT->Eta();
341           pt  [ic] = tpcT->Pt();
342
343           if (fH2REtaPhi) {
344             fH2REtaPhi->GetRandom2(etaR[ic], phiR[ic]);
345           } else {
346             phiR[ic] = 2. * TMath::Pi() * gRandom->Rndm();
347             etaR[ic] = 1.8 * gRandom->Rndm() - 0.9;
348           }
349           
350
351           if (pt[ic] > ptmax) {
352               ptmax = pt[ic];
353               icmax = ic;
354           }
355
356           fH2EtaPhi->Fill(eta[ic], phi[ic]);
357           ic++;
358       }
359   } //track loop 
360
361   for (Int_t i = 0; i < ic; i++) {
362     for (Int_t j = i+1; j < ic; j++) {
363       Double_t dphi = TMath::Abs(phi[i] - phi[j]);
364       fHDPhi->Fill(dphi);
365     }
366   }
367
368   //
369   // Cluster
370   if (ic < fNMin) {
371       PostData(1, fHists);
372       return;
373   }
374
375   //
376   Double_t rk0[10];
377   AliKMeansResult* res = 0;
378   AliKMeansResult best(10);
379   Float_t   rmaxG = -1.;
380
381   for (Int_t k = 0; k < 20; k++) {
382     Float_t   rmax   = -1.;
383     Int_t     imax   = 0;
384     for (Int_t i = 0; i < fK; i++) {
385       res = fA[i];
386       AliKMeansClustering::SoftKMeans2(i+1, ic, phi, eta, res->GetMx(),res->GetMy(), res->GetSigma2(), res->GetRk());
387       res->Sort(ic, phi, eta);
388       Int_t j = (res->GetInd())[0];
389       rk0[i]  = (res->GetTarget())[j];
390       if (rk0[i] > rmax)  {
391         rmax = rk0[i];
392         imax = i;
393       }
394     }
395     if (rmax > rmaxG) {
396       rmaxG = rmax;
397       best.CopyResults(fA[imax]);
398     }
399   }
400   Double_t* mPhi     = best.GetMx();
401   Double_t* mEta     = best.GetMy();
402   Double_t* sigma2   = best.GetSigma2();
403   Int_t     nk       = best.GetK();
404   Int_t     im       = (best.GetInd())[0];
405   Double_t  etaC     = mEta[im];
406
407   fHDensity->Fill(rmaxG / TMath::Pi());
408   fHCSize->Fill(2.28 * rmaxG * sigma2[im]);
409   fHNCluster->Fill(Float_t(nk));
410
411   Double_t dphic, detac;
412
413   if (rmaxG > 0. && TMath::Abs(etaC) < 0.4) {
414     // Analysis
415     //
416     // Cluster Multiplicity
417     Int_t mult[2] = {0, 0};
418     //
419     if (nk > 1) {
420       dphic = DeltaPhi(mPhi[0], mPhi[1]);
421       detac = TMath::Abs(mEta[0] - mEta[1]);
422       fH2DPhiEtaC->Fill(dphic, detac);  
423     }
424     //
425     // Random cluster position
426     Int_t ir = Int_t(Float_t(ic) * gRandom->Rndm());
427
428     Double_t crPhi = phi[ir];
429     Double_t crEta = eta[ir];
430     //
431     Double_t sumPt  = 0;
432     Double_t sumPtC = 0;
433
434     for (Int_t i = 0; i < 1; i++) {
435       fH1CEta->Fill(mEta[im]);
436       fH1CPhi->Fill(mPhi[im]);      
437       for (Int_t j = 0; j < ic; j++) {
438         Double_t r    = DeltaR(mPhi[im], mEta[im], phi[j], eta[j]);
439         Double_t dphi = DeltaPhi(mPhi[im], phi[j]);
440         Double_t deta = mEta[im] - eta[j];
441         Double_t rr   = DeltaR(crPhi, crEta, phi[j], eta[j]);
442
443         fH1DR->Fill(r);
444         fH1DPhi->Fill(dphi);
445         fH2DPhiEta->Fill(dphi, TMath::Abs(deta));
446         if (j == icmax) fH2DPhiEtaL->Fill(dphi, TMath::Abs(deta));
447         
448         if (r < 0.2) {
449           fH1PtC2->Fill(pt[j]);
450         }
451         if (r < 0.3) {
452           fH1PtC1->Fill(pt[j]);
453           fHPtDensity->Fill(rmaxG/TMath::Pi(), pt[j]);
454         }
455         if (rr < 0.3) {
456             fH1PtR->Fill(pt[j]);
457         }
458
459         if (r < 0.4)  {
460           sumPtC += pt[j];
461           mult[i]++;
462           fH1PtC->Fill(pt[j]);
463         } 
464         if (r > 0.7 && dphi < (TMath::Pi() - 0.3)) {
465           fH1Pt->Fill(pt[j]);
466         }
467         
468         if (r > 0.7 && r < 1.) {
469           sumPt += pt[j];
470         }
471         
472         if (dphi > (TMath::Pi() - 0.3)) {
473           fH1PtAS->Fill(pt[j], 1.);
474         }
475       }
476     }
477
478     fH2N1N2->Fill(Float_t(mult[0]), Float_t(mult[1]));
479     fH1SPt ->Fill(sumPt);
480     fH1SPtC->Fill(sumPtC);
481   }
482     //
483     // Randomized phi
484     //
485   rmaxG = -1.;
486   for (Int_t k = 0; k < 20; k++) {
487     Float_t rmax   = -1.;
488     Int_t   imax   =  0;
489     for (Int_t i = 0; i < fK; i++) {
490       res = fB[i];
491       AliKMeansClustering::SoftKMeans2(i+1, ic, phiR, etaR, res->GetMx(),res->GetMy(), res->GetSigma2(), res->GetRk());
492       res->Sort(ic, phiR, etaR);
493       Int_t j = (res->GetInd())[0];
494       rk0[i]  = (res->GetTarget())[j];
495       if (rk0[i] > rmax) {
496         rmax = rk0[i];
497         imax = i;
498       }
499     }
500     if (rmax > rmaxG) {
501       rmaxG = rmax;
502       best.CopyResults(fB[imax]);
503     }    
504   }
505   
506     mPhi    = best.GetMx();
507     mEta    = best.GetMy();
508     nk      = best.GetK();
509     im      = (best.GetInd())[0];
510     etaC    = mEta[im];    
511
512     //
513     // Cluster Multiplicity
514     if (rmaxG > 0. && TMath::Abs(etaC) < 0.4) {
515       for (Int_t i = 0; i < 1; i++) {
516         im = (best.GetInd())[i];
517         fH1CEtaR->Fill(mEta[im]);
518       }
519       
520       for (Int_t i = 0; i < 1; i++) {
521         im = (best.GetInd())[i];
522         for (Int_t j = 0; j < ic; j++) {
523           Double_t dphi = DeltaPhi(mPhi[im], phiR[j]);
524           Double_t deta = mEta[im] - etaR[j];
525           Double_t r = DeltaR(mPhi[im], mEta[im], phiR[j], etaR[j]);
526           fH1DRR->Fill(r);
527           fH2DPhiEtaR->Fill(dphi, TMath::Abs(deta));
528           if (j == icmax) fH2DPhiEtaLR->Fill(dphi, TMath::Abs(deta));
529         }
530       }
531       if (nk > 1) {
532         dphic = DeltaPhi(mPhi[0], mPhi[1]);
533         detac = TMath::Abs(mEta[0] - mEta[1]);
534         fH2DPhiEtaCR->Fill(dphic, detac);  
535       }
536     }
537     PostData(1, fHists);
538 }      
539
540 Double_t AliAnalysisTaskKMeans::DeltaPhi(Double_t phi1, Double_t phi2)
541 {
542     Double_t dphi = TMath::Abs(phi1 - phi2);
543     if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
544     return dphi;
545 }
546
547 Double_t AliAnalysisTaskKMeans::DeltaR(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
548 {
549     Double_t dphi = DeltaPhi(phi1, phi2);
550     Double_t deta = eta1 - eta2;
551     return (TMath::Sqrt(dphi * dphi + deta * deta));
552     
553 }
554
555 //________________________________________________________________________
556 void AliAnalysisTaskKMeans::Terminate(Option_t *) 
557 {
558 }  
559