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