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

// $Id$
//
// Calculation of rho from a collection of jets.
// If scale function is given the scaled rho will be exported
// with the name as "fRhoName".apppend("_Scaled").
//
// Authors: R.Reed, S.Aiola

#include "AliAnalysisTaskRho.h"

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

#include "AliAnalysisManager.h"
#include "AliCentrality.h"
#include "AliEmcalJet.h"
#include "AliLog.h"
#include "AliRhoParameter.h"
#include "AliVCluster.h"
#include "AliVEventHandler.h"

ClassImp(AliAnalysisTaskRho)

//________________________________________________________________________
AliAnalysisTaskRho::AliAnalysisTaskRho() : 
  AliAnalysisTaskRhoBase(),
  fTracksName(),
  fJetsName(),
  fRhoScaledName(""),
  fPhiMin(0),
  fPhiMax(0),
  fEtaMin(0),
  fEtaMax(0),
  fAreaCut(0),
  fAreaEmcCut(0),
  fNExclLeadJets(0),
  fScaleFunction(0),
  fCreateHisto(kFALSE),
  fTracks(0),
  fJets(0),
  fOutputList(0),
  fHistCentrality(0),
  fHistJetPt(0),
  fHistJetArea(0),
  fHistRhovsCent(0),
  fHistDeltaRhovsCent(0),
  fHistDeltaRhoScalevsCent(0),
  fHistJetPtvsCent(0),
  fHistJetAreavsCent(0),
  fHistNjetvsCent(0),
  fHistRhovsNtrack(0),
  fHistDeltaRhovsNtrack(0),
  fHistDeltaRhoScalevsNtrack(0),
  fHistJetPtvsNtrack(0),
  fHistJetAreavsNtrack(0),
  fHistNjetvsNtrack(0),
  fRhoScaled(0)
{
  // Constructor.
}

