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