[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 "AliGenPythiaEventHeader.h"
#include "AliMCEvent.h"

ClassImp(AliJetResponseMaker)

//________________________________________________________________________
AliJetResponseMaker::AliJetResponseMaker() : 
  AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
  fMCTracksName("MCParticles"),
  fMCJetsName("MCJets"),
  fMaxDistance(0.25),
  fDoWeighting(kFALSE),
  fVertexCut(10),
  fPythiaHeader(0),
  fEventWeight(0),
  fPtHardBin(0),
  fNTrials(0),
  fMCTracks(0),
  fMCJets(0),
  fHistNTrials(0),
  fHistAcceptedEvents(0),
  fHistEvents(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, kTRUE),
  fMCTracksName("MCParticles"),
  fMCJetsName("MCJets"),
  fMaxDistance(0.25),
  fDoWeighting(kFALSE),
  fVertexCut(10),
  fPythiaHeader(0),
  fEventWeight(0),
  fPtHardBin(0),
  fNTrials(0),
  fMCTracks(0),
  fMCJets(0),
  fHistNTrials(0),
  fHistAcceptedEvents(0),
  fHistEvents(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();

  const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
  const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};

  fHistNTrials = new TH1F("fHistNTrials", "fHistNTrials", 11, 0, 11);
  fHistNTrials->GetXaxis()->SetTitle("p_{T} hard bin");
  fHistNTrials->GetYaxis()->SetTitle("trials");
  fOutput->Add(fHistNTrials);

  fHistAcceptedEvents = new TH1F("fHistAcceptedEvents", "fHistAcceptedEvents", 11, 0, 11);
  fHistAcceptedEvents->GetXaxis()->SetTitle("p_{T} hard bin");
  fHistAcceptedEvents->GetYaxis()->SetTitle("accepted events");
  fOutput->Add(fHistAcceptedEvents);

  fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
  fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
  fHistEvents->GetYaxis()->SetTitle("total events");
  fOutput->Add(fHistEvents);

  for (Int_t i = 1; i < 12; i++) {
    fHistNTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
    fHistAcceptedEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
    fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
  }

  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
}

//________________________________________________________________________
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 = static_cast<AliEmcalJet*>(jets1->At(i));

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

    if (!AcceptJet(jet1, kTRUE, mc))
      continue;

    for (Int_t j = 0; j < nJets2; j++) {
      
      AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
      
      if (!jet2) {
	AliError(Form("Could not receive jet %d", j));
	continue;
      }  
      
      if (!AcceptJet(jet2, kTRUE, !mc))
	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::ExecOnce()
{
  // Execute once.

  if (!fMCJetsName.IsNull() && !fMCJets) {
    fMCJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCJetsName));
    if (!fMCJets) {
      AliError(Form("%s: Could not retrieve mc jets %s!", GetName(), fMCJetsName.Data()));
      return;
    }
    else if (!fMCJets->GetClass()->GetBaseClass("AliEmcalJet")) {
      AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fMCJetsName.Data())); 
      fMCJets = 0;
      return;
    }
  }

  if (!fMCTracksName.IsNull() && !fMCTracks) {
    fMCTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCTracksName));
    if (!fMCTracks) {
      AliError(Form("%s: Could not retrieve mc tracks %s!", GetName(), fMCTracksName.Data())); 
      return;
    }
    else {
      TClass *cl = fMCTracks->GetClass();
      if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
	AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fMCTracksName.Data())); 
	fMCTracks = 0;
	return;
      }
    }
  }

  AliAnalysisTaskEmcalJet::ExecOnce();
}

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

  const Int_t nMCJets = fMCJets->GetEntriesFast();

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

    AliEmcalJet* jet = static_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(), fEventWeight);
      }
      else {
	fHistClosestDistance->Fill(jet->ClosestJetDistance(), fEventWeight);

	Double_t deta = jet->MatchedJet()->Eta() - jet->Eta();
	fHistClosestDeltaEta->Fill(deta, fEventWeight);

	Double_t dphi = jet->MatchedJet()->Phi() - jet->Phi();
	fHistClosestDeltaPhi->Fill(dphi, fEventWeight);

	Double_t dpt = jet->MatchedJet()->Pt() - jet->Pt();
	fHistClosestDeltaPt->Fill(dpt, fEventWeight);

	fHistPartvsDetecPt->Fill(jet->MatchedJet()->Pt(), jet->Pt(), fEventWeight);
      }
    }
    else {
      fHistNonMatchedMCJetPt->Fill(jet->Pt(), fEventWeight);
      fHistMissedMCJets->Fill(jet->Pt(), fEventWeight);
    }

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

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

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

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

  const Int_t nJets = fJets->GetEntriesFast();

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

    AliEmcalJet* jet = static_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(), fEventWeight);
    }

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

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

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

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

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

  return kTRUE;
}

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

  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
    return kFALSE;
  
  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());

  if (!fPythiaHeader)
    return kFALSE;

  if (fDoWeighting)
    fEventWeight = fPythiaHeader->EventWeight();
  else
    fEventWeight = 1;

  const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
  const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};

  Double_t pthard = fPythiaHeader->GetPtHard();
  
  for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
    if (ptHardLo[fPtHardBin] < pthard && ptHardHi[fPtHardBin] > pthard)
      break;
  }

  fNTrials = fPythiaHeader->Trials();

  return kTRUE;
}

