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