PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalHJetMass.cxx

   1 //
   2 // Jet mass analysis task for jets recoiling from high pT hadron
   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 <TArrayI.h>
  19 #include <TArrayF.h>
  20 #include <TRandom3.h>
  21
  22 #include "AliVCluster.h"
  23 #include "AliVTrack.h"
  24 #include "AliEmcalJet.h"
  25 #include "AliRhoParameter.h"
  26 #include "AliLog.h"
  27 #include "AliEmcalParticle.h"
  28 #include "AliAnalysisManager.h"
  29 #include "AliJetContainer.h"
  30 #include "AliParticleContainer.h"
  31
  32 #include "AliAODEvent.h"
  33
  34 #include "AliAnalysisTaskEmcalHJetMass.h"
  35
  36 ClassImp(EmcalHJetMassAnalysis::AliAnalysisTaskEmcalHJetMass)
  37
  38 namespace EmcalHJetMassAnalysis {
  39
  40   //________________________________________________________________________
  41   AliAnalysisTaskEmcalHJetMass::AliAnalysisTaskEmcalHJetMass() :
  42     AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalHJetMass", kTRUE),
  43     fDoHJetAna(kTRUE),
  44     fDoNSHJetAna(kFALSE),
  45     fContainerBase(0),
  46     fContainerUnsub(1),
  47     fMinFractionShared(0),
  48     fUseUnsubJet(0),
  49     fJetMassType(kRaw),
  50     fDPhiHJetMax(0.6),
  51     fTriggerTrackType(kInclusive),
  52     fPtTTMin(0),
  53     fPtTTMax(0),
  54     fRandom(0),
  55     fEmbConstSel(0),
  56     fMarkMCLabel(-1),
  57     fh1PtHadron(0),
  58     fh1PtHadronMatch(0),
  59     fh3PtHPtJDPhi(0),
  60     fh3PtJet1VsMassVsHPtAllSel(0),
  61     fh3PtJet1VsMassVsHPtAllSelMatch(0),
  62     fh3PtJet1VsMassVsHPtTagged(0),
  63     fh3PtJet1VsMassVsHPtTaggedMatch(0),
  64     fh3PtJet1VsRatVsHPtAllSel(0),
  65     fh3PtJet1VsRatVsHPtAllSelMatch(0),
  66     fh3PtJet1VsRatVsHPtTagged(0),
  67     fh3PtJet1VsRatVsHPtTaggedMatch(0),
  68     fhnAllSel(0),
  69     fhnAllSelMatch(0),
  70     fhnTagged(0),
  71     fhnTaggedMatch(0)
  72   {
  73     // Default constructor.
  74
  75     fh1PtHadron                       = new TH1F*[fNcentBins];
  76     fh1PtHadronMatch                  = new TH1F*[fNcentBins];
  77     fh3PtHPtJDPhi                     = new TH3F*[fNcentBins];
  78     fh3PtJet1VsMassVsHPtAllSel        = new TH3F*[fNcentBins];
  79     fh3PtJet1VsMassVsHPtAllSelMatch   = new TH3F*[fNcentBins];
  80     fh3PtJet1VsMassVsHPtTagged        = new TH3F*[fNcentBins];
  81     fh3PtJet1VsMassVsHPtTaggedMatch   = new TH3F*[fNcentBins];
  82     fh3PtJet1VsRatVsHPtAllSel         = new TH3F*[fNcentBins];
  83     fh3PtJet1VsRatVsHPtAllSelMatch    = new TH3F*[fNcentBins];
  84     fh3PtJet1VsRatVsHPtTagged         = new TH3F*[fNcentBins];
  85     fh3PtJet1VsRatVsHPtTaggedMatch    = new TH3F*[fNcentBins];
  86     fhnAllSel                         = new THnSparse*[fNcentBins];
  87     fhnAllSelMatch                    = new THnSparse*[fNcentBins];
  88     fhnTagged                         = new THnSparse*[fNcentBins];
  89     fhnTaggedMatch                    = new THnSparse*[fNcentBins];
  90
  91     for (Int_t i = 0; i < fNcentBins; i++) {
  92       fh1PtHadron[i]                       = 0;
  93       fh1PtHadronMatch[i]                  = 0;
  94       fh3PtHPtJDPhi[i]                     = 0;
  95       fh3PtJet1VsMassVsHPtAllSel[i]        = 0;
  96       fh3PtJet1VsMassVsHPtAllSelMatch[i]   = 0;
  97       fh3PtJet1VsMassVsHPtTagged[i]        = 0;
  98       fh3PtJet1VsMassVsHPtTaggedMatch[i]   = 0;
  99       fh3PtJet1VsRatVsHPtAllSel[i]         = 0;
 100       fh3PtJet1VsRatVsHPtAllSelMatch[i]    = 0;
 101       fh3PtJet1VsRatVsHPtTagged[i]         = 0;
 102       fh3PtJet1VsRatVsHPtTaggedMatch[i]    = 0;
 103       fhnAllSel[i]                         = 0;
 104       fhnAllSelMatch[i]                    = 0;
 105       fhnTagged[i]                         = 0;
 106       fhnTaggedMatch[i]                    = 0;
 107     }
 108
 109     fPtTTMin = new TArrayF();
 110     fPtTTMax = new TArrayF();
 111
 112     SetMakeGeneralHistograms(kTRUE);
 113   }
 114
 115   //________________________________________________________________________
 116   AliAnalysisTaskEmcalHJetMass::AliAnalysisTaskEmcalHJetMass(const char *name) :
 117     AliAnalysisTaskEmcalJet(name, kTRUE),
 118     fDoHJetAna(kTRUE),
 119     fDoNSHJetAna(kFALSE),
 120     fContainerBase(0),
 121     fContainerUnsub(1),
 122     fMinFractionShared(0),
 123     fUseUnsubJet(0),
 124     fJetMassType(kRaw),
 125     fDPhiHJetMax(0.6),
 126     fTriggerTrackType(kInclusive),
 127     fPtTTMin(0),
 128     fPtTTMax(0),
 129     fRandom(0),
 130     fEmbConstSel(0),
 131     fMarkMCLabel(-1),
 132     fh1PtHadron(0),
 133     fh1PtHadronMatch(0),
 134     fh3PtHPtJDPhi(0),
 135     fh3PtJet1VsMassVsHPtAllSel(0),
 136     fh3PtJet1VsMassVsHPtAllSelMatch(0),
 137     fh3PtJet1VsMassVsHPtTagged(0),
 138     fh3PtJet1VsMassVsHPtTaggedMatch(0),
 139     fh3PtJet1VsRatVsHPtAllSel(0),
 140     fh3PtJet1VsRatVsHPtAllSelMatch(0),
 141     fh3PtJet1VsRatVsHPtTagged(0),
 142     fh3PtJet1VsRatVsHPtTaggedMatch(0),
 143     fhnAllSel(0),
 144     fhnAllSelMatch(0),
 145     fhnTagged(0),
 146     fhnTaggedMatch(0)
 147   {
 148     // Standard constructor.
 149
 150     fh1PtHadron                       = new TH1F*[fNcentBins];
 151     fh1PtHadronMatch                  = new TH1F*[fNcentBins];
 152     fh3PtHPtJDPhi                     = new TH3F*[fNcentBins];
 153     fh3PtJet1VsMassVsHPtAllSel        = new TH3F*[fNcentBins];
 154     fh3PtJet1VsMassVsHPtAllSelMatch   = new TH3F*[fNcentBins];
 155     fh3PtJet1VsMassVsHPtTagged        = new TH3F*[fNcentBins];
 156     fh3PtJet1VsMassVsHPtTaggedMatch   = new TH3F*[fNcentBins];
 157     fh3PtJet1VsRatVsHPtAllSel         = new TH3F*[fNcentBins];
 158     fh3PtJet1VsRatVsHPtAllSelMatch    = new TH3F*[fNcentBins];
 159     fh3PtJet1VsRatVsHPtTagged         = new TH3F*[fNcentBins];
 160     fh3PtJet1VsRatVsHPtTaggedMatch    = new TH3F*[fNcentBins];
 161     fhnAllSel                         = new THnSparse*[fNcentBins];
 162     fhnAllSelMatch                    = new THnSparse*[fNcentBins];
 163     fhnTagged                         = new THnSparse*[fNcentBins];
 164     fhnTaggedMatch                    = new THnSparse*[fNcentBins];
 165
 166     for (Int_t i = 0; i < fNcentBins; i++) {
 167       fh1PtHadron[i]                       = 0;
 168       fh1PtHadronMatch[i]                  = 0;
 169       fh3PtHPtJDPhi[i]                     = 0;
 170       fh3PtJet1VsMassVsHPtAllSel[i]        = 0;
 171       fh3PtJet1VsMassVsHPtAllSelMatch[i]   = 0;
 172       fh3PtJet1VsMassVsHPtTagged[i]        = 0;
 173       fh3PtJet1VsMassVsHPtTaggedMatch[i]   = 0;
 174       fh3PtJet1VsRatVsHPtAllSel[i]         = 0;
 175       fh3PtJet1VsRatVsHPtAllSelMatch[i]    = 0;
 176       fh3PtJet1VsRatVsHPtTagged[i]         = 0;
 177       fh3PtJet1VsRatVsHPtTaggedMatch[i]    = 0;
 178       fhnAllSel[i]                         = 0;
 179       fhnAllSelMatch[i]                    = 0;
 180       fhnTagged[i]                         = 0;
 181       fhnTaggedMatch[i]                    = 0;
 182     }
 183
 184     fPtTTMin = new TArrayF();
 185     fPtTTMax = new TArrayF();
 186
 187     SetMakeGeneralHistograms(kTRUE);
 188   }
 189
 190   //________________________________________________________________________
 191   AliAnalysisTaskEmcalHJetMass::~AliAnalysisTaskEmcalHJetMass()
 192   {
 193     // Destructor.
 194
 195     if(fRandom)      delete fRandom;
 196     if(fPtTTMin)     delete fPtTTMin;
 197     if(fPtTTMax)     delete fPtTTMax;
 198   }
 199
 200   //________________________________________________________________________
 201   void AliAnalysisTaskEmcalHJetMass::UserCreateOutputObjects()
 202   {
 203     // Create user output.
 204
 205     AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 206
 207     Bool_t oldStatus = TH1::AddDirectoryStatus();
 208     TH1::AddDirectory(kFALSE);
 209
 210     const Int_t nBinsPt  = 200;
 211     const Double_t minPt = -50.;
 212     const Double_t maxPt = 150.;
 213
 214     const Int_t nBinsM  = 90;
 215     const Double_t minM = -10.;
 216     const Double_t maxM = 80.;
 217
 218     const Int_t nBinsR  = 100;
 219     const Double_t minR = -0.2;
 220     const Double_t maxR = 0.8;
 221
 222     const Int_t nBinsPtH  = 100;
 223     const Double_t minPtH = 0.;
 224     const Double_t maxPtH = 100.;
 225
 226     const Int_t nBinsPhi  = 18*4;
 227     const Double_t minPhi = -0.5*TMath::Pi();
 228     const Double_t maxPhi = 1.5*TMath::Pi();
 229
 230     const Int_t nBinsSparse0 = 4; //PtJetAS,MJetAS,PtHNS,MJetNS
 231     const Int_t nBins0[nBinsSparse0] = {nBinsPt,nBinsM,nBinsPtH,nBinsM};
 232     const Double_t xmin0[nBinsSparse0]  = { minPt, minM, minPtH, minM};
 233     const Double_t xmax0[nBinsSparse0]  = { maxPt, maxM, maxPtH, maxM};
 234
 235     TString histName = "";
 236     TString histTitle = "";
 237     for (Int_t i = 0; i < fNcentBins; i++) {
 238       histName = TString::Format("fh1PtHadron_%d",i);
 239       histTitle = TString::Format("%s;#it{p}_{T,h}",histName.Data());
 240       fh1PtHadron[i] = new TH1F(histName.Data(),histTitle.Data(),200.,0.,200.);
 241       fOutput->Add(fh1PtHadron[i]);
 242
 243       histName = TString::Format("fh1PtHadronMatch_%d",i);
 244       histTitle = TString::Format("%s;#it{p}_{T,h}",histName.Data());
 245       fh1PtHadronMatch[i] = new TH1F(histName.Data(),histTitle.Data(),200.,0.,200.);
 246       fOutput->Add(fh1PtHadronMatch[i]);
 247
 248       histName = TString::Format("fh3PtHPtJDPhi_%d",i);
 249       histTitle = TString::Format("%s;#it{p}_{T,h};#it{p}_{T,jet};#Delta#varphi_{h,jet}",histName.Data());
 250       fh3PtHPtJDPhi[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPtH,minPtH,maxPtH,nBinsPt,minPt,maxPt,nBinsPhi,minPhi,maxPhi);
 251       fOutput->Add(fh3PtHPtJDPhi[i]);
 252
 253       if(fDoHJetAna) {
 254         histName = TString::Format("fh3PtJet1VsMassVsHPtAllSel_%d",i);
 255         histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
 256         fh3PtJet1VsMassVsHPtAllSel[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
 257         fOutput->Add(fh3PtJet1VsMassVsHPtAllSel[i]);
 258
 259         histName = TString::Format("fh3PtJet1VsMassVsHPtAllSelMatch_%d",i);
 260         histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
 261         fh3PtJet1VsMassVsHPtAllSelMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
 262         fOutput->Add(fh3PtJet1VsMassVsHPtAllSelMatch[i]);
 263
 264         histName = TString::Format("fh3PtJet1VsMassVsHPtTagged_%d",i);
 265         histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
 266         fh3PtJet1VsMassVsHPtTagged[i] =  new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
 267         fOutput->Add(fh3PtJet1VsMassVsHPtTagged[i]);
 268
 269         histName = TString::Format("fh3PtJet1VsMassVsHPtTaggedMatch_%d",i);
 270         histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
 271         fh3PtJet1VsMassVsHPtTaggedMatch[i] =  new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
 272         fOutput->Add(fh3PtJet1VsMassVsHPtTaggedMatch[i]);
 273
 274         //
 275         histName = TString::Format("fh3PtJet1VsRatVsHPtAllSel_%d",i);
 276         histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
 277         fh3PtJet1VsRatVsHPtAllSel[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
 278         fOutput->Add(fh3PtJet1VsRatVsHPtAllSel[i]);
 279
 280         histName = TString::Format("fh3PtJet1VsRatVsHPtAllSelMatch_%d",i);
 281         histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
 282         fh3PtJet1VsRatVsHPtAllSelMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
 283         fOutput->Add(fh3PtJet1VsRatVsHPtAllSelMatch[i]);
 284
 285         histName = TString::Format("fh3PtJet1VsRatVsHPtTagged_%d",i);
 286         histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
 287         fh3PtJet1VsRatVsHPtTagged[i] =  new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
 288         fOutput->Add(fh3PtJet1VsRatVsHPtTagged[i]);
 289
 290         histName = TString::Format("fh3PtJet1VsRatVsHPtTaggedMatch_%d",i);
 291         histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
 292         fh3PtJet1VsRatVsHPtTaggedMatch[i] =  new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
 293         fOutput->Add(fh3PtJet1VsRatVsHPtTaggedMatch[i]);
 294       }
 295
 296       if(fDoNSHJetAna) {
 297         histName = TString::Format("fhnAllSel_%d",i);
 298         histTitle = Form("%s;#it{p}_{T,jet}^{AS};#it{M}_{jet}^{AS};#it{p}_{T,h}^{NS};#it{M}_{jet}^{NS}",histName.Data());
 299         fhnAllSel[i] = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
 300         fOutput->Add(fhnAllSel[i]);
 301
 302         histName = TString::Format("fhnAllSelMatch_%d",i);
 303         histTitle = Form("%s;#it{p}_{T,jet}^{AS};#it{M}_{jet}^{AS};#it{p}_{T,h}^{NS};#it{M}_{jet}^{NS}",histName.Data());
 304         fhnAllSelMatch[i] = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
 305         fOutput->Add(fhnAllSelMatch[i]);
 306
 307         histName = TString::Format("fhnTagged_%d",i);
 308         histTitle = Form("%s;#it{p}_{T,jet}^{AS};#it{M}_{jet}^{AS};#it{p}_{T,h}^{NS};#it{M}_{jet}^{NS}",histName.Data());
 309         fhnTagged[i] = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
 310         fOutput->Add(fhnTagged[i]);
 311
 312         histName = TString::Format("fhnTaggedMatch_%d",i);
 313         histTitle = Form("%s;#it{p}_{T,jet}^{AS};#it{M}_{jet}^{AS};#it{p}_{T,h}^{NS};#it{M}_{jet}^{NS}",histName.Data());
 314         fhnTaggedMatch[i] = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
 315         fOutput->Add(fhnTaggedMatch[i]);
 316       }
 317     }
 318
 319     TH1::AddDirectory(oldStatus);
 320
 321     fRandom = new TRandom3(0);
 322
 323     PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 324   }
 325
 326   //________________________________________________________________________
 327   Bool_t AliAnalysisTaskEmcalHJetMass::Run()
 328   {
 329     // Run analysis code here, if needed. It will be executed before FillHistograms().
 330
 331     AliParticleContainer *pCont = GetParticleContainer(0);
 332     AliJetContainer      *jCont = GetJetContainer(fContainerBase);
 333     if(!pCont || !jCont) return kFALSE;
 334
 335     AliVParticle *vp = NULL;
 336
 337     if(fTriggerTrackType==kInclusive) {
 338       pCont->ResetCurrentID();
 339       while((vp = pCont->GetNextAcceptParticle())) {
 340         fh1PtHadron[fCentBin]->Fill(vp->Pt()); //all hadrons
 341         if(fMarkMCLabel>0 && TMath::Abs(vp->GetLabel()) >= fMarkMCLabel )
 342           fh1PtHadronMatch[fCentBin]->Fill(vp->Pt()); //hadrons matched to MC
 343         AliEmcalJet* jet = NULL;
 344         if(jCont) {
 345           jCont->ResetCurrentID();
 346           while((jet = jCont->GetNextAcceptJet())) {
 347             Double_t dphi = GetDeltaPhi(vp,jet)-TMath::Pi();
 348             fh3PtHPtJDPhi[fCentBin]->Fill(vp->Pt(),jet->Pt() - GetRhoVal(fContainerBase)*jet->Area(),dphi);
 349             if(dphi>fDPhiHJetMax) continue;
 350             FillHJetHistograms(vp,jet);
 351           }
 352         }
 353       }
 354     }
 355     else if(fTriggerTrackType==kSingleInclusive) {
 356      for(Int_t it = 0; it<fPtTTMin->GetSize(); it++) {
 357        vp = GetSingleInclusiveTT(pCont,fPtTTMin->At(it),fPtTTMax->At(it));
 358        if(!vp) continue;
 359        fh1PtHadron[fCentBin]->Fill(vp->Pt()); //all trigger tracks
 360         if(fMarkMCLabel>0 && TMath::Abs(vp->GetLabel()) >= fMarkMCLabel )
 361           fh1PtHadronMatch[fCentBin]->Fill(vp->Pt()); //hadrons matched to MC
 362        AliEmcalJet* jet = NULL;
 363        jCont->ResetCurrentID();
 364        while((jet = jCont->GetNextAcceptJet())) {
 365          Double_t dphi = GetDeltaPhi(vp,jet)-TMath::Pi();
 366          fh3PtHPtJDPhi[fCentBin]->Fill(vp->Pt(),jet->Pt() - GetRhoVal(fContainerBase)*jet->Area(),dphi);
 367          if(dphi>fDPhiHJetMax) continue;
 368          FillHJetHistograms(vp,jet);
 369        }
 370      }//trigger track types
 371     }
 372     return kTRUE;
 373   }
 374
 375   //________________________________________________________________________
 376   AliVParticle* AliAnalysisTaskEmcalHJetMass::GetSingleInclusiveTT(AliParticleContainer *pCont, Double_t ptmin, Double_t ptmax) const {
 377     AliVParticle *vp;
 378     TArrayI arr; arr.Set(pCont->GetNParticles());
 379     arr.Reset();
 380     Int_t counter = -1;
 381     pCont->ResetCurrentID();
 382     while((vp = pCont->GetNextAcceptParticle())) {
 383       if(vp->Pt()>=ptmin && vp->Pt()<ptmax ) {
 384         counter++;
 385         arr.SetAt(pCont->GetCurrentID(),counter);
 386       }
 387     }
 388     if(counter<0) return NULL;
 389     //select trigger track randomly
 390     fRandom->SetSeed(arr.At(0)); //random selection reproducible
 391     Double_t rnd = fRandom->Uniform() * counter;
 392     Int_t trigID = arr.At(TMath::FloorNint(rnd));
 393     vp = pCont->GetParticle(trigID);
 394     return vp;
 395   }
 396
 397   //________________________________________________________________________
 398   Bool_t AliAnalysisTaskEmcalHJetMass::FillHJetHistograms(const AliVParticle *vp, const AliEmcalJet *jet)
 399   {
 400     // Fill hadron-jet histograms.
 401     Double_t pt = vp->Pt();
 402     Double_t ptJet = jet->Pt() - GetRhoVal(fContainerBase)*jet->Area();
 403     Double_t mJet = GetJetMass(jet);
 404     Double_t rat = -1.;
 405     if(ptJet<0. || ptJet>0.) rat = mJet/ptJet;
 406
 407     Double_t var[4] = {ptJet,mJet,pt,-99.};
 408     if(fDoNSHJetAna) {
 409       AliEmcalJet *jetNS = FindNearSideJet(vp);
 410       if(jetNS) var[3] = GetJetMass(jetNS);
 411     }
 412
 413     if(fDoHJetAna) {
 414       fh3PtJet1VsMassVsHPtAllSel[fCentBin]->Fill(ptJet,mJet,pt);
 415       fh3PtJet1VsRatVsHPtAllSel[fCentBin]->Fill(ptJet,rat,pt);
 416     }
 417
 418     if(fDoNSHJetAna)
 419       fhnAllSel[fCentBin]->Fill(var);
 420
 421     Double_t fraction = 1.;
 422     if(fUseUnsubJet) {
 423       AliEmcalJet *jetUS = NULL;
 424       AliJetContainer *jetContUS = GetJetContainer(fContainerUnsub);
 425       Int_t ifound = 0;
 426       Int_t ilab = -1;
 427       for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
 428         jetUS = jetContUS->GetJet(i);
 429         if(jetUS->GetLabel()==jet->GetLabel()) {
 430           ifound++;
 431           if(ifound==1) ilab = i;
 432         }
 433       }
 434       if(ifound>1) AliDebug(2,Form("Found %d partners",ifound));
 435       if(ifound==0) jetUS = 0x0;
 436       else          jetUS = jetContUS->GetJet(ilab);
 437       fraction = jetContUS->GetFractionSharedPt(jetUS);
 438     } else {
 439       AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 440       fraction = jetCont->GetFractionSharedPt(jet);
 441     }
 442
 443     Bool_t mcMatch = kFALSE;
 444     if(fMarkMCLabel>0 && TMath::Abs(vp->GetLabel()) >= fMarkMCLabel ) mcMatch = kTRUE;
 445     if(fMinFractionShared>0. && fraction>fMinFractionShared) mcMatch = kTRUE;
 446     else mcMatch = kFALSE;
 447
 448     if(mcMatch) {
 449       if(fDoHJetAna) {
 450         fh3PtJet1VsMassVsHPtAllSelMatch[fCentBin]->Fill(ptJet,mJet,pt);
 451         fh3PtJet1VsRatVsHPtAllSelMatch[fCentBin]->Fill(ptJet,rat,pt);
 452       }
 453       if(fDoNSHJetAna)
 454         fhnAllSelMatch[fCentBin]->Fill(var);
 455     }
 456
 457     if(jet->GetTagStatus()<1 || !jet->GetTaggedJet())
 458       return kFALSE;
 459
 460     if(fDoHJetAna) {
 461       fh3PtJet1VsMassVsHPtTagged[fCentBin]->Fill(ptJet,mJet,pt);
 462       fh3PtJet1VsRatVsHPtTagged[fCentBin]->Fill(ptJet,rat,pt);
 463     }
 464     if(fDoNSHJetAna)
 465       fhnTagged[fCentBin]->Fill(var);
 466
 467     if(mcMatch) {
 468       if(fDoHJetAna) {
 469         fh3PtJet1VsMassVsHPtTaggedMatch[fCentBin]->Fill(ptJet,mJet,pt);
 470         fh3PtJet1VsRatVsHPtTaggedMatch[fCentBin]->Fill(ptJet,rat,pt);
 471       }
 472     if(fDoNSHJetAna)
 473       fhnTaggedMatch[fCentBin]->Fill(var);
 474     }
 475     return kTRUE;
 476   }
 477
 478   //________________________________________________________________________
 479   AliEmcalJet* AliAnalysisTaskEmcalHJetMass::FindNearSideJet(const AliVParticle *vp) {
 480     AliJetContainer      *jCont = GetJetContainer(fContainerBase);
 481     AliEmcalJet* jet = NULL;
 482     if(jCont) {
 483       jCont->ResetCurrentID();
 484       while((jet = jCont->GetNextAcceptJet())) {
 485         Int_t n = (Int_t)jet->GetNumberOfTracks();
 486         if (n < 1) continue;
 487         for (Int_t i = 0; i < n; i++) {
 488           AliVParticle *vp2 = static_cast<AliVParticle*>(jet->TrackAt(i, jCont->GetParticleContainer()->GetArray()));
 489           if(!vp2) continue;
 490           if(vp->Phi()==vp2->Phi()) return jet;
 491         }
 492       }
 493     }
 494     return jet;
 495   }
 496
 497   //________________________________________________________________________
 498   Double_t AliAnalysisTaskEmcalHJetMass::GetJetMass(const AliEmcalJet *jet) const {
 499     //calc subtracted jet mass
 500     if(fJetMassType==kRaw)
 501       return jet->M();
 502     else if(fJetMassType==kDeriv)
 503       return jet->GetSecondOrderSubtracted();
 504
 505     return 0;
 506   }
 507
 508   //________________________________________________________________________
 509   Double_t AliAnalysisTaskEmcalHJetMass::GetDeltaPhi(const AliVParticle *vp, const AliEmcalJet* jet) const {
 510     // Calculate azimuthal angle between particle and jet. range:[-0.5\pi,1.5\pi]
 511     return GetDeltaPhi(vp->Phi(),jet->Phi());
 512   }
 513
 514   //________________________________________________________________________
 515   Double_t AliAnalysisTaskEmcalHJetMass::GetDeltaPhi(Double_t phi1, Double_t phi2) const {
 516     // Calculate azimuthal angle between phi1 and phi2. range:[-0.5\pi,1.5\pi]
 517     Double_t dPhi = phi1-phi2;
 518     if(dPhi <-0.5*TMath::Pi())  dPhi += TMath::TwoPi();
 519     if(dPhi > 1.5*TMath::Pi())  dPhi -= TMath::TwoPi();
 520     return dPhi;
 521   }
 522
 523   //________________________________________________________________________
 524   Bool_t AliAnalysisTaskEmcalHJetMass::RetrieveEventObjects() {
 525     //
 526     // retrieve event objects
 527     //
 528     if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 529       return kFALSE;
 530
 531     return kTRUE;
 532   }
 533
 534   //________________________________________________________________________
 535   void AliAnalysisTaskEmcalHJetMass::AddTriggerTrackPtCuts(Float_t min, Float_t max) {
 536     if(!fPtTTMin) fPtTTMin = new TArrayF();
 537     if(!fPtTTMax) fPtTTMax = new TArrayF();
 538     Int_t newSize = fPtTTMin->GetSize()+1;
 539     fPtTTMin->Set(newSize);
 540     fPtTTMax->Set(newSize);
 541     fPtTTMin->AddAt(min,newSize-1);
 542     fPtTTMax->AddAt(max,newSize-1);
 543   }
 544
 545   //_______________________________________________________________________
 546   void AliAnalysisTaskEmcalHJetMass::Terminate(Option_t *)
 547   {
 548     // Called once at the end of the analysis.
 549   }
 550 }
 551