hjet mass ana update + bug fix
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetsTM.cxx
1 #include "TH1F.h"
2 #include "TH2F.h"
3 #include "TProfile.h"
4 #include "TClonesArray.h"
5 #include "TCanvas.h"
6 #include "TList.h"
7 #include "TMatrixD.h"
8 #include "TMatrixDSym.h"
9 #include "TMatrixDSymEigen.h"
10
11 #include "AliAODEvent.h"
12 #include "AliAODTrack.h"
13 #include "AliAODJet.h"
14 #include "AliAnalysisTaskJetsTM.h"
15
16 ClassImp(AliAnalysisTaskJetsTM)
17 //
18 // Thrust Major (TM) analysis of reconstructed jets.
19 // TM is the thrust in the plane perpendicular to the jet axis
20 // The present amalysis performs the following steps:
21 // (a) Construct to orthogonal unit vectors (e1, e2) in the plane perpendicular to the jet axis
22 // (b) Calculate the components of all particles with jT > 1 GeV with respect to e1, e2
23 // (c) Construct the sphericity matrix
24 // (d) Find the two orthogonal eigenvectors of the spericity matrix
25 // (e) Caluclate the components of all particles with jT < 1 GeV in the reference frame spanned by the eigenvectors
26 // (f) Calculate the azimuthal angle in this frame
27 //
28 //
29 // Author: andreas.morsch@cern.ch
30
31
32 AliAnalysisTaskJetsTM::AliAnalysisTaskJetsTM(const char *name) 
33     : AliAnalysisTaskSE(name) 
34       ,fHists(0)
35       ,fPtH(0)
36       ,fPtTH(0)
37       ,fPhiM(0)
38       ,fPhiMPt(0)
39       ,fPhiMPtJ(0)
40       ,fPtSum(0)
41 {
42   //
43   // Constructor
44   //
45   DefineOutput(1,  TList::Class());
46 }
47
48
49 //________________________________________________________________________
50 void AliAnalysisTaskJetsTM::UserCreateOutputObjects()
51 {
52   // Create histograms
53   // Called once
54     fHists    = new TList();
55     fPtH      = new TH1F("fPtH" ,     " p_{T}       distribution of jets"  ,            50, 0., 100.0);
56     fPtSum    = new TH2F("fPtSum" ,   " p_{T} sum in cone"                 ,            50, 0., 100.0, 50, 0, 100.);
57     fPtTH     = new TH1F("fPtTH",     " p_{T}       distribution of tracks",            50, 0., 100.0);
58     fPhiM     = new TH1F("fPhiM",     " phi^{major} distribution of tracks",            64, 0.,   3.2);
59     fPhiMPt   = new TH2F("fPhiMPt",   " phi^{major} distribution of tracks vs pt",      64, 0.,   3.2, 20, 0., 2.);
60     fPhiMPtJ  = new TH2F("fPhiMPtJ",  " phi^{major} distribution of tracks vs pt Jet",  64, 0.,   3.2, 40, 0., 100.);
61
62     fHists->SetOwner();
63     fHists->Add(fPtH);
64     fHists->Add(fPtSum);
65     fHists->Add(fPtTH);
66     fHists->Add(fPhiM);
67     fHists->Add(fPhiMPt);
68     fHists->Add(fPhiMPtJ);
69 }
70
71
72 //________________________________________________________________________
73 void AliAnalysisTaskJetsTM::UserExec(Option_t *) 
74 {
75   // Main loop
76   // Called for each event
77
78   if (!fInputEvent) {
79     Printf("ERROR: AOD not available");
80     return;
81   }
82
83   AliAODEvent* aodE  = dynamic_cast<AliAODEvent*>  (fInputEvent);
84
85   if (!aodE) {
86     Printf("ERROR: AOD not available");
87     return;
88   }
89
90   TClonesArray* jets = dynamic_cast<TClonesArray*> (aodE->FindListObject("jetsAOD_FASTKT04"));
91   if (!jets) {
92     Printf("ERROR: Jet branch not available");
93     return;
94   }
95
96   Int_t nJ = jets->GetEntries();
97
98   Float_t ptmax = 0.;
99   Int_t   imax  = -1;
100   
101 //
102 // Find highest pT jet with pt > 20 GeV
103 //
104   for (Int_t i = 0; i < nJ; i++) {
105       AliAODJet* jet = dynamic_cast<AliAODJet*> (jets->At(i));
106       if (!jet) continue;
107       Float_t ptJ  = jet->Pt();
108       Float_t etaJ = TMath::Abs(jet->Eta());
109       if ((ptJ > 20.) && (ptJ  > ptmax) && etaJ < 0.5) {
110           ptmax = ptJ;
111           imax = i;
112       }
113   }
114
115   if (imax == -1) return;
116   
117
118   AliAODJet* jet = dynamic_cast<AliAODJet*> (jets->At(imax));
119   Float_t ptJ  = jet->Pt();
120   Float_t phiJ = jet->Phi();
121   Float_t etaJ = jet->Eta();
122   
123   fPtH->Fill(ptJ);
124
125 //
126 // The transverse plane
127 //
128 // 2 normalized vectors perpendicular to the jet direction
129 //
130       Float_t pxJ = jet->Px();
131       Float_t pyJ = jet->Py();
132       Float_t pzJ = jet->Pz();
133       
134       TVector3  ppJ1(pxJ, pyJ, pzJ);
135       TVector3  ppJ3(- pxJ * pzJ, - pyJ * pzJ, pxJ * pxJ + pyJ * pyJ);
136       ppJ3.SetMag(1.);
137       TVector3  ppJ2(-pyJ, pxJ, 0);
138       ppJ2.SetMag(1.);
139
140 //
141 // 1st track loop to determine the sphericity matrix
142 //   
143       Int_t nT = aodE->GetNumberOfTracks();
144       // Initialize Shericity matrix
145       Float_t mxx    = 0.;
146       Float_t myy    = 0.;
147       Float_t mxy    = 0.;
148       Int_t   nc     = 0;
149       Float_t sump2  = 0.;
150       Float_t ptMax  = 0.;
151       Float_t etaMax = 0.;
152       Float_t phiMax = 0.;
153       Int_t   iMax   = -1;
154       Float_t ptsum  = 0.;
155       
156       for (Int_t i = 0; i < nT; i++) {
157           AliAODTrack* track = dynamic_cast<AliAODTrack*>(aodE->GetTrack(i));
158           if(!track) AliFatal("Not a standard AOD");
159 //
160 //    Track quality cuts
161           if (!track->TestFilterBit(1<<4)) continue;
162           fPtTH->Fill(track->Pt());
163
164 //
165 //
166           TVector3 pp(track->Px(), track->Py(), track->Pz());
167               
168           Float_t phi = track->Phi();
169           Float_t eta = track->Eta();
170           Float_t pt  = track->Pt();
171           Float_t jT  = pp.Perp(ppJ1);
172           
173           
174           // pT > 2.
175           
176             
177           Float_t deta = eta - etaJ;
178           Float_t dphi = phi - phiJ;
179
180           if (dphi >   TMath::Pi()) dphi =   2. * TMath::Pi() - dphi;
181           if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
182           Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
183
184           if (r < 0.4) ptsum += pt;
185           if (r < 0.5 && jT > 1.) {
186               TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
187               TVector3 pPerp = pp - pLong;
188
189               Float_t ppjX = pPerp.Dot(ppJ2);
190               Float_t ppjY = pPerp.Dot(ppJ3);
191               Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
192
193               mxx += (ppjX * ppjX / ppjT);
194               myy += (ppjY * ppjY / ppjT);
195               mxy += (ppjX * ppjY / ppjT);
196               nc++;
197               sump2 += ppjT;
198               // max pt
199               if (pt > ptMax) {
200                   ptMax  = pt;
201                   iMax   = i;
202                   etaMax = deta;
203                   phiMax = dphi;
204               }
205           } // R < 0.4
206       } // 1st Track Loop
207       fPtSum->Fill(ptJ, ptsum);
208       
209 //
210 // At this point we have mxx, myy, mxy
211       if (nc == 0) return;      
212 // Shericity Matrix     
213       const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};     
214       TMatrixDSym m0(2, ele);
215 // Find eigenvectors
216       TMatrixDSymEigen m(m0);
217       TVectorD eval(2);
218       TMatrixD evecm = m.GetEigenVectors();
219       eval  = m.GetEigenValues();
220 // Largest eigenvector
221       Int_t jev = 0;
222       if (eval[0] < eval[1]) jev = 1;
223       TVectorD evec0(2);
224 // Principle axis
225       evec0 = TMatrixDColumn(evecm, jev);
226       TVector2 evec(evec0[0], evec0[1]); 
227 // Principle axis from leading partice
228       Float_t phiM = TMath::ATan2(phiMax, etaMax);
229       TVector2 evecM(TMath::Cos(phiM), TMath::Sin(phiM)); 
230       Float_t phistM = evecM.DeltaPhi(evec);
231       if (TMath::Abs(phistM) > TMath::Pi()/2.) evec*=(-1.);
232
233 //
234 // 2nd correlation with principal axis
235 //   
236       for (Int_t i = 0; i < nT; i++) {
237           AliAODTrack* track = dynamic_cast<AliAODTrack*>(aodE->GetTrack(i));
238           if(!track) AliFatal("Not a standard AOD");
239 //
240 //    Track quality cuts
241           if (!track->TestFilterBit(1<<4)) continue;
242 //
243 //
244           TVector3 pp(track->Px(), track->Py(), track->Pz());
245           Float_t phi = track->Phi();
246           Float_t eta = track->Eta();
247           Float_t pt  = track->Pt();
248           Float_t jT  = pp.Perp(ppJ1);
249           // pT < 2.
250           if (jT > 1.) continue;            
251           Float_t deta = eta - etaJ;
252           Float_t dphi = phi - phiJ;
253           if (dphi >   TMath::Pi()) dphi =   2. * TMath::Pi() - dphi;
254           if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
255
256           Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
257
258           if (r > 0.7) continue;
259
260           TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
261           TVector3 pPerp = pp - pLong;
262           
263           Float_t ppjX = pPerp.Dot(ppJ2);
264           Float_t ppjY = pPerp.Dot(ppJ3);
265             
266           TVector2 vr(ppjX, ppjY) ;
267           Float_t phistr = evec.DeltaPhi(vr);
268           fPhiM->Fill(phistr);
269           fPhiMPt ->Fill(phistr, pt);
270           fPhiMPtJ->Fill(phistr, ptJ);
271           
272       } // 2nd Track loop
273       
274   
275   // Post output data.
276   PostData(1, fHists);
277 }      
278
279 //________________________________________________________________________
280 void AliAnalysisTaskJetsTM::Terminate(Option_t *) 
281 {
282 // Terminate
283 }