PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetMassResponse.cxx

   1 //
   2 // Jet mass response analysis task.
   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 "AliGenPythiaEventHeader.h"
  29 #include "AliAODMCHeader.h"
  30 #include "AliMCEvent.h"
  31 #include "AliAnalysisManager.h"
  32 #include "AliJetContainer.h"
  33 #include "AliParticleContainer.h"
  34
  35 #include "AliAODEvent.h"
  36
  37 #include "AliAnalysisTaskEmcalJetMassResponse.h"
  38
  39 ClassImp(AliAnalysisTaskEmcalJetMassResponse)
  40
  41 //________________________________________________________________________
  42 AliAnalysisTaskEmcalJetMassResponse::AliAnalysisTaskEmcalJetMassResponse() :
  43   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetMassResponse", kTRUE),
  44   fContainerBase(0),
  45   fMinFractionShared(0),
  46   f1JetMassAvg(0),
  47   fSingleTrackEmb(kFALSE),
  48   fSubtractMassless(kFALSE),
  49   fCreateTree(kFALSE),
  50   fh2PtJet1DeltaMNoSub(0),
  51   fh2PtJet2DeltaMNoSub(0),
  52   fh3PtJet1DeltaPtDeltaMCheat(0),
  53   fh3PtJet2DeltaPtDeltaMCheat(0),
  54   fh3PtJet1DeltaPtDeltaM(0),
  55   fh3PtJet2DeltaPtDeltaM(0),
  56   fh2PtJet1DeltaE(0),
  57   fh2PtJet2DeltaE(0),
  58   fh2PtJet1DeltaP(0),
  59   fh2PtJet2DeltaP(0),
  60   fh2PtJet2DeltaM(0),
  61   fh3PtJet1MJet1MJet2(0),
  62   fh3PtJet2MJet1MJet2(0),
  63   fh3PtJet1DeltaPtDeltaMRho(0),
  64   fh2PtJet1DeltaERho(0),
  65   fh3PtJet2DeltaPtDeltaMRho(0),
  66   fh2PtJet2DeltaPxRho(0),
  67   fh2PtJet2DeltaPyRho(0),
  68   fh2PtJet2DeltaPzRho(0),
  69   fh2PtJet2DeltaERho(0),
  70   fh2PtJet2DeltaMRho(0),
  71   fh2PtJet2DeltaPtRho(0),
  72   fh3PtJet2DeltaEDeltaMRho(0),
  73   fh3PtJet2DeltaPDeltaMRho(0),
  74   fh2PtJet1DeltaPtVecSub(0),
  75   fTreeJetBkg(),
  76   fJet1Vec(new TLorentzVector()),
  77   fJet2Vec(new TLorentzVector()),
  78   fBkgVec(new TLorentzVector()),
  79   fArea(0),
  80   fAreaPhi(0),
  81   fAreaEta(0),
  82   fRho(0),
  83   fRhoM(0),
  84   fNConst(0),
  85   fJetMassMassless(0)
  86 {
  87   // Default constructor.
  88
  89   fh2PtJet1DeltaMNoSub         = new TH2F*[fNcentBins];
  90   fh2PtJet2DeltaMNoSub         = new TH2F*[fNcentBins];
  91   fh3PtJet1DeltaPtDeltaMCheat  = new TH3F*[fNcentBins];
  92   fh3PtJet2DeltaPtDeltaMCheat  = new TH3F*[fNcentBins];
  93   fh3PtJet1DeltaPtDeltaM       = new TH3F*[fNcentBins];
  94   fh3PtJet2DeltaPtDeltaM       = new TH3F*[fNcentBins];
  95   fh2PtJet1DeltaE              = new TH2F*[fNcentBins];
  96   fh2PtJet2DeltaE              = new TH2F*[fNcentBins];
  97   fh2PtJet1DeltaP              = new TH2F*[fNcentBins];
  98   fh2PtJet2DeltaP              = new TH2F*[fNcentBins];
  99   fh2PtJet2DeltaM              = new TH2F*[fNcentBins];
 100   fh3PtJet1MJet1MJet2          = new TH3F*[fNcentBins];
 101   fh3PtJet2MJet1MJet2          = new TH3F*[fNcentBins];
 102   fh3PtJet1DeltaPtDeltaMRho    = new TH3F*[fNcentBins];
 103   fh2PtJet1DeltaERho           = new TH2F*[fNcentBins];
 104   fh3PtJet2DeltaPtDeltaMRho    = new TH3F*[fNcentBins];
 105   fh2PtJet2DeltaPxRho          = new TH2F*[fNcentBins];
 106   fh2PtJet2DeltaPyRho          = new TH2F*[fNcentBins];
 107   fh2PtJet2DeltaPzRho          = new TH2F*[fNcentBins];
 108   fh2PtJet2DeltaERho           = new TH2F*[fNcentBins];
 109   fh2PtJet2DeltaMRho           = new TH2F*[fNcentBins];
 110   fh2PtJet2DeltaPtRho          = new TH2F*[fNcentBins];
 111   fh3PtJet2DeltaEDeltaMRho     = new TH3F*[fNcentBins];
 112   fh3PtJet2DeltaPDeltaMRho     = new TH3F*[fNcentBins];
 113   fh2PtJet1DeltaPtVecSub       = new TH2F*[fNcentBins];
 114
 115   for (Int_t i = 0; i < fNcentBins; i++) {
 116     fh2PtJet1DeltaMNoSub[i]        = 0;
 117     fh2PtJet2DeltaMNoSub[i]        = 0;
 118     fh3PtJet1DeltaPtDeltaMCheat[i] = 0;
 119     fh3PtJet2DeltaPtDeltaMCheat[i] = 0;
 120     fh3PtJet1DeltaPtDeltaM[i]      = 0;
 121     fh3PtJet2DeltaPtDeltaM[i]      = 0;
 122     fh2PtJet1DeltaE[i]             = 0;
 123     fh2PtJet2DeltaE[i]             = 0;
 124     fh2PtJet1DeltaP[i]             = 0;
 125     fh2PtJet2DeltaP[i]             = 0;
 126     fh2PtJet2DeltaM[i]             = 0;
 127     fh3PtJet1MJet1MJet2[i]         = 0;
 128     fh3PtJet2MJet1MJet2[i]         = 0;
 129     fh3PtJet1DeltaPtDeltaMRho[i]   = 0;
 130     fh2PtJet1DeltaERho[i]          = 0;
 131     fh3PtJet2DeltaPtDeltaMRho[i]   = 0;
 132     fh2PtJet2DeltaPxRho[i]         = 0;
 133     fh2PtJet2DeltaPyRho[i]         = 0;
 134     fh2PtJet2DeltaPzRho[i]         = 0;
 135     fh2PtJet2DeltaERho[i]          = 0;
 136     fh2PtJet2DeltaMRho[i]          = 0;
 137     fh2PtJet2DeltaPtRho[i]         = 0;
 138     fh3PtJet2DeltaEDeltaMRho[i]    = 0;
 139     fh3PtJet2DeltaPDeltaMRho[i]    = 0;
 140     fh2PtJet1DeltaPtVecSub[i]      = 0;
 141   }
 142   SetMakeGeneralHistograms(kTRUE);
 143   DefineOutput(2, TTree::Class());
 144 }
 145
 146 //________________________________________________________________________
 147 AliAnalysisTaskEmcalJetMassResponse::AliAnalysisTaskEmcalJetMassResponse(const char *name) :
 148   AliAnalysisTaskEmcalJet(name, kTRUE),
 149   fContainerBase(0),
 150   fMinFractionShared(0),
 151   f1JetMassAvg(0),
 152   fSingleTrackEmb(kFALSE),
 153   fSubtractMassless(kFALSE),
 154   fCreateTree(kFALSE),
 155   fh2PtJet1DeltaMNoSub(0),
 156   fh2PtJet2DeltaMNoSub(0),
 157   fh3PtJet1DeltaPtDeltaMCheat(0),
 158   fh3PtJet2DeltaPtDeltaMCheat(0),
 159   fh3PtJet1DeltaPtDeltaM(0),
 160   fh3PtJet2DeltaPtDeltaM(0),
 161   fh2PtJet1DeltaE(0),
 162   fh2PtJet2DeltaE(0),
 163   fh2PtJet1DeltaP(0),
 164   fh2PtJet2DeltaP(0),
 165   fh2PtJet2DeltaM(0),
 166   fh3PtJet1MJet1MJet2(0),
 167   fh3PtJet2MJet1MJet2(0),
 168   fh3PtJet1DeltaPtDeltaMRho(0),
 169   fh2PtJet1DeltaERho(0),
 170   fh3PtJet2DeltaPtDeltaMRho(0),
 171   fh2PtJet2DeltaPxRho(0),
 172   fh2PtJet2DeltaPyRho(0),
 173   fh2PtJet2DeltaPzRho(0),
 174   fh2PtJet2DeltaERho(0),
 175   fh2PtJet2DeltaMRho(0),
 176   fh2PtJet2DeltaPtRho(0),
 177   fh3PtJet2DeltaEDeltaMRho(0),
 178   fh3PtJet2DeltaPDeltaMRho(0),
 179   fh2PtJet1DeltaPtVecSub(0),
 180   fTreeJetBkg(0),
 181   fJet1Vec(new TLorentzVector()),
 182   fJet2Vec(new TLorentzVector()),
 183   fBkgVec(new TLorentzVector()),
 184   fArea(0),
 185   fAreaPhi(0),
 186   fAreaEta(0),
 187   fRho(0),
 188   fRhoM(0),
 189   fNConst(0),
 190   fJetMassMassless(0)
 191 {
 192   // Standard constructor.
 193
 194   fh2PtJet1DeltaMNoSub         = new TH2F*[fNcentBins];
 195   fh2PtJet2DeltaMNoSub         = new TH2F*[fNcentBins];
 196   fh3PtJet1DeltaPtDeltaMCheat  = new TH3F*[fNcentBins];
 197   fh3PtJet2DeltaPtDeltaMCheat  = new TH3F*[fNcentBins];
 198   fh3PtJet1DeltaPtDeltaM       = new TH3F*[fNcentBins];
 199   fh3PtJet2DeltaPtDeltaM       = new TH3F*[fNcentBins];
 200   fh2PtJet1DeltaE              = new TH2F*[fNcentBins];
 201   fh2PtJet2DeltaE              = new TH2F*[fNcentBins];
 202   fh2PtJet1DeltaP              = new TH2F*[fNcentBins];
 203   fh2PtJet2DeltaP              = new TH2F*[fNcentBins];
 204   fh2PtJet2DeltaM              = new TH2F*[fNcentBins];
 205   fh3PtJet1MJet1MJet2          = new TH3F*[fNcentBins];
 206   fh3PtJet2MJet1MJet2          = new TH3F*[fNcentBins];
 207   fh3PtJet1DeltaPtDeltaMRho    = new TH3F*[fNcentBins];
 208   fh2PtJet1DeltaERho           = new TH2F*[fNcentBins];
 209   fh3PtJet2DeltaPtDeltaMRho    = new TH3F*[fNcentBins];
 210   fh2PtJet2DeltaPxRho          = new TH2F*[fNcentBins];
 211   fh2PtJet2DeltaPyRho          = new TH2F*[fNcentBins];
 212   fh2PtJet2DeltaPzRho          = new TH2F*[fNcentBins];
 213   fh2PtJet2DeltaERho           = new TH2F*[fNcentBins];
 214   fh2PtJet2DeltaMRho           = new TH2F*[fNcentBins];
 215   fh2PtJet2DeltaPtRho          = new TH2F*[fNcentBins];
 216   fh3PtJet2DeltaEDeltaMRho     = new TH3F*[fNcentBins];
 217   fh3PtJet2DeltaPDeltaMRho     = new TH3F*[fNcentBins];
 218   fh2PtJet1DeltaPtVecSub       = new TH2F*[fNcentBins];
 219
 220   for (Int_t i = 0; i < fNcentBins; i++) {
 221     fh2PtJet1DeltaMNoSub[i]        = 0;
 222     fh2PtJet2DeltaMNoSub[i]        = 0;
 223     fh3PtJet1DeltaPtDeltaMCheat[i] = 0;
 224     fh3PtJet2DeltaPtDeltaMCheat[i] = 0;
 225     fh3PtJet1DeltaPtDeltaM[i]      = 0;
 226     fh3PtJet2DeltaPtDeltaM[i]      = 0;
 227     fh2PtJet1DeltaE[i]             = 0;
 228     fh2PtJet2DeltaE[i]             = 0;
 229     fh2PtJet1DeltaP[i]             = 0;
 230     fh2PtJet2DeltaP[i]             = 0;
 231     fh2PtJet2DeltaM[i]             = 0;
 232     fh3PtJet1MJet1MJet2[i]         = 0;
 233     fh3PtJet2MJet1MJet2[i]         = 0;
 234     fh3PtJet1DeltaPtDeltaMRho[i]   = 0;
 235     fh2PtJet1DeltaERho[i]          = 0;
 236     fh3PtJet2DeltaPtDeltaMRho[i]   = 0;
 237     fh2PtJet2DeltaPxRho[i]         = 0;
 238     fh2PtJet2DeltaPyRho[i]         = 0;
 239     fh2PtJet2DeltaPzRho[i]         = 0;
 240     fh2PtJet2DeltaERho[i]          = 0;
 241     fh2PtJet2DeltaMRho[i]          = 0;
 242     fh2PtJet2DeltaPtRho[i]         = 0;
 243     fh3PtJet2DeltaEDeltaMRho[i]    = 0;
 244     fh3PtJet2DeltaPDeltaMRho[i]    = 0;
 245     fh2PtJet1DeltaPtVecSub[i]      = 0;
 246   }
 247
 248   SetMakeGeneralHistograms(kTRUE);
 249   DefineOutput(2, TTree::Class());
 250 }
 251
 252 //________________________________________________________________________
 253 AliAnalysisTaskEmcalJetMassResponse::~AliAnalysisTaskEmcalJetMassResponse()
 254 {
 255   // Destructor.
 256 }
 257
 258 //________________________________________________________________________
 259 void AliAnalysisTaskEmcalJetMassResponse::UserCreateOutputObjects()
 260 {
 261   // Create user output.
 262
 263   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 264
 265   Bool_t oldStatus = TH1::AddDirectoryStatus();
 266   TH1::AddDirectory(kFALSE);
 267
 268   const Int_t nBinsPt  = 200;
 269   const Double_t minPt = -50.;
 270   const Double_t maxPt = 150.;
 271
 272   const Int_t nBinsM  = 150;
 273   const Double_t minM = -50.;
 274   const Double_t maxM = 100.;
 275
 276   TString histName = "";
 277   TString histTitle = "";
 278   for (Int_t i = 0; i < fNcentBins; i++) {
 279     histName = TString::Format("fh2PtJet1DeltaMNoSub_%d",i);
 280     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{M}_{jet}",histName.Data());
 281     fh2PtJet1DeltaMNoSub[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 282     fOutput->Add(fh2PtJet1DeltaMNoSub[i]);
 283
 284     histName = TString::Format("fh2PtJet2DeltaMNoSub_%d",i);
 285     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{M}_{jet}",histName.Data());
 286     fh2PtJet2DeltaMNoSub[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 287     fOutput->Add(fh2PtJet2DeltaMNoSub[i]);
 288
 289     histName = TString::Format("fh3PtJet1DeltaPtDeltaMCheat_%d",i);
 290     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{p}_{T};#delta#it{M}_{jet}",histName.Data());
 291     fh3PtJet1DeltaPtDeltaMCheat[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 292     fOutput->Add(fh3PtJet1DeltaPtDeltaMCheat[i]);
 293
 294     histName = TString::Format("fh3PtJet2DeltaPtDeltaMCheat_%d",i);
 295     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{p}_{T};#delta#it{M}_{jet}",histName.Data());
 296     fh3PtJet2DeltaPtDeltaMCheat[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 297     fOutput->Add(fh3PtJet2DeltaPtDeltaMCheat[i]);
 298
 299     histName = TString::Format("fh3PtJet1DeltaPtDeltaM_%d",i);
 300     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{p}_{T};#delta#it{M}_{jet}",histName.Data());
 301     fh3PtJet1DeltaPtDeltaM[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 302     fOutput->Add(fh3PtJet1DeltaPtDeltaM[i]);
 303
 304     histName = TString::Format("fh3PtJet2DeltaPtDeltaM_%d",i);
 305     histTitle = TString::Format("%s;#it{p}_{T,jet2};#delta#it{p}_{T};#delta#it{M}_{jet}",histName.Data());
 306     fh3PtJet2DeltaPtDeltaM[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 307     fOutput->Add(fh3PtJet2DeltaPtDeltaM[i]);
 308
 309     histName = TString::Format("fh2PtJet1DeltaE_%d",i);
 310     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{E}",histName.Data());
 311     fh2PtJet1DeltaE[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 312     fOutput->Add(fh2PtJet1DeltaE[i]);
 313
 314     histName = TString::Format("fh2PtJet2DeltaE_%d",i);
 315     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{E}",histName.Data());
 316     fh2PtJet2DeltaE[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 317     fOutput->Add(fh2PtJet2DeltaE[i]);
 318
 319     histName = TString::Format("fh2PtJet1DeltaP_%d",i);
 320     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{p}",histName.Data());
 321     fh2PtJet1DeltaP[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 322     fOutput->Add(fh2PtJet1DeltaP[i]);
 323
 324     histName = TString::Format("fh2PtJet2DeltaP_%d",i);
 325     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{p}",histName.Data());
 326     fh2PtJet2DeltaP[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 327     fOutput->Add(fh2PtJet2DeltaP[i]);
 328
 329     histName = TString::Format("fh2PtJet2DeltaM_%d",i);
 330     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{M}",histName.Data());
 331     fh2PtJet2DeltaM[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 332     fOutput->Add(fh2PtJet2DeltaM[i]);
 333
 334     histName = TString::Format("fh3PtJet1MJet1MJet2_%d",i);
 335     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{M}_{jet2}",histName.Data());
 336     fh3PtJet1MJet1MJet2[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsM,minM,maxM);
 337     fOutput->Add(fh3PtJet1MJet1MJet2[i]);
 338
 339     histName = TString::Format("fh3PtJet2MJet1MJet2_%d",i);
 340     histTitle = TString::Format("%s;#it{p}_{T,jet2};#it{M}_{jet1};#it{M}_{jet2}",histName.Data());
 341     fh3PtJet2MJet1MJet2[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsM,minM,maxM);
 342     fOutput->Add(fh3PtJet2MJet1MJet2[i]);
 343
 344     histName = TString::Format("fh3PtJet1DeltaPtDeltaMRho_%d",i);
 345     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{p}_{T};#delta#it{M}_{jet}",histName.Data());
 346     fh3PtJet1DeltaPtDeltaMRho[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 347     fOutput->Add(fh3PtJet1DeltaPtDeltaMRho[i]);
 348
 349     histName = TString::Format("fh2PtJet1DeltaERho_%d",i);
 350     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{E}",histName.Data());
 351     fh2PtJet1DeltaERho[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 352     fOutput->Add(fh2PtJet1DeltaERho[i]);
 353
 354     histName = TString::Format("fh3PtJet2DeltaPtDeltaMRho_%d",i);
 355     histTitle = TString::Format("%s;#it{p}_{T,jet2};#delta#it{p}_{T};#delta#it{M}_{jet}",histName.Data());
 356     fh3PtJet2DeltaPtDeltaMRho[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 357     fOutput->Add(fh3PtJet2DeltaPtDeltaMRho[i]);
 358
 359     histName = TString::Format("fh2PtJet2DeltaPxRho_%d",i);
 360     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{p}_{x}",histName.Data());
 361     fh2PtJet2DeltaPxRho[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 362     fOutput->Add(fh2PtJet2DeltaPxRho[i]);
 363
 364     histName = TString::Format("fh2PtJet2DeltaPyRho_%d",i);
 365     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{p}_{y}",histName.Data());
 366     fh2PtJet2DeltaPyRho[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 367     fOutput->Add(fh2PtJet2DeltaPyRho[i]);
 368
 369     histName = TString::Format("fh2PtJet2DeltaPzRho_%d",i);
 370     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{p}_{z}",histName.Data());
 371     fh2PtJet2DeltaPzRho[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 372     fOutput->Add(fh2PtJet2DeltaPzRho[i]);
 373
 374     histName = TString::Format("fh2PtJet2DeltaERho_%d",i);
 375     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{E}",histName.Data());
 376     fh2PtJet2DeltaERho[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 377     fOutput->Add(fh2PtJet2DeltaERho[i]);
 378
 379     histName = TString::Format("fh2PtJet2DeltaMRho_%d",i);
 380     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{M}",histName.Data());
 381     fh2PtJet2DeltaMRho[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 382     fOutput->Add(fh2PtJet2DeltaMRho[i]);
 383
 384     histName = TString::Format("fh2PtJet2DeltaPtRho_%d",i);
 385     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{p}_{T}",histName.Data());
 386     fh2PtJet2DeltaPtRho[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 387     fOutput->Add(fh2PtJet2DeltaPtRho[i]);
 388
 389     histName = TString::Format("fh3PtJet2DeltaEDeltaMRho_%d",i);
 390     histTitle = TString::Format("%s;#it{p}_{T,jet2};#delta#it{E};#delta#it{M}_{jet}",histName.Data());
 391     fh3PtJet2DeltaEDeltaMRho[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 392     fOutput->Add(fh3PtJet2DeltaEDeltaMRho[i]);
 393
 394     histName = TString::Format("fh3PtJet2DeltaPDeltaMRho_%d",i);
 395     histTitle = TString::Format("%s;#it{p}_{T,jet2};#delta#it{p};#delta#it{M}_{jet}",histName.Data());
 396     fh3PtJet2DeltaPDeltaMRho[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 397     fOutput->Add(fh3PtJet2DeltaPDeltaMRho[i]);
 398
 399     histName = TString::Format("fh2PtJet1DeltaPtVecSub_%d",i);
 400     histTitle = TString::Format("%s;#it{p}_{T,jet1};#delta#it{p}_{T}",histName.Data());
 401     fh2PtJet1DeltaPtVecSub[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 402     fOutput->Add(fh2PtJet1DeltaPtVecSub[i]);
 403   }
 404
 405   // =========== Switch on Sumw2 for all histos ===========
 406   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 407     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 408     if (h1){
 409       h1->Sumw2();
 410       continue;
 411     }
 412     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 413     if(hn)hn->Sumw2();
 414   }
 415
 416   TH1::AddDirectory(oldStatus);
 417
 418   // Create a tree.
 419   if(fCreateTree) {
 420     fTreeJetBkg = new TTree("fTreeJetBkg", "fTreeJetBkg");
 421     fTreeJetBkg->Branch("fJet1Vec","TLorentzVector",&fJet1Vec);
 422     fTreeJetBkg->Branch("fJet2Vec","TLorentzVector",&fJet2Vec);
 423     fTreeJetBkg->Branch("fBkgVec","TLorentzVector",&fBkgVec);
 424     fTreeJetBkg->Branch("fArea",&fArea,"fArea/F");
 425     fTreeJetBkg->Branch("fAreaPhi",&fAreaPhi,"fAreaPhi/F");
 426     fTreeJetBkg->Branch("fAreaEta",&fAreaEta,"fAreaEta/F");
 427     fTreeJetBkg->Branch("fRho",&fRho,"fRho/F");
 428     fTreeJetBkg->Branch("fRhoM",&fRhoM,"fRhoM/F");
 429     fTreeJetBkg->Branch("fNConst",&fNConst,"fNConst/I");
 430     fTreeJetBkg->Branch("fJetMassMassless",&fJetMassMassless,"fJetMassMassless/F");
 431   }
 432
 433   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 434   if(fCreateTree) PostData(2, fTreeJetBkg);
 435 }
 436
 437 //________________________________________________________________________
 438 Bool_t AliAnalysisTaskEmcalJetMassResponse::Run()
 439 {
 440   // Run analysis code here, if needed. It will be executed before FillHistograms().
 441
 442   return kTRUE;
 443 }
 444
 445 //________________________________________________________________________
 446 Bool_t AliAnalysisTaskEmcalJetMassResponse::FillHistograms()
 447 {
 448   // Fill histograms.
 449
 450   AliEmcalJet* jet1 = NULL;
 451   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 452   if(jetCont) {
 453     jetCont->ResetCurrentID();
 454     while((jet1 = jetCont->GetNextAcceptJet())) {
 455       if(fSingleTrackEmb) {
 456         AliVParticle *vp = GetEmbeddedConstituent(jet1);
 457         if(!vp) continue;
 458         fJet2Vec->SetPxPyPzE(vp->Px(),vp->Py(),vp->Pz(),vp->E());
 459       } else {
 460         AliEmcalJet *jet2 = jet1->ClosestJet();
 461         if(!jet2) continue;
 462
 463         Double_t fraction = jetCont->GetFractionSharedPt(jet1);
 464         if(fMinFractionShared>0. && fraction<fMinFractionShared) continue;
 465         fJet2Vec->SetPxPyPzE(jet2->Px(),jet2->Py(),jet2->Pz(),jet2->E());
 466       }
 467
 468       Double_t ptJet1 = jet1->Pt() - GetRhoVal(fContainerBase)*jet1->Area();
 469       Double_t massJet1 = GetJetMass(jet1);//jet1->M();
 470
 471       Double_t deltaPt = ptJet1 - fJet2Vec->Pt();
 472       Double_t deltaM  = massJet1 - fJet2Vec->M();
 473       TLorentzVector vpS = GetSubtractedVector(jet1);
 474       Double_t deltaE  = vpS.E() - fJet2Vec->E();
 475       Double_t deltaP  = vpS.P() - fJet2Vec->P();
 476
 477       fh2PtJet1DeltaMNoSub[fCentBin]->Fill(ptJet1,jet1->M()-fJet2Vec->M());
 478       fh2PtJet2DeltaMNoSub[fCentBin]->Fill(fJet2Vec->Pt(),jet1->M()-fJet2Vec->M());
 479
 480       fh3PtJet1DeltaPtDeltaM[fCentBin]->Fill(ptJet1,deltaPt,deltaM);
 481       fh3PtJet2DeltaPtDeltaM[fCentBin]->Fill(fJet2Vec->Pt(),deltaPt,deltaM);
 482
 483       fh2PtJet1DeltaE[fCentBin]->Fill(ptJet1,deltaE);
 484       //  fh2PtJet2DeltaE[fCentBin]->Fill(ptJet2,deltaE);
 485
 486       fh2PtJet1DeltaP[fCentBin]->Fill(ptJet1,deltaP);
 487       //      fh2PtJet2DeltaP[fCentBin]->Fill(ptJet2,deltaP);
 488
 489       fh3PtJet1MJet1MJet2[fCentBin]->Fill(ptJet1,massJet1,fJet2Vec->M());
 490       fh3PtJet2MJet1MJet2[fCentBin]->Fill(fJet2Vec->Pt(),massJet1,fJet2Vec->M());
 491
 492       TLorentzVector vpBkgCheat = GetBkgVectorCheat(jet1);
 493       //      TLorentzVector vpBkgCheat; vpBkgCheat.SetPtEtaPhiM(vpBkgCheatB.Perp(),jet1->Eta(),jet1->Phi(),vpBkgCheatB.M());
 494       TLorentzVector vpBkg = GetBkgVector(jet1,jetCont);
 495       fh2PtJet2DeltaE[fCentBin]->Fill(fJet2Vec->Pt(),vpBkg.E()-vpBkgCheat.E());
 496       fh2PtJet2DeltaP[fCentBin]->Fill(fJet2Vec->Pt(),vpBkg.P()-vpBkgCheat.P());
 497       fh2PtJet2DeltaM[fCentBin]->Fill(fJet2Vec->Pt(),vpBkg.M()-vpBkgCheat.M());
 498
 499       // fh3PtJet1DeltaPtDeltaMRho[fCentBin]->Fill(jet1->Pt()-vpBkg.Perp()*jet1->Area(),jet1->Pt()-vpBkg.Perp()*jet1->Area()-ptJet2,jet1->M()-vpBkg.M()*jet1->Area()-jet2->M());
 500       // fh3PtJet2DeltaPtDeltaMRho[fCentBin]->Fill(ptJet2,jet1->Pt()-vpBkg.Perp()*jet1->Area()-ptJet2,jet1->M()-vpBkg.M()*jet1->Area()-jet2->M());
 501       Double_t px = jet1->Px()-vpBkg.Px();
 502       Double_t py = jet1->Py()-vpBkg.Py();
 503       Double_t pz = jet1->Pz()-vpBkg.Pz();
 504       Double_t E  = jet1->E()-vpBkg.E();
 505       Double_t p2 = px*px + py*py + pz*pz;
 506       Double_t msub = 0.;
 507       if((E*E)>p2) msub = TMath::Sqrt(E*E - p2);
 508       // fh3PtJet1DeltaPtDeltaMRho[fCentBin]->Fill(jet1->Pt()-vpBkg.Perp(),jet1->Pt()-vpBkg.Perp()-ptJet2,jet1->M()-vpBkg.M()-jet2->M());
 509       // fh3PtJet2DeltaPtDeltaMRho[fCentBin]->Fill(ptJet2,jet1->Pt()-vpBkg.Perp()-ptJet2,jet1->M()-vpBkg.M()-jet2->M());
 510
 511       fh3PtJet1DeltaPtDeltaMRho[fCentBin]->Fill(jet1->Pt()-vpBkg.Perp(),jet1->Pt()-vpBkg.Perp()-fJet2Vec->Pt(),msub-fJet2Vec->M());
 512       fh2PtJet1DeltaERho[fCentBin]->Fill(ptJet1,E - fJet2Vec->E());
 513
 514       fh3PtJet2DeltaPtDeltaMRho[fCentBin]->Fill(fJet2Vec->Pt(),jet1->Pt()-vpBkg.Perp()-fJet2Vec->Pt(),msub-fJet2Vec->M());
 515       fh2PtJet2DeltaPxRho[fCentBin]->Fill(fJet2Vec->Pt(),px - fJet2Vec->Px());
 516       fh2PtJet2DeltaPyRho[fCentBin]->Fill(fJet2Vec->Pt(),py - fJet2Vec->Py());
 517       fh2PtJet2DeltaPzRho[fCentBin]->Fill(fJet2Vec->Pt(),pz - fJet2Vec->Pz());
 518       fh2PtJet2DeltaERho[fCentBin]->Fill(fJet2Vec->Pt(),E - fJet2Vec->E());
 519       fh3PtJet2DeltaEDeltaMRho[fCentBin]->Fill(fJet2Vec->Pt(),E-fJet2Vec->E(),msub-fJet2Vec->M());
 520       fh3PtJet2DeltaPDeltaMRho[fCentBin]->Fill(fJet2Vec->Pt(),TMath::Sqrt(p2)-fJet2Vec->P(),msub-fJet2Vec->M());
 521
 522       fh2PtJet2DeltaMRho[fCentBin]->Fill(fJet2Vec->Pt(),msub - fJet2Vec->M());
 523       fh2PtJet2DeltaPtRho[fCentBin]->Fill(fJet2Vec->Pt(),TMath::Sqrt(px*px+py*py) - fJet2Vec->Pt());
 524
 525       TLorentzVector vpC = GetSubtractedVectorCheat(jet1);
 526       fh3PtJet1DeltaPtDeltaMCheat[fCentBin]->Fill(vpC.Perp(),vpC.Perp()-fJet2Vec->Pt(),vpC.M()-fJet2Vec->M());
 527       fh3PtJet2DeltaPtDeltaMCheat[fCentBin]->Fill(fJet2Vec->Pt(),vpC.Perp()-fJet2Vec->Pt(),vpC.M()-fJet2Vec->M());
 528
 529       if(fJet2Vec->Pt()>20. && fCreateTree) {
 530         fBkgVec->SetPxPyPzE(vpBkg.Px(),vpBkg.Py(),vpBkg.Pz(),vpBkg.E());
 531         fJet1Vec->SetPxPyPzE(px,py,pz,E);                                     //AA jet
 532         fArea = (Float_t)jet1->Area();
 533         fAreaPhi = (Float_t)jet1->AreaPhi();
 534         fAreaEta = (Float_t)jet1->AreaEta();
 535         fRho  = (Float_t)jetCont->GetRhoVal();
 536         fRhoM = (Float_t)jetCont->GetRhoMassVal();
 537         fNConst = (Int_t)jet1->GetNumberOfTracks();
 538         fJetMassMassless = (Float_t)GetJetMassMasslessConstituents(jet1);
 539         fTreeJetBkg->Fill();
 540       }
 541     }
 542   }
 543   return kTRUE;
 544 }
 545
 546 //________________________________________________________________________
 547 TLorentzVector AliAnalysisTaskEmcalJetMassResponse::GetBkgVector(AliEmcalJet *jet, AliJetContainer *cont) {
 548   //get background vector
 549
 550   Double_t rho  = cont->GetRhoVal();
 551   Double_t rhom = cont->GetRhoMassVal();
 552   TLorentzVector vpB;
 553   Double_t aphi = jet->AreaPhi();
 554   Double_t aeta = jet->AreaEta();
 555   //  vpB.SetPxPyPzE(rho*TMath::Cos(aphi)*jet->Area(),rho*TMath::Sin(aphi)*jet->Area(),(rho+rhom)*TMath::SinH(aeta)*jet->Area(),(rho+rhom)*TMath::CosH(aeta)*jet->Area());
 556   vpB.SetPxPyPzE(rho*TMath::Cos(aphi)*jet->Area(),rho*TMath::Sin(aphi)*jet->Area(),(rho+rhom)*TMath::SinH(aeta)*jet->Area(),(rho+rhom)*TMath::CosH(aeta)*jet->Area());
 557   return vpB;
 558 }
 559
 560 //________________________________________________________________________
 561 TLorentzVector AliAnalysisTaskEmcalJetMassResponse::GetSubtractedVector(AliEmcalJet *jet) {
 562   //get subtracted vector
 563   TLorentzVector vpS;
 564
 565   // if(f1JetMassAvg) {
 566   //   Double_t pt = jet->Pt() - GetRhoVal(fContainerBase)*jet->Area();
 567   //   TLorentzVector vpB; vpB.SetPtEtaPhiE(GetRhoVal(fContainerBase)*jet->Area(),0.,0.,f1JetMassAvg->Eval(pt));
 568   //   TLorentzVector vpAAboost; vpAAboost.SetPtEtaPhiM(jet->Pt(),0.,0.,jet->M());
 569   //   TLorentzVector vpSboost = vpAAboost - vpB;
 570   //   vpS.SetPtEtaPhiM(vpSboost.Perp(),jet->Eta(),jet->Phi(),vpSboost.M());
 571   // } else {
 572   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 573   TLorentzVector vpBkg = GetBkgVector(jet,jetCont);
 574   vpS.SetPxPyPzE(jet->Px()-vpBkg.Px(),jet->Py()-vpBkg.Py(),jet->Pz()-vpBkg.Pz(),jet->E()-vpBkg.E());
 575   // TLorentzVector vpAAboost; vpAAboost.SetPtEtaPhiE(jet->Pt(),0.,0.,jet->E());
 576   // TLorentzVector vpBkgboost; vpBkgboost.SetPtEtaPhiE(vpBkg.Perp(),0.,0.,vpBkg.E());
 577   // TLorentzVector vpSboost = vpAAboost - vpBkgboost;
 578   //  vpS.SetPtEtaPhiM(vpSboost.Perp(),jet->Eta(),jet->Phi(),vpSboost.M());
 579   // }
 580   return vpS;
 581 }
 582
 583 //________________________________________________________________________
 584 TLorentzVector AliAnalysisTaskEmcalJetMassResponse::GetBkgVectorCheat(AliEmcalJet *jet) {
 585   //get background vector with cheating
 586   TLorentzVector vpB;
 587   if(fJet2Vec) {
 588     TLorentzVector vpAAboost; vpAAboost.SetPtEtaPhiM(jet->Pt(),0.,0.,jet->M());
 589     TLorentzVector vpPPboost; vpPPboost.SetPtEtaPhiM(fJet2Vec->Pt(),0.,0.,fJet2Vec->M());
 590     Double_t dpt = vpAAboost.Perp()-vpPPboost.Perp();
 591     Double_t dE = vpAAboost.E()-vpPPboost.E();
 592     vpB.SetPtEtaPhiE(dpt,0.,0.,dE);
 593   }
 594   return vpB;
 595 }
 596
 597 //________________________________________________________________________
 598 TLorentzVector AliAnalysisTaskEmcalJetMassResponse::GetSubtractedVectorCheat(AliEmcalJet *jet) {
 599   //get subtracted vector taking pT and energy difference from MC match
 600   TLorentzVector vpS;
 601   TLorentzVector vpB = GetBkgVectorCheat(jet);
 602   TLorentzVector vpAAboost; vpAAboost.SetPtEtaPhiM(jet->Pt(),0.,0.,jet->M());
 603   TLorentzVector vpSboost = vpAAboost - vpB;
 604   vpS.SetPtEtaPhiM(vpSboost.Perp(),jet->Eta(),jet->Phi(),vpSboost.M());
 605   return vpS;
 606 }
 607
 608 //________________________________________________________________________
 609 Double_t AliAnalysisTaskEmcalJetMassResponse::GetJetMass(AliEmcalJet *jet) {
 610
 611   TLorentzVector vpS = GetSubtractedVector(jet);
 612
 613   AliEmcalJet *jet2 = jet->ClosestJet();
 614   if(jet2) fh2PtJet1DeltaPtVecSub[fCentBin]->Fill(vpS.Perp(),vpS.Perp()-jet2->Pt());
 615
 616   if(fSubtractMassless) {
 617     Double_t m = vpS.M()-GetJetMassMasslessConstituents(jet);
 618     return m;
 619   } else
 620     return vpS.M();
 621 }
 622
 623 //________________________________________________________________________
 624 Double_t AliAnalysisTaskEmcalJetMassResponse::GetJetMassMasslessConstituents(AliEmcalJet *jet) {
 625   //Get mass of jet assuming all particles are massless (E==P)
 626
 627   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 628   AliVParticle *vp = 0x0;
 629   Double_t px = 0.; Double_t py = 0.; Double_t pz = 0.; Double_t E = 0.;
 630   for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
 631     vp = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
 632     px+=vp->Px();
 633     py+=vp->Py();
 634     pz+=vp->Pz();
 635     E+=vp->P();
 636   }
 637
 638   Double_t m2 = E*E - px*px - py*py - pz*pz;
 639   if(m2>0.)
 640     return TMath::Sqrt(E*E - px*px - py*py - pz*pz);
 641   else
 642     return 0.;
 643 }
 644
 645 //________________________________________________________________________
 646 AliVParticle* AliAnalysisTaskEmcalJetMassResponse::GetEmbeddedConstituent(AliEmcalJet *jet) {
 647
 648   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 649   AliVParticle *vp = 0x0;
 650   AliVParticle *vpe = 0x0; //embedded particle
 651   Int_t nc = 0;
 652   for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
 653     vp = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); //check if fTracks is the correct track branch
 654     if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
 655     if(!vpe) vpe = vp;
 656     else if(vp->Pt()>vpe->Pt()) vpe = vp;
 657     nc++;
 658   }
 659   AliDebug(11,Form("Found %d embedded particles",nc));
 660   return vpe;
 661 }
 662
 663 //________________________________________________________________________
 664 Bool_t AliAnalysisTaskEmcalJetMassResponse::RetrieveEventObjects() {
 665   //
 666   // retrieve event objects
 667   //
 668
 669   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 670     return kFALSE;
 671
 672   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 673   jetCont->LoadRhoMass(InputEvent());
 674
 675   return kTRUE;
 676 }
 677
 678 //_______________________________________________________________________
 679 void AliAnalysisTaskEmcalJetMassResponse::Terminate(Option_t *)
 680 {
 681   // Called once at the end of the analysis.
 682 }
 683