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

// $Id$
//
// Emcal jet response matrix maker task.
//
// Author: S. Aiola

#include "AliJetResponseMaker.h"

#include <TChain.h>
#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TList.h>
#include <TLorentzVector.h>

#include "AliAnalysisManager.h"
#include "AliCentrality.h"
#include "AliFJWrapper.h"
#include "AliVCluster.h"
#include "AliVTrack.h"
#include "AliEmcalJet.h"
#include "AliMCEvent.h"

ClassImp(AliJetResponseMaker)

//________________________________________________________________________
AliJetResponseMaker::AliJetResponseMaker() : 
  AliAnalysisTaskEmcal("AliJetResponseMaker"),
  fMCTracksName("MCParticles"),
  fMCJetsName("MCJets"),
  fMaxDistance(0.2),
  fMCTracks(0),
  fMCJets(0),
  fHistMCJetPhiEta(0),
  fHistMCJetsPt(0),
  fHistMCJetsPtNonBias(0),
  fHistMCJetsNEFvsPt(0),
  fHistMCJetsZvsPt(0),
  fHistJetPhiEta(0),
  fHistJetsPt(0),
  fHistJetsPtNonBias(0),
  fHistJetsNEFvsPt(0),
  fHistJetsZvsPt(0),
  fHistClosestDistance(0),
  fHistClosestDeltaPhi(0),
  fHistClosestDeltaEta(0),
  fHistClosestDeltaPt(0),
  fHistPartvsDetecPt(0)
{
  // Default constructor.
}

//________________________________________________________________________
AliJetResponseMaker::AliJetResponseMaker(const char *name) : 
  AliAnalysisTaskEmcal(name),
  fMCTracksName("MCParticles"),
  fMCJetsName("MCJets"),
  fMaxDistance(0.2),
  fMCTracks(0),
  fMCJets(0),
  fHistMCJetPhiEta(0),
  fHistMCJetsPt(0),
  fHistMCJetsPtNonBias(0),
  fHistMCJetsNEFvsPt(0),
  fHistMCJetsZvsPt(0),
  fHistJetPhiEta(0),
  fHistJetsPt(0),
  fHistJetsPtNonBias(0),
  fHistJetsNEFvsPt(0),
  fHistJetsZvsPt(0),
  fHistClosestDistance(0),
  fHistClosestDeltaPhi(0),
  fHistClosestDeltaEta(0),
  fHistClosestDeltaPt(0),
  fHistPartvsDetecPt(0)
{
  // Standard constructor.

}

//________________________________________________________________________
AliJetResponseMaker::~AliJetResponseMaker()
{
  // Destructor
}

