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