]>
Commit | Line | Data |
---|---|---|
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 | ||
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); | |
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 | //________________________________________________________________________ | |
278 | void AliAnalysisTaskJetsTM::Terminate(Option_t *) | |
279 | { | |
280 | // Terminate | |
281 | } |