//________________________________________________________________________
void AliJetResponseMaker::UserCreateOutputObjects()
{
  // Create user objects.

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

  fHistJetPhiEta = new TH2F("fHistJetPhiEta", "fHistJetPhiEta", 20, -2, 2, 32, 0, 6.4);
  fHistJetPhiEta->GetXaxis()->SetTitle("#eta");
  fHistJetPhiEta->GetYaxis()->SetTitle("#phi");
  fOutput->Add(fHistJetPhiEta);
  
  fHistJetsPt = new TH1F("fHistJetsPt", "fHistJetsPt", fNbins, fMinPt, fMaxPt);
  fHistJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  fHistJetsPt->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistJetsPt);
  
  fHistJetsPtNonBias = new TH1F("fHistJetsPtNonBias", "fHistJetsPtNonBias", fNbins, fMinPt, fMaxPt);
  fHistJetsPtNonBias->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  fHistJetsPtNonBias->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistJetsPtNonBias);
  
  fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
  fHistJetsZvsPt->GetXaxis()->SetTitle("Z");
  fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
  fOutput->Add(fHistJetsZvsPt);
  
  if (fAnaType == kEMCAL) {
    fHistJetsNEFvsPt = new TH2F("fHistJetsNEFvsPt", "fHistJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
    fHistJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
    fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
    fOutput->Add(fHistJetsNEFvsPt);
  }

  fHistMCJetPhiEta = new TH2F("fHistMCJetPhiEta", "fHistMCJetPhiEta", 20, -2, 2, 32, 0, 6.4);
  fHistMCJetPhiEta->GetXaxis()->SetTitle("#eta");
  fHistMCJetPhiEta->GetYaxis()->SetTitle("#phi");
  fOutput->Add(fHistMCJetPhiEta);
  
  fHistMCJetsPt = new TH1F("fHistMCJetsPt", "fHistMCJetsPt", fNbins, fMinPt, fMaxPt);
  fHistMCJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  fHistMCJetsPt->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistMCJetsPt);
  
  fHistMCJetsPtNonBias = new TH1F("fHistMCJetsPtNonBias", "fHistMCJetsPtNonBias", fNbins, fMinPt, fMaxPt);
  fHistMCJetsPtNonBias->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  fHistMCJetsPtNonBias->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistMCJetsPtNonBias);
  
  fHistMCJetsZvsPt = new TH2F("fHistMCJetsZvsPt", "fHistMCJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
  fHistMCJetsZvsPt->GetXaxis()->SetTitle("Z");
  fHistMCJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
  fOutput->Add(fHistMCJetsZvsPt);
  
  if (fAnaType == kEMCAL) {
    fHistMCJetsNEFvsPt = new TH2F("fHistMCJetsNEFvsPt", "fHistMCJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
    fHistMCJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
    fHistMCJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
    fOutput->Add(fHistMCJetsNEFvsPt);
  }

  fHistClosestDistance = new TH1F("fHistClosestDistance", "fHistClosestDistance", 50, 0, fMaxDistance * 1.2);
  fHistClosestDistance->GetXaxis()->SetTitle("d");
  fHistClosestDistance->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistClosestDistance);

  fHistClosestDeltaPhi = new TH1F("fHistClosestDeltaPhi", "fHistClosestDeltaPhi", 64, -1.6, 4.8);
  fHistClosestDeltaPhi->GetXaxis()->SetTitle("#Delta#phi");
  fHistClosestDeltaPhi->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistClosestDeltaPhi);

  fHistClosestDeltaEta = new TH1F("fHistClosestDeltaEta", "fHistClosestDeltaEta", TMath::CeilNint(fMaxEta - fMinEta) * 20, fMinEta * 2, fMaxEta * 2);
  fHistClosestDeltaEta->GetXaxis()->SetTitle("#Delta#eta");
  fHistClosestDeltaEta->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistClosestDeltaEta);

  fHistClosestDeltaPt = new TH1F("fHistClosestDeltaPt", "fHistClosestDeltaPt", fNbins, -fMaxPt / 2, fMaxPt / 2);
  fHistClosestDeltaPt->GetXaxis()->SetTitle("#Delta p_{T}");
  fHistClosestDeltaPt->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistClosestDeltaPt);

  fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
  fHistPartvsDetecPt->GetXaxis()->SetTitle("p_{T}^{det}");
  fHistPartvsDetecPt->GetYaxis()->SetTitle("p_{T}^{rec}");
  fOutput->Add(fHistPartvsDetecPt);

  PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
}

