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"
15
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
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//________________________________________________________________________
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();
69}
70
71
72//________________________________________________________________________
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);
139a69f1 84
85 if (!aodE) {
86 Printf("ERROR: AOD not available");
87 return;
88 }
89
0eba9d77 90 TClonesArray* jets = dynamic_cast<TClonesArray*> (aodE->FindListObject("jetsAOD_FASTKT04"));
139a69f1 91 if (!jets) {
92 Printf("ERROR: Jet branch not available");
93 return;
94 }
95
0eba9d77 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));
781ef0e6 106 if (!jet) continue;
0eba9d77 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 = aodE->GetTrack(i);
158//
159// Track quality cuts
160 if (!track->TestFilterBit(1<<4)) continue;
161 fPtTH->Fill(track->Pt());
162
163//
164//
165 TVector3 pp(track->Px(), track->Py(), track->Pz());
166
167 Float_t phi = track->Phi();
168 Float_t eta = track->Eta();
169 Float_t pt = track->Pt();
170 Float_t jT = pp.Perp(ppJ1);
171
172
173 // pT > 2.
174
175
176 Float_t deta = eta - etaJ;
177 Float_t dphi = phi - phiJ;
178
179 if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
180 if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
181 Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
182
183 if (r < 0.4) ptsum += pt;
184 if (r < 0.5 && jT > 1.) {
185 TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
186 TVector3 pPerp = pp - pLong;
187
188 Float_t ppjX = pPerp.Dot(ppJ2);
189 Float_t ppjY = pPerp.Dot(ppJ3);
190 Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
191
192 mxx += (ppjX * ppjX / ppjT);
193 myy += (ppjY * ppjY / ppjT);
194 mxy += (ppjX * ppjY / ppjT);
195 nc++;
196 sump2 += ppjT;
197 // max pt
198 if (pt > ptMax) {
199 ptMax = pt;
200 iMax = i;
201 etaMax = deta;
202 phiMax = dphi;
203 }
204 } // R < 0.4
205 } // 1st Track Loop
206 fPtSum->Fill(ptJ, ptsum);
207
208//
209// At this point we have mxx, myy, mxy
210 if (nc == 0) return;
211// Shericity Matrix
212 const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
213 TMatrixDSym m0(2, ele);
214// Find eigenvectors
215 TMatrixDSymEigen m(m0);
216 TVectorD eval(2);
217 TMatrixD evecm = m.GetEigenVectors();
218 eval = m.GetEigenValues();
219// Largest eigenvector
220 Int_t jev = 0;
221 if (eval[0] < eval[1]) jev = 1;
222 TVectorD evec0(2);
223// Principle axis
224 evec0 = TMatrixDColumn(evecm, jev);
225 TVector2 evec(evec0[0], evec0[1]);
226// Principle axis from leading partice
227 Float_t phiM = TMath::ATan2(phiMax, etaMax);
228 TVector2 evecM(TMath::Cos(phiM), TMath::Sin(phiM));
229 Float_t phistM = evecM.DeltaPhi(evec);
230 if (TMath::Abs(phistM) > TMath::Pi()/2.) evec*=(-1.);
231
232//
233// 2nd correlation with principal axis
234//
235 for (Int_t i = 0; i < nT; i++) {
236 AliAODTrack* track = aodE->GetTrack(i);
237//
238// Track quality cuts
239 if (!track->TestFilterBit(1<<4)) continue;
240//
241//
242 TVector3 pp(track->Px(), track->Py(), track->Pz());
243 Float_t phi = track->Phi();
244 Float_t eta = track->Eta();
245 Float_t pt = track->Pt();
246 Float_t jT = pp.Perp(ppJ1);
247 // pT < 2.
248 if (jT > 1.) continue;
249 Float_t deta = eta - etaJ;
250 Float_t dphi = phi - phiJ;
251 if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi;
252 if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
253
254 Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
255
256 if (r > 0.7) continue;
257
258 TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
259 TVector3 pPerp = pp - pLong;
260
261 Float_t ppjX = pPerp.Dot(ppJ2);
262 Float_t ppjY = pPerp.Dot(ppJ3);
0eba9d77 263
264 TVector2 vr(ppjX, ppjY) ;
265 Float_t phistr = evec.DeltaPhi(vr);
266 fPhiM->Fill(phistr);
267 fPhiMPt ->Fill(phistr, pt);
268 fPhiMPtJ->Fill(phistr, ptJ);
269
270 } // 2nd Track loop
271
272
273 // Post output data.
274 PostData(1, fHists);
275}
276
277//________________________________________________________________________