[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;
  }
  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;
  }
  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]);

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