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