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

// $Id$
//
// Jet spectrum task.
//
// Author: R.Reed, M.Connors

#include "AliAnalysisTaskEmcalJetSpectra.h"

#include <TCanvas.h>
#include <TChain.h>
#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TParameter.h>
#include <TParticle.h>
#include <TTree.h>
#include <TVector3.h>

#include "AliAODEvent.h"
#include "AliAnalysisManager.h"
#include "AliAnalysisTask.h"
#include "AliCentrality.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliEmcalJet.h"
#include "AliVCluster.h"

ClassImp(AliAnalysisTaskEmcalJetSpectra)

//________________________________________________________________________
AliAnalysisTaskEmcalJetSpectra::AliAnalysisTaskEmcalJetSpectra() : 
  AliAnalysisTaskSE(), 
  fTracksName("tracks"),
  fJetsName("jets"),
  fClustersName("clusters"),
  fRhos1Name("fArrRhos"),
  fRhos2Name("fArrRhos"),
  fRhos3Name("fArrRhos"),
  fPhimin(-10), 
  fPhimax(10),
  fEtamin(-0.9), 
  fEtamax(0.9),
  fAreacut(0.0),
  fESD(0), 
  fOutputList(0), 
  fHistCentrality(0),  
  fHistDeltaRho12vsCent(0), 
  fHistDeltaRho13vsCent(0), 
  fHistDeltaRho23vsCent(0), 
  fHistDeltaJetPt12vsCent(0), 
  fHistDeltaJetPt13vsCent(0), 
  fHistDeltaJetPt23vsCent(0), 
  fHistRho1vsCent(0), 
  fHistRho2vsCent(0), 
  fHistRho3vsCent(0)
{
  // Default constructor.

  for (Int_t i = 0;i<6;++i) {
    fHistRawJetPt[i]       = 0;
    fHistAreavsRawPt[i]    = 0;
    for (Int_t j = 0;j<4;++j) {
      fHistNEFvsPt[i][j]   = 0;
       fHistZvsPt[i][j]     = 0;
       fHistZchvsPt[i][j]   = 0;
       fHistZemvsPt[i][j]   = 0;
       fHistJetPt[i][j]     = 0;
       fHistNconsvsPt[i][j] = 0;
       fHistJetPt5[i][j]    = 0;
       fHistJetPt6[i][j]    = 0;
       fHistJetPt7[i][j]    = 0;
       fHistJetPt8[i][j]    = 0;
    }
  }
}