//________________________________________________________________________
AliAnalysisTaskRho::AliAnalysisTaskRho(const char *name, Bool_t histo) :
  AliAnalysisTaskRhoBase(name),
  fTracksName("tracks"),
  fJetsName("KtJets"),
  fRhoScaledName(""),
  fPhiMin(0),
  fPhiMax(TMath::TwoPi()),
  fEtaMin(-0.5),
  fEtaMax(+0.5),
  fAreaCut(0.01),
  fAreaEmcCut(0),
  fNExclLeadJets(1),
  fScaleFunction(0),
  fCreateHisto(histo),
  fTracks(0),
  fJets(0),
  fOutputList(0),
  fHistCentrality(0),
  fHistJetPt(0),
  fHistJetArea(0),
  fHistRhovsCent(0),
  fHistDeltaRhovsCent(0),
  fHistDeltaRhoScalevsCent(0),
  fHistJetPtvsCent(0),
  fHistJetAreavsCent(0),
  fHistNjetvsCent(0),
  fHistRhovsNtrack(0),
  fHistDeltaRhovsNtrack(0),
  fHistDeltaRhoScalevsNtrack(0),
  fHistJetPtvsNtrack(0),
  fHistJetAreavsNtrack(0),
  fHistNjetvsNtrack(0),
  fRhoScaled(0)
{
  // Constructor.

  if (fCreateHisto)
    DefineOutput(1, TList::Class());
}

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

  AliAnalysisTaskRhoBase::UserCreateOutputObjects();

  if (fScaleFunction) {
    fRhoScaledName = fRhoName;
    fRhoScaledName += "_Scaled";
    fRhoScaled = new AliRhoParameter(fRhoScaledName, 0);  
  }

  if (!fCreateHisto)
    return;

  OpenFile(1);
  fOutputList = new TList();
  fOutputList->SetOwner();

  fHistCentrality             = new TH1F("Centrality",            "Centrality",            101, -1,  100);
  fHistRhovsCent              = new TH2F("RhovsCent",             "RhovsCent",             101, -1,  100,   500,  0,   500);
  fHistDeltaRhovsCent         = new TH2F("DeltaRhovsCent",        "DetlaRhovsCent",        101, -1,  100,   500, -250, 250);
  fHistDeltaRhoScalevsCent    = new TH2F("DeltaRhoScalevsCent",   "DeltaRhoScalevsCent",   101, -1,  100,   500, -250, 250);
  fHistJetPtvsCent            = new TH2F("JetPtvsCent",           "JetPtvsCent",           101, -1,  100,   200,  0,   500);
  fHistJetAreavsCent          = new TH2F("JetAreavsCent",         "JetAreavsCent",         101, -1,  100,   100,  0,   1.0);
  fHistNjetvsCent             = new TH2F("NjetvsCent",            "NjetvsCent",            101, -1,  100,   100,  0,   100);

  fHistRhovsNtrack            = new TH2F("RhovsNtrack",           "RhovsNtrack",           500,  0,  2500,  500,  0,   500);
  fHistDeltaRhovsNtrack       = new TH2F("DeltaRhovsNtrack",      "DeltaRhovsNtrack",      500,  0,  2500,  500, -250, 250);
  fHistDeltaRhoScalevsNtrack  = new TH2F("DeltaRhoScalevsNtrack", "DeltaRhoScalevsNtrack", 500,  0,  2500,  500, -250, 250);
  fHistJetPtvsNtrack          = new TH2F("JetPtvsNtrack",         "JetPtvsNtrack",         500,  0,  2500,  200,  0,   500);
  fHistJetAreavsNtrack        = new TH2F("JetAreavsNtrack",       "JetAreavsNtrack",       500,  0,  2500,  100,  0,   1.0);
  fHistNjetvsNtrack           = new TH2F("NjetvsNtrack",          "rNjetvsNtrack",         500,  0,  2500,  100,  0,   100);

  fHistJetPt                  = new TH1F("JetPt",                  "Jet Pt",               100,   0,    250);
  fHistJetArea                = new TH1F("JetArea",                "Jet Area",             100, 0.0,    1.0);
  
  fOutputList->Add(fHistCentrality);
  fOutputList->Add(fHistRhovsCent);
  fOutputList->Add(fHistJetPt);
  fOutputList->Add(fHistJetArea);
  fOutputList->Add(fHistDeltaRhovsCent);
  fOutputList->Add(fHistDeltaRhoScalevsCent);
  fOutputList->Add(fHistJetPtvsCent);
  fOutputList->Add(fHistJetAreavsCent);
  fOutputList->Add(fHistNjetvsCent);
  
  fOutputList->Add(fHistRhovsNtrack);
  fOutputList->Add(fHistDeltaRhovsNtrack);
  fOutputList->Add(fHistDeltaRhoScalevsNtrack);
  fOutputList->Add(fHistJetPtvsNtrack);
  fOutputList->Add(fHistJetAreavsNtrack);
  fOutputList->Add(fHistNjetvsNtrack);
  
  PostData(1, fOutputList);
}