//________________________________________________________________________
void AliJetResponseMaker::FillHistograms()
{
  // Fill histograms.

  // Find the closest jets
  DoJetLoop(fJets, fMCJets);
  DoJetLoop(fMCJets, fJets);

  Int_t nJets = fJets->GetEntriesFast();

  for (Int_t i = 0; i < nJets; i++) {

    AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJets->At(i));

    if (!jet) {
      AliError(Form("Could not receive jet %d", i));
      continue;
    }  

    if (!AcceptJet(jet))
      continue;

    fHistJetsPtNonBias->Fill(jet->Pt());

    if (jet->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC || jet->MaxClusterPt() < fPtBiasJetClus))
	continue;

    if (jet->ClosestJet() && jet->ClosestJet()->ClosestJet() == jet && jet->ClosestJetDistance() < fMaxDistance) {    // Matched jet found
      jet->SetMatchedToClosest();
      jet->ClosestJet()->SetMatchedToClosest();
      fHistClosestDistance->Fill(jet->ClosestJetDistance());
      Double_t deta = jet->Eta() - jet->MatchedJet()->Eta();
      fHistClosestDeltaEta->Fill(deta);
      Double_t dphi = jet->Phi() - jet->MatchedJet()->Phi();
      fHistClosestDeltaPhi->Fill(dphi);
      Double_t dpt = jet->Pt() - jet->MatchedJet()->Pt();
      fHistClosestDeltaPt->Fill(dpt);
      fHistPartvsDetecPt->Fill(jet->Pt(), jet->MatchedJet()->Pt());
    }

    fHistJetsPt->Fill(jet->Pt());

    fHistJetPhiEta->Fill(jet->Eta(), jet->Phi());

    if (fAnaType == kEMCAL)
      fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt());

    for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
      AliVParticle *track = jet->TrackAt(it, fTracks);
      if (track)
	fHistJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt());
    }

    for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
      AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
      if (cluster) {
	TLorentzVector nP;
	cluster->GetMomentum(nP, fVertex);
	fHistJetsZvsPt->Fill(nP.Pt() / jet->Pt(), jet->Pt());
      }
    }
  }

 Int_t nMCJets = fMCJets->GetEntriesFast();

  for (Int_t i = 0; i < nMCJets; i++) {

    AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fMCJets->At(i));

    if (!jet) {
      AliError(Form("Could not receive mc jet %d", i));
      continue;
    }  

    if (!AcceptJet(jet))
      continue;

    fHistMCJetsPtNonBias->Fill(jet->Pt());

    if (jet->MaxTrackPt() < fPtBiasJetTrack)
	continue;

    fHistMCJetsPt->Fill(jet->Pt());

    fHistMCJetPhiEta->Fill(jet->Eta(), jet->Phi());

    if (fAnaType == kEMCAL)
      fHistMCJetsNEFvsPt->Fill(jet->NEF(), jet->Pt());

    for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
      AliVParticle *track = jet->TrackAt(it, fMCTracks);
      if (track)
	fHistMCJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt());
    }
  }
}

//________________________________________________________________________
void AliJetResponseMaker::DoJetLoop(TClonesArray *jets1, TClonesArray *jets2)
{
  // Do the jet loop.

  Int_t nJets1 = jets1->GetEntriesFast();
  Int_t nJets2 = jets2->GetEntriesFast();

  for (Int_t i = 0; i < nJets1; i++) {

    AliEmcalJet* jet1 = dynamic_cast<AliEmcalJet*>(jets1->At(i));

    if (!jet1) {
      AliError(Form("Could not receive jet %d", i));
      continue;
    }  

    if (!AcceptJet(jet1))
      continue;
    
    if (jet1->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC || jet1->IsMC() || jet1->MaxClusterPt() < fPtBiasJetClus))
	continue;

    for (Int_t j = 0; j < nJets2; j++) {
      
      AliEmcalJet* jet2 = dynamic_cast<AliEmcalJet*>(jets2->At(j));
      
      if (!jet2) {
	AliError(Form("Could not receive jet %d", j));
	continue;
      }  
      
      if (!AcceptJet(jet2))
	continue;
      
      if (jet2->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC || jet2->IsMC() || jet2->MaxClusterPt() < fPtBiasJetClus))
	continue;
      
      Double_t deta = jet2->Eta() - jet1->Eta();
      Double_t dphi = jet2->Phi() - jet1->Phi();
      Double_t d = TMath::Sqrt(deta * deta + dphi * dphi);

      if (d < jet1->ClosestJetDistance()) {
	jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
	jet1->SetClosestJet(jet2, d);
      }
      else if (d < jet1->SecondClosestJetDistance()) {
	jet1->SetSecondClosestJet(jet2, d);
      }
    }
  }
}

//________________________________________________________________________
void AliJetResponseMaker::RetrieveEventObjects()
{
  // Retrieve event objects.

  AliAnalysisTaskEmcal::RetrieveEventObjects();
  
  if (!fMCJetsName.IsNull()) {
    fMCJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCJetsName));
    if (!fMCJets) {
      AliWarning(Form("Could not retrieve MC jets %s!", fMCJetsName.Data())); 
    }
  }

  if (!fMCTracksName.IsNull()) {
    fMCTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCTracksName));
    if (!fMCJets) {
      AliWarning(Form("Could not retrieve MC tracks %s!", fMCTracksName.Data())); 
    }
  }
}

