PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskDeltaPtJEmb.cxx

   1 // $Id$
   2 //
   3 // Jet embedding deltaPt task.
   4 //
   5 // Author: M. Verweij
   6
   7 #include <TClonesArray.h>
   8 #include <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TH3F.h>
  11 #include <TList.h>
  12 #include <TLorentzVector.h>
  13 #include <TRandom3.h>
  14
  15 #include "AliVCluster.h"
  16 #include "AliVParticle.h"
  17 #include "AliVTrack.h"
  18 #include "AliEmcalJet.h"
  19 #include "AliRhoParameter.h"
  20 #include "AliLog.h"
  21 #include "AliJetContainer.h"
  22 #include "AliParticleContainer.h"
  23 #include "AliClusterContainer.h"
  24
  25 #include "AliAnalysisTaskDeltaPtJEmb.h"
  26
  27 ClassImp(AliAnalysisTaskDeltaPtJEmb)
  28
  29 //________________________________________________________________________
  30 AliAnalysisTaskDeltaPtJEmb::AliAnalysisTaskDeltaPtJEmb() :
  31   AliAnalysisTaskEmcalJet("AliAnalysisTaskDeltaPtJEmb", kTRUE),
  32   fJetsCont(0),
  33   fTracksCont(0),
  34   fCaloClustersCont(0),
  35   fMinFractionShared(0.5),
  36   fHistEmbJetsPtArea(0),
  37   fHistEmbJetsCorrPtArea(0),
  38   fHistEmbPartPtvsJetPt(0),
  39   fHistEmbPartPtvsJetCorrPt(0),
  40   fHistJetPtvsJetCorrPt(0),
  41   fHistDistLeadPart2JetAxis(0),
  42   fHistEmbBkgArea(0),
  43   fHistRhoVSEmbBkg(0),
  44   fHistDeltaPtEmbArea(0),
  45   fHistDeltaPtEmbvsEP(0),
  46   fHistDeltaPtEmbPtProbe(0),
  47   fHistEmbJetsPhiEta(0),
  48   fHistLeadPartPhiEta(0)
  49 {
  50   // Default constructor.
  51
  52   fHistEmbJetsPtArea = new TH3*[fNcentBins];
  53   fHistEmbJetsCorrPtArea = new TH3*[fNcentBins];
  54   fHistEmbPartPtvsJetPt = new TH2*[fNcentBins];
  55   fHistEmbPartPtvsJetCorrPt = new TH2*[fNcentBins];
  56   fHistJetPtvsJetCorrPt = new TH2*[fNcentBins];
  57   fHistDistLeadPart2JetAxis = new TH1*[fNcentBins];
  58   fHistEmbBkgArea = new TH2*[fNcentBins];
  59   fHistRhoVSEmbBkg = new TH2*[fNcentBins];
  60   fHistDeltaPtEmbArea = new TH2*[fNcentBins];
  61   fHistDeltaPtEmbvsEP = new TH2*[fNcentBins];
  62   fHistDeltaPtEmbPtProbe = new TH2*[fNcentBins];
  63
  64   for (Int_t i = 0; i < fNcentBins; i++) {
  65     fHistEmbJetsPtArea[i] = 0;
  66     fHistEmbJetsCorrPtArea[i] = 0;
  67     fHistEmbPartPtvsJetPt[i] = 0;
  68     fHistEmbPartPtvsJetCorrPt[i] = 0;
  69     fHistJetPtvsJetCorrPt[i] = 0;
  70     fHistDistLeadPart2JetAxis[i] = 0;
  71     fHistEmbBkgArea[i] = 0;
  72     fHistRhoVSEmbBkg[i] = 0;
  73     fHistDeltaPtEmbArea[i] = 0;
  74     fHistDeltaPtEmbvsEP[i] = 0;
  75     fHistDeltaPtEmbPtProbe[i] = 0;
  76   }
  77
  78   SetMakeGeneralHistograms(kTRUE);
  79 }
  80
  81 //________________________________________________________________________
  82 AliAnalysisTaskDeltaPtJEmb::AliAnalysisTaskDeltaPtJEmb(const char *name) :
  83   AliAnalysisTaskEmcalJet(name, kTRUE),
  84   fJetsCont(0),
  85   fTracksCont(0),
  86   fCaloClustersCont(0),
  87   fMinFractionShared(0.5),
  88   fHistEmbJetsPtArea(0),
  89   fHistEmbJetsCorrPtArea(0),
  90   fHistEmbPartPtvsJetPt(0),
  91   fHistEmbPartPtvsJetCorrPt(0),
  92   fHistJetPtvsJetCorrPt(0),
  93   fHistDistLeadPart2JetAxis(0),
  94   fHistEmbBkgArea(0),
  95   fHistRhoVSEmbBkg(0),
  96   fHistDeltaPtEmbArea(0),
  97   fHistDeltaPtEmbvsEP(0),
  98   fHistDeltaPtEmbPtProbe(0),
  99   fHistEmbJetsPhiEta(0),
 100   fHistLeadPartPhiEta(0)
 101 {
 102   // Standard constructor.
 103
 104   fHistEmbJetsPtArea = new TH3*[fNcentBins];
 105   fHistEmbJetsCorrPtArea = new TH3*[fNcentBins];
 106   fHistEmbPartPtvsJetPt = new TH2*[fNcentBins];
 107   fHistEmbPartPtvsJetCorrPt = new TH2*[fNcentBins];
 108   fHistJetPtvsJetCorrPt = new TH2*[fNcentBins];
 109   fHistDistLeadPart2JetAxis = new TH1*[fNcentBins];
 110   fHistEmbBkgArea = new TH2*[fNcentBins];
 111   fHistRhoVSEmbBkg = new TH2*[fNcentBins];
 112   fHistDeltaPtEmbArea = new TH2*[fNcentBins];
 113   fHistDeltaPtEmbvsEP = new TH2*[fNcentBins];
 114   fHistDeltaPtEmbPtProbe = new TH2*[fNcentBins];
 115
 116   for (Int_t i = 0; i < fNcentBins; i++) {
 117     fHistEmbJetsPtArea[i] = 0;
 118     fHistEmbJetsCorrPtArea[i] = 0;
 119     fHistEmbPartPtvsJetPt[i] = 0;
 120     fHistEmbPartPtvsJetCorrPt[i] = 0;
 121     fHistJetPtvsJetCorrPt[i] = 0;
 122     fHistDistLeadPart2JetAxis[i] = 0;
 123     fHistEmbBkgArea[i] = 0;
 124     fHistRhoVSEmbBkg[i] = 0;
 125     fHistDeltaPtEmbArea[i] = 0;
 126     fHistDeltaPtEmbvsEP[i] = 0;
 127     fHistDeltaPtEmbPtProbe[i] = 0;
 128   }
 129
 130   SetMakeGeneralHistograms(kTRUE);
 131 }
 132
 133 //________________________________________________________________________
 134 void AliAnalysisTaskDeltaPtJEmb::UserCreateOutputObjects()
 135 {
 136   // Create user output.
 137
 138   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 139
 140   fJetsCont = GetJetContainer("Jets");
 141   fTracksCont = GetParticleContainer("Tracks");
 142   fCaloClustersCont = GetClusterContainer("CaloClusters");
 143
 144   fHistEmbJetsPhiEta = new TH2F("fHistEmbJetsPhiEta","fHistEmbJetsPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
 145   fHistEmbJetsPhiEta->GetXaxis()->SetTitle("#eta");
 146   fHistEmbJetsPhiEta->GetYaxis()->SetTitle("#phi");
 147   fOutput->Add(fHistEmbJetsPhiEta);
 148
 149   fHistLeadPartPhiEta = new TH2F("fHistLeadPartPhiEta","fHistLeadPartPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
 150   fHistLeadPartPhiEta->GetXaxis()->SetTitle("#eta");
 151   fHistLeadPartPhiEta->GetYaxis()->SetTitle("#phi");
 152   fOutput->Add(fHistLeadPartPhiEta);
 153
 154   TString histname;
 155
 156   const Int_t nbinsZ = 12;
 157   Double_t binsZ[nbinsZ+1] = {0,1,2,3,4,5,6,7,8,9,10,20,1000};
 158
 159   Double_t *binsPt       = GenerateFixedBinArray(fNbins, fMinBinPt, fMaxBinPt);
 160   Double_t *binsCorrPt   = GenerateFixedBinArray(fNbins*2, -fMaxBinPt, fMaxBinPt);
 161   Double_t *binsArea     = GenerateFixedBinArray(50, 0, 2);
 162
 163   for (Int_t i = 0; i < fNcentBins; i++) {
 164     histname = "fHistEmbJetsPtArea_";
 165     histname += i;
 166     fHistEmbJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins, binsPt, nbinsZ, binsZ);
 167     fHistEmbJetsPtArea[i]->GetXaxis()->SetTitle("area");
 168     fHistEmbJetsPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,raw} (GeV/#it{c})");
 169     fOutput->Add(fHistEmbJetsPtArea[i]);
 170
 171     histname = "fHistEmbJetsCorrPtArea_";
 172     histname += i;
 173     fHistEmbJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins * 2, binsCorrPt, nbinsZ, binsZ);
 174     fHistEmbJetsCorrPtArea[i]->GetXaxis()->SetTitle("area");
 175     fHistEmbJetsCorrPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,corr} (GeV/#it{c})");
 176     fOutput->Add(fHistEmbJetsCorrPtArea[i]);
 177
 178     histname = "fHistEmbPartPtvsJetPt_";
 179     histname += i;
 180     fHistEmbPartPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 181     fHistEmbPartPtvsJetPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
 182     fHistEmbPartPtvsJetPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
 183     fHistEmbPartPtvsJetPt[i]->GetZaxis()->SetTitle("counts");
 184     fOutput->Add(fHistEmbPartPtvsJetPt[i]);
 185
 186     histname = "fHistEmbPartPtvsJetCorrPt_";
 187     histname += i;
 188     fHistEmbPartPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
 189                                             fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
 190     fHistEmbPartPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
 191     fHistEmbPartPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
 192     fHistEmbPartPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
 193     fOutput->Add(fHistEmbPartPtvsJetCorrPt[i]);
 194
 195     histname = "fHistJetPtvsJetCorrPt_";
 196     histname += i;
 197     fHistJetPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
 198                                         fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
 199     fHistJetPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
 200     fHistJetPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
 201     fHistJetPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
 202     fOutput->Add(fHistJetPtvsJetCorrPt[i]);
 203
 204     histname = "fHistDistLeadPart2JetAxis_";
 205     histname += i;
 206     fHistDistLeadPart2JetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
 207     fHistDistLeadPart2JetAxis[i]->GetXaxis()->SetTitle("distance");
 208     fHistDistLeadPart2JetAxis[i]->GetYaxis()->SetTitle("counts");
 209     fOutput->Add(fHistDistLeadPart2JetAxis[i]);
 210
 211     histname = "fHistEmbBkgArea_";
 212     histname += i;
 213     fHistEmbBkgArea[i] = new TH2F(histname.Data(), histname.Data(), 50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
 214     fHistEmbBkgArea[i]->GetXaxis()->SetTitle("area");
 215     fHistEmbBkgArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
 216     fOutput->Add(fHistEmbBkgArea[i]);
 217
 218     histname = "fHistRhoVSEmbBkg_";
 219     histname += i;
 220     fHistRhoVSEmbBkg[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 221     fHistRhoVSEmbBkg[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
 222     fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
 223     fOutput->Add(fHistRhoVSEmbBkg[i]);
 224
 225     histname = "fHistDeltaPtEmbArea_";
 226     histname += i;
 227     fHistDeltaPtEmbArea[i] = new TH2F(histname.Data(), histname.Data(),
 228                                       50, 0, 2, fNbins * 2, -fMaxBinPt, fMaxBinPt);
 229     fHistDeltaPtEmbArea[i]->GetXaxis()->SetTitle("area");
 230     fHistDeltaPtEmbArea[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
 231     fHistDeltaPtEmbArea[i]->GetZaxis()->SetTitle("counts");
 232     fOutput->Add(fHistDeltaPtEmbArea[i]);
 233
 234     histname = "fHistDeltaPtEmbvsEP_";
 235     histname += i;
 236     fHistDeltaPtEmbvsEP[i] = new TH2F(histname.Data(), histname.Data(),
 237                                       402, -TMath::Pi()*1.01, TMath::Pi()*3.01, fNbins * 2, -fMaxBinPt, fMaxBinPt);
 238     fHistDeltaPtEmbvsEP[i]->GetXaxis()->SetTitle("#phi_{jet} - #Psi_{EP}");
 239     fHistDeltaPtEmbvsEP[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
 240     fHistDeltaPtEmbvsEP[i]->GetZaxis()->SetTitle("counts");
 241     fOutput->Add(fHistDeltaPtEmbvsEP[i]);
 242
 243     histname = "fHistDeltaPtEmbPtProbe_";
 244     histname += i;
 245     fHistDeltaPtEmbPtProbe[i] = new TH2F(histname.Data(), histname.Data(),
 246                                          fNbins, fMinBinPt, fMaxBinPt, fNbins * 2, -fMaxBinPt, fMaxBinPt);
 247     fHistDeltaPtEmbPtProbe[i]->GetXaxis()->SetTitle("#it{p}_{T,probe} (GeV/#it{c})");
 248     fHistDeltaPtEmbPtProbe[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
 249     fHistDeltaPtEmbPtProbe[i]->GetZaxis()->SetTitle("counts");
 250     fOutput->Add(fHistDeltaPtEmbPtProbe[i]);
 251   }
 252
 253   delete[] binsPt;
 254   delete[] binsCorrPt;
 255   delete[] binsArea;
 256
 257   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 258 }
 259
 260 //________________________________________________________________________
 261 Bool_t AliAnalysisTaskDeltaPtJEmb::FillHistograms()
 262 {
 263   // Fill histograms.
 264
 265   // ************
 266   // Embedding
 267   // Jets should already be matched, GetClosestJet gives matched partner. Cut on shared fraction to be applied
 268   // _________________________________
 269
 270   if (fJetsCont) {
 271     AliEmcalJet *jet1 = NULL;
 272     AliEmcalJet *jet2 = NULL;
 273     fJetsCont->ResetCurrentID();
 274     while((jet1 = fJetsCont->GetNextAcceptJet())) {
 275
 276       jet2 = jet1->ClosestJet();
 277       if(!jet2) continue;
 278
 279       Double_t fraction = fJetsCont->GetFractionSharedPt(jet1);
 280       if(fMinFractionShared>0. && fraction<fMinFractionShared) continue;
 281
 282       TLorentzVector mom;
 283       fJetsCont->GetLeadingHadronMomentum(mom,jet1);
 284
 285       Double_t distLeading2Jet = TMath::Sqrt((jet1->Eta() - mom.Eta()) * (jet1->Eta() - mom.Eta()) + (jet1->Phi() - mom.Phi()+TMath::Pi()) * (jet1->Phi() - mom.Phi()+TMath::Pi()));
 286
 287       fHistEmbPartPtvsJetPt[fCentBin]->Fill(jet2->Pt(), jet1->Pt());
 288       fHistEmbPartPtvsJetCorrPt[fCentBin]->Fill(jet2->Pt(), jet1->Pt() - jet1->Area() * fRhoVal);
 289       fHistLeadPartPhiEta->Fill(mom.Eta(), mom.Phi()+TMath::Pi());
 290       fHistDistLeadPart2JetAxis[fCentBin]->Fill(distLeading2Jet);
 291
 292       fHistEmbJetsPtArea[fCentBin]->Fill(jet1->Area(), jet1->Pt(), mom.Pt());
 293       fHistEmbJetsCorrPtArea[fCentBin]->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area(), mom.Pt());
 294       fHistEmbJetsPhiEta->Fill(jet1->Eta(), jet1->Phi());
 295       fHistJetPtvsJetCorrPt[fCentBin]->Fill(jet1->Pt(), jet1->Pt() - fRhoVal * jet1->Area());
 296
 297       fHistEmbBkgArea[fCentBin]->Fill(jet1->Area(), jet1->Pt() - jet2->Pt());
 298       fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal * jet1->Area(), jet1->Pt() - jet2->Pt());
 299       fHistDeltaPtEmbArea[fCentBin]->Fill(jet1->Area(), jet1->Pt() - jet1->Area() * fRhoVal - jet2->Pt());
 300       fHistDeltaPtEmbvsEP[fCentBin]->Fill(jet1->Phi() - fEPV0, jet1->Pt() - jet1->Area() * fRhoVal - jet2->Pt());
 301       fHistDeltaPtEmbPtProbe[fCentBin]->Fill(jet2->Pt(), jet1->Pt() - jet1->Area() * fRhoVal - jet2->Pt());
 302     }
 303   }
 304
 305   return kTRUE;
 306 }
 307
 308 //________________________________________________________________________
 309 void AliAnalysisTaskDeltaPtJEmb::ExecOnce()
 310 {
 311   // Initialize the analysis.
 312
 313   AliAnalysisTaskEmcalJet::ExecOnce();
 314
 315   if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
 316   if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
 317   if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
 318 }