//________________________________________________________________________
AliAnalysisTaskEmcalJetSpectra::AliAnalysisTaskEmcalJetSpectra(const char *name) :
  AliAnalysisTaskSE(name), 
  fTracksName("tracks"),
  fJetsName("jets"),
  fClustersName("clusters"),
  fRhos1Name("fArrRhos"),
  fRhos2Name("fArrRhos"),
  fRhos3Name("fArrRhos"),
  fPhimin(-10), 
  fPhimax(10),
  fEtamin(-0.9), 
  fEtamax(0.9),
  fAreacut(0.0),
  fESD(0), 
  fOutputList(0), 
  fHistCentrality(0),  
  fHistDeltaRho12vsCent(0), 
  fHistDeltaRho13vsCent(0), 
  fHistDeltaRho23vsCent(0), 
  fHistDeltaJetPt12vsCent(0), 
  fHistDeltaJetPt13vsCent(0), 
  fHistDeltaJetPt23vsCent(0), 
  fHistRho1vsCent(0), 
  fHistRho2vsCent(0), 
  fHistRho3vsCent(0)
{
  // Constructor

  for (Int_t i = 0;i<6;++i) {
    fHistRawJetPt[i]       = 0;
    fHistAreavsRawPt[i]    = 0;
    for (Int_t j = 0;j<4;++j) {
      fHistNEFvsPt[i][j]   = 0;
      fHistZvsPt[i][j]     = 0;
      fHistZchvsPt[i][j]   = 0;
      fHistZemvsPt[i][j]   = 0;
      fHistJetPt[i][j]     = 0;
      fHistNconsvsPt[i][j] = 0;
      fHistJetPt5[i][j]    = 0;
      fHistJetPt6[i][j]    = 0;
      fHistJetPt7[i][j]    = 0;
      fHistJetPt8[i][j]    = 0;
    }
  }
  
  DefineInput(0, TChain::Class());
  DefineOutput(1, TList::Class());
}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetSpectra::UserCreateOutputObjects()
{
  // Called once at the beginning of the analysis.

  AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
  if (!handler) {
    AliError("Input handler not available!");
    return;
  }

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

  fHistCentrality            = new TH1F("Centrality",             "Centrality",            101, -1,  100);
  fHistDeltaRho12vsCent      = new TH2F("DeltaRho12vsCent",       "DeltaRho12vsCent",      100, 0.0, 100.0, 500,-250,250);
  fHistDeltaRho13vsCent      = new TH2F("DeltaRho13vsCent",       "DeltaRho13vsCent",      100, 0.0, 100.0, 500,-250,250);
  fHistDeltaRho23vsCent      = new TH2F("DeltaRho23vsCent",       "DeltaRho23vsCent",      100, 0.0, 100.0, 500,-250,250);
  fHistDeltaJetPt12vsCent    = new TH2F("DeltaJetPt12vsCent",     "DeltaJetPt12vsCent",    100, 0.0, 100.0, 500,-250,250);
  fHistDeltaJetPt13vsCent    = new TH2F("DeltaJetPt13vsCent",     "DeltaJetPt13vsCent",    100, 0.0, 100.0, 500,-250,250);
  fHistDeltaJetPt23vsCent    = new TH2F("DeltaJetPt23vsCent",     "DeltaJetPt23vsCent",    100, 0.0, 100.0, 500,-250,250);
  fHistRho1vsCent            = new TH2F("Rho1vsCent",             "Rho1vsCent",            100, 0.0, 100.0, 500, 0, 500);
  fHistRho2vsCent            = new TH2F("Rho2vsCent",             "Rho2vsCent",            100, 0.0, 100.0, 500, 0, 500);
  fHistRho3vsCent            = new TH2F("Rho3vsCent",             "Rho3vsCent",            100, 0.0, 100.0, 500, 0, 500);

  for (Int_t i = 0;i<6;++i) {
    TString name00(Form("fHistRawJetPt_%i",i));
    fHistRawJetPt[i] = new TH1F(name00,name00,250,0,500);
    fOutputList->Add(fHistRawJetPt[i]);
    TString name01(Form("fHistAreavsRawPt_%i",i));
    fHistAreavsRawPt[i] = new TH2F(name01,name01,250,0,500,100,0,1);
    fOutputList->Add(fHistAreavsRawPt[i]);
    for (Int_t j = 0;j<4;j++) {
      TString name0(Form("fHistNEFvsPt_%i_%i",i,j));
      fHistNEFvsPt[i][j] = new TH2F(name0,name0,250,-250,250,200,0,2);
      fOutputList->Add(fHistNEFvsPt[i][j]);
      TString name1(Form("fHistZvsPt_%i_%i",i,j));
      fHistZvsPt[i][j] = new TH2F(name1,name1,250,-250,250,200,0,2);
      fOutputList->Add(fHistZvsPt[i][j]);
      TString name2(Form("fHistZchvsPt_%i_%i",i,j));
      fHistZchvsPt[i][j] = new TH2F(name2,name2,250,-250,250,200,0,2);
      fOutputList->Add(fHistZchvsPt[i][j]);
      TString name3(Form("fHistZemvsPt_%i_%i",i,j));
      fHistZemvsPt[i][j] = new TH2F(name3,name3,250,-250,250,200,0,2);
      fOutputList->Add(fHistZemvsPt[i][j]);
      TString name5(Form("fHistJetPt_%i_%i",i,j));
      fHistJetPt[i][j] = new TH1F(name5,name5,250,-250,250);
      fOutputList->Add(fHistJetPt[i][j]);
      TString name6(Form("fHistNconsvsPt_%i_%i",i,j));
      fHistNconsvsPt[i][j] = new TH2F(name6,name6,250,-250,250,500,0,500);
      fOutputList->Add(fHistNconsvsPt[i][j]);
      TString name7(Form("fHistJetPt5_%i_%i",i,j));
      fHistJetPt5[i][j] = new TH1F(name7,name7,250,-250,250);
      fOutputList->Add(fHistJetPt5[i][j]);
      TString name8(Form("fHistJetPt6_%i_%i",i,j));
      fHistJetPt6[i][j] = new TH1F(name8,name8,250,-250,250);
      fOutputList->Add(fHistJetPt6[i][j]);
      TString name9(Form("fHistJetPt7_%i_%i",i,j));
      fHistJetPt7[i][j] = new TH1F(name9,name9,250,-250,250);
      fOutputList->Add(fHistJetPt7[i][j]);
      TString name10(Form("fHistJetPt8_%i_%i",i,j));
      fHistJetPt8[i][j] = new TH1F(name10,name10,250,-250,250);
      fOutputList->Add(fHistJetPt8[i][j]);
    }
  }
  fOutputList->Add(fHistCentrality);
  fOutputList->Add(fHistDeltaRho12vsCent);
  fOutputList->Add(fHistDeltaRho13vsCent);
  fOutputList->Add(fHistDeltaRho23vsCent);
  fOutputList->Add(fHistDeltaJetPt12vsCent);
  fOutputList->Add(fHistDeltaJetPt13vsCent);
  fOutputList->Add(fHistDeltaJetPt23vsCent);
  fOutputList->Add(fHistRho1vsCent);
  fOutputList->Add(fHistRho2vsCent);
  fOutputList->Add(fHistRho3vsCent);
  
  PostData(1, fOutputList);
}

