[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)
{
  // 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)
{
  // 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);

  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 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;

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

 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, kTRUE, kFALSE))
      continue;

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

  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),
1b76c28f	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"),
99c67c1b	58	fMaxDistance(0.25),
7549c451	59	fMCTracks(0),
2949a2ec	60	fMCJets(0),
	61	fHistMCJetPhiEta(0),
	62	fHistMCJetsPt(0),
2949a2ec	63	fHistMCJetsNEFvsPt(0),
	64	fHistMCJetsZvsPt(0),
	65	fHistJetPhiEta(0),
	66	fHistJetsPt(0),
2949a2ec	67	fHistJetsNEFvsPt(0),
1b76c28f	68	fHistJetsZvsPt(0),
	69	fHistClosestDistance(0),
	70	fHistClosestDeltaPhi(0),
	71	fHistClosestDeltaEta(0),
	72	fHistClosestDeltaPt(0),
99c67c1b	73	fHistNonMatchedMCJetPt(0),
99c67c1b	74	fHistNonMatchedJetPt(0),
1b76c28f	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	fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
2949a2ec	106	fHistJetsZvsPt->GetXaxis()->SetTitle("Z");
	107	fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
	108	fOutput->Add(fHistJetsZvsPt);
	109
	110	if (fAnaType == kEMCAL) {
6fd5039f	111	fHistJetsNEFvsPt = new TH2F("fHistJetsNEFvsPt", "fHistJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
2949a2ec	112	fHistJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
	113	fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
	114	fOutput->Add(fHistJetsNEFvsPt);
	115	}
7549c451	116
	117	fHistMCJetPhiEta = new TH2F("fHistMCJetPhiEta", "fHistMCJetPhiEta", 20, -2, 2, 32, 0, 6.4);
	118	fHistMCJetPhiEta->GetXaxis()->SetTitle("#eta");
	119	fHistMCJetPhiEta->GetYaxis()->SetTitle("#phi");
	120	fOutput->Add(fHistMCJetPhiEta);
	121
6fd5039f	122	fHistMCJetsPt = new TH1F("fHistMCJetsPt", "fHistMCJetsPt", fNbins, fMinBinPt, fMaxBinPt);
7549c451	123	fHistMCJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	124	fHistMCJetsPt->GetYaxis()->SetTitle("counts");
	125	fOutput->Add(fHistMCJetsPt);
	126
6fd5039f	127	fHistMCJetsZvsPt = new TH2F("fHistMCJetsZvsPt", "fHistMCJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
7549c451	128	fHistMCJetsZvsPt->GetXaxis()->SetTitle("Z");
	129	fHistMCJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
	130	fOutput->Add(fHistMCJetsZvsPt);
	131
	132	if (fAnaType == kEMCAL) {
6fd5039f	133	fHistMCJetsNEFvsPt = new TH2F("fHistMCJetsNEFvsPt", "fHistMCJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
7549c451	134	fHistMCJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
	135	fHistMCJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
	136	fOutput->Add(fHistMCJetsNEFvsPt);
	137	}
	138
1b76c28f	139	fHistClosestDistance = new TH1F("fHistClosestDistance", "fHistClosestDistance", 50, 0, fMaxDistance * 1.2);
	140	fHistClosestDistance->GetXaxis()->SetTitle("d");
	141	fHistClosestDistance->GetYaxis()->SetTitle("counts");
	142	fOutput->Add(fHistClosestDistance);
	143
	144	fHistClosestDeltaPhi = new TH1F("fHistClosestDeltaPhi", "fHistClosestDeltaPhi", 64, -1.6, 4.8);
	145	fHistClosestDeltaPhi->GetXaxis()->SetTitle("#Delta#phi");
	146	fHistClosestDeltaPhi->GetYaxis()->SetTitle("counts");
	147	fOutput->Add(fHistClosestDeltaPhi);
	148
	149	fHistClosestDeltaEta = new TH1F("fHistClosestDeltaEta", "fHistClosestDeltaEta", TMath::CeilNint(fMaxEta - fMinEta) * 20, fMinEta * 2, fMaxEta * 2);
	150	fHistClosestDeltaEta->GetXaxis()->SetTitle("#Delta#eta");
	151	fHistClosestDeltaEta->GetYaxis()->SetTitle("counts");
	152	fOutput->Add(fHistClosestDeltaEta);
	153
6fd5039f	154	fHistClosestDeltaPt = new TH1F("fHistClosestDeltaPt", "fHistClosestDeltaPt", fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
1b76c28f	155	fHistClosestDeltaPt->GetXaxis()->SetTitle("#Delta p_{T}");
	156	fHistClosestDeltaPt->GetYaxis()->SetTitle("counts");
	157	fOutput->Add(fHistClosestDeltaPt);
	158
99c67c1b	159	fHistNonMatchedMCJetPt = new TH1F("fHistNonMatchedMCJetPt", "fHistNonMatchedMCJetPt", fNbins, fMinBinPt, fMaxBinPt);
	160	fHistNonMatchedMCJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	161	fHistNonMatchedMCJetPt->GetYaxis()->SetTitle("counts");
	162	fOutput->Add(fHistNonMatchedMCJetPt);
	163
	164	fHistNonMatchedJetPt = new TH1F("fHistNonMatchedJetPt", "fHistNonMatchedJetPt", fNbins, fMinBinPt, fMaxBinPt);
	165	fHistNonMatchedJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	166	fHistNonMatchedJetPt->GetYaxis()->SetTitle("counts");
	167	fOutput->Add(fHistNonMatchedJetPt);
	168
6fd5039f	169	fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
99c67c1b	170	fHistPartvsDetecPt->GetXaxis()->SetTitle("p_{T}^{rec}");
99c67c1b	171	fHistPartvsDetecPt->GetYaxis()->SetTitle("p_{T}^{gen}");
1b76c28f	172	fOutput->Add(fHistPartvsDetecPt);
1b76c28f	173
99c67c1b	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
99c67c1b	183	DoJetLoop(fJets, fMCJets, kFALSE);
99c67c1b	184	DoJetLoop(fMCJets, fJets, kTRUE);
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
1c59613f	200	if (jet->ClosestJet() && jet->ClosestJet()->ClosestJet() == jet &&
1c59613f	201	jet->ClosestJetDistance() < fMaxDistance) { // Matched jet found
1b76c28f	202	jet->SetMatchedToClosest();
	203	jet->ClosestJet()->SetMatchedToClosest();
	204	fHistClosestDistance->Fill(jet->ClosestJetDistance());
	205	Double_t deta = jet->Eta() - jet->MatchedJet()->Eta();
	206	fHistClosestDeltaEta->Fill(deta);
	207	Double_t dphi = jet->Phi() - jet->MatchedJet()->Phi();
	208	fHistClosestDeltaPhi->Fill(dphi);
	209	Double_t dpt = jet->Pt() - jet->MatchedJet()->Pt();
	210	fHistClosestDeltaPt->Fill(dpt);
	211	fHistPartvsDetecPt->Fill(jet->Pt(), jet->MatchedJet()->Pt());
2949a2ec	212	}
99c67c1b	213	else {
	214	fHistNonMatchedJetPt->Fill(jet->Pt());
	215	}
2949a2ec	216
	217	fHistJetsPt->Fill(jet->Pt());
	218
	219	fHistJetPhiEta->Fill(jet->Eta(), jet->Phi());
	220
	221	if (fAnaType == kEMCAL)
	222	fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt());
	223
	224	for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
1b76c28f	225	AliVParticle *track = jet->TrackAt(it, fTracks);
2949a2ec	226	if (track)
	227	fHistJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt());
	228	}
	229
1b76c28f	230	for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
	231	AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
	232	if (cluster) {
	233	TLorentzVector nP;
	234	cluster->GetMomentum(nP, fVertex);
	235	fHistJetsZvsPt->Fill(nP.Pt() / jet->Pt(), jet->Pt());
	236	}
	237	}
	238	}
	239
