[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.25),
  fMCTracks(0),
  fMCJets(0),
  fHistMCJetPhiEta(0),
  fHistMCJetsPt(0),
  fHistMCJetsNEFvsPt(0),
  fHistMCJetsZvsPt(0),
  fHistJetPhiEta(0),
  fHistJetsPt(0),
  fHistJetsNEFvsPt(0),
  fHistJetsZvsPt(0),
  fHistClosestDistance(0),
  fHistClosestDeltaPhi(0),
  fHistClosestDeltaEta(0),
  fHistClosestDeltaPt(0),
  fHistNonMatchedMCJetPt(0),
  fHistNonMatchedJetPt(0),
  fHistPartvsDetecPt(0),
  fHistMissedMCJets(0)
{
  // Default constructor.
}

//________________________________________________________________________
AliJetResponseMaker::AliJetResponseMaker(const char *name) : 
  AliAnalysisTaskEmcalJet(name),
  fMCTracksName("MCParticles"),
  fMCJetsName("MCJets"),
  fMaxDistance(0.25),
  fMCTracks(0),
  fMCJets(0),
  fHistMCJetPhiEta(0),
  fHistMCJetsPt(0),
  fHistMCJetsNEFvsPt(0),
  fHistMCJetsZvsPt(0),
  fHistJetPhiEta(0),
  fHistJetsPt(0),
  fHistJetsNEFvsPt(0),
  fHistJetsZvsPt(0),
  fHistClosestDistance(0),
  fHistClosestDeltaPhi(0),
  fHistClosestDeltaEta(0),
  fHistClosestDeltaPt(0),
  fHistNonMatchedMCJetPt(0),
  fHistNonMatchedJetPt(0),
  fHistPartvsDetecPt(0),
  fHistMissedMCJets(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);
  
  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);
  
  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);

  fHistNonMatchedMCJetPt = new TH1F("fHistNonMatchedMCJetPt", "fHistNonMatchedMCJetPt", fNbins, fMinBinPt, fMaxBinPt);
  fHistNonMatchedMCJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  fHistNonMatchedMCJetPt->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistNonMatchedMCJetPt);

  fHistNonMatchedJetPt = new TH1F("fHistNonMatchedJetPt", "fHistNonMatchedJetPt", fNbins, fMinBinPt, fMaxBinPt);
  fHistNonMatchedJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  fHistNonMatchedJetPt->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistNonMatchedJetPt);

  fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
  fHistPartvsDetecPt->GetXaxis()->SetTitle("p_{T}^{rec}");
  fHistPartvsDetecPt->GetYaxis()->SetTitle("p_{T}^{gen}");
  fOutput->Add(fHistPartvsDetecPt);

  fHistMissedMCJets = new TH1F("fHistMissedMCJets", "fHistMissedMCJets", fNbins, fMinBinPt, fMaxBinPt);
  fHistMissedMCJets->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  fHistMissedMCJets->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistMissedMCJets);

  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, kFALSE);
  DoJetLoop(fMCJets, fJets, kTRUE);

  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 jet %d", i));
      continue;
    }  

    if (!AcceptJet(jet, kTRUE, kFALSE))
      continue;

    if (jet->Pt() > fMaxBinPt)
      continue;

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

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

  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 mc jet %d", i));
      continue;
    }  
    
    if (!AcceptJet(jet))
      continue;

    if (!jet->MatchedJet()) {
      fHistNonMatchedJetPt->Fill(jet->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());
      }
    }
  }

  return kTRUE;
}

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