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

// $Id$
//
// Calculation of rho 
//
// Authors: R.Reed, S.Aiola

#include "AliAnalysisTaskRho.h"

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

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

ClassImp(AliAnalysisTaskRho)

//________________________________________________________________________
AliAnalysisTaskRho::AliAnalysisTaskRho() : 
  AliAnalysisTaskRhoBase(),
  fTracksName("tracks"),
  fJetsName("KtJets"),
  fRhoScaledName(""),
  fPhiMin(0),
  fPhiMax(0),
  fEtaMin(0),
  fEtaMax(0),
  fAreaCut(0),
  fNExclLeadJets(0),
  fScaleFunction(0),
  fCreateHisto(kFALSE),
  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) :
  AliAnalysisTaskRhoBase(name),
  fTracksName("tracks"),
  fJetsName("KtJets"),
  fRhoScaledName(""),
  fPhiMin(0),
  fPhiMax(0),
  fEtaMin(0),
  fEtaMax(0),
  fAreaCut(0),
  fNExclLeadJets(0),
  fScaleFunction(0),
  fCreateHisto(kFALSE),
  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(0),
  fEtaMin(0),
  fEtaMax(0),
  fAreaCut(0),
  fNExclLeadJets(0),
  fScaleFunction(0),
  fCreateHisto(histo),
  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();

  fRhoScaledName = fRhoName;
  fRhoScaledName += "_Scaled";
  fRhoScaled = new TParameter<Double_t>(fRhoScaledName, 0);  

  if (!fCreateHisto)
    return;

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

  if (!fCreateHisto) {
    PostData(1, fOutputList);
    return;
  }

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

//________________________________________________________________________
Double_t AliAnalysisTaskRho::GetScaleFactor(Double_t cent)
{
  Double_t scale = 1;
  if (fScaleFunction)
    scale = fScaleFunction->Eval(cent);
  return scale;
}


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

  fRho->SetVal(-1);

   // add rho to event if not yet there
  if (!(InputEvent()->FindListObject(fRhoScaledName))) {
    new(fRhoScaled) TParameter<Double_t>(fRhoScaledName, -1);
    InputEvent()->AddObject(fRhoScaled);
  }
  else {
    fRhoScaled->SetVal(-1);
  }

  // optimization in case autobranch loading is off
  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
  if (fTracksName == "Tracks")
    am->LoadBranch("Tracks");

  TClonesArray *jets = 0;
  TClonesArray *tracks = 0;

  tracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
  if (!tracks) {
  AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
   return;
  }
    
  jets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
  if (!jets) {
    AliError(Form("Pointer to jets %s == 0", fJetsName.Data() ));
    return;
  }
  
  if (fCreateHisto)
    fHistCentrality->Fill(fCent);

  const Int_t Ntracks = tracks->GetEntries();
  const Int_t Njets = jets->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*>(jets->At(ij));
      
      if (!jet) {
	AliError(Form("Could not receive jet %d", ij));
	continue;
      } 
      
      if (jet->Pt() > maxJetPts[0]) {
	maxJetPts[1] = maxJetPts[0];
	maxJetIds[1] = maxJetIds[0];
	maxJetPts[0] = jet->Pt();
	maxJetIds[0] = ij;
      }
      
      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*>(jets->At(iJets));

    if (!jet)
      continue; 

    // applying some other cuts
    if (jet->Area() < fAreaCut)
      continue;
    if ((jet->Phi() < fPhiMin) || (jet->Phi() > fPhiMax))
      continue;
    if ((jet->Eta() < fEtaMin) || (jet->Eta() > fEtaMax))
      continue;
    if (jet->Area() == 0)
      continue;

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

    NjetAcc++;

    if (fCreateHisto) {
      // filling histograms
      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 (fCreateHisto) {
    fHistNjetvsCent->Fill(fCent, NjetAcc);
    fHistNjetvsNtrack->Fill(Ntracks, NjetAcc);
  }

  Double_t scale = GetScaleFactor(fCent);
  Double_t rhochem = GetRhoFactor(fCent);

  Double_t rho0 = -1;
  
  if (NjetAcc > 0){
    //find median value
    rho0 = TMath::Median(NjetAcc, rhovec);

    Double_t rhoScaled = rho0 * scale;

    fRho->SetVal(rho0);
    fRhoScaled->SetVal(rhoScaled);

    if (fCreateHisto) {
      // filling other histograms
      fHistRhovsCent->Fill(fCent, rho0);
      fHistDeltaRhovsCent->Fill(fCent, rho0 - rhochem);
      fHistDeltaRhoScalevsCent->Fill(fCent, rhoScaled - rhochem);
      fHistRhovsNtrack->Fill(Ntracks, rho0);
      fHistDeltaRhovsNtrack->Fill(Ntracks, rho0 - rhochem);
      fHistDeltaRhoScalevsNtrack->Fill(Ntracks, rhoScaled - rhochem);
    }
  }

  if (fCreateHisto)
    PostData(1, fOutputList);
}      

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