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