PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetv2QA.cxx

   1 // Jet v2 task using QA method, based on jet sample task (S.Aiola).
   2 //
   3 // Authors: Jason Mueller (CERN summer student 2014) & Alice Ohlson
   4
   5
   6 #include <TClonesArray.h>
   7 #include <TH1F.h>
   8 #include <TH2F.h>
   9 #include <TH3F.h>
  10 #include <TList.h>
  11 #include <TLorentzVector.h>
  12 #include <TF1.h>
  13
  14 #include "AliVCluster.h"
  15 #include "AliAODCaloCluster.h"
  16 #include "AliESDCaloCluster.h"
  17 #include "AliVTrack.h"
  18 #include "AliEmcalJet.h"
  19 #include "AliRhoParameter.h"
  20 #include "AliLog.h"
  21 #include "AliJetContainer.h"
  22 #include "AliParticleContainer.h"
  23 #include "AliClusterContainer.h"
  24 #include "AliPicoTrack.h"
  25
  26 #include "AliAnalysisTaskEmcalJetv2QA.h"
  27
  28 ClassImp(AliAnalysisTaskEmcalJetv2QA)
  29
  30 //________________________________________________________________________
  31 AliAnalysisTaskEmcalJetv2QA::AliAnalysisTaskEmcalJetv2QA() :
  32 AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetv2QA", kTRUE),
  33   nCentBins(0),
  34   nJetPtBins(0),
  35   centBins(0x0),
  36   jetPtBins(0x0),
  37   fJetv2(0x0),
  38   doPtWeight(kFALSE),
  39   fHistTracksPt(0),
  40   fHistClustersPt(0),
  41   fHistLeadingJetPt(0),
  42   fHistLeadingJetPtCorr(0),
  43   fHistJetsPhiEta(0),
  44   fHistJetsPtArea(0),
  45   fHistJetsPtLeadHad(0),
  46   fHistJetsCorrPtArea(0),
  47   fHistPtDEtaDPhiTrackClus(0),
  48   fHistPtDEtaDPhiClusTrack(0),
  49   fDPhiJet(0),
  50   fDPhiJetPythia(0),
  51   fDPhiEP(0),
  52   hGx(0),
  53   hGy2(0),
  54   hGxGy2(0),
  55   hGy4(0),
  56   hGx2(0),
  57   hGx2Gy2(0),
  58   hGxGy4(0),
  59   hGy6(0),
  60   hGx2Gy4(0),
  61   hGxGy6(0),
  62   hGy8(0),
  63   hGy(0),
  64   hN(0),
  65   htv2std(0),
  66   htjv2std(0),
  67   htj2v2std(0),
  68   hV0jv2std(0),
  69   htdPsi(0),
  70   htjdPsi(0),
  71   htj2dPsi(0),
  72   hV0jdPsi(0),
  73   hAx(0),
  74   hAxDijet(0),
  75   hQx(0),
  76   hQy(0),
  77   hEventData(0),
  78   hNTracks(0),
  79   hNTracksCent(0),
  80   hGxTracks(0),
  81   hGyTracks(0),
  82   hGy2Tracks(0),
  83   hGxGy2Tracks(0),
  84   hGy4Tracks(0),
  85   htEPRes(0),
  86   htjEPRes(0),
  87   htj2EPRes(0),
  88   fJetsCont(0),
  89   fTracksCont(0),
  90   fCaloClustersCont(0)
  91 {
  92   // Default constructor.
  93
  94   cout << endl;
  95   cout << "*******************************************************" << endl;
  96   cout << "*** AliAnalysisTaskEmcalJetv2QA constructor! ***" << endl;
  97   cout << "*******************************************************" << endl;
  98   cout << endl;
  99
 100
 101
 102   SetMakeGeneralHistograms(kTRUE);
 103 }
 104
 105 //________________________________________________________________________
 106 AliAnalysisTaskEmcalJetv2QA::AliAnalysisTaskEmcalJetv2QA(const char *name) :
 107   AliAnalysisTaskEmcalJet(name, kTRUE),
 108   nCentBins(0),
 109   nJetPtBins(0),
 110   centBins(0x0),
 111   jetPtBins(0x0),
 112   fJetv2(0x0),
 113   doPtWeight(kFALSE),
 114   fHistTracksPt(0),
 115   fHistClustersPt(0),
 116   fHistLeadingJetPt(0),
 117   fHistLeadingJetPtCorr(0),
 118   fHistJetsPhiEta(0),
 119   fHistJetsPtArea(0),
 120   fHistJetsPtLeadHad(0),
 121   fHistJetsCorrPtArea(0),
 122   fHistPtDEtaDPhiTrackClus(0),
 123   fHistPtDEtaDPhiClusTrack(0),
 124   fDPhiJet(0),
 125   fDPhiJetPythia(0),
 126   fDPhiEP(0),
 127   hGx(0),
 128   hGy2(0),
 129   hGxGy2(0),
 130   hGy4(0),
 131   hGx2(0),
 132   hGx2Gy2(0),
 133   hGxGy4(0),
 134   hGy6(0),
 135   hGx2Gy4(0),
 136   hGxGy6(0),
 137   hGy8(0),
 138   hGy(0),
 139   hN(0),
 140   htv2std(0),
 141   htjv2std(0),
 142   htj2v2std(0),
 143   hV0jv2std(0),
 144   htdPsi(0),
 145   htjdPsi(0),
 146   htj2dPsi(0),
 147   hV0jdPsi(0),
 148   hAx(0),
 149   hAxDijet(0),
 150   hQx(0),
 151   hQy(0),
 152   hEventData(0),
 153   hNTracks(0),
 154   hNTracksCent(0),
 155   hGxTracks(0),
 156   hGyTracks(0),
 157   hGy2Tracks(0),
 158   hGxGy2Tracks(0),
 159   hGy4Tracks(0),
 160   htEPRes(0),
 161   htjEPRes(0),
 162   htj2EPRes(0),
 163   fJetsCont(0),
 164   fTracksCont(0),
 165   fCaloClustersCont(0)
 166 {
 167   // Standard constructor.
 168
 169
 170   SetMakeGeneralHistograms(kTRUE);
 171 }
 172
 173 //________________________________________________________________________
 174 AliAnalysisTaskEmcalJetv2QA::~AliAnalysisTaskEmcalJetv2QA()
 175 {
 176   // Destructor.
 177 }
 178
 179 //________________________________________________________________________
 180 void AliAnalysisTaskEmcalJetv2QA::UserCreateOutputObjects()
 181 {
 182   // Create user output.
 183
 184   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 185
 186   fJetsCont           = GetJetContainer(0);
 187   if(fJetsCont) { //get particles and clusters connected to jets
 188     fTracksCont       = fJetsCont->GetParticleContainer();
 189     fCaloClustersCont = fJetsCont->GetClusterContainer();
 190   } else {        //no jets, just analysis tracks and clusters
 191     fTracksCont       = GetParticleContainer(0);
 192     fCaloClustersCont = GetClusterContainer(0);
 193   }
 194
 195   nCentBins = 6; // set number of centrality bins
 196   nJetPtBins = 4; // set number of jetPtBins
 197   centBins = new Double_t[nCentBins+1];
 198   jetPtBins = new Double_t[nJetPtBins+1];
 199   centBins[0] = 0.; centBins[1] = 5.; centBins[2] = 10.; centBins[3] = 20.; centBins[4] = 30.; centBins[5] = 40.; centBins[6] = 50.; // set edges of bins
 200   jetPtBins[0] = 50.; jetPtBins[1] = 70.; jetPtBins[2] = 90.; jetPtBins[3] = 120.; jetPtBins[4] = 200.; // set edges of bins
 201
 202   fTracksCont->SetClassName("AliVTrack");
 203   fCaloClustersCont->SetClassName("AliVCluster");
 204
 205   fHistTracksPt = new TH1F("fHistTracksPt","fHistTracksPt", fNbins / 2, fMinBinPt, fMaxBinPt / 2);
 206   fHistTracksPt->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
 207   fHistTracksPt->GetYaxis()->SetTitle("counts");
 208   fOutput->Add(fHistTracksPt);
 209
 210   fHistClustersPt = new TH1F("fHistClustersPt", "fHistClustersPt", fNbins / 2, fMinBinPt, fMaxBinPt / 2);
 211   fHistClustersPt->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
 212   fHistClustersPt->GetYaxis()->SetTitle("counts");
 213   fOutput->Add(fHistClustersPt);
 214
 215   fHistLeadingJetPt = new TH1F("fHistLeadingJetPt", "fHistLeadingJetPt", fNbins, fMinBinPt, fMaxBinPt);
 216   fHistLeadingJetPt->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 217   fHistLeadingJetPt->GetYaxis()->SetTitle("counts");
 218   fOutput->Add(fHistLeadingJetPt);
 219
 220   fHistLeadingJetPtCorr = new TH1F("fHistLeadingJetPtCorr", "fHistLeadingJetPtCorr", fNbins, fMinBinPt, fMaxBinPt);
 221   fHistLeadingJetPtCorr->GetXaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
 222   fHistLeadingJetPtCorr->GetYaxis()->SetTitle("counts");
 223   fOutput->Add(fHistLeadingJetPtCorr);
 224
 225   fHistJetsPhiEta = new TH2F("fHistJetsPhiEta", "fHistJetsPhiEta", 50, -1, 1, 101, 0, TMath::Pi()*2 + TMath::Pi()/200);
 226   fHistJetsPhiEta->GetXaxis()->SetTitle("#eta");
 227   fHistJetsPhiEta->GetYaxis()->SetTitle("#phi");
 228   fOutput->Add(fHistJetsPhiEta);
 229
 230   fHistJetsPtArea = new TH2F("fHistJetsPtArea", "fHistJetsPtArea", fNbins, fMinBinPt, fMaxBinPt, 30, 0, 3);
 231   fHistJetsPtArea->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 232   fHistJetsPtArea->GetYaxis()->SetTitle("area");
 233   fOutput->Add(fHistJetsPtArea);
 234
 235   fHistJetsPtLeadHad = new TH2F("fHistJetsPtLeadHad", "fHistJetsPtLeadHad", fNbins, fMinBinPt, fMaxBinPt, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
 236   fHistJetsPtLeadHad->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 237   fHistJetsPtLeadHad->GetYaxis()->SetTitle("p_{T,lead} (GeV/c)");
 238   fHistJetsPtLeadHad->GetZaxis()->SetTitle("counts");
 239   fOutput->Add(fHistJetsPtLeadHad);
 240
 241   fHistJetsCorrPtArea = new TH2F("fHistJetsCorrPtArea", "fHistJetsCorrPtArea", fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
 242   fHistJetsCorrPtArea->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
 243   fHistJetsCorrPtArea->GetYaxis()->SetTitle("area");
 244   fOutput->Add(fHistJetsCorrPtArea);
 245
 246   fHistPtDEtaDPhiTrackClus = new TH3F("fHistPtDEtaDPhiTrackClus","fHistPtDEtaDPhiTrackClus;#it{p}_{T}^{track};#Delta#eta;#Delta#varphi",100,0.,100.,100,-0.1,0.1,100,-0.1,0.1);
 247   fOutput->Add(fHistPtDEtaDPhiTrackClus);
 248
 249   fHistPtDEtaDPhiClusTrack = new TH3F("fHistPtDEtaDPhiClusTrack","fHistPtDEtaDPhiClusTrack;#it{p}_{T}^{clus};#Delta#eta;#Delta#varphi",100,0.,100.,100,-0.1,0.1,100,-0.1,0.1);
 250   fOutput->Add(fHistPtDEtaDPhiClusTrack);
 251
 252   fDPhiJet = new TH1F("fDPhiJet","fDPhiJet",90, -TMath::Pi()/3, 5*TMath::Pi()/3);
 253   fDPhiJet->GetXaxis()->SetTitle("#Delta#varphi");
 254   fDPhiJet->GetYaxis()->SetTitle("counts");
 255   fOutput->Add(fDPhiJet);
 256
 257   fDPhiJetPythia = new TH1F("fDPhiJetPythia","fDPhiJetPythia",90, -TMath::Pi()/3, 5*TMath::Pi()/3);
 258   fDPhiJetPythia->GetXaxis()->SetTitle("#Delta#varphi");
 259   fDPhiJetPythia->GetYaxis()->SetTitle("counts");
 260   fOutput->Add(fDPhiJetPythia);
 261
 262   fDPhiEP = new TH1F("fDPhiEP","fDPhiEP",90, 0, 2*TMath::Pi());
 263   fDPhiEP->GetXaxis()->SetTitle("#Delta#varphi");
 264   fDPhiEP->GetYaxis()->SetTitle("counts");
 265   fOutput->Add(fDPhiEP);
 266
 267   hGx = new TH2D("hGx", "Gx v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
 268   hGx->GetXaxis()->SetTitle("Centrality (%)");
 269   hGx->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 270   fOutput->Add(hGx);
 271
 272   hGy2 = new TH2D("hGy2", "Gy2 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
 273   hGy2->GetXaxis()->SetTitle("Centrality (%)");
 274   hGy2->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 275   fOutput->Add(hGy2);
 276
 277   hGxGy2 = new TH2D("hGxGy2", "GxGy2 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
 278   hGxGy2->GetXaxis()->SetTitle("Centrality (%)");
 279   hGxGy2->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 280   fOutput->Add(hGxGy2);
 281
 282   hGy4 = new TH2D("hGy4", "Gy4 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
 283   hGy4->GetXaxis()->SetTitle("Centrality (%)");
 284   hGy4->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 285   fOutput->Add(hGy4);
 286
 287   hGx2 = new TH2D("hGx2", "Gx2 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
 288   hGx2->GetXaxis()->SetTitle("Centrality (%)");
 289   hGx2->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 290   fOutput->Add(hGx2);
 291
 292   hGx2Gy2 = new TH2D("hGx2Gy2", "Gx2Gy2 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
 293   hGx2Gy2->GetXaxis()->SetTitle("Centrality (%)");
 294   hGx2Gy2->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 295   fOutput->Add(hGx2Gy2);
 296
 297   hGxGy4 = new TH2D("hGxGy4", "GxGy4 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
 298   hGxGy4->GetXaxis()->SetTitle("Centrality (%)");
 299   hGxGy4->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 300   fOutput->Add(hGxGy4);
 301
 302   hGy6 = new TH2D("hGy6", "Gy6 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
 303   hGy6->GetXaxis()->SetTitle("Centrality (%)");
 304   hGy6->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 305   fOutput->Add(hGy6);
 306
 307   hGx2Gy4 = new TH2D("hGx2Gy4", "Gx2Gy4 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
 308   hGx2Gy4->GetXaxis()->SetTitle("Centrality (%)");
 309   hGx2Gy4->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 310   fOutput->Add(hGx2Gy4);
 311
 312   hGxGy6 = new TH2D("hGxGy6", "GxGy6 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
 313   hGxGy6->GetXaxis()->SetTitle("Centrality (%)");
 314   hGxGy6->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 315   fOutput->Add(hGxGy6);
 316
 317   hGy8 = new TH2D("hGy8", "Gy8 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
 318   hGy8->GetXaxis()->SetTitle("Centrality (%)");
 319   hGy8->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 320   fOutput->Add(hGy8);
 321
 322   hGy = new TH2D("hGy", "Gy v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
 323   hGy->GetXaxis()->SetTitle("Centrality (%)");
 324   hGy->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 325   fOutput->Add(hGy);
 326
 327   hN = new TH2D("hN", "N v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
 328   hN->GetXaxis()->SetTitle("Centrality (%)");
 329   hN->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 330   fOutput->Add(hN);
 331
 332   htv2std = new TH2D("htv2std", "v2std v Centrality v JetPt w/o jets", nCentBins, centBins, nJetPtBins, jetPtBins);
 333   htv2std->GetXaxis()->SetTitle("Centrality (%)");
 334   htv2std->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 335   fOutput->Add(htv2std);
 336
 337   htjv2std = new TH2D("htjv2std", "v2std v Centrality v JetPt w/ jets", nCentBins, centBins, nJetPtBins, jetPtBins);
 338   htjv2std->GetXaxis()->SetTitle("Centrality (%)");
 339   htjv2std->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 340   fOutput->Add(htjv2std);
 341
 342   htj2v2std = new TH2D("htj2v2std", "v2std v Centrality v JetPt w/ trackPt < 2 GeV", nCentBins, centBins, nJetPtBins, jetPtBins);
 343   htj2v2std->GetXaxis()->SetTitle("Centrality (%)");
 344   htj2v2std->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 345   fOutput->Add(htj2v2std);
 346
 347   hV0jv2std = new TH2D("hV0jv2std", "v2std v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins);
 348   hV0jv2std->GetXaxis()->SetTitle("Centrality (%)");
 349   hV0jv2std->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 350   fOutput->Add(hV0jv2std);
 351
 352   Int_t ndpsibins = 100;
 353   Double_t dpsibins[101];
 354   for(Int_t t = 0; t < 101; t++) dpsibins[t] = TMath::Pi()*t/50.;
 355
 356   htdPsi = new TH3D("htdPsi", "JetAxis - EventPlane w/o jets", nCentBins, centBins, nJetPtBins, jetPtBins, ndpsibins, dpsibins);
 357   htdPsi->GetZaxis()->SetTitle("#Psi_{jet} - #Psi_{EP}");
 358   fOutput->Add(htdPsi);
 359
 360   htjdPsi = new TH3D("htjdPsi", "JetAxis - EventPlane w/ jets", nCentBins, centBins, nJetPtBins, jetPtBins, ndpsibins, dpsibins);
 361   htjdPsi->GetZaxis()->SetTitle("#Psi_{jet} - #Psi_{EP}");
 362   fOutput->Add(htjdPsi);
 363
 364   htj2dPsi = new TH3D("htj2dPsi", "JetAxis - EventPlane w/ trackPt < 2 GeV", nCentBins, centBins, nJetPtBins, jetPtBins, ndpsibins, dpsibins);
 365   htj2dPsi->GetZaxis()->SetTitle("#Psi_{jet} - #Psi_{EP}");
 366   fOutput->Add(htj2dPsi);
 367
 368   hV0jdPsi = new TH3D("hV0jdPsi", "JetAxis - EventPlane", nCentBins, centBins, nJetPtBins, jetPtBins, ndpsibins, dpsibins);
 369   hV0jdPsi->GetZaxis()->SetTitle("#Psi_{jet} - #Psi_{EP}");
 370   fOutput->Add(hV0jdPsi);
 371
 372   hQx = new TH2D("hQx", "Qx Distribution in EP frame", 100, -0.3, 0.3, nCentBins, centBins);
 373   hQx->GetXaxis()->SetTitle("Qx");
 374   hQx->GetYaxis()->SetTitle("Centrality (%)");
 375   fOutput->Add(hQx);
 376
 377   hQy = new TH2D("hQy", "Qy Distribution in EP frame", 100, -0.3, 0.3, nCentBins, centBins);
 378   hQy->GetXaxis()->SetTitle("Qy");
 379   hQy->GetYaxis()->SetTitle("Centrality (%)");
 380   fOutput->Add(hQy);
 381
 382   hAx = new TH2D("hAx", "Ax Distribution in EP frame w/o Dijets", 100, -35, 70, nJetPtBins, jetPtBins);
 383   hAx->GetXaxis()->SetTitle("Ax");
 384   hAx->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 385   fOutput->Add(hAx);
 386
 387   hAxDijet = new TH2D("hAxDijet", "Ax Distribution in EP frame w/ Dijets", 100, -35, 70, nJetPtBins, jetPtBins);
 388   hAxDijet->GetXaxis()->SetTitle("Ax");
 389   hAxDijet->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 390   fOutput->Add(hAxDijet);
 391
 392   hEventData = new TH1F("hEventData","Events Kept and Discarded", 9, 0, 9);
 393   hEventData->GetYaxis()->SetTitle("counts");
 394   fOutput->Add(hEventData);
 395
 396   hNTracks = new TH2F("hNTracks","Number of Tracks Per Event", 100, 0, 3000, 3, 0, 3);
 397   hNTracks->GetXaxis()->SetTitle("# tracks");
 398   fOutput->Add(hNTracks);
 399
 400   hNTracksCent = new TH2D("hNTracksCent", "NTracks by centrality", 100, 0, 3000, (Int_t)centBins[nCentBins], centBins[0], centBins[nCentBins]);
 401   hNTracksCent->GetXaxis()->SetTitle("# tracks");
 402   hNTracksCent->GetYaxis()->SetTitle("Centrality (%)");
 403   fOutput->Add(hNTracksCent);
 404
 405   hGxTracks = new TH2D("hGxTracks", "Gx by NTracks", 200, -200, 200, 100, 0, 3000);
 406   hGxTracks->GetXaxis()->SetTitle("Gx");
 407   hGxTracks->GetYaxis()->SetTitle("# tracks");
 408   fOutput->Add(hGxTracks);
 409
 410   hGyTracks = new TH2D("hGyTracks", "Gy by NTracks", 200, -200, 200, 100, 0, 3000);
 411   hGyTracks->GetXaxis()->SetTitle("Gy");
 412   hGyTracks->GetYaxis()->SetTitle("# tracks");
 413   fOutput->Add(hGyTracks);
 414
 415   hGy2Tracks = new TH2D("hGy2Tracks", "Gy2 by NTracks", 100, 0, 20000, 100, 0, 3000);
 416   hGy2Tracks->GetXaxis()->SetTitle("Gy2");
 417   hGy2Tracks->GetYaxis()->SetTitle("# tracks");
 418   fOutput->Add(hGy2Tracks);
 419
 420   hGxGy2Tracks = new TH2D("hGxGy2Tracks", "GxGy2 by NTracks", 100, -100000, 100000, 100, 0, 3000);
 421   hGxGy2Tracks->GetXaxis()->SetTitle("GxGy2");
 422   hGxGy2Tracks->GetYaxis()->SetTitle("# tracks");
 423   fOutput->Add(hGxGy2Tracks);
 424
 425   hGy4Tracks = new TH2D("hGy4Tracks", "Gy4 by NTracks", 100, 0, 100000000, 100, 0, 3000);
 426   hGy4Tracks->GetXaxis()->SetTitle("Gy4");
 427   hGy4Tracks->GetYaxis()->SetTitle("# tracks");
 428   fOutput->Add(hGy4Tracks);
 429
 430   htEPRes = new TH2D("htEPRes", "EP Resolution w/o Jets", nCentBins, centBins, nJetPtBins, jetPtBins);
 431   htEPRes->GetXaxis()->SetTitle("Centrality (%)");
 432   htEPRes->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 433   fOutput->Add(htEPRes);
 434
 435   htjEPRes = new TH2D("htjEPRes", "EP Resolution w/ Jets", nCentBins, centBins, nJetPtBins, jetPtBins);
 436   htjEPRes->GetXaxis()->SetTitle("Centrality (%)");
 437   htjEPRes->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 438   fOutput->Add(htjEPRes);
 439
 440   htj2EPRes = new TH2D("htj2EPRes", "EP Resolution w/ trackPT < 2 GeV", nCentBins, centBins, nJetPtBins, jetPtBins);
 441   htj2EPRes->GetXaxis()->SetTitle("Centrality (%)");
 442   htj2EPRes->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)");
 443   fOutput->Add(htj2EPRes);
 444
 445   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 446 }
 447
 448 //________________________________________________________________________
 449 Bool_t AliAnalysisTaskEmcalJetv2QA::FillHistograms()
 450 {
 451   // Fill histograms.
 452
 453   return kTRUE;
 454 }
 455
 456
 457 //________________________________________________________________________
 458 void AliAnalysisTaskEmcalJetv2QA::ExecOnce() {
 459
 460   AliAnalysisTaskEmcalJet::ExecOnce();
 461
 462   if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
 463   if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
 464   if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
 465
 466 }
 467
 468 //________________________________________________________________________
 469 Bool_t AliAnalysisTaskEmcalJetv2QA::Run() // this part loops over each event
 470 {
 471   // Run analysis code here, if needed. It will be executed before FillHistograms().
 472
 473   Double_t jetPhi = -999; // angle of leading jet
 474   Double_t jetPt = -999; // pt of leading jet
 475   Double_t jetArea = -999;
 476   Double_t trackPt = -999;
 477   Double_t phi = -999; // track phi
 478   Double_t dPhi = -999; // track phi - jet axis
 479   Double_t dPhiQA = -999; // track phi - EP
 480   Double_t tSin = 0; // used for std EP calc
 481   Double_t tCos = 0; // used for std EP calc
 482   Double_t jSin = 0; // used for std EP calc
 483   Double_t jCos = 0; // used for std EP calc
 484   Double_t tSin2 = 0; // used for std EP calc with trackPt < 2 GeV
 485   Double_t tCos2 = 0; // used for std EP calc with trackPt < 2 GeV
 486   Double_t qx = 0; // used for Qx distribution
 487   Double_t qy = 0; // used for Qy distribution
 488   Double_t ax = 0; // used for Ax distribution
 489   Double_t tEP = -999; // EP w/o jets
 490   Double_t tjEP = -999; // EP w/ jets
 491   Double_t tjEP2 = -999; // EP w/ jets
 492   Double_t dPsi = -999; // jet axis - EP
 493   Double_t gx = 0; // used for G moment calc
 494   Double_t gy = 0; // used for G moment calc
 495   Int_t isDijet = 0; // if 0, no dijet. if 1, dijet.
 496   Int_t nTracksBkgnd = 0; // used to keep track of number of tracks in background
 497   Int_t nTracksJet =0; // used to keep track of number of tracks in jets
 498
 499   if(fJetsCont)
 500     {
 501
 502       AliEmcalJet *jettest = fJetsCont->GetNextAcceptJet(0);
 503       while(jettest) {
 504
 505         fHistJetsPtArea->Fill(jettest->Pt(), jettest->Area());
 506         fHistJetsPhiEta->Fill(jettest->Eta(), jettest->Phi());
 507
 508         Float_t ptLeading = fJetsCont->GetLeadingHadronPt(jettest);
 509         fHistJetsPtLeadHad->Fill(jettest->Pt(), ptLeading);
 510
 511         if (fHistJetsCorrPtArea) {
 512           Float_t corrPt = jettest->Pt() - fJetsCont->GetRhoVal() * jettest->Area();
 513           fHistJetsCorrPtArea->Fill(corrPt, jettest->Area());
 514         }
 515         jettest = fJetsCont->GetNextAcceptJet();
 516       }
 517
 518
 519       AliEmcalJet *jet = fJetsCont->GetLeadingJet();
 520       if(jet)
 521         {
 522           jetPhi = jet->Phi(); // get leading jet phi (jet axis)
 523           jetPt = jet->Pt(); // get leading jet pT to filter out low pT jet events
 524           jetArea = jet->Area();
 525         }
 526     }
 527
 528   // event cuts
 529   if(!fTracksCont)
 530     {
 531       hEventData->Fill("!fTracksCont",1);
 532       return kFALSE;
 533     }
 534   if(!fJetsCont)
 535     {
 536       hEventData->Fill("!fJetsCont",1);
 537       return kFALSE;
 538     }
 539   if(jetPt == -999)
 540     {
 541       hEventData->Fill("jetPt=-999",1);
 542       return kFALSE;
 543     }
 544   if(jetPt < jetPtBins[0])
 545     {
 546       hEventData->Fill("leadingJetPt<jetPtMin",1);
 547       return kFALSE;
 548     }
 549   if(jetPt > jetPtBins[nJetPtBins])
 550     {
 551       hEventData->Fill("leadingJetPt>jetPtMax",1);
 552       return kFALSE;
 553     }
 554   if(fCent < centBins[0])
 555     {
 556       hEventData->Fill("cent<centMin",1);
 557       return kFALSE;
 558     }
 559   if(fCent > centBins[nCentBins])
 560     {
 561       hEventData->Fill("cent>centMax",1);
 562       return kFALSE;
 563     }
 564   hEventData->Fill("good event",1);
 565
 566   AliEmcalJet *dijet = fJetsCont->GetNextAcceptJet(0); // check for dijet events
 567   while(dijet)
 568     {
 569       if(dijet->Pt() > 50 && fabs(jetPhi-dijet->Phi()-TMath::Pi()) < 0.4) // loop over jets with pT>50 and exclude leading jet and check that angular separation is < 0.4
 570         isDijet = 1;
 571       dijet = fJetsCont->GetNextAcceptJet();
 572     }
 573
 574   if (fCaloClustersCont)
 575     {
 576       AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
 577       while(cluster) {
 578         TLorentzVector nPart;
 579         cluster->GetMomentum(nPart, fVertex);
 580         fHistClustersPt->Fill(nPart.Pt());
 581         cluster = fCaloClustersCont->GetNextAcceptCluster();
 582       }
 583     }
 584
 585   fHistLeadingJetPt->Fill(jetPt);
 586   fHistLeadingJetPtCorr->Fill(jetPt-fJetsCont->GetRhoVal()*jetArea);
 587
 588   AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
 589   while(track)
 590     { // loop over all particles (including jet tracks)
 591       trackPt = track->Pt();
 592       fHistTracksPt->Fill(trackPt);
 593       phi = track->Phi(); // get track phi
 594
 595       dPhi = phi-jetPhi; // get track phi - jet axis
 596       if(dPhi < 0) dPhi += TMath::TwoPi();
 597       if(dPhi > TMath::TwoPi()) dPhi -= TMath::TwoPi();
 598
 599       dPhiQA = phi-fEPV0; // get track phi - EP
 600       if(dPhiQA < 0) dPhiQA += TMath::TwoPi();
 601       if(dPhiQA > TMath::TwoPi()) dPhiQA -= TMath::TwoPi();
 602       fDPhiEP->Fill(dPhiQA);
 603
 604       // fill jet-hadron correlation just to check if track labels make sense...
 605       if(track->Pt()>1.)
 606         {
 607           Double_t dphiJet = dPhi;
 608           if(dphiJet > 5*TMath::Pi()/3) dphiJet -= 2*TMath::Pi();
 609           if(dphiJet < -TMath::Pi()/3) dphiJet += 2*TMath::Pi();
 610           if(track->GetLabel() == 0)
 611             fDPhiJet->Fill(dphiJet);
 612           else
 613             fDPhiJetPythia->Fill(dphiJet);
 614         }
 615
 616       Double_t weight = 1.;
 617       if(doPtWeight) weight = trackPt;
 618
 619       gx += weight*cos(2*dPhi);
 620       gy += weight*sin(2*dPhi);
 621
 622       if(track->GetLabel() == 0)
 623         { // sum for std EP method
 624           tSin += weight*sin(2*phi); // bkgnd has label = 0
 625           tCos += weight*cos(2*phi);
 626           qx += weight*cos(2*dPhiQA);
 627           qy += weight*sin(2*dPhiQA);
 628           nTracksBkgnd += 1;
 629         }
 630       else
 631         {
 632           jSin += weight*sin(2*phi); // jets have label =/= 0
 633           jCos += weight*cos(2*phi);
 634           ax += weight*cos(2*dPhi);
 635           nTracksJet += 1;
 636         }
 637
 638       if(track->Pt() < 2.)
 639         { // sum for std EP method w/ trackPt < 2 GeV
 640           tSin2 += weight*sin(2*phi);
 641           tCos2 += weight*cos(2*phi);
 642         }
 643
 644
 645       track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()); // increment to next track
 646     } // close loop over particles
 647
 648   hNTracks->Fill(nTracksBkgnd,"Bkgnd Tracks",1);
 649   hNTracks->Fill(nTracksJet,"Jet Tracks",1);
 650   hNTracks->Fill(nTracksBkgnd+nTracksJet,"Total Tracks",1);
 651   hNTracksCent->Fill(nTracksBkgnd+nTracksJet,fCent,1);
 652
 653   if(nTracksBkgnd == 0)
 654     {
 655       hEventData->Fill("no tracks",1);
 656       hEventData->Fill("good event",-1);
 657       return kFALSE;
 658     }
 659
 660   Double_t v2weight = 1+2*fJetv2*cos(2*(jetPhi-fEPV0)); // set v2 weight for event
 661
 662   hGx->Fill(fCent, jetPt, v2weight*gx); // fill histograms for QA method
 663   hGy2->Fill(fCent, jetPt, v2weight*gy*gy);
 664   hGxGy2->Fill(fCent, jetPt, v2weight*gx*gy*gy);
 665   hGy4->Fill(fCent, jetPt, v2weight*gy*gy*gy*gy);
 666   hGx2->Fill(fCent, jetPt, v2weight*gx*gx);
 667   hGx2Gy2->Fill(fCent, jetPt, v2weight*gx*gx*gy*gy);
 668   hGxGy4->Fill(fCent, jetPt, v2weight*gx*gy*gy*gy*gy);
 669   hGy6->Fill(fCent, jetPt, v2weight*gy*gy*gy*gy*gy*gy);
 670   hGx2Gy4->Fill(fCent, jetPt, v2weight*gx*gx*gy*gy*gy*gy);
 671   hGxGy6->Fill(fCent, jetPt, v2weight*gx*gy*gy*gy*gy*gy*gy);
 672   hGy8->Fill(fCent, jetPt, v2weight*gy*gy*gy*gy*gy*gy*gy*gy);
 673   hGy->Fill(fCent, jetPt, v2weight*gy);
 674   hN->Fill(fCent, jetPt, v2weight);
 675
 676   hGxTracks->Fill(gx,nTracksBkgnd+nTracksJet);
 677   hGyTracks->Fill(gy,nTracksBkgnd+nTracksJet);
 678   hGy2Tracks->Fill(gy*gy,nTracksBkgnd+nTracksJet);
 679   hGxGy2Tracks->Fill(gx*gy*gy,nTracksBkgnd+nTracksJet);
 680   hGy4Tracks->Fill(gy*gy*gy*gy,nTracksBkgnd+nTracksJet);
 681   hQx->Fill(qx/nTracksBkgnd,fCent);
 682   hQy->Fill(qy/nTracksBkgnd,fCent);
 683
 684   if(isDijet == 0)
 685     hAx->Fill(ax,jetPt);
 686   if(isDijet == 1)
 687     hAxDijet->Fill(ax,jetPt);
 688
 689   tEP = 0.5*atan2(tSin,tCos); // calculate EP w/o jets
 690   tjEP = 0.5*atan2((tSin+jSin),(tCos+jCos)); // calculate EP w/ jets
 691   tjEP2 = 0.5*atan2(tSin2,tCos2); // calculate EP w/ trackPt < 2 GeV
 692
 693   htEPRes->Fill(fCent, jetPt, v2weight*cos(2*(tEP-fEPV0)));
 694   htjEPRes->Fill(fCent, jetPt, v2weight*cos(2*(tjEP-fEPV0)));
 695   htj2EPRes->Fill(fCent, jetPt, v2weight*cos(2*(tjEP2-fEPV0)));
 696
 697
 698   dPsi = jetPhi-tEP;
 699   if(dPsi < 0) dPsi += TMath::TwoPi();
 700   if(dPsi > TMath::TwoPi()) dPsi -= TMath::TwoPi();
 701
 702   htv2std->Fill(fCent, jetPt, v2weight*cos(2*dPsi)); // fill histogram with v2 data w/o jets
 703   htdPsi->Fill(fCent,jetPt,dPsi); // fill histogram with jet axis - EP w/o jets
 704
 705   dPsi = jetPhi-tjEP;
 706   if(dPsi < 0) dPsi += TMath::TwoPi();
 707   if(dPsi > TMath::TwoPi()) dPsi -= TMath::TwoPi();
 708
 709   htjv2std->Fill(fCent, jetPt, v2weight*cos(2*dPsi)); // fill histogram with v2 data w/ jets
 710   htjdPsi->Fill(fCent,jetPt,dPsi); // fill histogram with jet axis - EP w/ jets
 711
 712   dPsi = jetPhi-tjEP2;
 713   if(dPsi < 0) dPsi += TMath::TwoPi();
 714   if(dPsi > TMath::TwoPi()) dPsi -= TMath::TwoPi();
 715
 716   htj2v2std->Fill(fCent, jetPt, v2weight*cos(2*dPsi)); // fill histogram with v2 data w/ trackPt < 2 GeV
 717   htj2dPsi->Fill(fCent,jetPt,dPsi); // fill histogram with jet axis - EP w/ trackPt < 2 GeV
 718
 719   dPsi = jetPhi-fEPV0;
 720   if(dPsi < 0) dPsi += TMath::TwoPi();
 721   if(dPsi > TMath::TwoPi()) dPsi -= TMath::TwoPi();
 722
 723   hV0jv2std->Fill(fCent, jetPt, v2weight*cos(2*dPsi)); // fill histogram with v2 data
 724   hV0jdPsi->Fill(fCent,jetPt,dPsi); // fill histogram with jet axis - EPV0
 725
 726
 727   return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
 728 }
 729
 730 //________________________________________________________________________
 731 void AliAnalysisTaskEmcalJetv2QA::Terminate(Option_t *)
 732 {
 733   // Called once at the end of the analysis.
 734   if(centBins) delete [] centBins;
 735   if(jetPtBins) delete [] jetPtBins;
 736 }