PWGJE/EMCALJetTasks/AliAnalysisTaskRhoBase.cxx

   1 // $Id$
   2 //
   3 // Base class for rho calculation.
   4 // Calculates parameterized rho for given centrality independent of input.
   5 //
   6 // Author: S.Aiola
   7
   8 #include <TF1.h>
   9 #include <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TClonesArray.h>
  12
  13 #include "AliLog.h"
  14 #include "AliRhoParameter.h"
  15 #include "AliEmcalJet.h"
  16
  17 #include "AliAnalysisTaskRhoBase.h"
  18
  19 ClassImp(AliAnalysisTaskRhoBase)
  20
  21 //________________________________________________________________________
  22 AliAnalysisTaskRhoBase::AliAnalysisTaskRhoBase() :
  23   AliAnalysisTaskEmcalJetDev("AliAnalysisTaskRhoBase", kFALSE),
  24   fOutRhoName(),
  25   fOutRhoScaledName(),
  26   fCompareRhoName(),
  27   fCompareRhoScaledName(),
  28   fRhoFunction(0),
  29   fScaleFunction(0),
  30   fInEventSigmaRho(35.83),
  31   fAttachToEvent(kTRUE),
  32   fOutRho(0),
  33   fOutRhoScaled(0),
  34   fCompareRho(0),
  35   fCompareRhoScaled(0),
  36   fHistJetPtvsCent(0),
  37   fHistJetAreavsCent(0),
  38   fHistJetRhovsCent(0),
  39   fHistNjetvsCent(0),
  40   fHistJetPtvsNtrack(0),
  41   fHistJetAreavsNtrack(0),
  42   fHistNjetvsNtrack(0),
  43   fHistRhovsCent(0),
  44   fHistRhoScaledvsCent(0),
  45   fHistDeltaRhovsCent(0),
  46   fHistDeltaRhoScalevsCent(0),
  47   fHistRhovsNtrack(0),
  48   fHistRhoScaledvsNtrack(0),
  49   fHistDeltaRhovsNtrack(0),
  50   fHistDeltaRhoScalevsNtrack(0),
  51   fHistRhovsNcluster(0),
  52   fHistRhoScaledvsNcluster(0)
  53 {
  54   // Constructor.
  55
  56   for (Int_t i = 0; i < 4; i++) {
  57     fHistJetNconstVsPt[i] = 0;
  58     fHistJetRhovsEta[i] = 0;
  59   }
  60   for (Int_t i = 0; i < 12; i++) {
  61     fHistNjUEoverNjVsNj[i] = 0;
  62   }
  63 }
  64
  65 //________________________________________________________________________
  66 AliAnalysisTaskRhoBase::AliAnalysisTaskRhoBase(const char *name, Bool_t histo) :
  67   AliAnalysisTaskEmcalJetDev(name, histo),
  68   fOutRhoName(),
  69   fOutRhoScaledName(),
  70   fCompareRhoName(),
  71   fCompareRhoScaledName(),
  72   fRhoFunction(0),
  73   fScaleFunction(0),
  74   fInEventSigmaRho(35.83),
  75   fAttachToEvent(kTRUE),
  76   fOutRho(0),
  77   fOutRhoScaled(0),
  78   fCompareRho(0),
  79   fCompareRhoScaled(0),
  80   fHistJetPtvsCent(0),
  81   fHistJetAreavsCent(0),
  82   fHistJetRhovsCent(0),
  83   fHistNjetvsCent(0),
  84   fHistJetPtvsNtrack(0),
  85   fHistJetAreavsNtrack(0),
  86   fHistNjetvsNtrack(0),
  87   fHistRhovsCent(0),
  88   fHistRhoScaledvsCent(0),
  89   fHistDeltaRhovsCent(0),
  90   fHistDeltaRhoScalevsCent(0),
  91   fHistRhovsNtrack(0),
  92   fHistRhoScaledvsNtrack(0),
  93   fHistDeltaRhovsNtrack(0),
  94   fHistDeltaRhoScalevsNtrack(0),
  95   fHistRhovsNcluster(0),
  96   fHistRhoScaledvsNcluster(0)
  97 {
  98   // Constructor.
  99
 100   for (Int_t i = 0; i < 4; i++) {
 101     fHistJetNconstVsPt[i] = 0;
 102     fHistJetRhovsEta[i] = 0;
 103   }
 104   for (Int_t i = 0; i < 12; i++) {
 105     fHistNjUEoverNjVsNj[i] = 0;
 106   }
 107   SetMakeGeneralHistograms(histo);
 108 }
 109
 110 //________________________________________________________________________
 111 void AliAnalysisTaskRhoBase::UserCreateOutputObjects()
 112 {
 113   // User create output objects, called at the beginning of the analysis.
 114
 115   if (!fCreateHisto)
 116     return;
 117
 118   AliAnalysisTaskEmcalJetDev::UserCreateOutputObjects();
 119
 120   fHistRhovsCent = new TH2F("fHistRhovsCent", "fHistRhovsCent", 101, -1,  100, fNbins, fMinBinPt, fMaxBinPt*2);
 121   fHistRhovsCent->GetXaxis()->SetTitle("Centrality (%)");
 122   fHistRhovsCent->GetYaxis()->SetTitle("#rho (GeV/c * rad^{-1})");
 123   fOutput->Add(fHistRhovsCent);
 124
 125   if (fParticleCollArray.GetEntriesFast()>0) {
 126     fHistRhovsNtrack = new TH2F("fHistRhovsNtrack", "fHistRhovsNtrack", 150, 0, 6000, fNbins, fMinBinPt, fMaxBinPt*2);
 127     fHistRhovsNtrack->GetXaxis()->SetTitle("No. of tracks");
 128     fHistRhovsNtrack->GetYaxis()->SetTitle("#rho (GeV/c * rad^{-1})");
 129     fOutput->Add(fHistRhovsNtrack);
 130   }
 131
 132   if (fClusterCollArray.GetEntriesFast()>0) {
 133     fHistRhovsNcluster = new TH2F("fHistRhovsNcluster", "fHistRhovsNcluster", 50, 0, 1500, fNbins, fMinBinPt, fMaxBinPt*2);
 134     fHistRhovsNcluster->GetXaxis()->SetTitle("No. of tracks");
 135     fHistRhovsNcluster->GetYaxis()->SetTitle("#rho (GeV/c * rad^{-1})");
 136     fOutput->Add(fHistRhovsNcluster);
 137   }
 138
 139   if (fJetCollArray.GetEntriesFast()>0) {
 140     fHistJetPtvsCent = new TH2F("fHistJetPtvsCent", "fHistJetPtvsCent", 101, -1,  100, fNbins, fMinBinPt, fMaxBinPt);
 141     fHistJetPtvsCent->GetXaxis()->SetTitle("Centrality (%)");
 142     fHistJetPtvsCent->GetYaxis()->SetTitle("#it{p}_{T,jet} (GeV/c)");
 143     fOutput->Add(fHistJetPtvsCent);
 144
 145     fHistJetAreavsCent = new TH2F("fHistJetAreavsCent", "fHistJetAreavsCent", 101, -1, 100, 100, 0, 1);
 146     fHistJetAreavsCent->GetXaxis()->SetTitle("Centrality (%)");
 147     fHistJetAreavsCent->GetYaxis()->SetTitle("Jet area");
 148     fOutput->Add(fHistJetAreavsCent);
 149
 150     fHistJetRhovsCent = new TH2F("fHistJetRhovsCent", "fHistJetRhovsCent", 101, -1, 100, fNbins, fMinBinPt, fMaxBinPt*2);
 151     fHistJetRhovsCent->GetXaxis()->SetTitle("Centrality (%)");
 152     fHistJetRhovsCent->GetYaxis()->SetTitle("Jet #rho (GeV/c)");
 153     fOutput->Add(fHistJetRhovsCent);
 154
 155     fHistNjetvsCent = new TH2F("fHistNjetvsCent",  "fHistNjetvsCent", 101, -1, 100, 150, -0.5, 149.5);
 156     fHistNjetvsCent->GetXaxis()->SetTitle("Centrality (%)");
 157     fHistNjetvsCent->GetYaxis()->SetTitle("No. of jets");
 158     fOutput->Add(fHistNjetvsCent);
 159
 160
 161     if (fParticleCollArray.GetEntriesFast()>0) {
 162       fHistJetPtvsNtrack = new TH2F("fHistJetPtvsNtrack", "fHistJetPtvsNtrack", 150, 0, 6000, fNbins, fMinBinPt, fMaxBinPt);
 163       fHistJetPtvsNtrack->GetXaxis()->SetTitle("No. of tracks");
 164       fHistJetPtvsNtrack->GetYaxis()->SetTitle("#it{p}_{T,jet} (GeV/c)");
 165       fOutput->Add(fHistJetPtvsNtrack);
 166
 167       fHistJetAreavsNtrack = new TH2F("fHistJetAreavsNtrack", "fHistJetAreavsNtrack", 150, 0, 6000, 100, 0, 1);
 168       fHistJetAreavsNtrack->GetXaxis()->SetTitle("No. of tracks");
 169       fHistJetAreavsNtrack->GetYaxis()->SetTitle("Jet area");
 170       fOutput->Add(fHistJetAreavsNtrack);
 171
 172       fHistNjetvsNtrack = new TH2F("fHistNjetvsNtrack", "fHistNjetvsNtrack", 150, 0, 6000, 150, -0.5, 149.5);
 173       fHistJetAreavsNtrack->GetXaxis()->SetTitle("No. of tracks");
 174       fHistJetAreavsNtrack->GetYaxis()->SetTitle("Jet area");
 175       fOutput->Add(fHistNjetvsNtrack);
 176     }
 177
 178
 179     TString name;
 180     for (Int_t i = 0; i < 4; i++) {
 181       name = Form("fHistJetNconstVsPt_%d",i);
 182       fHistJetNconstVsPt[i] = new TH2F(name, name, 150, -0.5, 149.5, fNbins, fMinBinPt, fMaxBinPt);
 183       fHistJetNconstVsPt[i]->GetXaxis()->SetTitle("No. of constituents");
 184       fHistJetNconstVsPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
 185       fOutput->Add(fHistJetNconstVsPt[i]);
 186
 187       name = Form("fHistJetRhovsEta_%d",i);
 188       fHistJetRhovsEta[i] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt*2, 16, -0.8, 0.8);
 189       fHistJetRhovsEta[i]->GetXaxis()->SetTitle("#rho (GeV/c * rad^{-1})");
 190       fHistJetRhovsEta[i]->GetYaxis()->SetTitle("#eta");
 191       fOutput->Add(fHistJetRhovsEta[i]);
 192
 193       for (Int_t j = 0; j < 3; j++) {
 194         name = Form("NjUEoverNjVsNj_%d_%d",i,j+1);
 195         fHistNjUEoverNjVsNj[i*3+j] = new TH2F(name, name, 150, -0.5, 149.5, 120, 0.01, 1.21);
 196         fHistNjUEoverNjVsNj[i*3+j]->GetXaxis()->SetTitle("N_{jet}");
 197         fHistNjUEoverNjVsNj[i*3+j]->GetYaxis()->SetTitle("N_{jet_{UE}} / N_{jet}");
 198         fOutput->Add(fHistNjUEoverNjVsNj[i*3+j]);
 199       }
 200     }
 201   }
 202
 203   if (!fCompareRhoName.IsNull()) {
 204     fHistDeltaRhovsCent = new TH2F("fHistDeltaRhovsCent", "fHistDeltaRhovsCent", 101, -1, 100, fNbins, -fMaxBinPt, fMaxBinPt);
 205     fHistDeltaRhovsCent->GetXaxis()->SetTitle("Centrality (%)");
 206     fHistDeltaRhovsCent->GetYaxis()->SetTitle("#Delta#rho (GeV/c * rad^{-1})");
 207     fOutput->Add(fHistDeltaRhovsCent);
 208
 209     if (fParticleCollArray.GetEntriesFast()>0) {
 210       fHistDeltaRhovsNtrack = new TH2F("fHistDeltaRhovsNtrack", "fHistDeltaRhovsNtrack", 150, 0, 6000, fNbins, -fMaxBinPt, fMaxBinPt);
 211       fHistDeltaRhovsNtrack->GetXaxis()->SetTitle("No. of tracks");
 212       fHistDeltaRhovsNtrack->GetYaxis()->SetTitle("#Delta#rho (GeV/c * rad^{-1})");
 213       fOutput->Add(fHistDeltaRhovsNtrack);
 214     }
 215   }
 216
 217   if (fScaleFunction) {
 218     fHistRhoScaledvsCent = new TH2F("fHistRhoScaledvsCent", "fHistRhoScaledvsCent", 101, -1, 100, fNbins, fMinBinPt , fMaxBinPt*2);
 219     fHistRhoScaledvsCent->GetXaxis()->SetTitle("Centrality (%)");
 220     fHistRhoScaledvsCent->GetYaxis()->SetTitle("#rho_{scaled} (GeV/c * rad^{-1})");
 221     fOutput->Add(fHistRhoScaledvsCent);
 222
 223     if (fParticleCollArray.GetEntriesFast()>0) {
 224       fHistRhoScaledvsNtrack = new TH2F("fHistRhoScaledvsNtrack", "fHistRhoScaledvsNtrack", 150, 0, 6000, fNbins, fMinBinPt, fMaxBinPt*2);
 225       fHistRhoScaledvsNtrack->GetXaxis()->SetTitle("No. of tracks");
 226       fHistRhoScaledvsNtrack->GetYaxis()->SetTitle("#rho_{scaled} (GeV/c * rad^{-1})");
 227       fOutput->Add(fHistRhoScaledvsNtrack);
 228     }
 229
 230     if (fClusterCollArray.GetEntriesFast()>0) {
 231       fHistRhoScaledvsNcluster = new TH2F("fHistRhoScaledvsNcluster", "fHistRhoScaledvsNcluster", 50, 0, 1500, fNbins, fMinBinPt, fMaxBinPt*2);
 232       fHistRhoScaledvsNcluster->GetXaxis()->SetTitle("No. of clusters");
 233       fHistRhoScaledvsNcluster->GetYaxis()->SetTitle("#rho_{scaled} (GeV/c * rad^{-1})");
 234       fOutput->Add(fHistRhoScaledvsNcluster);
 235     }
 236
 237     if (!fCompareRhoScaledName.IsNull()) {
 238       fHistDeltaRhoScalevsCent = new TH2F("fHistDeltaRhoScalevsCent", "fHistDeltaRhoScalevsCent", 101, -1, 100, fNbins, -fMaxBinPt, fMaxBinPt);
 239       fHistDeltaRhoScalevsCent->GetXaxis()->SetTitle("Centrality (%)");
 240       fHistDeltaRhoScalevsCent->GetYaxis()->SetTitle("#Delta#rho_{scaled} (GeV/c * rad^{-1})");
 241       fOutput->Add(fHistDeltaRhoScalevsCent);
 242
 243       if (fParticleCollArray.GetEntriesFast()>0) {
 244         fHistDeltaRhoScalevsNtrack = new TH2F("fHistDeltaRhoScalevsNtrack", "fHistDeltaRhoScalevsNtrack", 150, 0, 6000, fNbins, -fMaxBinPt, fMaxBinPt);
 245         fHistDeltaRhoScalevsNtrack->GetXaxis()->SetTitle("No. of tracks");
 246         fHistDeltaRhoScalevsNtrack->GetYaxis()->SetTitle("#Delta#rho_{scaled} (GeV/c * rad^{-1})");
 247         fOutput->Add(fHistDeltaRhoScalevsNtrack);
 248       }
 249     }
 250   }
 251 }
 252
 253 //________________________________________________________________________
 254 Bool_t AliAnalysisTaskRhoBase::Run()
 255 {
 256   // Run the analysis.
 257
 258   Double_t rho = GetRhoFactor(fCent);
 259   fOutRho->SetVal(rho);
 260
 261   if (fScaleFunction) {
 262     Double_t rhoScaled = rho * GetScaleFactor(fCent);
 263     fOutRhoScaled->SetVal(rhoScaled);
 264   }
 265
 266   return kTRUE;
 267 }
 268
 269 //________________________________________________________________________
 270 Bool_t AliAnalysisTaskRhoBase::FillHistograms()
 271 {
 272   // Fill histograms.
 273
 274   Int_t Ntracks   = 0;
 275   Int_t Nclusters = 0;
 276
 277   if (fTracks)
 278     Ntracks = fTracks->GetEntries();
 279   if (fCaloClusters)
 280     Nclusters = fCaloClusters->GetEntries();
 281
 282   if (fJets) {
 283     Int_t    Njets         = fJets->GetEntries();
 284     Int_t    NjetAcc       = 0;
 285     Int_t    NjetUE1Sigma  = 0;
 286     Int_t    NjetUE2Sigma  = 0;
 287     Int_t    NjetUE3Sigma  = 0;
 288     Double_t rhoPlus1Sigma = fOutRho->GetVal() + fInEventSigmaRho;
 289     Double_t rhoPlus2Sigma = fOutRho->GetVal() + 2*fInEventSigmaRho;
 290     Double_t rhoPlus3Sigma = fOutRho->GetVal() + 3*fInEventSigmaRho;
 291
 292     for (Int_t i = 0; i < Njets; ++i) {
 293
 294       AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(i));
 295       if (!jet) {
 296         AliError(Form("%s: Could not receive jet %d", GetName(), i));
 297         continue;
 298       }
 299
 300       if (!AcceptJet(jet))
 301         continue;
 302
 303       fHistJetPtvsCent->Fill(fCent, jet->Pt());
 304       fHistJetAreavsCent->Fill(fCent, jet->Area());
 305       fHistJetRhovsCent->Fill(fCent, jet->Pt() / jet->Area());
 306       fHistJetRhovsEta[fCentBin]->Fill(jet->Pt() / jet->Area(), jet->Eta());
 307
 308       if (fTracks) {
 309         fHistJetPtvsNtrack->Fill(Ntracks, jet->Pt());
 310         fHistJetAreavsNtrack->Fill(Ntracks, jet->Area());
 311       }
 312
 313       fHistJetNconstVsPt[fCentBin]->Fill(jet->GetNumberOfConstituents(), jet->Pt());
 314
 315       if (jet->Pt() < rhoPlus1Sigma * jet->Area())
 316         NjetUE1Sigma++;
 317
 318       if (jet->Pt() < rhoPlus2Sigma * jet->Area())
 319         NjetUE2Sigma++;
 320
 321       if (jet->Pt() < rhoPlus3Sigma * jet->Area())
 322         NjetUE3Sigma++;
 323
 324       NjetAcc++;
 325     }
 326
 327     if (NjetAcc>0) {
 328       fHistNjUEoverNjVsNj[fCentBin*3  ]->Fill(NjetAcc,1.*NjetUE1Sigma/NjetAcc);
 329       fHistNjUEoverNjVsNj[fCentBin*3+1]->Fill(NjetAcc,1.*NjetUE2Sigma/NjetAcc);
 330       fHistNjUEoverNjVsNj[fCentBin*3+2]->Fill(NjetAcc,1.*NjetUE3Sigma/NjetAcc);
 331     }
 332
 333     fHistNjetvsCent->Fill(fCent, NjetAcc);
 334     if (fTracks)
 335       fHistNjetvsNtrack->Fill(Ntracks, NjetAcc);
 336   }
 337
 338   fHistRhovsCent->Fill(fCent, fOutRho->GetVal());
 339
 340   if (fTracks)
 341     fHistRhovsNtrack->Fill(Ntracks, fOutRho->GetVal());
 342   if (fCaloClusters)
 343     fHistRhovsNcluster->Fill(Nclusters, fOutRho->GetVal());
 344   if (fCompareRho) {
 345     fHistDeltaRhovsCent->Fill(fCent, fOutRho->GetVal() - fCompareRho->GetVal());
 346     if (fTracks)
 347       fHistDeltaRhovsNtrack->Fill(Ntracks, fOutRho->GetVal() - fCompareRho->GetVal());
 348   }
 349
 350   if (fOutRhoScaled) {
 351     fHistRhoScaledvsCent->Fill(fCent, fOutRhoScaled->GetVal());
 352     if (fTracks)
 353       fHistRhoScaledvsNtrack->Fill(Ntracks, fOutRhoScaled->GetVal());
 354     if (fCaloClusters)
 355       fHistRhoScaledvsNcluster->Fill(Nclusters,  fOutRhoScaled->GetVal());
 356     if (fCompareRhoScaled) {
 357       fHistDeltaRhoScalevsCent->Fill(fCent, fOutRhoScaled->GetVal() - fCompareRhoScaled->GetVal());
 358       if (fTracks)
 359         fHistDeltaRhoScalevsNtrack->Fill(Ntracks, fOutRhoScaled->GetVal() - fCompareRhoScaled->GetVal());
 360     }
 361   }
 362
 363   return kTRUE;
 364 }
 365
 366
 367 //________________________________________________________________________
 368 void AliAnalysisTaskRhoBase::ExecOnce()
 369 {
 370   // Init the analysis.
 371
 372   if (!fOutRho) {
 373     fOutRho = new AliRhoParameter(fOutRhoName, 0);
 374
 375     if (fAttachToEvent) {
 376       if (!(InputEvent()->FindListObject(fOutRhoName))) {
 377         InputEvent()->AddObject(fOutRho);
 378       } else {
 379         AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fOutRhoName.Data()));
 380         return;
 381       }
 382     }
 383   }
 384
 385   if (fScaleFunction && !fOutRhoScaled) {
 386     fOutRhoScaled = new AliRhoParameter(fOutRhoScaledName, 0);
 387
 388     if (fAttachToEvent) {
 389       if (!(InputEvent()->FindListObject(fOutRhoScaledName))) {
 390         InputEvent()->AddObject(fOutRhoScaled);
 391       } else {
 392         AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fOutRhoScaledName.Data()));
 393         return;
 394       }
 395     }
 396   }
 397
 398   if (!fCompareRhoName.IsNull() && !fCompareRho) {
 399     fCompareRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fCompareRhoName));
 400     if (!fCompareRho) {
 401       AliWarning(Form("%s: Could not retrieve rho %s!", GetName(), fCompareRhoName.Data()));
 402     }
 403   }
 404
 405   if (!fCompareRhoScaledName.IsNull() && !fCompareRhoScaled) {
 406     fCompareRhoScaled = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fCompareRhoScaledName));
 407     if (!fCompareRhoScaled) {
 408       AliWarning(Form("%s: Could not retrieve rho %s!", GetName(), fCompareRhoScaledName.Data()));
 409     }
 410   }
 411
 412   AliAnalysisTaskEmcalJetDev::ExecOnce();
 413 }
 414
 415 //________________________________________________________________________
 416 Double_t AliAnalysisTaskRhoBase::GetRhoFactor(Double_t cent)
 417 {
 418   // Return rho per centrality.
 419
 420   Double_t rho = 0;
 421   if (fRhoFunction)
 422     rho = fRhoFunction->Eval(cent);
 423   return rho;
 424 }
 425
 426 //________________________________________________________________________
 427 Double_t AliAnalysisTaskRhoBase::GetScaleFactor(Double_t cent)
 428 {
 429   // Get scale factor.
 430
 431   Double_t scale = 1;
 432   if (fScaleFunction)
 433     scale = fScaleFunction->Eval(cent);
 434   return scale;
 435 }