[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::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::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 (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin])
      break;
  }
  fPtHardBin--;

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

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


//________________________________________________________________________
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;
}
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	}
2949a2ec	226
b3ccdfe2	227	//________________________________________________________________________
	228	void AliJetResponseMaker::ExecOnce()
	229	{
	230	// Execute once.
	231
	232	if (!fMCJetsName.IsNull() && !fMCJets) {
	233	fMCJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCJetsName));
	234	if (!fMCJets) {
	235	AliError(Form("%s: Could not retrieve mc jets %s!", GetName(), fMCJetsName.Data()));
	236	return;
	237	}
	238	else if (!fMCJets->GetClass()->GetBaseClass("AliEmcalJet")) {
	239	AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fMCJetsName.Data()));
	240	fMCJets = 0;
	241	return;
	242	}
	243	}
	244
	245	if (!fMCTracksName.IsNull() && !fMCTracks) {
	246	fMCTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCTracksName));
	247	if (!fMCTracks) {
	248	AliError(Form("%s: Could not retrieve mc tracks %s!", GetName(), fMCTracksName.Data()));
	249	return;
	250	}
	251	else {
	252	TClass *cl = fMCTracks->GetClass();
	253	if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
	254	AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fMCTracksName.Data()));
	255	fMCTracks = 0;
	256	return;
	257	}
	258	}
	259	}
	260
	261	AliAnalysisTaskEmcalJet::ExecOnce();
	262	}
	263
	264	//________________________________________________________________________
	265	Bool_t AliJetResponseMaker::RetrieveEventObjects()
	266	{
	267	// Retrieve event objects.
	268
	269	if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
	270	return kFALSE;
	271
	272	fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
	273
	274	if (!fPythiaHeader)
	275	return kFALSE;
	276
	277	if (fDoWeighting)
	278	fEventWeight = fPythiaHeader->EventWeight();
	279	else
	280	fEventWeight = 1;
	281
	282	const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
	283	const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
	284
	285	Double_t pthard = fPythiaHeader->GetPtHard();
	286
	287	for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
	288	if (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin])
	289	break;
	290	}