//________________________________________________________________________
Bool_t AliJetResponseMaker::Run()
{
  // Find the closest jets

  fHistEvents->SetBinContent(fPtHardBin, fHistEvents->GetBinContent(fPtHardBin) + 1);

  if (TMath::Abs(fVertex[2]) > fVertexCut)
    return kFALSE;
  
  fHistAcceptedEvents->SetBinContent(fPtHardBin, fHistEvents->GetBinContent(fPtHardBin) + 1);
  fHistNTrials->SetBinContent(fPtHardBin, fHistNTrials->GetBinContent(fPtHardBin) + fNTrials);

  DoJetLoop(fJets, fMCJets, kFALSE);
  DoJetLoop(fMCJets, fJets, kTRUE);

  return kTRUE;
}
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"
1ee1b5b8	22	#include "AliGenPythiaEventHeader.h"
2949a2ec	23	#include "AliMCEvent.h"
	24
	25	ClassImp(AliJetResponseMaker)
	26
	27	//________________________________________________________________________
	28	AliJetResponseMaker::AliJetResponseMaker() :
e44e8726	29	AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
7549c451	30	fMCTracksName("MCParticles"),
7549c451	31	fMCJetsName("MCJets"),
99c67c1b	32	fMaxDistance(0.25),
1ee1b5b8	33	fDoWeighting(kFALSE),
	34	fVertexCut(10),
	35	fPythiaHeader(0),
	36	fEventWeight(0),
	37	fPtHardBin(0),
	38	fNTrials(0),
7549c451	39	fMCTracks(0),
2949a2ec	40	fMCJets(0),
1ee1b5b8	41	fHistNTrials(0),
	42	fHistAcceptedEvents(0),
	43	fHistEvents(0),
2949a2ec	44	fHistMCJetPhiEta(0),
2949a2ec	45	fHistMCJetsPt(0),
2949a2ec	46	fHistMCJetsNEFvsPt(0),
	47	fHistMCJetsZvsPt(0),
	48	fHistJetPhiEta(0),
	49	fHistJetsPt(0),
2949a2ec	50	fHistJetsNEFvsPt(0),
1b76c28f	51	fHistJetsZvsPt(0),
	52	fHistClosestDistance(0),
	53	fHistClosestDeltaPhi(0),
	54	fHistClosestDeltaEta(0),
	55	fHistClosestDeltaPt(0),
99c67c1b	56	fHistNonMatchedMCJetPt(0),
99c67c1b	57	fHistNonMatchedJetPt(0),
2bddb6ae	58	fHistPartvsDetecPt(0),
2bddb6ae	59	fHistMissedMCJets(0)
2949a2ec	60	{
	61	// Default constructor.
	62	}
	63
	64	//________________________________________________________________________
	65	AliJetResponseMaker::AliJetResponseMaker(const char *name) :