//________________________________________________________________________
void AliJetResponseMaker::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}
Commit	Line	Data
7549c451	1	// $Id$
2949a2ec	2	//
	3	// Emcal jet response matrix maker task.
	4	//
	5	// Author: S. Aiola
	6
	7	#include "AliJetResponseMaker.h"
	8
	9	#include <TChain.h>
	10	#include <TClonesArray.h>
	11	#include <TH1F.h>
	12	#include <TH2F.h>
	13	#include <TList.h>
	14	#include <TLorentzVector.h>
	15
	16	#include "AliAnalysisManager.h"
	17	#include "AliCentrality.h"
	18	#include "AliFJWrapper.h"
	19	#include "AliVCluster.h"
	20	#include "AliVTrack.h"
	21	#include "AliEmcalJet.h"
	22	#include "AliMCEvent.h"
	23
	24	ClassImp(AliJetResponseMaker)
	25
	26	//________________________________________________________________________
	27	AliJetResponseMaker::AliJetResponseMaker() :
	28	AliAnalysisTaskEmcal("AliJetResponseMaker"),
7549c451	29	fMCTracksName("MCParticles"),
7549c451	30	fMCJetsName("MCJets"),
1b76c28f	31	fMaxDistance(0.2),
7549c451	32	fMCTracks(0),
2949a2ec	33	fMCJets(0),
	34	fHistMCJetPhiEta(0),
	35	fHistMCJetsPt(0),
2949a2ec	36	fHistMCJetsPtNonBias(0),
2949a2ec	37	fHistMCJetsNEFvsPt(0),
	38	fHistMCJetsZvsPt(0),
	39	fHistJetPhiEta(0),
	40	fHistJetsPt(0),
2949a2ec	41	fHistJetsPtNonBias(0),
2949a2ec	42	fHistJetsNEFvsPt(0),
1b76c28f	43	fHistJetsZvsPt(0),
	44	fHistClosestDistance(0),
	45	fHistClosestDeltaPhi(0),
	46	fHistClosestDeltaEta(0),
	47	fHistClosestDeltaPt(0),
	48	fHistPartvsDetecPt(0)
2949a2ec	49	{
	50	// Default constructor.
	51	}
	52
	53	//________________________________________________________________________
	54	AliJetResponseMaker::AliJetResponseMaker(const char *name) :
	55	AliAnalysisTaskEmcal(name),
7549c451	56	fMCTracksName("MCParticles"),
7549c451	57	fMCJetsName("MCJets"),
1b76c28f	58	fMaxDistance(0.2),
7549c451	59	fMCTracks(0),
2949a2ec	60	fMCJets(0),
	61	fHistMCJetPhiEta(0),
	62	fHistMCJetsPt(0),
2949a2ec	63	fHistMCJetsPtNonBias(0),
2949a2ec	64	fHistMCJetsNEFvsPt(0),
	65	fHistMCJetsZvsPt(0),
	66	fHistJetPhiEta(0),
	67	fHistJetsPt(0),
2949a2ec	68	fHistJetsPtNonBias(0),
2949a2ec	69	fHistJetsNEFvsPt(0),
1b76c28f	70	fHistJetsZvsPt(0),
	71	fHistClosestDistance(0),
	72	fHistClosestDeltaPhi(0),
	73	fHistClosestDeltaEta(0),
	74	fHistClosestDeltaPt(0),
	75	fHistPartvsDetecPt(0)
