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