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   fh2PtJet1VsLeadPtTaggedMatch(0),
  47   fh2PtVsMassJet1All(0),
  48   fh2PtVsMassJet1Tagged(0),
  49   fh2PtVsMassJet1TaggedMatch(0),
  50   fpPtVsMassJet1All(0),
  51   fpPtVsMassJet1Tagged(0),
  52   fpPtVsMassJet1TaggedMatch(0),
  53   fh2MassVsAreaJet1All(0),
  54   fh2MassVsAreaJet1Tagged(0),
  55   fh2MassVsAreaJet1TaggedMatch(0),
  56   fh2MassVsNConstJet1All(0),
  57   fh2MassVsNConstJet1Tagged(0),
  58   fh2MassVsNConstJet1TaggedMatch(0),
  59   fh2EtMassOverEtRSq(0)
  60 {
  61   // Default constructor.
  62
  63   fh2PtJet1VsLeadPtAllSel        = new TH2F*[fNcentBins];
  64   fh2PtJet1VsLeadPtTagged        = new TH2F*[fNcentBins];
  65   fh2PtJet1VsLeadPtTaggedMatch   = new TH2F*[fNcentBins];
  66   fh2PtVsMassJet1All             = new TH2F*[fNcentBins];
  67   fh2PtVsMassJet1Tagged          = new TH2F*[fNcentBins];
  68   fh2PtVsMassJet1TaggedMatch     = new TH2F*[fNcentBins];
  69   fpPtVsMassJet1All              = new TProfile*[fNcentBins];
  70   fpPtVsMassJet1Tagged           = new TProfile*[fNcentBins];
  71   fpPtVsMassJet1TaggedMatch      = new TProfile*[fNcentBins];
  72   fh2MassVsAreaJet1All           = new TH2F*[fNcentBins];
  73   fh2MassVsAreaJet1Tagged        = new TH2F*[fNcentBins];
  74   fh2MassVsAreaJet1TaggedMatch   = new TH2F*[fNcentBins];
  75   fh2MassVsNConstJet1All         = new TH2F*[fNcentBins];
  76   fh2MassVsNConstJet1Tagged      = new TH2F*[fNcentBins];
  77   fh2MassVsNConstJet1TaggedMatch = new TH2F*[fNcentBins];
  78   fh2EtMassOverEtRSq             = new TH2F*[fNcentBins];
  79
  80   for (Int_t i = 0; i < fNcentBins; i++) {
  81     fh2PtJet1VsLeadPtAllSel[i]        = 0;
  82     fh2PtJet1VsLeadPtTagged[i]        = 0;
  83     fh2PtJet1VsLeadPtTaggedMatch[i]   = 0;
  84     fh2PtVsMassJet1All[i]             = 0;
  85     fh2PtVsMassJet1Tagged[i]          = 0;
  86     fh2PtVsMassJet1TaggedMatch[i]     = 0;
  87     fpPtVsMassJet1All[i]              = 0;
  88     fpPtVsMassJet1Tagged[i]           = 0;
  89     fpPtVsMassJet1TaggedMatch[i]      = 0;
  90     fh2MassVsAreaJet1All[i]           = 0;
  91     fh2MassVsAreaJet1Tagged[i]        = 0;
  92     fh2MassVsAreaJet1TaggedMatch[i]   = 0;
  93     fh2MassVsNConstJet1All[i]         = 0;
  94     fh2MassVsNConstJet1Tagged[i]      = 0;
  95     fh2MassVsNConstJet1TaggedMatch[i] = 0;
  96     fh2EtMassOverEtRSq[i]             = 0;
  97   }
  98
  99   SetMakeGeneralHistograms(kTRUE);
 100 }
 101
 102 //________________________________________________________________________
 103 AliAnalysisTaskEmcalJetMass::AliAnalysisTaskEmcalJetMass(const char *name) :
 104   AliAnalysisTaskEmcalJet(name, kTRUE),
 105   fContainerBase(0),
 106   fMinFractionShared(0),
 107   fJetMassAvg(0),
 108   fh2PtJet1VsLeadPtAllSel(0),
 109   fh2PtJet1VsLeadPtTagged(0),
 110   fh2PtJet1VsLeadPtTaggedMatch(0),
 111   fh2PtVsMassJet1All(0),
 112   fh2PtVsMassJet1Tagged(0),
 113   fh2PtVsMassJet1TaggedMatch(0),
 114   fpPtVsMassJet1All(0),
 115   fpPtVsMassJet1Tagged(0),
 116   fpPtVsMassJet1TaggedMatch(0),
 117   fh2MassVsAreaJet1All(0),
 118   fh2MassVsAreaJet1Tagged(0),
 119   fh2MassVsAreaJet1TaggedMatch(0),
 120   fh2MassVsNConstJet1All(0),
 121   fh2MassVsNConstJet1Tagged(0),
 122   fh2MassVsNConstJet1TaggedMatch(0),
 123   fh2EtMassOverEtRSq(0)
 124 {
 125   // Standard constructor.
 126
 127  fh2PtJet1VsLeadPtAllSel        = new TH2F*[fNcentBins];
 128   fh2PtJet1VsLeadPtTagged        = new TH2F*[fNcentBins];
 129   fh2PtJet1VsLeadPtTaggedMatch   = new TH2F*[fNcentBins];
 130   fh2PtVsMassJet1All             = new TH2F*[fNcentBins];
 131   fh2PtVsMassJet1Tagged          = new TH2F*[fNcentBins];
 132   fh2PtVsMassJet1TaggedMatch     = new TH2F*[fNcentBins];
 133   fpPtVsMassJet1All              = new TProfile*[fNcentBins];
 134   fpPtVsMassJet1Tagged           = new TProfile*[fNcentBins];
 135   fpPtVsMassJet1TaggedMatch      = new TProfile*[fNcentBins];
 136   fh2MassVsAreaJet1All           = new TH2F*[fNcentBins];
 137   fh2MassVsAreaJet1Tagged        = new TH2F*[fNcentBins];
 138   fh2MassVsAreaJet1TaggedMatch   = new TH2F*[fNcentBins];
 139   fh2MassVsNConstJet1All         = new TH2F*[fNcentBins];
 140   fh2MassVsNConstJet1Tagged      = new TH2F*[fNcentBins];
 141   fh2MassVsNConstJet1TaggedMatch = new TH2F*[fNcentBins];
 142   fh2EtMassOverEtRSq             = new TH2F*[fNcentBins];
 143
 144   for (Int_t i = 0; i < fNcentBins; i++) {
 145     fh2PtJet1VsLeadPtAllSel[i]        = 0;
 146     fh2PtJet1VsLeadPtTagged[i]        = 0;
 147     fh2PtJet1VsLeadPtTaggedMatch[i]   = 0;
 148     fh2PtVsMassJet1All[i]             = 0;
 149     fh2PtVsMassJet1Tagged[i]          = 0;
 150     fh2PtVsMassJet1TaggedMatch[i]     = 0;
 151     fpPtVsMassJet1All[i]              = 0;
 152     fpPtVsMassJet1Tagged[i]           = 0;
 153     fpPtVsMassJet1TaggedMatch[i]      = 0;
 154     fh2MassVsAreaJet1All[i]           = 0;
 155     fh2MassVsAreaJet1Tagged[i]        = 0;
 156     fh2MassVsAreaJet1TaggedMatch[i]   = 0;
 157     fh2MassVsNConstJet1All[i]         = 0;
 158     fh2MassVsNConstJet1Tagged[i]      = 0;
 159     fh2MassVsNConstJet1TaggedMatch[i] = 0;
 160     fh2EtMassOverEtRSq[i]             = 0;
 161   }
 162
 163   SetMakeGeneralHistograms(kTRUE);
 164 }
 165
 166 //________________________________________________________________________
 167 AliAnalysisTaskEmcalJetMass::~AliAnalysisTaskEmcalJetMass()
 168 {
 169   // Destructor.
 170 }
 171
 172 //________________________________________________________________________
 173 void AliAnalysisTaskEmcalJetMass::UserCreateOutputObjects()
 174 {
 175   // Create user output.
 176
 177   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 178
 179   Bool_t oldStatus = TH1::AddDirectoryStatus();
 180   TH1::AddDirectory(kFALSE);
 181
 182   const Int_t nBinsPt  = 250;
 183   const Double_t minPt = -50.;
 184   const Double_t maxPt = 200.;
 185
 186   const Int_t nBinsM  = 150;
 187   const Double_t minM = -50.;
 188   const Double_t maxM = 100.;
 189
 190   const Int_t nBinsArea = 100;
 191   const Double_t minArea = 0.;
 192   const Double_t maxArea = 1.;
 193
 194   const Int_t nBinsNConst = 100;
 195   const Double_t minNConst = 0.;
 196   const Double_t maxNConst = 500.;
 197
 198   TString histName = "";
 199   TString histTitle = "";
 200   for (Int_t i = 0; i < fNcentBins; i++) {
 201     histName = TString::Format("fh2PtJet1VsLeadPtAllSel_%d",i);
 202     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
 203     fh2PtJet1VsLeadPtAllSel[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
 204     fOutput->Add(fh2PtJet1VsLeadPtAllSel[i]);
 205
 206     histName = TString::Format("fh2PtJet1VsLeadPtTagged_%d",i);
 207     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
 208     fh2PtJet1VsLeadPtTagged[i] =  new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
 209     fOutput->Add(fh2PtJet1VsLeadPtTagged[i]);
 210
 211     histName = TString::Format("fh2PtJet1VsLeadPtTaggedMatch_%d",i);
 212     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
 213     fh2PtJet1VsLeadPtTaggedMatch[i] =  new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
 214     fOutput->Add(fh2PtJet1VsLeadPtTaggedMatch[i]);
 215
 216     histName = TString::Format("fh2PtVsMassJet1All_%d",i);
 217     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
 218     fh2PtVsMassJet1All[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 219     fOutput->Add(fh2PtVsMassJet1All[i]);
 220
 221     histName = TString::Format("fh2PtVsMassJet1Tagged_%d",i);
 222     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
 223     fh2PtVsMassJet1Tagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 224     fOutput->Add(fh2PtVsMassJet1Tagged[i]);
 225
 226     histName = TString::Format("fh2PtVsMassJet1TaggedMatch_%d",i);
 227     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
 228     fh2PtVsMassJet1TaggedMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 229     fOutput->Add(fh2PtVsMassJet1TaggedMatch[i]);
 230
 231     histName = TString::Format("fpPtVsMassJet1All_%d",i);
 232     histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
 233     fpPtVsMassJet1All[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsM,minM,maxM);
 234     fOutput->Add(fpPtVsMassJet1All[i]);
 235
 236     histName = TString::Format("fpPtVsMassJet1Tagged_%d",i);
 237     histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
 238     fpPtVsMassJet1Tagged[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsM,minM,maxM);
 239     fOutput->Add(fpPtVsMassJet1Tagged[i]);
 240
 241     histName = TString::Format("fpPtVsMassJet1TaggedMatch_%d",i);
 242     histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
 243     fpPtVsMassJet1TaggedMatch[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsM,minM,maxM);
 244     fOutput->Add(fpPtVsMassJet1TaggedMatch[i]);
 245
 246     histName = TString::Format("fh2MassVsAreaJet1All_%d",i);
 247     histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
 248     fh2MassVsAreaJet1All[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsArea,minArea,maxArea);
 249     fOutput->Add(fh2MassVsAreaJet1All[i]);
 250
 251     histName = TString::Format("fh2MassVsAreaJet1Tagged_%d",i);
 252     histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
 253     fh2MassVsAreaJet1Tagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsArea,minArea,maxArea);
 254     fOutput->Add(fh2MassVsAreaJet1Tagged[i]);
 255
 256     histName = TString::Format("fh2MassVsAreaJet1TaggedMatch_%d",i);
 257     histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
 258     fh2MassVsAreaJet1TaggedMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsArea,minArea,maxArea);
 259     fOutput->Add(fh2MassVsAreaJet1TaggedMatch[i]);
 260
 261     histName = TString::Format("fh2MassVsNConstJet1All_%d",i);
 262     histTitle = TString::Format("%s;#it{M}_{jet1};#it{N}_{constituents}",histName.Data());
 263     fh2MassVsNConstJet1All[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsNConst,minNConst,maxNConst);
 264     fOutput->Add(fh2MassVsNConstJet1All[i]);
 265
 266     histName = TString::Format("fh2MassVsNConstJet1Tagged_%d",i);
 267     histTitle = TString::Format("%s;#it{M}_{jet1};#it{N}_{constituents}",histName.Data());
 268     fh2MassVsNConstJet1Tagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsNConst,minNConst,maxNConst);
 269     fOutput->Add(fh2MassVsNConstJet1Tagged[i]);
 270
 271     histName = TString::Format("fh2MassVsNConstJet1TaggedMatch_%d",i);
 272     histTitle = TString::Format("%s;#it{M}_{jet1};#it{N}_{constituents}",histName.Data());
 273     fh2MassVsNConstJet1TaggedMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsNConst,minNConst,maxNConst);
 274     fOutput->Add(fh2MassVsNConstJet1TaggedMatch[i]);
 275
 276     histName = TString::Format("fh2EtMassOverEtRSq_%d",i);
 277     histTitle = TString::Format("%s;#it{E}_{T};(#it{M}/(#it{E}_{T}#it{R}))^{2}",histName.Data());
 278     fh2EtMassOverEtRSq[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,100,0.,1.);
 279     fOutput->Add(fh2EtMassOverEtRSq[i]);
 280   }
 281
 282   // =========== Switch on Sumw2 for all histos ===========
 283   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 284     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 285     if (h1){
 286       h1->Sumw2();
 287       continue;
 288     }
 289     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 290     if(hn)hn->Sumw2();
 291   }
 292
 293   TH1::AddDirectory(oldStatus);
 294
 295   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 296 }
 297
 298 //________________________________________________________________________
 299 Bool_t AliAnalysisTaskEmcalJetMass::Run()
 300 {
 301   // Run analysis code here, if needed. It will be executed before FillHistograms().
 302
 303   return kTRUE;
 304 }
 305
 306 //________________________________________________________________________
 307 Bool_t AliAnalysisTaskEmcalJetMass::FillHistograms()
 308 {
 309   // Fill histograms.
 310
 311   AliEmcalJet* jet1 = NULL;
 312   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 313   if(jetCont) {
 314     jetCont->ResetCurrentID();
 315     while((jet1 = jetCont->GetNextAcceptJet())) {
 316
 317       Double_t ptJet1 = jet1->Pt() - GetRhoVal(fContainerBase)*jet1->Area();
 318       Double_t mJet1 = GetJetMass(jet1);
 319       fh2PtJet1VsLeadPtAllSel[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
 320       fh2PtVsMassJet1All[fCentBin]->Fill(ptJet1,mJet1);
 321       fpPtVsMassJet1All[fCentBin]->Fill(ptJet1,mJet1);
 322       fh2MassVsAreaJet1All[fCentBin]->Fill(mJet1,jet1->Area());
 323       fh2MassVsNConstJet1All[fCentBin]->Fill(mJet1,jet1->GetNumberOfConstituents());
 324
 325       if(jet1->GetTagStatus()<1 || !jet1->GetTaggedJet())
 326         continue;
 327
 328       fh2PtJet1VsLeadPtTagged[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
 329       fh2PtVsMassJet1Tagged[fCentBin]->Fill(ptJet1,mJet1);
 330       fpPtVsMassJet1Tagged[fCentBin]->Fill(ptJet1,mJet1);
 331       fh2MassVsAreaJet1Tagged[fCentBin]->Fill(mJet1,jet1->Area());
 332       fh2MassVsNConstJet1Tagged[fCentBin]->Fill(mJet1,jet1->GetNumberOfConstituents());
 333
 334       Double_t fraction = jetCont->GetFractionSharedPt(jet1);
 335       if(fMinFractionShared>0. && fraction>fMinFractionShared) {
 336         fh2PtJet1VsLeadPtTaggedMatch[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
 337         fh2PtVsMassJet1TaggedMatch[fCentBin]->Fill(ptJet1,mJet1);
 338         fpPtVsMassJet1TaggedMatch[fCentBin]->Fill(ptJet1,mJet1);
 339         fh2MassVsAreaJet1TaggedMatch[fCentBin]->Fill(mJet1,jet1->Area());
 340         fh2MassVsNConstJet1TaggedMatch[fCentBin]->Fill(mJet1,jet1->GetNumberOfConstituents());
 341       }
 342
 343       Double_t Et2 = jet1->M()*jet1->M() + jet1->Pt()*jet1->Pt();
 344       Double_t Et = 0.;    Double_t massOverEtR = 0.;
 345       if(Et2>0.) Et = TMath::Sqrt(Et2);
 346       if((Et*jetCont->GetJetRadius())>0.)
 347         massOverEtR = jet1->M()/(Et*jetCont->GetJetRadius());
 348       fh2EtMassOverEtRSq[fCentBin]->Fill(Et,massOverEtR*massOverEtR);
 349     }
 350   }
 351
 352   return kTRUE;
 353 }
 354
 355 //________________________________________________________________________
 356 Double_t AliAnalysisTaskEmcalJetMass::GetJetMass(AliEmcalJet *jet) {
 357   //calc subtracted jet mass
 358   TLorentzVector vpS = GetSubtractedVector(jet);
 359   return vpS.M();
 360 }
 361
 362
 363 //________________________________________________________________________
 364 TLorentzVector AliAnalysisTaskEmcalJetMass::GetSubtractedVector(AliEmcalJet *jet) {
 365   //get subtracted vector
 366   TLorentzVector vpS;
 367   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 368   TLorentzVector vpBkg = GetBkgVector(jet,jetCont);
 369   vpS.SetPxPyPzE(jet->Px()-vpBkg.Px(),jet->Py()-vpBkg.Py(),jet->Pz()-vpBkg.Pz(),jet->E()-vpBkg.E());
 370   return vpS;
 371 }
 372
 373 //________________________________________________________________________
 374 TLorentzVector AliAnalysisTaskEmcalJetMass::GetBkgVector(AliEmcalJet *jet, AliJetContainer *cont) {
 375   //get background vector
 376
 377   Double_t rho  = cont->GetRhoVal();
 378   Double_t rhom = cont->GetRhoMassVal();
 379   TLorentzVector vpB;
 380   Double_t aphi = jet->AreaPhi();
 381   Double_t aeta = jet->AreaEta();
 382   vpB.SetPxPyPzE(rho*TMath::Cos(aphi)*jet->Area(),rho*TMath::Sin(aphi)*jet->Area(),(rho+rhom)*TMath::SinH(aeta)*jet->Area(),(rho+rhom)*TMath::CosH(aeta)*jet->Area());
 383   return vpB;
 384 }
 385
 386 //________________________________________________________________________
 387 Bool_t AliAnalysisTaskEmcalJetMass::RetrieveEventObjects() {
 388   //
 389   // retrieve event objects
 390   //
 391
 392   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 393     return kFALSE;
 394
 395   return kTRUE;
 396
 397 }
 398
 399 //_______________________________________________________________________
 400 void AliAnalysisTaskEmcalJetMass::Terminate(Option_t *)
 401 {
 402   // Called once at the end of the analysis.
 403 }
 404