[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliAnalysisTaskRhoBase.cxx

// $Id$
//
// Base class for rho calculation.
// Calculates parameterized rho for given centrality independent of input.
//
// Author: S.Aiola

#include <TF1.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TClonesArray.h>

#include "AliLog.h"
#include "AliRhoParameter.h"
#include "AliEmcalJet.h"

#include "AliAnalysisTaskRhoBase.h"

ClassImp(AliAnalysisTaskRhoBase)

//________________________________________________________________________
AliAnalysisTaskRhoBase::AliAnalysisTaskRhoBase() : 
  AliAnalysisTaskEmcalJet("AliAnalysisTaskRhoBase", kFALSE),
  fRhoScaledName(),
  fCompareRhoName(),
  fCompareRhoScaledName(),
  fRhoFunction(0),
  fScaleFunction(0),
  fInEventSigmaRho(35.83),
  fAttachToEvent(kTRUE),
  fRhoScaled(0),
  fCompareRho(0),
  fCompareRhoScaled(0),
  fHistJetPtvsCent(0),
  fHistJetAreavsCent(0),
  fHistJetRhovsCent(0),
  fHistNjetvsCent(0),
  fHistJetPtvsNtrack(0),
  fHistJetAreavsNtrack(0),
  fHistNjetvsNtrack(0),
  fHistRhovsCent(0),
  fHistRhoScaledvsCent(0),
  fHistDeltaRhovsCent(0),
  fHistDeltaRhoScalevsCent(0),
  fHistRhovsNtrack(0),
  fHistRhoScaledvsNtrack(0),
  fHistDeltaRhovsNtrack(0),
  fHistDeltaRhoScalevsNtrack(0),
  fHistRhovsNcluster(0),
  fHistRhoScaledvsNcluster(0)
{
  // Constructor.

  for (Int_t i = 0; i < 4; i++) {
    fHistJetNconstVsPt[i] = 0;
    fHistJetRhovsEta[i] = 0;
  }
  for (Int_t i = 0; i < 12; i++) {
    fHistNjUEoverNjVsNj[i] = 0;
  }
}

//________________________________________________________________________
AliAnalysisTaskRhoBase::AliAnalysisTaskRhoBase(const char *name, Bool_t histo) :
  AliAnalysisTaskEmcalJet(name, histo),
  fRhoScaledName(),
  fCompareRhoName(),
  fCompareRhoScaledName(),
  fRhoFunction(0),
  fScaleFunction(0),
  fInEventSigmaRho(35.83),
  fAttachToEvent(kTRUE),
  fRhoScaled(0),
  fCompareRho(0),
  fCompareRhoScaled(0),
  fHistJetPtvsCent(0),
  fHistJetAreavsCent(0),
  fHistJetRhovsCent(0),
  fHistNjetvsCent(0),
  fHistJetPtvsNtrack(0),
  fHistJetAreavsNtrack(0),
  fHistNjetvsNtrack(0),
  fHistRhovsCent(0),
  fHistRhoScaledvsCent(0),
  fHistDeltaRhovsCent(0),
  fHistDeltaRhoScalevsCent(0),
  fHistRhovsNtrack(0),
  fHistRhoScaledvsNtrack(0),
  fHistDeltaRhovsNtrack(0),
  fHistDeltaRhoScalevsNtrack(0),
  fHistRhovsNcluster(0),
  fHistRhoScaledvsNcluster(0)
{
  // Constructor.

  for (Int_t i = 0; i < 4; i++) {
    fHistJetNconstVsPt[i] = 0;
    fHistJetRhovsEta[i] = 0;
  }
  for (Int_t i = 0; i < 12; i++) {
    fHistNjUEoverNjVsNj[i] = 0;
  }
  SetMakeGeneralHistograms(histo);
}

//________________________________________________________________________
void AliAnalysisTaskRhoBase::UserCreateOutputObjects()
{
  // User create output objects, called at the beginning of the analysis.

  if (!fCreateHisto)
    return;

  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();

  fHistRhovsCent = new TH2F("RhovsCent", "RhovsCent", 101, -1,  100, fNbins, fMinBinPt, fMaxBinPt*2);
  fOutput->Add(fHistRhovsCent);

  if (!fTracksName.IsNull()) {
    fHistRhovsNtrack = new TH2F("RhovsNtrack", "RhovsNtrack", 125, 0, 4000, fNbins, fMinBinPt, fMaxBinPt*2);
    fOutput->Add(fHistRhovsNtrack);
  }

  if (!fCaloName.IsNull()) {
    fHistRhovsNcluster = new TH2F("RhovsNcluster", "RhovsNcluster", 50, 0, 1500, fNbins, fMinBinPt, fMaxBinPt*2);
    fOutput->Add(fHistRhovsNcluster);
  }

  if (!fJetsName.IsNull()) {
    fHistJetPtvsCent            = new TH2F("JetPtvsCent",           "JetPtvsCent",           101, -1,  100,   fNbins, fMinBinPt, fMaxBinPt);
    fHistJetAreavsCent          = new TH2F("JetAreavsCent",         "JetAreavsCent",         101, -1,  100,   30, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
    fHistJetRhovsCent           = new TH2F("fHistJetRhovsCent",     "fHistJetRhovsCent",     101, -1,  100,   fNbins, fMinBinPt, fMaxBinPt*2);
    fHistNjetvsCent             = new TH2F("NjetvsCent",            "NjetvsCent",            101, -1,  100,   150, -0.5, 149.5);

    fOutput->Add(fHistJetPtvsCent);
    fOutput->Add(fHistJetAreavsCent);
    fOutput->Add(fHistNjetvsCent);

    if (!fTracksName.IsNull()) {
      fHistJetPtvsNtrack        = new TH2F("JetPtvsNtrack",         "JetPtvsNtrack",         125,  0,  4000,  fNbins, fMinBinPt, fMaxBinPt);
      fHistJetAreavsNtrack      = new TH2F("JetAreavsNtrack",       "JetAreavsNtrack",       125,  0,  4000,  30, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
      fHistNjetvsNtrack         = new TH2F("NjetvsNtrack",          "rNjetvsNtrack",         125,  0,  4000,  150, -0.5, 149.5);

      fOutput->Add(fHistJetPtvsNtrack);
      fOutput->Add(fHistJetAreavsNtrack);
      fOutput->Add(fHistNjetvsNtrack);
    }


    TString name;
    for (Int_t i = 0; i < 4; i++) {
      name = Form("fHistJetNconstVsPt_%d",i);
      fHistJetNconstVsPt[i] = new TH2F(name, name, 150, -0.5, 149.5, fNbins, fMinBinPt, fMaxBinPt);
      fHistJetNconstVsPt[i]->GetXaxis()->SetTitle("# of constituents");
      fHistJetNconstVsPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
      fOutput->Add(fHistJetNconstVsPt[i]);

      name = Form("HistJetRhovsEta_%d",i);
      fHistJetRhovsEta[i] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt*2, 16, -0.8, 0.8);
      fHistJetRhovsEta[i]->GetXaxis()->SetTitle("Rho");
      fHistJetRhovsEta[i]->GetYaxis()->SetTitle("eta");
      fOutput->Add(fHistJetRhovsEta[i]);

      for (Int_t j = 0; j < 3; j++) {
	name = Form("NjUEoverNjVsNj_%d_%d",i,j+1);
	fHistNjUEoverNjVsNj[i*3+j] = new TH2F(name, name, 150, -0.5, 149.5, 120, 0.01, 1.21);
	fHistNjUEoverNjVsNj[i*3+j]->GetXaxis()->SetTitle("N_{jet}");
	fHistNjUEoverNjVsNj[i*3+j]->GetYaxis()->SetTitle("N_{jet_{UE}} / N_{jet}");
	fOutput->Add(fHistNjUEoverNjVsNj[i*3+j]);
      }
    }
  }
  
  if (!fCompareRhoName.IsNull()) {
    fHistDeltaRhovsCent = new TH2F("DeltaRhovsCent", "DetlaRhovsCent", 101, -1, 100, fNbins, -fMaxBinPt, fMaxBinPt);
    fOutput->Add(fHistDeltaRhovsCent);
    if (!fTracksName.IsNull()) {
      fHistDeltaRhovsNtrack = new TH2F("DeltaRhovsNtrack", "DeltaRhovsNtrack", 125, 0, 4000, fNbins, -fMaxBinPt, fMaxBinPt);
      fOutput->Add(fHistDeltaRhovsNtrack);
    }
  }

  if (fScaleFunction) {
    fHistRhoScaledvsCent = new TH2F("RhoScaledvsCent", "RhoScalevsCent", 101, -1, 100, fNbins, fMinBinPt , fMaxBinPt*2);
    fOutput->Add(fHistRhoScaledvsCent);

    if (!fTracksName.IsNull()) {
      fHistRhoScaledvsNtrack = new TH2F("RhoScaledvsNtrack", "RhoScaledvsNtrack", 125, 0, 4000, fNbins, fMinBinPt, fMaxBinPt*2);
      fOutput->Add(fHistRhoScaledvsNtrack);
    }

    if (!fCaloName.IsNull()) {
      fHistRhoScaledvsNcluster = new TH2F("RhoScaledvsNcluster", "RhoScaledvsNcluster", 50, 0, 1500, fNbins, fMinBinPt, fMaxBinPt*2);
      fOutput->Add(fHistRhoScaledvsNcluster);
    }

    if (!fCompareRhoScaledName.IsNull()) {
      fHistDeltaRhoScalevsCent = new TH2F("DeltaRhoScalevsCent", "DeltaRhoScalevsCent", 101, -1, 100, fNbins, -fMaxBinPt, fMaxBinPt);
      fOutput->Add(fHistDeltaRhoScalevsCent);
      
      if (!fTracksName.IsNull()) {
	fHistDeltaRhoScalevsNtrack = new TH2F("DeltaRhoScalevsNtrack", "DeltaRhoScalevsNtrack", 125, 0, 4000, fNbins, -fMaxBinPt, fMaxBinPt);
	fOutput->Add(fHistDeltaRhoScalevsNtrack);
      }
    }
  }
}

//________________________________________________________________________
Bool_t AliAnalysisTaskRhoBase::Run() 
{
  // Run the analysis.

  Double_t rho = GetRhoFactor(fCent);
  fRho->SetVal(rho);

  if (fScaleFunction) {
    Double_t rhoScaled = rho * GetScaleFactor(fCent);
    fRhoScaled->SetVal(rhoScaled);
  }

  return kTRUE;
}

//________________________________________________________________________
Bool_t AliAnalysisTaskRhoBase::FillHistograms() 
{
  // Fill histograms.

  Int_t Ntracks   = 0;
  Int_t Nclusters = 0;

  if (fTracks)
    Ntracks   = fTracks->GetEntriesFast();
  if (fCaloClusters)
    Nclusters = fCaloClusters->GetEntriesFast();

  if (fJets) {
    Int_t    Njets         = fJets->GetEntries();
    Int_t    NjetAcc       = 0;
    Int_t    NjetUE1Sigma  = 0;
    Int_t    NjetUE2Sigma  = 0;
    Int_t    NjetUE3Sigma  = 0;
    Double_t rhoPlus1Sigma = fRho->GetVal() + fInEventSigmaRho;
    Double_t rhoPlus2Sigma = fRho->GetVal() + 2*fInEventSigmaRho;
    Double_t rhoPlus3Sigma = fRho->GetVal() + 3*fInEventSigmaRho;

    for (Int_t i = 0; i < Njets; ++i) {
      
      AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(i));
      if (!jet) {
	AliError(Form("%s: Could not receive jet %d", GetName(), i));
	continue;
      } 
      
      if (!AcceptJet(jet))
	continue;
      
      fHistJetPtvsCent->Fill(fCent, jet->Pt());
      fHistJetAreavsCent->Fill(fCent, jet->Area());
      fHistJetRhovsCent->Fill(fCent, jet->Pt() / jet->Area());
      fHistJetRhovsEta[fCentBin]->Fill(jet->Pt() / jet->Area(), jet->Eta());

      if (fTracks) {
	fHistJetPtvsNtrack->Fill(Ntracks, jet->Pt());
	fHistJetAreavsNtrack->Fill(Ntracks, jet->Area());
      }

      fHistJetNconstVsPt[fCentBin]->Fill(jet->GetNumberOfConstituents(), jet->Pt());

      if (jet->Pt() < rhoPlus1Sigma * jet->Area())
	NjetUE1Sigma++;

      if (jet->Pt() < rhoPlus2Sigma * jet->Area())
	NjetUE2Sigma++;

      if (jet->Pt() < rhoPlus3Sigma * jet->Area())
	NjetUE3Sigma++;
      
      NjetAcc++;
    }
    
    if (NjetAcc>0) {
      fHistNjUEoverNjVsNj[fCentBin*3  ]->Fill(NjetAcc,1.*NjetUE1Sigma/NjetAcc);
      fHistNjUEoverNjVsNj[fCentBin*3+1]->Fill(NjetAcc,1.*NjetUE2Sigma/NjetAcc);
      fHistNjUEoverNjVsNj[fCentBin*3+2]->Fill(NjetAcc,1.*NjetUE3Sigma/NjetAcc);
    }

    fHistNjetvsCent->Fill(fCent, NjetAcc);
    if (fTracks)
      fHistNjetvsNtrack->Fill(Ntracks, NjetAcc);
  }
  
  fHistRhovsCent->Fill(fCent, fRho->GetVal());

  if (fTracks)
    fHistRhovsNtrack->Fill(Ntracks, fRho->GetVal());
  if (fCaloClusters)
    fHistRhovsNcluster->Fill(Nclusters, fRho->GetVal());
  if (fCompareRho) {
    fHistDeltaRhovsCent->Fill(fCent, fRho->GetVal() - fCompareRho->GetVal());
    if (fTracks)
      fHistDeltaRhovsNtrack->Fill(Ntracks, fRho->GetVal() - fCompareRho->GetVal());
  }

  if (fRhoScaled) {
    fHistRhoScaledvsCent->Fill(fCent, fRhoScaled->GetVal());
    if (fTracks)
      fHistRhoScaledvsNtrack->Fill(Ntracks, fRhoScaled->GetVal());
    if (fCaloClusters)
      fHistRhoScaledvsNcluster->Fill(Nclusters,  fRhoScaled->GetVal());
    if (fCompareRhoScaled) {
      fHistDeltaRhoScalevsCent->Fill(fCent, fRhoScaled->GetVal() - fCompareRhoScaled->GetVal());
      if (fTracks)
	fHistDeltaRhoScalevsNtrack->Fill(Ntracks, fRhoScaled->GetVal() - fCompareRhoScaled->GetVal());
    }
  }

  return kTRUE;
}      


//________________________________________________________________________
void AliAnalysisTaskRhoBase::ExecOnce() 
{
  // Init the analysis.

  if (!fRho) {
    fRho = new AliRhoParameter(fRhoName, 0);

    if (fAttachToEvent) {
      if (!(InputEvent()->FindListObject(fRhoName))) {
	InputEvent()->AddObject(fRho);
      } else {
	AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fRhoName.Data()));
	return;
      }
    }
    
    if (fScaleFunction && !fRhoScaled) {
      fRhoScaled = new AliRhoParameter(fRhoScaledName, 0);
      if (fAttachToEvent) {
	if (!(InputEvent()->FindListObject(fRhoScaledName))) {
	  InputEvent()->AddObject(fRhoScaled);
	} else {
	  AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fRhoScaledName.Data()));
	  return;
	}
      }
    }
  }

  if (!fCompareRhoName.IsNull() && !fCompareRho) {
    fCompareRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fCompareRhoName));
    if (!fCompareRho) {
      AliWarning(Form("%s: Could not retrieve rho %s!", GetName(), fCompareRhoName.Data()));
    }
  }

  if (!fCompareRhoScaledName.IsNull() && !fCompareRhoScaled) {
    fCompareRhoScaled = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fCompareRhoScaledName));
    if (!fCompareRhoScaled) {
      AliWarning(Form("%s: Could not retrieve rho %s!", GetName(), fCompareRhoScaledName.Data()));
    }
  }

  AliAnalysisTaskEmcalJet::ExecOnce();
}

//________________________________________________________________________
Double_t AliAnalysisTaskRhoBase::GetRhoFactor(Double_t cent)
{
  // Return rho per centrality.

  Double_t rho = 0;
  if (fRhoFunction)
    rho = fRhoFunction->Eval(cent);
  return rho;
}

//________________________________________________________________________
Double_t AliAnalysisTaskRhoBase::GetScaleFactor(Double_t cent)
{
  // Get scale factor.

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