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

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

#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"

#include "AliAnalysisTaskRho.h"

ClassImp(AliAnalysisTaskRho)

//________________________________________________________________________
AliAnalysisTaskRho::AliAnalysisTaskRho() : 
  AliAnalysisTaskRhoBase(),
  fTracksName("tracks"),
  fJetsName("KtJets"),
  fClustersName("caloClusters"),
  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"),
  fClustersName("caloClusters"),
  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

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

//________________________________________________________________________
void AliAnalysisTaskRho::UserCreateOutputObjects()
{
  AliAnalysisTaskRhoBase::UserCreateOutputObjects();

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

  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};
  Int_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] = maxJetPts[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 *) 
{

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