291	fPtHardBin--;
292
293	fNTrials = fPythiaHeader->Trials();
294
295	return kTRUE;
296	}
297
298	//________________________________________________________________________
299	Bool_t AliJetResponseMaker::Run()
300	{
301	// Find the closest jets
302
303	fHistEvents->SetBinContent(fPtHardBin, fHistEvents->GetBinContent(fPtHardBin) + 1);
304
305	if (TMath::Abs(fVertex[2]) > fVertexCut)
306	return kFALSE;
307
308	fHistAcceptedEvents->SetBinContent(fPtHardBin, fHistEvents->GetBinContent(fPtHardBin) + 1);
309	fHistNTrials->SetBinContent(fPtHardBin, fHistNTrials->GetBinContent(fPtHardBin) + fNTrials);
310
311	DoJetLoop(fJets, fMCJets, kFALSE);
312	DoJetLoop(fMCJets, fJets, kTRUE);
313
314	return kTRUE;
315	}
316
2949a2ec	317	//________________________________________________________________________
44476be7	318	void AliJetResponseMaker::DoJetLoop(TClonesArray jets1, TClonesArray jets2, Bool_t mc)
2949a2ec	319	{
44476be7	320	// Do the jet loop.
1ee1b5b8	321
44476be7	322	Int_t nJets1 = jets1->GetEntriesFast();
44476be7	323	Int_t nJets2 = jets2->GetEntriesFast();
1ee1b5b8	324
44476be7	325	for (Int_t i = 0; i < nJets1; i++) {
1ee1b5b8	326
44476be7	327	AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
2949a2ec	328
44476be7	329	if (!jet1) {
	330	AliError(Form("Could not receive jet %d", i));
	331	continue;
	332	}
	333
b3ccdfe2	334	if (!AcceptJet(jet1, kTRUE, !mc))
44476be7	335	continue;
	336
	337	for (Int_t j = 0; j < nJets2; j++) {
	338
	339	AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
	340
	341	if (!jet2) {
	342	AliError(Form("Could not receive jet %d", j));
	343	continue;
	344	}
	345
b3ccdfe2	346	if (!AcceptJet(jet2, kTRUE, mc))
44476be7	347	continue;
	348
	349	Double_t deta = jet2->Eta() - jet1->Eta();
	350	Double_t dphi = jet2->Phi() - jet1->Phi();
	351	Double_t d = TMath::Sqrt(deta * deta + dphi * dphi);
	352
	353	if (d < jet1->ClosestJetDistance()) {
	354	jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
	355	jet1->SetClosestJet(jet2, d);
	356	}
	357	else if (d < jet1->SecondClosestJetDistance()) {
	358	jet1->SetSecondClosestJet(jet2, d);
	359	}
	360	}
	361	}
1ee1b5b8	362	}
	363
	364
	365	//________________________________________________________________________
	366	Bool_t AliJetResponseMaker::FillHistograms()
	367	{
	368	// Fill histograms.
	369
2bddb6ae	370	const Int_t nMCJets = fMCJets->GetEntriesFast();
2949a2ec	371
2bddb6ae	372	for (Int_t i = 0; i < nMCJets; i++) {
2949a2ec	373
e44e8726	374	AliEmcalJet* jet = static_cast<AliEmcalJet*>(fMCJets->At(i));
2949a2ec	375
2949a2ec	376	if (!jet) {
1b76c28f	377	AliError(Form("Could not receive jet %d", i));
2949a2ec	378	continue;
	379	}
	380
2bddb6ae	381	if (!AcceptJet(jet, kTRUE, kFALSE))
	382	continue;
	383
	384	if (jet->Pt() > fMaxBinPt)
2949a2ec	385	continue;
2949a2ec	386
1c59613f	387	if (jet->ClosestJet() && jet->ClosestJet()->ClosestJet() == jet &&
1c59613f	388	jet->ClosestJetDistance() < fMaxDistance) { // Matched jet found
1b76c28f	389	jet->SetMatchedToClosest();
1b76c28f	390	jet->ClosestJet()->SetMatchedToClosest();
2bddb6ae	391	if (jet->MatchedJet()->Pt() > fMaxBinPt) {
1ee1b5b8	392	fHistMissedMCJets->Fill(jet->Pt(), fEventWeight);
2bddb6ae	393	}
2bddb6ae	394	else {
1ee1b5b8	395	fHistClosestDistance->Fill(jet->ClosestJetDistance(), fEventWeight);
1ee1b5b8	396
2bddb6ae	397	Double_t deta = jet->MatchedJet()->Eta() - jet->Eta();
1ee1b5b8	398	fHistClosestDeltaEta->Fill(deta, fEventWeight);
1ee1b5b8	399
2bddb6ae	400	Double_t dphi = jet->MatchedJet()->Phi() - jet->Phi();
1ee1b5b8	401	fHistClosestDeltaPhi->Fill(dphi, fEventWeight);
1ee1b5b8	402
2bddb6ae	403	Double_t dpt = jet->MatchedJet()->Pt() - jet->Pt();
1ee1b5b8	404	fHistClosestDeltaPt->Fill(dpt, fEventWeight);
	405
	406	fHistPartvsDetecPt->Fill(jet->MatchedJet()->Pt(), jet->Pt(), fEventWeight);
2bddb6ae	407	}
2949a2ec	408	}
99c67c1b	409	else {
1ee1b5b8	410	fHistNonMatchedMCJetPt->Fill(jet->Pt(), fEventWeight);
1ee1b5b8	411	fHistMissedMCJets->Fill(jet->Pt(), fEventWeight);
99c67c1b	412	}
2949a2ec	413
1ee1b5b8	414	fHistMCJetsPt->Fill(jet->Pt(), fEventWeight);
2949a2ec	415
1ee1b5b8	416	fHistMCJetPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
2949a2ec	417
2949a2ec	418	if (fAnaType == kEMCAL)
1ee1b5b8	419	fHistMCJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight);
2949a2ec	420
2949a2ec	421	for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
2bddb6ae	422	AliVParticle *track = jet->TrackAt(it, fMCTracks);
2949a2ec	423	if (track)
1ee1b5b8	424	fHistMCJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt(), fEventWeight);
1b76c28f	425	}
	426	}
	427
2bddb6ae	428	const Int_t nJets = fJets->GetEntriesFast();
1b76c28f	429
2bddb6ae	430	for (Int_t i = 0; i < nJets; i++) {
1b76c28f	431
e44e8726	432	AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
1b76c28f	433
	434	if (!jet) {
	435	AliError(Form("Could not receive mc jet %d", i));
2949a2ec	436	continue;
1b76c28f	437	}
99c67c1b	438
2bddb6ae	439	if (!AcceptJet(jet))
1b76c28f	440	continue;
2949a2ec	441
99c67c1b	442	if (!jet->MatchedJet()) {
1ee1b5b8	443	fHistNonMatchedJetPt->Fill(jet->Pt(), fEventWeight);
99c67c1b	444	}
1b76c28f	445
1ee1b5b8	446	fHistJetsPt->Fill(jet->Pt(), fEventWeight);
1b76c28f	447
1ee1b5b8	448	fHistJetPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
1b76c28f	449
1b76c28f	450	if (fAnaType == kEMCAL)
1ee1b5b8	451	fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight);
1b76c28f	452
1b76c28f	453	for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
2bddb6ae	454	AliVParticle *track = jet->TrackAt(it, fTracks);
1b76c28f	455	if (track)
1ee1b5b8	456	fHistJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt(), fEventWeight);
2bddb6ae	457	}
	458
	459	for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
	460	AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
	461	if (cluster) {
	462	TLorentzVector nP;
	463	cluster->GetMomentum(nP, fVertex);
1ee1b5b8	464	fHistJetsZvsPt->Fill(nP.Pt() / jet->Pt(), jet->Pt(), fEventWeight);
2bddb6ae	465	}
1b76c28f	466	}
1b76c28f	467	}
6fd5039f	468
6fd5039f	469	return kTRUE;
1b76c28f	470	}