//________________________________________________________________________
Int_t AliAnalysisTaskEmcalJetSpectra::GetCentBin(Double_t cent) const 
{
  // Get centrality bin.

  Int_t centbin = -1;
  if (cent>=0 && cent<10)
    centbin = 0;
  else if (cent>=10 && cent<20)
    centbin = 1;
  else if (cent>=20 && cent<30)
    centbin = 2;
  else if (cent>=30 && cent<40)
    centbin = 3;
  else if (cent>=40 && cent<50)
    centbin = 4;
  else if (cent>=50 && cent<90)
    centbin = 5;
  return centbin;
}

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

  // esd or aod mode
  Bool_t esdMode = kTRUE;
  if (dynamic_cast<AliAODEvent*>(InputEvent()))
    esdMode = kFALSE;

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

  // get centrality 
  Double_t fCent = -1; 
  TList *l = InputEvent()->GetList();
  AliCentrality *centrality = InputEvent()->GetCentrality() ;
  if (centrality)
    fCent = centrality->GetCentralityPercentile("V0M");
  else
    fCent=99; // probably pp data
  if (fCent<0) {
    AliError(Form("Centrality negative: %f", fCent));
    return;
  }

  Int_t centbin = GetCentBin(fCent);

  //don't want to analyze events 90-100
  if (centbin < 0)
    return;

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

  tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
  if (!tracks) {
    AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
    return;
  }
    
  jets = dynamic_cast<TClonesArray*>(l->FindObject(fJetsName));
  if (!jets) {
    AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
    return;
  }

  Double_t rho1 = -1;
  TParameter<Double_t> *Rho1Param = dynamic_cast<TParameter<Double_t>*>(InputEvent()->FindListObject(fRhos1Name));
  if (Rho1Param)
    rho1 = Rho1Param->GetVal();

  Double_t rho2 = -1;
  TParameter<Double_t> *Rho2Param = dynamic_cast<TParameter<Double_t>*>(InputEvent()->FindListObject(fRhos2Name));
  if (Rho2Param)
    rho2 = Rho2Param->GetVal();

  Double_t rho3 = -1;
  TParameter<Double_t> *Rho3Param = dynamic_cast<TParameter<Double_t>*>(InputEvent()->FindListObject(fRhos3Name));
  if (Rho3Param)
    rho3 = Rho3Param->GetVal();

  if (rho1>0)
    fHistRho1vsCent->Fill(fCent,rho1);
  if (rho2>0)
    fHistRho2vsCent->Fill(fCent,rho2);
  if (rho3>0)
    fHistRho3vsCent->Fill(fCent,rho3);

  if (( rho1>0 ) && ( rho2>0 )) 
    fHistDeltaRho12vsCent->Fill(fCent,rho1-rho2);
  if (( rho1>0 ) && ( rho3>0 )) 
    fHistDeltaRho13vsCent->Fill(fCent,rho1-rho3);
  if (( rho2>0 ) && ( rho3>0 )) 
    fHistDeltaRho23vsCent->Fill(fCent,rho2-rho3);
  fHistCentrality->Fill(fCent);

  Double_t lJetPt   = -500;
  Double_t lJetPtun = -500;
  const Int_t Njets = jets->GetEntries();
  Int_t NjetAcc = 0;
 
  for (Int_t iJets = 0; iJets < Njets; ++iJets) {
    AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(iJets));
    if (!jet)
      continue; 
    if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax))
      continue;
    if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax))
      continue;
    fHistAreavsRawPt[centbin]->Fill(jet->Pt(),jet->Area());
    if (jet->Area()<fAreacut)
      continue;
    //prevents 0 area jets from sneaking by when area cut == 0
    if (jet->Area()==0)
      continue;
    fHistRawJetPt[centbin]->Fill(jet->Pt());
    
    NjetAcc++;
    Double_t Z; 
    Double_t jetPt1 = -500;
    Double_t jetPt2 = -500;
    Double_t jetPt3 = -500;
    if (rho1 > 0)
      jetPt1 = jet->Pt()-jet->Area()*rho1;
    if (rho2 > 0)
      jetPt2 = jet->Pt()-jet->Area()*rho2;
    if (rho3 > 0)
      jetPt3 = jet->Pt()-jet->Area()*rho3;
    if (( rho1>0 ) && ( rho2>0 )) 
      fHistDeltaJetPt12vsCent->Fill(fCent,jetPt1-jetPt2);
    if (( rho1>0 ) && ( rho3>0 )) 
      fHistDeltaJetPt13vsCent->Fill(fCent,jetPt1-jetPt3);
    if (( rho2>0 ) && ( rho3>0 )) 
      fHistDeltaJetPt23vsCent->Fill(fCent,jetPt2-jetPt3);

    if (lJetPt < jetPt2) {
      lJetPt= jetPt2;
      lJetPtun = jet->Pt();
    }

    if (jet->MaxTrackPt()>jet->MaxClusterPt()) {
      Z = jet->MaxTrackPt();
    }
    else{
      Z = jet->MaxClusterPt();
    }

    fHistNEFvsPt[centbin][0]->Fill(jet->Pt(),jet->NEF());
    fHistZvsPt[centbin][0]->Fill(jet->Pt(),Z/jet->Pt());
    fHistZchvsPt[centbin][0]->Fill(jet->Pt(),jet->MaxTrackPt()/jet->Pt());
    fHistZemvsPt[centbin][0]->Fill(jet->Pt(),jet->MaxClusterPt()/jet->Pt());
    fHistJetPt[centbin][0]->Fill(jet->Pt());
    fHistNconsvsPt[centbin][0]->Fill(jet->Pt(),jet->N());
    if ((jet->MaxTrackPt()>5) || (jet->MaxClusterPt()>5))
      fHistJetPt5[centbin][0]->Fill(jet->Pt());
    if ((jet->MaxTrackPt()>6) || (jet->MaxClusterPt()>6))
      fHistJetPt6[centbin][0]->Fill(jet->Pt());
    if ((jet->MaxTrackPt()>7) || (jet->MaxClusterPt()>7))
      fHistJetPt7[centbin][0]->Fill(jet->Pt());
    if ((jet->MaxTrackPt()>8) || (jet->MaxClusterPt()>8))
      fHistJetPt8[centbin][0]->Fill(jet->Pt());

    fHistNEFvsPt[centbin][1]->Fill(jetPt1,jet->NEF());
    fHistZvsPt[centbin][1]->Fill(jetPt1,Z/jetPt1);
    fHistZchvsPt[centbin][1]->Fill(jetPt1,jet->MaxTrackPt()/jetPt1);
    fHistZemvsPt[centbin][1]->Fill(jetPt1,jet->MaxClusterPt()/jetPt1);
    fHistJetPt[centbin][1]->Fill(jetPt1);
    fHistNconsvsPt[centbin][1]->Fill(jetPt1,jet->N());
    if ((jet->MaxTrackPt()>5) || (jet->MaxClusterPt()>5))
      fHistJetPt5[centbin][1]->Fill(jetPt1);
    if ((jet->MaxTrackPt()>6) || (jet->MaxClusterPt()>6))
      fHistJetPt6[centbin][1]->Fill(jetPt1);
    if ((jet->MaxTrackPt()>7) || (jet->MaxClusterPt()>7))
      fHistJetPt7[centbin][1]->Fill(jetPt1);
    if ((jet->MaxTrackPt()>8) || (jet->MaxClusterPt()>8))
      fHistJetPt8[centbin][1]->Fill(jetPt1);

    fHistNEFvsPt[centbin][2]->Fill(jetPt2,jet->NEF());
    fHistZvsPt[centbin][2]->Fill(jetPt2,Z/jetPt2);
    fHistZchvsPt[centbin][2]->Fill(jetPt2,jet->MaxTrackPt()/jetPt2);
    fHistZemvsPt[centbin][2]->Fill(jetPt2,jet->MaxClusterPt()/jetPt2);
    fHistJetPt[centbin][2]->Fill(jetPt2);
    fHistNconsvsPt[centbin][2]->Fill(jetPt2,jet->N());
    if ((jet->MaxTrackPt()>5) || (jet->MaxClusterPt()>5))
      fHistJetPt5[centbin][2]->Fill(jetPt2);
    if ((jet->MaxTrackPt()>6) || (jet->MaxClusterPt()>6))
      fHistJetPt6[centbin][2]->Fill(jetPt2);
    if ((jet->MaxTrackPt()>7) || (jet->MaxClusterPt()>7))
      fHistJetPt7[centbin][2]->Fill(jetPt2);
    if ((jet->MaxTrackPt()>8) || (jet->MaxClusterPt()>8))
      fHistJetPt8[centbin][2]->Fill(jetPt2);

    fHistNEFvsPt[centbin][3]->Fill(jetPt3,jet->NEF());
    fHistZvsPt[centbin][3]->Fill(jetPt3,Z/jetPt3);
    fHistZchvsPt[centbin][3]->Fill(jetPt3,jet->MaxTrackPt()/jetPt3);
    fHistZemvsPt[centbin][3]->Fill(jetPt3,jet->MaxClusterPt()/jetPt3);
    fHistJetPt[centbin][3]->Fill(jetPt3);
    fHistNconsvsPt[centbin][3]->Fill(jetPt3,jet->N());
    if ((jet->MaxTrackPt()>5) || (jet->MaxClusterPt()>5))
      fHistJetPt5[centbin][3]->Fill(jetPt3);
    if ((jet->MaxTrackPt()>6) || (jet->MaxClusterPt()>6))
      fHistJetPt6[centbin][3]->Fill(jetPt3);
    if ((jet->MaxTrackPt()>7) || (jet->MaxClusterPt()>7))
      fHistJetPt7[centbin][3]->Fill(jetPt3);
    if ((jet->MaxTrackPt()>8) || (jet->MaxClusterPt()>8))
      fHistJetPt8[centbin][3]->Fill(jetPt3);
  }

  PostData(1, fOutputList);
}      