1c59613f	240	const Int_t nMCJets = fMCJets->GetEntriesFast();
1b76c28f	241
	242	for (Int_t i = 0; i < nMCJets; i++) {
	243
	244	AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fMCJets->At(i));
	245
	246	if (!jet) {
	247	AliError(Form("Could not receive mc jet %d", i));
2949a2ec	248	continue;
1b76c28f	249	}
99c67c1b	250
99c67c1b	251	if (!AcceptJet(jet, kTRUE, kFALSE))
1b76c28f	252	continue;
2949a2ec	253
99c67c1b	254	if (!jet->MatchedJet()) {
	255	fHistNonMatchedMCJetPt->Fill(jet->Pt());
	256	}
1b76c28f	257
	258	fHistMCJetsPt->Fill(jet->Pt());
	259
	260	fHistMCJetPhiEta->Fill(jet->Eta(), jet->Phi());
	261
	262	if (fAnaType == kEMCAL)
	263	fHistMCJetsNEFvsPt->Fill(jet->NEF(), jet->Pt());
	264
	265	for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
	266	AliVParticle *track = jet->TrackAt(it, fMCTracks);
	267	if (track)
	268	fHistMCJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt());
	269	}
	270	}
6fd5039f	271
6fd5039f	272	return kTRUE;
1b76c28f	273	}
	274
	275	//________________________________________________________________________