2949a2ec	76	{
	77	// Standard constructor.
	78
2949a2ec	79	}
	80
	81	//________________________________________________________________________
	82	AliJetResponseMaker::~AliJetResponseMaker()
	83	{
	84	// Destructor
	85	}
	86
	87	//________________________________________________________________________
	88	void AliJetResponseMaker::UserCreateOutputObjects()
	89	{
	90	// Create user objects.
	91
	92	OpenFile(1);
	93	fOutput = new TList();
	94	fOutput->SetOwner();
	95
7549c451	96	fHistJetPhiEta = new TH2F("fHistJetPhiEta", "fHistJetPhiEta", 20, -2, 2, 32, 0, 6.4);
	97	fHistJetPhiEta->GetXaxis()->SetTitle("#eta");
	98	fHistJetPhiEta->GetYaxis()->SetTitle("#phi");
	99	fOutput->Add(fHistJetPhiEta);
	100
2949a2ec	101	fHistJetsPt = new TH1F("fHistJetsPt", "fHistJetsPt", fNbins, fMinPt, fMaxPt);
	102	fHistJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	103	fHistJetsPt->GetYaxis()->SetTitle("counts");
	104	fOutput->Add(fHistJetsPt);
	105
2949a2ec	106	fHistJetsPtNonBias = new TH1F("fHistJetsPtNonBias", "fHistJetsPtNonBias", fNbins, fMinPt, fMaxPt);
	107	fHistJetsPtNonBias->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	108	fHistJetsPtNonBias->GetYaxis()->SetTitle("counts");
	109	fOutput->Add(fHistJetsPtNonBias);
	110
2949a2ec	111	fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
	112	fHistJetsZvsPt->GetXaxis()->SetTitle("Z");
	113	fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
	114	fOutput->Add(fHistJetsZvsPt);
	115
	116	if (fAnaType == kEMCAL) {
	117	fHistJetsNEFvsPt = new TH2F("fHistJetsNEFvsPt", "fHistJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
	118	fHistJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
	119	fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
	120	fOutput->Add(fHistJetsNEFvsPt);
	121	}
7549c451	122
	123	fHistMCJetPhiEta = new TH2F("fHistMCJetPhiEta", "fHistMCJetPhiEta", 20, -2, 2, 32, 0, 6.4);
	124	fHistMCJetPhiEta->GetXaxis()->SetTitle("#eta");
	125	fHistMCJetPhiEta->GetYaxis()->SetTitle("#phi");
	126	fOutput->Add(fHistMCJetPhiEta);
	127
	128	fHistMCJetsPt = new TH1F("fHistMCJetsPt", "fHistMCJetsPt", fNbins, fMinPt, fMaxPt);
	129	fHistMCJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	130	fHistMCJetsPt->GetYaxis()->SetTitle("counts");
	131	fOutput->Add(fHistMCJetsPt);
	132
7549c451	133	fHistMCJetsPtNonBias = new TH1F("fHistMCJetsPtNonBias", "fHistMCJetsPtNonBias", fNbins, fMinPt, fMaxPt);
	134	fHistMCJetsPtNonBias->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	135	fHistMCJetsPtNonBias->GetYaxis()->SetTitle("counts");
	136	fOutput->Add(fHistMCJetsPtNonBias);
	137
7549c451	138	fHistMCJetsZvsPt = new TH2F("fHistMCJetsZvsPt", "fHistMCJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
	139	fHistMCJetsZvsPt->GetXaxis()->SetTitle("Z");
	140	fHistMCJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
	141	fOutput->Add(fHistMCJetsZvsPt);
	142
	143	if (fAnaType == kEMCAL) {
	144	fHistMCJetsNEFvsPt = new TH2F("fHistMCJetsNEFvsPt", "fHistMCJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
	145	fHistMCJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
	146	fHistMCJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
	147	fOutput->Add(fHistMCJetsNEFvsPt);
	148	}
	149
1b76c28f	150	fHistClosestDistance = new TH1F("fHistClosestDistance", "fHistClosestDistance", 50, 0, fMaxDistance * 1.2);
	151	fHistClosestDistance->GetXaxis()->SetTitle("d");
	152	fHistClosestDistance->GetYaxis()->SetTitle("counts");
	153	fOutput->Add(fHistClosestDistance);
	154
	155	fHistClosestDeltaPhi = new TH1F("fHistClosestDeltaPhi", "fHistClosestDeltaPhi", 64, -1.6, 4.8);
	156	fHistClosestDeltaPhi->GetXaxis()->SetTitle("#Delta#phi");
	157	fHistClosestDeltaPhi->GetYaxis()->SetTitle("counts");
	158	fOutput->Add(fHistClosestDeltaPhi);
	159
	160	fHistClosestDeltaEta = new TH1F("fHistClosestDeltaEta", "fHistClosestDeltaEta", TMath::CeilNint(fMaxEta - fMinEta) * 20, fMinEta * 2, fMaxEta * 2);
	161	fHistClosestDeltaEta->GetXaxis()->SetTitle("#Delta#eta");
	162	fHistClosestDeltaEta->GetYaxis()->SetTitle("counts");
	163	fOutput->Add(fHistClosestDeltaEta);
	164
	165	fHistClosestDeltaPt = new TH1F("fHistClosestDeltaPt", "fHistClosestDeltaPt", fNbins, -fMaxPt / 2, fMaxPt / 2);
	166	fHistClosestDeltaPt->GetXaxis()->SetTitle("#Delta p_{T}");
	167	fHistClosestDeltaPt->GetYaxis()->SetTitle("counts");
	168	fOutput->Add(fHistClosestDeltaPt);
	169
	170	fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
	171	fHistPartvsDetecPt->GetXaxis()->SetTitle("p_{T}^{det}");
	172	fHistPartvsDetecPt->GetYaxis()->SetTitle("p_{T}^{rec}");
	173	fOutput->Add(fHistPartvsDetecPt);
	174
7549c451	175	PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
2949a2ec	176	}
	177
	178	//________________________________________________________________________
	179	void AliJetResponseMaker::FillHistograms()
	180	{
	181	// Fill histograms.
	182
1b76c28f	183	// Find the closest jets
	184	DoJetLoop(fJets, fMCJets);
	185	DoJetLoop(fMCJets, fJets);
2949a2ec	186
1b76c28f	187	Int_t nJets = fJets->GetEntriesFast();
2949a2ec	188
1b76c28f	189	for (Int_t i = 0; i < nJets; i++) {
2949a2ec	190
1b76c28f	191	AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJets->At(i));
2949a2ec	192
2949a2ec	193	if (!jet) {
1b76c28f	194	AliError(Form("Could not receive jet %d", i));
2949a2ec	195	continue;
	196	}
	197
	198	if (!AcceptJet(jet))
	199	continue;
	200
	201	fHistJetsPtNonBias->Fill(jet->Pt());
	202
2949a2ec	203	if (jet->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC \|\| jet->MaxClusterPt() < fPtBiasJetClus))
	204	continue;
	205
1b76c28f	206	if (jet->ClosestJet() && jet->ClosestJet()->ClosestJet() == jet && jet->ClosestJetDistance() < fMaxDistance) { // Matched jet found
	207	jet->SetMatchedToClosest();
	208	jet->ClosestJet()->SetMatchedToClosest();
	209	fHistClosestDistance->Fill(jet->ClosestJetDistance());
	210	Double_t deta = jet->Eta() - jet->MatchedJet()->Eta();
	211	fHistClosestDeltaEta->Fill(deta);
	212	Double_t dphi = jet->Phi() - jet->MatchedJet()->Phi();
	213	fHistClosestDeltaPhi->Fill(dphi);
	214	Double_t dpt = jet->Pt() - jet->MatchedJet()->Pt();
	215	fHistClosestDeltaPt->Fill(dpt);
	216	fHistPartvsDetecPt->Fill(jet->Pt(), jet->MatchedJet()->Pt());
2949a2ec	217	}
	218
	219	fHistJetsPt->Fill(jet->Pt());
	220
	221	fHistJetPhiEta->Fill(jet->Eta(), jet->Phi());
	222
	223	if (fAnaType == kEMCAL)
	224	fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt());
	225
	226	for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
1b76c28f	227	AliVParticle *track = jet->TrackAt(it, fTracks);
2949a2ec	228	if (track)
	229	fHistJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt());
	230	}
	231