e44e8726	66	AliAnalysisTaskEmcalJet(name, kTRUE),
7549c451	67	fMCTracksName("MCParticles"),
7549c451	68	fMCJetsName("MCJets"),
99c67c1b	69	fMaxDistance(0.25),
1ee1b5b8	70	fDoWeighting(kFALSE),
	71	fVertexCut(10),
	72	fPythiaHeader(0),
	73	fEventWeight(0),
	74	fPtHardBin(0),
	75	fNTrials(0),
7549c451	76	fMCTracks(0),
2949a2ec	77	fMCJets(0),
1ee1b5b8	78	fHistNTrials(0),
	79	fHistAcceptedEvents(0),
	80	fHistEvents(0),
2949a2ec	81	fHistMCJetPhiEta(0),
2949a2ec	82	fHistMCJetsPt(0),
2949a2ec	83	fHistMCJetsNEFvsPt(0),
	84	fHistMCJetsZvsPt(0),
	85	fHistJetPhiEta(0),
	86	fHistJetsPt(0),
2949a2ec	87	fHistJetsNEFvsPt(0),
1b76c28f	88	fHistJetsZvsPt(0),
	89	fHistClosestDistance(0),
	90	fHistClosestDeltaPhi(0),
	91	fHistClosestDeltaEta(0),
	92	fHistClosestDeltaPt(0),
