4 #include "TClonesArray.h"
8 #include "TMatrixDSym.h"
9 #include "TMatrixDSymEigen.h"
11 #include "AliAODEvent.h"
12 #include "AliAODTrack.h"
13 #include "AliAODJet.h"
14 #include "AliAnalysisTaskJetsTM.h"
16 ClassImp(AliAnalysisTaskJetsTM)
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
29 // Author: andreas.morsch@cern.ch
32 AliAnalysisTaskJetsTM::AliAnalysisTaskJetsTM(const char *name)
33 : AliAnalysisTaskSE(name)
45 DefineOutput(1, TList::Class());
49 //________________________________________________________________________
50 void AliAnalysisTaskJetsTM::UserCreateOutputObjects()
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.);
68 fHists->Add(fPhiMPtJ);
72 //________________________________________________________________________
73 void AliAnalysisTaskJetsTM::UserExec(Option_t *)
76 // Called for each event
79 Printf("ERROR: AOD not available");
83 AliAODEvent* aodE = dynamic_cast<AliAODEvent*> (fInputEvent);
86 Printf("ERROR: AOD not available");
90 TClonesArray* jets = dynamic_cast<TClonesArray*> (aodE->FindListObject("jetsAOD_FASTKT04"));
92 Printf("ERROR: Jet branch not available");
96 Int_t nJ = jets->GetEntries();
102 // Find highest pT jet with pt > 20 GeV
104 for (Int_t i = 0; i < nJ; i++) {
105 AliAODJet* jet = dynamic_cast<AliAODJet*> (jets->At(i));
107 Float_t ptJ = jet->Pt();
108 Float_t etaJ = TMath::Abs(jet->Eta());
109 if ((ptJ > 20.) && (ptJ > ptmax) && etaJ < 0.5) {
115 if (imax == -1) return;
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();
126 // The transverse plane
128 // 2 normalized vectors perpendicular to the jet direction
130 Float_t pxJ = jet->Px();
131 Float_t pyJ = jet->Py();
132 Float_t pzJ = jet->Pz();
134 TVector3 ppJ1(pxJ, pyJ, pzJ);
135 TVector3 ppJ3(- pxJ * pzJ, - pyJ * pzJ, pxJ * pxJ + pyJ * pyJ);
137 TVector3 ppJ2(-pyJ, pxJ, 0);
141 // 1st track loop to determine the sphericity matrix
143 Int_t nT = aodE->GetNumberOfTracks();
144 // Initialize Shericity matrix
156 for (Int_t i = 0; i < nT; i++) {
157 AliAODTrack* track = aodE->GetTrack(i);
159 // Track quality cuts
160 if (!track->TestFilterBit(1<<4)) continue;
161 fPtTH->Fill(track->Pt());
165 TVector3 pp(track->Px(), track->Py(), track->Pz());
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);
176 Float_t deta = eta - etaJ;
177 Float_t dphi = phi - phiJ;
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);
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;
188 Float_t ppjX = pPerp.Dot(ppJ2);
189 Float_t ppjY = pPerp.Dot(ppJ3);
190 Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
192 mxx += (ppjX * ppjX / ppjT);
193 myy += (ppjY * ppjY / ppjT);
194 mxy += (ppjX * ppjY / ppjT);
206 fPtSum->Fill(ptJ, ptsum);
209 // At this point we have mxx, myy, mxy
212 const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
213 TMatrixDSym m0(2, ele);
215 TMatrixDSymEigen m(m0);
217 TMatrixD evecm = m.GetEigenVectors();
218 eval = m.GetEigenValues();
219 // Largest eigenvector
221 if (eval[0] < eval[1]) jev = 1;
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.);
233 // 2nd correlation with principal axis
235 for (Int_t i = 0; i < nT; i++) {
236 AliAODTrack* track = aodE->GetTrack(i);
238 // Track quality cuts
239 if (!track->TestFilterBit(1<<4)) continue;
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);
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;
254 Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
256 if (r > 0.7) continue;
258 TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
259 TVector3 pPerp = pp - pLong;
261 Float_t ppjX = pPerp.Dot(ppJ2);
262 Float_t ppjY = pPerp.Dot(ppJ3);
264 TVector2 vr(ppjX, ppjY) ;
265 Float_t phistr = evec.DeltaPhi(vr);
267 fPhiMPt ->Fill(phistr, pt);
268 fPhiMPtJ->Fill(phistr, ptJ);
277 //________________________________________________________________________
278 void AliAnalysisTaskJetsTM::Terminate(Option_t *)