PWG4/JetTasks/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   TClonesArray* jets = dynamic_cast<TClonesArray*> (aodE->FindListObject("jetsAOD_FASTKT04"));
  85   Int_t nJ = jets->GetEntries();
  86
  87   Float_t ptmax = 0.;
  88   Int_t   imax  = -1;
  89
  90 //
  91 // Find highest pT jet with pt > 20 GeV
  92 //
  93   for (Int_t i = 0; i < nJ; i++) {
  94       AliAODJet* jet = dynamic_cast<AliAODJet*> (jets->At(i));
  95       Float_t ptJ  = jet->Pt();
  96       Float_t etaJ = TMath::Abs(jet->Eta());
  97       if ((ptJ > 20.) && (ptJ  > ptmax) && etaJ < 0.5) {
  98           ptmax = ptJ;
  99           imax = i;
 100       }
 101   }
 102
 103   if (imax == -1) return;
 104
 105
 106   AliAODJet* jet = dynamic_cast<AliAODJet*> (jets->At(imax));
 107   Float_t ptJ  = jet->Pt();
 108   Float_t phiJ = jet->Phi();
 109   Float_t etaJ = jet->Eta();
 110
 111   fPtH->Fill(ptJ);
 112
 113 //
 114 // The transverse plane
 115 //
 116 // 2 normalized vectors perpendicular to the jet direction
 117 //
 118       Float_t pxJ = jet->Px();
 119       Float_t pyJ = jet->Py();
 120       Float_t pzJ = jet->Pz();
 121
 122       TVector3  ppJ1(pxJ, pyJ, pzJ);
 123       TVector3  ppJ3(- pxJ * pzJ, - pyJ * pzJ, pxJ * pxJ + pyJ * pyJ);
 124       ppJ3.SetMag(1.);
 125       TVector3  ppJ2(-pyJ, pxJ, 0);
 126       ppJ2.SetMag(1.);
 127
 128 //
 129 // 1st track loop to determine the sphericity matrix
 130 //
 131       Int_t nT = aodE->GetNumberOfTracks();
 132       // Initialize Shericity matrix
 133       Float_t mxx    = 0.;
 134       Float_t myy    = 0.;
 135       Float_t mxy    = 0.;
 136       Int_t   nc     = 0;
 137       Float_t sump2  = 0.;
 138       Float_t ptMax  = 0.;
 139       Float_t etaMax = 0.;
 140       Float_t phiMax = 0.;
 141       Int_t   iMax   = -1;
 142       Float_t ptsum  = 0.;
 143
 144       for (Int_t i = 0; i < nT; i++) {
 145           AliAODTrack* track = aodE->GetTrack(i);
 146 //
 147 //    Track quality cuts
 148           if (!track->TestFilterBit(1<<4)) continue;
 149           fPtTH->Fill(track->Pt());
 150
 151 //
 152 //
 153           TVector3 pp(track->Px(), track->Py(), track->Pz());
 154
 155           Float_t phi = track->Phi();
 156           Float_t eta = track->Eta();
 157           Float_t pt  = track->Pt();
 158           Float_t jT  = pp.Perp(ppJ1);
 159
 160
 161           // pT > 2.
 162
 163
 164           Float_t deta = eta - etaJ;
 165           Float_t dphi = phi - phiJ;
 166
 167           if (dphi >   TMath::Pi()) dphi =   2. * TMath::Pi() - dphi;
 168           if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
 169           Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
 170
 171           if (r < 0.4) ptsum += pt;
 172           if (r < 0.5 && jT > 1.) {
 173               TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
 174               TVector3 pPerp = pp - pLong;
 175
 176               Float_t ppjX = pPerp.Dot(ppJ2);
 177               Float_t ppjY = pPerp.Dot(ppJ3);
 178               Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
 179
 180               mxx += (ppjX * ppjX / ppjT);
 181               myy += (ppjY * ppjY / ppjT);
 182               mxy += (ppjX * ppjY / ppjT);
 183               nc++;
 184               sump2 += ppjT;
 185               // max pt
 186               if (pt > ptMax) {
 187                   ptMax  = pt;
 188                   iMax   = i;
 189                   etaMax = deta;
 190                   phiMax = dphi;
 191               }
 192           } // R < 0.4
 193       } // 1st Track Loop
 194       fPtSum->Fill(ptJ, ptsum);
 195
 196 //
 197 // At this point we have mxx, myy, mxy
 198       if (nc == 0) return;
 199 // Shericity Matrix
 200       const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
 201       TMatrixDSym m0(2, ele);
 202 // Find eigenvectors
 203       TMatrixDSymEigen m(m0);
 204       TVectorD eval(2);
 205       TMatrixD evecm = m.GetEigenVectors();
 206       eval  = m.GetEigenValues();
 207 // Largest eigenvector
 208       Int_t jev = 0;
 209       if (eval[0] < eval[1]) jev = 1;
 210       TVectorD evec0(2);
 211 // Principle axis
 212       evec0 = TMatrixDColumn(evecm, jev);
 213       TVector2 evec(evec0[0], evec0[1]);
 214 // Principle axis from leading partice
 215       Float_t phiM = TMath::ATan2(phiMax, etaMax);
 216       TVector2 evecM(TMath::Cos(phiM), TMath::Sin(phiM));
 217       Float_t phistM = evecM.DeltaPhi(evec);
 218       if (TMath::Abs(phistM) > TMath::Pi()/2.) evec*=(-1.);
 219
 220 //
 221 // 2nd correlation with principal axis
 222 //
 223       for (Int_t i = 0; i < nT; i++) {
 224           AliAODTrack* track = aodE->GetTrack(i);
 225 //
 226 //    Track quality cuts
 227           if (!track->TestFilterBit(1<<4)) continue;
 228 //
 229 //
 230           TVector3 pp(track->Px(), track->Py(), track->Pz());
 231           Float_t phi = track->Phi();
 232           Float_t eta = track->Eta();
 233           Float_t pt  = track->Pt();
 234           Float_t jT  = pp.Perp(ppJ1);
 235           // pT < 2.
 236           if (jT > 1.) continue;
 237           Float_t deta = eta - etaJ;
 238           Float_t dphi = phi - phiJ;
 239           if (dphi >   TMath::Pi()) dphi =   2. * TMath::Pi() - dphi;
 240           if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
 241
 242           Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
 243
 244           if (r > 0.7) continue;
 245
 246           TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
 247           TVector3 pPerp = pp - pLong;
 248
 249           Float_t ppjX = pPerp.Dot(ppJ2);
 250           Float_t ppjY = pPerp.Dot(ppJ3);
 251           Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
 252
 253           TVector2 vr(ppjX, ppjY) ;
 254           Float_t phistr = evec.DeltaPhi(vr);
 255           fPhiM->Fill(phistr);
 256           fPhiMPt ->Fill(phistr, pt);
 257           fPhiMPtJ->Fill(phistr, ptJ);
 258
 259       } // 2nd Track loop
 260
 261
 262   // Post output data.
 263   PostData(1, fHists);
 264 }
 265
 266 //________________________________________________________________________
 267 void AliAnalysisTaskJetsTM::Terminate(Option_t *)
 268 {
 269 // Terminate
 270 }