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