//________________________________________________________________________
void AliAnalysisTaskEmcalJetSpectra::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}
Commit	Line	Data
9993af2d	1	// $Id$
	2	//
	3	// Jet spectrum task.
	4	//
	5	// Author: R.Reed, M.Connors
	6
	7	#include "AliAnalysisTaskEmcalJetSpectra.h"
	8
	9	#include <TCanvas.h>
74ec9ac9	10	#include <TChain.h>
9993af2d	11	#include <TClonesArray.h>
74ec9ac9	12	#include <TH1F.h>
74ec9ac9	13	#include <TH2F.h>
9993af2d	14	#include <TList.h>
020052e4	15	#include <TLorentzVector.h>
74ec9ac9	16	#include <TParameter.h>
9993af2d	17	#include <TParticle.h>
	18	#include <TTree.h>
	19	#include <TVector3.h>
020052e4	20
9993af2d	21	#include "AliAODEvent.h"
020052e4	22	#include "AliAnalysisManager.h"
9993af2d	23	#include "AliAnalysisTask.h"
9993af2d	24	#include "AliCentrality.h"
020052e4	25	#include "AliESDEvent.h"
020052e4	26	#include "AliESDInputHandler.h"
9993af2d	27	#include "AliEmcalJet.h"
9993af2d	28	#include "AliVCluster.h"
020052e4	29
	30	ClassImp(AliAnalysisTaskEmcalJetSpectra)
	31
	32	//________________________________________________________________________
