PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetShapeGR.cxx

   1 //
   2 // Analysis task for angular jet shape G(R) arXiv:1201.2688
   3 //
   4 // Author: M.Verweij
   5
   6 #include <TClonesArray.h>
   7 #include <TH1F.h>
   8 #include <TH2F.h>
   9 #include <TH3F.h>
  10 #include <THnSparse.h>
  11 #include <TF1.h>
  12 #include <TList.h>
  13 #include <TLorentzVector.h>
  14 #include <TProfile.h>
  15 #include <TChain.h>
  16 #include <TSystem.h>
  17 #include <TFile.h>
  18 #include <TKey.h>
  19 #include <TTree.h>
  20
  21 #include "AliVCluster.h"
  22 #include "AliVTrack.h"
  23 #include "AliEmcalJet.h"
  24 #include "AliRhoParameter.h"
  25 #include "AliLog.h"
  26 #include "AliEmcalParticle.h"
  27 #include "AliMCEvent.h"
  28 #include "AliAODEvent.h"
  29 #include "AliGenPythiaEventHeader.h"
  30 #include "AliAODMCHeader.h"
  31 #include "AliAnalysisManager.h"
  32 #include "AliJetContainer.h"
  33 #include "AliParticleContainer.h"
  34
  35 #include "AliAnalysisTaskJetShapeGR.h"
  36
  37 ClassImp(AliAnalysisTaskJetShapeGR)
  38
  39 //________________________________________________________________________
  40 AliAnalysisTaskJetShapeGR::AliAnalysisTaskJetShapeGR() :
  41   AliAnalysisTaskEmcalJet("AliAnalysisTaskJetShapeGR", kTRUE),
  42   fContainerBase(0),
  43   fContainerSub(1),
  44   fContainerTrue(2),
  45   fMinFractionShared(0),
  46   fSingleTrackEmb(kFALSE),
  47   fCreateTree(kFALSE),
  48   fTreeJetBkg(),
  49   fJet1Vec(new TLorentzVector()),
  50   fJet2Vec(new TLorentzVector()),
  51   fJetSubVec(new TLorentzVector()),
  52   fArea(0),
  53   fAreaPhi(0),
  54   fAreaEta(0),
  55   fRho(0),
  56   fRhoM(0),
  57   fNConst(0),
  58   fMatch(0),
  59   fDRStep(0.04),
  60   fMaxR(2.),
  61   fh2PtTrueDeltaGR(0x0),
  62   fh2PtTrueDeltaGRRel(0x0),
  63   fhnGRResponse(0x0),
  64   fh1PtTrue(0x0),
  65   fh3DeltaGRNumRPtTrue(0x0),
  66   fh3DeltaGRDenRPtTrue(0x0),
  67   fh2DeltaGRNumRPtTrue(0x0),
  68   fh2DeltaGRDenRPtTrue(0x0),
  69   fh1PtRaw(0x0),
  70   fh3DeltaGRNumRPtRaw(0x0),
  71   fh3DeltaGRDenRPtRaw(0x0),
  72   fh2DeltaGRNumRPtRaw(0x0),
  73   fh2DeltaGRDenRPtRaw(0x0),
  74   fh1PtRawMatch(0x0),
  75   fh3DeltaGRNumRPtRawMatch(0x0),
  76   fh3DeltaGRDenRPtRawMatch(0x0),
  77   fh2DeltaGRNumRPtRawMatch(0x0),
  78   fh2DeltaGRDenRPtRawMatch(0x0),
  79   fh1PtMatch(0x0),
  80   fh3DeltaGRNumRPtMatch(0x0),
  81   fh3DeltaGRDenRPtMatch(0x0),
  82   fh2DeltaGRNumRPtMatch(0x0),
  83   fh2DeltaGRDenRPtMatch(0x0),
  84   fh2DeltaGRNumRPtTrueMatch(0x0),
  85   fh2DeltaGRDenRPtTrueMatch(0x0)
  86 {
  87   // Default constructor.
  88
  89   fh2PtTrueDeltaGR     = new TH2F*[fNcentBins];
  90   fh2PtTrueDeltaGRRel  = new TH2F*[fNcentBins];
  91   fhnGRResponse        = new THnSparse*[fNcentBins];
  92   fh1PtTrue            = new TH1F*[fNcentBins];
  93   fh3DeltaGRNumRPtTrue = new TH3F*[fNcentBins];
  94   fh3DeltaGRDenRPtTrue = new TH3F*[fNcentBins];
  95   fh2DeltaGRNumRPtTrue = new TH2F*[fNcentBins];
  96   fh2DeltaGRDenRPtTrue = new TH2F*[fNcentBins];
  97   fh1PtRaw            = new TH1F*[fNcentBins];
  98   fh3DeltaGRNumRPtRaw = new TH3F*[fNcentBins];
  99   fh3DeltaGRDenRPtRaw = new TH3F*[fNcentBins];
 100   fh2DeltaGRNumRPtRaw = new TH2F*[fNcentBins];
 101   fh2DeltaGRDenRPtRaw = new TH2F*[fNcentBins];
 102   fh1PtRawMatch            = new TH1F*[fNcentBins];
 103   fh3DeltaGRNumRPtRawMatch = new TH3F*[fNcentBins];
 104   fh3DeltaGRDenRPtRawMatch = new TH3F*[fNcentBins];
 105   fh2DeltaGRNumRPtRawMatch = new TH2F*[fNcentBins];
 106   fh2DeltaGRDenRPtRawMatch = new TH2F*[fNcentBins];
 107   fh1PtMatch            = new TH1F*[fNcentBins];
 108   fh3DeltaGRNumRPtMatch = new TH3F*[fNcentBins];
 109   fh3DeltaGRDenRPtMatch = new TH3F*[fNcentBins];
 110   fh2DeltaGRNumRPtMatch = new TH2F*[fNcentBins];
 111   fh2DeltaGRDenRPtMatch = new TH2F*[fNcentBins];
 112   fh2DeltaGRNumRPtTrueMatch = new TH2F*[fNcentBins];
 113   fh2DeltaGRDenRPtTrueMatch = new TH2F*[fNcentBins];
 114
 115   for (Int_t i = 0; i < fNcentBins; i++) {
 116     fh2PtTrueDeltaGR[i]     = 0;
 117     fh2PtTrueDeltaGRRel[i]  = 0;
 118     fhnGRResponse[i]        = 0;
 119     fh1PtTrue[i]            = 0;
 120     fh3DeltaGRNumRPtTrue[i] = 0;
 121     fh3DeltaGRDenRPtTrue[i] = 0;
 122     fh2DeltaGRNumRPtTrue[i] = 0;
 123     fh2DeltaGRDenRPtTrue[i] = 0;
 124     fh1PtRaw[i]            = 0;
 125     fh3DeltaGRNumRPtRaw[i] = 0;
 126     fh3DeltaGRDenRPtRaw[i] = 0;
 127     fh2DeltaGRNumRPtRaw[i] = 0;
 128     fh2DeltaGRDenRPtRaw[i] = 0;
 129     fh1PtRawMatch[i]            = 0;
 130     fh3DeltaGRNumRPtRawMatch[i] = 0;
 131     fh3DeltaGRDenRPtRawMatch[i] = 0;
 132     fh2DeltaGRNumRPtRawMatch[i] = 0;
 133     fh2DeltaGRDenRPtRawMatch[i] = 0;
 134     fh1PtMatch[i]            = 0;
 135     fh3DeltaGRNumRPtMatch[i] = 0;
 136     fh3DeltaGRDenRPtMatch[i] = 0;
 137     fh2DeltaGRNumRPtMatch[i] = 0;
 138     fh2DeltaGRDenRPtMatch[i] = 0;
 139     fh2DeltaGRNumRPtTrueMatch[i] = 0;
 140     fh2DeltaGRDenRPtTrueMatch[i] = 0;
 141   }
 142
 143   SetMakeGeneralHistograms(kTRUE);
 144   DefineOutput(2, TTree::Class());
 145 }
 146
 147 //________________________________________________________________________
 148 AliAnalysisTaskJetShapeGR::AliAnalysisTaskJetShapeGR(const char *name) :
 149   AliAnalysisTaskEmcalJet(name, kTRUE),
 150   fContainerBase(0),
 151   fContainerSub(1),
 152   fContainerTrue(2),
 153   fMinFractionShared(0),
 154   fSingleTrackEmb(kFALSE),
 155   fCreateTree(kFALSE),
 156   fTreeJetBkg(0),
 157   fJet1Vec(new TLorentzVector()),
 158   fJet2Vec(new TLorentzVector()),
 159   fJetSubVec(new TLorentzVector()),
 160   fArea(0),
 161   fAreaPhi(0),
 162   fAreaEta(0),
 163   fRho(0),
 164   fRhoM(0),
 165   fNConst(0),
 166   fMatch(0),
 167   fDRStep(0.04),
 168   fMaxR(2.),
 169   fh2PtTrueDeltaGR(0x0),
 170   fh2PtTrueDeltaGRRel(0x0),
 171   fhnGRResponse(0x0),
 172   fh1PtTrue(0x0),
 173   fh3DeltaGRNumRPtTrue(0x0),
 174   fh3DeltaGRDenRPtTrue(0x0),
 175   fh2DeltaGRNumRPtTrue(0x0),
 176   fh2DeltaGRDenRPtTrue(0x0),
 177   fh1PtRaw(0x0),
 178   fh3DeltaGRNumRPtRaw(0x0),
 179   fh3DeltaGRDenRPtRaw(0x0),
 180   fh2DeltaGRNumRPtRaw(0x0),
 181   fh2DeltaGRDenRPtRaw(0x0),
 182   fh1PtRawMatch(0x0),
 183   fh3DeltaGRNumRPtRawMatch(0x0),
 184   fh3DeltaGRDenRPtRawMatch(0x0),
 185   fh2DeltaGRNumRPtRawMatch(0x0),
 186   fh2DeltaGRDenRPtRawMatch(0x0),
 187   fh1PtMatch(0x0),
 188   fh3DeltaGRNumRPtMatch(0x0),
 189   fh3DeltaGRDenRPtMatch(0x0),
 190   fh2DeltaGRNumRPtMatch(0x0),
 191   fh2DeltaGRDenRPtMatch(0x0),
 192   fh2DeltaGRNumRPtTrueMatch(0x0),
 193   fh2DeltaGRDenRPtTrueMatch(0x0)
 194 {
 195   // Standard constructor.
 196
 197   fh2PtTrueDeltaGR     = new TH2F*[fNcentBins];
 198   fh2PtTrueDeltaGRRel  = new TH2F*[fNcentBins];
 199   fhnGRResponse        = new THnSparse*[fNcentBins];
 200   fh1PtTrue            = new TH1F*[fNcentBins];
 201   fh3DeltaGRNumRPtTrue = new TH3F*[fNcentBins];
 202   fh3DeltaGRDenRPtTrue = new TH3F*[fNcentBins];
 203   fh2DeltaGRNumRPtTrue = new TH2F*[fNcentBins];
 204   fh2DeltaGRDenRPtTrue = new TH2F*[fNcentBins];
 205   fh1PtRaw            = new TH1F*[fNcentBins];
 206   fh3DeltaGRNumRPtRaw = new TH3F*[fNcentBins];
 207   fh3DeltaGRDenRPtRaw = new TH3F*[fNcentBins];
 208   fh2DeltaGRNumRPtRaw = new TH2F*[fNcentBins];
 209   fh2DeltaGRDenRPtRaw = new TH2F*[fNcentBins];
 210   fh1PtRawMatch            = new TH1F*[fNcentBins];
 211   fh3DeltaGRNumRPtRawMatch = new TH3F*[fNcentBins];
 212   fh3DeltaGRDenRPtRawMatch = new TH3F*[fNcentBins];
 213   fh2DeltaGRNumRPtRawMatch = new TH2F*[fNcentBins];
 214   fh2DeltaGRDenRPtRawMatch = new TH2F*[fNcentBins];
 215   fh1PtMatch            = new TH1F*[fNcentBins];
 216   fh3DeltaGRNumRPtMatch = new TH3F*[fNcentBins];
 217   fh3DeltaGRDenRPtMatch = new TH3F*[fNcentBins];
 218   fh2DeltaGRNumRPtMatch = new TH2F*[fNcentBins];
 219   fh2DeltaGRDenRPtMatch = new TH2F*[fNcentBins];
 220   fh2DeltaGRNumRPtTrueMatch = new TH2F*[fNcentBins];
 221   fh2DeltaGRDenRPtTrueMatch = new TH2F*[fNcentBins];
 222
 223   for (Int_t i = 0; i < fNcentBins; i++) {
 224     fh2PtTrueDeltaGR[i]     = 0;
 225     fh2PtTrueDeltaGRRel[i]  = 0;
 226     fhnGRResponse[i]        = 0;
 227     fh1PtTrue[i]            = 0;
 228     fh3DeltaGRNumRPtTrue[i] = 0;
 229     fh3DeltaGRDenRPtTrue[i] = 0;
 230     fh2DeltaGRNumRPtTrue[i] = 0;
 231     fh2DeltaGRDenRPtTrue[i] = 0;
 232     fh1PtRaw[i]            = 0;
 233     fh3DeltaGRNumRPtRaw[i] = 0;
 234     fh3DeltaGRDenRPtRaw[i] = 0;
 235     fh2DeltaGRNumRPtRaw[i] = 0;
 236     fh2DeltaGRDenRPtRaw[i] = 0;
 237     fh1PtRawMatch[i]            = 0;
 238     fh3DeltaGRNumRPtRawMatch[i] = 0;
 239     fh3DeltaGRDenRPtRawMatch[i] = 0;
 240     fh2DeltaGRNumRPtRawMatch[i] = 0;
 241     fh2DeltaGRDenRPtRawMatch[i] = 0;
 242     fh1PtMatch[i]            = 0;
 243     fh3DeltaGRNumRPtMatch[i] = 0;
 244     fh3DeltaGRDenRPtMatch[i] = 0;
 245     fh2DeltaGRNumRPtMatch[i] = 0;
 246     fh2DeltaGRDenRPtMatch[i] = 0;
 247     fh2DeltaGRNumRPtTrueMatch[i] = 0;
 248     fh2DeltaGRDenRPtTrueMatch[i] = 0;
 249   }
 250
 251   SetMakeGeneralHistograms(kTRUE);
 252   DefineOutput(2, TTree::Class());
 253 }
 254
 255 //________________________________________________________________________
 256 AliAnalysisTaskJetShapeGR::~AliAnalysisTaskJetShapeGR()
 257 {
 258   // Destructor.
 259 }
 260
 261 //________________________________________________________________________
 262 void AliAnalysisTaskJetShapeGR::UserCreateOutputObjects()
 263 {
 264   // Create user output.
 265
 266   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 267
 268   Bool_t oldStatus = TH1::AddDirectoryStatus();
 269   TH1::AddDirectory(kFALSE);
 270
 271   const Int_t nBinsPt  = 200;
 272   const Double_t minPt = -50.;
 273   const Double_t maxPt = 150.;
 274
 275   const Int_t nBinsM  = 150;
 276   const Double_t minM = -50.;
 277   const Double_t maxM = 100.;
 278
 279   const Int_t nBinsR  = 50;
 280   const Double_t minR = 0.;
 281   const Double_t maxR = 2.;
 282
 283   const Int_t nBinsDGR  = 100;
 284   const Double_t minDGR = 0.;
 285   const Double_t maxDGR = 10.;
 286
 287   //Binning for THnSparse
 288   const Int_t nBinsSparse0 = 4;
 289   const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt};
 290   const Double_t xmin0[nBinsSparse0]  = { minM, minM, minPt, minPt};
 291   const Double_t xmax0[nBinsSparse0]  = { maxM, maxM, maxPt, maxPt};
 292
 293   TString histName = "";
 294   TString histTitle = "";
 295   for (Int_t i = 0; i < fNcentBins; i++) {
 296     histName = Form("fh2PtTrueDeltaGR_%d",i);
 297     histTitle = Form("fh2PtTrueDeltaGR_%d;#it{p}_{T,true};#it{G(R)}_{sub}-#it{G(R)}_{true}",i);
 298     fh2PtTrueDeltaGR[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,100,-50.,50.);
 299     fOutput->Add(fh2PtTrueDeltaGR[i]);
 300
 301     histName = Form("fh2PtTrueDeltaGRRel_%d",i);
 302     histTitle = Form("fh2PtTrueDeltaGRRel_%d;#it{p}_{T,true};(#it{G(R)}_{sub}-#it{G(R)}_{true})/#it{G(R)}_{true}",i);
 303     fh2PtTrueDeltaGRRel[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,200,-1.,1.);
 304     fOutput->Add(fh2PtTrueDeltaGRRel[i]);
 305
 306     histName = Form("fhnGRResponse_%d",i);
 307     histTitle = Form("fhnGRResponse_%d;#it{G(R)}_{sub};#it{G(R)}_{true};#it{p}_{T,sub};#it{p}_{T,true}",i);
 308     fhnGRResponse[i] = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
 309     fOutput->Add(fhnGRResponse[i]);
 310
 311     //Histos for true jets
 312     histName = Form("fh1PtTrue_%d",i);
 313     histTitle = Form("%s;#it{p}_{T};#it{N}",histName.Data());
 314     fh1PtTrue[i] = new TH1F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
 315     fOutput->Add(fh1PtTrue[i]);
 316
 317     histName = Form("fh3DeltaGRNumRPtTrue_%d",i);
 318     histTitle = Form("%s;#Sigma#it{p}_{Tk,i}#it{p}_{Tk,j}#Delta#it{R}_{ij}^{2}#delta_{dR}(#it{r}-#it{R}_{ij});#it{r};#it{p}_{T,true}",histName.Data());
 319     fh3DeltaGRNumRPtTrue[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsDGR,minDGR,maxDGR,nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 320     fOutput->Add(fh3DeltaGRNumRPtTrue[i]);
 321
 322     histName = Form("fh3DeltaGRDenRPtTrue_%d",i);
 323     histTitle = Form("%s;#Sigma#it{p}_{Tk,i}#it{p}_{Tk,j}#Delta#it{R}_{ij}^{2}#Theta_{dR}(#it{r}-#it{R}_{ij});#it{r};#it{p}_{T,true}",histName.Data());
 324     fh3DeltaGRDenRPtTrue[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsDGR,minDGR,maxDGR,nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 325     fOutput->Add(fh3DeltaGRDenRPtTrue[i]);
 326
 327     histName = Form("fh2DeltaGRNumRPtTrue_%d",i);
 328     histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
 329     fh2DeltaGRNumRPtTrue[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 330     fOutput->Add(fh2DeltaGRNumRPtTrue[i]);
 331
 332     histName = Form("fh2DeltaGRDenRPtTrue_%d",i);
 333     histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
 334     fh2DeltaGRDenRPtTrue[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 335     fOutput->Add(fh2DeltaGRDenRPtTrue[i]);
 336
 337     //Histos for raw AA jets
 338     histName = Form("fh1PtRaw_%d",i);
 339     histTitle = Form("%s;#it{p}_{T};#it{N}",histName.Data());
 340     fh1PtRaw[i] = new TH1F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
 341     fOutput->Add(fh1PtRaw[i]);
 342
 343     histName = Form("fh3DeltaGRNumRPtRaw_%d",i);
 344     histTitle = Form("%s;#Sigma#it{p}_{Tk,i}#it{p}_{Tk,j}#Delta#it{R}_{ij}^{2}#delta_{dR}(#it{r}-#it{R}_{ij});#it{r};#it{p}_{T,true}",histName.Data());
 345     fh3DeltaGRNumRPtRaw[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsDGR,minDGR,maxDGR,nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 346     fOutput->Add(fh3DeltaGRNumRPtRaw[i]);
 347
 348     histName = Form("fh3DeltaGRDenRPtRaw_%d",i);
 349     histTitle = Form("%s;#Sigma#it{p}_{Tk,i}#it{p}_{Tk,j}#Delta#it{R}_{ij}^{2}#Theta_{dR}(#it{r}-#it{R}_{ij});#it{r};#it{p}_{T,true}",histName.Data());
 350     fh3DeltaGRDenRPtRaw[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsDGR,minDGR,maxDGR,nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 351     fOutput->Add(fh3DeltaGRDenRPtRaw[i]);
 352
 353     histName = Form("fh2DeltaGRNumRPtRaw_%d",i);
 354     histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
 355     fh2DeltaGRNumRPtRaw[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 356     fOutput->Add(fh2DeltaGRNumRPtRaw[i]);
 357
 358     histName = Form("fh2DeltaGRDenRPtRaw_%d",i);
 359     histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
 360     fh2DeltaGRDenRPtRaw[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 361     fOutput->Add(fh2DeltaGRDenRPtRaw[i]);
 362
 363     //Histos for raw AA jets matched to MC jet
 364     histName = Form("fh1PtRawMatch_%d",i);
 365     histTitle = Form("%s;#it{p}_{T};#it{N}",histName.Data());
 366     fh1PtRawMatch[i] = new TH1F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
 367     fOutput->Add(fh1PtRawMatch[i]);
 368
 369     histName = Form("fh3DeltaGRNumRPtRawMatch_%d",i);
 370     histTitle = Form("%s;#Sigma#it{p}_{Tk,i}#it{p}_{Tk,j}#Delta#it{R}_{ij}^{2}#delta_{dR}(#it{r}-#it{R}_{ij});#it{r};#it{p}_{T,true}",histName.Data());
 371     fh3DeltaGRNumRPtRawMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsDGR,minDGR,maxDGR,nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 372     fOutput->Add(fh3DeltaGRNumRPtRawMatch[i]);
 373
 374     histName = Form("fh3DeltaGRDenRPtRawMatch_%d",i);
 375     histTitle = Form("%s;#Sigma#it{p}_{Tk,i}#it{p}_{Tk,j}#Delta#it{R}_{ij}^{2}#Theta_{dR}(#it{r}-#it{R}_{ij});#it{r};#it{p}_{T,true}",histName.Data());
 376     fh3DeltaGRDenRPtRawMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsDGR,minDGR,maxDGR,nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 377     fOutput->Add(fh3DeltaGRDenRPtRawMatch[i]);
 378
 379     histName = Form("fh2DeltaGRNumRPtRawMatch_%d",i);
 380     histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
 381     fh2DeltaGRNumRPtRawMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 382     fOutput->Add(fh2DeltaGRNumRPtRawMatch[i]);
 383
 384     histName = Form("fh2DeltaGRDenRPtRawMatch_%d",i);
 385     histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
 386     fh2DeltaGRDenRPtRawMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 387     fOutput->Add(fh2DeltaGRDenRPtRawMatch[i]);
 388
 389     //Histos for matched jets
 390     histName = Form("fh1PtMatch_%d",i);
 391     histTitle = Form("%s;#it{p}_{T};#it{N}",histName.Data());
 392     fh1PtMatch[i] = new TH1F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
 393     fOutput->Add(fh1PtMatch[i]);
 394
 395     histName = Form("fh3DeltaGRNumRPtMatch_%d",i);
 396     histTitle = Form("%s;#Sigma#it{p}_{Tk,i}#it{p}_{Tk,j}#Delta#it{R}_{ij}^{2}#delta_{dR}(#it{r}-#it{R}_{ij});#it{r};#it{p}_{T,true}",histName.Data());
 397     fh3DeltaGRNumRPtMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsDGR,minDGR,maxDGR,nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 398     fOutput->Add(fh3DeltaGRNumRPtMatch[i]);
 399
 400     histName = Form("fh3DeltaGRDenRPtMatch_%d",i);
 401     histTitle = Form("%s;#Sigma#it{p}_{Tk,i}#it{p}_{Tk,j}#Delta#it{R}_{ij}^{2}#Theta_{dR}(#it{r}-#it{R}_{ij});#it{r};#it{p}_{T,true}",histName.Data());
 402     fh3DeltaGRDenRPtMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsDGR,minDGR,maxDGR,nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 403     fOutput->Add(fh3DeltaGRDenRPtMatch[i]);
 404
 405     histName = Form("fh2DeltaGRNumRPtMatch_%d",i);
 406     histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
 407     fh2DeltaGRNumRPtMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 408     fOutput->Add(fh2DeltaGRNumRPtMatch[i]);
 409
 410     histName = Form("fh2DeltaGRDenRPtMatch_%d",i);
 411     histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
 412     fh2DeltaGRDenRPtMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 413     fOutput->Add(fh2DeltaGRDenRPtMatch[i]);
 414
 415     histName = Form("fh2DeltaGRNumRPtTrueMatch_%d",i);
 416     histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
 417     fh2DeltaGRNumRPtTrueMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 418     fOutput->Add(fh2DeltaGRNumRPtTrueMatch[i]);
 419
 420     histName = Form("fh2DeltaGRDenRPtTrueMatch_%d",i);
 421     histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
 422     fh2DeltaGRDenRPtTrueMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
 423     fOutput->Add(fh2DeltaGRDenRPtTrueMatch[i]);
 424
 425   }
 426
 427   // =========== Switch on Sumw2 for all histos ===========
 428   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 429     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 430     if (h1){
 431       h1->Sumw2();
 432       continue;
 433     }
 434     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 435     if(hn)hn->Sumw2();
 436   }
 437
 438   TH1::AddDirectory(oldStatus);
 439
 440   // Create a tree.
 441   if(fCreateTree) {
 442     fTreeJetBkg = new TTree("fTreeJetSubGR", "fTreeJetSubGR");
 443     fTreeJetBkg->Branch("fJet1Vec","TLorentzVector",&fJet1Vec);
 444     fTreeJetBkg->Branch("fJet2Vec","TLorentzVector",&fJet2Vec);
 445     fTreeJetBkg->Branch("fJetSubVec","TLorentzVector",&fJetSubVec);
 446     fTreeJetBkg->Branch("fArea",&fArea,"fArea/F");
 447     fTreeJetBkg->Branch("fAreaPhi",&fAreaPhi,"fAreaPhi/F");
 448     fTreeJetBkg->Branch("fAreaEta",&fAreaEta,"fAreaEta/F");
 449     fTreeJetBkg->Branch("fRho",&fRho,"fRho/F");
 450     fTreeJetBkg->Branch("fRhoM",&fRhoM,"fRhoM/F");
 451     fTreeJetBkg->Branch("fMatch",&fMatch,"fMatch/I");
 452   }
 453
 454   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 455   if(fCreateTree) PostData(2, fTreeJetBkg);
 456 }
 457
 458 //________________________________________________________________________
 459 Bool_t AliAnalysisTaskJetShapeGR::Run()
 460 {
 461   // Run analysis code here, if needed. It will be executed before FillHistograms().
 462
 463   return kTRUE;
 464 }
 465
 466 //________________________________________________________________________
 467 Bool_t AliAnalysisTaskJetShapeGR::FillHistograms()
 468 {
 469   // Fill histograms.
 470   FillTrueJets();
 471
 472   AliEmcalJet* jet1 = NULL;
 473   AliEmcalJet *jet2 = NULL;
 474   AliEmcalJet *jetS = NULL;
 475   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 476   AliJetContainer *jetContS = GetJetContainer(fContainerSub);
 477   AliDebug(11,Form("NJets  Incl: %d  Csub: %d",jetCont->GetNJets(),jetContS->GetNJets()));
 478   if(jetCont && jetContS) {
 479     jetCont->ResetCurrentID();
 480     Double_t rmax = jetCont->GetJetRadius()+0.2;
 481     Double_t wr = 0.04;
 482     Int_t nr = TMath::CeilNint(rmax/wr);
 483     while((jet1 = jetCont->GetNextAcceptJet())) {
 484       Double_t fraction = 0.;
 485       fMatch = 0;
 486       fJet2Vec->SetPtEtaPhiM(0.,0.,0.,0.);
 487       if(fSingleTrackEmb) {
 488         AliVParticle *vp = GetEmbeddedConstituent(jet1);
 489         if(vp) {
 490           fJet2Vec->SetPxPyPzE(vp->Px(),vp->Py(),vp->Pz(),vp->E());
 491           fMatch = 1;
 492         }
 493       } else {
 494         jet2 = jet1->ClosestJet();
 495         if(!jet2) continue;
 496
 497         fraction = jetCont->GetFractionSharedPt(jet1);
 498         fMatch = 1;
 499         if(fMinFractionShared>0.) {
 500           if(fraction>fMinFractionShared) {
 501             fJet2Vec->SetPxPyPzE(jet2->Px(),jet2->Py(),jet2->Pz(),jet2->E());
 502             fMatch = 1;
 503           } else
 504             fMatch = 0;
 505         }
 506       }
 507
 508       //Fill histograms for all AA jets
 509       Double_t ptcorr = jet1->Pt()-jetCont->GetRhoVal()*jet1->Area();
 510       fh1PtRaw[fCentBin]->Fill(ptcorr);
 511       if(fMatch==1) fh1PtRawMatch[fCentBin]->Fill(ptcorr);
 512
 513       TArrayF numRaw = jet1->GetGRNumerator();
 514       TArrayF denRaw = jet1->GetGRDenominator();
 515       if(numRaw.GetSize()>0) {
 516         for(Int_t i = 0; i<nr; i++) {
 517           Double_t r = i*wr + 0.5*wr;
 518           fh3DeltaGRNumRPtRaw[fCentBin]->Fill(numRaw[i],r,ptcorr);
 519           fh3DeltaGRDenRPtRaw[fCentBin]->Fill(denRaw[i],r,ptcorr);
 520           fh2DeltaGRNumRPtRaw[fCentBin]->Fill(r,ptcorr,numRaw[i]);
 521           fh2DeltaGRDenRPtRaw[fCentBin]->Fill(r,ptcorr,denRaw[i]);
 522           if(fMatch==1) {
 523             fh3DeltaGRNumRPtRawMatch[fCentBin]->Fill(numRaw[i],r,ptcorr);
 524             fh3DeltaGRDenRPtRawMatch[fCentBin]->Fill(denRaw[i],r,ptcorr);
 525             fh2DeltaGRNumRPtRawMatch[fCentBin]->Fill(r,ptcorr,numRaw[i]);
 526             fh2DeltaGRDenRPtRawMatch[fCentBin]->Fill(r,ptcorr,denRaw[i]);
 527           }
 528         }
 529       }
 530
 531       //Fill histograms for matched jets
 532       if(fMatch==1) {
 533         fh1PtMatch[fCentBin]->Fill(ptcorr);
 534
 535         //now get second derivative vs R and do final calculation
 536         TArrayF num = jet1->GetGRNumeratorSub();
 537         TArrayF den = jet1->GetGRDenominatorSub();
 538         if(num.GetSize()>0) {
 539           for(Int_t i = 0; i<nr; i++) {
 540             Double_t r = i*wr + 0.5*wr;
 541             fh3DeltaGRNumRPtMatch[fCentBin]->Fill(num[i],r,ptcorr);
 542             fh3DeltaGRDenRPtMatch[fCentBin]->Fill(den[i],r,ptcorr);
 543             fh2DeltaGRNumRPtMatch[fCentBin]->Fill(r,ptcorr,num[i]);
 544             fh2DeltaGRDenRPtMatch[fCentBin]->Fill(r,ptcorr,den[i]);
 545             fh2DeltaGRNumRPtTrueMatch[fCentBin]->Fill(r,jet2->Pt(),num[i]);
 546             fh2DeltaGRDenRPtTrueMatch[fCentBin]->Fill(r,jet2->Pt(),den[i]);
 547
 548             Double_t dGR = 0.;
 549             fh2PtTrueDeltaGR[fCentBin]->Fill(jet2->Pt(),dGR);
 550           }
 551         }
 552       }
 553
 554       if(fCreateTree) {
 555         fJet1Vec->SetPxPyPzE(jet1->Px(),jet1->Py(),jet1->Pz(),jet1->E());
 556         if(jetS && jetS->Pt()>0.) fJetSubVec->SetPtEtaPhiM(jetS->Pt(),jetS->Eta(),jetS->Phi(),jetS->M());
 557         else fJetSubVec->SetPtEtaPhiM(0.,0.,0.,0.);
 558         fArea = (Float_t)jet1->Area();
 559         fAreaPhi = (Float_t)jet1->AreaPhi();
 560         fAreaEta = (Float_t)jet1->AreaEta();
 561         fRho  = (Float_t)jetCont->GetRhoVal();
 562         fRhoM = (Float_t)jetCont->GetRhoMassVal();
 563         fNConst = (Int_t)jet1->GetNumberOfTracks();
 564         fTreeJetBkg->Fill();
 565       }
 566     } //jet1 loop
 567   }//jetCont
 568
 569
 570   return kTRUE;
 571 }
 572
 573 //________________________________________________________________________
 574 Bool_t AliAnalysisTaskJetShapeGR::FillTrueJets() {
 575
 576   AliEmcalJet* jet1 = NULL;
 577   AliJetContainer *jetCont = GetJetContainer(fContainerTrue);
 578   if(!jetCont)
 579     return kFALSE;
 580
 581   AliDebug(11,Form("NJets True: %d",jetCont->GetNJets()));
 582
 583   //create arrays
 584   const Int_t nr = TMath::CeilNint(fMaxR/fDRStep);
 585   TArrayF *fNum = new TArrayF(nr);
 586   TArrayF *fDen = new TArrayF(nr);
 587
 588   if(jetCont) {
 589     jetCont->ResetCurrentID();
 590     while((jet1 = jetCont->GetNextAcceptJet())) {
 591       fh1PtTrue[fCentBin]->Fill(jet1->Pt());
 592
 593       //Double_t dev = CalcDeltaGR(jet1,fContainerTrue,fNum,fDen);//num,den);
 594       for(Int_t i = 0; i<nr; i++) {
 595         Double_t r = i*fDRStep + 0.5*fDRStep;
 596         fh3DeltaGRNumRPtTrue[fCentBin]->Fill(fNum->At(i),r,jet1->Pt());
 597         fh3DeltaGRDenRPtTrue[fCentBin]->Fill(fDen->At(i),r,jet1->Pt());
 598         fh2DeltaGRNumRPtTrue[fCentBin]->Fill(r,jet1->Pt(),fNum->At(i));
 599         fh2DeltaGRDenRPtTrue[fCentBin]->Fill(r,jet1->Pt(),fDen->At(i));
 600       }
 601     }
 602   }
 603   if(fNum)  delete fNum;
 604   if(fDen)  delete fDen;
 605   return kTRUE;
 606 }
 607
 608 //________________________________________________________________________
 609 Double_t AliAnalysisTaskJetShapeGR::CalcDeltaGR(AliEmcalJet *jet, Int_t ic, TArrayF *fNum, TArrayF *fDen) { //Double_t *num, Double_t *den) {
 610   //Calculate G(R)
 611
 612   //First clear the arrays
 613   const Int_t nr = TMath::CeilNint(fMaxR/fDRStep);
 614   for(Int_t i = 0; i<nr; i++) {
 615     fNum->SetAt(0.,i);
 616     fDen->SetAt(0.,i);
 617   }
 618
 619   AliJetContainer *jetCont = GetJetContainer(ic);
 620   AliVParticle *vp1 = 0x0;
 621   AliVParticle *vp2 = 0x0;
 622   Double_t A = 0.; Double_t B = 0.;
 623   if(jet->GetNumberOfTracks()<2) return 0.;
 624   for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
 625     vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
 626     if(!vp1) continue;
 627     for(Int_t j=i+1; j<jet->GetNumberOfTracks(); j++) {
 628       vp2 = static_cast<AliVParticle*>(jet->TrackAt(j, jetCont->GetParticleContainer()->GetArray()));
 629       if(!vp2) continue;
 630       Double_t dphi = GetDeltaPhi(vp1->Phi(),vp2->Phi());
 631       Double_t dr2 = (vp1->Eta()-vp2->Eta())*(vp1->Eta()-vp2->Eta()) + dphi*dphi;
 632       if(dr2>0.) {
 633         for(Int_t k = 0; k<nr; k++) {
 634           Double_t r = k*fDRStep + 0.5*fDRStep;
 635           //    Double_t x = jetCont->GetJetRadius()-TMath::Sqrt(dr2);
 636           Double_t dr = TMath::Sqrt(dr2);
 637           Double_t x = r-dr;
 638           //noisy function
 639           Double_t noise = TMath::Exp(-x*x/(2*fDRStep*fDRStep))/(TMath::Sqrt(2.*TMath::Pi())*fDRStep);
 640           //error function
 641           Double_t erf = 0.5*(1.+TMath::Erf(x/(TMath::Sqrt(2.)*fDRStep)));
 642
 643           A = vp1->Pt()*vp2->Pt()*dr2*noise;
 644           B = vp1->Pt()*vp2->Pt()*dr2*erf;
 645           fNum->AddAt(fNum->At(k)+A,k);
 646           fDen->AddAt(fDen->At(k)+B,k);
 647         }
 648       }
 649     }
 650   }
 651
 652   Double_t deltaGR = 0.;
 653   if(B>0.) deltaGR = A/B; //useless
 654   return deltaGR;
 655 }
 656
 657 //________________________________________________________________________
 658 Double_t AliAnalysisTaskJetShapeGR::CalcGR(AliEmcalJet *jet, Int_t ic) {
 659   //Calculate G(R)
 660   AliJetContainer *jetCont = GetJetContainer(ic);
 661   AliVParticle *vp1 = 0x0;
 662   AliVParticle *vp2 = 0x0;
 663   Double_t gR = 0.;
 664   Double_t wr = 0.04;
 665   const Int_t nr = TMath::CeilNint(jetCont->GetJetRadius()/wr);
 666   Double_t grArr[nr];
 667   for(Int_t i = 0; i<nr; i++)
 668     grArr[i] = 0.;
 669
 670   for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
 671     vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
 672     for(Int_t j=i; j<jet->GetNumberOfTracks(); j++) {
 673       vp2 = static_cast<AliVParticle*>(jet->TrackAt(j, jetCont->GetParticleContainer()->GetArray()));
 674       Double_t dphi = GetDeltaPhi(vp1->Phi(),vp2->Phi());
 675       Double_t dr2 = (vp1->Eta()-vp2->Eta())*(vp1->Eta()-vp2->Eta()) + dphi*dphi;
 676       Int_t bin = TMath::FloorNint(TMath::Sqrt(dr2)/wr);
 677       Double_t gr = vp1->Pt()*vp2->Pt()*dr2;
 678       if(bin<nr) grArr[bin]+=gr;
 679
 680       if(TMath::Sqrt(dr2)<jetCont->GetJetRadius())
 681         gR += gr;
 682     }
 683   }
 684   return gR;
 685 }
 686
 687 //________________________________________________________________________
 688 AliVParticle* AliAnalysisTaskJetShapeGR::GetEmbeddedConstituent(AliEmcalJet *jet) {
 689
 690   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 691   AliVParticle *vp = 0x0;
 692   AliVParticle *vpe = 0x0; //embedded particle
 693   Int_t nc = 0;
 694   for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
 695     vp = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); //check if fTracks is the correct track branch
 696     if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
 697     if(!vpe) vpe = vp;
 698     else if(vp->Pt()>vpe->Pt()) vpe = vp;
 699     nc++;
 700   }
 701   AliDebug(11,Form("Found %d embedded particles",nc));
 702   return vpe;
 703 }
 704
 705 //________________________________________________________________________
 706 Double_t AliAnalysisTaskJetShapeGR::GetDeltaPhi(Double_t phi1,Double_t phi2) {
 707   //
 708   // Calculate azimuthal angle difference
 709   //
 710
 711   Double_t dPhi = phi1-phi2;
 712
 713   if(dPhi <-1.*TMath::Pi())  dPhi += TMath::TwoPi();
 714   if(dPhi > 1.*TMath::Pi())  dPhi -= TMath::TwoPi();
 715
 716   return dPhi;
 717 }
 718
 719
 720 //________________________________________________________________________
 721 Bool_t AliAnalysisTaskJetShapeGR::RetrieveEventObjects() {
 722   //
 723   // retrieve event objects
 724   //
 725
 726   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 727     return kFALSE;
 728
 729   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 730   jetCont->LoadRhoMass(InputEvent());
 731
 732   return kTRUE;
 733 }
 734
 735 //_______________________________________________________________________
 736 void AliAnalysisTaskJetShapeGR::Terminate(Option_t *)
 737 {
 738   // Called once at the end of the analysis.
 739 }
 740