99c67c1b	93	fHistNonMatchedMCJetPt(0),
99c67c1b	94	fHistNonMatchedJetPt(0),
2bddb6ae	95	fHistPartvsDetecPt(0),
2bddb6ae	96	fHistMissedMCJets(0)
2949a2ec	97	{
2949a2ec	98	// Standard constructor.
2949a2ec	99	}
	100
	101	//________________________________________________________________________
	102	AliJetResponseMaker::~AliJetResponseMaker()
	103	{
	104	// Destructor
	105	}
	106
	107	//________________________________________________________________________
	108	void AliJetResponseMaker::UserCreateOutputObjects()
	109	{
	110	// Create user objects.
	111
	112	OpenFile(1);
	113	fOutput = new TList();
	114	fOutput->SetOwner();
	115
1ee1b5b8	116	const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
	117	const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
	118
	119	fHistNTrials = new TH1F("fHistNTrials", "fHistNTrials", 11, 0, 11);
	120	fHistNTrials->GetXaxis()->SetTitle("p_{T} hard bin");
	121	fHistNTrials->GetYaxis()->SetTitle("trials");
	122	fOutput->Add(fHistNTrials);
	123
	124	fHistAcceptedEvents = new TH1F("fHistAcceptedEvents", "fHistAcceptedEvents", 11, 0, 11);
	125	fHistAcceptedEvents->GetXaxis()->SetTitle("p_{T} hard bin");
	126	fHistAcceptedEvents->GetYaxis()->SetTitle("accepted events");
	127	fOutput->Add(fHistAcceptedEvents);
	128
	129	fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
	130	fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
	131	fHistEvents->GetYaxis()->SetTitle("total events");
	132	fOutput->Add(fHistEvents);
	133
	134	for (Int_t i = 1; i < 12; i++) {
	135	fHistNTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
	136	fHistAcceptedEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
	137	fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
	138	}
	139
7549c451	140	fHistJetPhiEta = new TH2F("fHistJetPhiEta", "fHistJetPhiEta", 20, -2, 2, 32, 0, 6.4);
	141	fHistJetPhiEta->GetXaxis()->SetTitle("#eta");
	142	fHistJetPhiEta->GetYaxis()->SetTitle("#phi");
	143	fOutput->Add(fHistJetPhiEta);
	144
6fd5039f	145	fHistJetsPt = new TH1F("fHistJetsPt", "fHistJetsPt", fNbins, fMinBinPt, fMaxBinPt);
2949a2ec	146	fHistJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	147	fHistJetsPt->GetYaxis()->SetTitle("counts");
	148	fOutput->Add(fHistJetsPt);
	149
6fd5039f	150	fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
2949a2ec	151	fHistJetsZvsPt->GetXaxis()->SetTitle("Z");
	152	fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
	153	fOutput->Add(fHistJetsZvsPt);
	154
	155	if (fAnaType == kEMCAL) {
6fd5039f	156	fHistJetsNEFvsPt = new TH2F("fHistJetsNEFvsPt", "fHistJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
2949a2ec	157	fHistJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
	158	fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
	159	fOutput->Add(fHistJetsNEFvsPt);
	160	}
7549c451	161
	162	fHistMCJetPhiEta = new TH2F("fHistMCJetPhiEta", "fHistMCJetPhiEta", 20, -2, 2, 32, 0, 6.4);
	163	fHistMCJetPhiEta->GetXaxis()->SetTitle("#eta");
	164	fHistMCJetPhiEta->GetYaxis()->SetTitle("#phi");
	165	fOutput->Add(fHistMCJetPhiEta);
	166
6fd5039f	167	fHistMCJetsPt = new TH1F("fHistMCJetsPt", "fHistMCJetsPt", fNbins, fMinBinPt, fMaxBinPt);
7549c451	168	fHistMCJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	169	fHistMCJetsPt->GetYaxis()->SetTitle("counts");
	170	fOutput->Add(fHistMCJetsPt);
	171
6fd5039f	172	fHistMCJetsZvsPt = new TH2F("fHistMCJetsZvsPt", "fHistMCJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
7549c451	173	fHistMCJetsZvsPt->GetXaxis()->SetTitle("Z");
	174	fHistMCJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
	175	fOutput->Add(fHistMCJetsZvsPt);
	176
	177	if (fAnaType == kEMCAL) {
6fd5039f	178	fHistMCJetsNEFvsPt = new TH2F("fHistMCJetsNEFvsPt", "fHistMCJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
7549c451	179	fHistMCJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
	180	fHistMCJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
	181	fOutput->Add(fHistMCJetsNEFvsPt);
	182	}
	183
1b76c28f	184	fHistClosestDistance = new TH1F("fHistClosestDistance", "fHistClosestDistance", 50, 0, fMaxDistance * 1.2);
	185	fHistClosestDistance->GetXaxis()->SetTitle("d");
	186	fHistClosestDistance->GetYaxis()->SetTitle("counts");
	187	fOutput->Add(fHistClosestDistance);
	188
	189	fHistClosestDeltaPhi = new TH1F("fHistClosestDeltaPhi", "fHistClosestDeltaPhi", 64, -1.6, 4.8);
	190	fHistClosestDeltaPhi->GetXaxis()->SetTitle("#Delta#phi");
	191	fHistClosestDeltaPhi->GetYaxis()->SetTitle("counts");
	192	fOutput->Add(fHistClosestDeltaPhi);
	193
	194	fHistClosestDeltaEta = new TH1F("fHistClosestDeltaEta", "fHistClosestDeltaEta", TMath::CeilNint(fMaxEta - fMinEta) * 20, fMinEta * 2, fMaxEta * 2);
	195	fHistClosestDeltaEta->GetXaxis()->SetTitle("#Delta#eta");
	196	fHistClosestDeltaEta->GetYaxis()->SetTitle("counts");
	197	fOutput->Add(fHistClosestDeltaEta);
	198
6fd5039f	199	fHistClosestDeltaPt = new TH1F("fHistClosestDeltaPt", "fHistClosestDeltaPt", fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
1b76c28f	200	fHistClosestDeltaPt->GetXaxis()->SetTitle("#Delta p_{T}");
	201	fHistClosestDeltaPt->GetYaxis()->SetTitle("counts");
	202	fOutput->Add(fHistClosestDeltaPt);
	203
99c67c1b	204	fHistNonMatchedMCJetPt = new TH1F("fHistNonMatchedMCJetPt", "fHistNonMatchedMCJetPt", fNbins, fMinBinPt, fMaxBinPt);
	205	fHistNonMatchedMCJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	206	fHistNonMatchedMCJetPt->GetYaxis()->SetTitle("counts");
	207	fOutput->Add(fHistNonMatchedMCJetPt);
	208
	209	fHistNonMatchedJetPt = new TH1F("fHistNonMatchedJetPt", "fHistNonMatchedJetPt", fNbins, fMinBinPt, fMaxBinPt);
	210	fHistNonMatchedJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	211	fHistNonMatchedJetPt->GetYaxis()->SetTitle("counts");
	212	fOutput->Add(fHistNonMatchedJetPt);
	213
6fd5039f	214	fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
99c67c1b	215	fHistPartvsDetecPt->GetXaxis()->SetTitle("p_{T}^{rec}");
99c67c1b	216	fHistPartvsDetecPt->GetYaxis()->SetTitle("p_{T}^{gen}");
1b76c28f	217	fOutput->Add(fHistPartvsDetecPt);
1b76c28f	218
2bddb6ae	219	fHistMissedMCJets = new TH1F("fHistMissedMCJets", "fHistMissedMCJets", fNbins, fMinBinPt, fMaxBinPt);
	220	fHistMissedMCJets->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	221	fHistMissedMCJets->GetYaxis()->SetTitle("counts");
	222	fOutput->Add(fHistMissedMCJets);
	223
99c67c1b	224	PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
2949a2ec	225	}
	226
	227	//________________________________________________________________________
44476be7	228	void AliJetResponseMaker::DoJetLoop(TClonesArray jets1, TClonesArray jets2, Bool_t mc)
2949a2ec	229	{
44476be7	230	// Do the jet loop.
1ee1b5b8	231
44476be7	232	Int_t nJets1 = jets1->GetEntriesFast();
44476be7	233	Int_t nJets2 = jets2->GetEntriesFast();
1ee1b5b8	234
44476be7	235	for (Int_t i = 0; i < nJets1; i++) {
1ee1b5b8	236
44476be7	237	AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
2949a2ec	238
44476be7	239	if (!jet1) {
	240	AliError(Form("Could not receive jet %d", i));
	241	continue;
	242	}
	243
	244	if (!AcceptJet(jet1, kTRUE, mc))
	245	continue;
	246
	247	for (Int_t j = 0; j < nJets2; j++) {
	248
	249	AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
	250
	251	if (!jet2) {
	252	AliError(Form("Could not receive jet %d", j));
	253	continue;
	254	}
	255
	256	if (!AcceptJet(jet2, kTRUE, !mc))
	257	continue;
	258
	259	Double_t deta = jet2->Eta() - jet1->Eta();
	260	Double_t dphi = jet2->Phi() - jet1->Phi();
	261	Double_t d = TMath::Sqrt(deta * deta + dphi * dphi);
	262
	263	if (d < jet1->ClosestJetDistance()) {
	264	jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
	265	jet1->SetClosestJet(jet2, d);
	266	}
	267	else if (d < jet1->SecondClosestJetDistance()) {
	268	jet1->SetSecondClosestJet(jet2, d);
	269	}
	270	}
	271	}
1ee1b5b8	272	}
1ee1b5b8	273
44476be7	274	//________________________________________________________________________
	275	void AliJetResponseMaker::ExecOnce()
	276	{
	277	// Execute once.
	278
	279	if (!fMCJetsName.IsNull() && !fMCJets) {
	280	fMCJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCJetsName));
	281	if (!fMCJets) {
	282	AliError(Form("%s: Could not retrieve mc jets %s!", GetName(), fMCJetsName.Data()));
	283	return;
	284	}
	285	else if (!fMCJets->GetClass()->GetBaseClass("AliEmcalJet")) {
	286	AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fMCJetsName.Data()));
	287	fMCJets = 0;
	288	return;
	289	}
	290	}
	291
	292	if (!fMCTracksName.IsNull() && !fMCTracks) {
	293	fMCTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCTracksName));
	294	if (!fMCTracks) {
	295	AliError(Form("%s: Could not retrieve mc tracks %s!", GetName(), fMCTracksName.Data()));
	296	return;
	297	}
	298	else {
	299	TClass *cl = fMCTracks->GetClass();
	300	if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
	301	AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fMCTracksName.Data()));
	302	fMCTracks = 0;
	303	return;
	304	}
	305	}
	306	}
	307
	308	AliAnalysisTaskEmcalJet::ExecOnce();
	309	}