1b76c28f	232	for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
	233	AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
	234	if (cluster) {
	235	TLorentzVector nP;
	236	cluster->GetMomentum(nP, fVertex);
	237	fHistJetsZvsPt->Fill(nP.Pt() / jet->Pt(), jet->Pt());
	238	}
	239	}
	240	}
	241
	242	Int_t nMCJets = fMCJets->GetEntriesFast();
	243
	244	for (Int_t i = 0; i < nMCJets; i++) {
	245
	246	AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fMCJets->At(i));
	247
	248	if (!jet) {
	249	AliError(Form("Could not receive mc jet %d", i));
2949a2ec	250	continue;
1b76c28f	251	}
2949a2ec	252
1b76c28f	253	if (!AcceptJet(jet))
1b76c28f	254	continue;
2949a2ec	255
1b76c28f	256	fHistMCJetsPtNonBias->Fill(jet->Pt());
	257
	258	if (jet->MaxTrackPt() < fPtBiasJetTrack)
	259	continue;
	260
	261	fHistMCJetsPt->Fill(jet->Pt());
	262
	263	fHistMCJetPhiEta->Fill(jet->Eta(), jet->Phi());
	264
	265	if (fAnaType == kEMCAL)
	266	fHistMCJetsNEFvsPt->Fill(jet->NEF(), jet->Pt());
	267
	268	for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
	269	AliVParticle *track = jet->TrackAt(it, fMCTracks);
	270	if (track)
	271	fHistMCJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt());
	272	}
	273	}
	274	}
	275
	276	//________________________________________________________________________
	277	void AliJetResponseMaker::DoJetLoop(TClonesArray jets1, TClonesArray jets2)
	278	{
	279	// Do the jet loop.
	280
	281	Int_t nJets1 = jets1->GetEntriesFast();
	282	Int_t nJets2 = jets2->GetEntriesFast();
	283
	284	for (Int_t i = 0; i < nJets1; i++) {
	285
	286	AliEmcalJet* jet1 = dynamic_cast<AliEmcalJet*>(jets1->At(i));
	287
	288	if (!jet1) {
	289	AliError(Form("Could not receive jet %d", i));
	290	continue;
	291	}
	292
	293	if (!AcceptJet(jet1))
	294	continue;
	295
	296	if (jet1->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC \|\| jet1->IsMC() \|\| jet1->MaxClusterPt() < fPtBiasJetClus))
	297	continue;
	298
	299	for (Int_t j = 0; j < nJets2; j++) {
	300
	301	AliEmcalJet* jet2 = dynamic_cast<AliEmcalJet*>(jets2->At(j));
	302
	303	if (!jet2) {
	304	AliError(Form("Could not receive jet %d", j));
	305	continue;
	306	}
	307
	308	if (!AcceptJet(jet2))
	309	continue;
	310
	311	if (jet2->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC \|\| jet2->IsMC() \|\| jet2->MaxClusterPt() < fPtBiasJetClus))
	312	continue;
	313
	314	Double_t deta = jet2->Eta() - jet1->Eta();
	315	Double_t dphi = jet2->Phi() - jet1->Phi();
	316	Double_t d = TMath::Sqrt(deta * deta + dphi * dphi);
	317
	318	if (d < jet1->ClosestJetDistance()) {
	319	jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
320	jet1->SetClosestJet(jet2, d);
321	}
322	else if (d < jet1->SecondClosestJetDistance()) {
323	jet1->SetSecondClosestJet(jet2, d);
2949a2ec	324	}
	325	}
	326	}
	327	}
	328
	329	//________________________________________________________________________
	330	void AliJetResponseMaker::RetrieveEventObjects()
	331	{
	332	// Retrieve event objects.
	333
	334	AliAnalysisTaskEmcal::RetrieveEventObjects();
	335
7549c451	336	if (!fMCJetsName.IsNull()) {
7549c451	337	fMCJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCJetsName));
2949a2ec	338	if (!fMCJets) {
	339	AliWarning(Form("Could not retrieve MC jets %s!", fMCJetsName.Data()));
	340	}
	341	}
7549c451	342
	343	if (!fMCTracksName.IsNull()) {
	344	fMCTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCTracksName));
	345	if (!fMCJets) {
	346	AliWarning(Form("Could not retrieve MC tracks %s!", fMCTracksName.Data()));
	347	}
	348	}
2949a2ec	349	}
	350
	351	//________________________________________________________________________
	352	void AliJetResponseMaker::Terminate(Option_t *)
	353	{
	354	// Called once at the end of the analysis.
	355	}