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     fContainerBase(0),
  44     fContainerUnsub(1),
  45     fMinFractionShared(0),
  46     fUseUnsubJet(0),
  47     fJetMassType(kRaw),
  48     fDPhiHJetMax(0.6),
  49     fTriggerTrackType(kInclusive),
  50     fPtTTMin(0),
  51     fPtTTMax(0),
  52     fRandom(0),
  53     fh1PtHadron(0),
  54     fh3PtHPtJDPhi(0),
  55     fh3PtJet1VsMassVsHPtAllSel(0),
  56     fh3PtJet1VsMassVsHPtTagged(0),
  57     fh3PtJet1VsMassVsHPtTaggedMatch(0),
  58     fh3PtJet1VsRatVsHPtAllSel(0),
  59     fh3PtJet1VsRatVsHPtTagged(0),
  60     fh3PtJet1VsRatVsHPtTaggedMatch(0)
  61   {
  62     // Default constructor.
  63
  64     fh1PtHadron                       = new TH1F*[fNcentBins];
  65     fh3PtHPtJDPhi                     = new TH3F*[fNcentBins];
  66     fh3PtJet1VsMassVsHPtAllSel        = new TH3F*[fNcentBins];
  67     fh3PtJet1VsMassVsHPtTagged        = new TH3F*[fNcentBins];
  68     fh3PtJet1VsMassVsHPtTaggedMatch   = new TH3F*[fNcentBins];
  69     fh3PtJet1VsRatVsHPtAllSel         = new TH3F*[fNcentBins];
  70     fh3PtJet1VsRatVsHPtTagged         = new TH3F*[fNcentBins];
  71     fh3PtJet1VsRatVsHPtTaggedMatch    = new TH3F*[fNcentBins];
  72
  73     for (Int_t i = 0; i < fNcentBins; i++) {
  74       fh1PtHadron[i]                       = 0;
  75       fh3PtHPtJDPhi[i]                     = 0;
  76       fh3PtJet1VsMassVsHPtAllSel[i]        = 0;
  77       fh3PtJet1VsMassVsHPtTagged[i]        = 0;
  78       fh3PtJet1VsMassVsHPtTaggedMatch[i]   = 0;
  79       fh3PtJet1VsRatVsHPtAllSel[i]         = 0;
  80       fh3PtJet1VsRatVsHPtTagged[i]         = 0;
  81       fh3PtJet1VsRatVsHPtTaggedMatch[i]    = 0;
  82     }
  83
  84     fPtTTMin = new TArrayF();
  85     fPtTTMax = new TArrayF();
  86
  87     SetMakeGeneralHistograms(kTRUE);
  88   }
  89
  90   //________________________________________________________________________
  91   AliAnalysisTaskEmcalHJetMass::AliAnalysisTaskEmcalHJetMass(const char *name) :
  92     AliAnalysisTaskEmcalJet(name, kTRUE),
  93     fContainerBase(0),
  94     fContainerUnsub(1),
  95     fMinFractionShared(0),
  96     fUseUnsubJet(0),
  97     fJetMassType(kRaw),
  98     fDPhiHJetMax(0.6),
  99     fTriggerTrackType(kInclusive),
 100     fPtTTMin(0),
 101     fPtTTMax(0),
 102     fRandom(0),
 103     fh1PtHadron(0),
 104     fh3PtHPtJDPhi(0),
 105     fh3PtJet1VsMassVsHPtAllSel(0),
 106     fh3PtJet1VsMassVsHPtTagged(0),
 107     fh3PtJet1VsMassVsHPtTaggedMatch(0),
 108     fh3PtJet1VsRatVsHPtAllSel(0),
 109     fh3PtJet1VsRatVsHPtTagged(0),
 110     fh3PtJet1VsRatVsHPtTaggedMatch(0)
 111   {
 112     // Standard constructor.
 113
 114     fh1PtHadron                       = new TH1F*[fNcentBins];
 115     fh3PtHPtJDPhi                     = new TH3F*[fNcentBins];
 116     fh3PtJet1VsMassVsHPtAllSel        = new TH3F*[fNcentBins];
 117     fh3PtJet1VsMassVsHPtTagged        = new TH3F*[fNcentBins];
 118     fh3PtJet1VsMassVsHPtTaggedMatch   = new TH3F*[fNcentBins];
 119     fh3PtJet1VsRatVsHPtAllSel         = new TH3F*[fNcentBins];
 120     fh3PtJet1VsRatVsHPtTagged         = new TH3F*[fNcentBins];
 121     fh3PtJet1VsRatVsHPtTaggedMatch    = new TH3F*[fNcentBins];
 122
 123     for (Int_t i = 0; i < fNcentBins; i++) {
 124       fh1PtHadron[i]                       = 0;
 125       fh3PtHPtJDPhi[i]                     = 0;
 126       fh3PtJet1VsMassVsHPtAllSel[i]        = 0;
 127       fh3PtJet1VsMassVsHPtTagged[i]        = 0;
 128       fh3PtJet1VsMassVsHPtTaggedMatch[i]   = 0;
 129       fh3PtJet1VsRatVsHPtAllSel[i]         = 0;
 130       fh3PtJet1VsRatVsHPtTagged[i]         = 0;
 131       fh3PtJet1VsRatVsHPtTaggedMatch[i]    = 0;
 132     }
 133
 134     fPtTTMin = new TArrayF();
 135     fPtTTMax = new TArrayF();
 136
 137     SetMakeGeneralHistograms(kTRUE);
 138   }
 139
 140   //________________________________________________________________________
 141   AliAnalysisTaskEmcalHJetMass::~AliAnalysisTaskEmcalHJetMass()
 142   {
 143     // Destructor.
 144
 145     if(fRandom)      delete fRandom;
 146     if(fPtTTMin)     delete fPtTTMin;
 147     if(fPtTTMax)     delete fPtTTMax;
 148   }
 149
 150   //________________________________________________________________________
 151   void AliAnalysisTaskEmcalHJetMass::UserCreateOutputObjects()
 152   {
 153     // Create user output.
 154
 155     AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 156
 157     Bool_t oldStatus = TH1::AddDirectoryStatus();
 158     TH1::AddDirectory(kFALSE);
 159
 160     const Int_t nBinsPt  = 200;
 161     const Double_t minPt = -50.;
 162     const Double_t maxPt = 150.;
 163
 164     const Int_t nBinsM  = 100;
 165     const Double_t minM = -20.;
 166     const Double_t maxM = 80.;
 167
 168     const Int_t nBinsR  = 100;
 169     const Double_t minR = -0.2;
 170     const Double_t maxR = 0.8;
 171
 172     const Int_t nBinsPtH  = 100;
 173     const Double_t minPtH = 0.;
 174     const Double_t maxPtH = 100.;
 175
 176     const Int_t nBinsPhi  = 18*4;
 177     const Double_t minPhi = -0.5*TMath::Pi();
 178     const Double_t maxPhi = 1.5*TMath::Pi();
 179
 180     TString histName = "";
 181     TString histTitle = "";
 182     for (Int_t i = 0; i < fNcentBins; i++) {
 183       histName = TString::Format("fh1PtHadron_%d",i);
 184       histTitle = TString::Format("%s;#it{p}_{T,h}",histName.Data());
 185       fh1PtHadron[i] = new TH1F(histName.Data(),histTitle.Data(),200.,0.,200.);
 186       fOutput->Add(fh1PtHadron[i]);
 187
 188       histName = TString::Format("fh3PtHPtJDPhi_%d",i);
 189       histTitle = TString::Format("%s;#it{p}_{T,h};#it{p}_{T,jet};#Delta#varphi_{h,jet}",histName.Data());
 190       fh3PtHPtJDPhi[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPtH,minPtH,maxPtH,nBinsPt,minPt,maxPt,nBinsPhi,minPhi,maxPhi);
 191       fOutput->Add(fh3PtHPtJDPhi[i]);
 192
 193       histName = TString::Format("fh3PtJet1VsMassVsHPtAllSel_%d",i);
 194       histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
 195       fh3PtJet1VsMassVsHPtAllSel[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
 196       fOutput->Add(fh3PtJet1VsMassVsHPtAllSel[i]);
 197
 198       histName = TString::Format("fh3PtJet1VsMassVsHPtTagged_%d",i);
 199       histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
 200       fh3PtJet1VsMassVsHPtTagged[i] =  new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
 201       fOutput->Add(fh3PtJet1VsMassVsHPtTagged[i]);
 202
 203       histName = TString::Format("fh3PtJet1VsMassVsHPtTaggedMatch_%d",i);
 204       histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,h}",histName.Data());
 205       fh3PtJet1VsMassVsHPtTaggedMatch[i] =  new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,nBinsPtH,minPtH,maxPtH);
 206       fOutput->Add(fh3PtJet1VsMassVsHPtTaggedMatch[i]);
 207
 208       //
 209       histName = TString::Format("fh3PtJet1VsRatVsHPtAllSel_%d",i);
 210       histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
 211       fh3PtJet1VsRatVsHPtAllSel[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
 212       fOutput->Add(fh3PtJet1VsRatVsHPtAllSel[i]);
 213
 214       histName = TString::Format("fh3PtJet1VsRatVsHPtTagged_%d",i);
 215       histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
 216       fh3PtJet1VsRatVsHPtTagged[i] =  new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
 217       fOutput->Add(fh3PtJet1VsRatVsHPtTagged[i]);
 218
 219       histName = TString::Format("fh3PtJet1VsRatVsHPtTaggedMatch_%d",i);
 220       histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}/#it{p}_{T,jet1};#it{p}_{T,h}",histName.Data());
 221       fh3PtJet1VsRatVsHPtTaggedMatch[i] =  new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,nBinsPtH,minPtH,maxPtH);
 222       fOutput->Add(fh3PtJet1VsRatVsHPtTaggedMatch[i]);
 223     }
 224
 225     TH1::AddDirectory(oldStatus);
 226
 227     fRandom = new TRandom3(0);
 228
 229     PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 230   }
 231
 232   //________________________________________________________________________
 233   Bool_t AliAnalysisTaskEmcalHJetMass::Run()
 234   {
 235     // Run analysis code here, if needed. It will be executed before FillHistograms().
 236
 237     AliParticleContainer *pCont = GetParticleContainer(0);
 238     AliJetContainer      *jCont = GetJetContainer(fContainerBase);
 239     if(!pCont || !jCont) return kFALSE;
 240
 241     AliVParticle *vp = NULL;
 242
 243     if(fTriggerTrackType==kInclusive) {
 244       pCont->ResetCurrentID();
 245       while((vp = pCont->GetNextAcceptParticle())) {
 246         fh1PtHadron[fCentBin]->Fill(vp->Pt()); //all hadrons
 247         AliEmcalJet* jet = NULL;
 248         if(jCont) {
 249           jCont->ResetCurrentID();
 250           while((jet = jCont->GetNextAcceptJet())) {
 251             Double_t dphi = GetDeltaPhi(vp,jet)-TMath::Pi();
 252             fh3PtHPtJDPhi[fCentBin]->Fill(vp->Pt(),jet->Pt() - GetRhoVal(fContainerBase)*jet->Area(),dphi);
 253             if(dphi>fDPhiHJetMax) continue;
 254             FillHJetHistograms(vp->Pt(),jet);
 255           }
 256         }
 257       }
 258     }
 259     else if(fTriggerTrackType==kSingleInclusive) {
 260       TArrayI arr; arr.Set(pCont->GetNParticles());
 261       for(Int_t it = 0; it<fPtTTMin->GetSize(); it++) {
 262         Int_t counter = -1;
 263         arr.Reset();
 264         pCont->ResetCurrentID();
 265         while((vp = pCont->GetNextAcceptParticle())) {
 266           if(vp->Pt()>fPtTTMin->At(it) && vp->Pt()<fPtTTMax->At(it) ) {
 267             counter++;
 268             arr.SetAt(pCont->GetCurrentID(),counter);
 269           }
 270         }
 271         if(counter<0) continue;
 272         //select trigger track randomly
 273         fRandom->SetSeed(arr.At(0)); //random selection reproducible
 274         Double_t rnd = fRandom->Uniform() * counter;
 275         Int_t trigID = arr.At(TMath::FloorNint(rnd));
 276         vp = pCont->GetParticle(trigID);
 277         fh1PtHadron[fCentBin]->Fill(vp->Pt()); //all hadrons
 278         AliEmcalJet* jet = NULL;
 279         jCont->ResetCurrentID();
 280         while((jet = jCont->GetNextAcceptJet())) {
 281           Double_t dphi = GetDeltaPhi(vp,jet)-TMath::Pi();
 282           fh3PtHPtJDPhi[fCentBin]->Fill(vp->Pt(),jet->Pt() - GetRhoVal(fContainerBase)*jet->Area(),dphi);
 283           if(dphi>fDPhiHJetMax) continue;
 284           FillHJetHistograms(vp->Pt(),jet);
 285         }
 286       }
 287     }
 288     return kTRUE;
 289   }
 290
 291   //________________________________________________________________________
 292   Bool_t AliAnalysisTaskEmcalHJetMass::FillHJetHistograms(Double_t pt, const AliEmcalJet *jet)
 293   {
 294     // Fill hadron-jet histograms.
 295     Double_t ptJet = jet->Pt() - GetRhoVal(fContainerBase)*jet->Area();
 296     Double_t mJet = GetJetMass(jet);
 297     Double_t rat = -1.;
 298     if(ptJet<0. || ptJet>0.) rat = mJet/ptJet;
 299
 300     fh3PtJet1VsMassVsHPtAllSel[fCentBin]->Fill(ptJet,mJet,pt);
 301     fh3PtJet1VsRatVsHPtAllSel[fCentBin]->Fill(ptJet,rat,pt);
 302
 303     if(jet->GetTagStatus()<1 || !jet->GetTaggedJet())
 304       return kFALSE;
 305     fh3PtJet1VsMassVsHPtTagged[fCentBin]->Fill(ptJet,mJet,pt);
 306     fh3PtJet1VsRatVsHPtTagged[fCentBin]->Fill(ptJet,rat,pt);
 307
 308     Double_t fraction = -1.;
 309     if(fUseUnsubJet) {
 310       AliEmcalJet *jetUS = NULL;
 311       AliJetContainer *jetContUS = GetJetContainer(fContainerUnsub);
 312       Int_t ifound = 0;
 313       Int_t ilab = -1;
 314       for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
 315         jetUS = jetContUS->GetJet(i);
 316         if(jetUS->GetLabel()==jet->GetLabel()) {
 317           ifound++;
 318           if(ifound==1) ilab = i;
 319         }
 320       }
 321       if(ifound>1) AliDebug(2,Form("Found %d partners",ifound));
 322       if(ifound==0) jetUS = 0x0;
 323       else          jetUS = jetContUS->GetJet(ilab);
 324       fraction = jetContUS->GetFractionSharedPt(jetUS);
 325     } else {
 326       AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 327       fraction = jetCont->GetFractionSharedPt(jet);
 328     }
 329
 330     if(fMinFractionShared>0. && fraction>fMinFractionShared) {
 331       fh3PtJet1VsMassVsHPtTaggedMatch[fCentBin]->Fill(ptJet,mJet,pt);
 332       fh3PtJet1VsRatVsHPtTaggedMatch[fCentBin]->Fill(ptJet,rat,pt);
 333     }
 334     return kTRUE;
 335   }
 336
 337   //________________________________________________________________________
 338   Double_t AliAnalysisTaskEmcalHJetMass::GetJetMass(const AliEmcalJet *jet) const {
 339     //calc subtracted jet mass
 340     if(fJetMassType==kRaw)
 341       return jet->M();
 342     else if(fJetMassType==kDeriv)
 343       return jet->GetSecondOrderSubtracted();
 344
 345     return 0;
 346   }
 347
 348   //________________________________________________________________________
 349   Double_t AliAnalysisTaskEmcalHJetMass::GetDeltaPhi(const AliVParticle *vp, const AliEmcalJet* jet) const {
 350     // Calculate azimuthal angle between particle and jet. range:[-0.5\pi,1.5\pi]
 351     return GetDeltaPhi(vp->Phi(),jet->Phi());
 352   }
 353
 354   //________________________________________________________________________
 355   Double_t AliAnalysisTaskEmcalHJetMass::GetDeltaPhi(Double_t phi1, Double_t phi2) const {
 356     // Calculate azimuthal angle between phi1 and phi2. range:[-0.5\pi,1.5\pi]
 357     Double_t dPhi = phi1-phi2;
 358     if(dPhi <-0.5*TMath::Pi())  dPhi += TMath::TwoPi();
 359     if(dPhi > 1.5*TMath::Pi())  dPhi -= TMath::TwoPi();
 360     return dPhi;
 361   }
 362
 363
 364   //________________________________________________________________________
 365   Bool_t AliAnalysisTaskEmcalHJetMass::RetrieveEventObjects() {
 366     //
 367     // retrieve event objects
 368     //
 369     if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 370       return kFALSE;
 371
 372     return kTRUE;
 373   }
 374
 375   //________________________________________________________________________
 376   void AliAnalysisTaskEmcalHJetMass::AddTriggerTrackPtCuts(Float_t min, Float_t max) {
 377     if(!fPtTTMin) fPtTTMin = new TArrayF();
 378     if(!fPtTTMax) fPtTTMax = new TArrayF();
 379     Int_t newSize = fPtTTMin->GetSize()+1;
 380     fPtTTMin->Set(newSize);
 381     fPtTTMax->Set(newSize);
 382     fPtTTMin->AddAt(min,newSize-1);
 383     fPtTTMax->AddAt(max,newSize-1);
 384   }
 385
 386   //_______________________________________________________________________
 387   void AliAnalysisTaskEmcalHJetMass::Terminate(Option_t *)
 388   {
 389     // Called once at the end of the analysis.
 390   }
 391 }
 392