//
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetsTM.cxx
CommitLineData
0eba9d77 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
16ClassImp(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
32AliAnalysisTaskJetsTM::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//________________________________________________________________________
50void 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//________________________________________________________________________
73void 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 TClonesArray* jets = dynamic_cast<TClonesArray*> (aodE->FindListObject("jetsAOD_FASTKT04"));
85 Int_t nJ = jets->GetEntries();
86
87 Float_t ptmax = 0.;
88 Int_t imax = -1;
89
90//
91// Find highest pT jet with pt > 20 GeV
92//
93 for (Int_t i = 0; i < nJ; i++) {
94 AliAODJet* jet = dynamic_cast<AliAODJet*> (jets->At(i));
95 Float_t ptJ = jet->Pt();
96 Float_t etaJ = TMath::Abs(jet->Eta());
97 if ((ptJ > 20.) && (ptJ > ptmax) && etaJ < 0.5) {
98 ptmax = ptJ;
99 imax = i;
100 }
101 }
102
103 if (imax == -1) return;
104
105
106 AliAODJet* jet = dynamic_cast<AliAODJet*> (jets->At(imax));
107 Float_t ptJ = jet->Pt();
108 Float_t phiJ = jet->Phi();
109 Float_t etaJ = jet->Eta();
110
111 fPtH->Fill(ptJ);
112
113//
114// The transverse plane
115//
116// 2 normalized vectors perpendicular to the jet direction
117//
118 Float_t pxJ = jet->Px();
119 Float_t pyJ = jet->Py();
120 Float_t pzJ = jet->Pz();
121
122 TVector3 ppJ1(pxJ, pyJ, pzJ);
123 TVector3 ppJ3(- pxJ * pzJ, - pyJ * pzJ, pxJ * pxJ + pyJ * pyJ);
124 ppJ3.SetMag(1.);
125 TVector3 ppJ2(-pyJ, pxJ, 0);
126 ppJ2.SetMag(1.);
127
128//
129// 1st track loop to determine the sphericity matrix
130//
131 Int_t nT = aodE->GetNumberOfTracks();
132 // Initialize Shericity matrix
133 Float_t mxx = 0.;
134 Float_t myy = 0.;
135 Float_t mxy = 0.;
136 Int_t nc = 0;
137 Float_t sump2 = 0.;
138 Float_t ptMax = 0.;
139 Float_t etaMax = 0.;
140 Float_t phiMax = 0.;
141 Int_t iMax = -1;
142 Float_t ptsum = 0.;
143
144 for (Int_t i = 0; i < nT; i++) {
145 AliAODTrack* track = aodE->GetTrack(i);
146//
147// Track quality cuts
148 if (!track->TestFilterBit(1<<4)) continue;
149 fPtTH->Fill(track->Pt());
150
151//
152//
153 TVector3 pp(track->Px(), track->Py(), track->Pz());
154
155 Float_t phi = track->Phi();
156 Float_t eta = track->Eta();
157 Float_t pt = track->Pt();
158 Float_t jT = pp.Perp(ppJ1);
159
160
161 // pT > 2.
162
163
164 Float_t deta = eta - etaJ;
165 Float_t dphi = phi - phiJ;
166
167 if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
168 if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
169 Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
170
171 if (r < 0.4) ptsum += pt;
172 if (r < 0.5 && jT > 1.) {
173 TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
174 TVector3 pPerp = pp - pLong;
175
176 Float_t ppjX = pPerp.Dot(ppJ2);
177 Float_t ppjY = pPerp.Dot(ppJ3);
178 Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
179
180 mxx += (ppjX * ppjX / ppjT);
181 myy += (ppjY * ppjY / ppjT);
182 mxy += (ppjX * ppjY / ppjT);
183 nc++;
184 sump2 += ppjT;
185 // max pt
186 if (pt > ptMax) {
187 ptMax = pt;
188 iMax = i;
189 etaMax = deta;
190 phiMax = dphi;
191 }
192 } // R < 0.4
193 } // 1st Track Loop
194 fPtSum->Fill(ptJ, ptsum);
195
196//
197// At this point we have mxx, myy, mxy
198 if (nc == 0) return;
199// Shericity Matrix
200 const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
201 TMatrixDSym m0(2, ele);
202// Find eigenvectors
203 TMatrixDSymEigen m(m0);
204 TVectorD eval(2);
205 TMatrixD evecm = m.GetEigenVectors();
206 eval = m.GetEigenValues();
207// Largest eigenvector
208 Int_t jev = 0;
209 if (eval[0] < eval[1]) jev = 1;
210 TVectorD evec0(2);
211// Principle axis
212 evec0 = TMatrixDColumn(evecm, jev);
213 TVector2 evec(evec0[0], evec0[1]);
214// Principle axis from leading partice
215 Float_t phiM = TMath::ATan2(phiMax, etaMax);
216 TVector2 evecM(TMath::Cos(phiM), TMath::Sin(phiM));
217 Float_t phistM = evecM.DeltaPhi(evec);
218 if (TMath::Abs(phistM) > TMath::Pi()/2.) evec*=(-1.);
219
220//
221// 2nd correlation with principal axis
222//
223 for (Int_t i = 0; i < nT; i++) {
224 AliAODTrack* track = aodE->GetTrack(i);
225//
226// Track quality cuts
227 if (!track->TestFilterBit(1<<4)) continue;
228//
229//
230 TVector3 pp(track->Px(), track->Py(), track->Pz());
231 Float_t phi = track->Phi();
232 Float_t eta = track->Eta();
233 Float_t pt = track->Pt();
234 Float_t jT = pp.Perp(ppJ1);
235 // pT < 2.
236 if (jT > 1.) continue;
237 Float_t deta = eta - etaJ;
238 Float_t dphi = phi - phiJ;
239 if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
240 if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
241
242 Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
243
244 if (r > 0.7) continue;
245
246 TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
247 TVector3 pPerp = pp - pLong;
248
249 Float_t ppjX = pPerp.Dot(ppJ2);
250 Float_t ppjY = pPerp.Dot(ppJ3);
251 Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
252
253 TVector2 vr(ppjX, ppjY) ;
254 Float_t phistr = evec.DeltaPhi(vr);
255 fPhiM->Fill(phistr);
256 fPhiMPt ->Fill(phistr, pt);
257 fPhiMPtJ->Fill(phistr, ptJ);
258
259 } // 2nd Track loop
260
261
262 // Post output data.
263 PostData(1, fHists);
264}
265
266//________________________________________________________________________
267void AliAnalysisTaskJetsTM::Terminate(Option_t *)
268{
269// Terminate
270}