//________________________________________________________________________
void AliAnalysisTaskRho::UserExec(Option_t *) 
{
  // Main loop, called for each event.

  if (!fIsInit) {
    ExecOnce();
    fIsInit = 1;
  }

  fRho->SetVal(0);
  if (fRhoScaled)
    fRhoScaled->SetVal(0);

  if (!fJets)
    return;

  DetermineCent();

  const Int_t Njets   = fJets->GetEntries();

  Int_t maxJetIds[]   = {-1, -1};
  Float_t maxJetPts[] = { 0,  0};

  if (fNExclLeadJets > 0) {
    for (Int_t ij = 0; ij < Njets; ++ij) {
      AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
      if (!jet) {
	AliError(Form("%s: Could not receive jet %d", GetName(), ij));
	continue;
      } 

      if (jet->Area() < fAreaCut)
        continue;
      if (jet->AreaEmc()<fAreaEmcCut)
        continue;
      if ((jet->Phi() < fPhiMin) || (jet->Phi() > fPhiMax))
        continue;
      if ((jet->Eta() < fEtaMin) || (jet->Eta() > fEtaMax))
        continue;

      if (jet->Pt() > maxJetPts[0]) {
	maxJetPts[1] = maxJetPts[0];
	maxJetIds[1] = maxJetIds[0];
	maxJetPts[0] = jet->Pt();
	maxJetIds[0] = ij;
      } else if (jet->Pt() > maxJetPts[1]) {
	maxJetPts[1] = jet->Pt();
	maxJetIds[1] = ij;
      }
    }
    if (fNExclLeadJets < 2) {
      maxJetIds[1] = -1;
      maxJetPts[1] = -1;
    }
  }

  static Double_t rhovec[999];
  Int_t NjetAcc = 0;

  // push all jets within selected acceptance into stack
  for (Int_t iJets = 0; iJets < Njets; ++iJets) {

    // exlcuding lead jets
    if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
      continue;

    AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
    if (!jet) {
      AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
      continue;
    } 

    if (jet->Area() < fAreaCut)
      continue;
    if (jet->AreaEmc()<fAreaEmcCut)
      continue;
    if ((jet->Phi() < fPhiMin) || (jet->Phi() > fPhiMax))
      continue;
    if ((jet->Eta() < fEtaMin) || (jet->Eta() > fEtaMax))
      continue;

    rhovec[NjetAcc] = jet->Pt() / jet->Area();
    ++NjetAcc;

    if (fCreateHisto) {
      // filling histograms
      const Int_t Ntracks = fTracks->GetEntriesFast();
      fHistJetPt->Fill(jet->Pt());
      fHistJetArea->Fill(jet->Area());
      fHistJetPtvsCent->Fill(fCent, jet->Pt());
      fHistJetPtvsNtrack->Fill(Ntracks, jet->Pt());
      fHistJetAreavsCent->Fill(fCent, jet->Area());
      fHistJetAreavsNtrack->Fill(Ntracks, jet->Area());
    }
  }

  if (NjetAcc > 0) {
    //find median value
    Double_t rho0      = TMath::Median(NjetAcc, rhovec);
    fRho->SetVal(rho0);
    Double_t rhoScaled = rho0;
    if (fRhoScaled) {
      rhoScaled *= GetScaleFactor(fCent);
      fRhoScaled->SetVal(rhoScaled);
    }

    if (fCreateHisto) {
      // filling other histograms
      Double_t rho        = GetRhoFactor(fCent);
      const Int_t Ntracks = fTracks->GetEntriesFast();
      fHistRhovsCent->Fill(fCent, rho0);
      fHistDeltaRhovsCent->Fill(fCent, rho0 - rho);
      fHistDeltaRhoScalevsCent->Fill(fCent, rhoScaled - rho);
      fHistRhovsNtrack->Fill(Ntracks, rho0);
      fHistDeltaRhovsNtrack->Fill(Ntracks, rho0 - rho);
      fHistDeltaRhoScalevsNtrack->Fill(Ntracks, rhoScaled - rho);
    }
  }

  if (fCreateHisto) {
    const Int_t Ntracks = fTracks->GetEntriesFast();
    fHistCentrality->Fill(fCent);
    fHistNjetvsCent->Fill(fCent, NjetAcc);
    fHistNjetvsNtrack->Fill(Ntracks, NjetAcc);
    PostData(1, fOutputList);
  }
}      

//________________________________________________________________________
void AliAnalysisTaskRho::ExecOnce() 
{
  // Initialize some settings that need to be determined in UserExec.

  AliAnalysisTaskRhoBase::ExecOnce();

  if (fRhoScaled) {
    // add rho to event if not yet there
    if (!(InputEvent()->FindListObject(fRhoScaledName))) {
      InputEvent()->AddObject(fRhoScaled);
    } else {
      AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fRhoScaledName.Data()));
      return;
    }
  }

  if (fCreateHisto) {
    fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
    if (!fTracks) {
      AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data() ));
      return;
    }
  }

  fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
  if (!fJets) {
    AliError(Form("%s: Pointer to jets %s == 0", GetName(), fJetsName.Data() ));
    return;
  }
}

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

  Double_t scale = 1;
  if (fScaleFunction)
    scale = fScaleFunction->Eval(cent);
  return scale;
}
Commit	Line	Data
5b5bad1d	1	// $Id$
192fc3f4	2	//
f09b22c5	3	// Calculation of rho from a collection of jets.
	4	// If scale function is given the scaled rho will be exported
	5	// with the name as "fRhoName".apppend("_Scaled").