9993af2d	33	AliAnalysisTaskEmcalJetSpectra::AliAnalysisTaskEmcalJetSpectra() :
	34	AliAnalysisTaskSE(),
	35	fTracksName("tracks"),
	36	fJetsName("jets"),
	37	fClustersName("clusters"),
	38	fRhos1Name("fArrRhos"),
	39	fRhos2Name("fArrRhos"),
	40	fRhos3Name("fArrRhos"),
	41	fPhimin(-10),
	42	fPhimax(10),
	43	fEtamin(-0.9),
	44	fEtamax(0.9),
	45	fAreacut(0.0),
	46	fESD(0),
	47	fOutputList(0),
	48	fHistCentrality(0),
	49	fHistDeltaRho12vsCent(0),
	50	fHistDeltaRho13vsCent(0),
	51	fHistDeltaRho23vsCent(0),
	52	fHistDeltaJetPt12vsCent(0),
	53	fHistDeltaJetPt13vsCent(0),
	54	fHistDeltaJetPt23vsCent(0),
	55	fHistRho1vsCent(0),
	56	fHistRho2vsCent(0),
	57	fHistRho3vsCent(0)
	58	{
	59	// Default constructor.
020052e4	60
9993af2d	61	for (Int_t i = 0;i<6;++i) {
	62	fHistRawJetPt[i] = 0;
	63	fHistAreavsRawPt[i] = 0;
	64	for (Int_t j = 0;j<4;++j) {
	65	fHistNEFvsPt[i][j] = 0;
020052e4	66	fHistZvsPt[i][j] = 0;
	67	fHistZchvsPt[i][j] = 0;
	68	fHistZemvsPt[i][j] = 0;
020052e4	69	fHistJetPt[i][j] = 0;
020052e4	70	fHistNconsvsPt[i][j] = 0;
6172360f	71	fHistJetPt5[i][j] = 0;
	72	fHistJetPt6[i][j] = 0;
	73	fHistJetPt7[i][j] = 0;
	74	fHistJetPt8[i][j] = 0;
9993af2d	75	}
9993af2d	76	}
020052e4	77	}
9993af2d	78
020052e4	79	//________________________________________________________________________
9993af2d	80	AliAnalysisTaskEmcalJetSpectra::AliAnalysisTaskEmcalJetSpectra(const char *name) :
	81	AliAnalysisTaskSE(name),
	82	fTracksName("tracks"),
	83	fJetsName("jets"),
	84	fClustersName("clusters"),
	85	fRhos1Name("fArrRhos"),
	86	fRhos2Name("fArrRhos"),
	87	fRhos3Name("fArrRhos"),
	88	fPhimin(-10),
	89	fPhimax(10),
	90	fEtamin(-0.9),
	91	fEtamax(0.9),
	92	fAreacut(0.0),
	93	fESD(0),
	94	fOutputList(0),
	95	fHistCentrality(0),
	96	fHistDeltaRho12vsCent(0),
	97	fHistDeltaRho13vsCent(0),
	98	fHistDeltaRho23vsCent(0),
	99	fHistDeltaJetPt12vsCent(0),
	100	fHistDeltaJetPt13vsCent(0),
	101	fHistDeltaJetPt23vsCent(0),
	102	fHistRho1vsCent(0),
	103	fHistRho2vsCent(0),
	104	fHistRho3vsCent(0)
