Adapt add macro and particle/cluster containers for the analysis on jets
[u/mrichter/AliRoot.git] / PWGJE / 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);
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++) {
f15c1f69 157 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aodE->GetTrack(i));
158 if(!track) AliFatal("Not a standard AOD");
0eba9d77 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++) {
f15c1f69 237 AliAODTrack* track = dynamic_cast<AliAODTrack*>(aodE->GetTrack(i));
238 if(!track) AliFatal("Not a standard AOD");
0eba9d77 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);
0eba9d77 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//________________________________________________________________________
280void AliAnalysisTaskJetsTM::Terminate(Option_t *)
281{
282// Terminate
283}