[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() : 
  AliAnalysisTaskEmcalJet("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) : 
  AliAnalysisTaskEmcalJet(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, fMinBinPt, fMaxBinPt);
  fHistJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  fHistJetsPt->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistJetsPt);
  
  fHistJetsPtNonBias = new TH1F("fHistJetsPtNonBias", "fHistJetsPtNonBias", fNbins, fMinBinPt, fMaxBinPt);
  fHistJetsPtNonBias->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  fHistJetsPtNonBias->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistJetsPtNonBias);
  
  fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
  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, fMinBinPt, fMaxBinPt);
    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, fMinBinPt, fMaxBinPt);
  fHistMCJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  fHistMCJetsPt->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistMCJetsPt);
  
  fHistMCJetsPtNonBias = new TH1F("fHistMCJetsPtNonBias", "fHistMCJetsPtNonBias", fNbins, fMinBinPt, fMaxBinPt);
  fHistMCJetsPtNonBias->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  fHistMCJetsPtNonBias->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistMCJetsPtNonBias);
  
  fHistMCJetsZvsPt = new TH2F("fHistMCJetsZvsPt", "fHistMCJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
  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, fMinBinPt, fMaxBinPt);
    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, -fMaxBinPt / 2, fMaxBinPt / 2);
  fHistClosestDeltaPt->GetXaxis()->SetTitle("#Delta p_{T}");
  fHistClosestDeltaPt->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistClosestDeltaPt);

  fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
  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
}

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

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

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

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

  return kTRUE;
}

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

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

  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
    return kFALSE;
  
  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())); 
    }
  }

  return kTRUE;
}