020052e4	105	{
020052e4	106	// Constructor
9993af2d	107
	108	for (Int_t i = 0;i<6;++i) {
	109	fHistRawJetPt[i] = 0;
	110	fHistAreavsRawPt[i] = 0;
	111	for (Int_t j = 0;j<4;++j) {
	112	fHistNEFvsPt[i][j] = 0;
	113	fHistZvsPt[i][j] = 0;
	114	fHistZchvsPt[i][j] = 0;
	115	fHistZemvsPt[i][j] = 0;
	116	fHistJetPt[i][j] = 0;
	117	fHistNconsvsPt[i][j] = 0;
	118	fHistJetPt5[i][j] = 0;
	119	fHistJetPt6[i][j] = 0;
	120	fHistJetPt7[i][j] = 0;
	121	fHistJetPt8[i][j] = 0;
	122	}
	123	}
020052e4	124
	125	DefineInput(0, TChain::Class());
	126	DefineOutput(1, TList::Class());
	127	}
	128
	129	//________________________________________________________________________
	130	void AliAnalysisTaskEmcalJetSpectra::UserCreateOutputObjects()
	131	{
9993af2d	132	// Called once at the beginning of the analysis.
020052e4	133
	134	AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
	135	if (!handler) {
	136	AliError("Input handler not available!");
	137	return;
	138	}
	139
020052e4	140	OpenFile(1);
	141	fOutputList = new TList();
	142	fOutputList->SetOwner();
	143
9993af2d	144	fHistCentrality = new TH1F("Centrality", "Centrality", 101, -1, 100);
	145	fHistDeltaRho12vsCent = new TH2F("DeltaRho12vsCent", "DeltaRho12vsCent", 100, 0.0, 100.0, 500,-250,250);
	146	fHistDeltaRho13vsCent = new TH2F("DeltaRho13vsCent", "DeltaRho13vsCent", 100, 0.0, 100.0, 500,-250,250);
	147	fHistDeltaRho23vsCent = new TH2F("DeltaRho23vsCent", "DeltaRho23vsCent", 100, 0.0, 100.0, 500,-250,250);
	148	fHistDeltaJetPt12vsCent = new TH2F("DeltaJetPt12vsCent", "DeltaJetPt12vsCent", 100, 0.0, 100.0, 500,-250,250);
	149	fHistDeltaJetPt13vsCent = new TH2F("DeltaJetPt13vsCent", "DeltaJetPt13vsCent", 100, 0.0, 100.0, 500,-250,250);
	150	fHistDeltaJetPt23vsCent = new TH2F("DeltaJetPt23vsCent", "DeltaJetPt23vsCent", 100, 0.0, 100.0, 500,-250,250);
020052e4	151	fHistRho1vsCent = new TH2F("Rho1vsCent", "Rho1vsCent", 100, 0.0, 100.0, 500, 0, 500);
	152	fHistRho2vsCent = new TH2F("Rho2vsCent", "Rho2vsCent", 100, 0.0, 100.0, 500, 0, 500);
	153	fHistRho3vsCent = new TH2F("Rho3vsCent", "Rho3vsCent", 100, 0.0, 100.0, 500, 0, 500);
	154
9993af2d	155	for (Int_t i = 0;i<6;++i) {
	156	TString name00(Form("fHistRawJetPt_%i",i));
	157	fHistRawJetPt[i] = new TH1F(name00,name00,250,0,500);
	158	fOutputList->Add(fHistRawJetPt[i]);
	159	TString name01(Form("fHistAreavsRawPt_%i",i));
	160	fHistAreavsRawPt[i] = new TH2F(name01,name01,250,0,500,100,0,1);
	161	fOutputList->Add(fHistAreavsRawPt[i]);
	162	for (Int_t j = 0;j<4;j++) {
	163	TString name0(Form("fHistNEFvsPt_%i_%i",i,j));
	164	fHistNEFvsPt[i][j] = new TH2F(name0,name0,250,-250,250,200,0,2);
	165	fOutputList->Add(fHistNEFvsPt[i][j]);
	166	TString name1(Form("fHistZvsPt_%i_%i",i,j));
	167	fHistZvsPt[i][j] = new TH2F(name1,name1,250,-250,250,200,0,2);
	168	fOutputList->Add(fHistZvsPt[i][j]);
	169	TString name2(Form("fHistZchvsPt_%i_%i",i,j));
	170	fHistZchvsPt[i][j] = new TH2F(name2,name2,250,-250,250,200,0,2);
	171	fOutputList->Add(fHistZchvsPt[i][j]);
	172	TString name3(Form("fHistZemvsPt_%i_%i",i,j));
	173	fHistZemvsPt[i][j] = new TH2F(name3,name3,250,-250,250,200,0,2);
	174	fOutputList->Add(fHistZemvsPt[i][j]);
	175	TString name5(Form("fHistJetPt_%i_%i",i,j));
	176	fHistJetPt[i][j] = new TH1F(name5,name5,250,-250,250);
	177	fOutputList->Add(fHistJetPt[i][j]);
	178	TString name6(Form("fHistNconsvsPt_%i_%i",i,j));
	179	fHistNconsvsPt[i][j] = new TH2F(name6,name6,250,-250,250,500,0,500);
	180	fOutputList->Add(fHistNconsvsPt[i][j]);
	181	TString name7(Form("fHistJetPt5_%i_%i",i,j));
	182	fHistJetPt5[i][j] = new TH1F(name7,name7,250,-250,250);
	183	fOutputList->Add(fHistJetPt5[i][j]);
	184	TString name8(Form("fHistJetPt6_%i_%i",i,j));
	185	fHistJetPt6[i][j] = new TH1F(name8,name8,250,-250,250);
	186	fOutputList->Add(fHistJetPt6[i][j]);
	187	TString name9(Form("fHistJetPt7_%i_%i",i,j));
	188	fHistJetPt7[i][j] = new TH1F(name9,name9,250,-250,250);
	189	fOutputList->Add(fHistJetPt7[i][j]);
	190	TString name10(Form("fHistJetPt8_%i_%i",i,j));
	191	fHistJetPt8[i][j] = new TH1F(name10,name10,250,-250,250);
	192	fOutputList->Add(fHistJetPt8[i][j]);
	193	}
020052e4	194	}
020052e4	195	fOutputList->Add(fHistCentrality);
020052e4	196	fOutputList->Add(fHistDeltaRho12vsCent);
	197	fOutputList->Add(fHistDeltaRho13vsCent);
	198	fOutputList->Add(fHistDeltaRho23vsCent);
6172360f	199	fOutputList->Add(fHistDeltaJetPt12vsCent);
	200	fOutputList->Add(fHistDeltaJetPt13vsCent);
	201	fOutputList->Add(fHistDeltaJetPt23vsCent);
020052e4	202	fOutputList->Add(fHistRho1vsCent);
	203	fOutputList->Add(fHistRho2vsCent);
	204	fOutputList->Add(fHistRho3vsCent);
	205
9993af2d	206	PostData(1, fOutputList);
020052e4	207	}
9993af2d	208
020052e4	209	//________________________________________________________________________
9993af2d	210	Int_t AliAnalysisTaskEmcalJetSpectra::GetCentBin(Double_t cent) const
	211	{
	212	// Get centrality bin.
020052e4	213
020052e4	214	Int_t centbin = -1;
	215	if (cent>=0 && cent<10)
	216	centbin = 0;
	217	else if (cent>=10 && cent<20)
	218	centbin = 1;
	219	else if (cent>=20 && cent<30)
	220	centbin = 2;
	221	else if (cent>=30 && cent<40)
	222	centbin = 3;
	223	else if (cent>=40 && cent<50)
	224	centbin = 4;
	225	else if (cent>=50 && cent<90)
	226	centbin = 5;
	227	return centbin;
	228	}
