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