99c67c1b	276	void AliJetResponseMaker::DoJetLoop(TClonesArray jets1, TClonesArray jets2, Bool_t mc)
1b76c28f	277	{
	278	// Do the jet loop.
	279
	280	Int_t nJets1 = jets1->GetEntriesFast();
	281	Int_t nJets2 = jets2->GetEntriesFast();
	282
	283	for (Int_t i = 0; i < nJets1; i++) {
	284
	285	AliEmcalJet* jet1 = dynamic_cast<AliEmcalJet*>(jets1->At(i));
	286
	287	if (!jet1) {
	288	AliError(Form("Could not receive jet %d", i));
	289	continue;
	290	}
	291
99c67c1b	292	if (!AcceptJet(jet1, kTRUE, mc))
1b76c28f	293	continue;
1b76c28f	294
1c59613f	295	if (jet1->MaxTrackPt() < fPtBiasJetTrack &&
1c59613f	296	(fAnaType == kTPC \|\| jet1->IsMC() \|\| jet1->MaxClusterPt() < fPtBiasJetClus))
1b76c28f	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
99c67c1b	308	if (!AcceptJet(jet2, kTRUE, !mc))
1b76c28f	309	continue;
1b76c28f	310
1c59613f	311	if (jet2->MaxTrackPt() < fPtBiasJetTrack &&
1c59613f	312	(fAnaType == kTPC \|\| jet2->IsMC() \|\| jet2->MaxClusterPt() < fPtBiasJetClus))
1b76c28f	313	continue;
	314
	315	Double_t deta = jet2->Eta() - jet1->Eta();
	316	Double_t dphi = jet2->Phi() - jet1->Phi();
	317	Double_t d = TMath::Sqrt(deta * deta + dphi * dphi);
	318
	319	if (d < jet1->ClosestJetDistance()) {
	320	jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
	321	jet1->SetClosestJet(jet2, d);
	322	}
	323	else if (d < jet1->SecondClosestJetDistance()) {
	324	jet1->SetSecondClosestJet(jet2, d);
2949a2ec	325	}
	326	}
	327	}
	328	}
	329
	330	//________________________________________________________________________
6fd5039f	331	Bool_t AliJetResponseMaker::RetrieveEventObjects()
2949a2ec	332	{
	333	// Retrieve event objects.
	334
1c59613f	335	if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
6fd5039f	336	return kFALSE;
2949a2ec	337
7549c451	338	if (!fMCJetsName.IsNull()) {
7549c451	339	fMCJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCJetsName));
2949a2ec	340	if (!fMCJets) {
	341	AliWarning(Form("Could not retrieve MC jets %s!", fMCJetsName.Data()));
	342	}
	343	}
7549c451	344
	345	if (!fMCTracksName.IsNull()) {
	346	fMCTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCTracksName));
	347	if (!fMCJets) {
	348	AliWarning(Form("Could not retrieve MC tracks %s!", fMCTracksName.Data()));
	349	}
	350	}
6fd5039f	351
6fd5039f	352	return kTRUE;
2949a2ec	353	}
	354
	355	//________________________________________________________________________
	356	void AliJetResponseMaker::Terminate(Option_t *)
	357	{
	358	// Called once at the end of the analysis.
	359	}