9993af2d	229
020052e4	230	//________________________________________________________________________
	231	void AliAnalysisTaskEmcalJetSpectra::UserExec(Option_t *)
	232	{
	233	// Main loop, called for each event.
	234
	235	// esd or aod mode
	236	Bool_t esdMode = kTRUE;
	237	if (dynamic_cast<AliAODEvent*>(InputEvent()))
	238	esdMode = kFALSE;
	239
9993af2d	240	if (esdMode) {
	241	// optimization in case autobranch loading is off
	242	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
	243	if (fTracksName == "Tracks")
	244	am->LoadBranch("Tracks");
	245	}
020052e4	246
	247	// get centrality
	248	Double_t fCent = -1;
	249	TList *l = InputEvent()->GetList();
	250	AliCentrality *centrality = InputEvent()->GetCentrality() ;
	251	if (centrality)
	252	fCent = centrality->GetCentralityPercentile("V0M");
	253	else
	254	fCent=99; // probably pp data
	255	if (fCent<0) {
	256	AliError(Form("Centrality negative: %f", fCent));
	257	return;
	258	}
	259
	260	Int_t centbin = GetCentBin(fCent);
	261
	262	//don't want to analyze events 90-100
	263	if (centbin < 0)
	264	return;
	265
	266	TClonesArray *jets = 0;
	267	TClonesArray *tracks = 0;
020052e4	268
	269	tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
	270	if (!tracks) {
	271	AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
	272	return;
	273	}
	274
	275	jets = dynamic_cast<TClonesArray*>(l->FindObject(fJetsName));
	276	if (!jets) {
	277	AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
	278	return;
	279	}
	280
6172360f	281	Double_t rho1 = -1;
	282	TParameter<Double_t> Rho1Param = dynamic_cast<TParameter<Double_t>>(InputEvent()->FindListObject(fRhos1Name));
	283	if (Rho1Param)
	284	rho1 = Rho1Param->GetVal();
74ec9ac9	285
6172360f	286	Double_t rho2 = -1;
	287	TParameter<Double_t> Rho2Param = dynamic_cast<TParameter<Double_t>>(InputEvent()->FindListObject(fRhos2Name));
	288	if (Rho2Param)
	289	rho2 = Rho2Param->GetVal();
020052e4	290
6172360f	291	Double_t rho3 = -1;
	292	TParameter<Double_t> Rho3Param = dynamic_cast<TParameter<Double_t>>(InputEvent()->FindListObject(fRhos3Name));
	293	if (Rho3Param)
	294	rho3 = Rho3Param->GetVal();
020052e4	295
6172360f	296	if (rho1>0)
020052e4	297	fHistRho1vsCent->Fill(fCent,rho1);
	298	if (rho2>0)
	299	fHistRho2vsCent->Fill(fCent,rho2);
	300	if (rho3>0)
	301	fHistRho3vsCent->Fill(fCent,rho3);
	302
	303	if (( rho1>0 ) && ( rho2>0 ))
	304	fHistDeltaRho12vsCent->Fill(fCent,rho1-rho2);
	305	if (( rho1>0 ) && ( rho3>0 ))
	306	fHistDeltaRho13vsCent->Fill(fCent,rho1-rho3);
	307	if (( rho2>0 ) && ( rho3>0 ))
	308	fHistDeltaRho23vsCent->Fill(fCent,rho2-rho3);
	309	fHistCentrality->Fill(fCent);
	310
9993af2d	311	Double_t lJetPt = -500;
020052e4	312	Double_t lJetPtun = -500;
020052e4	313	const Int_t Njets = jets->GetEntries();
020052e4	314	Int_t NjetAcc = 0;
020052e4	315
020052e4	316	for (Int_t iJets = 0; iJets < Njets; ++iJets) {
020052e4	317	AliEmcalJet jet = static_cast<AliEmcalJet>(jets->At(iJets));
	318	if (!jet)
	319	continue;
9993af2d	320	if ((jet->Phi()<fPhimin)\|\|(jet->Phi()>fPhimax))
020052e4	321	continue;
9993af2d	322	if ((jet->Eta()<fEtamin)\|\|(jet->Eta()>fEtamax))
020052e4	323	continue;
020052e4	324	fHistAreavsRawPt[centbin]->Fill(jet->Pt(),jet->Area());
9993af2d	325	if (jet->Area()<fAreacut)
020052e4	326	continue;
	327	//prevents 0 area jets from sneaking by when area cut == 0
	328	if (jet->Area()==0)
	329	continue;
	330	fHistRawJetPt[centbin]->Fill(jet->Pt());
	331
	332	NjetAcc++;
9993af2d	333	Double_t Z;
	334	Double_t jetPt1 = -500;
	335	Double_t jetPt2 = -500;
	336	Double_t jetPt3 = -500;
020052e4	337	if (rho1 > 0)
	338	jetPt1 = jet->Pt()-jet->Area()*rho1;
	339	if (rho2 > 0)
	340	jetPt2 = jet->Pt()-jet->Area()*rho2;
	341	if (rho3 > 0)
	342	jetPt3 = jet->Pt()-jet->Area()*rho3;
9993af2d	343	if (( rho1>0 ) && ( rho2>0 ))
	344	fHistDeltaJetPt12vsCent->Fill(fCent,jetPt1-jetPt2);
	345	if (( rho1>0 ) && ( rho3>0 ))
	346	fHistDeltaJetPt13vsCent->Fill(fCent,jetPt1-jetPt3);
	347	if (( rho2>0 ) && ( rho3>0 ))
	348	fHistDeltaJetPt23vsCent->Fill(fCent,jetPt2-jetPt3);
	349
	350	if (lJetPt < jetPt2) {
020052e4	351	lJetPt= jetPt2;
	352	lJetPtun = jet->Pt();
	353	}
	354
9993af2d	355	if (jet->MaxTrackPt()>jet->MaxClusterPt()) {
6172360f	356	Z = jet->MaxTrackPt();
020052e4	357	}
020052e4	358	else{
6172360f	359	Z = jet->MaxClusterPt();
020052e4	360	}
020052e4	361
74ec9ac9	362	fHistNEFvsPt[centbin][0]->Fill(jet->Pt(),jet->NEF());
6172360f	363	fHistZvsPt[centbin][0]->Fill(jet->Pt(),Z/jet->Pt());
	364	fHistZchvsPt[centbin][0]->Fill(jet->Pt(),jet->MaxTrackPt()/jet->Pt());
	365	fHistZemvsPt[centbin][0]->Fill(jet->Pt(),jet->MaxClusterPt()/jet->Pt());
74ec9ac9	366	fHistJetPt[centbin][0]->Fill(jet->Pt());
74ec9ac9	367	fHistNconsvsPt[centbin][0]->Fill(jet->Pt(),jet->N());
6172360f	368	if ((jet->MaxTrackPt()>5) \|\| (jet->MaxClusterPt()>5))
	369	fHistJetPt5[centbin][0]->Fill(jet->Pt());
	370	if ((jet->MaxTrackPt()>6) \|\| (jet->MaxClusterPt()>6))
	371	fHistJetPt6[centbin][0]->Fill(jet->Pt());
	372	if ((jet->MaxTrackPt()>7) \|\| (jet->MaxClusterPt()>7))
	373	fHistJetPt7[centbin][0]->Fill(jet->Pt());
	374	if ((jet->MaxTrackPt()>8) \|\| (jet->MaxClusterPt()>8))
	375	fHistJetPt8[centbin][0]->Fill(jet->Pt());
74ec9ac9	376
74ec9ac9	377	fHistNEFvsPt[centbin][1]->Fill(jetPt1,jet->NEF());
6172360f	378	fHistZvsPt[centbin][1]->Fill(jetPt1,Z/jetPt1);
	379	fHistZchvsPt[centbin][1]->Fill(jetPt1,jet->MaxTrackPt()/jetPt1);
	380	fHistZemvsPt[centbin][1]->Fill(jetPt1,jet->MaxClusterPt()/jetPt1);
74ec9ac9	381	fHistJetPt[centbin][1]->Fill(jetPt1);
74ec9ac9	382	fHistNconsvsPt[centbin][1]->Fill(jetPt1,jet->N());
6172360f	383	if ((jet->MaxTrackPt()>5) \|\| (jet->MaxClusterPt()>5))
	384	fHistJetPt5[centbin][1]->Fill(jetPt1);
	385	if ((jet->MaxTrackPt()>6) \|\| (jet->MaxClusterPt()>6))
	386	fHistJetPt6[centbin][1]->Fill(jetPt1);
	387	if ((jet->MaxTrackPt()>7) \|\| (jet->MaxClusterPt()>7))
	388	fHistJetPt7[centbin][1]->Fill(jetPt1);
	389	if ((jet->MaxTrackPt()>8) \|\| (jet->MaxClusterPt()>8))
	390	fHistJetPt8[centbin][1]->Fill(jetPt1);
74ec9ac9	391
74ec9ac9	392	fHistNEFvsPt[centbin][2]->Fill(jetPt2,jet->NEF());
6172360f	393	fHistZvsPt[centbin][2]->Fill(jetPt2,Z/jetPt2);
	394	fHistZchvsPt[centbin][2]->Fill(jetPt2,jet->MaxTrackPt()/jetPt2);
	395	fHistZemvsPt[centbin][2]->Fill(jetPt2,jet->MaxClusterPt()/jetPt2);
74ec9ac9	396	fHistJetPt[centbin][2]->Fill(jetPt2);
74ec9ac9	397	fHistNconsvsPt[centbin][2]->Fill(jetPt2,jet->N());
6172360f	398	if ((jet->MaxTrackPt()>5) \|\| (jet->MaxClusterPt()>5))
	399	fHistJetPt5[centbin][2]->Fill(jetPt2);
	400	if ((jet->MaxTrackPt()>6) \|\| (jet->MaxClusterPt()>6))
	401	fHistJetPt6[centbin][2]->Fill(jetPt2);
	402	if ((jet->MaxTrackPt()>7) \|\| (jet->MaxClusterPt()>7))
	403	fHistJetPt7[centbin][2]->Fill(jetPt2);
	404	if ((jet->MaxTrackPt()>8) \|\| (jet->MaxClusterPt()>8))
	405	fHistJetPt8[centbin][2]->Fill(jetPt2);
74ec9ac9	406
74ec9ac9	407	fHistNEFvsPt[centbin][3]->Fill(jetPt3,jet->NEF());
6172360f	408	fHistZvsPt[centbin][3]->Fill(jetPt3,Z/jetPt3);
	409	fHistZchvsPt[centbin][3]->Fill(jetPt3,jet->MaxTrackPt()/jetPt3);
	410	fHistZemvsPt[centbin][3]->Fill(jetPt3,jet->MaxClusterPt()/jetPt3);
74ec9ac9	411	fHistJetPt[centbin][3]->Fill(jetPt3);
74ec9ac9	412	fHistNconsvsPt[centbin][3]->Fill(jetPt3,jet->N());
6172360f	413	if ((jet->MaxTrackPt()>5) \|\| (jet->MaxClusterPt()>5))
	414	fHistJetPt5[centbin][3]->Fill(jetPt3);
	415	if ((jet->MaxTrackPt()>6) \|\| (jet->MaxClusterPt()>6))
	416	fHistJetPt6[centbin][3]->Fill(jetPt3);
	417	if ((jet->MaxTrackPt()>7) \|\| (jet->MaxClusterPt()>7))
	418	fHistJetPt7[centbin][3]->Fill(jetPt3);
	419	if ((jet->MaxTrackPt()>8) \|\| (jet->MaxClusterPt()>8))
	420	fHistJetPt8[centbin][3]->Fill(jetPt3);
020052e4	421	}
020052e4	422
020052e4	423	PostData(1, fOutputList);
	424	}
	425
	426	//________________________________________________________________________
	427	void AliAnalysisTaskEmcalJetSpectra::Terminate(Option_t *)
	428	{
9993af2d	429	// Called once at the end of the analysis.
020052e4	430	}
	431
	432