1ee1b5b8	310
	311	//________________________________________________________________________
	312	Bool_t AliJetResponseMaker::FillHistograms()
	313	{
	314	// Fill histograms.
	315
2bddb6ae	316	const Int_t nMCJets = fMCJets->GetEntriesFast();
2949a2ec	317
2bddb6ae	318	for (Int_t i = 0; i < nMCJets; i++) {
2949a2ec	319
e44e8726	320	AliEmcalJet* jet = static_cast<AliEmcalJet*>(fMCJets->At(i));
2949a2ec	321
2949a2ec	322	if (!jet) {
1b76c28f	323	AliError(Form("Could not receive jet %d", i));
2949a2ec	324	continue;
	325	}
	326
2bddb6ae	327	if (!AcceptJet(jet, kTRUE, kFALSE))
	328	continue;
	329
	330	if (jet->Pt() > fMaxBinPt)
2949a2ec	331	continue;
2949a2ec	332
1c59613f	333	if (jet->ClosestJet() && jet->ClosestJet()->ClosestJet() == jet &&
1c59613f	334	jet->ClosestJetDistance() < fMaxDistance) { // Matched jet found
1b76c28f	335	jet->SetMatchedToClosest();
1b76c28f	336	jet->ClosestJet()->SetMatchedToClosest();
2bddb6ae	337	if (jet->MatchedJet()->Pt() > fMaxBinPt) {
1ee1b5b8	338	fHistMissedMCJets->Fill(jet->Pt(), fEventWeight);
2bddb6ae	339	}
2bddb6ae	340	else {
1ee1b5b8	341	fHistClosestDistance->Fill(jet->ClosestJetDistance(), fEventWeight);
1ee1b5b8	342
2bddb6ae	343	Double_t deta = jet->MatchedJet()->Eta() - jet->Eta();
1ee1b5b8	344	fHistClosestDeltaEta->Fill(deta, fEventWeight);
1ee1b5b8	345
2bddb6ae	346	Double_t dphi = jet->MatchedJet()->Phi() - jet->Phi();
1ee1b5b8	347	fHistClosestDeltaPhi->Fill(dphi, fEventWeight);
1ee1b5b8	348
2bddb6ae	349	Double_t dpt = jet->MatchedJet()->Pt() - jet->Pt();
1ee1b5b8	350	fHistClosestDeltaPt->Fill(dpt, fEventWeight);
	351
	352	fHistPartvsDetecPt->Fill(jet->MatchedJet()->Pt(), jet->Pt(), fEventWeight);
2bddb6ae	353	}
2949a2ec	354	}
99c67c1b	355	else {
1ee1b5b8	356	fHistNonMatchedMCJetPt->Fill(jet->Pt(), fEventWeight);
1ee1b5b8	357	fHistMissedMCJets->Fill(jet->Pt(), fEventWeight);
99c67c1b	358	}
2949a2ec	359
1ee1b5b8	360	fHistMCJetsPt->Fill(jet->Pt(), fEventWeight);
2949a2ec	361
1ee1b5b8	362	fHistMCJetPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
2949a2ec	363
2949a2ec	364	if (fAnaType == kEMCAL)
1ee1b5b8	365	fHistMCJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight);
2949a2ec	366
2949a2ec	367	for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
2bddb6ae	368	AliVParticle *track = jet->TrackAt(it, fMCTracks);
2949a2ec	369	if (track)
1ee1b5b8	370	fHistMCJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt(), fEventWeight);
1b76c28f	371	}
	372	}
	373