//________________________________________________________________________
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() :
6fd5039f	28	AliAnalysisTaskEmcalJet("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) :
6fd5039f	55	AliAnalysisTaskEmcalJet(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	{
2949a2ec	77	// Standard constructor.
2949a2ec	78	}
	79
	80	//________________________________________________________________________
	81	AliJetResponseMaker::~AliJetResponseMaker()
	82	{
	83	// Destructor
	84	}
	85
	86	//________________________________________________________________________
	87	void AliJetResponseMaker::UserCreateOutputObjects()
	88	{
	89	// Create user objects.
	90
	91	OpenFile(1);
	92	fOutput = new TList();
	93	fOutput->SetOwner();
	94
7549c451	95	fHistJetPhiEta = new TH2F("fHistJetPhiEta", "fHistJetPhiEta", 20, -2, 2, 32, 0, 6.4);
	96	fHistJetPhiEta->GetXaxis()->SetTitle("#eta");
	97	fHistJetPhiEta->GetYaxis()->SetTitle("#phi");
	98	fOutput->Add(fHistJetPhiEta);
	99
6fd5039f	100	fHistJetsPt = new TH1F("fHistJetsPt", "fHistJetsPt", fNbins, fMinBinPt, fMaxBinPt);
2949a2ec	101	fHistJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	102	fHistJetsPt->GetYaxis()->SetTitle("counts");
	103	fOutput->Add(fHistJetsPt);
	104
6fd5039f	105	fHistJetsPtNonBias = new TH1F("fHistJetsPtNonBias", "fHistJetsPtNonBias", fNbins, fMinBinPt, fMaxBinPt);
2949a2ec	106	fHistJetsPtNonBias->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	107	fHistJetsPtNonBias->GetYaxis()->SetTitle("counts");
	108	fOutput->Add(fHistJetsPtNonBias);
	109
6fd5039f	110	fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
2949a2ec	111	fHistJetsZvsPt->GetXaxis()->SetTitle("Z");
	112	fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
	113	fOutput->Add(fHistJetsZvsPt);
	114
	115	if (fAnaType == kEMCAL) {
6fd5039f	116	fHistJetsNEFvsPt = new TH2F("fHistJetsNEFvsPt", "fHistJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
2949a2ec	117	fHistJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
	118	fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
	119	fOutput->Add(fHistJetsNEFvsPt);
	120	}
7549c451	121
	122	fHistMCJetPhiEta = new TH2F("fHistMCJetPhiEta", "fHistMCJetPhiEta", 20, -2, 2, 32, 0, 6.4);
	123	fHistMCJetPhiEta->GetXaxis()->SetTitle("#eta");
	124	fHistMCJetPhiEta->GetYaxis()->SetTitle("#phi");
	125	fOutput->Add(fHistMCJetPhiEta);
	126
6fd5039f	127	fHistMCJetsPt = new TH1F("fHistMCJetsPt", "fHistMCJetsPt", fNbins, fMinBinPt, fMaxBinPt);
7549c451	128	fHistMCJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	129	fHistMCJetsPt->GetYaxis()->SetTitle("counts");
	130	fOutput->Add(fHistMCJetsPt);
	131
6fd5039f	132	fHistMCJetsPtNonBias = new TH1F("fHistMCJetsPtNonBias", "fHistMCJetsPtNonBias", fNbins, fMinBinPt, fMaxBinPt);
7549c451	133	fHistMCJetsPtNonBias->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	134	fHistMCJetsPtNonBias->GetYaxis()->SetTitle("counts");
	135	fOutput->Add(fHistMCJetsPtNonBias);
	136
6fd5039f	137	fHistMCJetsZvsPt = new TH2F("fHistMCJetsZvsPt", "fHistMCJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
7549c451	138	fHistMCJetsZvsPt->GetXaxis()->SetTitle("Z");
	139	fHistMCJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
	140	fOutput->Add(fHistMCJetsZvsPt);
	141
	142	if (fAnaType == kEMCAL) {
6fd5039f	143	fHistMCJetsNEFvsPt = new TH2F("fHistMCJetsNEFvsPt", "fHistMCJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
7549c451	144	fHistMCJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
	145	fHistMCJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
	146	fOutput->Add(fHistMCJetsNEFvsPt);
	147	}
	148
1b76c28f	149	fHistClosestDistance = new TH1F("fHistClosestDistance", "fHistClosestDistance", 50, 0, fMaxDistance * 1.2);
	150	fHistClosestDistance->GetXaxis()->SetTitle("d");
	151	fHistClosestDistance->GetYaxis()->SetTitle("counts");
	152	fOutput->Add(fHistClosestDistance);
	153
	154	fHistClosestDeltaPhi = new TH1F("fHistClosestDeltaPhi", "fHistClosestDeltaPhi", 64, -1.6, 4.8);
	155	fHistClosestDeltaPhi->GetXaxis()->SetTitle("#Delta#phi");
	156	fHistClosestDeltaPhi->GetYaxis()->SetTitle("counts");
	157	fOutput->Add(fHistClosestDeltaPhi);
	158
	159	fHistClosestDeltaEta = new TH1F("fHistClosestDeltaEta", "fHistClosestDeltaEta", TMath::CeilNint(fMaxEta - fMinEta) * 20, fMinEta * 2, fMaxEta * 2);
	160	fHistClosestDeltaEta->GetXaxis()->SetTitle("#Delta#eta");
	161	fHistClosestDeltaEta->GetYaxis()->SetTitle("counts");
	162	fOutput->Add(fHistClosestDeltaEta);
	163
6fd5039f	164	fHistClosestDeltaPt = new TH1F("fHistClosestDeltaPt", "fHistClosestDeltaPt", fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
1b76c28f	165	fHistClosestDeltaPt->GetXaxis()->SetTitle("#Delta p_{T}");
	166	fHistClosestDeltaPt->GetYaxis()->SetTitle("counts");
	167	fOutput->Add(fHistClosestDeltaPt);
	168
6fd5039f	169	fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
1b76c28f	170	fHistPartvsDetecPt->GetXaxis()->SetTitle("p_{T}^{det}");
	171	fHistPartvsDetecPt->GetYaxis()->SetTitle("p_{T}^{rec}");
	172	fOutput->Add(fHistPartvsDetecPt);
	173
7549c451	174	PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
2949a2ec	175	}
	176
	177	//________________________________________________________________________
6fd5039f	178	Bool_t AliJetResponseMaker::FillHistograms()
2949a2ec	179	{
	180	// Fill histograms.
	181
1b76c28f	182	// Find the closest jets
	183	DoJetLoop(fJets, fMCJets);
	184	DoJetLoop(fMCJets, fJets);
2949a2ec	185
1c59613f	186	const Int_t nJets = fJets->GetEntriesFast();
2949a2ec	187
1b76c28f	188	for (Int_t i = 0; i < nJets; i++) {
2949a2ec	189
1b76c28f	190	AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJets->At(i));
2949a2ec	191
2949a2ec	192	if (!jet) {
1b76c28f	193	AliError(Form("Could not receive jet %d", i));
2949a2ec	194	continue;
	195	}
	196
	197	if (!AcceptJet(jet))
	198	continue;
	199
	200	fHistJetsPtNonBias->Fill(jet->Pt());
	201
2949a2ec	202	if (jet->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC \|\| jet->MaxClusterPt() < fPtBiasJetClus))
	203	continue;
	204
1c59613f	205	if (jet->ClosestJet() && jet->ClosestJet()->ClosestJet() == jet &&
1c59613f	206	jet->ClosestJetDistance() < fMaxDistance) { // Matched jet found
1b76c28f	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
1c59613f	242	const Int_t nMCJets = fMCJets->GetEntriesFast();
1b76c28f	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	}
6fd5039f	274
6fd5039f	275	return kTRUE;
1b76c28f	276	}
	277
	278	//________________________________________________________________________
	279	void AliJetResponseMaker::DoJetLoop(TClonesArray jets1, TClonesArray jets2)
	280	{
	281	// Do the jet loop.
	282
	283	Int_t nJets1 = jets1->GetEntriesFast();
	284	Int_t nJets2 = jets2->GetEntriesFast();
	285
	286	for (Int_t i = 0; i < nJets1; i++) {
	287
	288	AliEmcalJet* jet1 = dynamic_cast<AliEmcalJet*>(jets1->At(i));
	289
	290	if (!jet1) {
	291	AliError(Form("Could not receive jet %d", i));
	292	continue;
	293	}
	294
	295	if (!AcceptJet(jet1))
	296	continue;
	297
1c59613f	298	if (jet1->MaxTrackPt() < fPtBiasJetTrack &&
1c59613f	299	(fAnaType == kTPC \|\| jet1->IsMC() \|\| jet1->MaxClusterPt() < fPtBiasJetClus))
1b76c28f	300	continue;
	301
	302	for (Int_t j = 0; j < nJets2; j++) {
	303
	304	AliEmcalJet* jet2 = dynamic_cast<AliEmcalJet*>(jets2->At(j));
	305
	306	if (!jet2) {
	307	AliError(Form("Could not receive jet %d", j));
	308	continue;
	309	}
	310
	311	if (!AcceptJet(jet2))
	312	continue;
	313
1c59613f	314	if (jet2->MaxTrackPt() < fPtBiasJetTrack &&
1c59613f	315	(fAnaType == kTPC \|\| jet2->IsMC() \|\| jet2->MaxClusterPt() < fPtBiasJetClus))
1b76c28f	316	continue;
	317
	318	Double_t deta = jet2->Eta() - jet1->Eta();
	319	Double_t dphi = jet2->Phi() - jet1->Phi();
	320	Double_t d = TMath::Sqrt(deta * deta + dphi * dphi);
	321
	322	if (d < jet1->ClosestJetDistance()) {
	323	jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
	324	jet1->SetClosestJet(jet2, d);
	325	}
	326	else if (d < jet1->SecondClosestJetDistance()) {
	327	jet1->SetSecondClosestJet(jet2, d);
2949a2ec	328	}
	329	}
	330	}
	331	}
	332
	333	//________________________________________________________________________
6fd5039f	334	Bool_t AliJetResponseMaker::RetrieveEventObjects()
2949a2ec	335	{
	336	// Retrieve event objects.
	337
1c59613f	338	if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
6fd5039f	339	return kFALSE;
2949a2ec	340
7549c451	341	if (!fMCJetsName.IsNull()) {
7549c451	342	fMCJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCJetsName));
2949a2ec	343	if (!fMCJets) {
	344	AliWarning(Form("Could not retrieve MC jets %s!", fMCJetsName.Data()));
	345	}
	346	}
7549c451	347
	348	if (!fMCTracksName.IsNull()) {
	349	fMCTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCTracksName));
	350	if (!fMCJets) {
	351	AliWarning(Form("Could not retrieve MC tracks %s!", fMCTracksName.Data()));
	352	}
	353	}
6fd5039f	354
6fd5039f	355	return kTRUE;
2949a2ec	356	}
	357
	358	//________________________________________________________________________
	359	void AliJetResponseMaker::Terminate(Option_t *)
	360	{
	361	// Called once at the end of the analysis.
	362	}