PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetMass.cxx

   1 //
   2 // Jet mass 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
  19 #include "AliVCluster.h"
  20 #include "AliVTrack.h"
  21 #include "AliEmcalJet.h"
  22 #include "AliRhoParameter.h"
  23 #include "AliLog.h"
  24 #include "AliEmcalParticle.h"
  25 #include "AliMCEvent.h"
  26 #include "AliGenPythiaEventHeader.h"
  27 #include "AliAODMCHeader.h"
  28 #include "AliMCEvent.h"
  29 #include "AliAnalysisManager.h"
  30 #include "AliJetContainer.h"
  31
  32 #include "AliAODEvent.h"
  33
  34 #include "AliAnalysisTaskEmcalJetMass.h"
  35
  36 ClassImp(AliAnalysisTaskEmcalJetMass)
  37
  38 //________________________________________________________________________
  39 AliAnalysisTaskEmcalJetMass::AliAnalysisTaskEmcalJetMass() :
  40   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetMass", kTRUE),
  41   fContainerBase(0),
  42   fMinFractionShared(0),
  43   fJetMassType(kRaw),
  44   fh3PtJet1VsMassVsLeadPtAllSel(0),
  45   fh3PtJet1VsMassVsLeadPtTagged(0),
  46   fh3PtJet1VsMassVsLeadPtTaggedMatch(0),
  47   fpPtVsMassJet1All(0),
  48   fpPtVsMassJet1Tagged(0),
  49   fpPtVsMassJet1TaggedMatch(0),
  50   fh2MassVsAreaJet1All(0),
  51   fh2MassVsAreaJet1Tagged(0),
  52   fh2MassVsAreaJet1TaggedMatch(0),
  53   fh2MassVsNConstJet1All(0),
  54   fh2MassVsNConstJet1Tagged(0),
  55   fh2MassVsNConstJet1TaggedMatch(0),
  56   fh2EtMassOverEtRSq(0)
  57 {
  58   // Default constructor.
  59
  60   fh3PtJet1VsMassVsLeadPtAllSel        = new TH3F*[fNcentBins];
  61   fh3PtJet1VsMassVsLeadPtTagged        = new TH3F*[fNcentBins];
  62   fh3PtJet1VsMassVsLeadPtTaggedMatch   = new TH3F*[fNcentBins];
  63   fpPtVsMassJet1All                    = new TProfile*[fNcentBins];
  64   fpPtVsMassJet1Tagged                 = new TProfile*[fNcentBins];
  65   fpPtVsMassJet1TaggedMatch            = new TProfile*[fNcentBins];
  66   fh2MassVsAreaJet1All                 = new TH2F*[fNcentBins];
  67   fh2MassVsAreaJet1Tagged              = new TH2F*[fNcentBins];
  68   fh2MassVsAreaJet1TaggedMatch         = new TH2F*[fNcentBins];
  69   fh2MassVsNConstJet1All               = new TH2F*[fNcentBins];
  70   fh2MassVsNConstJet1Tagged            = new TH2F*[fNcentBins];
  71   fh2MassVsNConstJet1TaggedMatch       = new TH2F*[fNcentBins];
  72   fh2EtMassOverEtRSq                   = new TH2F*[fNcentBins];
  73
  74   for (Int_t i = 0; i < fNcentBins; i++) {
  75     fh3PtJet1VsMassVsLeadPtAllSel[i]        = 0;
  76     fh3PtJet1VsMassVsLeadPtTagged[i]        = 0;
  77     fh3PtJet1VsMassVsLeadPtTaggedMatch[i]   = 0;
  78     fpPtVsMassJet1All[i]                    = 0;
  79     fpPtVsMassJet1Tagged[i]                 = 0;
  80     fpPtVsMassJet1TaggedMatch[i]            = 0;
  81     fh2MassVsAreaJet1All[i]                 = 0;
  82     fh2MassVsAreaJet1Tagged[i]              = 0;
  83     fh2MassVsAreaJet1TaggedMatch[i]         = 0;
  84     fh2MassVsNConstJet1All[i]               = 0;
  85     fh2MassVsNConstJet1Tagged[i]            = 0;
  86     fh2MassVsNConstJet1TaggedMatch[i]       = 0;
  87     fh2EtMassOverEtRSq[i]                   = 0;
  88   }
  89
  90   SetMakeGeneralHistograms(kTRUE);
  91 }
  92
  93 //________________________________________________________________________
  94 AliAnalysisTaskEmcalJetMass::AliAnalysisTaskEmcalJetMass(const char *name) :
  95   AliAnalysisTaskEmcalJet(name, kTRUE),
  96   fContainerBase(0),
  97   fMinFractionShared(0),
  98   fJetMassType(kRaw),
  99   fh3PtJet1VsMassVsLeadPtAllSel(0),
 100   fh3PtJet1VsMassVsLeadPtTagged(0),
 101   fh3PtJet1VsMassVsLeadPtTaggedMatch(0),
 102   fpPtVsMassJet1All(0),
 103   fpPtVsMassJet1Tagged(0),
 104   fpPtVsMassJet1TaggedMatch(0),
 105   fh2MassVsAreaJet1All(0),
 106   fh2MassVsAreaJet1Tagged(0),
 107   fh2MassVsAreaJet1TaggedMatch(0),
 108   fh2MassVsNConstJet1All(0),
 109   fh2MassVsNConstJet1Tagged(0),
 110   fh2MassVsNConstJet1TaggedMatch(0),
 111   fh2EtMassOverEtRSq(0)
 112 {
 113   // Standard constructor.
 114
 115   fh3PtJet1VsMassVsLeadPtAllSel        = new TH3F*[fNcentBins];
 116   fh3PtJet1VsMassVsLeadPtTagged        = new TH3F*[fNcentBins];
 117   fh3PtJet1VsMassVsLeadPtTaggedMatch   = new TH3F*[fNcentBins];
 118   fpPtVsMassJet1All                    = new TProfile*[fNcentBins];
 119   fpPtVsMassJet1Tagged                 = new TProfile*[fNcentBins];
 120   fpPtVsMassJet1TaggedMatch            = new TProfile*[fNcentBins];
 121   fh2MassVsAreaJet1All                 = new TH2F*[fNcentBins];
 122   fh2MassVsAreaJet1Tagged              = new TH2F*[fNcentBins];
 123   fh2MassVsAreaJet1TaggedMatch         = new TH2F*[fNcentBins];
 124   fh2MassVsNConstJet1All               = new TH2F*[fNcentBins];
 125   fh2MassVsNConstJet1Tagged            = new TH2F*[fNcentBins];
 126   fh2MassVsNConstJet1TaggedMatch       = new TH2F*[fNcentBins];
 127   fh2EtMassOverEtRSq                   = new TH2F*[fNcentBins];
 128
 129   for (Int_t i = 0; i < fNcentBins; i++) {
 130     fh3PtJet1VsMassVsLeadPtAllSel[i]        = 0;
 131     fh3PtJet1VsMassVsLeadPtTagged[i]        = 0;
 132     fh3PtJet1VsMassVsLeadPtTaggedMatch[i]   = 0;
 133     fpPtVsMassJet1All[i]                    = 0;
 134     fpPtVsMassJet1Tagged[i]                 = 0;
 135     fpPtVsMassJet1TaggedMatch[i]            = 0;
 136     fh2MassVsAreaJet1All[i]                 = 0;
 137     fh2MassVsAreaJet1Tagged[i]              = 0;
 138     fh2MassVsAreaJet1TaggedMatch[i]         = 0;
 139     fh2MassVsNConstJet1All[i]               = 0;
 140     fh2MassVsNConstJet1Tagged[i]            = 0;
 141     fh2MassVsNConstJet1TaggedMatch[i]       = 0;
 142     fh2EtMassOverEtRSq[i]                   = 0;
 143   }
 144
 145   SetMakeGeneralHistograms(kTRUE);
 146 }
 147
 148 //________________________________________________________________________
 149 AliAnalysisTaskEmcalJetMass::~AliAnalysisTaskEmcalJetMass()
 150 {
 151   // Destructor.
 152 }
 153
 154 //________________________________________________________________________
 155 void AliAnalysisTaskEmcalJetMass::UserCreateOutputObjects()
 156 {
 157   // Create user output.
 158
 159   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 160
 161   Bool_t oldStatus = TH1::AddDirectoryStatus();
 162   TH1::AddDirectory(kFALSE);
 163
 164   const Int_t nBinsPt  = 200;
 165   const Double_t minPt = -50.;
 166   const Double_t maxPt = 150.;
 167
 168   const Int_t nBinsM  = 100;
 169   const Double_t minM = -20.;
 170   const Double_t maxM = 80.;
 171
 172   const Int_t nBinsArea = 50;
 173   const Double_t minArea = 0.;
 174   const Double_t maxArea = 1.;
 175
 176   const Int_t nBinsNConst = 100;
 177   const Double_t minNConst = 0.;
 178   const Double_t maxNConst = 500.;
 179
 180   TString histName = "";
 181   TString histTitle = "";
 182   for (Int_t i = 0; i < fNcentBins; i++) {
 183     histName = TString::Format("fh3PtJet1VsMassVsLeadPtAllSel_%d",i);
 184     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,lead trk}",histName.Data());
 185     fh3PtJet1VsMassVsLeadPtAllSel[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,20,0.,20.);
 186     fOutput->Add(fh3PtJet1VsMassVsLeadPtAllSel[i]);
 187
 188     histName = TString::Format("fh3PtJet1VsMassVsLeadPtTagged_%d",i);
 189     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,lead trk}",histName.Data());
 190     fh3PtJet1VsMassVsLeadPtTagged[i] =  new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,20,0.,20.);
 191     fOutput->Add(fh3PtJet1VsMassVsLeadPtTagged[i]);
 192
 193     histName = TString::Format("fh3PtJet1VsMassVsLeadPtTaggedMatch_%d",i);
 194     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1};#it{p}_{T,lead trk}",histName.Data());
 195     fh3PtJet1VsMassVsLeadPtTaggedMatch[i] =  new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM,20,0.,20.);
 196     fOutput->Add(fh3PtJet1VsMassVsLeadPtTaggedMatch[i]);
 197
 198     histName = TString::Format("fpPtVsMassJet1All_%d",i);
 199     histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
 200     fpPtVsMassJet1All[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
 201     fOutput->Add(fpPtVsMassJet1All[i]);
 202
 203     histName = TString::Format("fpPtVsMassJet1Tagged_%d",i);
 204     histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
 205     fpPtVsMassJet1Tagged[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
 206     fOutput->Add(fpPtVsMassJet1Tagged[i]);
 207
 208     histName = TString::Format("fpPtVsMassJet1TaggedMatch_%d",i);
 209     histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
 210     fpPtVsMassJet1TaggedMatch[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
 211     fOutput->Add(fpPtVsMassJet1TaggedMatch[i]);
 212
 213     histName = TString::Format("fh2MassVsAreaJet1All_%d",i);
 214     histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
 215     fh2MassVsAreaJet1All[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsArea,minArea,maxArea);
 216     fOutput->Add(fh2MassVsAreaJet1All[i]);
 217
 218     histName = TString::Format("fh2MassVsAreaJet1Tagged_%d",i);
 219     histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
 220     fh2MassVsAreaJet1Tagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsArea,minArea,maxArea);
 221     fOutput->Add(fh2MassVsAreaJet1Tagged[i]);
 222
 223     histName = TString::Format("fh2MassVsAreaJet1TaggedMatch_%d",i);
 224     histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
 225     fh2MassVsAreaJet1TaggedMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsArea,minArea,maxArea);
 226     fOutput->Add(fh2MassVsAreaJet1TaggedMatch[i]);
 227
 228     histName = TString::Format("fh2MassVsNConstJet1All_%d",i);
 229     histTitle = TString::Format("%s;#it{M}_{jet1};#it{N}_{constituents}",histName.Data());
 230     fh2MassVsNConstJet1All[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsNConst,minNConst,maxNConst);
 231     fOutput->Add(fh2MassVsNConstJet1All[i]);
 232
 233     histName = TString::Format("fh2MassVsNConstJet1Tagged_%d",i);
 234     histTitle = TString::Format("%s;#it{M}_{jet1};#it{N}_{constituents}",histName.Data());
 235     fh2MassVsNConstJet1Tagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsNConst,minNConst,maxNConst);
 236     fOutput->Add(fh2MassVsNConstJet1Tagged[i]);
 237
 238     histName = TString::Format("fh2MassVsNConstJet1TaggedMatch_%d",i);
 239     histTitle = TString::Format("%s;#it{M}_{jet1};#it{N}_{constituents}",histName.Data());
 240     fh2MassVsNConstJet1TaggedMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsNConst,minNConst,maxNConst);
 241     fOutput->Add(fh2MassVsNConstJet1TaggedMatch[i]);
 242
 243     histName = TString::Format("fh2EtMassOverEtRSq_%d",i);
 244     histTitle = TString::Format("%s;#it{E}_{T};(#it{M}/(#it{E}_{T}#it{R}))^{2}",histName.Data());
 245     fh2EtMassOverEtRSq[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,100,0.,1.);
 246     fOutput->Add(fh2EtMassOverEtRSq[i]);
 247   }
 248
 249   // =========== Switch on Sumw2 for all histos ===========
 250   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 251     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 252     if (h1){
 253       h1->Sumw2();
 254       continue;
 255     }
 256     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 257     if(hn)hn->Sumw2();
 258   }
 259
 260   TH1::AddDirectory(oldStatus);
 261
 262   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 263 }
 264
 265 //________________________________________________________________________
 266 Bool_t AliAnalysisTaskEmcalJetMass::Run()
 267 {
 268   // Run analysis code here, if needed. It will be executed before FillHistograms().
 269
 270   return kTRUE;
 271 }
 272
 273 //________________________________________________________________________
 274 Bool_t AliAnalysisTaskEmcalJetMass::FillHistograms()
 275 {
 276   // Fill histograms.
 277
 278   AliEmcalJet* jet1 = NULL;
 279   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 280   if(jetCont) {
 281     jetCont->ResetCurrentID();
 282     while((jet1 = jetCont->GetNextAcceptJet())) {
 283
 284       Double_t ptJet1 = jet1->Pt() - GetRhoVal(fContainerBase)*jet1->Area();
 285       Double_t mJet1 = GetJetMass(jet1);
 286       fh3PtJet1VsMassVsLeadPtAllSel[fCentBin]->Fill(ptJet1,mJet1,jet1->MaxTrackPt());
 287       fpPtVsMassJet1All[fCentBin]->Fill(ptJet1,mJet1);
 288       fh2MassVsAreaJet1All[fCentBin]->Fill(mJet1,jet1->Area());
 289       fh2MassVsNConstJet1All[fCentBin]->Fill(mJet1,jet1->GetNumberOfConstituents());
 290
 291       if(jet1->GetTagStatus()<1 || !jet1->GetTaggedJet())
 292         continue;
 293
 294       fh3PtJet1VsMassVsLeadPtTagged[fCentBin]->Fill(ptJet1,mJet1,jet1->MaxTrackPt());
 295       fpPtVsMassJet1Tagged[fCentBin]->Fill(ptJet1,mJet1);
 296       fh2MassVsAreaJet1Tagged[fCentBin]->Fill(mJet1,jet1->Area());
 297       fh2MassVsNConstJet1Tagged[fCentBin]->Fill(mJet1,jet1->GetNumberOfConstituents());
 298
 299       Double_t fraction = jetCont->GetFractionSharedPt(jet1);
 300       if(fMinFractionShared>0. && fraction>fMinFractionShared) {
 301         fh3PtJet1VsMassVsLeadPtTaggedMatch[fCentBin]->Fill(ptJet1,mJet1,jet1->MaxTrackPt());
 302         fpPtVsMassJet1TaggedMatch[fCentBin]->Fill(ptJet1,mJet1);
 303         fh2MassVsAreaJet1TaggedMatch[fCentBin]->Fill(mJet1,jet1->Area());
 304         fh2MassVsNConstJet1TaggedMatch[fCentBin]->Fill(mJet1,jet1->GetNumberOfConstituents());
 305       }
 306
 307       Double_t Et2 = mJet1*mJet1 + ptJet1*ptJet1;
 308       Double_t Et = 0.;    Double_t massOverEtR = 0.;
 309       if(Et2>0.) Et = TMath::Sqrt(Et2);
 310       if((Et*jetCont->GetJetRadius())>0.)
 311         massOverEtR = mJet1/(Et*jetCont->GetJetRadius());
 312       fh2EtMassOverEtRSq[fCentBin]->Fill(Et,massOverEtR*massOverEtR);
 313     }
 314   }
 315
 316   return kTRUE;
 317 }
 318
 319 //________________________________________________________________________
 320 Double_t AliAnalysisTaskEmcalJetMass::GetJetMass(AliEmcalJet *jet) {
 321   //calc subtracted jet mass
 322   if(fJetMassType==kRaw)
 323     return jet->M();
 324   else if(fJetMassType==kDeriv)
 325     return jet->GetSecondOrderSubtracted();
 326
 327   return 0;
 328 }
 329
 330 //________________________________________________________________________
 331 Bool_t AliAnalysisTaskEmcalJetMass::RetrieveEventObjects() {
 332   //
 333   // retrieve event objects
 334   //
 335   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 336     return kFALSE;
 337
 338   return kTRUE;
 339 }
 340
 341 //_______________________________________________________________________
 342 void AliAnalysisTaskEmcalJetMass::Terminate(Option_t *)
 343 {
 344   // Called once at the end of the analysis.
 345 }
 346