2bddb6ae	374	const Int_t nJets = fJets->GetEntriesFast();
1b76c28f	375
2bddb6ae	376	for (Int_t i = 0; i < nJets; i++) {
1b76c28f	377
e44e8726	378	AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
1b76c28f	379
	380	if (!jet) {
	381	AliError(Form("Could not receive mc jet %d", i));
2949a2ec	382	continue;
1b76c28f	383	}
99c67c1b	384
2bddb6ae	385	if (!AcceptJet(jet))
1b76c28f	386	continue;
2949a2ec	387
99c67c1b	388	if (!jet->MatchedJet()) {
1ee1b5b8	389	fHistNonMatchedJetPt->Fill(jet->Pt(), fEventWeight);
99c67c1b	390	}
1b76c28f	391
1ee1b5b8	392	fHistJetsPt->Fill(jet->Pt(), fEventWeight);
1b76c28f	393
1ee1b5b8	394	fHistJetPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
1b76c28f	395
1b76c28f	396	if (fAnaType == kEMCAL)
1ee1b5b8	397	fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight);
1b76c28f	398
1b76c28f	399	for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
2bddb6ae	400	AliVParticle *track = jet->TrackAt(it, fTracks);
1b76c28f	401	if (track)
1ee1b5b8	402	fHistJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt(), fEventWeight);
2bddb6ae	403	}
	404
	405	for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
	406	AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
	407	if (cluster) {
	408	TLorentzVector nP;
	409	cluster->GetMomentum(nP, fVertex);
1ee1b5b8	410	fHistJetsZvsPt->Fill(nP.Pt() / jet->Pt(), jet->Pt(), fEventWeight);
2bddb6ae	411	}
1b76c28f	412	}
1b76c28f	413	}
6fd5039f	414
6fd5039f	415	return kTRUE;
1b76c28f	416	}
1b76c28f	417
2949a2ec	418	//________________________________________________________________________
6fd5039f	419	Bool_t AliJetResponseMaker::RetrieveEventObjects()
2949a2ec	420	{
44476be7	421	// Retrieve event objects.
44476be7	422
1c59613f	423	if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
6fd5039f	424	return kFALSE;
2949a2ec	425
1ee1b5b8	426	fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
	427
	428	if (!fPythiaHeader)
	429	return kFALSE;
	430
	431	if (fDoWeighting)
	432	fEventWeight = fPythiaHeader->EventWeight();
	433	else
	434	fEventWeight = 1;
	435
	436	const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
	437	const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
	438
	439	Double_t pthard = fPythiaHeader->GetPtHard();
	440
	441	for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
	442	if (ptHardLo[fPtHardBin] < pthard && ptHardHi[fPtHardBin] > pthard)
	443	break;
	444	}
	445
	446	fNTrials = fPythiaHeader->Trials();
	447
	448	return kTRUE;
	449	}
	450
	451	//________________________________________________________________________
44476be7	452	Bool_t AliJetResponseMaker::Run()
1ee1b5b8	453	{
44476be7	454	// Find the closest jets
1ee1b5b8	455
44476be7	456	fHistEvents->SetBinContent(fPtHardBin, fHistEvents->GetBinContent(fPtHardBin) + 1);
7549c451	457
44476be7	458	if (TMath::Abs(fVertex[2]) > fVertexCut)
	459	return kFALSE;
	460
	461	fHistAcceptedEvents->SetBinContent(fPtHardBin, fHistEvents->GetBinContent(fPtHardBin) + 1);
94d22b76	462	fHistNTrials->SetBinContent(fPtHardBin, fHistNTrials->GetBinContent(fPtHardBin) + fNTrials);
6fd5039f	463
44476be7	464	DoJetLoop(fJets, fMCJets, kFALSE);
44476be7	465	DoJetLoop(fMCJets, fJets, kTRUE);
2949a2ec	466
44476be7	467	return kTRUE;
2949a2ec	468	}
44476be7	469