PWGJE/AliAnalysisTaskJetsTM.cxx

   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);
  84
  85   if (!aodE) {
  86     Printf("ERROR: AOD not available");
  87     return;
  88   }
  89
  90   TClonesArray* jets = dynamic_cast<TClonesArray*> (aodE->FindListObject("jetsAOD_FASTKT04"));
  91   if (!jets) {
  92     Printf("ERROR: Jet branch not available");
  93     return;
  94   }
  95
  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));
 106       if (!jet) continue;
 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);
 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 }