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   fh2PtJet1VsLeadPtAllSel(0),
  43   fh2PtJet1VsLeadPtTagged(0),
  44   fh2PtVsMassJet1All(0),
  45   fh2PtVsMassJet1Tagged(0),
  46   fpPtVsMassJet1All(0),
  47   fpPtVsMassJet1Tagged(0),
  48   fh2MassVsAreaJet1All(0),
  49   fh2MassVsAreaJet1Tagged(0)
  50 {
  51   // Default constructor.
  52
  53   fh2PtJet1VsLeadPtAllSel      = new TH2F*[fNcentBins];
  54   fh2PtJet1VsLeadPtTagged      = new TH2F*[fNcentBins];
  55   fh2PtVsMassJet1All           = new TH2F*[fNcentBins];
  56   fh2PtVsMassJet1Tagged        = new TH2F*[fNcentBins];
  57   fpPtVsMassJet1All            = new TProfile*[fNcentBins];
  58   fpPtVsMassJet1Tagged         = new TProfile*[fNcentBins];
  59   fh2MassVsAreaJet1All         = new TH2F*[fNcentBins];
  60   fh2MassVsAreaJet1Tagged      = new TH2F*[fNcentBins];
  61
  62   for (Int_t i = 0; i < fNcentBins; i++) {
  63     fh2PtJet1VsLeadPtAllSel[i]     = 0;
  64     fh2PtJet1VsLeadPtTagged[i]     = 0;
  65     fh2PtVsMassJet1All[i]          = 0;
  66     fh2PtVsMassJet1Tagged[i]       = 0;
  67     fpPtVsMassJet1All[i]           = 0;
  68     fpPtVsMassJet1Tagged[i]        = 0;
  69     fh2MassVsAreaJet1All[i]        = 0;
  70     fh2MassVsAreaJet1Tagged[i]     = 0;
  71   }
  72
  73   SetMakeGeneralHistograms(kTRUE);
  74
  75 }
  76
  77 //________________________________________________________________________
  78 AliAnalysisTaskEmcalJetMass::AliAnalysisTaskEmcalJetMass(const char *name) :
  79   AliAnalysisTaskEmcalJet(name, kTRUE),
  80   fContainerBase(0),
  81   fh2PtJet1VsLeadPtAllSel(0),
  82   fh2PtJet1VsLeadPtTagged(0),
  83   fh2PtVsMassJet1All(0),
  84   fh2PtVsMassJet1Tagged(0),
  85   fpPtVsMassJet1All(0),
  86   fpPtVsMassJet1Tagged(0),
  87   fh2MassVsAreaJet1All(0),
  88   fh2MassVsAreaJet1Tagged(0)
  89 {
  90   // Standard constructor.
  91
  92   fh2PtJet1VsLeadPtAllSel      = new TH2F*[fNcentBins];
  93   fh2PtJet1VsLeadPtTagged      = new TH2F*[fNcentBins];
  94   fh2PtVsMassJet1All           = new TH2F*[fNcentBins];
  95   fh2PtVsMassJet1Tagged        = new TH2F*[fNcentBins];
  96   fpPtVsMassJet1All            = new TProfile*[fNcentBins];
  97   fpPtVsMassJet1Tagged         = new TProfile*[fNcentBins];
  98   fh2MassVsAreaJet1All         = new TH2F*[fNcentBins];
  99   fh2MassVsAreaJet1Tagged      = new TH2F*[fNcentBins];
 100
 101   for (Int_t i = 0; i < fNcentBins; i++) {
 102     fh2PtJet1VsLeadPtAllSel[i]     = 0;
 103     fh2PtJet1VsLeadPtTagged[i]     = 0;
 104     fh2PtVsMassJet1All[i]          = 0;
 105     fh2PtVsMassJet1Tagged[i]       = 0;
 106     fpPtVsMassJet1All[i]           = 0;
 107     fpPtVsMassJet1Tagged[i]        = 0;
 108     fh2MassVsAreaJet1All[i]        = 0;
 109     fh2MassVsAreaJet1Tagged[i]     = 0;
 110   }
 111
 112   SetMakeGeneralHistograms(kTRUE);
 113 }
 114
 115 //________________________________________________________________________
 116 AliAnalysisTaskEmcalJetMass::~AliAnalysisTaskEmcalJetMass()
 117 {
 118   // Destructor.
 119 }
 120
 121 //________________________________________________________________________
 122 void AliAnalysisTaskEmcalJetMass::UserCreateOutputObjects()
 123 {
 124   // Create user output.
 125
 126   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 127
 128   Bool_t oldStatus = TH1::AddDirectoryStatus();
 129   TH1::AddDirectory(kFALSE);
 130
 131   const Int_t nBinsPt  = 250;
 132   const Double_t minPt = -50.;
 133   const Double_t maxPt = 200.;
 134
 135   const Int_t nBinsArea = 100;
 136   const Double_t minArea = 0.;
 137   const Double_t maxArea = 1.;
 138
 139   TString histName = "";
 140   TString histTitle = "";
 141   for (Int_t i = 0; i < fNcentBins; i++) {
 142     histName = TString::Format("fh2PtJet1VsLeadPtAllSel_%d",i);
 143     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
 144     fh2PtJet1VsLeadPtAllSel[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
 145     fOutput->Add(fh2PtJet1VsLeadPtAllSel[i]);
 146
 147     histName = TString::Format("fh2PtJet1VsLeadPtTagged_%d",i);
 148     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
 149     fh2PtJet1VsLeadPtTagged[i] =  new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
 150     fOutput->Add(fh2PtJet1VsLeadPtTagged[i]);
 151
 152     histName = TString::Format("fh2PtVsMassJet1All_%d",i);
 153     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
 154     fh2PtVsMassJet1All[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 155     fOutput->Add(fh2PtVsMassJet1All[i]);
 156
 157     histName = TString::Format("fh2PtVsMassJet1Tagged_%d",i);
 158     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
 159     fh2PtVsMassJet1Tagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 160     fOutput->Add(fh2PtVsMassJet1Tagged[i]);
 161
 162     histName = TString::Format("fpPtVsMassJet1All_%d",i);
 163     histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
 164     fpPtVsMassJet1All[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
 165     fOutput->Add(fpPtVsMassJet1All[i]);
 166
 167     histName = TString::Format("fpPtVsMassJet1Tagged_%d",i);
 168     histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
 169     fpPtVsMassJet1Tagged[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
 170     fOutput->Add(fpPtVsMassJet1Tagged[i]);
 171
 172     histName = TString::Format("fh2MassVsAreaJet1All_%d",i);
 173     histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
 174     fh2MassVsAreaJet1All[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsArea,minArea,maxArea);
 175     fOutput->Add(fh2MassVsAreaJet1All[i]);
 176
 177     histName = TString::Format("fh2MassVsAreaJet1Tagged_%d",i);
 178     histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
 179     fh2MassVsAreaJet1Tagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsArea,minArea,maxArea);
 180     fOutput->Add(fh2MassVsAreaJet1Tagged[i]);
 181
 182   }
 183
 184   // =========== Switch on Sumw2 for all histos ===========
 185   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 186     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 187     if (h1){
 188       h1->Sumw2();
 189       continue;
 190     }
 191     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 192     if(hn)hn->Sumw2();
 193   }
 194
 195   TH1::AddDirectory(oldStatus);
 196
 197   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 198 }
 199
 200 //________________________________________________________________________
 201 Bool_t AliAnalysisTaskEmcalJetMass::Run()
 202 {
 203   // Run analysis code here, if needed. It will be executed before FillHistograms().
 204
 205   return kTRUE;
 206 }
 207
 208 //________________________________________________________________________
 209 Bool_t AliAnalysisTaskEmcalJetMass::FillHistograms()
 210 {
 211   // Fill histograms.
 212
 213   for(int i = 0; i < GetNJets(fContainerBase);++i) {
 214     AliEmcalJet *jet1 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, fContainerBase));
 215     if(!jet1) continue;
 216
 217     Double_t ptJet1 =  jet1->Pt() - GetRhoVal(fContainerBase)*jet1->Area();
 218     fh2PtJet1VsLeadPtAllSel[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
 219     fh2PtVsMassJet1All[fCentBin]->Fill(ptJet1,jet1->M());
 220     fpPtVsMassJet1All[fCentBin]->Fill(ptJet1,jet1->M());
 221     fh2MassVsAreaJet1All[fCentBin]->Fill(jet1->M(),jet1->Area());
 222
 223     if(jet1->GetTagStatus()<1)
 224       continue;
 225
 226     AliEmcalJet *jet2 = jet1->GetTaggedJet();
 227     if(!jet2) continue;
 228
 229     fh2PtJet1VsLeadPtTagged[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
 230     fh2PtVsMassJet1Tagged[fCentBin]->Fill(ptJet1,jet1->M());
 231     fpPtVsMassJet1Tagged[fCentBin]->Fill(ptJet1,jet1->M());
 232     fh2MassVsAreaJet1Tagged[fCentBin]->Fill(jet1->M(),jet1->Area());
 233
 234   }
 235
 236   return kTRUE;
 237
 238 }
 239
 240
 241 //________________________________________________________________________
 242 Bool_t AliAnalysisTaskEmcalJetMass::RetrieveEventObjects() {
 243   //
 244   // retrieve event objects
 245   //
 246
 247   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 248     return kFALSE;
 249
 250   return kTRUE;
 251
 252 }
 253
 254 //_______________________________________________________________________
 255 void AliAnalysisTaskEmcalJetMass::Terminate(Option_t *)
 256 {
 257   // Called once at the end of the analysis.
 258 }
 259