PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetMassBkg.cxx

   1 //
   2 // Jet mass background 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 <TList.h>
  12 #include <TLorentzVector.h>
  13 #include <TProfile.h>
  14 #include <TChain.h>
  15 #include <TSystem.h>
  16 #include <TFile.h>
  17 #include <TKey.h>
  18 #include <TRandom3.h>
  19
  20 #include "AliVCluster.h"
  21 #include "AliVTrack.h"
  22 #include "AliEmcalJet.h"
  23 #include "AliRhoParameter.h"
  24 #include "AliLog.h"
  25 #include "AliEmcalParticle.h"
  26 #include "AliMCEvent.h"
  27 #include "AliGenPythiaEventHeader.h"
  28 #include "AliAODMCHeader.h"
  29 #include "AliMCEvent.h"
  30 #include "AliAnalysisManager.h"
  31 #include "AliJetContainer.h"
  32 #include "AliClusterContainer.h"
  33 #include "AliParticleContainer.h"
  34
  35 #include "AliAnalysisTaskEmcalJetMassBkg.h"
  36
  37 ClassImp(AliAnalysisTaskEmcalJetMassBkg)
  38
  39 //________________________________________________________________________
  40 AliAnalysisTaskEmcalJetMassBkg::AliAnalysisTaskEmcalJetMassBkg() :
  41   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetMassBkg", kTRUE),
  42   fContainerBase(0),
  43   fMinRC2LJ(-1),
  44   fRCperEvent(10),
  45   fConeRadius(0.2),
  46   fConeMinEta(-0.9),
  47   fConeMaxEta(0.9),
  48   fConeMinPhi(0),
  49   fConeMaxPhi(TMath::Pi()*2),
  50   fJetsCont(0),
  51   fTracksCont(0),
  52   fCaloClustersCont(0),
  53   fh2PtVsMassRC(0),
  54   fh2PtVsMassRCExLJDPhi(0),
  55   fh2PtVsMassPerpConeLJ(0),
  56   fh2PtVsMassPerpConeTJ(0),
  57   fh2PtVsERC(0),
  58   fh2PtVsERCExLJDPhi(0),
  59   fh2PtVsEPerpConeLJ(0),
  60   fh2PtVsEPerpConeTJ(0),
  61   fpPtVsMassRC(0),
  62   fpPtVsMassRCExLJ(0),
  63   fpPtVsMassPerpConeLJ(0),
  64   fpPtVsMassPerpConeTJ(0),
  65   fh2EtaVsMassRC(0),
  66   fh2EtaVsMassRCExLJ(0),
  67   fh2EtaVsMassPerpConeLJ(0),
  68   fh2EtaVsMassPerpConeTJ(0),
  69   fh2CentVsMassRC(0),
  70   fh2CentVsMassRCExLJ(0),
  71   fh2CentVsMassPerpConeLJ(0),
  72   fh2CentVsMassPerpConeTJ(0),
  73   fh2MultVsMassRC(0),
  74   fh2MultVsMassRCExLJ(0),
  75   fh2MultVsMassPerpConeLJ(0),
  76   fh2MultVsMassPerpConeTJ(0),
  77   fh2CentVsMedianMassRC(0),
  78   fh2CentVsMedianMassRCExLJ(0),
  79   fh2MultVsMedianMassRC(0),
  80   fh2MultVsMedianMassRCExLJ(0),
  81   fh2CentVsMeanMassRC(0),
  82   fh2CentVsMeanMassRCExLJ(0),
  83   fh2MultVsMeanMassRC(0),
  84   fh2MultVsMeanMassRCExLJ(0),
  85   fh2CentVsMedianMassPerAreaRC(0),
  86   fh2CentVsMedianMassPerAreaRCExLJ(0),
  87   fh2MultVsMedianMassPerAreaRC(0),
  88   fh2MultVsMedianMassPerAreaRCExLJ(0)
  89 {
  90   // Default constructor.
  91
  92   fh2PtVsMassRC            = new TH2F*[fNcentBins];
  93   fh2PtVsMassRCExLJDPhi    = new TH3F*[fNcentBins];
  94   fh2PtVsMassPerpConeLJ    = new TH2F*[fNcentBins];
  95   fh2PtVsMassPerpConeTJ    = new TH2F*[fNcentBins];
  96
  97   fh2PtVsERC               = new TH2F*[fNcentBins];
  98   fh2PtVsERCExLJDPhi       = new TH3F*[fNcentBins];
  99   fh2PtVsEPerpConeLJ       = new TH2F*[fNcentBins];
 100   fh2PtVsEPerpConeTJ       = new TH2F*[fNcentBins];
 101
 102   fpPtVsMassRC             = new TProfile*[fNcentBins];
 103   fpPtVsMassRCExLJ         = new TProfile*[fNcentBins];
 104   fpPtVsMassPerpConeLJ     = new TProfile*[fNcentBins];
 105   fpPtVsMassPerpConeTJ     = new TProfile*[fNcentBins];
 106
 107   fh2EtaVsMassRC            = new TH2F*[fNcentBins];
 108   fh2EtaVsMassRCExLJ        = new TH2F*[fNcentBins];
 109   fh2EtaVsMassPerpConeLJ    = new TH2F*[fNcentBins];
 110   fh2EtaVsMassPerpConeTJ    = new TH2F*[fNcentBins];
 111
 112   for (Int_t i = 0; i < fNcentBins; i++) {
 113     fh2PtVsMassRC[i]             = 0;
 114     fh2PtVsMassRCExLJDPhi[i]     = 0;
 115     fh2PtVsMassPerpConeLJ[i]     = 0;
 116     fh2PtVsMassPerpConeTJ[i]     = 0;
 117
 118     fh2PtVsERC[i]                = 0;
 119     fh2PtVsERCExLJDPhi[i]        = 0;
 120     fh2PtVsEPerpConeLJ[i]        = 0;
 121     fh2PtVsEPerpConeTJ[i]        = 0;
 122
 123     fpPtVsMassRC[i]              = 0;
 124     fpPtVsMassRCExLJ[i]          = 0;
 125     fpPtVsMassPerpConeLJ[i]      = 0;
 126     fpPtVsMassPerpConeTJ[i]      = 0;
 127
 128     fh2EtaVsMassRC[i]             = 0;
 129     fh2EtaVsMassRCExLJ[i]         = 0;
 130     fh2EtaVsMassPerpConeLJ[i]     = 0;
 131     fh2EtaVsMassPerpConeTJ[i]     = 0;
 132   }
 133
 134   SetMakeGeneralHistograms(kTRUE);
 135
 136 }
 137
 138 //________________________________________________________________________
 139 AliAnalysisTaskEmcalJetMassBkg::AliAnalysisTaskEmcalJetMassBkg(const char *name) :
 140   AliAnalysisTaskEmcalJet(name, kTRUE),
 141   fContainerBase(0),
 142   fMinRC2LJ(-1),
 143   fRCperEvent(10),
 144   fConeRadius(0.2),
 145   fConeMinEta(-0.9),
 146   fConeMaxEta(0.9),
 147   fConeMinPhi(0),
 148   fConeMaxPhi(TMath::Pi()*2),
 149   fJetsCont(0),
 150   fTracksCont(0),
 151   fCaloClustersCont(0),
 152   fh2PtVsMassRC(0),
 153   fh2PtVsMassRCExLJDPhi(0),
 154   fh2PtVsMassPerpConeLJ(0),
 155   fh2PtVsMassPerpConeTJ(0),
 156   fh2PtVsERC(0),
 157   fh2PtVsERCExLJDPhi(0),
 158   fh2PtVsEPerpConeLJ(0),
 159   fh2PtVsEPerpConeTJ(0),
 160   fpPtVsMassRC(0),
 161   fpPtVsMassRCExLJ(0),
 162   fpPtVsMassPerpConeLJ(0),
 163   fpPtVsMassPerpConeTJ(0),
 164   fh2EtaVsMassRC(0),
 165   fh2EtaVsMassRCExLJ(0),
 166   fh2EtaVsMassPerpConeLJ(0),
 167   fh2EtaVsMassPerpConeTJ(0),
 168   fh2CentVsMassRC(0),
 169   fh2CentVsMassRCExLJ(0),
 170   fh2CentVsMassPerpConeLJ(0),
 171   fh2CentVsMassPerpConeTJ(0),
 172   fh2MultVsMassRC(0),
 173   fh2MultVsMassRCExLJ(0),
 174   fh2MultVsMassPerpConeLJ(0),
 175   fh2MultVsMassPerpConeTJ(0),
 176   fh2CentVsMedianMassRC(0),
 177   fh2CentVsMedianMassRCExLJ(0),
 178   fh2MultVsMedianMassRC(0),
 179   fh2MultVsMedianMassRCExLJ(0),
 180   fh2CentVsMeanMassRC(0),
 181   fh2CentVsMeanMassRCExLJ(0),
 182   fh2MultVsMeanMassRC(0),
 183   fh2MultVsMeanMassRCExLJ(0),
 184   fh2CentVsMedianMassPerAreaRC(0),
 185   fh2CentVsMedianMassPerAreaRCExLJ(0),
 186   fh2MultVsMedianMassPerAreaRC(0),
 187   fh2MultVsMedianMassPerAreaRCExLJ(0)
 188 {
 189   // Standard constructor.
 190
 191   fh2PtVsMassRC            = new TH2F*[fNcentBins];
 192   fh2PtVsMassRCExLJDPhi    = new TH3F*[fNcentBins];
 193   fh2PtVsMassPerpConeLJ    = new TH2F*[fNcentBins];
 194   fh2PtVsMassPerpConeTJ    = new TH2F*[fNcentBins];
 195
 196   fh2PtVsERC               = new TH2F*[fNcentBins];
 197   fh2PtVsERCExLJDPhi       = new TH3F*[fNcentBins];
 198   fh2PtVsEPerpConeLJ       = new TH2F*[fNcentBins];
 199   fh2PtVsEPerpConeTJ       = new TH2F*[fNcentBins];
 200
 201   fpPtVsMassRC             = new TProfile*[fNcentBins];
 202   fpPtVsMassRCExLJ         = new TProfile*[fNcentBins];
 203   fpPtVsMassPerpConeLJ     = new TProfile*[fNcentBins];
 204   fpPtVsMassPerpConeTJ     = new TProfile*[fNcentBins];
 205
 206   fh2EtaVsMassRC            = new TH2F*[fNcentBins];
 207   fh2EtaVsMassRCExLJ        = new TH2F*[fNcentBins];
 208   fh2EtaVsMassPerpConeLJ    = new TH2F*[fNcentBins];
 209   fh2EtaVsMassPerpConeTJ    = new TH2F*[fNcentBins];
 210
 211   for (Int_t i = 0; i < fNcentBins; i++) {
 212     fh2PtVsMassRC[i]             = 0;
 213     fh2PtVsMassRCExLJDPhi[i]     = 0;
 214     fh2PtVsMassPerpConeLJ[i]     = 0;
 215     fh2PtVsMassPerpConeTJ[i]     = 0;
 216
 217     fh2PtVsERC[i]                = 0;
 218     fh2PtVsERCExLJDPhi[i]        = 0;
 219     fh2PtVsEPerpConeLJ[i]        = 0;
 220     fh2PtVsEPerpConeTJ[i]        = 0;
 221
 222     fpPtVsMassRC[i]              = 0;
 223     fpPtVsMassRCExLJ[i]          = 0;
 224     fpPtVsMassPerpConeLJ[i]      = 0;
 225     fpPtVsMassPerpConeTJ[i]      = 0;
 226
 227     fh2EtaVsMassRC[i]             = 0;
 228     fh2EtaVsMassRCExLJ[i]         = 0;
 229     fh2EtaVsMassPerpConeLJ[i]     = 0;
 230     fh2EtaVsMassPerpConeTJ[i]     = 0;
 231   }
 232
 233   SetMakeGeneralHistograms(kTRUE);
 234 }
 235
 236 //________________________________________________________________________
 237 AliAnalysisTaskEmcalJetMassBkg::~AliAnalysisTaskEmcalJetMassBkg()
 238 {
 239   // Destructor.
 240 }
 241
 242 //________________________________________________________________________
 243 void AliAnalysisTaskEmcalJetMassBkg::UserCreateOutputObjects()
 244 {
 245   // Create user output.
 246
 247   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 248
 249   fJetsCont         = GetJetContainer(fContainerBase);
 250   fTracksCont       = fJetsCont->GetParticleContainer();
 251   fCaloClustersCont = fJetsCont->GetClusterContainer();
 252
 253   Bool_t oldStatus = TH1::AddDirectoryStatus();
 254   TH1::AddDirectory(kFALSE);
 255
 256   const Int_t nBinsPt  = 250;
 257   const Double_t minPt = -50.;
 258   const Double_t maxPt = 200.;
 259
 260   const Int_t nBinsE  = 250;
 261   const Double_t minE = -50.;
 262   const Double_t maxE = 200.;
 263
 264   const Int_t nBinsM  = 150;
 265   const Double_t minM = -50.;
 266   const Double_t maxM = 100.;
 267
 268   const Int_t nBinsEta  = 100;
 269   const Double_t minEta = -1.;
 270   const Double_t maxEta =  1.;
 271
 272   const Int_t nBinsCent  = 100;
 273   const Double_t minCent = 0.;
 274   const Double_t maxCent = 100.;
 275
 276   const Int_t nBinsMult  = 400;
 277   const Double_t minMult = 0.;
 278   const Double_t maxMult = 4000.;
 279
 280   fh2CentVsMassRC = new TH2F("fh2CentVsMassRC","fh2CentVsMassRC;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsM,minM,maxM);
 281   fOutput->Add(fh2CentVsMassRC);
 282
 283   fh2CentVsMassRCExLJ = new TH2F("fh2CentVsMassRCExLJ","fh2CentVsMassRCExLJ;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsM,minM,maxM);
 284   fOutput->Add(fh2CentVsMassRCExLJ);
 285
 286   fh2CentVsMassPerpConeLJ = new TH2F("fh2CentVsMassPerpConeLJ","fh2CentVsMassPerpConeLJ;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsM,minM,maxM);
 287   fOutput->Add(fh2CentVsMassPerpConeLJ);
 288
 289   fh2CentVsMassPerpConeTJ = new TH2F("fh2CentVsMassPerpConeTJ","fh2CentVsMassPerpConeTJ;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsM,minM,maxM);
 290   fOutput->Add(fh2CentVsMassPerpConeTJ);
 291
 292   fh2MultVsMassRC = new TH2F("fh2MultVsMassRC","fh2MultVsMassRC;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsM,minM,maxM);
 293   fOutput->Add(fh2MultVsMassRC);
 294
 295   fh2MultVsMassRCExLJ = new TH2F("fh2MultVsMassRCExLJ","fh2MultVsMassRCExLJ;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsM,minM,maxM);
 296   fOutput->Add(fh2MultVsMassRCExLJ);
 297
 298   fh2MultVsMassPerpConeLJ = new TH2F("fh2MultVsMassPerpConeLJ","fh2MultVsMassPerpConeLJ;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsM,minM,maxM);
 299   fOutput->Add(fh2MultVsMassPerpConeLJ);
 300
 301   fh2MultVsMassPerpConeTJ = new TH2F("fh2MultVsMassPerpConeTJ","fh2MultVsMassPerpConeTJ;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsM,minM,maxM);
 302   fOutput->Add(fh2MultVsMassPerpConeTJ);
 303
 304   fh2CentVsMedianMassRC = new TH2F("fh2CentVsMedianMassRC","fh2CentVsMedianMassRC;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsM,minM,maxM);
 305   fOutput->Add(fh2CentVsMedianMassRC);
 306
 307   fh2CentVsMedianMassRCExLJ = new TH2F("fh2CentVsMedianMassRCExLJ","fh2CentVsMedianMassRCExLJ;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsM,minM,maxM);
 308   fOutput->Add(fh2CentVsMedianMassRCExLJ);
 309
 310   fh2MultVsMedianMassRC = new TH2F("fh2MultVsMedianMassRC","fh2MultVsMedianMassRC;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsM,minM,maxM);
 311   fOutput->Add(fh2MultVsMedianMassRC);
 312
 313   fh2MultVsMedianMassRCExLJ = new TH2F("fh2MultVsMedianMassRCExLJ","fh2MultVsMedianMassRCExLJ;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsM,minM,maxM);
 314   fOutput->Add(fh2MultVsMedianMassRCExLJ);
 315
 316   fh2CentVsMeanMassRC = new TH2F("fh2CentVsMeanMassRC","fh2CentVsMeanMassRC;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsM,minM,maxM);
 317   fOutput->Add(fh2CentVsMeanMassRC);
 318
 319   fh2CentVsMeanMassRCExLJ = new TH2F("fh2CentVsMeanMassRCExLJ","fh2CentVsMeanMassRCExLJ;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsM,minM,maxM);
 320   fOutput->Add(fh2CentVsMeanMassRCExLJ);
 321
 322   fh2MultVsMeanMassRC = new TH2F("fh2MultVsMeanMassRC","fh2MultVsMeanMassRC;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsM,minM,maxM);
 323   fOutput->Add(fh2MultVsMeanMassRC);
 324
 325   fh2MultVsMeanMassRCExLJ = new TH2F("fh2MultVsMeanMassRCExLJ","fh2MultVsMeanMassRCExLJ;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsM,minM,maxM);
 326   fOutput->Add(fh2MultVsMeanMassRCExLJ);
 327
 328   fh2CentVsMedianMassPerAreaRC = new TH2F("fh2CentVsMedianMassPerAreaRC","fh2CentVsMedianMassPerAreaRC;cent;#it{M}_{RC}/A",nBinsCent,minCent,maxCent,nBinsM,minM,maxM);
 329   fOutput->Add(fh2CentVsMedianMassPerAreaRC);
 330
 331   fh2CentVsMedianMassPerAreaRCExLJ = new TH2F("fh2CentVsMedianMassPerAreaRCExLJ","fh2CentVsMedianMassPerAreaRCExLJ;cent;#it{M}_{RC}/A",nBinsCent,minCent,maxCent,nBinsM,minM,maxM);
 332   fOutput->Add(fh2CentVsMedianMassPerAreaRCExLJ);
 333
 334   fh2MultVsMedianMassPerAreaRC = new TH2F("fh2MultVsMedianMassPerAreaRC","fh2MultVsMedianMassPerAreaRC;#it{N}_{track};#it{M}_{RC}/A",nBinsMult,minMult,maxMult,nBinsM,minM,maxM);
 335   fOutput->Add(fh2MultVsMedianMassPerAreaRC);
 336
 337   fh2MultVsMedianMassPerAreaRCExLJ = new TH2F("fh2MultVsMedianMassPerAreaRCExLJ","fh2MultVsMedianMassPerAreaRCExLJ;#it{N}_{track};#it{M}_{RC}/A",nBinsMult,minMult,maxMult,nBinsM,minM,maxM);
 338   fOutput->Add(fh2MultVsMedianMassPerAreaRCExLJ);
 339
 340   TString histName = "";
 341   TString histTitle = "";
 342   for (Int_t i = 0; i < fNcentBins; i++) {
 343     histName = TString::Format("fh2PtVsMassRC_%d",i);
 344     histTitle = TString::Format("%s;#it{p}_{T,RC};#it{M}_{RC}",histName.Data());
 345     fh2PtVsMassRC[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 346     fOutput->Add(fh2PtVsMassRC[i]);
 347
 348     histName = TString::Format("fh2PtVsMassRCExLJDPhi_%d",i);
 349     histTitle = TString::Format("%s;#it{p}_{T,RC};#it{M}_{RC}",histName.Data());
 350     fh2PtVsMassRCExLJDPhi[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,72,-0.5*TMath::Pi(),1.5*TMath::Pi());
 351     fOutput->Add(fh2PtVsMassRCExLJDPhi[i]);
 352
 353     histName = TString::Format("fh2PtVsMassPerpConeLJ_%d",i);
 354     histTitle = TString::Format("%s;#it{p}_{T,PerpConeLJ};#it{M}_{PerpConeLJ}",histName.Data());
 355     fh2PtVsMassPerpConeLJ[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 356     fOutput->Add(fh2PtVsMassPerpConeLJ[i]);
 357
 358     histName = TString::Format("fh2PtVsMassPerpConeTJ_%d",i);
 359     histTitle = TString::Format("%s;#it{p}_{T,PerpConeTJ};#it{M}_{PerpConeTJ}",histName.Data());
 360     fh2PtVsMassPerpConeTJ[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 361     fOutput->Add(fh2PtVsMassPerpConeTJ[i]);
 362     //
 363     histName = TString::Format("fh2PtVsERC_%d",i);
 364     histTitle = TString::Format("%s;#it{p}_{T,RC};#it{M}_{RC}",histName.Data());
 365     fh2PtVsERC[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsE,minE,maxE);
 366     fOutput->Add(fh2PtVsERC[i]);
 367
 368     histName = TString::Format("fh2PtVsERCExLJDPhi_%d",i);
 369     histTitle = TString::Format("%s;#it{p}_{T,RC};#it{M}_{RC}",histName.Data());
 370     fh2PtVsERCExLJDPhi[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,72,-0.5*TMath::Pi(),1.5*TMath::Pi());
 371     fOutput->Add(fh2PtVsERCExLJDPhi[i]);
 372
 373     histName = TString::Format("fh2PtVsEPerpConeLJ_%d",i);
 374     histTitle = TString::Format("%s;#it{p}_{T,PerpConeLJ};#it{M}_{PerpConeLJ}",histName.Data());
 375     fh2PtVsEPerpConeLJ[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsE,minE,maxE);
 376     fOutput->Add(fh2PtVsEPerpConeLJ[i]);
 377
 378     histName = TString::Format("fh2PtVsEPerpConeTJ_%d",i);
 379     histTitle = TString::Format("%s;#it{p}_{T,PerpConeTJ};#it{M}_{PerpConeTJ}",histName.Data());
 380     fh2PtVsEPerpConeTJ[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsE,minE,maxE);
 381     fOutput->Add(fh2PtVsEPerpConeTJ[i]);
 382     //
 383     histName = TString::Format("fh2EtaVsMassRC_%d",i);
 384     histTitle = TString::Format("%s;#eta_{RC};#it{M}_{RC}",histName.Data());
 385     fh2EtaVsMassRC[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsEta,minEta,maxEta,nBinsM,minM,maxM);
 386     fOutput->Add(fh2EtaVsMassRC[i]);
 387
 388     histName = TString::Format("fh2EtaVsMassRCExLJ_%d",i);
 389     histTitle = TString::Format("%s;#eta_{RC};#it{M}_{RC}",histName.Data());
 390     fh2EtaVsMassRCExLJ[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsEta,minEta,maxEta,nBinsM,minM,maxM);
 391     fOutput->Add(fh2EtaVsMassRCExLJ[i]);
 392
 393     histName = TString::Format("fh2EtaVsMassPerpConeLJ_%d",i);
 394     histTitle = TString::Format("%s;#eta_{PerpConeLJ};#it{M}_{PerpConeLJ}",histName.Data());
 395     fh2EtaVsMassPerpConeLJ[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsEta,minEta,maxEta,nBinsM,minM,maxM);
 396     fOutput->Add(fh2EtaVsMassPerpConeLJ[i]);
 397
 398     histName = TString::Format("fh2EtaVsMassPerpConeTJ_%d",i);
 399     histTitle = TString::Format("%s;#eta_{PerpConeTJ};#it{M}_{PerpConeTJ}",histName.Data());
 400     fh2EtaVsMassPerpConeTJ[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsEta,minEta,maxEta,nBinsM,minM,maxM);
 401     fOutput->Add(fh2EtaVsMassPerpConeTJ[i]);
 402
 403     histName = TString::Format("fpPtVsMassRC_%d",i);
 404     histTitle = TString::Format("%s;#it{p}_{T,RC};Avg #it{M}_{RC}",histName.Data());
 405     fpPtVsMassRC[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
 406     fOutput->Add(fpPtVsMassRC[i]);
 407
 408     histName = TString::Format("fpPtVsMassRCExLJ_%d",i);
 409     histTitle = TString::Format("%s;#it{p}_{T,RC};Avg #it{M}_{RC}",histName.Data());
 410     fpPtVsMassRCExLJ[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
 411     fOutput->Add(fpPtVsMassRCExLJ[i]);
 412
 413     histName = TString::Format("fpPtVsMassPerpConeLJ_%d",i);
 414     histTitle = TString::Format("%s;#it{p}_{T,RC};Avg #it{M}_{RC}",histName.Data());
 415     fpPtVsMassPerpConeLJ[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
 416     fOutput->Add(fpPtVsMassPerpConeLJ[i]);
 417
 418     histName = TString::Format("fpPtVsMassPerpConeTJ_%d",i);
 419     histTitle = TString::Format("%s;#it{p}_{T,RC};Avg #it{M}_{RC}",histName.Data());
 420     fpPtVsMassPerpConeTJ[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
 421     fOutput->Add(fpPtVsMassPerpConeTJ[i]);
 422   }
 423
 424   // =========== Switch on Sumw2 for all histos ===========
 425   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 426     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 427     if (h1){
 428       h1->Sumw2();
 429       continue;
 430     }
 431     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 432     if(hn)hn->Sumw2();
 433   }
 434
 435   TH1::AddDirectory(oldStatus);
 436
 437   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 438 }
 439
 440 //________________________________________________________________________
 441 Bool_t AliAnalysisTaskEmcalJetMassBkg::Run()
 442 {
 443   // Run analysis code here, if needed. It will be executed before FillHistograms().
 444
 445   return kTRUE;
 446 }
 447
 448 //________________________________________________________________________
 449 Bool_t AliAnalysisTaskEmcalJetMassBkg::FillHistograms()
 450 {
 451   // Fill histograms.
 452
 453   const Float_t rcArea = fConeRadius * fConeRadius * TMath::Pi();
 454
 455   //Event properties
 456   Double_t rho = GetRhoVal(fContainerBase);
 457   Int_t trackMult = fTracksCont->GetNAcceptedParticles();
 458
 459   //Leading jet
 460   AliEmcalJet* jet = NULL;
 461   if (fJetsCont)
 462     jet = fJetsCont->GetLeadingJet("rho");
 463
 464   TLorentzVector lvRC(0.,0.,0.,0.);
 465   Float_t RCpt = 0;
 466   Float_t RCeta = 0;
 467   Float_t RCphi = 0;
 468   Float_t RCmass = 0.;
 469   Float_t RCE = 0.;
 470
 471   static Double_t massvecRC[999];
 472   static Double_t massPerAreavecRC[999];
 473
 474   static Double_t massvecRCExLJ[999];
 475   static Double_t massPerAreavecRCExLJ[999];
 476
 477   Int_t nRCAcc = 0;
 478   Int_t nRCExLJAcc = 0;
 479
 480   for (Int_t i = 0; i < fRCperEvent; i++) {
 481     // Simple random cones
 482     lvRC.SetPxPyPzE(0.,0.,0.,0.);
 483     RCpt = 0;
 484     RCeta = 0;
 485     RCphi = 0;
 486     GetRandomCone(lvRC,RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, 0);
 487     RCmass = lvRC.M();
 488     RCE = lvRC.E();
 489     if (RCpt > 0) {
 490       fh2PtVsMassRC[fCentBin]->Fill(RCpt - rho*rcArea,RCmass);
 491       fh2PtVsERC[fCentBin]->Fill(RCpt - rho*rcArea,RCE);
 492       fpPtVsMassRC[fCentBin]->Fill(RCpt - rho*rcArea,RCmass);
 493       fh2EtaVsMassRC[fCentBin]->Fill(RCeta,RCmass);
 494       fh2CentVsMassRC->Fill(fCent,RCmass);
 495       fh2MultVsMassRC->Fill(trackMult,RCmass);
 496
 497       massvecRC[nRCAcc] = RCmass;
 498       massPerAreavecRC[nRCAcc] = RCmass/rcArea;
 499       ++nRCAcc;
 500     }
 501
 502     if (fJetsCont && jet) {
 503       // Random cones far away from leading jet(s)
 504       lvRC.SetPxPyPzE(0.,0.,0.,0.);
 505       RCpt = 0;
 506       RCeta = 0;
 507       RCphi = 0;
 508       GetRandomCone(lvRC,RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet);
 509       RCmass = lvRC.M();
 510       RCE = lvRC.E();
 511       if (RCpt > 0 && jet) {
 512         Float_t dphi = RCphi - jet->Phi();
 513         if (dphi > 1.5*TMath::Pi()) dphi -= TMath::Pi() * 2;
 514         if (dphi < -0.5*TMath::Pi()) dphi += TMath::Pi() * 2;
 515         fh2PtVsMassRCExLJDPhi[fCentBin]->Fill(RCpt - rho*rcArea,RCmass,dphi);
 516         fh2PtVsERCExLJDPhi[fCentBin]->Fill(RCpt - rho*rcArea,RCE,dphi);
 517         fpPtVsMassRCExLJ[fCentBin]->Fill(RCpt - rho*rcArea,RCmass);
 518         fh2EtaVsMassRCExLJ[fCentBin]->Fill(RCeta,RCmass);
 519         fh2CentVsMassRCExLJ->Fill(fCent,RCmass);
 520         fh2MultVsMassRCExLJ->Fill(trackMult,RCmass);
 521
 522         massvecRCExLJ[nRCExLJAcc] = RCmass;
 523         massPerAreavecRCExLJ[nRCExLJAcc] = RCmass/rcArea;
 524         ++nRCExLJAcc;
 525       }
 526     }
 527   }//RC loop
 528
 529   Double_t medianRC, medianRCExLJ = 0.;
 530   medianRC = TMath::Median(nRCAcc,massvecRC);
 531   medianRCExLJ = TMath::Median(nRCExLJAcc,massvecRCExLJ);
 532
 533   fh2CentVsMedianMassRC->Fill(fCent,medianRC);
 534   fh2CentVsMedianMassRCExLJ->Fill(fCent,medianRCExLJ);
 535
 536   fh2MultVsMedianMassRC->Fill(trackMult,medianRC);
 537   fh2MultVsMedianMassRCExLJ->Fill(trackMult,medianRCExLJ);
 538
 539   Double_t meanRC = 0.; Double_t meanRCExLJ = 0.;
 540   if(nRCAcc>0)   meanRC = TMath::Mean(nRCAcc,massvecRC);
 541   if(nRCExLJAcc) meanRCExLJ = TMath::Mean(nRCExLJAcc,massvecRCExLJ);
 542
 543   fh2CentVsMeanMassRC->Fill(fCent,meanRC);
 544   fh2CentVsMeanMassRCExLJ->Fill(fCent,meanRCExLJ);
 545
 546   fh2MultVsMeanMassRC->Fill(trackMult,meanRC);
 547   fh2MultVsMeanMassRCExLJ->Fill(trackMult,meanRCExLJ);
 548
 549   Double_t medianPerAreaRC, medianPerAreaRCExLJ = 0.;
 550   medianPerAreaRC = TMath::Median(nRCAcc,massPerAreavecRC);
 551   medianPerAreaRCExLJ = TMath::Median(nRCExLJAcc,massPerAreavecRCExLJ);
 552
 553   fh2CentVsMedianMassPerAreaRC->Fill(fCent,medianPerAreaRC);
 554   fh2CentVsMedianMassPerAreaRCExLJ->Fill(fCent,medianPerAreaRCExLJ);
 555
 556   fh2MultVsMedianMassPerAreaRC->Fill(trackMult,medianPerAreaRC);
 557   fh2MultVsMedianMassPerAreaRCExLJ->Fill(trackMult,medianPerAreaRCExLJ);
 558
 559
 560   if(fJetsCont && jet) {
 561     //cone perpendicular to leading jet
 562     TLorentzVector lvPC(0.,0.,0.,0.);
 563     Float_t PCpt = 0;
 564     Float_t PCeta = 0;
 565     Float_t PCphi = 0;
 566     Float_t PCmass = 0.;
 567     Float_t PCE = 0.;
 568     if(jet) {
 569       GetPerpCone(lvPC,PCpt, PCeta, PCphi, fTracksCont, fCaloClustersCont, jet);
 570       PCmass = lvPC.M();
 571       PCE = lvPC.E();
 572       if(PCpt>0.) {
 573         fh2PtVsMassPerpConeLJ[fCentBin]->Fill(PCpt-rho*rcArea,PCmass);
 574         fh2PtVsEPerpConeLJ[fCentBin]->Fill(PCpt-rho*rcArea,PCE);
 575         fpPtVsMassPerpConeLJ[fCentBin]->Fill(PCpt-rho*rcArea,PCmass);
 576         fh2EtaVsMassPerpConeLJ[fCentBin]->Fill(PCeta,PCmass);
 577         fh2CentVsMassPerpConeLJ->Fill(fCent,PCmass);
 578         fh2MultVsMassPerpConeLJ->Fill(trackMult,PCmass);
 579       }
 580     }
 581     //cone perpendicular to all tagged jets
 582     for(int i = 0; i < fJetsCont->GetNJets();++i) {
 583       jet = static_cast<AliEmcalJet*>(fJetsCont->GetAcceptJet(i));
 584       if(!jet) continue;
 585
 586       if(jet->GetTagStatus()<1)
 587         continue;
 588
 589       lvPC.SetPxPyPzE(0.,0.,0.,0.);
 590       PCpt = 0;
 591       PCeta = 0;
 592       PCphi = 0;
 593       GetPerpCone(lvPC,PCpt, PCeta, PCphi, fTracksCont, fCaloClustersCont, jet);
 594       PCmass = lvPC.M();
 595       PCE = lvPC.E();
 596       if(PCpt>0.) {
 597         fh2PtVsMassPerpConeTJ[fCentBin]->Fill(PCpt-rho*rcArea,PCmass);
 598         fh2PtVsEPerpConeTJ[fCentBin]->Fill(PCpt-rho*rcArea,PCE);
 599         fpPtVsMassPerpConeTJ[fCentBin]->Fill(PCpt-rho*rcArea,PCmass);
 600         fh2EtaVsMassPerpConeTJ[fCentBin]->Fill(PCeta,PCmass);
 601         fh2CentVsMassPerpConeTJ->Fill(fCent,PCmass);
 602         fh2MultVsMassPerpConeTJ->Fill(trackMult,PCmass);
 603       }
 604     }//jet loop
 605   }
 606
 607   return kTRUE;
 608
 609 }
 610
 611 //________________________________________________________________________
 612 void AliAnalysisTaskEmcalJetMassBkg::GetRandomCone(TLorentzVector& lvRC,Float_t &pt, Float_t &eta, Float_t &phi,
 613                                                    AliParticleContainer* tracks, AliClusterContainer* clusters,
 614                                                    AliEmcalJet *jet) const
 615 {
 616   // Get rigid cone.
 617   lvRC.SetPxPyPzE(0.,0.,0.,0.);
 618
 619   eta = -999;
 620   phi = -999;
 621   pt = 0;
 622
 623   if (!tracks && !clusters)
 624     return;
 625
 626   Float_t LJeta = 999;
 627   Float_t LJphi = 999;
 628
 629   if (jet) {
 630     LJeta = jet->Eta();
 631     LJphi = jet->Phi();
 632   }
 633
 634   Float_t maxEta = fConeMaxEta;
 635   Float_t minEta = fConeMinEta;
 636   Float_t maxPhi = fConeMaxPhi;
 637   Float_t minPhi = fConeMinPhi;
 638
 639   if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
 640   if (minPhi < 0) minPhi = 0;
 641
 642   Float_t dLJ = 0;
 643   Int_t repeats = 0;
 644   Bool_t reject = kTRUE;
 645   do {
 646     eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
 647     phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
 648     dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
 649
 650     repeats++;
 651   } while (dLJ < fMinRC2LJ && repeats < 999 && reject);
 652
 653   if (repeats == 999) {
 654     AliWarning(Form("%s: Could not get random cone!", GetName()));
 655     return;
 656   }
 657
 658   GetCone(lvRC,pt,eta,phi,tracks,clusters);
 659
 660
 661 }
 662
 663 //________________________________________________________________________
 664 void AliAnalysisTaskEmcalJetMassBkg::GetCone(TLorentzVector& lvRC,Float_t &pt, Float_t eta, Float_t phi, AliParticleContainer* tracks, AliClusterContainer* clusters) const
 665 {
 666
 667   pt = 0.;
 668   lvRC.SetPxPyPzE(0.,0.,0.,0.);
 669
 670   if (clusters) {
 671     AliVCluster* cluster = clusters->GetNextAcceptCluster(0);
 672     while (cluster) {
 673       TLorentzVector nPart;
 674       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
 675
 676       Float_t cluseta = nPart.Eta();
 677       Float_t clusphi = nPart.Phi();
 678
 679       if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi + 2 * TMath::Pi()))
 680         clusphi += 2 * TMath::Pi();
 681       if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi - 2 * TMath::Pi()))
 682         clusphi -= 2 * TMath::Pi();
 683
 684       Float_t d = TMath::Sqrt((cluseta - eta) * (cluseta - eta) + (clusphi - phi) * (clusphi - phi));
 685       if (d <= fConeRadius) {
 686         pt += nPart.Pt();
 687         TLorentzVector lvcl(nPart.Px(),nPart.Py(),nPart.Pz(),nPart.E());
 688         lvRC += lvcl;
 689       }
 690
 691       cluster = clusters->GetNextAcceptCluster();
 692     }
 693   }
 694
 695   if (tracks) {
 696     AliVParticle* track = tracks->GetNextAcceptParticle(0);
 697     while(track) {
 698       Float_t tracketa = track->Eta();
 699       Float_t trackphi = track->Phi();
 700
 701       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
 702         trackphi += 2 * TMath::Pi();
 703       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
 704         trackphi -= 2 * TMath::Pi();
 705
 706       Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
 707       if (d <= fConeRadius) {
 708         pt += track->Pt();
 709         TLorentzVector lvtr(track->Px(),track->Py(),track->Pz(),track->E());
 710         lvRC += lvtr;
 711       }
 712
 713       track = tracks->GetNextAcceptParticle();
 714     }
 715   }
 716
 717 }
 718
 719 //________________________________________________________________________
 720 void AliAnalysisTaskEmcalJetMassBkg::GetPerpCone(TLorentzVector& lvRC,Float_t &pt, Float_t &eta, Float_t &phi, AliParticleContainer* tracks, AliClusterContainer* clusters, AliEmcalJet *jet) const
 721 {
 722   // Get rigid cone.
 723   lvRC.SetPxPyPzE(0.,0.,0.,0.);
 724
 725   eta = -999;
 726   phi = -999;
 727   pt = 0;
 728
 729   if (!tracks && !clusters)
 730     return;
 731
 732   if(!jet)
 733     return;
 734
 735   Float_t LJeta = jet->Eta();
 736   Float_t LJphi = jet->Phi();
 737
 738   eta = LJeta;
 739   phi = LJphi + 0.5*TMath::Pi();
 740   if(phi>TMath::TwoPi()) phi-=TMath::TwoPi();
 741
 742   GetCone(lvRC,pt,eta,phi,tracks,clusters);
 743 }
 744
 745 //________________________________________________________________________
 746 void AliAnalysisTaskEmcalJetMassBkg::SetConeEtaPhiEMCAL()
 747 {
 748   // Set default cuts for full cones
 749
 750   SetConeEtaLimits(-0.7+fConeRadius,0.7-fConeRadius);
 751   SetConePhiLimits(1.405+fConeRadius,3.135-fConeRadius);
 752 }
 753
 754 //________________________________________________________________________
 755 void AliAnalysisTaskEmcalJetMassBkg::SetConeEtaPhiTPC()
 756 {
 757   // Set default cuts for charged cones
 758
 759   SetConeEtaLimits(-0.9+fConeRadius, 0.9-fConeRadius);
 760   SetConePhiLimits(-10, 10);
 761 }
 762
 763 //________________________________________________________________________
 764 void AliAnalysisTaskEmcalJetMassBkg::ExecOnce() {
 765
 766   AliAnalysisTaskEmcalJet::ExecOnce();
 767
 768   if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
 769   if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
 770
 771   if (fRCperEvent < 0) {
 772     Double_t area = (fConeMaxEta - fConeMinEta) * (fConeMaxPhi - fConeMinPhi);
 773     Double_t rcArea = TMath::Pi() * fConeRadius * fConeRadius;
 774     fRCperEvent = TMath::FloorNint(area / rcArea - 0.5);
 775     if (fRCperEvent == 0)
 776       fRCperEvent = 1;
 777   }
 778
 779   if (fMinRC2LJ < 0)
 780     fMinRC2LJ = fConeRadius * 1.5;
 781
 782   const Float_t maxDist = TMath::Max(fConeMaxPhi - fConeMinPhi, fConeMaxEta - fConeMinEta) / 2;
 783   if (fMinRC2LJ > maxDist) {
 784     AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
 785                     "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
 786     fMinRC2LJ = maxDist;
 787   }
 788
 789 }
 790
 791 //________________________________________________________________________
 792 Bool_t AliAnalysisTaskEmcalJetMassBkg::RetrieveEventObjects() {
 793   //
 794   // retrieve event objects
 795   //
 796
 797   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 798     return kFALSE;
 799
 800   return kTRUE;
 801
 802 }
 803
 804 //_______________________________________________________________________
 805 void AliAnalysisTaskEmcalJetMassBkg::Terminate(Option_t *)
 806 {
 807   // Called once at the end of the analysis.
 808 }
 809