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