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