192fc3f4	6	//
	7	// Authors: R.Reed, S.Aiola
	8
5b5bad1d	9	#include "AliAnalysisTaskRho.h"
5b5bad1d	10
020052e4	11	#include <TClonesArray.h>
192fc3f4	12	#include <TF1.h>
297edd60	13	#include <TH1F.h>
	14	#include <TH2F.h>
	15	#include <TList.h>
c60e0a21	16	#include <TMath.h>
020052e4	17
020052e4	18	#include "AliAnalysisManager.h"
020052e4	19	#include "AliCentrality.h"
020052e4	20	#include "AliEmcalJet.h"
297edd60	21	#include "AliLog.h"
297edd60	22	#include "AliRhoParameter.h"
020052e4	23	#include "AliVCluster.h"
297edd60	24	#include "AliVEventHandler.h"
020052e4	25
020052e4	26	ClassImp(AliAnalysisTaskRho)
	27
	28	//________________________________________________________________________
	29	AliAnalysisTaskRho::AliAnalysisTaskRho() :
192fc3f4	30	AliAnalysisTaskRhoBase(),
6a12db52	31	fTracksName(),
6a12db52	32	fJetsName(),
192fc3f4	33	fRhoScaledName(""),
	34	fPhiMin(0),
	35	fPhiMax(0),
	36	fEtaMin(0),
	37	fEtaMax(0),
	38	fAreaCut(0),
6a12db52	39	fAreaEmcCut(0),
c60e0a21	40	fNExclLeadJets(0),
192fc3f4	41	fScaleFunction(0),
c60e0a21	42	fCreateHisto(kFALSE),
f09b22c5	43	fTracks(0),
f09b22c5	44	fJets(0),
020052e4	45	fOutputList(0),
	46	fHistCentrality(0),
	47	fHistJetPt(0),
192fc3f4	48	fHistJetArea(0),
020052e4	49	fHistRhovsCent(0),
192fc3f4	50	fHistDeltaRhovsCent(0),
192fc3f4	51	fHistDeltaRhoScalevsCent(0),
020052e4	52	fHistJetPtvsCent(0),
192fc3f4	53	fHistJetAreavsCent(0),
020052e4	54	fHistNjetvsCent(0),
020052e4	55	fHistRhovsNtrack(0),
192fc3f4	56	fHistDeltaRhovsNtrack(0),
192fc3f4	57	fHistDeltaRhoScalevsNtrack(0),
020052e4	58	fHistJetPtvsNtrack(0),
	59	fHistJetAreavsNtrack(0),
	60	fHistNjetvsNtrack(0),
c60e0a21	61	fRhoScaled(0)
020052e4	62	{
6a12db52	63	// Constructor.
020052e4	64	}
020052e4	65
5d845887	66	//________________________________________________________________________
	67	AliAnalysisTaskRho::AliAnalysisTaskRho(const char *name, Bool_t histo) :
	68	AliAnalysisTaskRhoBase(name),
	69	fTracksName("tracks"),
	70	fJetsName("KtJets"),
5d845887	71	fRhoScaledName(""),
5d845887	72	fPhiMin(0),
f09b22c5	73	fPhiMax(TMath::TwoPi()),
	74	fEtaMin(-0.5),
	75	fEtaMax(+0.5),
	76	fAreaCut(0.01),
6a12db52	77	fAreaEmcCut(0),
f09b22c5	78	fNExclLeadJets(1),
5d845887	79	fScaleFunction(0),
5d845887	80	fCreateHisto(histo),
f09b22c5	81	fTracks(0),
f09b22c5	82	fJets(0),
5d845887	83	fOutputList(0),
	84	fHistCentrality(0),
	85	fHistJetPt(0),
	86	fHistJetArea(0),
	87	fHistRhovsCent(0),
	88	fHistDeltaRhovsCent(0),
	89	fHistDeltaRhoScalevsCent(0),
	90	fHistJetPtvsCent(0),
	91	fHistJetAreavsCent(0),
	92	fHistNjetvsCent(0),
	93	fHistRhovsNtrack(0),
	94	fHistDeltaRhovsNtrack(0),
	95	fHistDeltaRhoScalevsNtrack(0),
	96	fHistJetPtvsNtrack(0),
	97	fHistJetAreavsNtrack(0),
	98	fHistNjetvsNtrack(0),
	99	fRhoScaled(0)
	100	{
6a12db52	101	// Constructor.
5d845887	102
	103	if (fCreateHisto)
	104	DefineOutput(1, TList::Class());
	105	}
	106
