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   TLorentzVector lvRCS(0.,0.,0.,0.);
 466   Float_t RCpt = 0;
 467   Float_t RCeta = 0;
 468   Float_t RCphi = 0;
 469
 470   static Double_t massvecRC[999];
 471   static Double_t massPerAreavecRC[999];
 472
 473   static Double_t massvecRCExLJ[999];
 474   static Double_t massPerAreavecRCExLJ[999];
 475
 476   Int_t nRCAcc = 0;
 477   Int_t nRCExLJAcc = 0;
 478
 479   for (Int_t i = 0; i < fRCperEvent; i++) {
 480     // Simple random cones
 481     lvRC.SetPxPyPzE(0.,0.,0.,0.);
 482     RCpt = 0;
 483     RCeta = 0;
 484     RCphi = 0;
 485     GetRandomCone(lvRC,RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, 0);
 486     if (RCpt > 0) {
 487       lvRCS = GetSubtractedVector(lvRC.Pt(),lvRC.Eta(),lvRC.Phi(),lvRC.E());
 488       fh2PtVsMassRC[fCentBin]->Fill(RCpt-rho*rcArea,lvRCS.M());
 489       fh2PtVsERC[fCentBin]->Fill(RCpt-rho*rcArea,lvRCS.E());
 490       fpPtVsMassRC[fCentBin]->Fill(RCpt-rho*rcArea,lvRCS.M());
 491       fh2EtaVsMassRC[fCentBin]->Fill(lvRCS.Eta(),lvRCS.M());
 492       fh2CentVsMassRC->Fill(fCent,lvRCS.M());
 493       fh2MultVsMassRC->Fill(trackMult,lvRCS.M());
 494
 495       massvecRC[nRCAcc] = lvRCS.M();
 496       massPerAreavecRC[nRCAcc] = massvecRC[nRCAcc]/rcArea;
 497       ++nRCAcc;
 498     }
 499
 500     if (fJetsCont && jet) {
 501       // Random cones far away from leading jet(s)
 502       lvRC.SetPxPyPzE(0.,0.,0.,0.);
 503       RCpt = 0;
 504       RCeta = 0;
 505       RCphi = 0;
 506       GetRandomCone(lvRC,RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet);
 507       if (RCpt > 0 && jet) {
 508         lvRCS = GetSubtractedVector(lvRC.Pt(),lvRC.Eta(),lvRC.Phi(),lvRC.E());
 509         Float_t dphi = RCphi - jet->Phi();
 510         if (dphi > 1.5*TMath::Pi()) dphi -= TMath::Pi() * 2;
 511         if (dphi < -0.5*TMath::Pi()) dphi += TMath::Pi() * 2;
 512         fh2PtVsMassRCExLJDPhi[fCentBin]->Fill(RCpt-rho*rcArea,lvRCS.M(),dphi);
 513         fh2PtVsERCExLJDPhi[fCentBin]->Fill(RCpt-rho*rcArea,lvRCS.E(),dphi);
 514         fpPtVsMassRCExLJ[fCentBin]->Fill(RCpt-rho*rcArea,lvRCS.M());
 515         fh2EtaVsMassRCExLJ[fCentBin]->Fill(lvRCS.Eta(),lvRCS.M());
 516         fh2CentVsMassRCExLJ->Fill(fCent,lvRCS.M());
 517         fh2MultVsMassRCExLJ->Fill(trackMult,lvRCS.M());
 518
 519         massvecRCExLJ[nRCExLJAcc] = lvRCS.M();
 520         massPerAreavecRCExLJ[nRCExLJAcc] = lvRCS.M()/rcArea;
 521         ++nRCExLJAcc;
 522       }
 523     }
 524   }//RC loop
 525
 526   Double_t medianRC, medianRCExLJ = 0.;
 527   medianRC = TMath::Median(nRCAcc,massvecRC);
 528   medianRCExLJ = TMath::Median(nRCExLJAcc,massvecRCExLJ);
 529
 530   fh2CentVsMedianMassRC->Fill(fCent,medianRC);
 531   fh2CentVsMedianMassRCExLJ->Fill(fCent,medianRCExLJ);
 532
 533   fh2MultVsMedianMassRC->Fill(trackMult,medianRC);
 534   fh2MultVsMedianMassRCExLJ->Fill(trackMult,medianRCExLJ);
 535
 536   Double_t meanRC = 0.; Double_t meanRCExLJ = 0.;
 537   if(nRCAcc>0)   meanRC = TMath::Mean(nRCAcc,massvecRC);
 538   if(nRCExLJAcc) meanRCExLJ = TMath::Mean(nRCExLJAcc,massvecRCExLJ);
 539
 540   fh2CentVsMeanMassRC->Fill(fCent,meanRC);
 541   fh2CentVsMeanMassRCExLJ->Fill(fCent,meanRCExLJ);
 542
 543   fh2MultVsMeanMassRC->Fill(trackMult,meanRC);
 544   fh2MultVsMeanMassRCExLJ->Fill(trackMult,meanRCExLJ);
 545
 546   Double_t medianPerAreaRC, medianPerAreaRCExLJ = 0.;
 547   medianPerAreaRC = TMath::Median(nRCAcc,massPerAreavecRC);
 548   medianPerAreaRCExLJ = TMath::Median(nRCExLJAcc,massPerAreavecRCExLJ);
 549
 550   fh2CentVsMedianMassPerAreaRC->Fill(fCent,medianPerAreaRC);
 551   fh2CentVsMedianMassPerAreaRCExLJ->Fill(fCent,medianPerAreaRCExLJ);
 552
 553   fh2MultVsMedianMassPerAreaRC->Fill(trackMult,medianPerAreaRC);
 554   fh2MultVsMedianMassPerAreaRCExLJ->Fill(trackMult,medianPerAreaRCExLJ);
 555
 556
 557   if(fJetsCont && jet) {
 558     //cone perpendicular to leading jet
 559     TLorentzVector lvPC(0.,0.,0.,0.);
 560     TLorentzVector lvPCS(0.,0.,0.,0.);
 561     Float_t PCpt = 0;
 562     Float_t PCeta = 0;
 563     Float_t PCphi = 0;
 564     if(jet) {
 565       GetPerpCone(lvPC,PCpt, PCeta, PCphi, fTracksCont, fCaloClustersCont, jet);
 566       if(PCpt>0.) {
 567         lvPCS = GetSubtractedVector(lvPC.Pt(),lvPC.Eta(),lvPC.Phi(),lvPC.E());
 568         fh2PtVsMassPerpConeLJ[fCentBin]->Fill(PCpt-rho*rcArea,lvPCS.M());
 569         fh2PtVsEPerpConeLJ[fCentBin]->Fill(PCpt-rho*rcArea,lvPCS.E());
 570         fpPtVsMassPerpConeLJ[fCentBin]->Fill(PCpt-rho*rcArea,lvPCS.M());
 571         fh2EtaVsMassPerpConeLJ[fCentBin]->Fill(lvPCS.Eta(),lvPCS.M());
 572         fh2CentVsMassPerpConeLJ->Fill(fCent,lvPCS.M());
 573         fh2MultVsMassPerpConeLJ->Fill(trackMult,lvPCS.M());
 574       }
 575     }
 576     //cone perpendicular to all tagged jets
 577     for(int i = 0; i < fJetsCont->GetNJets();++i) {
 578       jet = static_cast<AliEmcalJet*>(fJetsCont->GetAcceptJet(i));
 579       if(!jet) continue;
 580
 581       if(jet->GetTagStatus()<1)
 582         continue;
 583
 584       lvPC.SetPxPyPzE(0.,0.,0.,0.);
 585       PCpt = 0;
 586       PCeta = 0;
 587       PCphi = 0;
 588       GetPerpCone(lvPC,PCpt, PCeta, PCphi, fTracksCont, fCaloClustersCont, jet);
 589       if(PCpt>0.) {
 590         lvPCS = GetSubtractedVector(lvPC.Pt(),lvPC.Eta(),lvPC.Phi(),lvPC.E());
 591         fh2PtVsMassPerpConeTJ[fCentBin]->Fill(PCpt-rho*rcArea,lvPCS.M());
 592         fh2PtVsEPerpConeTJ[fCentBin]->Fill(PCpt-rho*rcArea,lvPCS.E());
 593         fpPtVsMassPerpConeTJ[fCentBin]->Fill(PCpt-rho*rcArea,lvPCS.M());
 594         fh2EtaVsMassPerpConeTJ[fCentBin]->Fill(lvPCS.Eta(),lvPCS.M());
 595         fh2CentVsMassPerpConeTJ->Fill(fCent,lvPCS.M());
 596         fh2MultVsMassPerpConeTJ->Fill(trackMult,lvPCS.M());
 597       }
 598     }//jet loop
 599   }
 600
 601   return kTRUE;
 602
 603 }
 604
 605 //________________________________________________________________________
 606 void AliAnalysisTaskEmcalJetMassBkg::GetRandomCone(TLorentzVector& lvRC,Float_t &pt, Float_t &eta, Float_t &phi,
 607                                                    AliParticleContainer* tracks, AliClusterContainer* clusters,
 608                                                    AliEmcalJet *jet) const
 609 {
 610   // Get rigid cone.
 611   lvRC.SetPxPyPzE(0.,0.,0.,0.);
 612
 613   eta = -999;
 614   phi = -999;
 615   pt = 0;
 616
 617   if (!tracks && !clusters)
 618     return;
 619
 620   Float_t LJeta = 999;
 621   Float_t LJphi = 999;
 622
 623   if (jet) {
 624     LJeta = jet->Eta();
 625     LJphi = jet->Phi();
 626   }
 627
 628   Float_t maxEta = fConeMaxEta;
 629   Float_t minEta = fConeMinEta;
 630   Float_t maxPhi = fConeMaxPhi;
 631   Float_t minPhi = fConeMinPhi;
 632
 633   if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
 634   if (minPhi < 0) minPhi = 0;
 635
 636   Float_t dLJ = 0;
 637   Int_t repeats = 0;
 638   Bool_t reject = kTRUE;
 639   do {
 640     eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
 641     phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
 642     dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
 643
 644     repeats++;
 645   } while (dLJ < fMinRC2LJ && repeats < 999 && reject);
 646
 647   if (repeats == 999) {
 648     AliWarning(Form("%s: Could not get random cone!", GetName()));
 649     return;
 650   }
 651
 652   GetCone(lvRC,pt,eta,phi,tracks,clusters);
 653
 654
 655 }
 656
 657 //________________________________________________________________________
 658 void AliAnalysisTaskEmcalJetMassBkg::GetCone(TLorentzVector& lvRC,Float_t &pt, Float_t eta, Float_t phi, AliParticleContainer* tracks, AliClusterContainer* clusters) const
 659 {
 660
 661   pt = 0.;
 662   lvRC.SetPxPyPzE(0.,0.,0.,0.);
 663
 664   if (clusters) {
 665     AliVCluster* cluster = clusters->GetNextAcceptCluster(0);
 666     while (cluster) {
 667       TLorentzVector nPart;
 668       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
 669
 670       Float_t cluseta = nPart.Eta();
 671       Float_t clusphi = nPart.Phi();
 672
 673       if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi + 2 * TMath::Pi()))
 674         clusphi += 2 * TMath::Pi();
 675       if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi - 2 * TMath::Pi()))
 676         clusphi -= 2 * TMath::Pi();
 677
 678       Float_t d = TMath::Sqrt((cluseta - eta) * (cluseta - eta) + (clusphi - phi) * (clusphi - phi));
 679       if (d <= fConeRadius) {
 680         pt += nPart.Pt();
 681         TLorentzVector lvcl(nPart.Px(),nPart.Py(),nPart.Pz(),nPart.E());
 682         lvRC += lvcl;
 683       }
 684
 685       cluster = clusters->GetNextAcceptCluster();
 686     }
 687   }
 688
 689   if (tracks) {
 690     AliVParticle* track = tracks->GetNextAcceptParticle(0);
 691     while(track) {
 692       Float_t tracketa = track->Eta();
 693       Float_t trackphi = track->Phi();
 694
 695       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
 696         trackphi += 2 * TMath::Pi();
 697       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
 698         trackphi -= 2 * TMath::Pi();
 699
 700       Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
 701       if (d <= fConeRadius) {
 702         pt += track->Pt();
 703         TLorentzVector lvtr(track->Px(),track->Py(),track->Pz(),track->E());
 704         lvRC += lvtr;
 705       }
 706
 707       track = tracks->GetNextAcceptParticle();
 708     }
 709   }
 710
 711 }
 712
 713 //________________________________________________________________________
 714 void AliAnalysisTaskEmcalJetMassBkg::GetPerpCone(TLorentzVector& lvRC,Float_t &pt, Float_t &eta, Float_t &phi, AliParticleContainer* tracks, AliClusterContainer* clusters, AliEmcalJet *jet) const
 715 {
 716   // Get rigid cone.
 717   lvRC.SetPxPyPzE(0.,0.,0.,0.);
 718
 719   eta = -999;
 720   phi = -999;
 721   pt = 0;
 722
 723   if (!tracks && !clusters)
 724     return;
 725
 726   if(!jet)
 727     return;
 728
 729   Float_t LJeta = jet->Eta();
 730   Float_t LJphi = jet->Phi();
 731
 732   eta = LJeta;
 733   phi = LJphi + 0.5*TMath::Pi();
 734   if(phi>TMath::TwoPi()) phi-=TMath::TwoPi();
 735
 736   GetCone(lvRC,pt,eta,phi,tracks,clusters);
 737 }
 738
 739 //________________________________________________________________________
 740 void AliAnalysisTaskEmcalJetMassBkg::SetConeEtaPhiEMCAL()
 741 {
 742   // Set default cuts for full cones
 743
 744   SetConeEtaLimits(-0.7+fConeRadius,0.7-fConeRadius);
 745   SetConePhiLimits(1.405+fConeRadius,3.135-fConeRadius);
 746 }
 747
 748 //________________________________________________________________________
 749 void AliAnalysisTaskEmcalJetMassBkg::SetConeEtaPhiTPC()
 750 {
 751   // Set default cuts for charged cones
 752
 753   SetConeEtaLimits(-0.9+fConeRadius, 0.9-fConeRadius);
 754   SetConePhiLimits(-10, 10);
 755 }
 756
 757 //________________________________________________________________________
 758 void AliAnalysisTaskEmcalJetMassBkg::ExecOnce() {
 759
 760   AliAnalysisTaskEmcalJet::ExecOnce();
 761
 762   if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
 763   if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
 764
 765   if (fRCperEvent < 0) {
 766     Double_t area = (fConeMaxEta - fConeMinEta) * (fConeMaxPhi - fConeMinPhi);
 767     Double_t rcArea = TMath::Pi() * fConeRadius * fConeRadius;
 768     fRCperEvent = TMath::FloorNint(area / rcArea - 0.5);
 769     if (fRCperEvent == 0)
 770       fRCperEvent = 1;
 771   }
 772
 773   if (fMinRC2LJ < 0)
 774     fMinRC2LJ = fConeRadius * 1.5;
 775
 776   const Float_t maxDist = TMath::Max(fConeMaxPhi - fConeMinPhi, fConeMaxEta - fConeMinEta) / 2;
 777   if (fMinRC2LJ > maxDist) {
 778     AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
 779                     "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
 780     fMinRC2LJ = maxDist;
 781   }
 782
 783 }
 784
 785   //________________________________________________________________________
 786 TLorentzVector AliAnalysisTaskEmcalJetMassBkg::GetBkgVector(Double_t eta, Double_t phi, AliJetContainer *cont) {
 787   //get background vector
 788
 789   Double_t rho  = cont->GetRhoVal();
 790   Double_t rhom = cont->GetRhoMassVal();
 791   TLorentzVector vpB(0.,0.,0.,0.);
 792   Double_t aRC = TMath::Pi()*fConeRadius*fConeRadius;
 793   vpB.SetPxPyPzE(rho*TMath::Cos(phi)*aRC,rho*TMath::Sin(phi)*aRC,(rho+rhom)*TMath::SinH(eta)*aRC,(rho+rhom)*TMath::CosH(eta)*aRC);
 794   return vpB;
 795 }
 796
 797 //________________________________________________________________________
 798 TLorentzVector AliAnalysisTaskEmcalJetMassBkg::GetSubtractedVector(Double_t pt, Double_t eta, Double_t phi, Double_t e) {
 799   //get subtracted vector
 800   TLorentzVector vpS;
 801   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 802   TLorentzVector vpBkg = GetBkgVector(eta,phi,jetCont);
 803   vpS.SetPxPyPzE(pt*TMath::Cos(phi)-vpBkg.Px(),pt*TMath::Sin(phi)-vpBkg.Py(),pt*TMath::SinH(eta)-vpBkg.Pz(),e-vpBkg.E());
 804   return vpS;
 805 }
 806
 807 //________________________________________________________________________
 808 Bool_t AliAnalysisTaskEmcalJetMassBkg::RetrieveEventObjects() {
 809   //
 810   // retrieve event objects
 811   //
 812
 813   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 814     return kFALSE;
 815
 816   return kTRUE;
 817
 818 }
 819
 820 //_______________________________________________________________________
 821 void AliAnalysisTaskEmcalJetMassBkg::Terminate(Option_t *)
 822 {
 823   // Called once at the end of the analysis.
 824 }
 825