020052e4	107	//________________________________________________________________________
	108	void AliAnalysisTaskRho::UserCreateOutputObjects()
	109	{
5b5bad1d	110	// User create output objects, called at the beginning of the analysis.
5b5bad1d	111
192fc3f4	112	AliAnalysisTaskRhoBase::UserCreateOutputObjects();
020052e4	113
f09b22c5	114	if (fScaleFunction) {
	115	fRhoScaledName = fRhoName;
	116	fRhoScaledName += "_Scaled";
	117	fRhoScaled = new AliRhoParameter(fRhoScaledName, 0);
	118	}
020052e4	119
5d845887	120	if (!fCreateHisto)
	121	return;
	122
020052e4	123	OpenFile(1);
	124	fOutputList = new TList();
	125	fOutputList->SetOwner();
	126
	127	fHistCentrality = new TH1F("Centrality", "Centrality", 101, -1, 100);
	128	fHistRhovsCent = new TH2F("RhovsCent", "RhovsCent", 101, -1, 100, 500, 0, 500);
	129	fHistDeltaRhovsCent = new TH2F("DeltaRhovsCent", "DetlaRhovsCent", 101, -1, 100, 500, -250, 250);
	130	fHistDeltaRhoScalevsCent = new TH2F("DeltaRhoScalevsCent", "DeltaRhoScalevsCent", 101, -1, 100, 500, -250, 250);
	131	fHistJetPtvsCent = new TH2F("JetPtvsCent", "JetPtvsCent", 101, -1, 100, 200, 0, 500);
	132	fHistJetAreavsCent = new TH2F("JetAreavsCent", "JetAreavsCent", 101, -1, 100, 100, 0, 1.0);
	133	fHistNjetvsCent = new TH2F("NjetvsCent", "NjetvsCent", 101, -1, 100, 100, 0, 100);
	134
	135	fHistRhovsNtrack = new TH2F("RhovsNtrack", "RhovsNtrack", 500, 0, 2500, 500, 0, 500);
	136	fHistDeltaRhovsNtrack = new TH2F("DeltaRhovsNtrack", "DeltaRhovsNtrack", 500, 0, 2500, 500, -250, 250);
	137	fHistDeltaRhoScalevsNtrack = new TH2F("DeltaRhoScalevsNtrack", "DeltaRhoScalevsNtrack", 500, 0, 2500, 500, -250, 250);
	138	fHistJetPtvsNtrack = new TH2F("JetPtvsNtrack", "JetPtvsNtrack", 500, 0, 2500, 200, 0, 500);
	139	fHistJetAreavsNtrack = new TH2F("JetAreavsNtrack", "JetAreavsNtrack", 500, 0, 2500, 100, 0, 1.0);
	140	fHistNjetvsNtrack = new TH2F("NjetvsNtrack", "rNjetvsNtrack", 500, 0, 2500, 100, 0, 100);
	141
	142	fHistJetPt = new TH1F("JetPt", "Jet Pt", 100, 0, 250);
	143	fHistJetArea = new TH1F("JetArea", "Jet Area", 100, 0.0, 1.0);
	144
	145	fOutputList->Add(fHistCentrality);
	146	fOutputList->Add(fHistRhovsCent);
	147	fOutputList->Add(fHistJetPt);
	148	fOutputList->Add(fHistJetArea);
	149	fOutputList->Add(fHistDeltaRhovsCent);
	150	fOutputList->Add(fHistDeltaRhoScalevsCent);
	151	fOutputList->Add(fHistJetPtvsCent);
	152	fOutputList->Add(fHistJetAreavsCent);
	153	fOutputList->Add(fHistNjetvsCent);
	154
	155	fOutputList->Add(fHistRhovsNtrack);
	156	fOutputList->Add(fHistDeltaRhovsNtrack);
	157	fOutputList->Add(fHistDeltaRhoScalevsNtrack);
	158	fOutputList->Add(fHistJetPtvsNtrack);
	159	fOutputList->Add(fHistJetAreavsNtrack);
	160	fOutputList->Add(fHistNjetvsNtrack);
	161
	162	PostData(1, fOutputList);
020052e4	163	}
020052e4	164
020052e4	165	//________________________________________________________________________
	166	void AliAnalysisTaskRho::UserExec(Option_t *)
	167	{
	168	// Main loop, called for each event.
e44e8726	169
f09b22c5	170	if (!fIsInit) {
	171	ExecOnce();
	172	fIsInit = 1;
	173	}
020052e4	174
a5621834	175	fRho->SetVal(0);
	176	if (fRhoScaled)
	177	fRhoScaled->SetVal(0);
	178
e44e8726	179	if (!fJets)
	180	return;
	181
f09b22c5	182	DetermineCent();
020052e4	183
f09b22c5	184	const Int_t Njets = fJets->GetEntries();
020052e4	185
f09b22c5	186	Int_t maxJetIds[] = {-1, -1};
f09b22c5	187	Float_t maxJetPts[] = { 0, 0};
c60e0a21	188
c60e0a21	189	if (fNExclLeadJets > 0) {
f09b22c5	190	for (Int_t ij = 0; ij < Njets; ++ij) {
f09b22c5	191	AliEmcalJet jet = static_cast<AliEmcalJet>(fJets->At(ij));
c60e0a21	192	if (!jet) {
f09b22c5	193	AliError(Form("%s: Could not receive jet %d", GetName(), ij));
c60e0a21	194	continue;
c60e0a21	195	}
ed777bb8	196
	197	if (jet->Area() < fAreaCut)
	198	continue;
	199	if (jet->AreaEmc()<fAreaEmcCut)
	200	continue;
	201	if ((jet->Phi() < fPhiMin) \|\| (jet->Phi() > fPhiMax))
	202	continue;
	203	if ((jet->Eta() < fEtaMin) \|\| (jet->Eta() > fEtaMax))
	204	continue;
	205
c60e0a21	206	if (jet->Pt() > maxJetPts[0]) {
c60e0a21	207	maxJetPts[1] = maxJetPts[0];
5b5bad1d	208	maxJetIds[1] = maxJetIds[0];
c60e0a21	209	maxJetPts[0] = jet->Pt();
c60e0a21	210	maxJetIds[0] = ij;
f09b22c5	211	} else if (jet->Pt() > maxJetPts[1]) {
c60e0a21	212	maxJetPts[1] = jet->Pt();
	213	maxJetIds[1] = ij;
	214	}
	215	}
	216	if (fNExclLeadJets < 2) {
	217	maxJetIds[1] = -1;
	218	maxJetPts[1] = -1;
	219	}
	220	}
	221
	222	static Double_t rhovec[999];
020052e4	223	Int_t NjetAcc = 0;
020052e4	224
c60e0a21	225	// push all jets within selected acceptance into stack
020052e4	226	for (Int_t iJets = 0; iJets < Njets; ++iJets) {
f09b22c5	227
c60e0a21	228	// exlcuding lead jets
	229	if (iJets == maxJetIds[0] \|\| iJets == maxJetIds[1])
	230	continue;
	231
f09b22c5	232	AliEmcalJet jet = static_cast<AliEmcalJet>(fJets->At(iJets));
ed777bb8	233	if (!jet) {
	234	AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
	235	continue;
	236	}
c60e0a21	237
020052e4	238	if (jet->Area() < fAreaCut)
020052e4	239	continue;
6a12db52	240	if (jet->AreaEmc()<fAreaEmcCut)
6a12db52	241	continue;
020052e4	242	if ((jet->Phi() < fPhiMin) \|\| (jet->Phi() > fPhiMax))
	243	continue;
	244	if ((jet->Eta() < fEtaMin) \|\| (jet->Eta() > fEtaMax))
	245	continue;
c60e0a21	246
c60e0a21	247	rhovec[NjetAcc] = jet->Pt() / jet->Area();
f09b22c5	248	++NjetAcc;
c60e0a21	249
	250	if (fCreateHisto) {
	251	// filling histograms
f09b22c5	252	const Int_t Ntracks = fTracks->GetEntriesFast();
c60e0a21	253	fHistJetPt->Fill(jet->Pt());
	254	fHistJetArea->Fill(jet->Area());
	255	fHistJetPtvsCent->Fill(fCent, jet->Pt());
	256	fHistJetPtvsNtrack->Fill(Ntracks, jet->Pt());
	257	fHistJetAreavsCent->Fill(fCent, jet->Area());
	258	fHistJetAreavsNtrack->Fill(Ntracks, jet->Area());
	259	}
020052e4	260	}
192fc3f4	261
f09b22c5	262	if (NjetAcc > 0) {
020052e4	263	//find median value
f09b22c5	264	Double_t rho0 = TMath::Median(NjetAcc, rhovec);
c60e0a21	265	fRho->SetVal(rho0);
f09b22c5	266	Double_t rhoScaled = rho0;
	267	if (fRhoScaled) {
	268	rhoScaled *= GetScaleFactor(fCent);
	269	fRhoScaled->SetVal(rhoScaled);
	270	}
020052e4	271
c60e0a21	272	if (fCreateHisto) {
c60e0a21	273	// filling other histograms
f09b22c5	274	Double_t rho = GetRhoFactor(fCent);
f09b22c5	275	const Int_t Ntracks = fTracks->GetEntriesFast();
c60e0a21	276	fHistRhovsCent->Fill(fCent, rho0);
f09b22c5	277	fHistDeltaRhovsCent->Fill(fCent, rho0 - rho);
f09b22c5	278	fHistDeltaRhoScalevsCent->Fill(fCent, rhoScaled - rho);
c60e0a21	279	fHistRhovsNtrack->Fill(Ntracks, rho0);
f09b22c5	280	fHistDeltaRhovsNtrack->Fill(Ntracks, rho0 - rho);
f09b22c5	281	fHistDeltaRhoScalevsNtrack->Fill(Ntracks, rhoScaled - rho);
020052e4	282	}
020052e4	283	}
020052e4	284
f09b22c5	285	if (fCreateHisto) {
	286	const Int_t Ntracks = fTracks->GetEntriesFast();
	287	fHistCentrality->Fill(fCent);
	288	fHistNjetvsCent->Fill(fCent, NjetAcc);
	289	fHistNjetvsNtrack->Fill(Ntracks, NjetAcc);
c60e0a21	290	PostData(1, fOutputList);
f09b22c5	291	}
c60e0a21	292	}
020052e4	293
f09b22c5	294	//________________________________________________________________________
	295	void AliAnalysisTaskRho::ExecOnce()
	296	{
	297	// Initialize some settings that need to be determined in UserExec.
	298
	299	AliAnalysisTaskRhoBase::ExecOnce();
	300
	301	if (fRhoScaled) {
	302	// add rho to event if not yet there
	303	if (!(InputEvent()->FindListObject(fRhoScaledName))) {
	304	InputEvent()->AddObject(fRhoScaled);
	305	} else {
	306	AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fRhoScaledName.Data()));
	307	return;
	308	}
	309	}
	310
	311	if (fCreateHisto) {
	312	fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
	313	if (!fTracks) {
	314	AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data() ));
	315	return;
	316	}
	317	}
	318
	319	fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
	320	if (!fJets) {
	321	AliError(Form("%s: Pointer to jets %s == 0", GetName(), fJetsName.Data() ));
	322	return;
	323	}
	324	}
	325
	326	//________________________________________________________________________
	327	Double_t AliAnalysisTaskRho::GetScaleFactor(Double_t cent)
	328	{
	329	// Get scale factor.
	330
	331	Double_t scale = 1;
	332	if (fScaleFunction)
	333	scale = fScaleFunction->Eval(cent);
	334	return scale;
	335	}